Laurence Beaude
email: laurence.beaude@unice.fr

Konstantin Brenner
email: konstantin.brenner@unice.fr

Simon Lopez
email: s.lopez@brgm.fr

Roland Masson
email: roland.masson@unice.fr

Farid F Smaï

Farid Smai
email: f.smai@brgm.fr

Non-isothermal compositional liquid gas Darcy flow: formulation, soil-atmosphere boundary condition and application to high energy geothermal simulations

Keywords: Non-isothermal compositional Darcy flow, geothermal energy, soil-atmosphere boundary condition, outflow boundary condition, porous medium drying, finite volume scheme

come

the world current final energy consumption [START_REF] Limberger | Geothermal energy in deep aquifers: A global assessment of the resource base for direct heat utilization[END_REF].

The quantitative understanding of the shallow parts of geothermal systems is challenging both for the exploration and exploitation of high energy geothermal resources. First, the unsaturated zone and/or cooler superficial water flows can considerably alter evidences of the presence of a deeper geothermal resource. In some cases, the resource may be totally hidden. In terms of exploitation, as some systems underlay urbanized areas (e.g. Rotuara in New Zealand or Bouillante in the French Caribbean), resource exploitation much be carefully monitored and controlled in order to avoid unwanted induced surface manifestations or risks. Moreover, several features such as geysers, have a major cultural significance for indigenous populations and must therefore be protected and kept unaltered [START_REF] O'sullivan | Using Reservoir Models for Monitoring Geothermal Surface Features[END_REF].

Numerical modelling has become essential in all phases of geothermal operations. It is used in the exploration phases to assess the geothermal potential, validate conceptual hypothesis and help well siting. Field development and resource management need quantitative estimation to prevent resource exhaustion and achieve its sustainable exploitation (production/injection scenarios). Finally, numerical modelling is also helpful in studying exploitation related to industrial risks such as the interaction with shallow water levels (drinking water resources, hydrothermal vents or eruption) (e.g. [START_REF] O'sullivan | Using Reservoir Models for Monitoring Geothermal Surface Features[END_REF]). The physics embedded in the numerical model should correctly handle non-linear evolution of saturation transients with water table fluctuations, high temperature gradients and phase change processes in the shallow levels of the simulation domain.

Current software suffers several limitations in terms of boundary conditions which are known to play a major role in geothermal flows [START_REF] O'sullivan | Geothermal reservoir simulation: the state-of-practice and emerging trends[END_REF][START_REF] O'sullivan | Future directions in geothermal modelling[END_REF]. Mixed-type transient boundary conditions are not supported which impedes the convenient modelling of natural processes such as recharge or seepage or water table fluctuations. Workarounds may exist (e.g. [START_REF] Hurwitz | Groundwater flow, heat transport, and water table position within volcanic edifices: Implications for volcanic processes in the Cascade Range[END_REF]) but are relatively tedious to implement and are not formulated in a generic way. Transient complex upper or lower boundary conditions are mandatory to take into account some crucial processes. In volcanic island settings, the inland water table may be excessively deep and the interactions between the vadose zone and the fresh water recharge may hide geothermal resources ( [START_REF] Byrdina | Influence of the regional topography on the remote emplacement of hydrothermal systems with examples of Ticsani and Ubinas volcanoes. Southern Peru[END_REF], [START_REF] Hurwitz | Groundwater flow, heat transport, and water table position within volcanic edifices: Implications for volcanic processes in the Cascade Range[END_REF]). In sedimentary basins the interactions with the topography and recharge areas must be correctly taken into account to reproduce head distributions at basin scale [START_REF] Dentzer | Réévaluation du potentiel géothermique en Ile de France[END_REF]. Though many groundwater simulation software programs can deal with the vadose zone, they are rarely designed to study multiphasic hydrothermal processes. Conversely, some geothermal reservoir simulators propose to take into account a gas/air component [START_REF] O'sullivan | Improved Convergence for Air-Water and CO2-Water TOUGH2 Simulations[END_REF] but they are still restricted to rather simple boundary conditions and most of the time the alternative is between fixed value/Dirichlet type for all primary variables or fixed fluxes/Neumann type for all conserved quantities (e.g. [START_REF] Pruess | TOUGH2 user's guide, version 2[END_REF]).

The objective of this work is to investigate a new formulation for non-isothermal compositional gas liquid Darcy flows and to couple it with an advanced soil-atmosphere boundary condition. The compositional model should typically account for the water component which can vaporize into the gas phase and for a set of gaseous components which can dissolve into the liquid phase. The soil-atmosphere boundary condition, based on mole and energy balance equations set at the interface, should take into account the vaporization of the liquid phase in the atmosphere, the convective molar and energy transfer, a liquid outflow condition as well as the precipitation recharge and the heat radiation.

Different formulations have been studied for isothermal and non-isothermal compositional Darcy flows. They basically differ by their choice of the principal unknowns and equations and by the way they deal with phase transitions, which is one of the main difficulty of this type of models. The objective of such a choice is usually to reduce the non-linearity of the successive non-linear systems that typically arise when solving a transient problem with an Euler fully implicit time integration scheme.

Let us distinguish between variable switch and persistent variable formulations. The first ones adapt the set of principal unknowns and equations to the set of present phases which can vary in space and time.

The most well known formulation in this family is the so called natural variable or Coats' formulation widely used in reservoir simulations [START_REF] Coats | Implicit compositional simulation of single-porosity and dual-porosity reservoirs[END_REF][START_REF] Class | Numerical simulation of non-isothermal multiphase multicomponent processes in porous media.: 1. An efficient solution technique[END_REF][START_REF] Eymard | Vertex-centred discretization of multiphase compositional Darcy flows on general meshes[END_REF][START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal Darcy flows in fractured porous media[END_REF]. It has the advantage to use the physical variables of the thermodynamic and hydrodynamic laws as the set of principal unknowns. Its main drawback is to require a cumbersome switch of this set of variables depending on the set of present phases at each point of the space time domain. On the other hand, persistent variable formulations are based either on natural physical quantities such as overall component molar fractions or total specific enthalpy (see [START_REF] Voskov | Comparison of nonlinear formulations for two-phase multi-component EoS based simulation[END_REF]), or alternatively on nonstandard principal unknowns such as in [START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media: application to gas migration in a nuclear waste repository[END_REF][START_REF] Marchand | Fully coupled generalized hybrid-mixed finite element approximation of two-phase twocomponent flow in porous media. Part I: formulation and properties of the mathematical model[END_REF]. Another strategy to avoid the switch of variables is based on the extension of some physical quantities such as the phase molar fractions like in [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] using component fugacities, or the phase molar fractions and pressures like in [START_REF] Angelini | Finite volume approximation of a diffusion-dissolution model and application to nuclear waste storage[END_REF][START_REF] Masson | Formulations of two phase liquid gas compositional Darcy flows with phase transitions[END_REF]. Let us also mention the negative saturation formulations [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF][START_REF] Salimi | Negative Saturation Approach for Non-Isothermal Compositional Two-Phase Flow Simulations[END_REF] belonging to this family. A comparison between some of these formulations can be found in [START_REF] Voskov | Comparison of nonlinear formulations for two-phase multi-component EoS based simulation[END_REF][START_REF] Masson | Formulations of two phase liquid gas compositional Darcy flows with phase transitions[END_REF] in the case of isothermal compositional Darcy flows.

Our choice of the formulation is based, like in the Coats' formulation, on the phase pressures, phase saturations, temperature and phase molar fractions as set of principal unknowns. This is a convenient choice since all the physical laws can be directly expressed using subsets of this set of variables. It is also a very natural choice in single phase regions which are usually dominant in geothermal applications. In order to avoid the switch of variables, this choice of the principal unknowns is combined with an extension of the phase molar fractions of an absent phase by the molar fractions at thermodynamic equilibrium with the present phase. It results that the set of principal unknowns does not depend on the set of present phases. Moreover, the phase transitions can be expressed as complementarity constraints as in [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] which means that the non-linear systems can be solved using semi-smoothed Newton techniques such as the Newton-min algorithm [START_REF] Kräutle | The semi-smooth Newton method for multicomponent reactive transport with minerals[END_REF][START_REF] Ben Gharbia | Gas phase appearance and disappearance as a problem with complementarity constraints[END_REF]. This formulation has been preferred to the formulation proposed in [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] since the use of the component fugacities as principal unknowns rather than the phase molar fractions results in additional non-linear couplings between the molar fractions and the temperature which are not desirable for non-isothermal Darcy flows. In this work, the efficiency of this formulation combined with different improvements of the Newton-min non-linear solver will be investigated on several test cases. It will also be compared in terms of non-linear convergence with the formulation proposed in [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] over one simulation.

As mentioned in [START_REF] O'sullivan | Geothermal reservoir simulation: the state-of-practice and emerging trends[END_REF][START_REF] O'sullivan | Future directions in geothermal modelling[END_REF], the interaction between the flow in the porous medium and the atmosphere plays an important role in geothermal flows. Since the coupling between the porous medium and surface flows is not computationally realistic at the space and time scales of a geothermal flow, our objective is rather to model the soil-atmosphere interaction using an advanced boundary condition accounting for the matter (mole) and energy balance at the interface between the porous medium and the atmosphere. Such model should consider the vaporization of the liquid phase, the convective molar and energy transfer, a liquid outflow condition at seepage surfaces, as well as the heat radiation and the precipitation influx.

Assuming the vaporization of the liquid phase at the soil-atmosphere interface, the molar and energy normal fluxes at the interface on the atmosphere side are frequently approximated in hydrogeology by two-point fluxes between the gas phase at the interface and the atmosphere at a reference height [START_REF] Cui | Simulating the water content and temperature changes in an experimental embankment using meteorological data[END_REF][START_REF] Chalhoub | A simple heat and moisture transfer model to predict ground temperature for shallow ground heat exchangers[END_REF]. The transmissivities of these two-point fluxes are based on convective molar and energy transfer coefficients. Such approximation basically assumes that the lateral variations in wind, air temperature and humidity can be neglected [START_REF] Vanderborght | Heat and water transport in soils and across the soilatmosphere interface: 1. Theory and different model concepts[END_REF]. Let us refer to the textbook [START_REF] Monteith | Principles of Environmental Physics[END_REF] on meteorology for the computation of convective molar and energy transfer coefficients at the soil-atmosphere interface depending on the roughness of the soil surface including the effect of the vegetation, on the wind velocity, on the eddy diffusivity in the air stream and stability of the air above the heated soil surface. The radiation which is absorbed by and emitted from the soil surface as well as the precipitation recharge can also be incorporated in such models [START_REF] Cui | Simulating the water content and temperature changes in an experimental embankment using meteorological data[END_REF][START_REF] Chalhoub | A simple heat and moisture transfer model to predict ground temperature for shallow ground heat exchangers[END_REF].

Outflow boundary conditions are frequently used in hydrogeology at seepage surfaces allowing the groundwater discharge to occur where the water table intercepts a sloping land surface. They have already been used for geothermal applications as in [START_REF] Hurwitz | Groundwater flow, heat transport, and water table position within volcanic edifices: Implications for volcanic processes in the Cascade Range[END_REF] for a single component liquid gas Darcy flow model. For the Richards equation, outflow boundary conditions are modelled by complementarity constraints between the non negative liquid normal flux and the non negative capillary pressure (see [START_REF] Schweiser | Regularization of outflow problems in unsaturated porous media with dry regions[END_REF]). For liquid gas Darcy flow models, they are combined with a Dirichlet boundary condition for the gas pressure [START_REF] Lenzinger | Two-phase flow equations with outflow boundary conditions in the hydrophobic-hydrophilic case[END_REF]. To our knowledge, their extension to general non-isothermal compositional liquid gas Darcy flows has not yet been derived.

