D Castanon Quiroz

L Jeannin

S Lopez
email: s.lopez@brgm.fr

R Masson

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

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

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 an 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. It 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.

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 the place of two-phase flows and phase changes. During apraisal 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 of the fluids produced at the surface. A thorough understanding of wells and their modeling allow a better understanding of the reservoir behaviour.

Physical phenomena and exchanges with the reservoir are complex. For example, the startup 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 produce 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 behaviour 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 occuring 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, [START_REF] Shippen | Steady-state Multiphase Flow: Past, Present, and Future with Flow Assurance Perspective[END_REF] for steady flows, [START_REF] Danielson | Transient Multiphase Flow: Past, Present, and Future with Flow Assurance Perspective[END_REF] for transient flows or [START_REF] Jerez-Carrizales | Prediction of Multiphase Flow in Pipelines: Literature Review[END_REF]). The nuclear industry has also a strong interest in multiphase thermal flows in pipes ( [START_REF] Delhaye | Thermohydraulique des réacteurs nucléaires[END_REF] 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 [START_REF] Tonkin | A review of mathematical models for geothermal wellbore simulation[END_REF], 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 [START_REF] Bjornsson | A multi-feedzone geothermal wellbore simulator[END_REF] 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, [START_REF] Gudmundsdottir | The Wellbore Simulator FloWell -Model Enhancement and Verification[END_REF] for matching field data or [START_REF] Battistelli | Modeling of coupled wellborereservoir flow in steam-like supercritical geothermal systems[END_REF] for the modeling of the IDDP-1 supercritical well).

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 [START_REF] Zuber | Average Volumetric Concentration in Two-Phase Flow Systems[END_REF][START_REF] Wallis | One-dimensional two-phase flow[END_REF]. 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 [START_REF] Wallis | One-dimensional two-phase flow[END_REF][START_REF] Delhaye | Thermohydraulique des réacteurs nucléaires[END_REF][START_REF] Tonkin | A review of mathematical models for geothermal wellbore simulation[END_REF]. 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 [START_REF] Tonkin | A review of mathematical models for geothermal wellbore simulation[END_REF][START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF].

It is usually assumed that the thermodynamical equilibrium assumption applies in the well, but some authors are interested in non-equilibrium phenomena [START_REF] Pereira | Exsolving two-phase flow in oil wells[END_REF][START_REF] Aunzo | Wellbore Models GWELL, GWNACL, and HOLA, user's guide[END_REF] to describe scaling or kinetics phenomena. [START_REF] Satman | The effect of calcite deposition on geothermal well in flow performance[END_REF][START_REF] Lopez | 40 years of Dogger aquifer management in Ile-de-France[END_REF][START_REF] Khasani | Numerical study of the effects of CO 2 gas in geothermal water on the fluid-flow characteristics in production wells[END_REF][START_REF] Zolfagharroshan | A rigorous approach to scale formation and deposition modelling in geothermal wellbores[END_REF]. 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 [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF]. 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 [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF], 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 [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] for the slip law. Moreover, following [START_REF] Tran | The No Pressure Wave (NPW) Model: Application to Oil and Gas Transport[END_REF][START_REF] Livescu | A fully-coupled thermal multiphase wellbore flow model for use in reservoir simulation[END_REF][START_REF] Masella | Transient simulation of two-phase flows in pipes[END_REF], 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 assumption remains valid to describe pressure transients (see [START_REF] Tonkin | A review of mathematical models for geothermal wellbore simulation[END_REF] for a numerical validation).

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 [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF] 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 velocites is obtained on the full range of saturations [0, 1] using an hybrid upwinding approach exploiting the mathematical properties of the drift velocity and profile parameter of the DFM model proposed in [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF]. Note that the methodology developped in this work could be applied to other slip laws sharing the same types of mathematical properties.

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 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 [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF] 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, density and enthalpy. The detailed specific construction of the monotone superficial velocity flux function based on an hybrid upwinding approach for the DFM model proposed in [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] is reported to Appendix 7.1. The well monitoring conditions are introduced in Subsection 3.4 in the case of a production well and Appendix 7.2 recall 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.

Well physical model

