Jérôme Droniou
email: jerome.droniou@umontpellier.fr

Mohamed Laaziri
email: mohamed.laaziri@univ-cotedazur.fr

Roland Masson
email: roland.masson@univ-cotedazur.fr

Discretisations of mixed-dimensional Thermo-Hydro-Mechanical models preserving energy estimates

Keywords: Thermo-Hydro-Mechanical (THM) model, Mixed-dimensional model, Discrete Fracture Matrix model, Contact-mechanics, Thermodynamically consistent discretisation, Finite Volume, Mixed Finite Element formulation

 based on first and second principles of thermodynamics. Two formulations of the energy equation will be considered based either on energy conservation or on the entropy balance, assuming a vanishing thermo-poro-elastic dissipation. Our focus is on space time discretisations preserving energy estimates for both types of formulations and for a general single phase fluid thermodynamical model. This is achieved by a Finite Volume discretisation of the non-isothermal flow based on coercive fluxes and a tailored discretisation of the nonconservative convective terms. It is combined with a mixed Finite Element formulation of the contact-mechanical model with face-wise constant Lagrange multipliers accounting for the surface tractions, which preserves the dissipative properties of the contact terms. The discretisations of both THM formulations are investigated and compared in terms of convergence, accuracy and robustness on 2D test cases. It includes a Discrete Fracture Matrix model with a convection dominated thermal regime, and either a weakly compressible liquid or a highly compressible gas thermodynamical model.

Introduction

Thermo-Hydro-Mechanical (THM) models in fractured/faulted porous rocks play an important role in addressing the challenges of a sustainable exploitation of subsurface resources. This is for example typically the case of deep geothermal energy production and CO 2 geological storage. THM models offer valuable insights into the complex interactions between temperature changes, fluid flow, rock deformation and fracture/fault mechanical behavior within the subsurface. They play a pivotal role in assessing potential risks, informing mitigation strategies and managing hydraulic stimulation of geothermal systems. Similarly, in CO 2 sequestration projects, THM models are instrumental in predicting and preventing issues like fault reactivation which can potentially induce CO 2 leakage or seismicity.

The THM models considered in this work initially integrate a mixed-dimensional approach, coupling a non-isothermal Poiseuille flow within the network of fractures/faults represented as co-dimension one surfaces, to the non-isothermal Darcy flow in the surrounding porous rock, known as the matrix. Let us refer to [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF][START_REF] Flauraud | Domain decomposition for an asymptotic geological fault modeling[END_REF][START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF][START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF][START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF][START_REF] Brenner | Gradient Discretization of Hybrid Dimensional Darcy Flows in Fractured Porous Media with discontinuous pressure at matrix fracture interfaces[END_REF][START_REF] Flemisch | Benchmarks for single-phase flow in fractured porous media[END_REF][START_REF] Nordbotten | Unified approach to discretization of flow in fractured porous media[END_REF][START_REF] Berre | Verification benchmarks for single-phase flow in three-dimensional fractured porous media[END_REF] and the references there-in for the derivation and for various discretisations of such mixed-dimensional models in the context of single phase Darcy flows. The second key component is the thermo-poromechanical model coupling the deformation of the rock with the non-isothermal Darcy flow within the matrix domain. Let us refer to the monograph [START_REF] Coussy | Poromechanics[END_REF] as a key reference textbook on a thermodynamically consistent derivation of such models. In this work, we will follow this framework further assuming small strains and porosity variations as well as a linear thermo-poro-elastic behavior of the porous rock. The third ingredient is related to the mechanical behavior of the fractures/faults typically based on contact-mechanics taking into account frictional contact at matrix-fracture interfaces. This type of mixeddimensional poromechanical models have been the object of many recent works both in the isothermal case [START_REF] Nejati | A finite element framework for modeling internal frictional contact in three-dimensional fractured media using unstructured tetrahedral meshes[END_REF][START_REF] Garipov | Discrete fracture model for coupled flow and geomechanics[END_REF][START_REF] Berge | Finite volume discretization for poroelastic media with fractures modeled by contact mechanics[END_REF][START_REF] Franceschini | Algebraically stabilized Lagrange multiplier method for frictional contact mechanics with hydraulically active fractures[END_REF][START_REF] Bonaldi | Gradient discretization of twophase flows coupled with mechanical deformation in fractured porous media[END_REF][START_REF] Bonaldi | Gradient discretization of two-phase poro-mechanical models with discontinuous pressures at matrix fracture interfaces[END_REF][START_REF] Bonaldi | Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces[END_REF][START_REF] Boon | Mixed-dimensional poromechanical models of fractured porous media[END_REF] and in the non-isothermal case [START_REF] Salimzadeh | A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems[END_REF][START_REF] Garipov | Discrete fracture modeling approach for simulating coupled thermo-hydro-mechanical effects in fractured reservoirs[END_REF][START_REF] Stefansson | A fully coupled numerical model of thermo-hydro-mechanical processes and fracture contact mechanics in porous media[END_REF].

Building upon existing research, our work focuses on discretisations of mixed-dimensional THM models preserving energy estimates for a general single phase fluid thermodynamical model, an aspect not thoroughly explored in previous studies. We first extend to the mixed-dimensional framework the thermodynamically consistent THM models derived in [START_REF] Coussy | Poromechanics[END_REF] based on first and second principles of thermodynamics. Two formulations of the energy equation will be considered based either on energy conservation or on the entropy balance, assuming a vanishing thermo-poro-elastic dissipation. The entropy balance formulation is frequently combined with a small Darcy velocity assumption and a linearisation of the Fourier term based on a small temperature variation assumption [START_REF] Kirkesaether Brun | Monolithic and splitting solution schemes for fully coupled quasistatic thermo-poroelasticity with nonlinear convective transport[END_REF][START_REF] Jakub | The gradient flow structures of thermo-poro-visco-elastic processes in porous media[END_REF]. It leads to an approximate entropy equation which will be investigated in this work both theoretically and numerically from the points of view of energy estimates and accuracy of the solution. One difficulty that needs to be dealt with when starting from such entropy balance equation is the ability to return back to the energy equation at the discrete level. This will be addressed in this work by preserving at the discrete level the links between the energy conservation and entropy balance equations.

To achieve these goals, the space and time discretisations of the non-isothermal flow and of the contact-mechanics must be selected carefully. Regarding the flow, our framework is based on coercive fluxes and incorporate a possible upwind approximation of the convection terms in order to deal with convection dominated regimes. The discretisation of the nonconservative convective terms in the entropy balance formulation is designed in order to preserve the link with the energy formulation. Although, this setting accounts for a large class of coercive Finite Volume schemes, we focus in the following on the Hybrid Finite Volume (HFV) discretisation [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF][START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF] in order to simplify the presentation. The time integration scheme is based on a forward Euler discretisation including a semi-implicit variation in the entropy balance approach. Regarding the discretisation of the mechanical model, it will be based, for the sake of simplicity, on a conforming Finite Element approximation of the displacement field accounting for its discontinuity at matrix fracture interfaces. The discretisation of the frictional contact along the fracture network is a key feature both for the robustness of the scheme and for the derivation of energy estimates. Different formulations such as mixed or stabilized mixed formulations [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF][START_REF] Lleras | A stabilized Lagrange multiplier method for the finite element approximation of frictional contact problems in elastostatics[END_REF], augmented Lagrangian [START_REF] Burman | The augmented Lagrangian method as a framework for stabilised methods in computational mechanics[END_REF] and Nitsche methods [START_REF] Chouly | An overview of recent results on nitsche's method for contact problems[END_REF][START_REF] Chouly | Nitsche method for contact with coulomb friction: existence results for the static and dynamic finite element formulations[END_REF][START_REF] Beaude | Mixed and Nitsche's discretizations of coulomb frictional contact-mechanics for mixed dimensional poromechanical models[END_REF] have been developed to discretise Coulomb frictional contact. In line with [START_REF] Drouet | An accurate local average contact method for nonmatching meshes[END_REF][START_REF] Franceschini | Algebraically stabilized Lagrange multiplier method for frictional contact mechanics with hydraulically active fractures[END_REF][START_REF] Bonaldi | Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces[END_REF], a mixed formulation with face-wise constant Lagrange multipliers is selected. This choice maintains at the discrete level the dissipative properties of the contact terms and the lower bound on the fracture aperture. It enables the straightforward handling of complex fracture networks, including corners, tips, and intersections, and lead to local expressions of the contact equations, making possible the use of efficient semi-smooth Newton nonlinear solvers.

