# Multi-segmented non-isothermal compositional liquid gas well model for geothermal processes

Daniel Castanon Quiroz, Laurent Jeannin, Simon Lopez, Roland Masson

► To cite this version:

Daniel Castanon Quiroz, Laurent Jeannin, Simon Lopez, Roland Masson. Multi-segmented non-isothermal compositional liquid gas well model for geothermal processes. SMAI Journal of Computational Mathematics, 2025, 11, pp.203 - 231. 10.5802/smai-jcm.123. hal-05454736

HAL Id: hal-05454736

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

Submitted on 12 Jan 2026

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 CC BY 4.0 - Attribution - International License
Multi-segmented non-isothermal
compositional liquid gas well model
for geothermal processes

DANIEL CASTANON QUIROZ, LAURENT JEANNIN,
SIMON LOPEZ & ROLAND MASSON

Volume 11 (2025), p. 203-231.

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

© The authors, 2025.

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

CENTRE
MERSENNE

The SMAI Journal of Computational Mathematics is a member
of the Centre Mersenne for Open Scientific Publishing

http://www.centre-mersenne.org/

Submissions at https://smai-jcm.centre-mersenne.org/ojs/submission

e-ISSN: 2426-8399
<figure><img src="image_3.png" /></figure>

# Multi-segmented non-isothermal compositional liquid gas well model for geothermal processes

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

LAURENT JEANNIN<fnref n="2" />

SIMON LOPEZ<fnref n="3" />

ROLAND MASSON<fnref n="4" />

**Abstract.** We consider a non-isothermal compositional gas liquid model for the simulation of well operations in geothermal processes. The model accounts for phase transitions assumed to be at thermodynamical equilibrium and is based on a hydrodynamical Drift Flux Model (DFM) combined with a No Pressure Wave approximation of the momentum equation. The focus of this work is on the design of a robust discretization accounting for slanted and multibranch wells with the ability to simulate both transient behavior such as well opening as well as coupled simulations at the time scale of the reservoir. It is based on a staggered finite volume scheme in space combined with a fully implicit Euler time integration. The construction of consistent and stable numerical fluxes is a key feature for a robust numerical method. This is achieved by combining a monotone flux approximation for the phase superficial velocities with an upwind approximation of the phase molar fractions, density and enthalpy. In order to facilitate the coupling of the well and reservoir models, the Newton linearization accounts for the elimination of the hydrodynamical unknowns leading to Jacobian systems using the same primary unknowns than those of the reservoir model. The efficiency of our approach is investigated on both stand alone well test cases without and with cross flow, and on a fully coupled well-reservoir simulation.

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

**Keywords.** thermal well model, drift flux model, liquid gas compositional model, geothermal system, multi-segmented well, finite volume scheme, staggered finite volume, monotone flux, Coats formulation.

## 1. Introduction

Wells are key objects during the operation of geothermal reservoirs: they provide the connection between the surface and the geothermal reservoir at depth, enable fluid and heat exchanges with the geological structure, and are among the places of two-phase flows and phase changes. During appraisal phases well tests provide valuable information on the target and during geothermal field operation, it is very important to correctly describe the behavior of wellbore flows in order to predict the quality

<footnote n="1">¹ 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)</footnote>

<footnote>E-mail address: daniel.castanon@iimas.unam.mx</footnote>

<footnote n="2">² STORENGY, 12 rue Raoul Nordling - Djinn - CS 70001 92274 Bois Colombes Cedex, France</footnote>

<footnote>E-mail address: laurent.jeannin@storengy.com</footnote>

<footnote n="3">³ BRGM, 3 avenue Claude-Guillemin, BP 36009, 45060 Orléans Cedex 2, France</footnote>

<footnote>E-mail address: s.lopez@brgm.fr</footnote>

<footnote n="4">⁴ Université Côte d'Azur, Inria, CNRS, LJAD, UMR 7351 CNRS, team GALETS, Parc Valrose 06108 Nice Cedex 02, France</footnote>

<footnote>E-mail address: roland.masson@univ-cotedazur.fr.</footnote>

<footnote>The work of Daniel Castanon Quiroz was partially supported by UNAM-PAPIIT grant IA-101723. The authors also acknowledge the support of the French Geological Survey (BRGM) and Storengy.</footnote>

<footnote>https://doi.org/10.5802/smai-jcm.123</footnote>

<footnote>© The authors, 2025</footnote>
of the fluids produced at the surface. A thorough understanding of wells and their modeling allow a better understanding of the reservoir behavior.

Physical phenomena and exchanges with the reservoir are complex. For example, the start-up of a high-temperature geothermal well, fed by a hot liquid, gives rise to transient flows with flash in the well. Geothermal wells may also be connected by several feed zones at different pressures, and some feed zones will be produced in the well, while others may receive fluids from the well, even when the latter is in production. Finally, in addition to describing wells over transient short periods, well modeling can also be used to predict the operational behavior of a reservoir over a period of several years. In this case, the well model has to be coupled with a flow model in the geological structures.

Flows occurring in geothermal wells are generally described as two-phase thermal pipe flows. Given the well's radius/length aspect ratio, the usual assumption is to rely on 1D models. This type of flow is relevant not only to the geothermal industry, but also to the oil and gas industry (for example, [25] for steady flows, [8] for transient flows or [13]). The nuclear industry has also a strong interest in multiphase thermal flows in pipes ([9] for a review). We present below a few bibliographical elements relevant to our work for geothermal processes, but readers interested in more details may find further information in [28], where Tonkin et al. recently proposed an in-depth review of mathematical models of flows in geothermal boreholes.

The specifics of geothermal flows lie in the description of the interaction between the well and the reservoir, and the consideration of feed zones and cross-flow phenomena. Bjornsson [6] was one of the first to propose a multi-feedzone well simulator. These well models are either transient or steady-state; transient well models are particularly useful for interpreting well tests or, for example, well start-up, while steady-state models are more dedicated to operations follow-up (for example, [11] for matching field data or [3] for the modeling of the IDDP-1 supercritical well). Let us also refer to [15, 16, 19, 27] for the development of numerical non-isothermal multi-segmented well models and their coupling with the reservoir.

One of the greatest difficulties in modeling multiphase flows in geothermal wells arises from the complexity and variety of flow regimes encountered. The simplest homogeneous flow model, in which all phases flow at the same velocity, is unsuitable because it cannot satisfactorily reproduce the in-situ well volume fraction and flow rate of each phase. Drift Flux Models (DFM) are more complex models in which a drift velocity (or slip law) is introduced to describe the relative motion between phases [30, 33]. These DFM require a number of empirical correlations derived from experiments for different flow regimes in order to correctly describe phase slip and in situ phase volume fractions [9, 28, 30]. Given the diversity of flow regimes encountered in geothermal wells, preference is given to DFM capable of describing all flow regimes encountered during production well operations [16, 19, 24, 26, 28].

It is usually assumed that the thermodynamical equilibrium assumption applies in the well, but some authors are interested in non-equilibrium phenomena [2, 22] to describe scaling or kinetics phenomena [14, 17, 23, 32]. Nevertheless, these reactive phenomena will not be considered in our work. Our objective is to set-up an extensible framework for the modeling of transient multi-component two-phase geothermal wellbore flows with a focus on robust numerical schemes and the ability to couple reservoir and wellbore flow models.

The considered well thermodynamical model accounts for a two-phase liquid gas compositional non-isothermal system as for the reservoir model described in [31]. This two-phase compositional system uses a Coats-type formulation in order to account for a large class of equilibrium thermodynamical models ranging from single phase, two-phase immiscible to two-phase fully compositional with phase appearance and disappearance. Choosing the same formulation for both the well and the reservoir ensures robust coupling between these two systems. Then, following [4], to be able to consider complex well geometries with slanted or multi-branched paths, the geometry of each well is defined as a set of edges of the reservoir mesh assumed to define a rooted tree oriented away from the root.
The hydrodynamical model is based on the DFM introduced in [24] for the slip law. Moreover, following [16, 18, 29], we consider a “No Pressure Wave” (NPW) approximation, where the mixture momentum equation is reduced to a static balance accounting for gravity and wall friction pressure losses. This simplified NPW model has been shown to provide good results in pipe flow simulations even in the highly transient case of severe slugging [18, 29] as well as in typical reservoir simulations as stated in [16]. It is possible, however, that for some flow regimes the kinetic terms will not be negligible compared with the gravity and wall friction pressure drop terms. In that case our approach could be extended to incorporate these additional terms, but at the price of giving up the elimination of the edge mixture velocities from the Jacobian system at each nonlinear iteration. This investigation goes beyond the present work.

The discretization is based on a staggered finite volume scheme with edge center superficial velocities and nodal pressure, temperature, saturation and molar fraction unknowns. This choice is consistent with the Vertex Approximate Gradient (VAG) nodal discretization [31] used in the reservoir. It is combined with a two-point monotone flux for the phase superficial velocities and a phase based upwind approximation of the phase molar fractions, molar density and molar enthalpy. The monotonicity of the superficial velocity fluxes is a key property for the stability of the scheme. Moreover, it guarantees the consistent definition of the phase based upwinding according to the sign of the phase superficial velocity. The clear understanding of this interplay between the upwinding of the phase superficial velocities and of the upwinding of the phase transported variables is one of the main contribution of this work. The monotonicity property of the superficial velocities is obtained on the full range of saturations [0, 1] using a hybrid upwinding approach exploiting the mathematical properties of the drift velocity, and profile parameters of the DFM model proposed in [24]. Note that this hybrid upwinding methodology could be applied to other slip laws of generalized Zuber-Findlay type.

The time integration is fully implicit to account for the large reservoir time scale and for the stiffness of the system at the gas phase appearance. The Jacobian system is reduced by the elimination of the phase superficial velocities and of the well flow rates using the slip law, the momentum NPW equation and the well monitoring conditions. This elimination allows the use of the same primary unknowns as for the reservoir model [31] which facilitates the fully implicit coupling with the reservoir system.

The remainder of this paper is organized as follows. Section 2 presents the well physical model, and its staggered finite volume discretization is described in Section 3. It starts, in Subsection 3.1, by the definition of the well mesh and discrete unknowns, while the staggered finite volume scheme is detailed in Subsection 3.2. Subsection 3.3 defines the monotonicity properties required for the phase superficial velocities and underlines the interplay between these properties and the upwind approximation of the phase molar fractions, the density and the enthalpy. The detailed specific construction of the monotone superficial velocity flux function based on a hybrid upwinding approach for the DFM model proposed in [24] is postponed to Appendix A. The well monitoring conditions are introduced in Subsection 3.4 in the case of a production well and Appendix B recalls the definition of the reservoir-well fluxes based on the Peaceman indexes. Then, the algorithm used to solve the coupled nonlinear system at each time step of the simulation is addressed in Section 4. It is based on a Newton-Raphson algorithm combined with an active set method for the phase appearance/disappearance and the well monitoring conditions. Its main originality is related to the elimination in the linearization process of the thermodynamical, hydrodynamical and monitoring closure laws leading to the same set of primary unknowns as for the reservoir model in the Jacobian system. This procedure is described respectively in Subsections 4.1, 4.2 and 4.3. The objective of the numerical Section 5 is to validate the discrete model and to investigate the efficiency of our approach on both stand alone well test cases in Subsection 5.1 and for a fully coupled well-reservoir model in Subsection 5.2. The first test case in Subsection 5.1.1 considers the validation of the numerical model on a simple two-phase incompressible immiscible flow which reduces to a Buckley Leverett scalar hyperbolic equation. The second and third test cases consider a high energy liquid vapor single component thermal flow along a multi-branch production well without and
with cross flow in respectively Subsections 5.1.2 and 5.1.3. Finally, a single well liquid vapor geothermal test case is considered in Subsection 5.2 to validate the fully coupled model by comparison to a simpler well model based on a single implicit well unknown.

## 2. Well physical model