In this work, both the evaporation and liquid outflow models are combined in a single boundary condition which automatically switches from evaporation to evaporation and liquid outflow boundary condition. It assumes that the liquid phase does not accumulate at the surface on the atmosphere side, considering that standing water condition such as lake or sea can easily be expressed in the form of Dirichlet conditions. Alternatively, one would need to model the flow of the liquid phase at the soilatmosphere surface which has not been considered here and might induce a coupling with a system with different time scales than the underlying geothermal system. Our boundary condition is coupled with the general non-isothermal compositional liquid gas Darcy flow model. The previous formulation of the Darcy flow model is adapted to account for the new unknowns and equations at the soil-atmosphere interface. The derivation of the two-point molar and energy fluxes in the atmosphere is obtained starting from the transmission conditions proposed in [START_REF] Mosthaf | A coupling concept for two-phase compositional porous-medium and single-phase compositional free flow[END_REF] (see also [START_REF] Masson | Coupling compositional liquid gas Darcy and free gas flows at porous and free-flow domains interface[END_REF][START_REF] Vanderborght | Heat and water transport in soils and across the soilatmosphere interface: 1. Theory and different model concepts[END_REF]) for the coupling of non-isothermal compositional liquid gas Darcy and free gas flows. The complementarity constraint for the liquid outflow is extended to non-isothermal compositional flows using a switching criterion based on the thermodynamic equilibrium between the gas and liquid phases at the interface on the free-flow side.

The structure of the rest of this article is as follows. The non-isothermal compositional two-phase Darcy flow model and its formulation are presented in Section 2. The soil-atmosphere boundary condition is derived in Section 3. Then, the finite volume two-point flux discretization of the model as well as the Newton-min algorithms used to solve the fully coupled systems at each time step of the simulation are introduced in Section 4. The convergence of the model is then studied over a semi-analytical 1D stationary solution in Section 5. In Section 6, the soil-atmosphere evaporation boundary condition is compared with a full-dimensional free gas flow and transport model coupled to the liquid gas Darcy flow. In Section 7, the model and its formulation are studied numerically in terms of solution and convergence of the Newton-min non-linear solvers on several geothermal test cases including 1D test cases and 2D cuts of the Bouillante geothermal field in Guadeloupe.

Non-isothermal compositional two-phase Darcy flow model

We consider a non-isothermal compositional liquid gas Darcy flow model with P = {g, l} denoting the set of gas and liquid phases. Each phase α ∈ P is a mixture of a set of components denoted by C including typically a water component, denoted w, which can vaporize in the gas phase and a set of gaseous components which can dissolve in the liquid phase. The thermodynamic properties of each phase α ∈ P depend on its pressure P α , the local equilibrium temperature of the system T and its molar fractions

C α = (C α i ) i∈C .
For each phase α ∈ P, we denote by ζ α (P α , T,C α ) its molar density, by ρ α (P α , T,C α ) its mass 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. For the gas phase, assuming an ideal mixture, the molar enthalpy is defined by

h g (P g , T,C g ) = ∑ i∈C C g i h g i (P g , T ), (2.1)
where h g i (P g , T ) is the molar enthalpy of the component i in the gas phase. Thermodynamic equilibrium between the gas and liquid phases is assumed for each component and governed by the phase fugacities denoted by f α (P α , T,C α ) = ( f α i (P α , T,C α )) i∈C , α ∈ P. The rock porosity is denoted by ϕ (x) and the rock permeability tensor by Λ Λ Λ (x) where x denotes the spatial coordinates. The hydrodynamic Darcy laws are characterized by the relative permeability k α r (x, S α ) of each phase α ∈ P, as a function of the phase saturation S α , and by the capillary pressure P c (x, S g ) = P g -P l . The relative permeabilities and capillary pressure are piecewise constant for each rocktype, thus their dependence on x is omitted in the following for the sake of simplicity.

Our formulation of the model (called T-PSC in the following) is based on the fixed set of unknowns defined by X =

( P α , T, S α ,C α , α ∈ P ) . (2.2)
Note that, as opposed to the Coats' variable switch formulation [START_REF] Coats | Implicit compositional simulation of single-porosity and dual-porosity reservoirs[END_REF][START_REF] Class | Numerical simulation of non-isothermal multiphase multicomponent processes in porous media.: 1. An efficient solution technique[END_REF][START_REF] Eymard | Vertex-centred discretization of multiphase compositional Darcy flows on general meshes[END_REF][START_REF] Xing | Parallel numerical modeling of hybrid-dimensional compositional non-isothermal Darcy flows in fractured porous media[END_REF], the molar fractions C α of an absent phase α are extended by the ones at equilibrium with the present phase in the sense that the equality of the gas and liquid fugacities f g (P g , T,C g ) = f l (P l , T,C l ) always holds.

Let n i (X) be the number of moles of the component i ∈ C per unit pore volume defined by

n i (X) = ∑ α∈P ζ α S α C α i , i ∈ C .
We introduce the rock energy per unit rock volume defined by E r (T ) and the fluid energy per unit pore volume defined by

E(X) = ∑ α∈P ζ α S α e α .
Let us denote by g the gravitational acceleration vector. The generalized Darcy velocity of the phase α ∈ P is given by

V α = - k α r µ α Λ Λ Λ (x) ( ∇P α -ρ α g ) . ( 2.3)
The total molar flux of the component i ∈ C is denoted by q i and the energy flux by q e , with

q i = ∑ α∈P C α i ζ α V α , q e = ∑ α∈P h α ζ α V α -λ ∇T, (2.4)
where λ stands for the bulk thermal conductivity of the fluid and rock mixture.

The system of equations accounts for the molar conservation of each component i ∈ C together with the energy conservation

ϕ (x)∂ t n i + div(q i ) = 0, i ∈ C , ϕ (x)∂ t E + (1 -ϕ (x))∂ t E r + div(q e ) = 0. (2.5)
It is complemented by the following capillary relation between the two phase pressures and the pore volume balance