Let W denote the set of wells. As in [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF], each multi-branch well ω ∈ 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 ( ) 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 α , and molar fractions c α = (c α i ) i∈C , α ∈ P = { , 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 α ⊂ C the set of components in phase α and by P i the set of phases containing the component i ∈ C.

The thermodynamical model uses a Coats' formulation already developped for the reservoir model (see e.g. [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF]). It is based on the natural variables p, T , s α , c α , α ∈ Q with a set of unknowns and equations depending on the additional unknown Q ⊂ 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 α, we denote by ζ α (p, T, c α ) its molar density, by ρ α (p, T, c α ) its specific density, by µ α (p, T, c α ) its dynamic viscosity, by e α (p, T, c α ) its molar internal energy, and by h α (p, T, c α ) its molar enthalpy. The fugacities are denoted by f α i (p, T, c α ).

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 [START_REF] Zuber | Average volumetric concentration in two-phase flow systems[END_REF]. 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 [START_REF] Tran | The No Pressure Wave (NPW) Model: Application to Oil and Gas Transport[END_REF][START_REF] Livescu | A fully-coupled thermal multiphase wellbore flow model for use in reservoir simulation[END_REF].

Let τ denote the curvilinear coordinate along the well and e τ the unit tangential vector along the well assumed to be oriented toward the well root node. We denote by u α = u α e τ the superficial velocity of each phase α ∈ P. The mixture velocity is defined by

u m = u + u g = u m e τ .
The tangential divergence along the well is denoted by div τ and the tangential gradient by ∇ τ . The section along the well is denoted by S ω (τ ). We also set g τ = (g • e τ )e τ where g is the acceleration of gravity vector.

Let ω (τ ) = (∇τ z(τ ))•eτ ) |∇τ z(τ )|
∈ {-1, 1} define the orientation along the well, with arbitrary value for ∇ τ z = 0 along an horizontal branch. Following [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF], let us introduce the function

o ω (τ ) = ω (τ ) cos(θ ω (τ )) 1/2 (1 + sin(θ ω (τ ))) 2 ,
where θ ω (τ ) ∈ [0, π/2] is the edge angle w.r.t. the vertical direction. Note that o ω (τ ) = 0 along an horizontal branch. The slip closure law introduced in [START_REF] Zuber | Average volumetric concentration in two-phase flow systems[END_REF] 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 , H)o ω + s g C 0 (s g , H)u m , (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, H) = 1, U d (s g = 1, H) = 0, s g C 0 (s g , H) ≤ 1 and s g C 0 (s g , H) non decreasing w.r.t. s g . The drift velocity and profile parameter depend on a set of thermodynamical variables denoted by 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 ω .

We denote by

T f = T f e τ ,
the wall friction law depending on the mixture velocity u m and on thermodynamical quantities.

To fix ideas, we will use in the following a Darcy-Forchheimer type law:

T f = - 8µ m r 2 ω + f q ρ m 1 4r ω |u m | u m , (2)
where f q is the friction coefficient, µ m = s g µ g + s µ the mean viscosity and ρ m = s g ρ g + s ρ the mean specific density.

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

α∈P i ∂ t (S ω c α i ζ α s α ) + div τ (S ω c α i ζ α u α ) = q r→ω i , i ∈ C, ( 3a
) α∈P ∂ t (S ω e α ζ α s α ) + div τ (S ω h α ζ α u α ) + div τ (-S ω λ∇ τ T ) = q r→ω e , (3b)
∇ τ p = T f + ρ m g τ , (3c)
f g i (p, T, c g ) = f l i (p, T, c l ) if Q ∩ P i = P, i∈C α c α i = 1, α ∈ Q, α∈Q s α = 1, s α = 0 if α ∈ Q, (3d)
u g = s g U d (s g , H)o ω + s g C 0 (s g , H)u m if Q = P, (3e)
where q r→ω i and q r→ω e 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.

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 7.1 based on the model proposed in [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF]. Thanks to the monotonicity and consistency properties of the superficial velocities, it is 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.

Well discretization and notations

The set of nodes of a well ω ∈ W is denoted by V ω and its root node is denoted by v ω . A partial ordering is defined on the set of vertices V ω with v < ω v if and only if the unique path from the root v ω to v passes through v. The set of edges of the well ω is denoted by E ω and for each edge a ∈ E ω we set a = vv with v < ω v (i.e. v is the parent node of v , see Figure 1). It is

assumed that V ω 1 ∩ V ω 2 = ∅ for any ω 1 , ω 2 ∈ W such that ω 1 = ω 2 .
Let |S ω a | and |S ω vω | denote the well section at respectively the center of the edge a and at the head node v ω . For a ∈ E ω , |a| denotes the lengh of the edge a. Let us define E ω v ⊂ E ω as the set of well edges sharing the node

v ∈ V ω . For all vv = a ∈ E ω , let us set κ a,v = -1 and κ a,v = 1. well a = vv u α a p v , T v , s α v , c α v u α ω v v ω v q r →ω v,i q r →ω v,e
Figure 1: Example of multi-branch well ω with its root node v ω , one edge a = vv (v parent node of v ) and the main physical quantities: the head node superficial velocities u α ω (non negative for production wells and non positive for injection wells), the molar and energy flow rates between the reservoir and the well q r→ω v,i , q r→ω v,e , the well node pressure, temperature, saturations and phase molar fractions p v , T v , s α v , c α v , and the phase superficial velocity u α a at the edge a oriented positively from v to v.