The rest of this article is organised as follows. Section 2 presents the mixed-dimensional THM models with energy equations based either on an energy conservation formulation or on an approximate entropy balance formulation. The energy estimates satisfied by both models are formally derived in Section 2.3. Their discretisations are described in Section 3 including the Finite Volume discretisations of the non-isothermal flow in Section 3.2 and the mixed formulation of the contact-mechanics in Section 3.3 . The discrete energy estimates are derived in Section 3.4 for the discretisations of both models which are compared numerically in Section 4 in terms of convergence, accuracy and robustness. We first consider in Section 4.1 a manufactured solution for an incompressible fluid and different Peclet numbers on a 2D square domain without fractures. Then, a 2D Discrete Fracture Matrix (DFM) model with a six fracture network is investigated in Section 4.2 both for the case of a weakly compressible liquid and for the case of a highly compressible perfect gas.

2 Mixed-dimensional Thermo-Hydro-Mechanical models

Mixed-dimensional geometry and function spaces

In what follows, scalar fields are represented by lightface letters, vector fields by boldface letters. We let Ω Ă R d , d P t2, 3u, denote a bounded polytopal domain, partitioned into a fracture domain Γ and a matrix domain ΩzΓ. The network of fractures is defined by

Γ " ď iPI Γ i ,
where each fracture Γ i Ă Ω, i P I, is a planar and simply connected polygonal domain, which is relatively open in the hyperplane it spans. Without restriction of generality, we will assume that the fractures may only intersect at their boundaries (Figure 1), that is, for any i, j P I, i ‰ j it holds Γ i X Γ j " H, but not necessarily Γ i X Γ j " H.

The two sides of a given fracture of Γ are denoted by ˘in the matrix domain, with unit normal vectors n ˘oriented outward from the sides ˘. We denote by γ a the trace operators on the side a P t`, ´u of Γ for functions in H 1 pΩzΓq and by γ BΩ the trace operator for the same functions on BΩ. The jump operator on Γ for functions u in pH 1 pΩzΓqq d is defined by u " γ `u ´γ´u , and we denote by its normal and tangential components. The notation q m n will also be used to denote the normal jump of flux functions q m P H div pΩzΓq defined by q m n " γ ǹ q m `γń q m with γ a n the normal trace operator on the side a of Γ oriented outward to the side a P t`, ´u. The notation γ a n will also be applied to tensor fields. The tangential gradient and divergence along the fractures are respectively denoted by ∇ τ and div τ . The symmetric gradient operator ε is defined such that εpvq " 1 2 p∇v `p∇vq t q for a given vector field v P H 1 pΩzΓq d . To fix ideas, the pressure p and temperature T are assumed to be continuous at the matrixfracture interface, that is, these functions belong to H 1 pΩq. As a consequence, letting γ : H 1 pΩq Ñ L 2 pΓq be the trace operator on the fracture, the fracture pressure and temperature are given by p f " γp and T f " γT .

u n " u ¨n`a nd u τ " u ´ u n n 3
The space for the displacement field is defined by

U 0 " tv P H 1 pΩzΓq d : γ BΩ v " 0u.

Mixed-dimensional models

We consider a Thermo-Hydro-Mechanical (THM) model under the hypothesis of small perturbations for the skeleton accounting for small strain, displacement and variations of porosity [START_REF] Coussy | Poromechanics[END_REF]. Linear isotropic thermo-poro-elastic constitutive laws are considered for the skeleton assuming small variations of temperature around the reference temperature T ref .

The Darcy law is used for the fluid velocity and the Fourier law for the thermal conduction. Thermal equilibrium is assumed between the fluid and the skeleton, and the mechanical inertial term is modelled using the frozen specific average fluid-rock density m 0 .

The fluid thermodynamical properties, depending on its pressure p and temperature T , are

• ϱpp, T q: specific density,

• epp, T q: specific internal energy,

• hpp, T q " epp, T q `p ϱpp,T q : specific enthalpy, and

• ηpp, T q: dynamic viscosity.

For a fluid property Ξ " ϱ, e, h, η, we use the short notations Ξ m " Ξpp, T q in the matrix and Ξ f " Ξpp f , T f q " γΞ m along the fracture network.

The primary unknowns of the model are the fluid pressure p, the fluid temperature T and the skeleton displacement field u. They are solutions of the nonlinear system of PDEs coupling the mixed-dimensional fluid mass conservation equation, the mixed-dimensional total energy conservation equation, the skeleton momentum balance equation, and the frictional contact conditions at matrix fracture interfaces. The model is detailed below starting with the equations in the matrix followed by the equations along the fracture network. Two different formulations of the energy equation will be considered.

Matrix model

In the matrix, the model accounts for the mass and energy conservation equations coupled to the skeleton momentum balance equation. Letting t F be the final simulation time, we therefore consider

B t pϱ m ϕq `divpϱ m V m q " G m in p0, t F q ˆΩ, (1a)
T B t S s `pB t ϕ `Bt pϱ m ϕe m q `divpϱ m h m V m `qm q " H m in p0, t F q ˆΩ, (1b)
m 0 B 2 t u ´divσ " F in p0, t F q ˆΩ, (1c)
with the Darcy and Fourier laws governing respectively the matrix fluid velocity and the matrix conductive thermal flux:

V m " ´Kpϕq η m ∇p, q m " ´Λm pϕq∇T, (1d)
where K and Λ m are respectively the rock permeability and the fluid rock average thermal conductivity, both possibly depending on the matrix porosity ϕ. The energy equation (1b) has been obtained under the hypothesis of reversible mechanical deformation in the sense of zero thermo-poro-mechanical dissipation:

T B t S s `pB t ϕ `σ : B t εpuq ´Bt E s " 0. ( 2
)
Assuming a linear isotropic thermo-poro-elastic behavior of the skeleton, this gives the following constitutive laws

B t ϕ " b divB t u ´αϕ B t T `1 N B t p , (3a)
B t S s " α s K s divB t u ´αϕ B t p `Cs T ref B t T, (3b)
σ " σ e puq ´b pI ´αs K s pT ´Tref qI, (3c)

σ e puq " 2µ εpuq `λ divu I, (3d)
which derive from the following volumetric skeleton internal energy [START_REF] Coussy | Poromechanics[END_REF]:

E s " µ|εpuq| 2 `λ 2 pdivpuqq 2 `1 2 " p T ı M « p T ff `αs K s T ref divu, (4)
where the matrix M , defined below, is assumed to be definite positive:

M :" « 1 N ´αϕ ´αϕ Cs T ref ff .
In the equations above, S s is the volumetric skeleton entropy, σ the total stress tensor and σ e the effective stress tensor. The parameters µ and λ are the effective Lame coefficients, N is the Biot modulus, b the Biot coefficient, K s is the bulk modulus, α s is the volumetric skeleton thermal dilation coefficient, α ϕ is the volumetric thermal dilation coefficient related to the porosity, and C s is the skeleton volumetric heat capacity.

Considering the combination 1 T ˆ(1b) ´hm T ˆ(1a) leads to the following alternative formulation of the energy equation

B t S s `ϱm ϕ T B t e m `pϱ m T ϕB t 1 ϱ m `1 T ϱ m V m ¨∇h m `1 T divq m " H m T ´hm T G m . (5)
An approximation of ( 5) is obtained based on the classical assumptions of small Darcy velocity V m and small variations of temperature around T ref .

Then, using the approximation

ϱ m V m ¨∇h m " ϱ m V m ¨∇e m `pϱ m V m ¨∇ 1 ϱ m `Vm ¨∇p " ϱ m V m ¨∇e m `pϱ m V m ¨∇ 1 ϱ m ,
and the linearisation 1 T divq m " 1 T ref divq m , we obtain the following approximate equation:

B t S s `ϱm ϕ T B t e m `pϱ m T ϕB t 1 ϱ m `1 T ϱ m V m ¨∇e m `p T ϱ m V m ¨∇ 1 ϱ m `1 T ref divq m " H m T ´hm T G m . (6)
Note that, using T ds " de `pd 1 ϱ with s the fluid specific entropy, equation (6) becomes

B t S s `ϱm ϕB t s m `ϱm V m ¨∇s m `1 T ref divq m " H m T ´hm T G m .
This equation is a classical approximate non conservative formulation of the entropy equation [START_REF] Coussy | Poromechanics[END_REF], which motivates the terminology adopted in the following of approximate entropy equation for (6) and of entropy equation for (5).

Fracture model

The mass and energy conservation equations of the reduced fracture model are obtained by integration along the fracture width of the equi-dimensional equations taking into account the mass and energy normal flux continuity at matrix fracture interfaces. This process leads to

B t pϱ f d f q `div τ pϱ f V f q ´ ϱ m V m n " G f in p0, t F q ˆΓ, (7a)
p f B t d f `Bt pϱ f d f e f q `div τ pϱ f h f V f `qf q ´ ϱ m h m V m `qm n " H f in p0, t F q ˆΩ, (7b)
with the fluid tangential velocity and thermal conductive flux integrated along the fracture width defined by