Let $W$ denote the set of wells. As in [4], each multi-branch well $\omega \in W$ is defined by a set of oriented edges of the reservoir mesh assumed to define a rooted tree oriented away from the root. This orientation corresponds to the drilling direction of the well. We consider a two-phase liquid gas, compositional, and non-isothermal flow model. The liquid ($\ell$) and gas (g) phases are described by their pressure $p$ (both liquid and gas pressures are assumed to match along the well), temperature $T$, volume fraction or saturation $s^\alpha$, and molar fractions $c^\alpha = (c_i^\alpha)_{i \in C}$, $\alpha \in P = \{\ell, g\}$, where $C$ denotes the set of components and $P$ the set of liquid and gas phases. The components are not necessarily present in both phases but at least in one of them. We consequently denote by $C^\alpha \subset C$ the set of components in phase $\alpha$ and by $P_i$ the set of phases containing the component $i \in C$.

The thermodynamical model uses a Coats formulation already developed for the reservoir model (see e.g. [31]). It is based on the natural variables $p, T, s^\alpha, c^\alpha, \alpha \in Q$ with a set of unknowns and equations depending on the additional unknown $Q \subset P$ representing the subset of present phases at each point of the space time domain. This formulation has the advantage to account for an arbitrary number of phases and components and also to allow the components to be either present or absent of each given phase. This allows to account in the same framework for models ranging from single phase liquid, single phase gas to two-phase gas liquid, and from immiscible to fully compositional.

For each phase $\alpha$, we denote by $\zeta^\alpha(p, T, c^\alpha)$ its molar density, by $\rho^\alpha(p, T, c^\alpha)$ its specific density, by $\mu^\alpha(p, T, c^\alpha)$ its dynamic viscosity, by $e^\alpha(p, T, c^\alpha)$ its molar internal energy, and by $h^\alpha(p, T, c^\alpha)$ its molar enthalpy. The fugacities are denoted by $f_i^\alpha(p, T, c^\alpha)$.

The hydrodynamical model is based on the Drift Flux Model (DFM) with slip closure laws expressing the phase superficial velocities as functions of the mixture velocity and drift velocity terms [34]. It is combined with the No-Pressure-Wave (NPW) approximation of the momentum equation relating the pressure gradient to the friction and gravity pressure losses [16, 29].

Let $\tau$ denote the curvilinear coordinate along the well and $e_\tau$ the unit tangential vector along the well assumed to be oriented toward the well root node. We denote by $\mathbf{u}^\alpha = u^\alpha e_\tau$ the superficial velocity of each phase $\alpha \in P$. The mixture velocity is defined by

$$ \mathbf{u}^m = \mathbf{u}^\ell + \mathbf{u}^g = u^m e_\tau. $$

The tangential divergence along the well is denoted by $\mathrm{div}_\tau$ and the tangential gradient by $\nabla_\tau$. The section along the well is denoted by $S^\omega(\tau)$. We also set $g_\tau = (g \cdot e_\tau)e_\tau$ where $g$ is the acceleration of gravity vector.

Let $\epsilon_\omega(\tau) = \frac{(\nabla_\tau z(\tau)) \cdot e_\tau}{|\nabla_\tau z(\tau)|} \in \{-1, 1\}$ define the orientation along the well, with arbitrary value for $\nabla_\tau z = 0$ along a horizontal branch. Following [24], let us introduce the function

$$ o_\omega(\tau) = \epsilon_\omega(\tau) \cos(\theta_\omega(\tau))^{1/2} (1 + \sin(\theta_\omega(\tau)))^2, $$

where $\theta_\omega(\tau) \in [0, \pi/2]$ is the edge angle w.r.t. the vertical direction. Note that $o_\omega(\tau) = 0$ along a horizontal branch. The slip closure law introduced in [34] expresses the gas superficial velocity $u^g$ as a function of the mixture velocity, the drift velocity, and the profile parameter:

$$ u^g = s^g U_d(s^g, \bar{H}) o_\omega + s^g C_0(s^g, \bar{H}) u^m, \quad (2.1) $$

with the non-negative drift velocity $U_d$ (for vertical wells) and the profile parameter $C_0$ such that $C_0(s^g = 1, \bar{H}) = 1$, $U_d(s^g = 1, \bar{H}) = 0$, $s^g C_0(s^g, \bar{H}) \le 1$ and $s^g C_0(s^g, \bar{H})$ non-decreasing w.r.t. $s^g$. The drift velocity and profile parameter depend on a set of thermodynamical variables denoted by $\bar{H}$.
and typically comprising the densities of the phases and the liquid gas interfacial tension. They also depend on other fixed parameters such as the well radius $r_\omega$.

We denote by

$$ \mathbf{T}^f = T^f \mathbf{e}_\tau, $$

the wall friction law depending on the mixture velocity $\mathbf{u}^m$ and on thermodynamical quantities. As a typical example, we will use in the following a Darcy–Forchheimer type law:

$$ T^f = - \left( \frac{8\mu^m}{r_\omega^2} + f_q \rho^m \frac{1}{4r_\omega} |u^m| \right) u^m, \qquad (2.2) $$

where $f_q$ is the friction coefficient, $\mu^m = s^g\mu^g + s^l\mu^l$ the mean viscosity and $\rho^m = s^g\rho^g + s^l\rho^l$ the mean specific density.

The well model is based on the molar conservation of each component (2.3a), the total energy conservation (2.3b) and the mixture momentum equation based on the NPW model (2.3c). It is combined with the DFM slip law (2.3e) together with the thermodynamical equilibrium and the sum to one of the saturations (2.3d):

$$ \sum_{\alpha \in \mathcal{P}_i} \left( \partial_t (S^\omega c_i^\alpha \zeta^\alpha s^\alpha) + \mathrm{div}_\tau (S^\omega c_i^\alpha \zeta^\alpha \mathbf{u}^\alpha) \right) = q_i^{r \to \omega}, \quad i \in \mathcal{C}, \qquad (2.3a) $$

$$ \sum_{\alpha \in \mathcal{P}} \left( \partial_t (S^\omega e^\alpha \zeta^\alpha s^\alpha) + \mathrm{div}_\tau (S^\omega h^\alpha \zeta^\alpha \mathbf{u}^\alpha) \right) + \mathrm{div}_\tau (-S^\omega \lambda \nabla_\tau T) = q_e^{r \to \omega}, \qquad (2.3b) $$

$$ \nabla_\tau p = \mathbf{T}^f + \rho^m \mathbf{g}_\tau, \qquad (2.3c) $$

$$ \begin{gathered} f_i^g(p, T, c^g) = f_i^l(p, T, c^l) \quad \text{if } Q \cap \mathcal{P}_i = \mathcal{P}, \\ \sum_{i \in \mathcal{C}^\alpha} c_i^\alpha = 1, \quad \alpha \in Q, \\ s^\alpha = 0 \quad \text{if } \alpha \notin Q, \end{gathered} \qquad (2.3d) $$

$$ \sum_{\alpha \in Q} s^\alpha = 1 $$

$$ u^g = s^g U_d(s^g, \bar{H}) o_\omega + s^g C_0(s^g, \bar{H}) u^m \quad \text{if } Q = \mathcal{P}, \qquad (2.3e) $$

where $q_i^{r \to \omega}$ and $q_e^{r \to \omega}$ are the molar and energy exchange terms with the reservoir. The system is closed with monitoring conditions at the root node and with a flash computation to determine the set of present phases $Q$. The temperature and pressure continuity is also assumed at the well junctions. Note that here the term flash is used in a broad sense of an algorithm to manage phase appearance and disappearance which implementation depends on the considered thermodynamical laws. It can be typically a negative flash or a stability analysis for phase appearance, and the sign of the saturation for phase disappearance.

### 3. Staggered Finite Volume discretization of the well model

Notations about the well mesh and the discrete variables are introduced in Subsection 3.1. The scheme presented in Subsection 3.2 is based on a fully implicit time integration to cope with large time steps at the reservoir time scale. It is combined with a staggered finite volume discretization in space using node centred control volumes for the molar and energy conservation equations and edge control volumes for the momentum equation. A key ingredient is the discretization of the convective fluxes based on a monotone two-point flux for the approximation of the superficial velocities. This framework is presented in Subsection 3.3 and an example is detailed in Appendix A based on the model proposed in [24]. Thanks to the monotonicity and consistency properties of the superficial velocities, they can
be combined with an upwind approximation of the phase molar fractions, density and enthalpy w.r.t. the sign of the phase superficial velocity. The monitoring conditions at the head node of the well are described in Subsection 3.4 in the case of a production well considered in the numerical Section.

### 3.1. Well discretization and notations

The set of nodes of a well $\omega \in W$ is denoted by $V_\omega$ and its root node is denoted by $v_\omega$. A partial ordering is defined on the set of vertices $V_\omega$ with $v <_\omega v'$ if and only if the unique path from the root $v_\omega$ to $v'$ passes through $v$. The set of edges of the well $\omega$ is denoted by $\mathcal{E}_\omega$ and for each edge $a \in \mathcal{E}_\omega$ we set $a = vv'$ with $v <_\omega v'$ (i.e. $v$ is the parent node of $v'$, see Figure 3.1). It is assumed that $V_{\omega_1} \cap V_{\omega_2} = \emptyset$ for any $\omega_1, \omega_2 \in W$ such that $\omega_1 \neq \omega_2$.

Let $|S_a^\omega|$ and $|S_{v_\omega}^\omega|$ denote the well section at respectively the center of the edge $a$ and at the head node $v_\omega$. For $a \in \mathcal{E}_\omega$, $|a|$ denotes the length of the edge $a$. Let us define $\mathcal{E}_V^\omega \subset \mathcal{E}_\omega$ as the set of well edges sharing the node $v \in V_\omega$. For all $vv' = a \in \mathcal{E}_\omega$, let us set $\kappa_{a,v'} = -1$ and $\kappa_{a,v} = 1$.

<figure><img src="image_4.png" /><figcaption>FIGURE 3.1. Example of multi-branch well $\omega$ with its root node $v_\omega$, one edge $a = vv'$ ($v$ parent node of $v'$) and the main physical quantities: the head node superficial velocities $u_\omega^\alpha$ (non-negative for production wells and non-positive for injection wells), the molar and energy flow rates between the reservoir and the well $q_{v,i}^{r\to\omega}$, $q_{v,e}^{r\to\omega}$, the well node pressure, temperature, saturations and phase molar fractions $p_v, T_v, s_v^\alpha, c_v^\alpha$, and the phase superficial velocity $u_a^\alpha$ at the edge $a$ oriented positively from $v'$ to $v$.</figcaption></figure>

For each edge $a = vv' \in \mathcal{E}_\omega$, let us denote by $u_a^\alpha$ the superficial velocity of phase $\alpha$ along the edge $a$ oriented positively from $v'$ to $v$. Using the Coats formulation, the set of well unknowns is defined at each node $v \in V_\omega$ by

* the set of present phases $Q_v$,

* the well pressure $p_v$,
• the well temperature $T_v$,

• the well saturations $s_v^\alpha$ for $\alpha \in \mathcal{P}$,

• the well molar fractions $c_v^\alpha$ for $\alpha \in Q_v$,

• the number of moles $\tilde{n}_{v,i}$ for $i \in \tilde{C}_{Q_v}$,

and by

• the edge superficial velocities $u_a^\alpha$ for $\alpha \in \mathcal{P}$ at each well edge $a \in \mathcal{E}_\omega$,

• the head node superficial velocities $u_\omega^\alpha$ for $\alpha \in \mathcal{P}$ at the well head node $v_\omega$.

In the above definition, $\tilde{C}_Q$ denotes the set of components not contained in $\cup_{\alpha \in Q} C^\alpha$. For $\tilde{C}_{Q_v} \neq \emptyset$, the additional unknowns $\tilde{n}_{v,i}$ corresponding to the number of moles of the absent components $i \in \tilde{C}_{Q_v}$ are needed first to avoid the singularity of the Jacobian and second to track the appearance of the missing phase containing these absent components (see [31] for details). Note that the saturation $s_v^\alpha$ of an absent phase $\alpha \notin Q_v$ vanishes. In the following, we assume that the monitoring conditions at the well head node correspond to the case of a production well in line with the numerical section.