For each edge a = vv ∈ E ω , let us denote by u α a the superficial velocity of phase α 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 ∈ V ω by

• the set of present phases Q v ,
• the well pressure p v ,

• the well temperature T v ,

• the well saturations s α v for α ∈ P, • the well molar fractions c α v for α ∈ Q v , • the number of moles ñv,i for i ∈ C Qv , and by

• the edge superficial velocities u α a for α ∈ P at each well edge a ∈ E ω ,

• the head node superficial velocities u α ω for α ∈ P at the well head node v ω . In the above definition, C Q denotes the set of components not contained in α∈Q C α . For C Qv = ∅, the additional unknowns ñv,i corresponding to the number of moles of the absent components i ∈ C Qv 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 [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF] for details). Note that the saturation s α v of an absent phase α ∈ Q v vanishes. In the following, to fix ideas, we assume that the monitoring conditions at the well head node correspond to the case of a production well consistently with the numerical section.

For a given thermodynamical phase property ξ α at a node v, we will use the notation

ξ α v = ξ α (p v , T v , c α v ), with typically ξ α = ζ α , ρ α , e α , h α , f α i , µ α .
The molar and energy flow rates between the reservoir and the well ω at a given node v ∈ V ω are defined by a two point flux approximation between the reservoir and well properties at node v based on the Peaceman approach combined with a phase potential upwinding. They are denoted respectively by q r→ω v,i for each component i ∈ C and q r→ω v,e in what follows (see Figure 1 and Appendix 7.2 for their definitions).

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 =        ( a∈E ω v |S ω a | |a| 2 ) α∈Qv∩P i c α v,i s α v ζ α v if i ∈ C Qv , ñv,i else, n v,e = ( a∈E ω v |S ω a | |a| 2 ) α∈Qv s α v ζ α v e α v ,
where the case of an absent component i ∈ C Qv is accounted for by the introduction of the additional unknown ñv,i in order to avoid the singularity of the discrete system.

For each edge a = vv ∈ E ω , and each phase α, 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:

ξ α a = ξ α v if u α a ≥ 0, ξ α v if u α a < 0, (4)
with

ξ α = c α , ζ α , h α .
We denote by (0, t F ) the time interval and consider its discretization given by t n , n = 0, • • • , N with t 0 = 0, t N = t F and ∆t n = t nt n-1 > 0 for all n = 1, • • • , N . The time discretization of the system (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 (3a)-(3b)-(3c)-(3d) at each time step n = 1, • • • , N is defined as follows:

n v,i -n n-1 v,i ∆t n + a∈E ω v α∈P i -κ a,v |S ω a |c α a,i ζ α a u α a = q r→ω v,i -δ vω v α∈P i |S ω vω |c α v,i ζ α v u α ω , i ∈ C, v ∈ V ω , (5a)
n v,e -n n-1 v,e ∆t n + a∈E ω v |S ω a | λ a |a| (T v -T v ) + α∈P -κ a,v h α a ζ α a u α a = q r→ω v,e -δ vω v α∈P |S ω vω |h α v ζ α v u α ω , v ∈ V ω , (5b)
p v -p v + ρm a g(z v -z v ) = T f a |a|, vv = a ∈ E ω , (5c)
f g i,v = f i,v if Q v ∩ P i = P, v ∈ V ω , i∈C α c α v,i = 1, α ∈ Q v , v ∈ V ω , α∈Qv s α v = 1, s α v = 0 if α ∈ Q v , v ∈ V ω . ( 5d
)
where δ stands for the Kronecker symbol with 5), the edge thermal conductivity of the mixture is defined by