V f " ´Cf pd f q η f ∇ τ p f , q f " ´Λf pd f q∇ τ T f , (7c)
and where

d f " d 0 ´ u n (7d)
is the fracture aperture with d 0 the aperture at contact state as illustrated in Figure 2. In (7c), C f is the fracture hydraulic conductivity typically given by the Poiseuille law

C f " pd f q 3
12 , and Λ f is the fracture thermal conductivity also possibly depending on d f . At matrix fracture interfaces, a contact Coulomb frictional model is considered defined as follows

$ ' ' ' ' ' ' & ' ' ' ' ' ' % T ``T ´" 0 on p0, t F q ˆΓ, T n ď 0, u n ď 0, u n T n " 0 on p0, t F q ˆΓ, |T τ | ď ´F T n on p0, t F q ˆΓ, pB t u τ q ¨Tτ ´F T n |B t u τ | " 0 on p0, t F q ˆΓ, (7e)
where the vectorial surface tractions and their normal and tangential components are defined by

$ ' ' ' & ' ' ' %
T a " γ a n σpuq `pf n a on p0, t F q ˆΓ, a P t`, ´u,

T n " T `¨n `on p0, t F q ˆΓ, T τ " T `´pT `¨n `qn `on p0, t F q ˆΓ.

In a similar way as in the matrix, we also consider the entropy equation defined by the

combination 1 T ˆ(7b) ´hm T ˆ(7a) leading to ϱ f d f T f B t e f `pf ϱ f T f d f B t 1 ϱ f `1 T f ϱ f V f ¨∇h f `1 T f pdivq f ´ q m n q ´1 T f ´ ϱ m h m V m n ´hf ϱ m V m n ¯" H f T f ´hf T f G f . (8)
By continuity of the pressure and temperature (and thus of the fluid specific enthalpy), ϱ m h m V m n ´hf ϱ m V m n " 0 and the last term in the left-hand side of (8) therefore vanishes. However, at the discrete level, as a result of a possible upwinding this compensation may not necessarily occur. To keep that in mind in the continuous model, and justify the discretisation chosen in Section 3, we therefore write this vanishing term in the form of a (fictitious) jump of the enthalpy between the matrix and the fracture.

ϱ m h m V m n ´hf ϱ m V m n " ÿ aPt`,´u γ a n pϱ m V m qpγ a h m ´hf q.
Therefore, (8) can be recast as

ϱ f d f T f B t e f `pf ϱ f T f d f B t 1 ϱ f `1 T f ϱ f V f ¨∇h f `1 T f pdivq f ´ q m n q ´1 T f ÿ aPt`,´u γ a n pϱ m V m qpγ a h m ´hf q " H f T f ´hf T f G f . (9)
Its approximation based on small Darcy velocity and temperature variations assumptions writes The model is closed by considering no-flow and no-energy flux boundary conditions both for the matrix and fractures. Let us also recall that the pressure and temperature are assumed to be continuous at matrix pressure interfaces. At fracture intersections, the fracture pressure p f and temperature T f are classically also continuous, and mass and thermal flux conservation is imposed.

ϱ f d f T f B t e f `pf ϱ f T f d f B t 1 ϱ f `1 T f ϱ f V f ¨∇e f `pf T f ϱ f V f ¨∇ 1 ϱ f `1 T ref pdivq f ´ q m n q ´1 T f ÿ aPt`,´u γ a n pϱ m V m q ´γa e m ´ef `pf p 1 γ a ϱ m ´1 ϱ f q ¯" H f T f ´hf T f G f . (10)
In the sequel, we use the following terminology to refer to the two mixed-dimensional THM systems:

(1)-( 2)-( 3)-( 7): "enthalpy-based " THM model, (H-model)

(1)-( 2)-( 3)-( 7) with (1b) replaced by the approximation ( 6) and (7b) replaced by the approximation (10): "entropy-based " THM model. (S-model)

Energy estimates

Let us first consider the enthalpy-based THM model (H-model). Summing the integral over Ω of (1b) `Bt u ¨(1c) and the integral over Γ of (7b), taking into account the contact and boundary conditions, the following contact persistence condition (which results from (7e))

ż Γ T n B t u n dσ " 0,
and the zero dissipation equation ( 2), we obtain formally the following energy estimate:

B t ż Ω pE s `ϱm ϕe m qdx `Bt ż Γ ϱ f d f e f dσ`ż Γ ´F T n | B t u τ |dσ " ż Ω pH m `F ¨Bt uqdx `żΓ H f dσ. (11)
Let us now consider the entropy-based THM model (S-model). By construction, the energy estimate for this model can be deduced from the previous estimate [START_REF] Brenner | Gradient Discretization of Hybrid Dimensional Darcy Flows in Fractured Porous Media with discontinuous pressure at matrix fracture interfaces[END_REF] just by adding to the left hand side of (11) the term

ż Ω ´´V m ¨∇p `p T T ref ´1qdivq m ¯dx `żΓ ´´V f ¨∇τ p f `p T f T ref ´1qpdivq f ´ q m n q ¯dσ,
which corresponds to the difference between the exact terms and their approximations.

Recalling the definitions (1d) and (7c) of V m , q m , V f and q f , this leads to the following stronger form of the energy estimate

B t ż Ω pE s `ϱm ϕe m qdx `Bt ż Γ ϱ f d f e f dσ `żΓ ´F T n | B t u τ |dσ `żΩ ´K η m ∇p ¨∇p `Λm T ref |∇T | 2 ¯dx `żΓ ´Cf η f |∇ τ p f | 2 `Λf T ref |∇ τ T f | 2 ¯dσ " ż Ω pH m `F ¨Bt uqdx `żΓ H f dσ. ( 12
)
3 Discretisation

The objective is to design a discretisation preserving the energy estimates ( 11) and ( 12) for respectively the enthalpy-based and entropy-based THM models. Both types of formulations will be compared in the numerical section. The schemes must account for both diffusive and convective dominated energy transport. This motivates the choice of a Finite Volume formulation of the mass and energy equations with possible upwinding of the mobilities. To fix ideas, we will consider in this work the mixed-dimensional Hybrid Finite Volume (HFV) discretisation introduced in [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF] accounting for the continuity of the pressure and temperature at matrix fracture interfaces. It will be combined with a mixed variational formulation of the contact-mechanics with a conforming Finite Element discretisation of the displacement field and a facewise constant approximation of the Lagrange multiplier introduced to represent the surface tractions along the fractures. This choice handles fracture networks with corners and intersections, and leads to a local formulation of the contact conditions. It also ensures the discrete persistence property of the contact term and a dissipative frictional term which are key conditions to obtain an energy estimate for the coupled THM systems.

Space and time discretisations

Let M denote the set of polytopal cells, and F the set of faces of the mesh, with internal faces gathered in F int and boundary faces in F ext . The subset F K Ă F denotes the set of faces of the cell K P M. We denote by σ " K|L the internal face shared by the two cells K, L gathered in the subset M σ " tK, Lu. The notations σ " K|¨and M σ " tKu are used for a face σ P F K X F ext . The mesh is assumed conforming to the fracture network Γ in the sense that there exists a subset

F Γ of F such that Γ " ď σPF Γ σ.
The subset of edges of a face σ is denoted by E σ and we define the set of edges of Γ by The HFV method is a Gradient Discretisation (GD) defined by the vector space of discrete unknowns

E Γ " Ť σPF Γ E σ .
X D " tv D " ppv K q KPM , pv σ q σPF , pv ζ q ζPE Γ qu " R MYF YE Γ ,
and the reconstruction operators Π m D , ∇ m D in the matrix and Π f D , ∇ f D along the fractures (see [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF] for their detailed definition). We assume that Π m D and Π f D are piecewise constant reconstructions [17, Definition 2.12] and, more specifically, pΠ m D v D q |K " v K for all K P M, and pΠ f D v D q |σ " v σ for all σ P F Γ . From the gradient reconstruction operator ∇ m D in the matrix, and given any K P M and any symmetric positive definite d ˆd matrix D K , fluxes F K,σ pD K ; ¨q : X D Ñ R (for σ P F K ) are defined such that, for all v D , w D P X D and all