{ P c (S g ) = P g -P l , ∑ α∈P S α = 1. (2.6)
In the spirit of [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF][START_REF] Masson | Formulations of two phase liquid gas compositional Darcy flows with phase transitions[END_REF], the liquid gas thermodynamic equilibrium can be expressed as the following complementarity constraint for each phase α ∈ P combined with the equality of the gas and liquid fugacities of each component

{ S α ≥ 0, 1 -∑ i∈C C α i ≥ 0, S α (1 -∑ i∈C C α i ) = 0, α ∈ P, f g i (P g , T,C g ) = f l i (P l , T,C l ), i ∈ C . (2.7)
Note that our formulation of the model leads, independently on the set of present phases, to the fix sets of 2#C + 5 unknowns (2.2) and of 2#C + 5 equations (2.5)-(2.6)-(2.7) including the #C + 1 conservation equations (2.5) and the remaining #C + 4 local closure laws (2.6)-(2.7). It also has the advantage to express the thermodynamic equilibrium as complementarity constraints which will allow the use of semi-smooth Newton methods [START_REF] Kräutle | The semi-smooth Newton method for multicomponent reactive transport with minerals[END_REF][START_REF] Ben Gharbia | Gas phase appearance and disappearance as a problem with complementarity constraints[END_REF] to solve the non-linear systems at each time step of the simulation as specified in Section 4.1.

Soil-atmosphere boundary condition for non-isothermal compositional liquid gas Darcy flow

The fluid and energy transport in high energy geothermal systems is deeply governed by the conditions set at the boundary of the computational domain. In particular, it is well known that the modelling of the interaction between the porous medium model and the atmosphere plays an important role [START_REF] O'sullivan | Geothermal reservoir simulation: the state-of-practice and emerging trends[END_REF][START_REF] O'sullivan | Future directions in geothermal modelling[END_REF]. In this section, a boundary condition is derived to model the soil-atmosphere interaction based on mole and energy balance equations set at the interface. The model takes into account two coupling processes: on the one hand, the vaporization of the liquid phase and the convective molar and energy transfer in the atmosphere described in Subsection 3.1, on the other hand, a liquid outflow condition described in Subsection 3.2. Both coupling processes will be combined in a single boundary condition assuming that the liquid phase does not accumulate at the surface. The radiation and the precipitation recharge are also considered.

3.1 Convective molar and energy transfer in the atmosphere

Transmission conditions at the interface between a non-isothermal compositional liquid gas Darcy flow and a gas free flow

The derivation of the boundary condition accounting for convective molar and energy transfer in the atmosphere can be explained starting from the transmission conditions introduced in [START_REF] Mosthaf | A coupling concept for two-phase compositional porous-medium and single-phase compositional free flow[END_REF] (see also [START_REF] Masson | Coupling compositional liquid gas Darcy and free gas flows at porous and free-flow domains interface[END_REF][START_REF] Vanderborght | Heat and water transport in soils and across the soilatmosphere interface: 1. Theory and different model concepts[END_REF]) to couple a non-isothermal compositional liquid gas Darcy flow with a gas free flow. These conditions state the continuity of the component molar and energy normal fluxes, assuming instantaneous vaporization of the liquid phase, as well as the continuity of the gas molar fractions, of the temperature and of the gas pressure, neglecting the gas pressure jump. It is complemented by the thermodynamic equilibrium between the liquid and gas phases and by the Beavers-Joseph condition. On the free-flow side, the component molar and energy fluxes are defined by

w i = ζ g (P, T,C) ( C i u -D t ∇C i ) , i ∈ C , w e = ζ g (P, T,C)h g (P, T,C)u -λ t ∇T -∑ i∈C ζ g (P, T,C)h g i (P, T )D t ∇C i , (3.1)
where u denotes the gas velocity, P the pressure, C = (C i ) i∈C the gas molar fractions, T the temperature, D t the turbulent diffusivity and λ t the turbulent thermal conductivity. The continuity of the component molar normal fluxes states that

w i • n = q i • n, i ∈ C , ( 3.2)
where the unit normal vector n at the interface is oriented outward from the porous medium domain. The last term in the free-flow energy flux in (3.1) introduces a strong non-linear coupling between the component molar and energy fluxes which raises an additional difficulty in the two-point approximation of the normal fluxes. This can be addressed in a simple and efficient way if the dissolution of the gaseous components in the liquid phase is small which corresponds to the usual case. In such a case, using

ζ g << ζ l , C w << 1, ∑ i∈C C i = 1 and (3.
2), we can derive that

|q i • n| << |q w • n| ∼ ζ g |u • n| ∼ ζ g D t |∇C w • n|, (3.3)
for all i ∈ C \ {w}. Using that ∑ i∈C C i = 1, one has , temperature T atm ∞ and pressure P atm , which fixes the far field atmospheric specific gas enthalpy of the water component h g,atm w,∞ = h g w (P atm , T atm ∞ ). From the transmission conditions stated above, the temperature, the gas molar fractions and the gas pressure defined at the interface on the atmospheric side, match with their values on the porous medium side and consequently they are denoted respectively by T , C g and P g . The twopoint flux approximations account for the turbulent boundary layers of the gas flow and transport in the atmosphere using the molar and energy transfer coefficients H m and H T . These coefficients are usually obtained from correlations used for environmental gas flows depending on the roughness of the soil surface including the effect of the vegetation, on the wind velocity, on the eddy diffusivity in the air stream and stability of the air above the heated soil surface [START_REF] Monteith | Principles of Environmental Physics[END_REF]. The two-point fluxes also take into account the convective normal fluxes using, as additional unknown, the gas molar flow rate q g,atm at the interface on the atmosphere side oriented outward from the porous medium domain. It is combined with an upwinding of the gas molar fractions and of the gas enthalpy of the water component between the interface and the far field atmospheric conditions. This leads us to define the following two-point fluxes oriented outward from the porous medium domain

(w e + λ t ∇T ) • n = ζ g (P, T,C) ( h g (P, T,C)u • n -∑ i∈C h g i (P, T )D t ∇C i • n ) = ζ g (P, T,C) ( ∑ i∈C h g i (P, T )C i u • n -∑ i∈C \{w} (h g i (P, T ) -h g w (P, T ))D t ∇C i • n ) = h g w (P, T )ζ g (P, T,C)u • n + ∑ i∈C \{w} (h g i (P, T ) -h g w (P, T ))w i • n. ( 3
q g,atm i = (q g,atm ) + C g i -(q g,atm ) -C g,atm i,∞ + H m ( C g i -C g,atm i,∞
) , i ∈ C , q g,atm e = (q g,atm ) + h g w (P g , T ) -(q g,atm

) -h g,atm w,∞ + H T (T -T atm ∞ ), (3.6)
where for any real u, we have set (u) + = max(0, u) and (u) -= max(0, -u).

Neglecting the variations of pressure in the atmosphere leads to the following continuity equation for the gas pressure

P g = P atm . (3.7)
Thermodynamic equilibrium is always assumed at the interface in the sense that the gas molar fractions and pressure at the interface on the porous medium side are extended by the one at equilibrium with the liquid phase in the absence of the gas phase. On the other hand, the liquid phase can appear or disappear according to the liquid phase complementarity constraint. It results that the following equations hold at the interface

                 f g i (P g , T,C g ) = f l i (P l , T,C l ), i ∈ C , ∑ i∈C C g i = 1, S l ≥ 0, 1 -∑ i∈C C l i ≥ 0, S l (1 -∑ i∈C C l i ) = 0, S g = S g (P g -P l ), ∑ α∈P S α = 1, (3.8)
where P l is the liquid pressure, C l the liquid molar fractions and S α , α ∈ P the saturations at the interface on the porous medium side and S g (P c ) denotes the inverse of the monotone graph extension of the capillary pressure function P c (S g ). As detailed in [START_REF] Brenner | Improving Newton's Method Performance by Parametrization: The Case of the Richards Equation[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] and in Subsection 7.1.2, a switch of variables between S g and P c could also be used in order to account for non invertible capillary functions.

Regarding the interface energy balance, the model can also account for the solar and long wave radiation that is absorbed by and emitted from the soil surface defined by the following net radiation R n ( W. m -2 )

R n = (1 -a)R s + R a -εσ SB T 4 , ( 3.9)
where R a ( W. m -2 ) is the incoming long-wave radiation emitted by the atmosphere, R s ( W. m -2 ) is the net short-wave radiation, a is the surface albedo, σ SB ( W. m -2 . K -4 ) is the Stephan-Boltzmann constant and ε the soil emissivity.

Liquid outflow complementarity constraint

The liquid phase is assumed to vaporize instantaneously when leaving the porous medium as long as the atmosphere is not saturated with water vapour. As soon as the atmosphere is vapour saturated at the interface, the component molar and energy normal fluxes in the liquid phase defined by

q l,atm i = C l,atm i q l,atm , i ∈ C , q l,atm e = h l (P l , T,C l,atm )q l,atm , (3.10)
are allowed to exit the porous medium, where q l,atm ≥ 0 is an additional unknown corresponding to the total liquid molar flow rate oriented positively outward to the porous medium domain. It is assumed that the liquid phase does not accumulate at the surface on the atmosphere side, because modelling the flow of the liquid phase at the soil-atmosphere surface might induce a coupling with a system with different time scales than the underlying geothermal system. Moreover, the accumulation of water such as lake or sea can easily be expressed in the form of Dirichlet conditions. In (3.10), the liquid molar fractions C l,atm = (C l,atm i

) i∈C at the interface on the atmosphere side are those at thermodynamic equilibrium with the gas phase thus are such that f l (P atm , T,C l,atm ) = f g (P g , T,C g ).

(3.11)

Note that, due to the jump of the capillary pressure which vanishes on the atmosphere side, C l,atm does not match in general with the liquid molar fractions on the porous medium side C l which satisfies

f l (P l , T,C l ) = f g (P g , T,C g ). (3.12)
The liquid molar outflow rate q l,atm is determined by the following complementarity constraint accounting for the thermodynamic equilibrium between the liquid and gas phases at the interface on the atmosphere side

   (1 -∑ i∈C C l,atm i ) q l,atm = 0, 1 -∑ i∈C C l,atm i ⩾ 0, q l,atm ⩾ 0. (3.13)
It remains to eliminate the liquid molar fractions C l,atm from (3.10) and (3.13). Let us consider for f ∈ R C the function C l ( f , P l , T ) ∈ R C defined as the unique solution of the equation f l (P l , T,C l ) = f . From f g (P g , T,C g ) = f l (P g , T,C l,atm ) = f l (P l , T,C l ) := f given by the equations (3.11) and (3.12), it results that

C l,atm = C l ( f , P g , T ).
On the one hand, if S l > 0, it follows that

1 -∑ i∈C C l,atm i = ∑ i∈C ( C l i -C l,atm i ) = ∑ i∈C ( C l i ( f , P l , T ) -C l i ( f , P g , T )
) .

(3.14)

Following [START_REF] Masson | Formulations of two phase liquid gas compositional Darcy flows with phase transitions[END_REF], since the function ∑ i∈C C l i ( f , P, T ) is strictly decreasing with respect to P, it results that the complementarity constraint (3.13) is equivalent to

{ (P g -P l ) q l,atm = 0, P g -P l ⩾ 0, q l,atm ⩾ 0. (3.15)
On the other hand, if S l = 0 then one has P g -P l = P c (1) > 0 and ∑ i∈C C l,atm i < 1. It results that both conditions (3.15) and (3.13) imply that q l,atm = 0. Finally, let us remark that if (3.15) holds, the liquid outflow fluxes in (3.10) rewrite as follows

q l,atm i = C l i q l,atm , i ∈ C , q l,atm e = h l (P l , T,C l )q l,atm . (3.16)
The model also takes into account the following component molar and energy flow rates which represent the precipitation recharge

q l,rain i = C l,rain i q l,rain , i ∈ C , q l,rain e = h l (P atm , T atm ∞ ,C l,rain )q l,rain , (3.17)
with the rain molar fractions denoted by C l,rain = (C l,rain i ) i∈C , a temperature assumed at equilibrium with the far field atmosphere, the rain molar enthalpy denoted by h l,rain = h l (P atm , T atm ∞ ,C l,rain ) and a non positive rain molar flow rate q l,rain ≤ 0.

Evaporation-outflow boundary condition

Both the liquid outflow and evaporation models are combined in a single boundary condition, assuming that the liquid does not accumulate at the surface. Gathering the equations (3.6), (3.7), (3.8), (3.15), (3.16), (3.17) together with the component molar an energy balance equations, the evaporation-outflow boundary condition at the interface is defined by the sets of 7 + 2#C unknowns

X Γ = ( q g,atm , q l,atm , T, P α , S α ,C α , α ∈ P ) ,
and equations

                                                         q i • n = (q g,atm ) + C g i -(q g,atm ) -C g,atm i,∞ + H m ( C g i -C g,atm i,∞
)

+C l i q l,atm +C l,rain i q l,rain , i ∈ C , q e • n = (q g,atm ) + h g w (P g , T ) -(q g,atm ) -h g,atm w,∞ + H T (T -T atm ∞ ) -R n + h l (P l , T,C l )q l,
atm + h l,rain q l,rain ,

P g = P atm , S g = S g (P g -P l ), ∑ α∈P S α = 1, ∑ i∈C C g i = 1, min ( S l , 1 -∑ i∈C C l i ) = 0, f g i (P g , T,C g ) = f l i (P l , T,C l ), i ∈ C min ( P g -P l , q l,atm ) = 0. (3.18)
4 Discretization and non-linear solvers

The system of equations (2.5)-(2.6)-(2.7)-(3.18) is discretized using a finite volume discretization in space with a Two-Point Flux Approximation (TPFA) of the Darcy and Fourier fluxes [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF] combined with a phase based upwind scheme for the approximation of the mobilities [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF][START_REF] Eymard | Finite volume methods[END_REF][START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. A mesh satisfying the admissibility condition of TPFA schemes at both inner and boundary faces is used [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF]. It can be typically a triangular mesh with acute angles and isotropic permeability, a Voronoi mesh with isotropic permeability or a Cartesian mesh with anisotropic permeability aligned with the axes. In all cases the permeability is assumed cellwise constant.

Let Ω denote a bounded polytopal domain of R d (polygonal for d = 2 or polyhedral for d = 3) and Γ ⊂ ∂ Ω the boundary on which the soil-atmosphere evaporation-outflow boundary condition is imposed. Let us denote by M the set of cells that are disjoint open polytopal subsets of Ω such that

∪ K∈M K = Ω .
It is assumed that there exists a subset F Γ of the set of faces such that

Γ = ∪ σ ∈F Γ σ .
For the sake of simplicity only the boundary faces σ ∈ F Γ are considered in the following. The time integration is based on a fully implicit Euler scheme to avoid severe restrictions on the time steps. For N ∈ N * , let us consider the time discretization

t 0 = 0 < t 1 < • • • < t n-1 < t n < • • • < t N = t f of the time interval [0,t f ]. We denote the time steps by ∆t n = t n -t n-1 for all n = 1, • • • , N.
Let us denote the set of physical unknowns by

X K = ( P α K , T K , S α K ,C α K , α ∈ P ) , ( 4.1)
for each cell K ∈ M and by

X σ = ( q g,atm σ , q l,atm σ , P α σ , T σ , S α σ ,C α σ , α ∈ P ) , ( 4.2)
for each boundary face σ ∈ F Γ . The full set of unknowns is denoted by

X D = {X K , X σ , K ∈ M , σ ∈ F Γ }. (4.3) For each degree of freedom ν ∈ M ∪ F Γ , we denote by R ν (X D ) the residual vector ( R ν,i (X D ), i ∈ C ∪{e} )
of the component and energy conservation equations and by L ν (X ν ) the residual vector of the local closure laws. It defines the following non-linear system at each time step n = 1, ..., N

0 = R(X D ) =        ( R K (X D ) L K (X K ) ) K ∈ M , ( R σ (X D ) L σ (X σ ) ) σ ∈ F Γ , (4.4)
where the current time step superscript n has been dropped.

Newton-min non-linear solver

The non-linear system R(X D ) = 0 is solved using a Newton-min solver [START_REF] Kräutle | The semi-smooth Newton method for multicomponent reactive transport with minerals[END_REF][START_REF] Ben Gharbia | Gas phase appearance and disappearance as a problem with complementarity constraints[END_REF] as detailed below. In order to reduce the size of the linear systems to be solved at each Newton-min iteration to #C + 1 equations and unknowns for each degree of freedom ν ∈ M ∪ F Γ , the set of unknowns X ν is split into #C + 1 primary unknowns X p ν and the remaining secondary unknowns X s ν . This splitting is done for each degree of freedom in such a way that the Jacobian of the local closure laws ∂ L ν ∂ X s ν (X ν ) with respect to the secondary unknowns is non singular.

The Newton-min algorithm is initialized with an initial guess X (0) D usually given by the previous time step solution and iterates on the following steps for r = 0, • • • , until the following stopping criterion is satisfied

max i∈C ∪{e}     ∑ ν∈M ∪F Γ |R ν,i (X (r) D )| ∑ ν∈M ∪F Γ |R ν,i (X (0) D )|     ≤ ε R or dim(X) ∑ i=1    max ν∈M |dX (r) ν,i | ∆ X i    + dim(X Γ ) ∑ i=1    max ν∈F Γ |dX (r) ν,i | ∆ X i    ≤ ε X with dX (r) ν,i = { dX p (r) ν,i if i is a primary unknown, dX s (r) ν,i if i is a secondary unknown. and ε R = 10 -8 , ε X = 10 -6 , given ∆ X i > 0, i = 1, • • • , dim(X ν )
and with the Newton's steps dX

s (r) ν,i , dX p (r) ν,i
defined respectively in (4.5) and (4.6). If the convergence it not met after r max Newton iterations, the time step is chopped.

1. Computation of the residual R(X (r) D ) and of the Jacobian matrix with elimination of the secondary unknowns. It starts with the choice of the primary and secondary unknowns for each degree of freedom ν ∈ M ∪ F Γ depending only on the active complementarity constraints, choice specified in Table (4.1). Then the matrices A sp (r) ν and the vectors B s (r)

ν , ν ∈ M ∪ F Γ defined by A sp (r) ν = - ( ∂ L ν ∂ X s (r) ν ) -1 ∂ L ν ∂ X p (r) ν , B s (r) ν = ( ∂ L ν ∂ X s (r) ν ) -1 L ν ,
and such that dX

s (r) ν = A sp (r) ν dX p (r) ν + B s (r) ν , ( 4.5)
are computed to obtain the reduced Jacobian J

(r) = (J (r) ν,ν ′ ) (ν,ν ′ )∈(M ∪F Γ ) 2 defined by the square ma- trices of size #C + 1 J (r) ν,ν ′ = ∂ R ν ∂ X p (r) ν ′ (X (r) D ) + ∂ R ν ∂ X s (r) ν ′ (X (r) D )A sp (r) ν ′
B (r) ν = -R ν (X (r) D ) -∑ ν ′ ∈M ∪F Γ ∂ R ν ∂ X s (r) ν ′ (X (r) D )B s (r) ν ′ .
2. Computation of the solution of the reduced linear system

J (r) dX p (r) D = B (r) . (4.6) 3. Update of the unknowns X (r) ν , ν ∈ M ∪ F Γ with a possible relaxation θ (r) ∈ (0, 1] { X p (r+1) ν = X p (r) ν + θ (r) dX p (r) ν , X s (r+1) ν = X s (r) ν + θ (r) ( A sp (r) ν dX p (r) ν + B s (r) ν ) . (4.7)
4. Additional updates of some unknowns in order to satisfy exactly some non-linear closure laws to be specified.

The step 4 of the above algorithm allows proposing different improvements of the basic Newton-min algorithm that are detailed in the following paragraphs. Table 4.1 Choice of the primary unknowns of the degree of freedom ν ∈ M ∪ F Γ depending on the active complementarity constraints of the Newton-min algorithm.

ν ∈ F Γ ν ∈ M q l,atm ν = 0 q g,atm ν , P c,ν , (C l i,ν ) i=1,#C -1 1 -∑ i∈C C g i,ν = 0 P g ν , S g ν , (C l i,ν ) i=1,#C -1 1 -∑ i∈C C l i,ν = 0 1 -∑ i∈C C l i,ν = 0 P g ν -P l ν = 0 q g,atm ν , q l,atm ν , T ν , S g ν = 0 P g ν , T ν , (C l i,ν ) i=1,#C -1 1 -∑ i∈C C l i,ν = 0 (C l i,ν ) i=1,#C -2 1 -∑ i∈C C l i,ν = 0 q l,atm ν = 0 q g,atm ν , T ν , (C g i,ν ) i=1,#C -1 1 -∑ i∈C C g i,ν = 0 P g ν , T ν , (C g i,ν ) i=1,#C -1 S l ν = 0 S l ν = 0 4.1.

Basic Newton-min algorithm

The basic version of the Newton-min algorithm only enforces at each iterate the following non-linear closure law for ν ∈ M ∪ F Γ P g ν -P l ν = P c (S g ν ).

Our objective for this basic Newton-min algorithm was to use no projections of the physical unknowns onto their physical bounds. However, in order to obtain the convergence of the Newton-min algorithm, it was necessary to project the molar fractions of a present phase within the range say [-0.2; 1.2] at each Newton iterate.

Newton-min with projection on the complementarity constraints

In order to obtain a better convergence of the Newton-min algorithm, all the complementarity constraints of type min(X 1 , X 2 ) = 0 are enforced at the initial guess and at each Newton iterate. In addition, P g ν -P l ν = P c (S g ν ) is also enforced at each Newton iterate and the following physical ranges are imposed on the molar fractions of a present phase and on the saturations

     if S α ν > 0 then 0 ≤ C α i,ν ≤ 1, i ∈ C , α ∈ P, S α ν ⩾ 0, α ∈ P, ∑ α∈P S α ν = 1.
An additional improvement also studied in the numerical section is to test the appearance of a missing phase using the molar fractions at equilibrium with the present phase rather than their linear Newton updates.

Newton-min with projection on the complementarity constraints and thermodynamic equilibrium

In addition to the previous updates, the molar fractions C ᾱ which are secondary unknowns (see Table 4.1), complemented by the temperature if both phases are present, are updated in order to verify the following closure laws at each Newton iterate

   min ( S ᾱ ν , 1 -∑ i∈C C ᾱ i,ν ) = 0, f g i (P g ν , T ν ,C g ν ) = f l i (P l ν , T ν ,C l ν ), i ∈ C
for all ν ∈ M ∪ F Γ . Note that the first equation is already satisfied if only one phase is present as S ᾱ ν = 0.

5 Convergence of the Finite Volume scheme to a semi-analytical 1D stationary solution

In this section, a simplified 1D non-isothermal liquid gas model with a single water component is used to analyse the convergence of the TPFA scheme with the T-PSC formulation to a semi-analytical stationary solution.

Let us consider a single water component (C = {w}) liquid gas non-isothermal Darcy flow model on the 1D domain Ω = (z 0 , z 3 ), z 3 > z 0 . Let us assume constant liquid and gas molar densities ζ l and ζ g , and viscosities µ l and µ g . The liquid and gas molar enthalpies are defined by h α (T ) = c α p T , α ∈ P, with the molar heat capacities c g p > c l p . Let us denote by P sat (T ) the saturated vapour pressure and by T sat (P) its inverse. The thermal conductivity is assumed constant λ > 0. The relative permeabilities k α r (S α ), α = l, g are such that k α r (0) = 0 and k α r (1) = 1. The capillary pressure is neglected. The absolute permeability of the porous medium Λ Λ Λ is also assumed constant.

Let us recall the definition of the water (or total) molar flow rate

q w = ζ l V l + ζ g V g , with V α = -k α r (S α )
µ α Λ Λ Λ ∂ z (P + ρ α |g|z) and ρ α = ζ α m w where m w denotes the water molar mass. We also set S = S g = 1 -S l . The energy flow rate is

q e = c l p T ζ l V l + c g p T ζ g V g -λ ∂ z T.
The stationary solution P(z), T (z), S(z), z ∈ Ω is such that ∂ z q w = 0, ∂ z q e = 0 and

   P g = P sat (T ), if 0 < S < 1, P g ≥ P sat (T ), if S = 0, P g ≤ P sat (T ), if S = 1.
The boundary conditions specify at the bottom boundary z = z 0 the input saturation S = 0, the input molar flow rate q w,0 > 0 and the input energy flow rate q e,0 . The temperature at z = z 0 is denoted by T 0 and the pressure by P 0 . It is assumed that P 0 > P sat (T 0 ) meaning that only the liquid phase is present at the bottom boundary. At the top boundary we impose a gaseous state of temperature T 3 and pressure P 3 which are such that P 3 < P sat (T 3 ).

We assume that the solution has three zones defined by z 0 < z 1 < z 2 < z 3 and such that S = 0 on (z 0 , z 1 ), S = 1 on (z 2 , z 3 ) and P = P sat (T ) on (z 1 , z 2 ).

Let us denote by P i , T i the pressure and temperature at z = z i for i = 0, 1, 2, 3. The saturation S 1 at point z 1 is continuous and equal to S 1 = 0. At point z 2 the saturation may be discontinuous with S + 2 = 1 and S - 2 to be determined. We define the following set of 8 unknowns U = (P 0 , T 0 , P 1 , T 1 , P 2 , T 2 , z 1 , z 2 ) for which we are going to define 8 equations.

Since S = 0 on [z 0 , z 1 ],

q w,0 = - ζ l µ l Λ Λ Λ ∂ z (P + ρ l |g|z).
Integrating this equation between z 0 and z 1 , we obtain

P 1 -P 0 = -(ρ l |g| + µ l q w,0 ζ l Λ Λ Λ )(z 1 -z 0 ).
Moreover, ∂ z q e = 0 gives

a l ∂ z T = ∂ z 2 T
with a l = c l p q w,0 λ , from which we deduce that T (z) = E l + D l e a l z , with

α l i = e a l z i , i = 0, 1, E l = α l 1 T 0 -α l 0 T 1 α l 1 -α l 0 , D l = T 1 -T 0 α l 1 -α l 0 .
Introducing this formula into q e,0 = c l p T q w,0 -λ ∂ z T, we get q e,0 = c l p E l q w,0 + D l e a l z (c l p q w,0 -λ a l ) = c l p q w,0 E l . Likewise on ]z 2 , z 3 ], we have S = 1 from which we deduce that

P 3 -P 2 = -(ρ g |g| + µ g q w,0 ζ g Λ Λ Λ )(z 3 -z 2 ) and T (z) = E g + D g e a g z , with a g = c g p q w,0 λ , α g i = e a g z i , i = 2, 3, E g = α g 3 T 2 -α g 2 T 3 α g 3 -α g 2 , D g = T 3 -T 2 α g 3 -α g 2 .
It results that q e,0 = c g p q w,0 E g . This provides the following 6 equations

                           P 1 = P sat (T 1 ), P 2 = P sat (T 2 ), P 1 -P 0 = -(ρ l |g| + µ l q w,0 ζ l Λ Λ Λ )(z 1 -z 0 ), P 3 -P 2 = -(ρ g |g| + µ g q w,0 ζ g Λ Λ Λ )(z 3 -z 2 ),
c l p q w,0 E l = q e,0 , c g p q w,0 E g = q e,0 .

The two missing equations are obtained by integration of the following system on (z

1 , z 2 )                    q α = - k α r (S α ) µ α Λ Λ Λ ∂ z (P + ρ α |g|z), ζ l q l + ζ g q g = q w,0 , c l p T ζ l q l + c g p T ζ g q g -λ ∂ z T = q e,0 , P = P sat (T ), P(z 2 ) = P 2 , (5.1)
and by setting S(z 1 ) = 0 and P(z 1 ) = P 1 .

We also obtain an algebraic equation using that ζ l q l = q w,0 , q g = 0 at z = z 1 as well as S(z 1 ) = 0, T (z 1 ) = T 1 . It leads by elimination of ∂ z T and using ∂ z P = P ′ sat (T )∂ z T to the following equation which allows computing directly T 1

- ζ l Λ Λ Λ µ l ( P ′ sat (T 1 )
1 λ (c l p T 1 q w,0q e,0 ) + ρ l |g| ) = q w,0 .

To simplify the computation of the solution we assume that e -a l (z 1 -z 0 ) and e -a g (z 3 -z 2 ) are of the order of the machine round off error. This will be checked a posteriori but it basically holds when λ is small enough. In that case we have E l = T 0 and E g = T 2 , which gives T 0 = q e,0 c l p q w,0

, T 2 = q e,0 C g p q w,0

, P 2 = P sat (T 2 ) and z 2 = z 3 -P 2 -P 3

ρ g |g| + µ g q w,0 ζ g Λ Λ Λ
.

Then, the temperature T 1 can be computed as the solution of the equation

- ζ l Λ Λ Λ µ l ( P ′ sat (T 1 ) 1 λ (c l p T 1 q w,0 -q e,0 ) + ρ l |g| ) = q w,0 ,
as well as P 1 = P sat (T 1 ).

The position z 1 is computed by dichotomy in order to solve T (z 1 ) = T 1 (or equivalently S(z 1 ) = 0) using a numerical integration of (5.1) on (z 1 , z 2 ). This integration is done using an Euler implicit scheme. It amounts to solve a non-linear system for T (z), S(z) at each node z of a uniform discretization of the interval (z 1 , z 2 ) in decreasing z order starting from z = z 2 . Once z 1 is obtained, we can compute P 0 as follows

P 0 = P 1 + (ρ l |g| + µ l q w,0 ζ l Λ Λ Λ )(z 1 -z 0 ).
Let us consider the following data set

z 0 = 0 m, z 3 = 3000 m, T 3 = 400 K, P 3 = Pa, ζ g = 10 m w , ζ l = 1000 m w , µ g = 10 -4 Pa. s -1 , µ l = -3 Pa. s -1 , λ = 1 W. m -1 . K -1 , k α r (S α ) = (S α ) 2 , c l p = 2000m w , c g p = 2800m w , q w,0 = 6 10 -4
m w , q e,0 = 600 q w,0 c l p , |g| = 9.81 m. s -2 , m w = 0.018 kg. mol -1 ,

Λ Λ Λ = 10 -12 m 2 .
The saturated vapour pressure is defined by the following Rankine formula

P sat = 10 5 exp(13.7 - 5120 T ).
The corresponding stationary solution computed as described above is exhibited in Figure 5.1. Then the convergence of the finite volume approximation of the transient liquid gas Darcy flow model to this stationary solution is tested. The simulation on a given mesh is stopped once a stationary solution is obtained. Figure 5.2 exhibits the good convergence of the errors obtained for the pressure, temperature and gas saturation on three uniform meshes of sizes N = 50, 150, 450. porous medium domains are those introduced in [START_REF] Mosthaf | A coupling concept for two-phase compositional porous-medium and single-phase compositional free flow[END_REF]. They assume the vaporization of the liquid phase in the free-flow domain and account for the gas molar fraction and molar and energy normal flux continuity, the liquid gas thermodynamic equilibrium, the no slip condition and the normal component of the normal stress continuity.

We consider a 2D test case from [START_REF] Birgle | A domain decomposition method to couple non-isothermal compositional gas liquid Darcy and free gas flows[END_REF] which simulates the mass and energy exchanges occurring within deep geological radioactive waste disposal at the interface between a geological formation with low permeable porous medium and a ventilation excavated gallery. The data set is derived from lab experiments and in accordance with the deep disposal center for French radioactive waste project. In this test case, the porous medium initially saturated with the liquid phase is dried by suction in the neighbourhood of the interface between the porous and free-flow domains. The gas phase penetrates the porous domain and the liquid phase is vaporized in the free-flow domain.

As exhibited in Figure 6.1, the porous medium domain is defined by Ω pm = (0, l) × (h f f , h pm ) with l = 100 m, h f f = 5 m, h pm = 20 m. It corresponds to the computational domain of the Darcy flow model coupled with the soil-atmosphere evaporation-outflow boundary condition at the interface Γ = (0, l) × {h f f }. The computational domain (0, l) × (0, h pm ) of the coupled Darcy and full-dimensional free-flow models is the union of the porous medium domain Ω pm and of the free-flow domain Ω f f = (0, l) × (0, h f f ).

A single rocktype defined by the Callovo Oxfordian clay (Cox) is considered in the porous medium with the homogeneous porosity ϕ (x) = 0.15 and isotropic permeability Λ Λ Λ (x) = K × I with K = 5 • 10 -20 m 2 . The relative permeabilities and capillary pressure are given by the following Van Genuchten laws with the parameters n = 1.49, m = 1 -1 n , Pr = 15 • 10 6 Pa and the residual liquid and gas saturations S l r = 0.4, S g r = 0 (see Figure 6.3) The liquid and gas phases are a mixture of two components, the water denoted by w and the air denoted by a. The gas molar density is defined by the perfect gas law ζ g = P g RT , with R = 8.314 J. K -1 . mol -1 and the liquid molar density is fixed to ζ l = 55555 mol. m -3 . The phases viscosities are fixed to µ g = 18.51 • 10 -6 Pa. s -1 and µ l = 10 -3 Pa. s -1 . The gas fugacities are given by Daltons law for an ideal mixture of perfect gas f g i = C g i P g , i = a, w. The fugacity of the air component in the liquid phase is given by Henry's law f l a = C l a H a (T ) with the temperature dependent Henry constant H a (T ) = Ha 1 + (Ha 2 -Ha 1 ) T -T 1 T 2 -T 1 where Ha 1 = 6 • 10 9 Pa, Ha 2 = 10 10 Pa, T 1 = 293 K and T 2 = 353 K. For the water component in the liquid phase, the fugacity is taken from [START_REF] Schmidt | Properties of water and steam in S.I. units[END_REF] f l w = C l w P sat (T ) exp

k l r (S l ) =        0 if S l < S l r , 1 if S l > 1 -S g r , √ Sl ( 1 -(1 -( Sl ) 1 m ) m ) 2 if S l r ≤ S l ≤ 1 -S g
k g r (S g ) =        0 if S g < S g r , 1 if S g > 1 -S l r , √ 1 -Sl ( 1 -( Sl ) 1 m ) 2m if S g r ≤ S g ≤ 1 -S l r , P c (S l ) = Pr(( Sl ) -1 m -1) 1 n if 0 < Sl ≤ 1, with Sl = S l -S l
( - P sat (T ) -P l ζ l RT
) ,

where P sat (T ) is the vapour pressure of the pure water given by the Rankine formula

P sat (T ) = 1.013 • 10 5 exp ( 13.7 - 5120 T
) .

The liquid molar enthalpy h l and the gas molar enthalpy of the water component h g w are taken from [START_REF] Schmidt | Properties of water and steam in S.I. units[END_REF]. The gas molar enthalpy of the air component is given by h g a (T ) = c g p,a m a T where c g p,a = 1000 J. K -1 . kg -1 is the specific heat capacity of pure air and m a = 29 • 10 -3 kg. mol -1 is the air molar mass. The gas molar enthalpy is then defined by (2.1). The bulk rock thermal conductivity is fixed to λ = 10 W. m -1 . K -1 and the rock energy per unit volume is given by E r (T ) = 2 • 10 6 T in J. m -3 with T in K.

The initial and top boundary conditions in the porous medium are defined by a liquid phase S l = 1 with pure water C l w = 1, C l a = 0, a temperature T 0 pm fixed either to 303 or 333 K and a hydrostatic pressure with 38 atm at the top boundary. The lateral boundaries are considered thermally isolated (no Fourier flux) and impervious (no Darcy flux). The soil-atmosphere evaporation-outflow boundary condition at the interface Γ is set with no radiation nor precipitation recharge and the outflow liquid flux is vanishing during the overall simulation in the following test cases.

At the output boundary Γ out = {l} × (0, h f f ) of the free-flow domain, the pressure P out = 10 5 Pa is the atmospheric pressure. The velocity at the input boundary Γ in = {0} × (0, h f f ) is defined by the uncoupled turbulent velocity profile u 0 (y) = ( u 0 (y) 0

) , computed from the Prandtl algebraic turbulent model (see [START_REF] Masson | Coupling compositional liquid gas Darcy and free gas flows at porous and free-flow domains interface[END_REF]), with an average velocity

u in = 1 h ff ∫ h ff 0 u 0 (y)dy = 0.5 m.s -1 ,
and such that u 0 (h f f ) = u 0 (0) = 0. The temperature at the input boundary Γ in is fixed to T in = 303 K and the input molar fractions to

C in = ( C w,in 1 -C w,in )
with the relative humidity

P out C w,in P sat (T in ) = 0.5.
An homogeneous Dirichlet boundary condition for the velocity and homogeneous Neumann boundary conditions for the molar and energy transport are used at the boundary (0, l) × {0}. The turbulent viscosity µ t used in the free-flow domain to define the RANS stress tensor is given by the Prandtl algebraic turbulent model as in [START_REF] Masson | Coupling compositional liquid gas Darcy and free gas flows at porous and free-flow domains interface[END_REF] and computed once and for all from the uncoupled solution in the free-flow model. The turbulent diffusivity

D t = D g + µ t -µ g ρ g S c (6.1)
is deduced using the gas Fickian diffusion D g = 2 • 10 -5 m 2 . s -1 and the Schmidt number S c = 1. The turbulent thermal conductivity is similarly defined by λ t = λ g + c g p,a (µ t -µ g ) with the gas thermal con- ductivity λ g = 0.026 W. m -1 . K -1 .

The simulation is run over the time interval [0,t f ] with t f = 200 years, using an adaptive time step starting with an initial time stepping of 1 s and a maximum time step of 10 years. The Cartesian mesh is uniform in the x direction with N x = 100 edges and refined exponentially in the vertical y direction on both sides of the interface Γ to account for the turbulent boundary layer and for the high gradient of the liquid pressure. More precisely, the porous medium mesh is defined by N e > 0, N y > N e , r > 1, ∆ y r > 0, such that h f f + ∆ y r (r N e -1) < h pm . Numbering the y-edges (y i , y i+1 ), i = 1, • • • , N y + 1 from bottom to top, we set

y i =    h f f + ∆ y r (r i-1 -1), 1 ≤ i ≤ N e + 1,
y N e +1 + (i -N e -1) h pm -y N e +1 N y -N e , N e + 2 ≤ i ≤ N y + 1.
The numerical performances of the Darcy flow model coupled with the soil-atmosphere boundary condition are assessed on the following meshes

( cu 0 -D t ∇c ) = 0 on Ω f f , c = 1 on Γ , c = 0 on Γ in , ∇c • n = 0 on Γ out , ∇c • n = 0 on (0, l) × {0}, (6.2) and div ( ζ g (P out , T in ) ∂ h g w ∂ T (P out , T in )T u 0 -λ t ∇T ) = 0 on Ω f f , T = 1 on Γ , T = 0 on Γ in , ∇T • n = 0 on Γ out , ∇T • n = 0 on (0, l) × {0}. (6.3) Then, we set H m (x) = ζ g (P out , T in )D g ∇c • n f f | Γ , H T (x) = λ g ∇T • n f f | Γ , (6.4)
with n f f the normal at Γ oriented outward to the free-flow domain.

Tables 6.1 and 6.2 compare, respectively for T 0 pm = 303 K and T 0 pm = 333 K and for the three meshes, the numerical efficiency of the Newton-min non-linear solvers with their different improvements introduced in Section 4.1. Each table contains the number of successful time steps, the number of time step chops (i.e. the number of Newton-min convergence failures), the average number of Newton iterations per time step and the CPU time (in seconds on 2.9 GHz Intel Core i5 processor and 8Go RAM). It is shown that the basic Newton-min algorithm fails to converge in these test cases and that the use of the equilibrium phase molar fractions for the phase appearance criterion is necessary to obtain the convergence for the finest mesh at T 0 pm = 303 K. On the other hand, imposing the thermodynamic equilibrium at each Newton iterate improves only slightly the convergence.

The solutions of the liquid gas Darcy flow coupled either with the soil-atmosphere boundary condition or with the full-dimensional gas free flow are compared using the finest mesh with N y = 90. Figures 6.4 and 6.5 exhibit, respectively for T 0 pm = 303 K and T 0 pm = 333 K, the evolution in time of the mean relative humidity, the mean temperature and the mean molar flow rate of the water component at the interface Γ for both models. The two stages, typical of drying processes, are clearly identified. The first stage corresponds to a high liquid evaporation rate combined with a relative humidity at the interface close to one. This stage is mainly governed by the free turbulent flow as long as the interface is water saturated. The second stage, triggered by the desaturation of the porous medium, corresponds to the drop of both the evaporation rate and the relative humidity reaching their stationary state say at time 200 years. The cooling effect of the liquid evaporation at the interface is also clearly observed in the temperature plot of Figure 6.4. Figures 6.4 and 6.5 show that the soil-atmosphere boundary condition combined with the convective molar and energy transfer coefficients (6.4) provides a very good approximation of the coupled nonisothermal liquid gas Darcy and full-dimensional gas free flow model. The mismatch is larger for T 0 pm = 333 K than for T 0 pm = 303 K on the evaporation rate due to larger variations in time of the convective molar and energy transfer coefficients not captured by H m and H T . Nevertheless, the temperature and relative humidity at the interface remains very well approximated in both cases.

1D and 2D geothermal test cases

In these simulations, the porous medium is homogeneous of porosity ϕ (x) = 0.35 and of isotropic permeability Λ Λ Λ (x) = K × I with K = 1 D. The relative permeabilities are defined by k α r (S α ) = (S α ) 2 for each phase α ∈ P. The capillary pressure function is given by the Corey law with b = 2 • 10 5 Pa and S 1 = 0.99. It is regularized for S g ∈ (S 1 ; 1] to allow for the disappearance of the liquid phase (see Figure 7.1). Since there is no entry capillary pressure (in the sense that P c (0) = 0), the complementarity constraint min(P c (S g ), q l,atm ) = 0 from (3.15) is equivalent to min(S g , q l,atm ) = 0. It results that the gas saturation can be used in the following test cases as primary unknown at the interface Γ rather than the capillary pressure P c (refer to Table 4.1). Different choices, including P c and a variable switch between P c and S g are compared in Paragraph 7.1.2.

P c (S g ) = { -b ln(1 -S g ) if 0 ≤ S g ≤ S 1 -b ln(1 -S 1 ) + b 1-S 1 (S g -S 1 ) if S 1 < S g ≤ 1
The liquid and gas phases are a mixture of two components, the water denoted by w and the air denoted by a. The gas thermodynamic laws are defined by the perfect gas molar density ζ g = P g RT with R = 8.314 J. K -1 . mol -1 and the viscosity µ g = (0.361T -10.2)•10 -7 in Pa. s. The liquid molar enthalpy h l and the gas molar enthalpies of each component h g a , h g w are taken from [START_REF] Schmidt | Properties of water and steam in S.I. units[END_REF]. The gas molar enthalpy is then defined by (2.1). The liquid molar density and viscosity are also from [START_REF] Schmidt | Properties of water and steam in S.I. units[END_REF] and defined by with the salinity fixed to C s = 35 • 10 -3 kg. kg -1 . The mass density is defined by ρ α = ∑ i∈C C α i m i ζ α with the molar masses m a = 0.029 and m w = 0.018 kg. mol -1 . The vapour pressure P sat (T ) is given by the Clausius-Clapeyron equation

ζ l = (780.83795 + 1.62692 T -3.06354 • 10 -3 T 2 )(1 + 0.651 C s ) 0.018 , ( 7
P sat (T ) = 100 exp ( 46.784 - 6435 T -3.868 log(T ) ) ,
and the Henry constant of the air component is set to H a = 10 8 Pa. The molar internal energy of each phase is considered to be equal to its enthalpy. Finally, the fugacities are given by     

f g i = C g i P g , i = a, w, f l a = C l a H a , f l w = C l w P sat (T ) exp ( -P sat (T )-P l 1000RT /0.018
) .

The thermal conductivity is fixed to λ = 3 W. m -1 . K -1 and the rock energy per unit volume is given by E r (T ) = 2 • 10 6 T in J. m -3 with T in K.

1D geothermal test cases

The aim of the 1D test cases is to study the soil-atmosphere evaporation-outflow boundary condition introduced in Section 3 and in particular the appearance and disappearance of the liquid outflow. The domain is a box of length (0 m, 11000 m) and height (-3000 m, 1000 m) meshed with 1000 cells in the vertical direction. The initial condition is defined by a liquid phase S l = 1 composed of pure water C l w = 1, C l a = 0 at hydrostatic pressure with P l = 1 atm at the top boundary and a linear temperature between 300 K at the top boundary and 550 K at the bottom boundary. The lateral boundaries of the domain are thermally isolated (no Fourier flux) and impervious (no Darcy flux). In addition to the fixed temperature 550 K, we impose at the bottom boundary an input molar flow rate q w,in < 0 composed of pure liquid water.

The soil-atmosphere evaporation-outflow boundary condition is imposed at the top boundary. The short and long wave radiation coming from the atmosphere and reaching the soil surface is fixed to (1a)R s + R a = 340 W. m -2 and the soil emissivity to ε = 0.97. The convective molar and energy transfer coefficients are set to H m = 0.69 mol. m -2 . s -1 and H T = 29 × H m = 20 W. m -2 . K -1 . The far field atmospheric conditions are set to C g,atm a,∞ = 0.99, C g,atm w,∞ = 10 -2 , T atm ∞ = 300 K and P atm = 1 atm, corresponding roughly to a relative humidity of 0.5. The precipitation recharge is not considered in these 1D test cases.

The simulations are run with t f = 1200 years using an adaptive time stepping starting with an initial time step of 6 days and with a maximum time step of 700 days.

1D geothermal test case with appearance and disappearance of the outflow

In this test case, as exhibited in Figure 7.2, the time-dependent input molar flow rate q w,in (t) =    -2.9 • 10 -2 for 0 < t ≤ 300 years, 0 for 300 < t ≤ 900 years, -1.45 • 10 -2 for 900 < t ≤ t f , in mol. m -2 . s -1 is imposed at the bottom boundary in order to test the appearance and disappearance of the liquid outflow. Table 7.1 shows the numerical behaviour of the simulation with the different versions of the Newtonmin algorithm described in Section 4.1, including the number of successful time steps N t f , the number of time step chops N chops , the average number of Newton-min iterations per successful time step N Newton and the CPU time (in seconds on 2.9 GHz Intel Core i5 processor and 8Go RAM). It shows that the numerical performances are drastically increased when enforcing the complementarity constraints of type min(X 1 , X 2 ) = 0 to hold exactly at each Newton iteration. The non-linear appearance criterion slightly increases the number of time step chops compared with the linear appearance criterion. Enforcing the thermodynamic equilibrium at each Newton iterate reduces the number of time step chops but deteriorates the number of Newton iterations so it remains worse than the convergence obtained with the linear appearance criterion. Figures 7.3 and 7.4 show the time evolution of the phase pressures, temperature and gas saturation both at the top cell and at the top boundary, and the relative humidity and the phase molar flow rates at the interface. Figure 7.5 exhibits the solutions obtained at the end of each period at times t = 300, 900, 1200 years. The simulation is governed by the time-dependent input molar flow rate exhibited in Figure 7.2 and by the boundary condition at the top of the geothermal column. Typically, the hot liquid front rises by buoyancy and viscous forces to the top of the reservoir. As soon as the temperature is larger than the saturated vapour temperature, a bubble of vapour grows and rises by gravity. The connection with the top of the reservoir crucially depends on the boundary condition imposed at the top. In our case, as expected, the soil-atmosphere evaporation-outflow boundary condition automatically switches from a liquid-outflow to a vanishing liquid phase boundary condition when the input bottom molar flow rate vanishes at time t = 300 years. It switches back to a liquid outflow at time t = 900 years when the input molar flow rate becomes again strictly negative. Since there is no entry capillary pressure in this test case, the liquid outflow complementarity constraint (3.15) is equivalent to min ( S g , q l,atm ) = 0 which explains why, in Figure 7.3, the gas saturation at the interface is vanishing as soon as the liquid molar flow rate q l,atm is strictly positive (see Figure 7.4). This can also be observed in Figure 7.5 in the gas saturation blue and red plots. The solutions exhibited respectively at times t = 300 year and at final time correspond to the stationary states obtained with their respective input molar flow rate at the bottom boundary. During the time interval 300 years < t ≤ 900 years, since there is no input flow rate at the bottom boundary, the liquid outflow and the gas molar flow rate vanish rapidly. The solutions at the end of this time interval, exhibited in green in Figure 7.5, show that the liquid phase drops, air penetrates at the top of the geothermal column and the water vapour rises by gravity. The liquid phase vanishes from the top of the geothermal column for say z ∈ (978 m, 1000 m). Note also that the stationary linear conductive solution for the temperature is far from being reached after 600 years at time t = 900 years. From 630 m to say 845 m, a step in the temperature curve can be noticed. This step corresponds to the domain where the fluid is diphasic and the air molar fraction in the gas phase is null. Neglecting the Kelvin correction and the dissolution of air in the liquid phase, the equality of the liquid and gas water fugacities imposes that f g w = C g w P g ∼ P sat (T ). Since f g w = C g w P g is roughly equal to 1 atm, the temperature is also roughly constant. Tables 7.2 and 7.3 exhibit a significant gain in terms of non-linear convergence obtained with the variable switch parametrization compared with the scaled capillary pressure.

The plots of Figures 7.7 and 7.8 highlight that the liquid outflow appears at the first time step and remains throughout the simulation. This is due to a high input molar flux at the bottom boundary which saturates the atmosphere at the top boundary. Figure 7.8 exhibits that the evaporation rate increases with the temperature at the top boundary while the liquid outflow decreases. The stationary state is reached after say 200 years for this test case and is very similar to the one obtained at the end of the first period of the previous test case. This is expected since the outflow regime is reached at the first time step of both simulations and since both capillary pressure curves are quite similar.

2D Bouillante geothermal test cases

The two dimensional test case illustrated in Figure 7.9 represents a simplified 2D cut of the Bouillante geothermal reservoir. It is run with two different upper boundary conditions to compare the solutions obtained with the evaporation-outflow boundary condition introduced in Section 3 and with Dirichlet boundary conditions. The initial and left side conditions are defined by a pure water liquid phase (S l = 1, C l w = 1, C l a = 0) at hydrostatic pressure and by a linear temperature between the fixed top and bottom temperatures. The bottom boundary is impervious (no Darcy flux) with a fixed temperature of 400 K except The simulations are run over the time interval [0,t f ], t f = 1000 years, with an adaptive time stepping starting with an initial time step of 6 days in the Dirichlet case and of 1 day with the evaporation-outflow boundary condition. The maximum time step is fixed to 700 days in both cases.

Two Voronoi meshes, a coarse and a fine one, satisfying the admissibility condition of TPFA schemes at both inner and boundary faces have been generated. The coarse mesh contains approximatively 1500 cells (about 1700 degrees of freedom) and is refined at the neighbourhood of the top boundary with a volume ratio of 22 between the smallest and the largest cells of the mesh. The fine mesh contains approximatively 3500 cells (around 4000 degrees of freedom) and the refinement at the interface is characterized by a volume ratio of 115 between the smallest and the largest cells of the mesh. The coarse mesh is exhibited in Figure 7.10 as well as a zoom of the top right zone for both the coarse and the fine meshes.

2D geothermal test case with Dirichlet top boundary conditions

In this test case, the upper boundary is composed of three parts corresponding to the seabed (z ≤ 0 m and 0 ≤ x ≤ 5000 m) described above, a sunny plain zone (0 < z ≤ 500 m and 5000 m < x ≤ 8450 m) and a rainy mountain zone (z > 500 m and 8450 m < x ≤ 11000 m). The sunny plain zone is defined with the same parameters than the far field atmospheric conditions used in Subsection 7.2.2, which means that the relative humidity is fixed to 0.5, the temperature to 300 K and the gas pressure to P g = 1 atm from which we deduce that only the gas phase with the water and air molar fractions of about C g a ≃ 0.99, C g w ≃ 10 -2 . The rainy mountain zone is characterized by a two-phase flow at thermodynamic equilibrium which is fitted in such a way that the liquid flux entering the domain is similar to the one obtained in Subsection 7.2.2 with the evaporation-outflow top boundary condition including the precipitation recharge. Then, the Dirichlet boundary condition for z > 500 m (which corresponds to x > 8450 m) is defined by a fixed temperature, gas pressure and relative humidity corresponding to the following physical values S g ≃ 0.72, S l ≃ 0.28,

P g = 1 atm, P l ≃ -153671 Pa, C g a ≃ 0.97, C g w ≃ 0.03, C l a ≃ 10 -3 , C l w ≃ 0.999, T = 300 K.
The efficiency of the T-PSC formulation is compared to an alternative persistent variable formulation which is a basic extension to the non-isothermal case of the formulation introduced in [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] in the isothermal case. This formulation, denoted by T-PSF, is based on the phase pressures, temperature, phase saturations and component fugacities f = ( f i ) i∈C as set of principal unknowns. The T-PSF formulation defines the molar fractions C α as the function C α ( f , P α , T ) of the component fugacities, phase pressure and temperature obtained by inversion of the non-linear system f α (P α , T,C α ) = f . This allows extending the phase molar fractions of an absent phase by those at equilibrium with the present phase and to express the thermodynamic equilibrium as the following complementarity constraints, in a similar way as in the T-PSC formulation

S α ≥ 0, 1 -∑ i∈C C α i (P α , T, f ) ≥ 0, S α (1 -∑ i∈C C α i (P α , T, f )) = 0, α ∈ P.
Let us refer to [START_REF] Lauser | A new approach for phase transitions in miscible multi-phase flow in porous media[END_REF] for a more detailed description. Two versions of Newton updates have been implemented for the T-PSF formulation which can be viewed as variants of the Newton-min algorithm. In the first version, (called "T-PSF with projection of the saturations"), the phase disappearance is governed by the negative sign of the updated phase saturation for which we enforce S α ≥ 0 and ∑ αP S α = 1 at each Newton iterate. The negative sign of the second constraint 1 -∑ i∈C C α i (P α , T, f ) governs the phase appearance but we do not enforce 1 -∑ i∈C C α i (P α , T, f ) to vanish at each Newton iterate even if the phase is present since it would involve the solution of an additional non-linear system. Note that, for the T-PSC formulation, this equality can be easily enforced if the phase is present since the phase molar fractions are principal unknowns. This is a major advantage of the T-PSC formulation compared with the T-PSF formulation. In the second version (called "T-PSF with non-linear update of the fugacities"), the constraint ∑ i∈C C l i (P l , T, f ) = 1 is additionally enforced at each Newton iterate as soon as the liquid phase is present, taking advantage of the particular expressions of the fugacities, in the following way. After linear update of the unknowns (in particular of the gas pressure P g ν , the temperature T ν and the air fugacity f a,ν ), we set C l a,ν = min(1, max(0,

f a,ν H a )), C l w,ν = 1 -C l a,ν and f w,ν = f l w (P l ν , T ν ,C l ν )
. This value of the water fugacity is taken as Newton update and also used to test the appearance of the gas phase.

The numerical performances of both implementations of the T-PSF formulation are exhibited in Table 7.4 whereas Table 7.5 summarizes the convergence behaviour of the Newton-min algorithms with the T-PSC formulation. The comparison between Tables 7.4 and 7.5 shows that the Newton-min algorithm with projection of the T-PSC formulation outperforms on the fine mesh the Newton-min algorithm with projection of the saturations of the T-PSF formulation. It is also a much simpler algorithm since it does not involve the inversion of the systems f α (P α , T,C α ) = f which will be costly in combination with an Equation of State thermodynamic system. The non-linear update of the fugacities in the T-PSF formulation appears necessary to obtain a similar efficiency of the Newton-min algorithm on the fine mesh compared with the Newton-min algorithms of the T-PSC formulation, excluding the basic Newton algorithm which should clearly not be used. Still, the number of time step chops remains significantly larger on the fine mesh for the T-PSF formulation with non-linear update of the fugacities. Considering that this non-linear update of the fugacities seems difficult to generalize to more general thermodynamic systems, this motivates our choice of the T-PSC formulation rather than the T-PSF formulation in this work.

2D geothermal test case with the soil-atmosphere evaporation-outflow boundary condition

In this subsection, the Dirichlet conditions on the sunny plain and rainy mountain zones are replaced by the evaporation-outflow boundary condition developed in Section 3. The radiation, the convective molar and energy transfer coefficients and the far field atmospheric conditions are those defined in the 1D geothermal test cases in Subsection 7.1. The precipitation recharge is null on the sunny plain zone and fixed to q l,rain = -3.2 • 10 -2 mol. m -2 . s -1 on the rainy mountain zone with C l,rain w = 0.999 and C l,rain a = 10 -3 . This precipitation recharge corresponds to roughly twice the observed rainfall of 9 m in 2016. It has been doubled since the reservoir 2D cut is assumed to be along a fault plane which favours the water intrusion.

Figure 7.12 exhibits the temperature and the gas saturation in the whole domain at final time. When compared with Figure 7.11, one can notice that at final time the evaporation-outflow boundary condition It also results that the temperature drop near the shoreline does no longer appear. The gas saturation remains null below the seabed and the desaturated zone is shifted to x > 5000 m. It can also be noticed that the desaturated zone is deeper with the evaporation-outflow than with the Dirichlet boundary condition.

2D geothermal test case with a water-air-salt thermodynamic system

In this subsection, the previous test case of Subsection 7.2.2 is extended to take into account the dissolution of the salt component in the liquid phase. Since our model assumes all components to be present in both phases, the liquid and gas phases are now a mixture of three components, the water denoted by w, the air denoted by a and the salt denoted by s, setting C = {w, a, s}. The liquid molar density (7.1) and viscosity (7.2) are functions of the salinity C s in kg. kg -1 which is now related to the liquid molar fractions by Table 7.7 summarizes the convergence behaviour of the different versions of the Newton-min algorithm. Figure 7.13 exhibits the temperature, the gas saturation and the salt mass fraction in the liquid phase in the reservoir at final time. It is clear from the comparison between Figure 7.13 and Figures 7.11 -7.12 that the sea water intrusion prevents as expected the development of the convective thermal instabilities from the left side of the reservoir. This is due to the higher salinity of the sea water compared with the left side and bottom salinity. It also explains why the high temperature zone is shifted to the right in this simulation compared with the previous one. This shift is responsible for the vaporisation of the liquid water component near the top boundary which can be observed in the gas saturation in Figure 7.13. Indeed, the high temperature zone is closer to the shoreline where the pressure is lower which favours the vaporisation of the liquid phase. The plot of the salt molar fraction in the liquid phase at final time in Figure 7.13 clearly shows that the reservoir is split in 3 zones depending on the source of the water flux, the sea water zone on the left, the rain water zone on the right and the high temperature water zone in between. A high salt molar fraction in the liquid phase can also be noticed at the top boundary due to the liquid vaporization, it could induce the precipitation of the salt not taken into account in this model. Tables 7.5, 7.6 and 7.7 confirm that enforcing the complementarity constraints to hold at each Newton iterate considerably improves the convergence compared with the basic Newton-min version.

In both simulations with the evaporation-outflow boundary condition (Tables 7.6 and 7.7) the basic Newton-min algorithm fails to converge while the Newton-min algorithm with projection on the complementarity constraints exhibits a good non-linear convergence. In most of the cases, the non-linear phase appearance criterion also improves the non-linear convergence and its combination with the thermodynamic equilibrium update gives mixed results with sometimes a small reduction of the CPU time.

The non-linear convergence of the test case using Dirichlet boundary conditions is easier to achieve, hence the basic Newton-min algorithm succeeds in converging (refer to Table 7.5) but is about twice longer than the other versions of the Newton-min algorithm. In this simulation, the Newton-min algorithm with projection and thermodynamic equilibrium is the most efficient.

Conclusion

In this work, a new formulation for non-isothermal compositional gas liquid Darcy flows based on natural variables and using extended phase molar fractions has been introduced. The non-isothermal compositional model is coupled with an advanced soil-atmosphere boundary condition accounting for the vaporization of the liquid phase in the atmosphere, the convective molar and energy transfer, a liquid outflow condition as well as the precipitation recharge and the radiation. Newton-min algorithms with various improvements have been investigated to solve the non-linear systems obtained at each time step after an Euler implicit time integration. The numerical efficiency of the formulation and the soil-atmosphere evaporation-outflow boundary condition have been studied on several 1D and 2D test cases including in particular a 2D cut of the Bouillante high energy geothermal field in Guadeloupe with both air-water and air-water-salt thermodynamic systems. The importance for geothermal simulations of the top boundary condition taking into account the seabed, the sunny plain and the rainy mountain zones is enlightened by comparison with a fitted Dirichlet boundary condition. Regarding the non-linear solver efficiency, it is shown that enforcing the complementarity constraints to hold at each Newton iterate considerably improves the non-linear convergence. A more robust convergence is also obtained when using the nonlinear phase appearance criterion rather than its linearized version.

,

  and the reduced right hand side B (r) = (B (r) ν ) ν∈M ∪F Γ defined by the vectors of size #C + 1

Fig. 5 . 1 Fig. 5 . 2

 5152 Fig. 5.1 Pressure, temperature and gas saturation stationary solutions.

= h f f 20 mFig. 6 . 1 Fig. 6 . 2

 206162 Fig. 6.1 Computational domain of the coupled Darcy and full-dimensional free-flow models.

Fig. 6 . 3

 63 Fig. 6.3 Relative permeabilities of both phases k α r , α = g, l and capillary pressure P c as functions of the liquid saturation S l of the Callovo Oxfordian clay.

N y = 30

 30 with N e = 10, r = 1.58, ∆ y r = 1.43 • 10 -2 , N y = 60 with N e = 20, r = 1.28, ∆ y r = 1.02 • 10 -2 , N y = 90 with N e = 30, r = 1.19, ∆ y r = 8.40 • 10 -3 .The convective molar and energy transfer coefficients are computed from the following low frequency diagonal approximations of the Dirichlet to Neumann operators related to the uncoupled convection diffusion equations in the free-flow domain. Let us define the solutions c and T of the following linear convection diffusion equations by ζ g (P out , T in )div

Fig. 6 . 4

 64 Fig. 6.4 Mean relative humidity, mean temperature (in K) and mean molar flow rate of the water component (in mol. m -2 . s -1 ) at the interface Γ as functions of time (in years) for both models with T 0 pm = 303 K and N y = 90.

Fig. 6 . 5

 65 Fig. 6.5 Mean relative humidity, mean temperature (in K) and mean molar flow rate of the water component (in mol. m -2 . s -1 ) at the interface Γ as functions of time (in years) for both models with T 0 pm = 333 K and N y = 90.

Fig. 7 . 1

 71 Fig. 7.1 Relative permeabilities of both phases k α r , α = g, l and capillary pressure P c as functions of the liquid saturation S l .

Fig. 7 . 2 1

 721 Fig. 7.2 Input molar flow rate q w,in (t) at the bottom boundary as a function of time.

Fig. 7 . 3

 73 Fig. 7.3 Gas and liquid pressures (in MPa), temperature (in K) and gas saturation as functions of time both at the top boundary and at the top cell, obtained for the time-dependent input molar flow rate.

Fig. 7 . 4

 74 Fig. 7.4 Relative humidity and gas and liquid molar flow rates (in mol. m -2 . s -1 ) at the top boundary as functions of time, obtained for the time-dependent input molar flow rate.

Fig. 7 . 5

 75 Fig. 7.5 Gas and liquid pressures (in MPa), temperature (in K) and gas saturation as functions of depth (in m) and air molar fraction in the gas phase obtained at times t = 300 years, 900 years and 1200 years for the time-dependent input molar flow rate.

  of successful time steps N t f , of time step chops N chops , average number of Newton iterations per successful time step N newton and CPU time for the different versions of the Newton-min algorithm using the variable switch parametrization.

Fig. 7 . 7

 77 Fig. 7.7 Relative humidity and gas and liquid molar flow rates (in mol. m -2 . s -1 ) at the top boundary as functions of time with entry capillary pressure.

Fig. 7 . 8

 78 Fig. 7.8 Gas and liquid pressures (in MPa), temperature (in K) and gas saturation at the top boundary and at the top cell as functions of time with entry capillary pressure.

T 1 T 1 Fig. 7 . 9

 1179 Fig. 7.9 Illustration of the two dimensional domain and its boundary conditions.

Fig. 7 .

 7 Fig. 7.10 Coarse mesh (on the top) of the 2D cut of the Bouillante geothermal reservoir and zoom of the top right zone of the coarse mesh (in the bottom left) and of the fine mesh (on the bottom right).

Fig. 7 .

 7 Fig. 7.11 Temperature (in Celsius) and gas saturation above the threshold of 10 -2 at final time (1000 years) obtained with the fine mesh and the Dirichlet top boundary conditions.

Figure 7 .

 7 Figure 7.11 exhibits the temperature and the gas saturation in the reservoir at final time. The hot liquid plume rises by buoyancy and viscous forces from the bottom injection boundary to the top of the reservoir. During the simulation, convective thermal instabilities are initially observed, then the hot liquid plume is stabilized and reaches a stationary state in between the cold water intrusion from the rainy mountain boundary on the right side of the reservoir and the sea water intrusion on the left side of the reservoir. The desaturation deepens by gravity at the right top side of the reservoir down to a stationary state at the end of the simulation. A small amount of water vapour can also be observed at the top of the hot liquid plume close to the surface due to the high temperature combined with the low pressure. Convective thermal instabilities are still observed at final time in Figure 7.11 corresponding to water intrusion from the left side boundary. Taking into account the salinity dependent liquid dynamic viscosity and mass density of the sea water (see paragraph 7.2.3) will be shown to remove these thermal instabilities induced by the left side Dirichlet boundary condition. From Figure 7.11, let us remark that the hot liquid plume goes out of the reservoir at the top boundary on both sides of the shoreline approximately in the interval 3575 m ≤ x ≤ 5550 m. Inside this interval, we can observe a temperature drop in the interval 4800 m < x < 5200 m. It is explained by the vaporisation of the liquid phase which cools down the surface neighbourhood.

  4.50/201 593/11/4.93/1025 with projection and non-linear phase appearance criterion 577/1/4.40/207 614/15/4.74/1074 Newton-min 577/1/4.41/203 615/15/4.71/1071 with projection and thermodynamic equilibrium Table 7.6 Number of successful time steps, of time step chops, average number of Newton iterations per successful time step and CPU time obtained for the different versions of the Newton-min algorithm for both meshes and with the evaporation-outflow boundary condition.

Fig. 7 .

 7 Fig. 7.12 Temperature (in Celsius) and gas saturation above the threshold of 10 -2 at final time (1000 years) obtained with the fine mesh and the evaporation-outflow boundary condition.

, 7

 7 with m s = 58.44 • 10 -3 , m w = 18 • 10 -3 , m a = 29 • 10 -3 kg. mol -1 . The air and water fugacities in both phases are still given by (7.3) and the fugacities of the salt component are defined by{ f g s = C g s P g , f l s = C l s H s ,with a very low Henry constant H s = 10 -1 Pa in order to keep the vaporization of the salt component in the gas phase negligible. The Dirichlet boundary condition at the interface between the sea and the reservoir now uses the input salinity C s = 35 • 10 -3 kg. kg -1 of the sea water. The input salinity at the left side of the reservoir as well as at the bottom boundary is fixed to the lower value C s = 20 • 10 -3 kg. kg -1 . The remaining boundary and initial conditions are unchanged compared with the previous test case, considering that the initial water in the reservoir and the precipitation recharge contain no salt. Number of successful time steps, of time step chops, average number of Newton iterations per successful time step and CPU time obtained with the different versions of the Newton-min algorithm for both meshes and the air-water-salt test case.

Fig. 7 .

 7 Fig. 7.13 Temperature (in Celsius), gas saturation above the threshold of 10 -2 and salinity of the liquid phase (in g.Kg -1 ) at final time (1000 years) obtained with the fine mesh and the air-water-salt test case.

  The boundary conditions are obtained by two-point flux approximations of the component normal fluxes w i • n, i ∈ C and of the energy normal flux w e • n. These two-point fluxes are computed between the interface on the atmosphere side and the far field atmospheric conditions at a given reference height. The far field atmospheric conditions are defined by the constant gas molar fractions C g,atm ∞

	From (3.3) and (3.2), it results that
	(w e + λ t ∇T ) • n ∼ h g w (P, T )ζ g (P, T,C)u • n,
	allowing the following approximation of w e • n
	w e • n = h g w (P, T )ζ g (P, T,C)u • n -λ t ∇T • n.	(3.5)
	3.1.2 Two-point flux approximation
		.4)

Table 6 . 1

 61 Number of successful time steps, of time step chops, average number of Newton iterations per successful time step and CPU time for the three Newton methods obtained with N y = 30, 60, 90 and T 0 pm = 303 K.

	N x × N y	100 × 30	100 × 60	100 × 90
	Basic Newton-min	×	×	×
	Newton-min with projection	157/0/3.71/147 157/1/3.97/552	×
	and non-linear phase appearance criterion	157/0/3.44/147 157/0/3.80/502 157/0/3.92/1012
	Newton-min
	with projection and	157/0/3.46/140 157/0/3.78/487	157/0/3.97/988
	thermodynamic equilibrium
	N x × N y	100 × 30	100 × 60	100 × 90
	Basic Newton-min	×	×	×
	Newton-min with projection	157/0/3.52/142 157/0/3.90/497 180/2/4.09/1360
	and non-linear phase appearance criterion	157/0/3.44/145 157/0/3.80/467	157/0/3.91/999
	Newton-min
	with projection and	157/0/3.39/138 157/0/3.73/523	157/0/3.90/972
	thermodynamic equilibrium

Table 6 . 2

 62 Number of successful time steps, of time step chops, average number of Newton iterations per successful time step and CPU time for the three Newton methods obtained with N y = 30, 60, 90 and T 0 pm = 333 K.

Table 7 . 2

 72 Number of successful time steps N t f , of time step chops N chops , average number of Newton iterations per successful time step N newton and CPU time for the different versions of the Newton-min algorithm using the scaled capillary pressure as parameter.

		N t f	N chops	N newton	CPU(s)
	Basic Newton-min	738	11	2.02	135
	Newton-min with projection	718	4	1.88	115
	and non-linear phase appearance criterion	718	4	1.88	116
	Newton-min
	with projection and	714	5	1.88	114
	thermodynamic equilibrium
		N t f	N chops	N newton	CPU(s)
	Basic Newton-min	727	4	1.95	117
	Newton-min with projection	716	2	1.91	113
	and non-linear phase appearance criterion	719	1	1.94	115
	Newton-min
	with projection and	716	2	1.91	109
	thermodynamic equilibrium

Table 7 . 3

 73 Number

Table 7 . 4

 74 Number of successful time steps, of time step chops, average number of Newton iterations per successful time step and CPU time obtained for the different versions of the T-PSF formulation for both meshes with the Dirichlet top boundary condition.

		Coarse mesh	Fine mesh
	T-PSF with projection of the saturations	561/5/4.74/220 619/26/5.23/1349
	T-PSF with non-linear update of the fugacities	561/5/3.36/163	594/13/3.70/885
		Coarse mesh	Fine mesh
	Basic Newton-min	614/24/3.86/238 835/98/4.11/1919
	Newton-min with projection	570/7/3.72/181	581/7/4.28/894
	and non-linear phase appearance criterion	570/4/3.69/177	581/5/4.25/863
	Newton-min
	with projection and	566/2/3.75/174	574/1/4.32/852
	thermodynamic equilibrium

Table 7 . 5

 75 Number of successful time steps, of time step chops, average number of Newton iterations per successful time step and CPU time obtained for the different versions of the Newton-min algorithm for both meshes with the Dirichlet top boundary condition.

Acknowledgements We would like to thank the BRGM and the Provence-Alpes-Côte d'Azur Region for the co-funding of the PhD of Laurence Beaude as well as the support of the CHARMS ANR project (ANR-16-CE06-0009).

1D geothermal test case with entry capillary pressure

In this test case, the input molar flow rate is fixed to q w,in = -2.9 • 10 -2 mol. m -2 . s -1 during the overall simulation and the capillary pressure curve includes an entry pressure P e > 0 defined by the following regularization of the Corey law (see Figure 7.6)

with P e = 10 5 Pa, b = 2 • 10 5 Pa and S 1 = 0.99. The Corey law is again regularized for S g ∈ (S 1 ; 1] to allow for the disappearance of the liquid phase. The soil-atmosphere evaporation-outflow boundary condition (3.18) should account for capillary pressures P g -P l in the interval [0, P e ]. It results that, for non zero entry pressures P e , the gas saturation cannot be used anymore as primary unknown at the top boundary. Following [START_REF] Brenner | Improving Newton's Method Performance by Parametrization: The Case of the Richards Equation[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF], let us introduce a parameter τ and two continuously differentiable non-decreasing functions S : R → [0, 1] and P c : R → R, chosen such that P c (S(τ)) = P c (τ).

Then, τ is defined as an additional unknown at the top boundary and equation S g = S g (P g -P l ) in (3.18) is replaced by S g = S(τ) and P g -P l = P c (τ).

Two choices of parametrization S(τ) and P c (τ) are compared in terms of convergence of the different versions of the Newton-min algorithm in Tables 7.2 and 7.3. The first choice uses the capillary pressure scaled by the entry pressure P e as parameter τ and is defined by

The second choice is based on a variable switch between the capillary pressure and the gas saturation which is shown in [START_REF] Brenner | Improving Newton's Method Performance by Parametrization: The Case of the Richards Equation[END_REF][START_REF] Brenner | Immiscible two-phase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media[END_REF] to improve the non-linear convergence and also allows accounting for non invertible capillary functions. It is defined by

P c (τ) =

{ P e τ if τ ∈ [0, 1), P c (τ -1) if τ ∈ [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF][START_REF]Decarbonization Wedges[END_REF].