δ vω v = 1 if v = v ω , else 0. In (
λ a = α∈P s α v + s α v 2 λ α ,
with λ α the thermal conductivity of the phase α assumed constant for simplicity. The wall friction term T f a is given by the Darcy-Forchheimer law (6) leading to

T f a = - 8μ m a r 2 a + f q ρm a 1 4r a |u m a | u m a , (6)
with the following arithmetic means of the mixture specific density and dynamic viscosity at the edge a

ρm a = 1 2 α∈Qv s α v ρ α v + α∈Q v s α v ρ α v , μm a = 1 2 α∈Qv s α v µ α v + α∈Q v s α v µ α v . (7)
It remains to define the discretization of the edge superficial velocities u α a accounting for the slip law (3e) at given mixture velocity u m = u g + u . This framework is detailed in the next subsection based on a monotone two-point flux.

Monotone two-point flux for the edge superficial velocities

Let us define the arithmetic average of a thermodynamical variable ξ α for a given phase α at a well edge a = vv with α

∈ Q v ∪ Q v by ξα a =              ξ α (p v , T v , c α v ) if α ∈ Q v and α ∈ Q v , ξ α (p v , T v , c α v ) else if α ∈ Q v and α ∈ Q v , s α v ξ α (p v , T v , c α v ) + s α v ξ α (p v , T v , c α v ) s α v + s α v else if s α v + s α v > , 1 2 ξ α (p v , T v , c α v ) + ξ α (p v , T v , c α v ) else.
This arithmetic averaging is used to compute the thermodynamical variables Ha at the well edge a = vv which enter in the definition of the slip law (3e). Then, the superficial velocities u α a are obtained at given mixture velocity u m a using a numerical two-point monotone flux denoted by

F g a (s g v , s g v , Ha ),
for the continuous flux function

f a (s g ) = s g U d (s g , Ha )o a + s g C 0 (s g , Ha )u m a ,
where o a is the value of the function o ω along the edge a. The numerical flux function F g a must satisfy the following properties:

• consistency property:

F g a (s g , s g , Ha ) = f a (s g ) for all s g ∈ [0, 1],
• monotonicity property which specifies that, for all (

s g v , s g v ) ∈ [0, 1]×[0, 1], F g a (s g v , s g v , Ha
) is non decreasing w.r.t. its first argument s g v and non increasing w.r.t. its second argument s g v .

Then the superficial velocities u α a are defined by

u g a = F g a (s g v , s g v , Ha ), u a = u m a -u g a . (8)
From the consistency and monotonicity properties, it results that

• s g v = 0 ⇒ u g a = F g (0, s g v ) ≤ 0 for all s g v ∈ [0, 1], • s g v = 0 ⇒ u g a = F g (s g v , 0) ≥ 0 for all s g v ∈ [0, 1], • s g v = 1 ⇒ u a = u m a -F g (1, s g v ) ≤ 0 for all s g v ∈ [0, 1], • s g v = 1 ⇒ u a = u m a -F g (s g v , 1
) ≥ 0 for all s g v ∈ [0, 1]. These properties are key for the definition (4) 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 ξ α a u α a is always properly defined, with

ξ α = ζ α , c α , h α . Note also that if Q v ∪ Q v = {α} (
single phase α present at both nodes), the properties of the numerical flux function imply that u α a = u m a and u β a = 0 for the absent phase β. It means that, in that cases, the thermodynamical variables Ha 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 g a depends on the laws C 0 and U d . We propose in Appendix 7.1 a monotone two-point flux for the gas liquid DFM model introduced in [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] on the full range of gas saturation.

Monitoring conditions

Assuming the case of a production well to fix ideas, the monitoring conditions prescribe a minimum head node pressure pω and a maximum well molar flow rate qω ≥ 0. Setting

p ω = p vω and q ω = α∈P |S ω vω |ζ α vω u α ω , (9)
the boundary condition at the head node v ω combines the complementary constraints on the pair qω -q ω , p ω -pω with the slip law:

q ω ≤ qω , p ω ≥ pω , (q ω -qω )(p ω -pω ) = 0, u g ω = F g (s g vω , s g vω , H vω ) = s g vω C 0 (s g vω , H vω )u m ω + s g vω U d (s g vω , H vω ), u m ω = u g ω + u ω . ( 10
)

Nonlinear solver

Given the subset of present phases Q v , let us define the set of well unknowns at node v ∈ V ω

X v = p v , T v , s g v , s v , c α v , α ∈ Q v , ñv,i , i ∈ C Qv .
The discrete nonlinear system to be solved at each time step is defined for each well by the set of unknowns (X v ) v∈Vω , (Q v ) v∈Vω , (u g a , u a ) a∈Eω , (u g ω , u ω ), and the set of equations ( 5)-( 8)-( 10) complemented by the flash equations of type Q v = flash(X v ) for each node v ∈ V ω . It is coupled fully implicitly to the reservoir system through the source terms q r→ω v,i , q r→ω v,e for v ∈ V ω . Let us refer to Appendix 7.2 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. [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF] and [START_REF] Beaude | Non-isothermal compositional liquid gas darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations[END_REF]).

The Jacobian system is assembled 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 reduction of the Jacobian system by elimination of the local thermodynamical closure laws (5d), of the hydrodynamical equations (5c)-( 8) and of the well monitoring conditions [START_REF] Gudmundsdottir | The Wellbore Simulator FloWell -Model Enhancement and Verification[END_REF]. The elimination of the thermodynamical 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.

Elimination of the thermodynamical closure laws

Let us denote by C v (X v ) the thermodynamical closure equations (5d) at node v. The elimination of the thermodynamical closure laws is based on a splitting of the unknowns X v = (X p v , X s v ) into #C + 1 primary unknowns X p v and the remaining secondary unknowns X s v . This splitting classically depends on the set of present phases Q v and is such that