K P M, ż K D K ∇ m D v D ¨∇m D w D dx " ÿ σPF K F K,σ pD K ; v D qpw K ´wσ q. ( 13
)
Likewise, from the fracture tangential gradient reconstruction operator ∇ f D , for any σ P F Γ and any symmetric positive definite pd ´1q ˆpd ´1q matrix D σ , fracture fluxes F σ,ζ pD σ ; ¨q :

X D Ñ R are defined (for ζ P E σ ) such that, for all v D , w D P X D , ż σ D σ ∇ f D v D ¨∇f D w D dσ " ÿ ζPEσ F σ,ζ pD σ ; v D qpw σ ´wζ q. ( 14
)
These fluxes are local to each cell K P M (resp. each fracture face σ P F Γ ) in the sense that they only depend on pv K , pv σ q σPF K q (resp. pv σ , pv ζ q ζPEσ q) [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF].

In order to account for convective dominated energy transport, an upwind approximation needs to be introduced according to the sign of the Darcy flux F K,σ pK K ; p D q using the following notations. Let the index "Kσ, `" denote either σ if no upwinding is used, or an upwind choice between K and L if σ " K|L P F int zF Γ , or an upwind choice between K and σ if σ " K|¨P F ext Y F Γ (note that, if σ " K|L is an internal face and not a fracture face, this upwinding must not depend on K or L; for example, ϱ Kσ,`" ϱ Lσ,`) .

Similarly, for a given fracture edge ζ, the index "σζ, `" denote either ζ if no upwinding is used, or an upwind value between σ and σ 1 if F Γ,ζ " tσ, σ 1 u and ζ is not a boundary edge, or an upwind value between σ and ζ in the other cases. The choice of the upwind value is done according to the sign of the Darcy flux F σ,ζ pC f,σ ; p D q.

For any fluid thermodynamical property Ξ " h, e, η, ϱ, and p D , T D P X D , we define Ξ D P X D component by component by setting Ξ ν " Ξpp ν , T ν q for all ν P M Y F Y E Γ .

The displacement field u is discretised using a finite dimensional subspace U D of U 0 , typically given by a Finite Element Method.

We consider a time discretisation pt n q n"0,...,N of the time interval p0, t F q with t 0 " 0 and t N " t F , and denote by ∆t n " t n ´tn´1 the time step n. If f " pf n q n"0,...,N is a family of functions or vectors of X D , the discrete time derivative of f is defined as

δ n t f " f n ´f n´1 ∆t n .
The second time derivative of the discrete displacement field is defined by

δ n t 9 u D " 2 9 u n D ´9 u n´1 D t n ´tn´2 with 9 u n D " δ n t u D .
It corresponds to the generalisation to non-uniform time steps of the standard centered discretisation of B 2 t u. To alleviate notations, whenever w D P X N `1

D

and ' " tm, f u, we set δ n t,' w D " δ n t pΠ ' D w D q. We also consider the discrete time derivative of products of functions and elements of X D , or elements of X D , by setting:

δ n t,' pf v D w D q " δ n t pf Π ' D v D Π ' D w D q
for all families of functions f " pf n q n"0,...,N and elements v D , w D P X N `1 D . The discrete porosities pϕ n D q n"0,...,N and skeleton entropies pS n s,D q n"0,...,N are families of piecewise constant functions Ω Ñ R on M such that ϕ 0 D , S 0 s,D are given (e.g., as projections of the continuous initial porosity and entropy) and, for all n " 1, . . . , N ,

δ n t ϕ D " bπ M pδ n t divu D q ´αϕ δ n t,m T D `1 N δ n t,m p D , (15)
δ n t S s,D " α s K s π M pδ n t divu D q ´αϕ δ n t,m p D `Cs T ref δ n t,m T D , (16)
where π M is the projection on piecewise constant functions on M (that is,

pπ M f q |K " 1 |K| ş K f for all K P M).
In the fracture, we define the discrete apertures as pd n f,D q n"0,...,N , with

d n f,D : Γ Ñ R given by d n f,D " π F ,Γ pd 0 ´ u n D n q. ( 17
)
Here, π F ,Γ is the L 2 -projection on facewise constant functions on Γ.

We also define, for n " 1, . . . , N , the function p G n m : Ω Ñ R which is piecewise constant on M equal on K P M to the average of G m on pt n´1 , t n q ˆK. The piecewise constant functions p

H n m and p F n on M, as well as p G n f and p H n f on F Γ are similarly defined. In the following, we use the notations

V n K,σ " F K,σ ´Kpϕ n´1 K q η n´1 K ; p n D ¯and V n σ,ζ " F σ,ζ ´Cf pd n´1 f,σ q η n´1 σ ; p n D ¯, (18)
for the matrix and fracture Darcy flux, as well as

Q n K,σ " F K,σ ´Λm pϕ n´1 K q; T n D ¯and Q n σ,ζ " F σ,ζ ´Λf pd n´1 f,σ q; T n D ¯( 19
)
for the matrix and fracture Fourier fluxes. Note that K, Λ m , C f and Λ f could also explicitly depend on x. This dependence is omitted for the sake of simplicity.

Discretisation of the mass and energy equations

Let us first consider the HFV discretisation of the mass and energy conservation equations of the enthalpy-based THM model (H-model). For all time step n " 1, . . . , N , it couples the discrete mass conservation equations

|K|δ n t pϱ K ϕ K q `ÿ σPF K ϱ n Kσ,`V n K,σ " |K| p G n m,K @K P M, (20a)
|σ|δ n t pϱ σ d f,σ q `ÿ ζPEσ ϱ n σζ,`V n σ,ζ ´ÿ KPMσ ϱ n Kσ,`V n K,σ " |σ| p G n f,σ @σ P F Γ , (20b)
ÿ KPMσ V n K,σ " 0 @σ P FzF Γ , (20c)
´ÿ σPF Γ,ζ ϱ n σζ,`V n σ,ζ " 0 @ζ P E Γ , (20d)
to the discrete total energy conservation equations

|K| ´T n K δ n t S s,K `pn K δ n t ϕ K `δn t pϱ K ϕ K e K q ÿ σPF K pϱ n Kσ,`h n Kσ,`V n K,σ `Qn K,σ q " |K| p H n m,K @K P M, (21a)
|σ| ´pn σ δ n t d f,σ `δn t pϱ σ d f,σ e σ q ¯`ÿ ζPEσ pϱ n σζ,`h n σζ,`V n σ,ζ `Qn σ,ζ q ´ÿ KPMσ pϱ n Kσ,`h n Kσ,`V n K,σ `Qn K,σ q " |σ| p H n f,σ @σ P F Γ , (21b)
ÿ KPMσ Q n K,σ " 0 @σ P FzF Γ , (21c)
´ÿ σPF Γ,ζ pϱ n σζ,`h n σζ,`V n σ,ζ `Qn σ,ζ q " 0 @ζ P E Γ , (21d)
Note that (20c), (21c) embed both the conservation of fluxes across internal faces, as well as the zero-flux conditions on boundary faces. The same remark holds for the flux conservation equations (20d), (21d) at fracture edges but keeping the full nonlinear fluxes due to the dependence on σ of the upwind density and enthalpy when F Γ,ζ contains more than two faces (which happens at internal crossing lines between three or more fractures). Similar considerations hold for (23c) and (23d) below.

Temporarily dropping the time index n for simplicity, the HFV discretisation of the approximate entropy equations ( 6)-( 10) is based on the following discretisation of ϱ m V m ¨∇w m " divpϱ m w m V m q ´wm divpϱ m V m q on a given cell K, which includes a possible upwinding: for

w D P X D , ÿ σPF K w Kσ,`ϱKσ,`VK,σ ´wK ÿ σPF K ϱ Kσ,`VK,σ " ÿ σPF K ϱ Kσ,`VK,σ pw Kσ,`´wK q.
This is a key choice which preserves the link between the non conservative entropy and the conservative energy formulations in the sense that it is designed to satisfy the following discrete version of

w m divpϱ m V m q `ϱm V m ¨∇w m " divpϱ m w m V m q: w K ÿ σPF K ϱ Kσ,`VK,σ `ÿ σPF K ϱ Kσ,`VK,σ pw Kσ,`´wK q " ÿ σPF K
w Kσ,`ϱKσ,`VK,σ .

Likewise, in the fractures, we use the following discretisations of ϱ f V f ¨∇τ w f ÿ ζPEσ ϱ σζ,`Vσ,ζ pw σζ,`´wσ q, on σ P F Γ , and

ÿ σPF Γ,ζ
´ϱσζ,`Vσ,ζ pw σζ,`´wζ q, on ζ P E Γ . In the same spirit, the terms such as ϱ m ϕB t w m in the matrix are discretised by

Π m D ϱ n´1 ϕ n´1 D δ n t,m w D to ensure that Π m D ϱ n´1 D ϕ n´1 D δ n t,m w D `Πm D w n D ˆδn t,m pϱ D ϕ D q " δ n t,m pϱ D ϕ D w D q. ( 22
)
The same discretisation is applied for the terms such as ϱ f d f B t w f in the fractures. Following this methodology, the approximate entropy equations are discretised for all time step n " 1, . . . , N by

|K|δ n t S s,K `|K| T n K ϱ n´1 K ϕ n´1 K ´δn t e K `pn K δ n t 1 ϱ K ¯`1 T ref ÿ σPF K Q n K,σ `1 T n K ÿ σPF K ϱ n Kσ,`V n K,σ ´en Kσ,`´e n K `pn K p 1 ϱ n Kσ,`´1 ϱ n K q " |K| T n K p p H n m,K ´hn K p G n m,K q @K P M, ( 23a
) |σ| T n σ ϱ n´1 σ d n´1 f,σ ´δn t e σ `pn σ δ n t 1 ϱ σ ¯`1 T ref ´ÿ ζPEσ Q n σ,ζ ´ÿ KPMσ Q n K,σ 1 T n σ ÿ ζPEσ ϱ n σζ,`V n σ,ζ ´en σζ,`´e n σ `pn σ p 1 ϱ n σζ,`´1 ϱ n σ q 1 T n σ ÿ KPMσ ϱ n Kσ,`V n K,σ ´en Kσ,`´e n σ `pn σ p 1 ϱ n Kσ,`´1 ϱ n σ q " |σ| T n σ p p G n f,σ ´hn σ p G n f,σ q @σ P F Γ , (23b)
ÿ KPMσ Q n K,σ " 0 @σ P FzF Γ , (23c)
´1 T n ζ ÿ σPF Γ,ζ ϱ n σζ,`V n σ,ζ ´en σζ,`´e n ζ `pn ζ p 1 ϱ n σζ,`´1 ϱ n ζ q ¯´1 T ref ÿ σPF Γ,ζ Q n σ,ζ " 0 @ζ P E Γ , (23d)
By construction, applying this space time discretisation to the entropy equations ( 5)-( 9) would lead to an equivalence between the entropy and energy formulations at the discrete level when combined with the discrete mass equations. It results that the following Lemma states the equivalence up to the terms which have been neglected or linearised in the approximate entropy equations.

Lemma 3.1. For all cell K P M, T n K ˆ(23a) `hn K ˆ(20a) is equivalent to (21a) upon correcting the left-hand side by adding

ˆT n K T ref ´1˙ÿ σPF K Q n K,σ ´ÿ σPF K V n K,σ pp n Kσ,`´p n K q.
For all fracture face σ P F Γ , T n σ ˆ(23b) `hn σ ˆ(20b) is equivalent to (21b) upon correcting the left-hand side by adding

ˆT n σ T ref ´1˙´ÿ ζPEσ Q n σ,ζ ´ÿ KPMσ Q n K,σ ¯´ÿ ζPEσ V n σ,ζ pp n σζ,`´p n σ q `ÿ KPMσ V n K,σ pp n Kσ,`´p n σ q.
For all fracture edge ζ P E Γ , T n ζ ˆ(23d) `hn ζ ˆ(20d) is equivalent to (21d) upon correcting the left-hand side by adding

´ˆT n ζ T ref ´1˙ÿ σPF Γ,ζ Q n σ,ζ `ÿ σPF Γ,ζ V n σ,ζ pp n σζ,`´p n ζ q.
Proof. Let us detail only the proof for the cell terms since fracture face and edge terms are similar. Gathering the terms of T n K ˆ(23a) `hn K ˆ(20a) leads to the equation:

|K| " T n K δ n t S s,K `hn K δ n t pϱ K ϕ K q `ϱn´1 K ϕ n´1 K ´δn t e K `pn K δ n t 1 ϱ K ¯ȷ `hn K ÿ σPF K ϱ n Kσ,`V n K,σ `ÿ σPF K ϱ n Kσ,`V n K,σ ´en Kσ,`´e n K `pn K p 1 ϱ n Kσ,`´1 ϱ n K q T n K T ref ÿ σPF K Q n K,σ " |K| p H n m,K .
Recalling that h " e `p{ϱ, we obtain

|K| " T n K δ n t S s,K `´e n K δ n t pϱ K ϕ K q `ϱn´1 K ϕ n´1 K δ n t e K pn K ´ϱn´1 K ϕ n´1 K δ n t 1 ϱ K `1 ϱ n K δ n t pϱ K ϕ K q ¯ȷ ÿ σPF K ϱ n Kσ,`V n K,σ ´hn Kσ,``p p n K ´pn Kσ,`q 1 ϱ n Kσ,`¯`T n K T ref ÿ σPF K Q n K,σ " |K| p H n m,K .
Invoking [START_REF] Flemisch | Benchmarks for single-phase flow in fractured porous media[END_REF], this relation reduces to

|K| ´T n K δ n t S s,K `pn K δ n t ϕ K `δn t pϱ K ϕ K e K q ¯`ÿ σPF K ´ϱn Kσ,`h n Kσ,`V n K,σ `Qn K,σ ÿ σPF K V n K,σ pp n Kσ,`´p n K q `ˆT n K T ref ´1˙ÿ σPF K Q n K,σ " |K| p H n m,K ,
which is the expected result for the cell term.

Mixed variational discretisation of the contact-mechanical model

The subspace M D Ă L 2 pΓq denotes the set of piecewise constant functions on the partition F Γ and we set M D " pM D q d . For λ D in M D , we use the decomposition λ D " pλ D,n , λ D,τ q with λ D,n " λ D ¨n`, λ D,τ " λ D ´λD,n n `. We denote by λ σ the constant value of λ D P M D on the face σ P F Γ and by λ n,σ and λ τ,σ its normal and tangential components.

We consider the mixed formulation based on facewise constant Lagrange multipliers, which seamlessly deals with fracture intersections, corners and tips and leads to local expressions of the discrete contact conditions and efficient semi-smooth Newton solvers. On the other hand, this approach requires to assume the following uniform inf-sup condition between the space U D of displacement fields and the space M D of Lagrange multipliers: there exists c ‹ independent on the mesh such that inf

µ D PM D sup v D PU D ż Γ µ D ¨ v D dσ }v D } U 0 }µ D } H ´1 2 pΓq d ě c ‹ ą 0. ( 24
)
Let us define the discrete dual cone of normal Lagrange multipliers as

Λ D " tλ D,n P M D | λ D,n ě 0 on Γu,
and the discrete dual cone of vectorial Lagrange multipliers given λ D,n P Λ D as

Λ D pλ D,n q " tµ D " pµ D,n , µ D,τ q P M D | µ D,n ě 0, |µ D,τ | ď F λ D,n on Γu.
Note that the friction coefficient F is assumed to be facewise constant on the partition F Γ .

The mixed discretisation of the quasi static contact-mechanical model reads: find pu D , λ D q P U D ˆΛD pλ D,n q such that, for all pv D , µ D q P U D ˆΛD pλ D,n q,

ż Ω ´m0 δ n t 9 u D ¨vD `σpu n D q : εpv D q ´rb Π m D p n D `αs K s pΠ m D T n D ´Tref qs divpv D q ¯dx `żΓ λ n D ¨ v D dσ `żΓ Π f D p n D v D n dσ " ż Ω p F n ¨vD dx, (25a) ż
Γ ´pµ D,n ´λn D,n q u n D n `pµ D,τ ´λn D,τ q ¨ δ n t u D τ ¯dσ ď 0, (25b)
The variational inequality in (25b) is equivalent to the contact conditions between λ σ and the face average of u D denoted by u D σ (see e.g. Lemma 4.1 of [START_REF] Bonaldi | Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces[END_REF]). This is a key property to obtain the following discrete persistence property of the normal contact term

ż Γ λ n D,n δ n t u D n dσ ě 0, (26)
leading to the dissipative property of the contact term

ż Γ λ n D ¨ δ n t u D dσ ě ż Γ F λ n D,n | δ n t u D τ |dσ ě 0. ( 27
)
We refer the interested reader to [START_REF] Bonaldi | Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces[END_REF] for the proof.

Summary of the schemes and discrete energy estimates

To summarise, the discretisations of the enthalpy-based and entropy-based formulations are as follows:

(20)-( 21)-( 25)-( 15)-( 16)-( 17) for the enthalpy-based model (H-model).

(H-scheme) ( 20)-( 23)-( 25)-( 15)-( 16)-( 17) for the entropy-based model (S-model).

(S-scheme)

Let us define the discrete skeleton internal energy E n s : Ω Ñ R by

E n s " 1 2 " Π m D p n D Π m D T n D ı M « Π m D p n D Π m D T n D ff `αs K s T ref divu n D `µ|εpu n D q| 2 `λ 2 pdivu n D q 2 .
Then, the following propositions state the discrete energy estimates for both formulations of the energy equation.

Proposition 3.2. The scheme (H-scheme) based on the conservative enthalpy formulation of the energy equation satisfies, for all n " 1, . . . , N , the following energy estimate

ż Ω m 0 2 δ n t | 9 u D | 2 dx `żΩ δ n t E s `żΩ δ n t,m pϱ D ϕ D e D q `żΓ δ n t,f pϱ D d f,D e D q `żΓ F λ n D,n | δ n t u D τ |dσ ď ż Ω ´p H n m `p F n ¨δn t u D ¯dx `żΓ p H n f dσ. (28)
Proof. Adding (25a) with v D " δ n t u D " 9 u n D to the sums over M of (21a), over F Γ of (21b) and over E Γ of (21d), and taking into account the flux conservativity and the homogeneous Neumann boundary conditions (20c)-(21c), together with the definitions of the discrete fracture aperture [START_REF] Droniou | The Gradient Discretisation Method[END_REF], porosity [START_REF] Chouly | Nitsche method for contact with coulomb friction: existence results for the static and dynamic finite element formulations[END_REF] and skeleton entropy [START_REF] Coussy | Poromechanics[END_REF], we obtain the following equality

ż Ω m 0 δ n t 9 u D ¨9 u n D dx `żΩ δ n t,m pϱ D ϕ D e D qdx `żΓ δ n t,f pϱ D d f,D e D qdσ `żΩ ´"Π m D p n D Π m D T n D ı M δ n t,m « p D T D ff dx `σe pu n D q : δ n t εpu D qdx `αs K s T ref divδ n t u D ¯dx `żΓ λ n D ¨ δ n t u D dσ " ż Ω ´p H n m `p F n ¨δn t u D ¯dx `żΓ p H n f dσ.
Then, using the diffusive properties of the implicit Euler time stepping (that is e.g., δ n t 9

u D 9 u n D ě 1 2 | 9 u n D | 2 ´1 2 | 9 u n´1 D | 2
) and the dissipative property of the discrete contact term [START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF], we obtain the energy estimate [START_REF] Lleras | A stabilized Lagrange multiplier method for the finite element approximation of frictional contact problems in elastostatics[END_REF].

Proposition 3.3. The scheme (S-scheme) based on the non-conservative approximate entropy equation satisfies, for all n " 1, . . . , N , the following energy estimate

ż Ω m 0 2 δ n t | 9 u D | 2 dx `żΩ δ n t E s `żΩ δ n t,m pϱ D ϕ D e D q `żΓ δ n t,f pϱ D d f,D e D q `żΩ ´Kpϕ n´1 D q Π m D η n´1 D ∇ m D p n D ¨∇m D p n D `Λm pϕ n´1 D q T ref |∇ m D T n D | 2 ¯dx `żΓ ´Cf pd n´1 f,D q Π f D η n´1 D |∇ f D p n D | 2 `Λf pd n´1 f,D q T ref |∇ f D T n D | 2 ¯dσ `żΓ F λ n D,n | δ n t u D τ |dσ ď ż Ω ´p H n m `p F n ¨δn t u D ¯dx `żΓ p H n f dσ. ( 29
)
Proof. The estimate ( 29) is a consequence of the proof of Proposition 3.2 combined with Lemma 3.1 and the definitions of the Darcy [START_REF] Droniou | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF] and Fourier [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF] fluxes from the coercive matrix ( 13) and fracture ( 14) fluxes.

Numerical experiments

Two test cases are considered in this Section. The first is based on a manufactured solution for both the entropy and enthalpy-based models on a square domain with no fracture. The objective is to assess and compare the convergence of the schemes for both models and both for centred and upwind approximations of the thermal convection. Two values of the permeability are tested to induce either convection dominated or equilibrated diffusion and convection thermal regimes.

The second example considers a Discrete Fracture Matrix (DFM) model introduced in [START_REF] Berge | Finite volume discretization for poroelastic media with fractures modeled by contact mechanics[END_REF] including a six fracture network. The setting of the simulation follows the test case proposed in [START_REF] Stefansson | A fully coupled numerical model of thermo-hydro-mechanical processes and fracture contact mechanics in porous media[END_REF] and the physical and numerical behavior of the discrete models are investigated both for the case of a slightly compressible liquid and for the case of a perfect gas.

In all these numerical experiments, the HFV discretisations of the non-isothermal flow are combined with the P 2 conforming Finite Element discretisation of the displacement field. This choice ensures the inf-sup condition [START_REF] Garipov | Discrete fracture model for coupled flow and geomechanics[END_REF] for the P 2 ´P0 mixed formulation of the contact-mechanics [START_REF] Garipov | Discrete fracture modeling approach for simulating coupled thermo-hydro-mechanical effects in fractured reservoirs[END_REF]. It also satisfies the inf-sup condition between the displacement and pressure discrete spaces which prevents potential oscillations of the pressure field at short times in the undrained regime.

The coupled nonlinear system is solved at each time step using a fixed-point method on the function:

g p,T : pp, T q ÝÑ Contact Mechanics Solve u h ÝÑ Flow Solve pp, T q,
accelerated by a Newton-Krylov algorithm -which have proven to be efficient within this context in the isothermal case [START_REF] Bonaldi | Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces[END_REF]. The stopping criteria is set to 10 ´10 on the relative residual. We refer to [START_REF] Kirkesaether Brun | Monolithic and splitting solution schemes for fully coupled quasistatic thermo-poroelasticity with nonlinear convective transport[END_REF] for an investigation of various fixed-point algorithms for THM models.

At each iteration of the Newton-Krylov algorithm, the p, T sub-system is solved using a Newton-Raphson algorithm and the contact-mechanics is solved using a semi-smooth Newton method. In both cases, the stopping criteria is defined by a relative norm of the residual set to 10 ´10 or a scaled maximum Newton increment of 10 ´10 . At each Newton iteration, the linear system is solved using the sparse direct solver SuperLU version 4.3 both for the non-isothermal flow and for the contact-mechanics.

Manufactured solution without fractures

We investigate in this section the numerical convergence of the schemes (H-scheme) and (S-scheme) for respectively the (H-model) and (S-model) models using the following analytical solution upx, tq " 10 ´1e ´t ˜x2 y 2 ´x2 y 2 ¸, ppx, tq " e ´t sin pxq sin pyq , T px, tq " e ´t p2 ´cos pxq cos pyqq , on the domain Ω " p0, 1q 2 m 2 and time interval p0, t F q with t F " 1 s. The fluid is assumed incompressible and its specific internal energy is defined by epT q " T . Two values K " 100 I m 2 and K " I m 2 (with I the identity matrix) of the homogeneous isotropic permeability are considered. They are chosen such that the resulting Darcy velocity V m corresponds to a thermal convection dominated regime in the first case, and to similar orders of magnitude for thermal convection and diffusion in the second case. Dirichlet boundary conditions are imposed for p, T and u on p0, t F q ˆBΩ and the source terms G m , H m and F are computed from the analytical solution based on the data set defined in Table 1. The domain Ω is discretised using the first family of triangular meshes from [START_REF] Droniou | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF] as illustrated in Figure 3. Each mesh indexed by m P t1, 2, 3, 4u includes #M " 56 ˆ4m´1 triangles. We consider a uniform time stepping of p0, t F q with time step ∆t " 2.10 ´5 s chosen small enough to reduce the error due to the time discretisation and focus on the convergence in space.