For a given thermodynamical phase property $\xi^\alpha$ at a node $v$, we will use the notation

$$ \xi_v^\alpha = \xi^\alpha(p_v, T_v, c_v^\alpha), $$

with typically $\xi^\alpha = \zeta^\alpha, \rho^\alpha, e^\alpha, h^\alpha, f_i^\alpha, \mu^\alpha$.

The molar and energy flow rates between the reservoir and the well $\omega$ at a given node $v \in \mathcal{V}_\omega$ are defined by a two point flux approximation between the reservoir, and the well properties at node $v$ based on the Peaceman approach combined with a phase potential upwinding. They are denoted respectively by $q_{v,i}^{r \to \omega}$ for each component $i \in C$ and $q_{v,e}^{r \to \omega}$ in what follows (see Figure 3.1 and Appendix B for their definitions).

## 3.2. Fully implicit Finite Volume scheme

Let us define at each nodal control volume $v$, the total number of mole of each component and the total energy as follows

$$ n_{v,i} = \begin{cases} \left( \sum_{a \in \mathcal{E}_v^\omega} |S_a^\omega| \frac{|a|}{2} \right) \sum_{\alpha \in Q_v \cap \mathcal{P}_i} c_{v,i}^\alpha s_v^\alpha \zeta_v^\alpha & \text{if } i \notin \tilde{C}_{Q_v}, \\ \tilde{n}_{v,i} & \text{otherwise,} \end{cases} $$

$$ n_{v,e} = \left( \sum_{a \in \mathcal{E}_v^\omega} |S_a^\omega| \frac{|a|}{2} \right) \sum_{\alpha \in Q_v} s_v^\alpha \zeta_v^\alpha e_v^\alpha, $$

where the case of an absent component $i \in \tilde{C}_{Q_v}$ is accounted for by the introduction of the additional unknown $\tilde{n}_{v,i}$ in order to avoid the singularity of the discrete system.

For each edge $a = vv' \in \mathcal{E}_\omega$, and each phase $\alpha$, let us define the following upwind approximation of the phase molar fractions, density and enthalpy w.r.t. the sign of the phase superficial velocity:

$$ \xi_a^\alpha = \begin{cases} \xi_{v'}^\alpha & \text{if } u_a^\alpha \ge 0, \\ \xi_v^\alpha & \text{if } u_a^\alpha < 0, \end{cases} \qquad (3.1) $$

with $\xi^\alpha = c^\alpha, \zeta^\alpha, h^\alpha$.
We denote by $(0, t_F)$ the time interval and consider its discretization given by $t^n$, $n = 0, \dots, N$ with $t^0 = 0$, $t^N = t_F$ and $\Delta t^n = t^n - t^{n-1} > 0$ for all $n = 1, \dots, N$. The time discretization of the system (2.3) is based on an implicit Euler integration scheme. To simplify the notations only the accumulation terms $n_{v,i}$, $n_{v,e}$ at time $t^{n-1}$ will be specified by the $n-1$ superscript. For all other quantities considered at the current time $t^n$, we will drop the $n$ superscript for simplicity.

Then, the discretization of the well equations (2.3a)-(2.3b)-(2.3c)-(2.3d) at each time step $n = 1, \dots, N$ are respectively defined by the equations (3.2a)-(3.2b)-(3.2c)-(3.2d) as follows:

$$ \frac{n_{v,i} - n_{v,i}^{n-1}}{\Delta t^n} + \sum_{a \in \mathcal{E}_{v}^{\omega}} \sum_{\alpha \in \mathcal{P}_i} -\kappa_{a,v} |S_{a}^{\omega}| c_{a,i}^{\alpha} \zeta_{a}^{\alpha} u_{a}^{\alpha} \\ = q_{v,i}^{r \to \omega} - \delta_{v}^{\nu_{\omega}} \sum_{\alpha \in \mathcal{P}_i} |S_{v_{\omega}}^{\omega}| c_{v,i}^{\alpha} \zeta_{v}^{\alpha} u_{\omega}^{\alpha}, \quad i \in \mathcal{C}, \quad v \in \mathcal{V}_{\omega}, \qquad (3.2a) $$

$$ \frac{n_{v,e} - n_{v,e}^{n-1}}{\Delta t^n} + \sum_{a \in \mathcal{E}_{v}^{\omega}} |S_{a}^{\omega}| \left( \frac{\lambda_{a}}{|a|} (T_{v} - T_{v'}) + \sum_{\alpha \in \mathcal{P}} -\kappa_{a,v} h_{a}^{\alpha} \zeta_{a}^{\alpha} u_{a}^{\alpha} \right) \\ = q_{v,e}^{r \to \omega} - \delta_{v}^{\nu_{\omega}} \sum_{\alpha \in \mathcal{P}} |S_{v_{\omega}}^{\omega}| h_{v}^{\alpha} \zeta_{v}^{\alpha} u_{\omega}^{\alpha}, \quad v \in \mathcal{V}_{\omega}, \qquad (3.2b) $$

$$ p_{v} - p_{v'} = -\bar{\rho}_{a}^{m} g(z_{v} - z_{v'}) + T_{a}^{f} |a|, \qquad v v' = a \in \mathcal{E}_{\omega}, \qquad (3.2c) $$

$$ \begin{aligned} f_{i,v}^{g} &= f_{i,v}^{l} \quad \text{if } Q_{v} \cap \mathcal{P}_{i} = \mathcal{P}, \quad v \in \mathcal{V}_{\omega}, \\ \sum_{i \in \mathcal{C}^{\alpha}} c_{v,i}^{\alpha} &= 1, \quad \alpha \in Q_{v}, \quad v \in \mathcal{V}_{\omega}, \\ s_{v}^{\alpha} &= 0 \quad \text{if } \alpha \notin Q_{v}, \quad v \in \mathcal{V}_{\omega}, \\ \sum_{\alpha \in Q_{v}} s_{v}^{\alpha} &= 1, \quad v \in \mathcal{V}_{\omega}. \end{aligned} \qquad (3.2d) $$

where $\delta$ stands for the Kronecker symbol with $\delta_{v}^{\nu_{\omega}} = 1$ if $v = v_{\omega}$, else 0. In (3.2), the edge thermal conductivity of the mixture is defined by

$$ \lambda_a = \sum_{\alpha \in \mathcal{P}} \frac{s_{v}^{\alpha} + s_{v'}^{\alpha}}{2} \lambda^{\alpha}, $$

with $\lambda^\alpha$ the thermal conductivity of the phase $\alpha$ assumed constant for simplicity. The wall friction term $T_a^f$ is given by the Darcy–Forchheimer law (2.2) leading to

$$ T_a^f = - \left( \frac{8 \bar{\mu}_a^m}{r_a^2} + f_q \bar{\rho}_a^m \frac{1}{4 r_a} |u_a^m| \right) u_a^m, \qquad (3.3) $$

with the following arithmetic means of the mixture specific density and dynamic viscosity at the edge a

$$ \bar{\rho}_a^m = \frac{1}{2} \left( \sum_{\alpha \in Q_v} s_v^\alpha \rho_v^\alpha + \sum_{\alpha \in Q_{v'}} s_{v'}^\alpha \rho_{v'}^\alpha \right), \quad \bar{\mu}_a^m = \frac{1}{2} \left( \sum_{\alpha \in Q_v} s_v^\alpha \mu_v^\alpha + \sum_{\alpha \in Q_{v'}} s_{v'}^\alpha \mu_{v'}^\alpha \right). \qquad (3.4) $$

It remains to define the discretization of the edge superficial velocities $u_a^\alpha$ accounting for the slip law (2.3e) at given mixture velocity $u_a^m = u_a^g + u_a^l$. This framework is detailed in the next subsection based on a monotone two-point flux. Note that the edge mixture velocity $u_a^m$ will be derived from the momentum equation (3.2c) combined with the wall friction law (2.2), as detailed in Section 4.2.1.
### 3.3. Monotone two-point flux for the edge superficial velocities

Let us define the arithmetic average of a thermodynamical variable $\xi^\alpha$ for a given phase $\alpha$ at a well edge $\mathbf{a} = \mathbf{vv}'$ with $\alpha \in Q_\mathbf{v} \cup Q_{\mathbf{v}'}$ by

$$ \bar{\xi}_{\mathbf{a}}^{\alpha} = \begin{cases} \xi^{\alpha}(p_{\mathbf{v}}, T_{\mathbf{v}}, c_{\mathbf{v}}^{\alpha}) & \text{if } \alpha \in Q_{\mathbf{v}} \text{ and } \alpha \notin Q_{\mathbf{v}'}, \\ \xi^{\alpha}(p_{\mathbf{v}'}, T_{\mathbf{v}'}, c_{\mathbf{v}'}^{\alpha}) & \text{else if } \alpha \notin Q_{\mathbf{v}} \text{ and } \alpha \in Q_{\mathbf{v}'}, \\ \frac{s_{\mathbf{v}}^{\alpha} \xi^{\alpha}(p_{\mathbf{v}}, T_{\mathbf{v}}, c_{\mathbf{v}}^{\alpha}) + s_{\mathbf{v}'}^{\alpha} \xi^{\alpha}(p_{\mathbf{v}'}, T_{\mathbf{v}'}, c_{\mathbf{v}'}^{\alpha})}{s_{\mathbf{v}}^{\alpha} + s_{\mathbf{v}'}^{\alpha}} & \text{else if } s_{\mathbf{v}}^{\alpha} + s_{\mathbf{v}'}^{\alpha} > \epsilon, \\ \frac{1}{2} \left( \xi^{\alpha}(p_{\mathbf{v}}, T_{\mathbf{v}}, c_{\mathbf{v}}^{\alpha}) + \xi^{\alpha}(p_{\mathbf{v}'}, T_{\mathbf{v}'}, c_{\mathbf{v}'}^{\alpha}) \right) & \text{otherwise.} \end{cases} $$

This arithmetic averaging is used to compute the thermodynamical variables $\bar{H}_\mathbf{a}$ at the well edge $\mathbf{a} = \mathbf{vv}'$ which enter in the definition of the slip law (2.3e).

Then, the superficial velocities $u_\mathbf{a}^\alpha$ are obtained at given mixture velocity $u_\mathbf{a}^m$ using a numerical two-point monotone flux denoted by

$$ F_{\mathbf{a}}^{g}(s_{\mathbf{v}'}^{g}, s_{\mathbf{v}}^{g}, \bar{H}_{\mathbf{a}}), $$

for the continuous flux function

$$ f_{\mathbf{a}}(s^g) = s^g U_d(s^g, \bar{H}_{\mathbf{a}}) o_{\mathbf{a}} + s^g C_0(s^g, \bar{H}_{\mathbf{a}}) u_{\mathbf{a}}^m, $$

where $o_\mathbf{a}$ is the value of the function $o_\omega$ along the edge $\mathbf{a}$. Note that the flux function $f_\mathbf{a}(s^g)$ is related to the Buckley Leverett scalar hyperbolic equation obtained assuming immiscible liquid and gas phases with constant molar densities, a constant well section and no source terms, leading to the gas phase volume conservation equation

$$ \partial_t s^g + \operatorname{div}_\tau (f_\mathbf{a}(s^g) e_\tau) = 0, \quad (3.5) $$

combined with the incompressibility condition on the mixture velocity $\operatorname{div}_\tau(u^m e_\tau) = 0$. This motivates our choice of the numerical flux $F_\mathbf{a}^g$ which, for such scalar hyperbolic equation (3.5), must satisfy the following consistency and stability properties:

* consistency property:

$$ F_{\mathbf{a}}^{g}(s^{g}, s^{g}, \bar{H}_{\mathbf{a}}) = f_{\mathbf{a}}(s^{g}) \text{ for all } s^{g} \in [0, 1], $$

* monotonicity property which specifies that, for all $(s_{\mathbf{v}'}^g, s_{\mathbf{v}}^g) \in [0, 1] \times [0, 1]$, $F_\mathbf{a}^g(s_{\mathbf{v}'}^g, s_{\mathbf{v}}^g, \bar{H}_\mathbf{a})$ is non-decreasing w.r.t. its first argument $s_{\mathbf{v}'}^g$ and non-increasing w.r.t. its second argument $s_{\mathbf{v}}^g$.

Then the superficial velocities $u_\mathbf{a}^\alpha$ are defined by

$$ \left\{ \begin{array}{l} u_{\mathbf{a}}^{g} = F_{\mathbf{a}}^{g}(s_{\mathbf{v}'}^{g}, s_{\mathbf{v}}^{g}, \bar{H}_{\mathbf{a}}), \\ u_{\mathbf{a}}^{l} = u_{\mathbf{a}}^{m} - u_{\mathbf{a}}^{g}. \end{array} \right. \qquad (3.6) $$

From the consistency and monotonicity properties, it results that

* $s_{\mathbf{v}'}^g = 0 \Rightarrow u_\mathbf{a}^g = F^g(0, s_{\mathbf{v}}^g) \le 0$ for all $s_{\mathbf{v}}^g \in [0, 1]$,

* $s_{\mathbf{v}}^g = 0 \Rightarrow u_\mathbf{a}^g = F^g(s_{\mathbf{v}'}^g, 0) \ge 0$ for all $s_{\mathbf{v}'}^g \in [0, 1]$,

* $s_{\mathbf{v}'}^g = 1 \Rightarrow u_\mathbf{a}^l = u_\mathbf{a}^m - F^g(1, s_{\mathbf{v}}^g) \le 0$ for all $s_{\mathbf{v}}^g \in [0, 1]$,

* $s_{\mathbf{v}}^g = 1 \Rightarrow u_\mathbf{a}^l = u_\mathbf{a}^m - F^g(s_{\mathbf{v}'}^g, 1) \ge 0$ for all $s_{\mathbf{v}'}^g \in [0, 1]$,
which numerically accounts for the fact that gas cannot flow out of a liquid only node control volume, and that liquid cannot flow out of a gas only node control volume. These properties are key for the definition (3.1) of the upwind approximations of the phase molar fractions, phase molar density and phase molar enthalpy w.r.t. the sign of the phase superficial velocity. Indeed, they ensure that either the phase is present at the upwind node or that its superficial velocity vanishes ensuring that the product $\xi_a^\alpha u_a^\alpha$ is always properly defined, with $\xi^\alpha = \zeta^\alpha, c^\alpha, h^\alpha$.

Note also that if $Q_v \cup Q_{v'} = \{\alpha\}$ (single phase $\alpha$ present at both nodes), the properties of the numerical flux function imply that $u_a^\alpha = u_a^m$ and $u_a^\beta = 0$ for the absent phase $\beta$. It means that, in that cases, the thermodynamical variables $\bar{H}_a$ do not need to be computed which would have required a cumbersome extension of the molar fractions of the absent phase.

The specific construction of the two-point monotone flux function $F_a^g$ depends on the profile parameter $C_0$ and drift velocity $U_d$ laws. We propose in Appendix A a monotone two-point flux for the gas liquid DFM model introduced in [24] on the full range of gas saturation. Our methodology is based on a hybrid upwinding strategy [10, 12] using a physical splitting of the flux function $f_a(s^g) = f_a^m(s^g) + f_a^d(s^g)$ as the sum of the mixture $f_a^m(s^g) = s^g C_0(s^g, \bar{H}_a)u_a^m$ and drift $f_a^d(s^g) = s^g U_d(s^g, \bar{H}_a)u_a^d$ velocity terms. Then, a consistent monotone numerical flux $F_a^g$ is obtained by summing consistent monotone fluxes for each continuous flux functions $f_a^m$ and $f_a^d$. Since $s^g C_0(s^g, \bar{H}_a)$ must be physically a non-decreasing function, the numerical flux for $f_a^m(s^g)$ is just based on an upwinding of the gas saturation w.r.t. the sign of $u^m$ corresponding to the Godunov scheme. On the other hand the function $f_a^d(s^g)$ is never monotone on the full range of gas saturations since it vanishes both for $s^g = 0$ and $s^g = 1$. As a result, the design of a monotone flux for $f_a^d(s^g)$ must be adapted to the specific definition of $U_d$ and will also be dependent on the orientation of the pipe corresponding to the sign of $o_a$. The specific construction described in Appendix A uses the fact that $s^g U_d(s^g, \bar{H}_a)$ is the product of two non-negative functions which are respectively non-decreasing and non-increasing w.r.t. $s^g$.

### 3.4. Monitoring conditions

Assuming the case of a production well, the monitoring conditions prescribe a minimum head node pressure $\bar{p}_\omega$ and a maximum well molar flow rate $\bar{q}_\omega \ge 0$. Setting

$$p_\omega = p_{v_\omega} \quad \text{and} \quad q_\omega = \sum_{\alpha \in \mathcal{P}} |S_{v_\omega}^\omega| \zeta_{v_\omega}^\alpha u_\omega^\alpha, \qquad (3.7)$$

the boundary condition at the head node $v_\omega$ combines the complementary constraints on the pair $(\bar{q}_\omega - q_\omega, p_\omega - \bar{p}_\omega)$ with the slip law:

$$\begin{align}
& q_\omega \le \bar{q}_\omega, \quad p_\omega \ge \bar{p}_\omega, \quad (q_\omega - \bar{q}_\omega)(p_\omega - \bar{p}_\omega) = 0, \nonumber \\
& u_\omega^g = F^g(s_{v_\omega}^g, s_{v_\omega}^g, H_{v_\omega}) = s_{v_\omega}^g C_0(s_{v_\omega}^g, H_{v_\omega}) u_\omega^m + s_{v_\omega}^g U_d(s_{v_\omega}^g, H_{v_\omega}), \tag{3.8} \\
& u_\omega^m = u_\omega^g + u_\omega^l. \nonumber
\end{align}$$

### 4. Nonlinear solver

Given the subset of present phases $Q_v$, let us define the set of well unknowns at node $v \in V_\omega$

$$X_v = (p_v, T_v, s_v^g, s_v^l, c_v^\alpha, \alpha \in Q_v, \tilde{n}_{v,i}, i \in \tilde{C}_{Q_v}).$$

The discrete nonlinear system to be solved at each time step is defined for each well by the set of unknowns

$$(X_v)_{v \in V_\omega}, (Q_v)_{v \in V_\omega}, (u_a^g, u_a^l)_{a \in \mathcal{E}_\omega}, (u_\omega^g, u_\omega^l),$$
and the set of equations (3.2)-(3.6)-(3.8) complemented by the flash equations of type $Q_v = \text{flash}(X_v)$ for each node $v \in V_\omega$. It is coupled fully implicitly to the reservoir system through the source terms $q_{v,i}^{r \to \omega}$, $q_{v,e}^{r \to \omega}$ for $v \in V_\omega$. Let us refer to Appendix B for their expressions depending on the reservoir and well pressure, temperature, molar fractions and saturations at node $v$.

This system is solved using a Newton-Raphson algorithm based on an active set method for the subsets of present phases and for the well active constraints. This type of algorithm can be viewed as an active set formulation of a Newton-Min semi-smooth Newton algorithm (see e.g. [31] and [5]).

The assembly of the Jacobian system is performed at each Newton iteration given the sets of present phases and the well active constraints which are updated after each Newton update. An important feature of the implementation is the Jacobian system reduction by elimination of the local thermodynamical and sum to one of the saturations closure laws (3.2d), of the hydrodynamical equations (3.2c)-(3.6) and of the well monitoring conditions (3.8). The elimination of the local thermodynamical and sum to one of the saturations closure laws is also applied to the reservoir equations. The reduced linear system couples the component and energy conservation equations with a block structure of $\#C + 1$ equations and primary unknowns at each well and reservoir nodes. This elimination procedure detailed in the next subsections facilitates the assembly of the coupled linear system and reduce the cost of its resolution.

## 4.1. Elimination of the thermodynamical and sum to one of the saturations closure laws

Let us denote by $C_v(X_v)$ the thermodynamical and sum to one of the saturations closure equations (3.2d) at node $v$. The elimination of these closure laws is based on a splitting of the unknowns $X_v = (X_v^p, X_v^s)$ into $\#C + 1$ primary unknowns $X_v^p$ and the remaining secondary unknowns $X_v^s$. This splitting classically depends on the set of present phases $Q_v$ and is such that $\frac{\partial C_v(X_v^p, X_v^s)}{\partial X_v^s}$ is non-singular. It results that, after Newton linearization, the secondary unknowns $dX_v^s$ can be expressed as a linear function of $dX_v^p$ and of the residual $C_v(X_v)$ (see [31] for details in the case of the Coats formulation of the reservoir model).

## 4.2. Elimination of the hydrodynamical equations

For all well edge $a = v v'$, given $X_v, X_{v'}, Q_v, Q_{v'}$, the hydrodynamical module computes the phase superficial velocities

$$u_a^\alpha \text{ for } \alpha \in Q_v \cup Q_{v'}$$

and their Newton linearization accounting for the elimination of the thermodynamical and sum to one of the saturations closure laws, such that

$$du_a^\alpha = A_{a,v}^\alpha dX_v^p + A_{a,v'}^\alpha dX_{v'}^p + B_a^\alpha \text{ for } \alpha \in Q_v \cup Q_{v'} \quad (4.1)$$

This computation is detailed in the next two paragraphs. Note that, from the properties of the numerical flux function $F_a^g$, we have $u_a^\alpha = 0$ for $\alpha \notin Q_v \cup Q_{v'}$, hence this case does not need to be considered.

### 4.2.1. Mixture velocity

Let us define the difference of potential at the well edge $a = v v'$ by

$$\Delta_a \Phi = p_v - p_{v'} + \bar{\rho}_a^m g(z_v - z_{v'})$$

Considering the wall friction law (2.2), the momentum equation (3.2c) can be solved for the mixture velocity as a function of $\Delta_a \Phi$ depending on the sign of $\Delta_a \Phi$ as follows

$$u_a^m = -\text{sign}(\Delta_a \Phi) \frac{(\alpha_c - \alpha_b)}{2\alpha_a}, \quad (4.2)$$
with

$$ \alpha_a = |\mathbf{a}| f_q \frac{\bar{\rho}^m}{4r_a}, \quad \alpha_b = 8|\mathbf{a}| \frac{\bar{\mu}^m}{(r_a)^2}, \quad \alpha_c = \sqrt{(\alpha_b)^2 + 4|\Delta_a \Phi|\alpha_a}, $$

and the definition (3.4) of the mean specific density $\bar{\rho}_a^m$ and viscosity $\bar{\mu}_a^m$.

Note that for more general wall friction laws, the elimination of the mixture velocity is not always possible at the nonlinear level. In that case it is done at the linear level after Newton linearization taking into account the momentum equation residual. The mixture velocity $u_a^m$ is computed from (4.2) as a function of the entries

$$ \Delta_a \Phi, \bar{\rho}_a^m, \bar{\mu}_a^m, $$

as well as its derivative w.r.t. the 3 entries.

Then, using the two-point monotone flux function $F_a^g$, we can compute both superficial velocities as detailed below.

### 4.2.2. Superficial velocities

For $Q_v \cup Q_{v'} = \{\alpha\} \in \mathcal{P}, \{\beta\} = \mathcal{P} \setminus \{\alpha\}$, we simply have

$$ u_a^\alpha = u_a^m, \quad u_a^\beta = 0, $$

with $u_a^m$ defined by (4.2).

For $Q_v \cup Q_{v'} = \mathcal{P}$, we first compute $u_a^m$ from (4.2), and then, the gas superficial velocity is obtained from the two-point monotone flux $F_a^g$. Using the model detailed in Appendix A, the entries of $F_a^g$ are

$$ s_{v'}^g, s_v^g, \bar{\rho}_a^\ell, \bar{\rho}_a^g, \bar{\sigma}_{g\ell,a}, u_a^m. $$

The flux function $F_a^g$ computes $u_a^g$ as well as its derivatives w.r.t. the 6 entries. Then the liquid superficial velocity is given by

$$ u_a^\ell = u_a^m - u_a^g. $$

Combining the derivatives of the mixture velocities and superficial velocities w.r.t. to the above entries with the derivatives of the thermodynamical variables, we can compute (4.1).

## 4.3. Elimination of the monitoring conditions

The complementary constraints are solved using an active set method meaning that we impose either a fixed pressure $p_\omega = \bar{p}_\omega$ with $q_\omega \le \bar{q}_\omega$ or a fixed molar flow rate $q_\omega = \bar{q}_\omega$ with $p_\omega \ge \bar{p}_\omega$. The active constraint is updated at each Newton iteration according to the remaining inequality constraint.

The elimination of the head node superficial velocities is performed as follows. In case of active flow rate constraint, the well total molar flow rate is prescribed with

$$ q_\omega = \bar{q}_\omega, \tag{4.3} $$

In case of active pressure constraint, it is obtained by the sum over $i$ of the molar conservation equations at node $v_\omega$

$$ q_\omega = - \sum_{i \in \mathcal{C}} \frac{n_{v_\omega, i} - n_{v_\omega, i}^{n-1}}{\Delta t^n} + \sum_{i \in \mathcal{C}} q_{v_\omega, i}^{r \to \omega} + \sum_{a \in \mathcal{E}_{v_\omega}^\omega} \sum_{\alpha \in \mathcal{P}} \kappa_{a, v_\omega} |S_a^\omega| \zeta_a^\alpha u_a^\alpha. \tag{4.4} $$

Note that this equation must hence be substituted by the equation $p_\omega = \bar{p}_\omega$ in the system of primary equations and unknowns. Using the total molar flow rate $q_\omega$ computed from (4.3) or (4.4), we can eliminate the head node mixture velocity $u_\omega^m$ as follows. Using (3.7) and (3.8), we set

$$ u_\omega^m = \frac{q_\omega}{|S_{v_\omega}^\omega| \zeta_{v_\omega}^\alpha} \quad \text{if } Q_{v_\omega} = \{\alpha\}, $$
and, if $Q_{v_\omega} = \mathcal{P}$, setting $C_{0,v_\omega} = C_0(s_{v_\omega}^g, H_{v_\omega})$, $U_{d,v_\omega} = U_d(s_{v_\omega}^g, H_{v_\omega})$, we obtain the equation

$$ \frac{q_\omega}{|S_{v_\omega}^g|} = \zeta_{v_\omega}^g (s_{v_\omega}^g C_{0,v_\omega} u_\omega^m + s_{v_\omega}^g U_{d,v_\omega}) + \zeta_{v_\omega}^l ((1 - s_{v_\omega}^g C_{0,v_\omega}) u_\omega^m - s_{v_\omega}^g U_{d,v_\omega}). $$

providing

$$ u_{\omega}^{m} = \frac{\frac{q_{\omega}}{|S_{v_{\omega}}^{g}|} + (\zeta_{v_{\omega}}^{l} - \zeta_{v_{\omega}}^{g}) s_{v_{\omega}}^{g} U_{d,v_{\omega}}}{\zeta_{v_{\omega}}^{g} s_{v_{\omega}}^{g} C_{0,v_{\omega}} + \zeta_{v_{\omega}}^{l} (1 - s_{v_{\omega}}^{g} C_{0,v_{\omega}})}. $$

Once $u_\omega^m$ is known we can compute the head node superficial velocities $u_\omega^\alpha$ from (3.8).

# 5. Numerical experiments

The objectives of this numerical Section is to validate the numerical model and investigate its ability to simulate both well opening transient test cases without and with cross flow and fully coupled test cases at the reservoir time scale. The well stand alone test cases are presented in Subsection 5.1 and a fully coupled test case is considered in Subsection 5.2 where it is compared for validation to a simpler single implicit unknown well model.

## 5.1. Stand alone well model

### 5.1.1. Validation on a Buckley Leverett solution

Let us consider a two component $C = \{1, 2\}$, two-phase $\mathcal{P} = \{\ell, g\}$, immiscible, incompressible and isothermal model obtained by setting in the general compositional framework $C^\ell = \{1\}$ and $C^g = \{2\}$. The well is vertical of diameter 0.1 m and centerline $(0, H)$ with $z = 0$ corresponding to the leaf node and $z = H = 100$ m to the head node. The molar and specific densities as well as the dynamic viscosities are considered constant.

Let us set

$$ u^g(s^g) = s^g U_d(s^g) + s^g C_0(s^g) u^m, $$

and

$$ u^{\ell}(s^{g}) = u^{m} - u^{g}(s^{g}). $$

The volume conservation equations for each phase write

$$ \begin{cases} \partial_t s^g + \partial_z u^g(s^g) = 0, \\ \partial_t s^\ell + \partial_z u^\ell(s^g) = 0. \end{cases} $$

Summing the above equations we obtain $\partial_z u_m = 0$ meaning that the mixture velocity depends only on time. It will be fixed by the input boundary condition at the leaf node fixing the mixture velocity to $u^m = 0.5 \text{ m} \cdot \text{s}^{-1}$. Thus, we obtain the scalar hyperbolic Buckely Leverett equation

$$ \partial_t s^g + \partial_z u^g(s^g) = 0, $$

with the flux function $u^g(s^g)$ at fixed mixture velocity. In this test case we compare the solution of the model developed in this work to the one obtained by solving numerically the scalar hyperbolic equation. The same Euler implicit time discretization is used for both implementations. The scalar hyperbolic discretization is based on a cell centered Finite Volume (FV) scheme while the compositional model uses the node centered FV scheme described in this work. The data set is defined by $\rho^g = 4 \text{ kg} \cdot \text{m}^{-3}$, $\rho^\ell = 1000 \text{ kg} \cdot \text{m}^{-3}$, by the input gas superficial velocity at the leaf node $u^g = 0.55 \text{ m} \cdot \text{s}^{-1}$ and the output liquid velocity at the leaf node $u^\ell = -0.05 \text{ m} \cdot \text{s}^{-1}$. The DFM model [24] (see Appendix A) is used with the parameters $A = 1.2$, $B = 0.3$, $a_1 = 0.2$, $a_2 = 0.4$, $K_u = 1.5$ and $\sigma_{g\ell}$ is fixed to $71.97 \times 10^{-3}$ (see the resulting flux function on the left Figure 5.1). The solutions are computed with both models using a uniform mesh of 200 cells, a uniform time stepping with 200 time steps and the
simulation time $t_F = 50$ s. The right Figure 5.1 shows that both models provide basically the same gas saturation plot at final time $t_F$. We also plot in Figure 5.2 the fine mesh gas saturation (Left) and superficial velocities (Right) at final time obtained with the compositional model using 1000 cells and 1000 time steps. The gas moves up at a higher velocity than the mixture velocity as a result of the drift velocity allowing the liquid to go out of the well at the leaf node.

<figure><img src="image_5.png" /><figcaption>FIGURE 5.1. (Left): gas superficial velocity as a function of the gas saturation (flux function) at fixed mixture velocity $u^m = 0.5 \text{ m} \cdot \text{s}^{-1}$ for the Buckley Leverett test case. (Right): Comparison between the Buckley Leverett (BL) and compositional (Comp) gas saturations at final time on the Buckley Leverett test case</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5.2. (Left): Compositional gas saturation at final time using the fine mesh with 1000 cells and 1000 time steps on the Buckley Leverett test case. (Right): liquid and gas superficial velocities obtained with the compositional model at final time using the fine mesh with 1000 cells and 1000 time steps on the Buckley Leverett test case.</figcaption></figure>

### 5.1.2. Thermal test case with a chair shaped production well

We consider a single H₂O component liquid vapor thermal flow along the chair shaped production well illustrated in Figure 5.3. The internal energy, mass density and viscosity of H₂O in the liquid and gas phases are defined by analytical laws as functions of the pressure and temperature (refer [1, Section 4.1] for details). The vapour pressure $P_{sat}(T)$ is given in Pa by

$$p_{sat}(T) = 10^{-3}(T - 273)^4.$$
The DFM model [24] is used with the parameters $A = 1.2$, $B = 0.3$, $a_1 = 0.2$, $a_2 = 0.4$, $K_u = 1.5$ and $\sigma_{g\ell}$ is fixed to $71.97 \times 10^{-3}$ (see Appendix A). The Darcy–Forchheimer friction law (2.2) is set up with the friction parameter $f_q = 6 \times 10^{-2}$. The well radius is fixed to 0.05 m, the thermal conductivity is considered constant for simplicity with $\lambda = 2$.

The well is monitored with the minimum well head pressure set to $\bar{p}_\omega = 5 \times 10^5$ Pa and the maximum well mass flow rate set to $\bar{q}_\omega = 15 \text{ kg} \cdot \text{s}^{-1}$. The initial temperature is set to 350 K along the well and the initial pressure is hydrostatic at liquid state with pressure set to $5 \times 10^5$ Pa at the head node.

Feed zones are modeled at each of the two leaf nodes at high temperature and high pressure in liquid state. At both leaves the reservoir pressure is fixed to $1.1 \times 10^7$ Pa, the temperature to 520 K, the Darcy well index to $WI^D = 10^{-12}$ m and the Fourier well index to $WI^F = 100 \text{ J} \cdot \text{s}^{-1} \cdot \text{K}^{-1}$.

The well is meshed using a uniform discretization of each of the four branches with 40 edges for both bottom vertical branches, 40 edges for the horizontal branch and 60 edges for the top vertical branch. The simulation time is fixed to $t_F = 2000$ s at which the stationary state is basically reached, and the time stepping is set up with an initial time step of 10 s and a maximum time step of 40 s. The stopping criterion of the Newton nonlinear solver is fixed to either $10^{-8}$ on the relative residual $l^1$ norm or to $10^{-10}$ on the $l^\infty$ norm of the Newton increment $ds^g + \frac{dp}{10^5} + \frac{dT}{100}$. The time step is multiplied by the factor 1.1 until it reaches the maximum time step in case of Newton convergence in less than 50 iterations and restarted with a twice smaller time step otherwise. Using this setting, the simulation runs in 58 time steps with no time step failure and a total number of 438 Newton iterations. Figures 5.4 and 5.5 exhibit the final pressure, temperature, gas saturation and superficial velocities. The rise of the hot temperature front induces the appearance of the gas phase at the top of the well starting at roughly $t = 1200$ s (see Figure 5.6). Figure 5.6 plots the time histories of the well pressure $\bar{p}_\omega$, the mass flow rate $\bar{q}_\omega$, the leaf pressures and the gas volume inside the well. The well starts to produce at the fixed minimum well head pressure of $5 \times 10^5$ Pa with a rising mass flow rate. This rising flow rate is induced by the temperature increase along the well which reduces the weight of the liquid column. This can be checked in the leaf pressures plot showing the decrease of the leaf pressures in the first part of the simulation at fixed minimum well pressure. Then, the well flow rate reaches its maximum value of $15 \text{ kg} \cdot \text{s}^{-1}$ and the well operates at fixed maximum mass flow rate with a rising well head pressure which speeds up when the gas phase appears. At around $t = 1400$ s, the well head pressure starts to decrease rapidly as a result of the pressure drop induces by the high gas velocity, until the well operates again at fixed minimum well head pressure of $5 \times 10^5$ Pa. We can also notice, in the leaf pressures plot, the gap between the two leaf pressures which results from the additional wall friction along the horizontal branch.

<figure><img src="image_7.png" /><figcaption>FIGURE 5.3. Chair shaped well with one junction and four branches of sizes 600 m (upper vertical branch), 400 m (lower left vertical branch), 400 m (lower vertical branch) and 400 m (horizontal branch).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 5.4. Pressure (Pa), temperature (K), gas saturation solutions at final time for the chair shaped thermal well test case.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 5.5. Liquid and gas superficial velocities (m·s⁻¹) at final time for the chair shaped thermal well test case.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 5.6. Well mass flow rate, well head pressure, leaf pressures, and gas volume in the well as a function of time for the chair shaped thermal well test case.</figcaption></figure>

5.1.3. *Thermal test case with cross flow*

We consider the same liquid vapor thermal model as in the previous test case along the T shaped
production well illustrated in Figure 5.7. The DFM model [24] is again used with the same parameters
as in the previous test case. In order to obtain cross flow, i.e. a well acting as a producer at some
points along the well and as an injector at others, the reservoir properties at the two leaf nodes are
set up as follows. The bottom left leaf will act as an injection point and is set up with the reservoir
pressure 9 × 10⁶ Pa, the reservoir temperature 520 K, the Darcy well index WID = 10⁻¹² m and the
Fourier well index WIF = 100 J·s⁻¹·K⁻¹. The right leaf will act as a feed zone and is set up with the
reservoir pressure 7 × 10⁶ Pa, the reservoir temperature 500 K, the Darcy well index WID = 10⁻¹² m
and the Fourier well index WIF = 100 J·s⁻¹·K⁻¹.

The well is monitored as in the previous test case with the minimum well head pressure set to $\bar{p}_{\omega} = 5 \times 10^5$ Pa and the maximum well mass flow rate set to $\bar{q}_{\omega} = 15 \text{ kg} \cdot \text{s}^{-1}$. The initial temperature is set to 350 K along the well and the initial pressure is hydrostatic at liquid state with pressure set to $5 \times 10^5$ Pa at the head node.

The well is meshed using a uniform discretization of each of the three branches with 40 edges for the bottom vertical branch, 40 edges for the horizontal branch and 60 edges for the top vertical branch. The simulation time is fixed to $t_F = 2500$ s to reach the stationary state and the time stepping is set up with an initial time step of 0.1 s and a maximum time step of 40 s. Using the same nonlinear solver setting as in the previous test case, the simulation runs in 147 time steps with 5 time step failures and a total number of 1056 Newton iterations.
Figure 5.9 exhibits the final superficial liquid and gas velocities and Figure 5.12 plots the liquid and gas superficial velocities along the vertical (b1) and horizontal (b2) parts of the well at times 1000 s and 2500 s. These figures clearly exhibit the cross flow between the right leaf acting as a feed zone and the bottom leaf acting as an injection point.

Figure 5.8 exhibits the final pressure, temperature, and gas saturation. The high temperature front propagates from the right leaf to both the top and the bottom sides of the vertical part of the well. The rise of the hot temperature front at the low pressure top side of the well induces the appearance of the gas phase starting at roughly *t* = 1800 s (see also Figure 5.10). Figure 5.10 plots the time histories of the well pressure *p̅*<sub>ω</sub>, the mass flow rate *q̅*<sub>ω</sub>, the leaf pressures and the gas volume inside the well. It is shown that the well is monitored at the minimum well head pressure until the gas phase appears reducing the weight of the column which reduces the pressure at the right leaf node and allows to reach the maximum mass flow rate. The well rapidly switches back to the minimum pressure monitoring due to the increase of the wall friction induced by the high gas velocity.

Figure 5.11 exhibits the plots of the pressure, temperature, gas saturation along the vertical (b1) and horizontal (b2) parts of the well at times 1000 s and 2500 s showing the propagation of the temperature front on both sides of the vertical part of the well, as well as the decrease of the pressure along the well as a result of the temperature increase.

<figure><img src="image_11.png" /><figcaption>FIGURE 5.7. T shaped well with one junction and three branches of sizes 600 m (upper vertical branch), 400 m (lower left vertical branch), and 400 m (horizontal branch).</figcaption></figure>
MULTI-SEGMENTED NON-ISOTHERMAL COMPOSITIONAL WELL MODEL

<figure><img src="image_12.png" /><figcaption>FIGURE 5.8. Pressure (Pa), temperature (K), gas saturation solutions at final time for the cross flow thermal well test case.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 5.9. Liquid and gas superficial velocities (m·s⁻¹) at final time for the cross flow thermal well test case.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 5.10. Well mass flow rate, well head pressure, leaf pressures, and gas volume in the well as a function of time for the cross flow thermal well test case.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 5.11. Pressure, temperature and gas saturation at different times along the left vertical left branches (b1) and the right horizontal branch (b2) for the cross flow thermal well test case.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 5.12. Liquid and gas superficial velocities at different times along the left vertical left branches (b1) and the right horizontal branch (b2) for the cross flow thermal well test case.</figcaption></figure>

## 5.2. One production well coupled with a reservoir

The objective of this test case is to validate the coupled reservoir multi-segmented well model by comparison with the results obtained using a simpler well model based on a single implicit unknown. We consider single component H₂O two-phase (liquid and vapor), non-isothermal flow. The DFM model [24] parameters are the same as the ones used in the test case 5.1.2. However, the Darcy-Forchheimer friction law (2.2) is set up with the friction parameter $f_q \in \{0.001, 0.06\}$ where the case $f_q = f_{q1} = 0.001$ corresponds to a very small friction coefficient. This choice will be useful for comparison with previous published results in the literature without friction as explained below. The case $f_q = f_{q2} = 0.06$ corresponds a typical order of magnitude. The geothermal reservoir is defined by the domain $\Omega = (-H, H)^2 \times (0, H_z)$ where $H = 1000$ m and $H_z = 200$ m, and we consider 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 = k\mathbb{I}, k = 5 \times 10^{-14} \text{m}^2$ and porosity $\phi = 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 defined by analytical laws as functions of the pressure and temperature (refer [1, Section 4.1] for details). The vapour pressure $P_{sat}(T)$ is given in Pa by the Clausius–Clapeyron equation

$$p_{sat}(T) = 100 \exp \left( 46.784 - \frac{6435}{T} - 3.868 \log(T) \right).$$

The reservoir thermal conductivity is fixed to $\lambda^r = 2 \text{ W} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$, and the rock volumetric heat capacity is given by $C_s = 1.6 \text{ MJ} \cdot \text{K}^{-1} \cdot \text{m}^{-3}$ with $E^s(p, T) = C_s T$. The relative permeabilities are set to $k_r^\alpha(s^\alpha) = (s^\alpha)^2$ for both phases $\alpha \in \{\ell, g\}$. The gravity vector is as usual $g = (0, 0, -g_z)$ with $g_z = 9.81 \text{ m} \cdot \text{s}^{-2}$. The simulation consists in two stages. At the first one, the well is closed and we impose a Dirichlet boundary condition at the top of the domain prescribing the reservoir pressure and the temperature equal to $p^r = 4 \text{ MPa}$ and $T^r = (p_{sat})^{-1}(p^r) - 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^r$ 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. The first stage is run until the simulation reaches a stationary state with the liquid phase only, a constant temperature and a 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$ bar (never reached in practice) and the maximum total mass flow rate $\bar{q}_{\omega} = 200 \text{ ton} \cdot \text{h}^{-1}$. The second stage is run on the time interval $(0, t_F)$ with $t_F = 30$ days. Moreover, at the top of the reservoir and at the root of the well, there is no coupling between them during both stages. The entire simulation runs on uniform Cartesian mesh of size $n_x \times n_y \times n_z$ to discretize the domain $\Omega$ with $(n_x, n_y, n_z) = (20, 20, 10)$. The well indexes $WI^D$ and $WI^F$ introduced in Appendix B are computed at each node of the well following [4].

In order to validate the coupled model, we compare the numerical results against the ones published in [1, Section 4.1] using the same thermodynamic parameters, the same initialization and stages, and the same mesh. The model presented in [4] is based on the following assumptions on the well model:

(i) the wall friction, the thermal conduction and the transient terms are all neglected all along the well,

(ii) the thermal conduction is neglected between the reservoir and the well,

(iii) there is no cross flow and both the gas and liquid velocities are oriented in the same direction all along the well,

(iv) the pressure drop along the well is computed based on an explicit approximation of the mean density $\rho^m$.

In that case, the well model can be shown (see e.g. [4]) to reduce to a single implicit unknown, the well head node pressure $p_{\omega}$, and to a single equation, the monitoring conditions, fully coupled to the reservoir system. The computation of the mean density $\rho^m$ along the well is based on a gas liquid flash computation providing the well temperature and gas saturation at each well node. This flash computation is based on the lagged in time values of the well pressures and molar and energy flow rates $q_{V,i}^{r \to \omega}$, $q_{V,e}^{r \to \omega}$. It usually assumes a zero slip law $u^g = s^g u^m$ as detailed in [4], but it can also easily account for a general slip law (2.1) by a simple modification of the gas liquid flash computation. In order to investigate further the effects of the DFM slip law, we consider in this test case the model proposed in [4] with both a zero slip law and its enhancement using the same slip law as the one used in the multi-segmented well model. For the sake of completeness, the modified flash computations implied by the non-zero slip law is detailed in Appendix C.

To make easier the comparison of the results obtained by each model, we label by MSwell-fq₁ and MSwell-fq₂ the results obtained by the proposed Multi-Segmented well model (3.2) with $f_{q1} = 0.001$ and $f_{q2} = 0.06$, respectively. We label by SIUwell the ones obtained by the Single Implicit Unknown well model proposed in [4], and we label by SIUwell-DFM its enhancement including the DFM slip law (2.1) in the gas liquid flash computation. Figure 5.13 compares the pressure, the temperature and the gas saturation along the well; respectively, at final time $t_F$ using this three different models. Figure 5.14 show the total volume of gas inside the well and the reservoir as functions of time for the three models. It can be noticed from those figures that the results obtained with the MSwell model (3.2) with $f_{q1} = 0.001$ are very close to the ones obtained using the well model of [4] provided that the DFM slip law is taken into account. This is expected for such configuration with no cross flow, no significant thermal conduction losses between the well and the reservoir, no significant wall friction, and no significant transient effects at the reservoir time scale. On the other hand, the more realistic value $f_{q2}$ of the friction coefficient induces a larger pressure drop which results in significant variations between the MSwell-fq₂ and the SIUwell-DFM models.
Table 5.1 shows the numerical efficiency of the nonlinear solver for all models using the same mesh for the second stage of the simulation. We denote by $N_{\Delta t}$ the number of time steps and by $N_{\text{Newton}}$ the average number of Newton iterations per time step. It exhibits as expected the higher number of Newton iterations for the multi-segmented well models compared with the single implicit unknown well models, as a result of stronger nonlinearities. We note also that the nonlinear convergence of the multi-segmented well model is sensitive to very small (not physical) friction coefficients due to an increased stiffness of the mixture velocity as a function of $\Delta_a \Phi$ in such cases (see (4.2)).

<figure><img src="image_17.png" /><figcaption>FIGURE 5.13. Pressure in Pa, temperature in °C and gas saturation along the well at final time using the three different models.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 5.14. Total gas volume inside the well (left) and inside the reservoir (right) as a function of time using the three different models.</figcaption></figure>
<table><caption>TABLE 5.1. Numerical behavior of the second stage of the simulation for different models using the same mesh. $N_{\Delta t}$ is the number of time steps, and $N_{\text{Newton}}$ the average number of Newton iterations per time step.</caption><thead><tr><th>Model</th><th>#M</th><th>N<sub>Δt</sub></th><th>N<sub>Newton</sub></th></tr></thead><tbody><tr><td>MSwell-fq<sub>1</sub></td><td>32000</td><td>121</td><td>10.12</td></tr><tr><td>MSwell-fq<sub>2</sub></td><td>32000</td><td>120</td><td>6.15</td></tr><tr><td>SIUwell</td><td>32000</td><td>133</td><td>1.56</td></tr><tr><td>SIUwell-DFM</td><td>32000</td><td>133</td><td>1.58</td></tr></tbody></table>

## 6. Conclusion

A numerical model for two-phase compositional non-isothermal flow in geothermal multi-branch wells is developed to simulate the flow and transport along the wells in geothermal field operations. The model combines the flexible thermodynamical Coats formulation with the drift flux hydrodynamical model leading to a system coupling the conservation equations for each component, momentum and energy.

The discretization is fully implicit in time and based on a staggered finite volume scheme in space, with node centred control volumes for components and energy conservations and edge control volumes for the momentum conservation. The numerical fluxes combine a monotone flux approximation for the phase superficial velocities on the full range of gas saturation with an upwind approximation of the phase molar fractions, density and enthalpy. The interplay between these two parts of the fluxes is a key ingredient of the stability of the scheme.

The nonlinear solver benefits from the elimination of the well superficial velocities and flow rates in the Newton linearization process leading to a Jacobian system with the same primary unknowns both for the well and the reservoir, which makes the coupling between both models much easier.

This numerically robust well model can be used to simulate complex transient flows such as those occurring during the start-up of a geothermal well, as well as cross-flow configurations along the wells. It can be used either by modeling the reservoir with source terms using well indexes, or fully coupled with a subsurface flow model. Different configurations (transient, stationary, coupled or not with a reservoir flow model) have been tested to validate the model performance.

Future work will focus on the implementation of alternative drift flux models in particular to take downward flows into account, and on coupled well-reservoir modeling for large scale industrial case studies.

## Appendix A. Example of monotone two-point flux for the DFM model from [24]

The gas liquid DFM model introduced in [24] is based on the following choices of the drift velocity and profile parameter.

Let $\sigma_{g\ell}$ be the gas liquid interfacial tension, and let us define the following characteristic velocity for the rise of a gas bubble in a liquid:

$$U_c = \left( \sigma_{g\ell} g \frac{(\rho^\ell - \rho^g)}{(\rho^\ell)^2} \right)^{\frac{1}{4}}.$$
The profile parameter is defined by

$$ egin{cases} C_0(s^g) = \frac{A}{1 + (A - 1)\gamma^2}, \\ \gamma = P_{[0,1]}\left(\frac{\beta - B}{1 - B}\right), \\ \beta = \max\left(s^g, F_\nu s^g \frac{|u^m|}{V_{sgf}}\right), \\ V_{sgf} = K_u \sqrt{\frac{\rho^\ell}{\rho^g}} U_c, \end{cases} \quad (A.1) $$

with constant parameters A, B which must be such that $B < (2 - A)/A$ to ensure that $s^gC_0(s^g) \le 1$ and $s^gC_0(s^g)$ non-decreasing. Typically we set $A = 1.2$, $B = 0.3$. The Critical Kutateladze number $K_u \in [0, 3.5]$ depends on the dimensionless diameter of the pipe defined by

$$ \hat{D} = \left( g \frac{(\rho^\ell - \rho^g)}{\sigma_{g\ell}} \right)^{\frac{1}{2}} D, $$

where D is the diameter of the pipe. A graphical representation of the function $K_u$ is provided in [24]. The constant parameter $F_\nu$ is typically set to 1 and the projection $P_{[0,1]}$ is defined by

$$ P_{[0,1]}(x) = \begin{cases} 0 & \text{if } x \le 0, \\ 1 & \text{if } x \ge 1, \\ x & \text{if } x \in (0, 1). \end{cases} $$

The drift velocity times $s^g$ is defined by

$$ s^g U_d(s^g) = G(s^g) \tilde{K}(s^g) U_c, \quad (A.2) $$

with

$$ G(s^g) = \frac{(1 - s^g C_0(s^g))}{s^g C_0(s^g) \sqrt{\frac{\rho^g}{\rho^\ell}} + 1 - s^g C_0(s^g)}, $$

and

$$ \tilde{K}(s^g) = \begin{cases} 1.53 \, s^g & \text{if } s^g \le a_1, \\ K_u C_0(s^g) s^g & \text{if } s^g \ge a_2, \\ 1.53 \, a_1 \frac{s^g - a_2}{a_1 - a_2} + K_u C_0(a_2) \frac{s^g - a_1}{a_2 - a_1} & \text{if } s^g \in (a_1, a_2). \end{cases} $$

Note that the linear interpolation between $a_1$ and $a_2$ is done in [24] on $\frac{\tilde{K}(s^g)}{s^g C_0(s^g)}$ rather than on $\tilde{K}(s^g)$. It has been modified here to simplify the design of a monotone flux. The constant parameters $0 < a_1 < a_2 < 1$ are typically set to $a_1 = 0.2$ and $a_2 = 0.4$. We can check that $G(s^g)$ is a non-increasing function w.r.t. $s^g$. We also assume in the following that $\tilde{K}(s^g)$ is a non-decreasing function w.r.t. $s^g$ which is the case provided that the condition

$$ 1.53 \, a_1 \le a_2 K_u C_0(a_2), $$

is satisfied. In practice, it seems that this condition is not very restrictive. This condition guarantees that the following flux function $F_a^g$ is a monotone two-point flux:

$$ \left\{ \begin{aligned} F_a^g(u, v) = & u C_0(u) (u_a^m)^+ + v C_0(v) (u_a^m)^- \\ & + G(v) \tilde{K}(u) (U_c o_a)^+ + G(u) \tilde{K}(v) (U_c o_a)^-. \end{aligned} \right. \quad (A.3) $$
In the case $1.53 a_1 > a_2 K_u C_0(a_2)$, a Godunov numerical flux $\tilde{K}_g(u, v)$ could still be computed for the continuous flux function $\tilde{K}(s^g)$ since $\tilde{K}(s^g)$ is in that case a piecewise monotone function non decreasing on $(0, a_1) \cup (a_2, 1)$ and non-increasing on $(a_1, a_2)$. Then, the second line in (A.3) will be replaced by $G(v) \tilde{K}_g(u, v)(U_c o_a)^+ + G(u) \tilde{K}_g(v, u)(U_c o_a)^-$.

## Appendix B. Reservoir source terms

It is assumed that the radius of the wells are 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 $v \in V_\omega$ is obtained using the Two Point Flux Approximation

$$V_v^\alpha = W_v^D(p_v^{r,\alpha} - p_v),$$

where $p_v^{r,\alpha}$ is the reservoir phase pressure at node $v$. Fourier fluxes between the reservoir and the well are discretized in the same way using the Two Point Flux Approximation

$$F_v = W_v^F (T_v^r - T_v),$$

with $T_v^r$ denoting the reservoir temperature at node $v$. The Well Indexes $W_v^D$ and $W_v^F$ are typically computed using Peaceman's approach (see [4, 7, 20, 21]) and take into account the unresolved singularity of respectively the pressure and temperature solutions in the neighborhood of the well. Let us denote by $k_{r,v}^\alpha(s^\alpha)$ the phase relative permeability at node $v$ as a function of the phase saturation $s^\alpha$, by $s_v^{r,\alpha}$ the reservoir saturation of phase $\alpha$, and by $c_v^{r,\alpha}$ the reservoir molar fractions of phase $\alpha$.

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

$$q_{v,\alpha,i}^{r \to \omega} = c_{v,i}^{\alpha} \frac{\zeta_v^{\alpha}}{\mu_v^{\alpha}} k_{r,v}^{\alpha}(s_v^{\alpha})(V_v^{\alpha})^{-} + c_{v,i}^{\alpha} \frac{\zeta^{\alpha}(p_v^r, T_v^r, c_v^{r,\alpha})}{\mu^{\alpha}(p_v^r, T_v^r, c_v^{r,\alpha})} k_{r,v}^{\alpha}(s_v^{r,\alpha})(V_v^{\alpha})^{+}, \quad (B.1)$$

$$q_{v,i}^{r \to \omega} = \sum_{\alpha \in \mathcal{P}_i} q_{v,\alpha,i}^{r \to \omega},$$

and the energy flow rate is defined similarly by

$$q_{v,e}^{r \to \omega} = \sum_{\alpha \in \mathcal{P}} \left( h_{v}^{\alpha} (q_{v,\alpha}^{r \to \omega})^{-} + h^{\alpha} (p_{v}^{r}, T_{v}, c_{v}^{r,\alpha}) (q_{v,\alpha}^{r \to \omega})^{+} \right) + F_{v}, \quad (B.2)$$

$$\text{with } q_{v,\alpha}^{r \to \omega} = \sum_{i \in \mathcal{C}^{\alpha}} q_{v,\alpha,i}^{r \to \omega}.$$

## Appendix C. Computations of the well temperatures and saturations for the single implicit unknown well model

We detail in this Appendix the computations of the production well temperatures and saturations for the single implicit unknown well model of Subsection 5.2 at given well pressures $p_v$, molar $q_{v,h_{20}}^{r \to \omega}$ and energy $q_{v,e}^{r \to \omega}$ flow rates for $v \in V_\omega$. Using the assumptions (i) and (ii) stated in Subsection 5.2 for the SIUwell-DFM well model, and considering the case of a single component H$_2$O, the conservation equations (3.2a)-(3.2b) reduce to

$$\sum_{a \in \mathcal{E}_v^\omega} \sum_{\alpha \in \mathcal{P}} -\kappa_{a,v} |S_a^\omega| \zeta_a^\alpha u_a^\alpha = q_{v,h_{20}}^{r \to \omega} - \delta_v^{v,\omega} \sum_{\alpha \in \mathcal{P}} |S_{v,\omega}^\omega| \zeta_v^\alpha u_\omega^\alpha,$$

$$\sum_{a \in \mathcal{E}_v^\omega} \sum_{\alpha \in \mathcal{P}} -\kappa_{a,v} |S_a^\omega| h_a^\alpha \zeta_a^\alpha u_a^\alpha = q_{v,e}^{r \to \omega} - \delta_v^{v,\omega} \sum_{\alpha \in \mathcal{P}} |S_{v,\omega}^\omega| h_v^\alpha \zeta_v^\alpha u_\omega^\alpha,$$
for $\mathbf{v} \in \mathcal{V}_\omega$. Summing these molar and energy conservation equations over all nodes $\mathbf{v}'' \ge \mathbf{v}$, and using assumption (iii) for the enthalpy and molar density upwind values, we obtain for all $\mathbf{a} = \mathbf{v}'\mathbf{v} \in \mathcal{E}_\omega$ that

$$ \sum_{\alpha \in \mathcal{P}} |S_{\mathbf{a}}^{\omega}| \zeta^{\alpha}(p_{\mathbf{v}}, T_{\mathbf{v}}) u_{\mathbf{a}}^{\alpha} = \sum_{\mathbf{v}'' \in \mathcal{V}_{\omega} | \mathbf{v}'' \ge \mathbf{v}^{\omega}} q_{\mathbf{v}'' , \mathbf{h}_{20}}^{r \to \omega} := Q_{\mathbf{v} , \mathbf{h}_{20}}^{\omega}, $$

$$ \sum_{\alpha \in \mathcal{P}} |S_{\mathbf{a}}^{\omega}| h^{\alpha}(p_{\mathbf{v}}, T_{\mathbf{v}}) \zeta^{\alpha}(p_{\mathbf{v}}, T_{\mathbf{v}}) u_{\mathbf{a}}^{\alpha} = \sum_{\mathbf{v}'' \in \mathcal{V}_{\omega} | \mathbf{v}'' \ge \mathbf{v}^{\omega}} q_{\mathbf{v}'' , e}^{r \to \omega} := Q_{\mathbf{v}, e}^{\omega}. $$

It results that the thermodynamical equilibrium at fixed well pressure $p_{\mathbf{v}}$, molar $Q_{\mathbf{v},h20}^{\omega}$ and energy $Q_{\mathbf{v},e}^{\omega}$ provides the well temperature $T_{\mathbf{v}}$ and the well saturations $s_{\mathbf{v}}^{\alpha}$ at node $\mathbf{v}$ as follows. Let us define the phase molar fractions $c_{\mathbf{v}}^{\alpha}$, $\alpha \in \mathcal{P}$ such that

$$ c_{\mathbf{v}}^{\alpha} Q_{\mathbf{v}, h_{20}}^{\omega} = |S_{\mathbf{a}}^{\omega}| \zeta^{\alpha}(p_{\mathbf{v}}, T_{\mathbf{v}}) u_{\mathbf{a}}^{\alpha}. $$

We first assume that both phases are present which implies that $T_{\text{sat}} = (p_{\text{sat}})^{-1}(p_{\mathbf{v}})$ and that the liquid molar fraction is given by

$$ c_{\mathbf{v}}^{\ell} = \frac{h^{\mathrm{g}}(p_{\mathbf{v}}, T_{\mathrm{sat}}) - \frac{Q_{\mathbf{v},e}^{\omega}}{Q_{\mathbf{v},h_{20}}^{\omega}}}{h^{\mathrm{g}}(p_{\mathbf{v}}, T_{\mathrm{sat}}) - h^{\ell}(p_{\mathbf{v}}, T_{\mathrm{sat}})}, $$

from which we can compute $c_{\mathbf{v}}^{\mathrm{g}} = 1 - c_{\mathbf{v}}^{\ell}$ and the superficial velocities $u_{\mathbf{a}}^{\ell}$ and $u_{\mathbf{a}}^{\mathrm{g}}$. The following alternatives are checked:

**Two-phase state:** if $0 < c_{\mathbf{v}}^{\ell} < 1$, the two-phase state is confirmed. Using the slip law $u^{\mathrm{g}} = s^{\mathrm{g}}U_{\mathrm{d}}(s^{\mathrm{g}}) + s^{\mathrm{g}}C^{0}(s^{\mathrm{g}})(u^{\mathrm{g}} + u^{\ell})$, then $T_{\mathbf{v}} = T_{\text{sat}}$ and the gas saturation $s_{\mathbf{v}}^{\mathrm{g}}$ is solution of the equation

$$ s_{\mathbf{v}}^{\mathrm{g}} U_{\mathrm{d}}(s_{\mathbf{v}}^{\mathrm{g}}) + s_{\mathbf{v}}^{\mathrm{g}} C^{0}(s_{\mathbf{v}}^{\mathrm{g}})(u_{\mathbf{a}}^{\mathrm{g}} + u_{\mathbf{a}}^{\ell}) - u_{\mathbf{a}}^{\mathrm{g}} = 0. $$

In the no slip case, corresponding to $U_{\mathrm{d}} = 0$ and $C^{0} = 1$, it reduces to $s_{\mathbf{v}}^{\mathrm{g}} = \frac{u_{\mathbf{a}}^{\mathrm{g}}}{u_{\mathbf{a}}^{\mathrm{g}} + u_{\mathbf{a}}^{\ell}}$.

**Liquid state:** if $c_{\mathbf{v}}^{\ell} \ge 1$, then only the liquid phase is present, we set $s_{\mathbf{v}}^{\ell} = 1$, $s_{\mathbf{v}}^{\mathrm{g}} = 0$, and $T_{\mathbf{v}}$ is the solution of

$$ h^{\ell}(p_{\mathbf{v}}, T_{\mathbf{v}}) = \frac{Q_{\mathbf{v},e}^{\omega}}{Q_{\mathbf{v},h_{20}}^{\omega}}. $$

**Gas state:** if $c_{\mathbf{v}}^{\ell} \le 0$, then only the gas phase is present, we set $s_{\mathbf{v}}^{\ell} = 0$, $s_{\mathbf{v}}^{\mathrm{g}} = 1$, and $T_{\mathbf{v}}$ is the solution of

$$ h^{\mathrm{g}}(p_{\mathbf{v}}, T_{\mathbf{v}}) = \frac{Q_{\mathbf{v},e}^{\omega}}{Q_{\mathbf{v},h_{20}}^{\omega}}. $$

We note that the same computations are done at the head node $\mathbf{v}_{\omega}$ based on the equations

$$ \sum_{\alpha \in \mathcal{P}} |S_{\mathbf{v}_{\omega}}^{\omega}| \zeta^{\alpha}(p_{\mathbf{v}_{\omega}}, T_{\mathbf{v}_{\omega}}) u_{\omega}^{\alpha} = \sum_{\mathbf{v}'' \in \mathcal{V}_{\omega}} q_{\mathbf{v}'' , \mathbf{h}_{20}}^{r \to \omega} := Q_{\mathbf{v}_{\omega}, \mathbf{h}_{20}}^{\omega}, $$

$$ \sum_{\alpha \in \mathcal{P}} |S_{\mathbf{v}_{\omega}}^{\omega}| h^{\alpha}(p_{\mathbf{v}_{\omega}}, T_{\mathbf{v}_{\omega}}) \zeta^{\alpha}(p_{\mathbf{v}_{\omega}}, T_{\mathbf{v}_{\omega}}) u_{\omega}^{\alpha} = \sum_{\mathbf{v}'' \in \mathcal{V}_{\omega}} q_{\mathbf{v}'' , e}^{r \to \omega} := Q_{\mathbf{v}_{\omega}, e}^{\omega}. $$

Then, the pressures $p_{\mathbf{v}}$, temperatures $T_{\mathbf{v}}$ and saturations $s_{\mathbf{v}}^{\alpha}$, $\mathbf{v} \in \mathcal{V}_{\omega}$ are used to compute the edge mean densities $\bar{\rho}_{\mathbf{a}}^{\mathrm{m}}$ defined by (3.4). These mean densities are frozen in the momentum equation for the computation of the next time step. Together with the zero wall friction assumption (i) (which could be relaxed by an explicit approximation), it results that the next time step well pressures depend only on the well head node pressure $p_{\omega}$. This is the only well unknown implicitly coupled to the reservoir system combined with the single well equation reducing to the complementary constraints on the pair $(\bar{q}_{\omega} - q_{\omega}, p_{\omega} - \bar{p}_{\omega})$ (first equation of (3.8)). Let us refer to see [4] for more details about this type of single implicit unknown well model.
References

[1] A. Armandine Les Landes, D. Castanon Quiroz, L. Jeannin, S. Lopez, and R. Masson. Two-phase geothermal model with fracture network and multi-branch wells. SMAI J. Comput. Math., 9:121-149, 2023.

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

[3] A. Battistelli, S. Finsterle, M. Marcolini, and L. Pan. Modeling of coupled wellbore-reservoir flow in steam-like supercritical geothermal systems. Geothermics, 86: article no. 101793, 2020.

[4] L. Beaude, T. Beltzung, K. Brenner, S. Lopez, R. Masson, F. Smai, J.-F. Thebault, and F. Xing. Parallel Geothermal Numerical Model with Fractures and Multi-Branch Wells. ESAIM: ProcS, 63:109-134, 2018.

[5] L. Beaude, K. Brenner, S. Lopez, R. Masson, and F. 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.

[6] G. Bjornsson. A multi-feedzone geothermal wellbore simulator. Technical report, 1987.

[7] Z. Chen and Y. Zhang. Well flow models for various numerical methods. Int. J. Numer. Anal. Model., 6:375-388, 2009.

[8] T. J. Danielson. Transient Multiphase Flow: Past, Present, and Future with Flow Assurance Perspective. Energy Fuels, 26(7):4137-4144, 2012.

[9] J.-M. Delhaye. Thermohydraulique des réacteurs nucléaires. Génie atomique. EDP Sciences, 2008.

[10] T. Gallouët and A. Pfertzel. On some upstream weighting schemes for oil recovery simulation. Les Annales de l'Enit, 4(2), 1990.

[11] H. Gudmundsdottir and M. T. Jonsson. The Wellbore Simulator FloWell - Model Enhancement and Verification. In World Geothermal Congress 2015, 2015.

[12] François P. Hamon, Bradley T. Mallison, and Hamdi A. Tchelepi. Implicit Hybrid Upwind scheme for coupled multiphase flow and transport with buoyancy. Comput. Methods Appl. Mech. Eng., 311:599-624, 2016.

[13] M. Jerez-Carrizales, J. E. Jaramillo, and D. Fuentes. Prediction of Multiphase Flow in Pipelines: Literature Review. Ingenieria y Ciencia, 11(22):213-233, 2015.

[14] Khasani, Deendarlianto, and R. Itoi. Numerical study of the effects of CO₂ gas in geothermal water on the fluid-flow characteristics in production wells. Eng. Appl. Comput. Fluid Mech., 15(1):111-129, 2021.

[15] B. Lamy-Chappuis, A. Yapparova, and T. Driesner. Advanced well model for superhot and saline geothermal reservoirs. Geothermics, 105: article no. 102529, 2022.

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

[17] S. Lopez, V. Hamm, M. Le Brun, L. Schaper, F. Boissier, C. Cotiche, and E. Giuglaris. 40 years of Dogger aquifer management in Ile-de-France, Paris Basin, France. Geothermics, 39(4):339-356, 2010.

[18] J. M. Masella, Q. H. Tran, D. Ferre, and C. Pauchon. Transient simulation of two-phase flows in pipes. Int. J. Multiphase Flow, 24(5):739-755, 1998.

[19] L. Pan and C. M. Oldenburg. T2Well-An integrated wellbore-reservoir simulator. Comput. Geosci., 65:46-55, 2014.

[20] D. W. Peaceman. Interpretation of Well-Block Pressures in Numerical. SPE J., 18(03):183-194, 1978.

[21] D. W. Peaceman. Interpretation of Well-Block Pressures in Numerical Reservoir Simulation with Nonsquare Grid Blocks and Anisotropic Permeability. SPE J., 23(03):531-543, 1983.
[22] V. E. Pereira and A. C. Fowler. Exsolving two-phase flow in oil wells. *Geophys. Astrophys. Fluid Dyn.*, 114(3):283-305, 2020.

[23] A. Batman, Z. Ugur, and M. Onur. The effect of calcite deposition on geothermal well in flow performance. *Geothermics*, 28(3):425-444, 1999.

[24] H. Shi, J. A. Holmes, L. J. Durlofsky, K. Aziz, L. Diaz, B. Alkaya, and G. Oddie. Drift-Flux Modeling of Two-Phase Flow in Wellbores. *SPE J.*, 10(01):24–33, 2005.

[25] M. Shippen and W. J. Bailey. Steady-state Multiphase Flow: Past, Present, and Future with Flow Assurance Perspective. *Energy Fuels*, 26(7):4145–4157, 2012.

[26] R. A. Tonkin, J. P. O'Sullivan, M. Gravatt, and M. J. O'Sullivan. Simulation of complex transient flows in geothermal wells. *Geothermics*, 108: article no. 102631, 2023.

[27] R. A. Tonkin, J. P. O'Sullivan, M. Gravatt, and M. J. O'Sullivan. A transient geothermal wellbore simulator. *Geothermics*, 110: article no. 102653, 2023.

[28] R. A. Tonkin, M. J. O'Sullivan, and J. P. O'Sullivan. A review of mathematical models for geothermal wellbore simulation. *Geothermics*, 97: article no. 102255, 2021.

[29] Q. H. Tran, I. Faille, C. Pauchon, and F. Willien. The No Pressure Wave (NPW) Model: Application to Oil and Gas Transport. In *Third International Symposium on Finite Volume for Complex Applications Porquerolles France*, 2002.

[30] G. B. Wallis. *One-dimensional two-phase flow*. McGraw-Hill, 1969.

[31] 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.

[32] M. Zolfagharroshan and E. Khamechi. A rigorous approach to scale formation and deposition modelling in geothermal wellbores. *Geothermics*, 87: article no. 101841, 2020.

[33] N. Zuber and J. A. Findlay. Average Volumetric Concentration in Two-Phase Flow Systems. *J. Heat Transfer*, 87(4):453–468, 1965.

[34] N. Zuber and J. A. Findlay. Average volumetric concentration in two-phase flow systems. *J. Heat Transfer*, 87:453–458, 1965.