∂C $ (X p v ,X s v ) ∂X s v is non singular.
It results that, after Newton linearization, the secondary unknowns dX s v can be expressed as a linear function of dX p v and of the residual C v (X v ) (see [START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal darcy flows in fractured porous media[END_REF] for details in the case of the Coats' formulation of the reservoir model).

Elimination of the hydrodynamical equations

For all well edge a = vv , given X v , X v , Q v , Q v , the hydrodynamical module computes the phase superficial velocities

u α a for α ∈ Q v ∪ Q v ,
and their Newton linearization accounting for the elimination of the thermodynamical closure laws, such that

du α a = A α a,v dX p v + A α a,v dX p v + B α a for α ∈ Q v ∪ Q v . ( 11
)
This computation is detailed in the next two paragraphs. Note that, from the properties of the numerical flux function F g a , we have u α a = 0 for α ∈ Q v ∪ Q v , hence this case does not need to be considered.

Mixture velocity

Let us define the difference of potential at the well edge a = vv by

∆ a Φ = p v -p v + ρm a g(z v -z v ).
Considering the wall friction law (6), the momentum equation (5c) can be solved for the mixture velocity as a function of ∆ a Φ depending on the sign of ∆ a Φ as follows

u m a = -sign(∆ a Φ) (α c -α b ) 2α a , (12)
with

α a = |a|f q ρm 4r a , α b = 8|a| μm (r a ) 2 , α c = (α b ) 2 + 4|∆ a Φ|α a ,
and the definition [START_REF] Chen | Well flow models for various numerical methods[END_REF] of the mean specific density ρm a and viscosity μm a . Note that for more general wall friction laws, the elimination of the mixture velocity is not always possible at the non linear 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 m a is computed from [START_REF] Khasani | Numerical study of the effects of CO 2 gas in geothermal water on the fluid-flow characteristics in production wells[END_REF] as a function of the entries ∆ a Φ, ρm a , μm a , as well as its derivative w.r.t. the 3 entries. Then, using the two-point monotone flux function F g a , we can compute both superficial velocities as detailed below.

Superficial velocities

For Q v ∪ Q v = {α} ∈ P, {β} = P \ {α}, we simply have u α a = u m a , u β a = 0,
with u m a defined by [START_REF] Khasani | Numerical study of the effects of CO 2 gas in geothermal water on the fluid-flow characteristics in production wells[END_REF]. [START_REF] Khasani | Numerical study of the effects of CO 2 gas in geothermal water on the fluid-flow characteristics in production wells[END_REF], and then, the gas superficial velocity is obtained from the two-point monotone flux F g a . Using the model detailed in Appendix 7.1, the entries of F g a are s g v , s g v , ρ a , ρg a , σg ,a , u m a . The flux function F g a computes u g a as well as its derivatives w.r.t. the 6 entries. Then the liquid superficial velocity is given by u a = u m a -u g a . 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 [START_REF] Jerez-Carrizales | Prediction of Multiphase Flow in Pipelines: Literature Review[END_REF].

For Q v ∪ Q v = P, we first compute u m a from

Elimination of the monitoring conditions

The complementary constraints are solved using an active set method meaning that we impose either a fixed pressure p ω = pω with q ω ≤ qω or a fixed molar flow rate q ω = qω with p ω ≥ pω . 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 ω = qω , (13)
In case of active pressure constraint, it is obtained by the sum over i of the molar conservation equations at node v ω

q ω = - i∈C n vω,i -n n-1 vω,i ∆t n + i∈C q r→ω vω,i + a∈E ω vω α∈P κ a,vω |S ω a |ζ α a u α a . (14)
Note that this equation must hence be substituted by the equation p ω = pω in the system of primary equations and unknowns. Using the total molar flow rate q ω computed from ( 13) or ( 14), we can eliminate the head node mixture velocity u m ω as follows. Using ( 9) and ( 10), we set

u m ω = q ω |S ω vω |ζ α vω if Q vω = {α},
and, if

Q vω = P, setting C 0,vω = C 0 (s g vω , H vω ), U d,vω = U d (s g vω , H vω ), we obtain the equation q ω |S ω vω | = ζ g vω s g vω C 0,vω u m ω + s g vω U d,vω + ζ vω (1 -s g vω C 0,vω )u m ω -s g vω U d,vω .
providing

u m ω = qω |S ω vω | + (ζ vω -ζ g vω )s g vω U d,vω ζ g vω s g vω C 0,vω + ζ vω (1 -s g vω C 0,vω )
.

Once u m ω is known we can compute the head node superficial velocities u α ω from (10).

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. Let us set u g (s g ) = s g U d (s g ) + s g C 0 (s g )u m , and u (s g ) = u mu g (s g ).

Stand alone well model

The volume conservation equations for each phase write

∂ t s g + ∂ z u g (s g ) = 0, ∂ t s + ∂ z u (s g ) = 0.
Summing the above equations we obtain ∂ 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 m.s -1 . Thus, we obtain the scalar hyperbolic Buckely Leverett equation

∂ t s g + ∂ 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 developped 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 ρ g = 4 Kg.m -3 , ρ l = 1000 Kg.m -3 , by the input gas superficial velocity at the leaf node u g = 0.55 m.s -1 and the output liquid velocity at the leaf node u = -0.05 m.s -1 . The DFM model [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] (see Appendix 7.1) is used with the parameters A = 1.2, B = 0.3, a 1 = 0.2, a 2 = 0.4, K u = 1.5 and σ g is fixed to 71.97 10 -3 (see the resulting flux function on the left figure 2). 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 2 shows that both models provide basically the same gas saturation plot at final time t F . We also plot in Figure 3 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. The DFM model [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] is used with the parameters A = 1.2, B = 0.3, a 1 = 0.2, a 2 = 0.4, K u = 1.5 and σ g is fixed to 71.97 10 -3 (see Appendix 7.1). The Darcy-Forccheimer friction law ( 6) is set up with the friction parameter f q = 6 10 -2 . The well radius is fixed to 0.05 m, the thermal conductivity is considered constant for simplicity with λ = 2.

The well is monitored with the minimum well head pressure set to pω = 5 10 5 Pa and the maximum well mass flow rate set to qω = 15 Kg.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 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 10 7 Pa, the temperature to 520 K, the Darcy well index to W I D = 10 -12 m and the Fourier well index to W I F = 100 J.s -1 .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 stationnary 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 criteria 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 ∞ norm of the Newton increment ds g + dp 10 5 + 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 and6 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 7). Figure 7 plots the time histories of the well pressure pω , the mass flow rate qω , the leaf pressures and the gas volume inside the well. The well starts to produce at the fixed minimum well head pressure of 5 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 Kg.s -1 and the well operates at fixed maximun 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 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.

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 8. The DFM model [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] is again used with the same parameters as in the previous test case. In order to obtain cross flow, 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 6 Pa, the reservoir temperature 520 K, the Darcy well index W I D = 10 -12 m the Fourier well index W I F = 100 J.s -1 .K -1 . The right leaf will act as a feed zone and is set up with the reservoir pressure 7 10 6 Pa, the reservoir temperature 500 K, the Darcy well index W I D = 10 -12 m and the Fourier well index W I F = 100 J.s -1 .K -1 .

The well is monitored as in the previous test case with the minimum well head pressure set to pω = 5 10 5 Pa and the maximum well mass flow rate set to qω = 15 Kg.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 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 stationnary 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 10 exhibits the final superficial liquid and gas velocities and Figure 13 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 9 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 11). Figure 11 plots the time histories of the well pressure pω , the mass flow rate qω , 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 12 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.

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 2 O two-phase (liquid and vapor), non-isothermal flow. The DFM model [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] parameters are the same as the ones used in the test case 5.1.2. However, the Darcy-Forccheimer friction law ( 6) is set up with the friction parameter f q ∈ {0.001, 0.06} where the case f q = q 1 = 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 q 2 = 0.06 corresponds a typical order of magnitude. The geothermal reservoir is defined by the domain Ω = (-H, H) The reservoir thermal conductivity is fixed to λ r = 2 W.m -1 .K -1 , and the rock volumetric heat capacity is given by C s = 1.6 MJ.K -1 .m -3 with E s (p, T ) = C s T . The relative permeabilities are set to k α r (s α ) = (s α ) 2 for both phases α ∈ { , g}. The gravity vector is as usual g = (0, 0, -g z ) with g z = 9.81 m.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 MPa and T r = (p sat ) -1 (p r ) -1 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 an hydrostatic pressure depending only on the vertical coordinate.