The convergence of the L 2 space time errors for p, T, u and their gradients are exhibited for both schemes and both permeabilities in Figures 45as functions of the mesh step. Both the upwind and centred schemes are considered for the thermal convection. Note that the centred schemes are not presented for K " 100 I since they fail to provide a solution due to their instability in the convection dominated regime. From Figures 45, we first observe that the discretisations of both models provide a very similar convergence behavior for a given choice of the permeability and of the centred or upwind approximation of the thermal convection terms. Regarding the displacement field, second and first order convergence rates are observed in all cases for respectively u and ∇u. This is in accordance with the cellwise constant reconstruction Π m D of the pressure and temperature in the displacement field variational formulation (25a). The convergence rates for p and ∇p are respectively roughly 2 and 1 in all cases as could be expected. On the other hand, the converge rates for T and ∇T depend on the approximation of the convection term and on the convectiondiffusion regime. For equilibrated convection and diffusion, the centred scheme provides a higher convergence rate of order roughly 2 for T than the upwind scheme of order between 1 and 2. An order 1 is observed on ∇T for both the centred and upwind schemes. In the convection dominated regime, the upwind scheme exhibits a convergence rate slightly better than 1 for T , and an order between 0.5 and 1 for ∇T slightly better for the entropy-based than for the enthalpy-based model. Figure 3. Square domain Ω with its triangular mesh m " 2 using 56 ˆ4 cells.

2D Discrete Fracture Matrix (DFM) model

We consider the DFM model exhibited in Figure 6 and introduced in [START_REF] Berge | Finite volume discretization for poroelastic media with fractures modeled by contact mechanics[END_REF]. It is defined on the domain Ω " p0, 1q ˆp0, 2q m 2 and includes a network Γ of six fractures, among which fracture 1 made up of two sub-fractures forming a corner, nearly intersecting fractures and the non-immersed fracture 5 with one tip lying on the left boundary. The setting of the simulation follows the test case proposed in [START_REF] Stefansson | A fully coupled numerical model of thermo-hydro-mechanical processes and fracture contact mechanics in porous media[END_REF] and consists in three stages on the time intervals I 1 " p0, t p1q q, I 2 " rt p1q , t p2q q and I 3 " rt p2q , t F s. These three stages are monitored by the boundary conditions triggering respectively mechanical, hydraulic and thermal driving forces. As exhibited in Figure 7, the displacement u " t p5, ´2q 10 ´4 m is imposed at the top boundary for t ą 0. At the left boundary, a high pressure p " 8. 10 6 Pa is set for t ě t p1q , and a low temperature T " 285 K is prescribed for t ě t p2q .

The initial pressure and temperature are fixed to p 0 " 10 5 Pa and T 0 " 300 K. The initial displacement field u 0 is computed by solving the discrete contact-mechanical problem (25a) at given p 0 , T 0 with a zero displacement u " 0 at the bottom and top boundaries and a free boundary condition at the left and right sides. Throughout the simulation, the top and bottom boundaries are assumed impervious with zero heat flux while the right boundary condition imposes a fixed pressure p " 10 5 Pa and temperature T " 300 K. A zero displacement u " 0 is prescribed at the bottom boundary and the left and right sides are kept with homogeneous total stress conditions.