For the second stage, homogeneous Neumann boundary conditions are prescribed at the bottom and at the top of the domain Ω, 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 pω = 1 bar (never reached in practice) and the maximum total mass flow rate qω = 200 ton.hour -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 × n y × n z to discretize the domain Ω with (n x , n y , n z ) = [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF][START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF][START_REF] Gudmundsdottir | The Wellbore Simulator FloWell -Model Enhancement and Verification[END_REF]. The well indexes W I D and W I F introduced in Appendix 7.2 are computed at each node of the well following [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF].

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 [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF] 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 ρ m .

In that case, the well model can be shown (see e.g. [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF]) to reduce to a single implicit unknown, the well head node pressure p ω , and to a single equation, the monitoring conditions, fully coupled to the reservoir system. The computation of the mean density ρ 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 r→ω v,i , q r→ω v,e . It usually assumes a zero slip law u g = s g u m as detailed in [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF], but it can also easily account for a general slip law (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 [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF] 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 7.3.

To make easier the comparison of the results obtained by each model, we label by MSwell-fq 1 and MSwell-fq 2 the results obtained by the proposed Multi-Segmented well model ( 5) with f q 1 = 0.001 and f q 2 = 0.06, respectively. We label by SIUwell the ones obtained by the Single Implicit Unknown well model proposed in [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF], and we label by SIUwell-DFM its enhancement including the DFM slip law (1) in the gas liquid flash computation. Figure 14 compares the pressure, the temperature and the gas saturation along the well; respectively, at final time t F using this three different models. Figure 15 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 ( 5) with f q 1 = 0.001 are very close to the ones obtained using the well model of [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF] 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 q 2 of the friction coefficient induces a larger pressure drop which results in significant variations between the MSwell-fq 2 and the SIUwell-DFM models.

Table 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 ∆t the number of time steps and by N 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 ∆ a Φ in such cases (see ( 12)).

2.8x10 6 3.2x10 6 3.6x10 6 Table 1: Numerical behavior of the second stage of the simulation for different models using the same mesh. N ∆t is the number of time steps, and N Newton the average number of Newton iterations per time step.

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 a 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.

Appendices

7.1 Appendix I: example of monotone two-point flux for the DFM model from [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] The gas liquid DFM model introduced in [START_REF] Shi | Drift-flux modeling of two-phase flow in wellbores[END_REF] is based on the following choices of the drift velocity and profile parameter.

Let σ g 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 = σ g g (ρ -ρ g ) (ρ ) 2

Appendix II: 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 ∈ V ω is obtained using the Two Point Flux Approximation

V α v = WI D v (p r,α v -p v )
, where p r,α v 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 = WI F v (T r v -T v
), with T r v denoting the reservoir temperature at node v. The Well Indexes WI D v and WI F v are typically computed using Peaceman's approach (see [START_REF] Peaceman | Interpretation of Well-Block Pressures in Numerical[END_REF][START_REF] Peaceman | Interpretation of Well-Block Pressures in Numerical Reservoir Simulation with Nonsquare Grid Blocks and Anisotropic Permeability[END_REF][START_REF] Chen | Well flow models for various numerical methods[END_REF][START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF]) 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 (s α ) the phase relative permeability at node v as a function of the phase saturation s α , by s r,α v the reservoir saturation of phase α, and by c r,α v the reservoir molar fractions of phase α.

For any a ∈ R, let us define a + = max(a, 0) and a -= min(a, 0). The molar flow rates between the reservoir and the well ω at a given node v ∈ V ω are defined by the following phase based upwind approximation of the mobilities:

q r→ω v,α,i = c α v,i ζ α v µ α v k α r,v (s α v )(V α v ) -+ c α v,i ζ α (p r v , T r v , c r,α v ) µ α (p r v , T r v , c r,α v ) k α r,v (s r,α v )(V α v ) + , q r→ω v,i = α∈P i q r→ω v,α,i , (18)
and the energy flow rate is defined similarly by q r→ω v,e = α∈P h α v (q r→ω v,α ) -+ h α (p r v , T r v , c r,α v )(q r→ω v,α )

+ + F v , (19)
with q r→ω v,α = i∈C α q r→ω v,α,i .

Appendix III: 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 r→ω v,h 2 o and energy q r→ω v,e flow rates for v ∈ V ω . 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 (5a)-(5b) reduce to

a∈E ω v α∈P -κ a,v |S ω a |ζ α a u α a = q r→ω v,h 2 o -δ vω v α∈P |S ω vω |ζ α v u α ω , a∈E ω v α∈P -κ a,v |S ω a |h α a ζ α a u α a = q r→ω v,e -δ vω v α∈P |S ω vω |h α v ζ α v u α ω ,
for v ∈ V ω . Summing these molar and energy conservation equations over all nodes v ≥ ω v, and using assumption (iii) for the enthalpy and molar density upwind values, we obtain for all

a = v v ∈ E ω that α∈P |S ω a |ζ α (p v , T v )u α a = v ∈Vω|v ≥ ω v q r→ω v ,h 2 o := Q ω v,h 2 o , α∈P |S ω a |h α (p v , T v )ζ α (p v , T v )u α a = v ∈Vω|v ≥ ω v q r→ω v ,e := Q ω v,e .
It results that the thermodynamical equilibrium at fixed well pressure p v , molar Q ω v,h 2 o and energy Q ω v,e provides the well temperature T v and the well saturations s α v at node v as follows. Let us define the phase molar fractions c α v , α ∈ P such that

c α v Q ω v,h 2 o = |S ω a |ζ α (p v , T v )u α a .
We first assume that both phases are present which implies that T sat = (p sat ) -1 (p v ) and that the liquid molar fraction is given by

c v = h g (p v , T sat ) - Q ω v,e Q ω v,h 2 o h g (p v , T sat ) -h (p v , T sat ) ,
from which we can compute c g v = 1 -c v and the superficial velocities u a and u g a . The following alternatives are checked: Two-phase state: if 0 < c v < 1, the two-phase state is confirmed. Using the slip law u g = s g U d (s g ) + s g C 0 (s g )(u g + u ), then T v = T sat and the gas saturation s g v is solution of the equation s g v U d (s g v ) + s g v C 0 (s g v )(u g a + u a ) -u g a = 0.

In the no slip case, corresponding to U d = 0 and C 0 = 1, it reduces to s g v = u g a u g a +u a

.

Liquid state: if c v ≥ 1, then only the liquid phase is present, we set s v = 1, s g v = 0, and T v is the solution of

h (p v , T v ) = Q ω v,e Q ω v,h 2 o
.

Gas state: if c v ≤ 0, then only the gas phase is present, we set s v = 0, s g v = 1, and T v is the solution of

h g (p v , T v ) = Q ω v,e Q ω v,h 2 o
.

We note that the same computations are done at the head node v ω based on the equations Then, the pressures p v , temperatures T v and saturations s α v , v ∈ V ω are used to compute the edge mean densities ρm a defined by [START_REF] Chen | Well flow models for various numerical methods[END_REF]. 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 ω . 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 qω -q ω , p ω -pω (first equation of [START_REF] Gudmundsdottir | The Wellbore Simulator FloWell -Model Enhancement and Verification[END_REF]). Let us refer to see [START_REF] Beaude | Parallel geothermal numerical model with fractures and multi-branch wells[END_REF] for more details about this type of single implicit unknown well model.

  1 and 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.

5. 1 . 1

 11 Validation on a Buckley Leverett solutionLet us consider a two component C = {1, 2}, two-phase P = { , g}, immiscible, incompressible and isothermal model obtained by setting in the general compositional framework C = {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.

Figure 2 :Figure 3 :

 23 Figure 2: (Left): gas superficial velocity as a function of the gas saturation (flux function) at fixed mixture velocity u m = 0.5 m.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

Figure 4 :

 4 Figure 4: 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).

Figure 5 :

 5 Figure 5: Pressure (Pa), temperature (K), gas saturation solutions at final time for the chair shaped thermal well test case.

Figure 6 :

 6 Figure 6: Liquid and gas superficial velocities (m.s -1 ) at final time for the chair shaped thermal well test case.

Figure 7 :

 7 Figure 7: 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.

Figure 8 :

 8 Figure 8: 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).