In order to compare the (H-model) and (S-model) models and their discretisations on different fluid behaviors, we consider in the following two cases corresponding first to a slightly compressible liquid and second to a perfect gas. The simulation parameters, not depending on the fluid thermodynamical model, are reported in Table 2 together with α ϕ " pb ´ϕ0 qα s , K s " λ `µ, N " Ks pb´ϕ 0 qp1´bq . The times t p1q , t p2q , t F will be fixed according to the fluid thermodynamical model in such a way that a stationary state is roughly reached at the end of each stage.

The simulations are performed on a family of 4 uniformly refined meshes indexed by m P t0, ¨¨¨, 3u with 2855ˆ4 m triangular cells and 88ˆ2 m fracture faces. The finest mesh m " 3 is used for the reference solution in the numerical convergence investigations. Due to the thermal convection dominated regime during stage 3, the convection terms are upwinded for the discretisation of both models. . Set up of the simulation in terms of initial condition at time t " 0 and of boundary conditions during each of the three stages corresponding to the time intervals I 1 " p0, t p1q q, I 2 " rt p1q , t p2q q and I 3 " rt p2q , t F s. Here 1 bar corresponds to 10 5 Pa.

Weakly compressible liquid case

We consider a weakly compressible liquid with thermodynamical constitutive laws deriving from a free enthalpy potential. It is characterised by its specific density ϱpp, T q such that ϱ ref ϱpp, T q " 1 ´pp ´pref q K f `αf pT ´Tref q, and its specific internal energy

epp, T q " C f T ´αf ϱ ref ´pp ´pref qT ref `ppT ´Tref q ¯`pp 2 ´p2 ref q 2ϱ ref K f ,
given the parameters T ref " 300 K, p ref " 10 5 Pa, ϱ ref " 10 3 Kg.m ´3, K f " 2.18 GPa, α f " 2.07 10 ´4K ´1, and C f " 4180 J.Kg ´1.K ´1. The fluid viscosity is set to η " 10 ´3 Pa s. The time intervals for each of the three stages are given by t p1q " 100 s, t p2q " 200 s and t F " 5 days. The time stepping is defined by a small initial time step of 0.1 s in order to capture the undrained regime time scale at the beginning of stage 1 and by the maximum time steps of 5 s for stage 1 and 2 and of 0.1 day for stage 3. A very good match is observed on all variables with only small differences during stage 3 which has been checked to result from the neglected terms V ¨∇p in the approximate entropy equations of the (S-model) model. When these terms are added back, the discrete solutions of both models fully match. Note that the linearisation of the Fourier term is not a significant source of discrepancy for this test case due to the thermal convection dominated regime.

Figure 11 compares, for the discrete solutions of both models, the convergence of the L 2 errors in time of mean in space variables as functions of the mesh step h (1/total number of fracture faces). The errors are computed w.r.t. the fine mesh m " 3 reference solutions obtained using the same time stepping and the same model. The observed convergence behaviors of the discretisations of the two models are very similar, with a rate of convergence roughly equal to 1.5.

Perfect gas case

We now consider the case of a perfect gas characterised by the following specific enthalpy and density hpp, T q " C f T, ϱpp, T q " M f p RT , with gas specific heat capacity C f " 1000 J.Kg ´1.K ´1, molar mass M f " 28.9645 10

´3

Kg.mol ´1 and the perfect gas constant R " 8.3149 J.mol ´1.K ´1. The gas viscosity is fixed to η " 1.72 10 ´5 Pa s.

The time intervals for each of the three stages are given by t p1q " 100 s, t p2q " 500 s and t F " 3.005 10 5 s. The time stepping is defined with a single time step during stage 1 since the coupling of the displacement with the pressure is very small during this stage due to the high gas compressibility (see the pressure solution during stage 1 in Figure 12 (a)). Stage 2 is initialised with a time step of 10 s and the maximum time steps are set to 100 s for stage 2 and 8000 s for stage 3.

In Figure 13 we present the L 2 -errors in time, vs. the mesh size, for the spacial means of relevant variables. Both schemes provide roughly a convergence order of 1.5. It can also be noticed that the enthalpy-based discrete solution is not very well captured by the coarser mesh, which explains the strong decrease of the errors between the first two meshes. The significant differences between the convergence plots of the two models result from the high discrepancy between the enthalpy-based and entropy-based discrete solutions.

Figure 14 better illustrates the difference between the solutions of both models; we see in particular that the entropy-based model exhibits a temperature variation far below ´15 K during stage 3 while it should physically be close to 15 K as it is the case for the enthalpy-based model. We demonstrate in Figure 15 that this discrepancy is due to the neglected terms V m ¨∇p in the approximate entropy equations of the (S-model) model. Once these terms are added, the discrete solutions of both models fully match. This can be easily explained by comparing during stage 3 the order of magnitude of the neglected term V m ¨∇p with that of ϱ m V m ¨∇h m . The ratio between both terms is of the order of roughly 10 for the gas test case while it is roughly 0.1 for the liquid test case. This is explained by the rather high pressure gradient combined with the low density and heat capacity in the gas case compared with the liquid case. We can conclude that the terms V m ¨∇p cannot be neglected in the entropy-based model in the gas case with low specific density.

Performances of the nonlinear solver

Figures 16 and 17 exhibit the total numbers of Newton iterations against time for the Thermo-Hydro and Mechanical models. The total number of time steps is 110 for the liquid case in Figure 16 and 87 for the gas case in Figure 17. We can notice the robustness in both cases of the nonlinear solvers w.r.t. the mesh size. Remarkably, in the liquid case, the entropy and enthalpy-based discrete models provide similar numbers of iterations through time. In the gas case, the number of Thermo-Hydro Newton iterations gets moderately larger during stage 3 for the entropy-based than for the enthalpy-based discrete models as a result of a much larger temperature variation for the entropy-based simulation. It has been checked that adding the V m ¨∇p terms to the entropy-based discrete model gives back essentially the same Newton behavior as for the enthalpy-based model.

Conclusion

This work focuses on discretisations of mixed-dimensional THM models preserving energy estimates for a general single phase fluid thermodynamical model. Our approach uses a Finite Volume scheme for the mass and energy equations with a possible upwinding of the convection terms in order to account for convection dominated regimes. It is combined with a mixed Finite Element discretisation of the contact-mechanics with face-wise con-stant Lagrange multipliers, keeping the dissipative property of the contact terms at matrix fracture interfaces. Two formulations of the energy equation are considered and compared numerically. It is built either directly from the energy conservation or obtained from an approximate entropy balance equation based on a small Darcy velocity and small temperature variation assumptions. The Finite Volume discretisation of the entropy-based model, and in particular of the non-conservative convection terms, is carefully designed in order to preserve the link between both formulations, which in turns guarantees that it satisfies an energy estimate. Both discrete models are assessed and compared in terms of convergence, accuracy and robustness on 2D test cases including a convective dominated regime, and either a weakly compressible liquid or highly compressible gas. It is shown that both approaches provides similar results in terms of spatial convergence and robustness of the nonlinear solver. On the other hand, in the gas case, for low specific density and high pressure gradient, the terms V ¨∇p that are typically neglected in the entropy balance approach must be accounted for in order to provide the physical solution.

Figure 1 .

 1 Figure1. Illustration of the dimension reduction in the fracture aperture for a 2D domain Ω with three intersecting fractures Γ i , i P t1, 2, 3u, with the equi-dimensional geometry on the left and the mixed-dimensional geometry on the right.

Figure 2 .

 2 Figure 2. Conceptual fracture model with contact at asperities, d 0 being the fracture aperture at contact state.

  For a given edge ζ P E Γ , let us denote by F Γ,ζ the subset of fracture faces sharing the edge ζ. We denote by |K| the d-dimensional measure of the cell K P M, and by |σ| the pd ´1q-dimensional measure of the face σ P F.

Figure 4 .Figure 5 .

 45 Figure 4. Convergence of the relative L 2 errors for the temperature T , pressure p, and displacement u and their gradients for the discretisation of the (S-model) model and both the centred (for K " I) and upwind (for K " I, 100 I) schemes.

Figure 6 .

 6 Figure 6. Two-dimensional domain Ω " p0, 1q ˆp0, 2q m 2 including the six fracture network Γ.

Figure 7

 7 Figure 7. Set up of the simulation in terms of initial condition at time t " 0 and of boundary conditions during each of the three stages corresponding to the time intervals I 1 " p0, t p1q q, I 2 " rt p1q , t p2q q and I 3 " rt p2q , t F s. Here 1 bar corresponds to 10 5 Pa.

Figure 9 Figure 8 .

 98 Figure 9 exhibits the evolution of the contact state (open, contact stick or contact slip) along the fractures at different times during the three stages. At time t ą 0, due to the imposed displacement at the top, most of the fractures switch from open to contact. During the undrained regime, at the very beginning of stage 1, the high increase of the pressure (see Figure 8 (a)) induces a slip state for most of the fractures in contact, as a result of the reduction of the normal surface traction. Toward the end of stage 1, these fractures switch back to stick state due to the pressure relaxation (see the evolution of the mean pressure in Figure 10 (a)). During stage 2, fractures 1 and 3 switch back to slip state as a result of the high pressure front propagation (see Figure 8 (b) and 10 (a)) while they partially

Figure 10

 10 Figure10compares the evolution in time of mean in space variables (p, p f , T , T f , and scaled d f ´d0 , u τ , and ϕ´ϕ 0 ) for the entropy-based and enthalpy-based discrete solutions obtained on the finest mesh m " 3. A very good match is observed on all variables with only small differences during stage 3 which has been checked to result from the neglected terms V ¨∇p in the approximate entropy equations of the (S-model) model. When these terms are added back, the discrete solutions of both models fully match. Note that the linearisation of the Fourier term is not a significant source of discrepancy for this test case due to the thermal convection dominated regime.

Figure 9 .Figure 11 .´1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 ´1Figure 10 .

 911610 Figure 9. Contact state (open, contact stick, contact slip) along the fractures at different times for the discrete solution on the finest mesh m " 3 of the enthalpy-based model.

Figure 12 .Figure 13 .Figure 14 .Figure 15 .

 12131415 Figure 12. Discrete solution obtained with the enthalpy-based model on the finest mesh m " 3: (a) pressure p at the end of stage 1, (b) pressure p at t " 170 s during stage 2, (c) temperature T at t " 42730 s during stage 3.

H´1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 0´1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 0Figure 16 .Figure 17 .

 661617 Figure 16. Total numbers of Newton iterations for the Thermo-Hydro and Mechanical models as a function of time, for the liquid case and both enthalpy-based (H) and entropybased (S) schemes. (Left) mesh m " 3, (right) mesh m " 1, with a total number of 110 time steps in all cases.

Table 1 .

 1 Material Properties

	Symbol	Quantity	Value	Unit

Table 2 .

 2 Material Properties common to the liquid and gas test cases.

	Symbol	Quantity	Value	Unit
	E	Young modulus	40	GPa
	ν	Poisson coefficient	0.15	F
		Friction coefficient	0.5	b
		Biot coefficient	0.65	Ḱ
		Permeability coefficient	˜1 0 0 0.5 ¸ˆ10 ´15 m 2
	ϕ 0	Initial porosity	0.1	d0
		Contact aperture	5. 10 ´4	m
	Λ m	Effective thermal conductivity	2	W m ´1 K ´1

Acknowledgements: the authors are grateful to Andra and BRGM for partially funding this work and to Laurence Beaude, Marc Leconte, Simon Lopez, Antoine Pasteau and Farid Smai for fruitful discussions during the elaboration of this work.

J. Droniou would like to acknowledge a partial funding by the European Union (ERC Synergy, NEMESIS, project number 101115663). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.


Resolution:

Writing a multiauthored manual on the topic of Dark Earth Target group: field archaeologists

Publication plan:

Facing the problem:

Meeting of specialists and discussing the phenomenon aim: make field archaeologists aware of the Earth phenomenon and to enable them to recognise, preserve and sample these soils

Earth Manual

  case studies with cultural interpretations/implications 5) Field guide (description and short graphical version) 6) Bibliography, information about the authors and useful addresses

  by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany´s Excellence Strategy -EXC 2150 -390870439. And Doing research: Plant pot trial to determine the real yield potential of Nordic DE Planned publication date: 2nd half 2024 Chapter 3: Methods Overview in Tables, short descriptions and graphic illustrations Example Example

Table 1 : Basic physical, chemical, and biological analysis of Dark Earth and their interpretation.

 1

	Analysis	Interpretation
	Particle size	Soil texture, parent material similar to surrounding soils, water holding
	distribution	capacity, nutrient holding capacity, water infiltration
	Sedimentology	Origin, deposition, transformation of parent material
	Munsell colour	Soil characteristics
	pH value	Soil reaction, plant growth, Aluminum toxicity, nutrient availability
	Total organic	Soil organic matter content, nutrient holding capacity, rooting capability, plant
	carbon	growth and yield, water holding capacity
	Total nitrogen	Soil organic matter quality, natural nitrogen availability, plant growth and yield
	Phosphorus	Past human activities such as bone and excrement deposition
	(fractions)
	Plant-available	Soil fertility and plant growth and yield
	nutrients

Table 2 : Sophisticated physical, chemical, and biological analysis of Dark Earth and their interpretation.

 2

	Analysis	Interpretation
	a. Micromorphology	Spatial distribution of human remains such as charcoal, bones etc.
		Natural soil processes such as clay movement
	b. Charcoal / Black	Quantity and quality of charring residues even if not visible as particulate
	carbon	charcoal
	d. Stable isotopes	Type of vegetation (woody vs. grass), meat consumption, fish
		consumption
	e. Urease	Presence of feces or manure
	f. Microorganisms	Presence of feces or manure
	g. Ancient soil DNA	Presence of humans, animals, plants
	(aDNA)
	h. Archaeobotany	Vegetation history
	i. Phytoliths	Vegetation history
	j. Radiocarbon	Determination of age / dating

c. Biomarkers Type of vegetation (phenols, terpenoids, cannabinoids), plant-vs. aquatic origin (n-alkanes), plant-vs. microbial origin (sugars), bacterial vs. fungal origin (amino sugars), microbial diversity (phospholipid fatty acids), human vs. animal excrements (steroids and bile acids) k. OSL Determination of age / dating