Figure 9 :

 9 Figure 9: Pressure (Pa), temperature (K), gas saturation solutions at final time for the cross flow thermal well test case.

Figure 10 :

 10 Figure 10: Liquid and gas superficial velocities (m.s -1 ) at final time for the cross flow thermal well test case.

Figure 11 :

 11 Figure 11: 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.

Figure 12 :Figure 13 :

 1213 Figure 12: 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.

2 ×

 2 (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) ∈ Ω | x = y = 0} of radius r ω = 0.1 m. The reservoir is assumed homogeneous with isotropic permeability K = kI, k = 5×10 -14 m 2 and porosity φ = 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 46.784 -6435 T -3.868 log(T ) .

Figure 14 :Figure 15 :

 1415 Figure 14: Pressure in Pa, temperature in • C and gas saturation along the well at final time using the three different models.

  α∈P |S ω vω |ζ α (p vω , T vω )u α ω = v ∈Vω q r→ω v ,h 2 o := Q ω vω,h 2 o , α∈P |S ω vω |h α (p vω , T vω )ζ α (p vω , T vω )u α ω = v ∈Vω q r→ω v ,e := Q ω vω,e .

The profile parameter is defined by

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

where D is the diameter of the pipe. The constant parameter F ν is typically set to 1 and the projection P [0,1] is defined by

The drift velocity times s g is defined by

with G(s g ) = (1 -s g C 0 (s g ))

Note that the linear interpolation between a 1 and a 2 is done in [20] on K(s g ) s g C 0 (s g ) rather than on 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 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 ≤ 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 g a is a monotone two-point flux:

If 1.53 a 1 > a 2 K u C 0 (a 2 ), a Godunov scheme could still be computed for K(s g ) since it is in that case a piecewise monotone function non decreasing on (0, a 1 ) ∪ (a 2 , 1) and non increasing on (a 1 , a 2 ).