Francesco Bonaldi
email: francesco.bonaldi@univ-cotedazur.fr

Konstantin Brenner
email: konstantin.brenner@univ-cotedazur.fr

Jérôme Droniou
email: jerome.droniou@monash.edu

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

Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media

Keywords: 65M12, 76S05, 74B10 poromechanics, discrete fracture matrix models, two-phase Darcy flows, Gradient Discretization Method, convergence analysis

come

Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media

Introduction

Many real-life applications in geosciences involve processes like multi-phase flow and hydromechanical coupling in heterogeneous porous media. Such mathematical models are coupled systems of partial differential equations, including non-linear and degenerate parabolic ones. Besides the inherent difficulties posed by such equations, further complexities stem from the heterogeneity of the medium and the presence of discontinuities like fractures. This has a strong impact on the complexity of the models, challenging their mathematical and numerical analysis and the development of efficient simulation tools. This work focuses on the so called hybrid-dimensional matrix fracture models obtained by averaging both the unknowns and the equations across the fracture width and by imposing appropriate transmission conditions at the matrix fracture interfaces. Given the high geometrical complexity of real-life fracture networks, the main advantages of these hybriddimensional compared to full-dimensional models are to facilitate the mesh generation and the discretization of the model, and to reduce the computational cost of the resulting schemes. This type of hybrid-dimensional models has been the object of intensive researches over the last twenty years due to the ubiquity of fractures in geology and their large impact on flow, transport and mechanical behavior of rocks. For the derivation and analysis of such models, 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] Angot | Asymptotic and numerical modelling of flows in fractured 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] Nordbotten | Unified approach to discretization of flow in fractured porous media[END_REF] for single-phase Darcy flows, [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for two-phase flow in fractured porous media[END_REF][START_REF] Monteagudo | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF][START_REF] Jaffré | A discrete fracture model for two-phase flow with matrix-fracture interaction[END_REF][START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF][START_REF] Brenner | Hybrid dimensional modelling of two-phase flow through fractured with enhanced matrix fracture transmission conditions[END_REF][START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF] for two-phase Darcy flows, and [START_REF] Khoei | Numerical modeling of two-phase fluid flow in deformable fractured porous media using the extended finite element method and an equivalent continuum model[END_REF][START_REF] Kim | Rigorous coupling of geomechanics and multiphase flow with strong capillarity[END_REF][START_REF] Jha | Coupled Modeling of Multiphase Flow and Fault Poromechanics during geologic CO 2 storage[END_REF][START_REF] Girault | A lubrication fracture model in a poro-elastic medium[END_REF][START_REF] Hanowski | The hydromechanical equilibrium state of poroelastic media with a static fracture: A dimension-reduced model with existence results in weighted Sobolev spaces and simulations with an XFEM discretization[END_REF][START_REF] Garipov | Discrete fracture model for coupled flow and geomechanics[END_REF][START_REF] Jin | Fully coupled nonlinear fluid flow and poroelasticity in arbitrarily fractured porous media: A hybrid-dimensional computational model[END_REF][START_REF] Giovanardi | Unfitted FEM for modelling the interaction of multiple fractures in a poroelastic medium[END_REF][START_REF] Ucar | A finite-volume discretization for deformation of fractured media[END_REF] for poroelastic models.

In this article, we consider the two-phase Darcy flow in a network of pre-existing fractures represented as pd ´1qdimensional planar surfaces coupled with the surrounding d-dimensional matrix. The fractures are assumed to be open and filled by the fluids. Both phase pressures are assumed continuous across the fractures. This is a classical assumption for open fractures given the low pressure drop in the width of the fractures [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for two-phase flow in fractured porous media[END_REF][START_REF] Monteagudo | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF][START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF]. For single-phase flows, Poiseuille's law is classically used to model the flow along the fractures. This leads to a Darcy-like tangential flow with conductivity equal to

d 3 f 12
, where d f is the fracture aperture [START_REF] Girault | A lubrication fracture model in a poro-elastic medium[END_REF][START_REF] Hanowski | The hydromechanical equilibrium state of poroelastic media with a static fracture: A dimension-reduced model with existence results in weighted Sobolev spaces and simulations with an XFEM discretization[END_REF]. Following [START_REF] Khoei | Numerical modeling of two-phase fluid flow in deformable fractured porous media using the extended finite element method and an equivalent continuum model[END_REF], the extension to a two-phase flow is based on the generalized Darcy laws involving appropriate relative permeabilities and the capillary pressure-saturation relation. This hybrid-dimensional two-phase Darcy flow model is coupled with the matrix mechanical deformation assuming small strains and a linear poroelastic behavior [START_REF] Khoei | Numerical modeling of two-phase fluid flow in deformable fractured porous media using the extended finite element method and an equivalent continuum model[END_REF][START_REF] Kim | Rigorous coupling of geomechanics and multiphase flow with strong capillarity[END_REF][START_REF] Jha | Coupled Modeling of Multiphase Flow and Fault Poromechanics during geologic CO 2 storage[END_REF]. The extension of the single-phase poromechanical coupling [START_REF] Girault | A lubrication fracture model in a poro-elastic medium[END_REF][START_REF] Hanowski | The hydromechanical equilibrium state of poroelastic media with a static fracture: A dimension-reduced model with existence results in weighted Sobolev spaces and simulations with an XFEM discretization[END_REF][START_REF] Garipov | Discrete fracture model for coupled flow and geomechanics[END_REF][START_REF] Jin | Fully coupled nonlinear fluid flow and poroelasticity in arbitrarily fractured porous media: A hybrid-dimensional computational model[END_REF][START_REF] Ucar | A finite-volume discretization for deformation of fractured media[END_REF] to two-phase Darcy flows is based on the so-called equivalent pressure used both in the matrix for the total stress and at both sides of the fractures as boundary condition for the mechanics. Typically, the equivalent pressure is defined as a convex combination of the phase pressures and several different combinations have been proposed in the literature [START_REF] Nuth | Effective stress concept in unsaturated soils: Clarification and validation of a unified framework[END_REF]. Our choice of the equivalent pressure follows the pioneer monograph by Coussy [START_REF] Coussy | Poromechanics[END_REF] and involves the capillary energy which, as already noticed in [START_REF] Kim | Rigorous coupling of geomechanics and multiphase flow with strong capillarity[END_REF][START_REF] Jha | Coupled Modeling of Multiphase Flow and Fault Poromechanics during geologic CO 2 storage[END_REF], plays a key role to obtain energy estimates for the coupled system. From the open fracture assumption, the fracture mechanical behavior reduces to the continuity of the normal stresses at both sides of the fracture matching with the fracture equivalent pressure times the unit normal vector. To our best knowledge, no theoretical or numerical analysis of the complete poromechanical model, with all non-linear coupling, has been carried out so far.

In this work, the hybrid-dimensional coupled model is discretized using the gradient discretization method (GDM) [START_REF] Droniou | The Gradient Discretisation Method[END_REF]. This framework is based on abstract vector spaces of discrete unknowns combined with reconstruction operators. The gradient scheme is then obtained by substitution of the continuous operators by their discrete counterparts in the weak formulation of the coupled model. The main asset of this framework is to enable a generic convergence analysis based on general properties of the discrete operators that hold for a large class of conforming and non conforming discretizations. Two essential ingredients to discretize the coupled model are the discretizations of the hybrid-dimensional two-phase Darcy flow and the discretization of the mechanics. Let us briefly mention, in both cases, a few families of discretizations typically satisfying the gradient discretization properties. For the discretization of the Darcy flow, the gradient discretization framework typically covers the case of cell-centered finite volume schemes with Two-Point Flux Approximation on strongly admissible meshes [START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF][START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF][START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF], or some symmetric Multi-Point Flux Approximations [START_REF] Tunc | A model for conductive faults with non matching grids[END_REF][START_REF] Sandve | An efficient multi-point flux approximation method for discrete fracture-matrix simulations[END_REF][START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF] on tetrahedral or hexahedral meshes. It also accounts for the families of Mixed Hybrid Mimetic and Mixed or Mixed Hybrid Finite Element discretizations such as in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[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] Antonietti | Mimetic finite difference approximation of flows in fractured porous media[END_REF]. The case of vertexbased discretizations such as Control Volume Finite Element approaches (i.e. conforming finite element with mass lumping) [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF][START_REF] Reichenberger | A mixed-dimensional finite volume method for two-phase flow in fractured porous media[END_REF][START_REF] Monteagudo | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF] or the Vertex Approximate Gradient scheme [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] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF][START_REF] Brenner | Hybrid dimensional modelling of two-phase flow through fractured with enhanced matrix fracture transmission conditions[END_REF] is also accounted for. For the discretization of the elastic mechanical model, the gradient discretization framework covers conforming finite element methods such as in [START_REF] Girault | A lubrication fracture model in a poro-elastic medium[END_REF], as well as the Crouzeix-Raviart discretization [START_REF] Hansbo | Discontinuous Galerkin and the Crouzeix-Raviart element: Application to elasticity[END_REF][START_REF] Pietro | An extension of the Crouzeix-Raviart space to general meshes with application to quasi-incompressible linear elasticity and Stokes flow[END_REF], Discontinuous Galerkin methods [START_REF] Eymard | Discontinuous Galerkin gradient discretisations for the approximation of secondorder differential operators in divergence form[END_REF], the Hybrid High Order discretization [START_REF] Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF], and the Virtual Element Method [START_REF] Beirão Da | Virtual elements for linear elasticity problems[END_REF]. Note that many of these methods are actually applicable to both the flow and the mechanical component of the model. Without taking into account the poromechanical coupling, convergence results have been obtained in [START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF][START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[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] for hybrid-dimensional single-phase Darcy flow models, and in [START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF] for hybrid-dimensional two-phase Darcy flow models. The well-posedness and convergence analysis of single-phase poromechanical models is studied in [START_REF] Girault | A lubrication fracture model in a poro-elastic medium[END_REF][START_REF] Hanowski | The hydromechanical equilibrium state of poroelastic media with a static fracture: A dimension-reduced model with existence results in weighted Sobolev spaces and simulations with an XFEM discretization[END_REF]. Nevertheless those analyses consider a linear approximation of the coupled model obtained by freezing the fracture conductivity d 3 f 12 , and hence eliminating the non-linear coupling between the fracture aperture and the Darcy flow. Let us also mention the related recent work [START_REF] Both | Global existence of a weak solution to unsaturated poroelasticity[END_REF] on unsaturated poroelasticity based on the Richards approximation of the two-phase flow model, using partial linearizations, non-degeneracy conditions and Kirchhoff transformation (which is made possible by assuming that the saturation-capillary pressure law is uniform across the domain). Note that fractures are not considered in this work.

Our main result is the proof of convergence, in the GDM setting, of the approximate solutions to the weak solution of the non-linear coupled model with two-phase flows. To our best knowledge, this is the first convergence result for this type of hybrid-dimensional model taking into account the full non-linear poromechanical coupling. Since it is based on discrete compactness techniques, the convergence is that of a subsequence of approximate solutions (precisely, we prove that sequences of approximate solutions are compact, and that any of their limit points is a weak solution of the continuous model). To establish this result, we make the following main assumptions. It is first assumed that the approximate matrix porosity remains bounded below by a strictly positive constant and that the approximate fracture aperture remains larger than some given aperture vanishing only at the tips. Let us point out that these assumptions are due to the limitations of the model itself rather than to the shortcomings of the numerical analysis. They cannot be avoided since the continuous model does not ensure the positivity of the porosity nor of the fracture aperture, properties needed to guarantee existence of solutions. We note that previous works on similar models circumvent these limitations by linearization processes (complete or partial freezing of the matrix porosity and fracture apertures). Regarding the assumption on the fracture aperture, it could possibly be overcome by introducing contact mechanics in the model [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]. This direction will be investigated in a future work. It is also assumed in the numerical analysis that the mobility functions are bounded below by strictly positive constants. Independently of the poromechanical coupling, this is a classical assumption to enable the stability and convergence analysis of two-phase Darcy flows with spatial discontinuity of the capillary pressure functions, as it is always the case in the presence of fractures (see [START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF][START_REF] Brenner | Vertex Approximate Gradient Scheme for Hybrid Dimensional Two-Phase Darcy Flows in Fractured Porous Media[END_REF][START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF]). To our knowledge, the only convergence analyses covering both the degeneracy of the mobilities and discontinuous capillary pressures are limited to Two-Point Flux Approximations (see [START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF][START_REF] Brenner | Total-velocity-based finite volume discretization of twophase darcy flow in highly heterogeneous media with discontinuous capillary pressure[END_REF]). Extending such analyses, even considering only the TPFA method for the flow, to the poromechanical model considered here is far from straightforward and seems to bring additional challenges; given that our analysis is already quite technical, we postpone this extension to degenerate mobility functions to a future work.

The rest of the article is organized as follows. Section 2 introduces the continuous hybrid-dimensional coupled model. Section 3 describes the gradient discretization method for the coupled model including the definition of the reconstruction operators, the discrete variational formulation and the properties of the gradient discretization needed for the subsequent convergence analysis. Section 4 proceeds with the convergence analysis. The a priori estimates are established in Subsection 4.1, the compactness properties in Subsection 4.2 and the convergence to a weak solution is proved in Subsection 4.3. In Section 5, numerical experiments based on the Two-Point Flux Approximation finite volume scheme for the flows and second-order finite elements for the mechanical deformation are carried out for a cross-shaped fracture network in a two-dimensional porous medium, and illustrate the numerical convergence of the solution. Appendices A.1 and A.2 state some technical results used in the convergence analysis.

Continuous model

We consider a bounded polytopal domain Ω of R d , d P t2, 3u, 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 polygonal simply connected open domain with angles strictly lower than 2π. Without restriction of generality, we will assume that the fractures may intersect exclusively at their boundaries, that is for any i, j P I, i ‰ j one has Γ i X Γ j " H, but not necessarily Γ i X Γ j " H. Since one can split a general (non-simply connected) planar polygon into several simply connected pieces intersecting only at their boundaries (see Figure 1) our assumptions on the fracture network are in fact quite general. Roughly speaking we only exclude the non-planar fractures.

Since the fractures are assumed open with no contact, we also have to assume in the following that the boundary of each connected component of ΩzΓ has a non zero measure intersection with BΩ.

The two sides of a given fracture of Γ are denoted by ˘in the matrix domain, with unit normal vectors n ˘oriented outward of the sides ˘. We denote by γ the trace operator on Γ for functions in H 1 pΩq, by γ BΩ the trace operator for the same functions on BΩ, and by ¨ the normal trace jump operator on Γ for functions in H div pΩzΓq, defined by ū " ū`¨n``ū´¨n´f or all ū P H div pΩzΓq.

We denote by ∇ τ the tangential gradient and by div τ the tangential divergence on the fracture network Γ. The symmetric gradient operator ε is defined such that εpvq " 1 2 p∇v `t p∇vqq for a given vector field v P H 1 pΩzΓq d . The fracture aperture, denoted by df , is defined by df " ´ ū for a displacement field ū P H 1 pΩzΓq d .

Let us fix a continuous function d 0 : Γ Ñ p0, `8q vanishing at BΓzpBΓ X BΩq (i.e. at the tips of Γ) and taking strictly positive values at BΓ X BΩ. The discrete fracture aperture will be assumed to be greater than or equal to d 0 almost everywhere (by the established convergence result, the same will hold for its limit). We note that the assumptions on d 0 are minimal, allowing for very general behavior of the fracture aperture at the tips.

Let us introduce some relevant function spaces:

U 0 " tv P pH 1 pΩzΓqq d | γ BΩ v " 0u (1)
Figure 2: Example of a 2D domain Ω with its fracture network Γ, the unit normal vectors n ˘to Γ, the phase pressures pα in the matrix and γ pα in the fracture network, the displacement vector field ū, the matrix Darcy velocities q α m and the fracture tangential Darcy velocities q α f integrated along the fracture width.

for the displacement vector, and

V 0 " tv P H 1 0 pΩq | γv P H 1 d0 pΓqu (2)
for each phase pressure, where the space

H 1 d0 pΓq is made of functions v Γ in L 2 pΓq, such that d 3 {2 0 ∇ τ v Γ is in L 2 pΓq d´1
, and whose traces are continuous at fracture intersections BΓ i X BΓ j , pi, jq P I ˆI (i ‰ j) and vanish on the boundary BΓ X BΩ.

The matrix and fracture rock types are denoted by the indices rt " m and rt " f , respectively, and the non-wetting and wetting phases by the superscripts α " nw and α " w, respectively. Each rock type rt P tm, f u is characterized by its own set of mobility functions pη α rt q αPtnw,wu and capillary pressure-saturation relation pS α rt q αPtnw,wu .

The PDEs model reads: find the phase pressures pα , α P tnw, wu, and the displacement vector field ū, both satisfying homogeneous Dirichlet boundary conditions on BΩ, such that pc " pnw ´p w and, for α P tnw, wu,

$ ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' % B t `φ m S α m pp c q ˘`div pq α m q " h α m on p0, T q ˆΩzΓ, q α m " ´ηα m pS α m pp c qqK m ∇p α on p0, T q ˆΩzΓ, B t ´d f S α f pγ pc q ¯`div τ pq α f q ´ q α m " h α f on p0, T q ˆΓ, q α f " ´ηα f pS α f pγ pc qqp 1 12 d3
f q∇ τ γ pα on p0, T q ˆΓ, ´div ´σpūq ´b pE m I ¯" f on p0, T q ˆΩzΓ σpūq " 2µ εpūq `λ divpūq I on p0, T q ˆΩzΓ,

$ ' ' & ' ' % B t φm " b divB t ū `1 M B t pE m on p0, T q ˆΩzΓ, pσpūq ´b pE m Iqn ˘" ´p E f n ˘on p0, T q ˆΓ, df " ´ ū on p0, T q ˆΓ, (3) with
and the initial conditions pα | t"0 " pα 0 , φm | t"0 " φ0 m . Here, we have denoted by pc the capillary pressure, and the equivalent pressures pE m and pE f are defined, following [START_REF] Coussy | Poromechanics[END_REF], by pE m " ÿ αPtnw,wu pα S α m pp c q ´Um pp c q, pE f " ÿ αPtnw,wu γ pα S α f pγ pc q ´Uf pγ pc q, where U rt pp c q "

ż pc 0 z pS nw rt q 1 pzq dz (5)
is the capillary energy density function of the rock type rt P tm, f u. As already noticed in [START_REF] Kim | Rigorous coupling of geomechanics and multiphase flow with strong capillarity[END_REF][START_REF] Jha | Coupled Modeling of Multiphase Flow and Fault Poromechanics during geologic CO 2 storage[END_REF], this is a key choice to obtain the energy estimates that are the starting point for the convergence analysis.

We make the following main assumptions on the data:

(H1) For each phase α P tnw, wu and rock type rt P tm, f u, the mobility function η α rt is continuous, non-decreasing, and there exist 0 ă η α rt,min ď η α rt,max ă `8 such that η α rt,min ď η α rt psq ď η α rt,max for all s P r0, 1s.

(H2) For each rock type rt P tm, f u, the non-wetting phase saturation function S nw rt is a non-decreasing Lipschitz continuous function with values in r0, 1s, and S w rt " 1 ´Snw rt . (H3) b P r0, 1s is the Biot coefficient, M ą 0 is the Biot modulus, and λ ą 0, µ ą 0 are the Lamé coefficients. These coefficients are assumed to be constant for simplicity.

(H4) The initial pressures are such that pα 0 P V 0 X L 8 pΩq and γ pα 0 P L 8 pΓq, α P tnw, wu; the initial porosity is such that φ0 m P L 8 pΩq. (H5) The source terms satisfy f P L 2 pΩq d , h α m P L 2 pp0, T q ˆΩq, and h α f P L 2 pp0, T q ˆΓq. (H6) The matrix permeability tensor K m is symmetric and uniformly elliptic on Ω. Note that the variation of the matrix permeability with the porosity is neglected.

The notion of weak solution for (3)-( 4) is classically obtained multiplying each flow equation and the mechanical equation by a separate test function, integrating by parts and, for each phase, adding together the equations resulting from the flows in the matrix and the fractures. When the capillary pressure has continuous first temporal and second spatial derivatives in p0, T q ˆpΩzΓq, its trace has continuous first temporal and second tangential derivatives in p0, T qˆΓ, and the displacement has continuous second spatial derivatives, the following weak formulation is equivalent to the PDE model.

Definition 2.1 (Weak solution of the model). A weak solution of the model is given by pα P L 2 p0, T ; V 0 q, α P tnw, wu, and ū P L 8 p0, T ; U 0 q, such that, for any α P tnw, wu, d 3 {2 f ∇ τ γ pα P L 2 pp0, T q ˆΓq d´1 and, for all φα P C 8 c pr0, T q ˆΩq and all smooth functions v : r0, T s ˆpΩzΓq Ñ R d vanishing on BΩ and whose derivatives of any order admit finite limits on each side of Γ,

ż T 0 ż Ω ´´φ m S α m pp c qB t φα `ηα m pS α m pp c qqK m ∇p α ¨∇ φα ¯dxdt `ż T 0 ż Γ ´´d f S α f pγ pc qB t γ φα `ηα f pS α f pγ pc qq d 3 f 12 ∇ τ γ pα ¨∇τ γ φα ¯dσpxqdt ´żΩ φ0 m S α m pp 0 c q φα p0, ¨qdx ´żΓ d0 f S α f pγ p0 c qγ φα p0, ¨qdσpxq " ż T 0 ż Ω h α m φα dxdt `ż T 0 ż Γ h α f γ φα dσpxqdt, (6)
ż T 0 ż Ω ´σpūq : εpvq ´b pE m divpvq ¯dxdt `ż T 0 ż Γ pE f v dσpxqdt " ż T 0 ż Ω f ¨v dxdt, (7)
with pc " pnw ´p w , df " ´ ū , φm ´φ 0 m " b divpū ´ū 0 q `1 M pp E m ´p E,0 m q, d0 f " ´ ū0 , where ū0 is the solution of (7) without the time integral and using the initial equivalent pressures pE,0 m and pE,0 f obtained from the initial pressures pα 0 and γ pα 0 , α P tnw, wu. Remark 2.2 (Regularity of the fracture aperture). Notice that, by the Sobolev-trace embeddings [1, Theorem 4.12], ū P L 8 p0, T ; U 0 q implies that df " ´ ū P L 8 p0, T ; L 4 pΓqq. All the integrals above are thus well-defined.

The gradient discretization method

The gradient discretization (GD) for the Darcy continuous pressure model, introduced in [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF], is defined by a finitedimensional vector space of discrete unknowns X 0 Dp and • two discrete gradient linear operators on the matrix and fracture domains

∇ m Dp : X 0 Dp Ñ L 8 pΩq d , ∇ f Dp : X 0 Dp Ñ L 8 pΓq d´1 ;
• two function reconstruction linear operators on the matrix and fracture domains

Π m Dp : X 0 Dp Ñ L 8 pΩq, Π f Dp : X 0 Dp Ñ L 8 pΓq,
which are piecewise constant [START_REF] Droniou | The Gradient Discretisation Method[END_REF]Definition 2.12].

A consequence of the piecewise-constant property is the following: there is a basis pe i q iPI of X 0 Dp such that, if v " ř iPI v i e i and if, for a mapping g : R Ñ R with gp0q " 0, we define gpvq " ř iPI gpv i qe i P X 0 Dp by applying g component-wise, then Π rt Dp gpvq " gpΠ rt Dp vq for rt P tm, f u. Note that the basis pe i q iPI is usually canonical and chosen in the design of X 0 Dp . The vector space X 0 Dp is endowed with

}v} Dp -}∇ m Dp v} L 2 pΩq `}d 3 {2 0 ∇ f Dp v} L 2 pΓq ,
assumed to define a norm on X 0 Dp .

The gradient discretization for the mechanics is defined by a finite-dimensional vector space of discrete unknowns X Remark 3.1 (On the boundary conditions). The exponent 0 in the spaces means that homogeneous Dirichlet boundary conditions are encoded in these spaces. We restrict our analysis to these boundary conditions for simplicity but, as shown in [START_REF] Droniou | The Gradient Discretisation Method[END_REF], the GDM analysis can easily be adapted to other types of boundary conditions -in particular to mixed Dirichlet/Neumann boundary conditions (with non-homogeneous Dirichlet values) as used for the flow part of the model in the numerical tests of Section 5.

A spatial GD can be extended into a space-time GD by complementing it with

• a discretization 0 " t 0 ă t 1 ă ¨¨¨ă t N " T of the time interval r0, T s;

• interpolators I Dp : V 0 Ñ X 0 Dp and I m Dp : L 2 pΩq Ñ X 0 Dp of initial conditions.

For n P t0, . . . , N u, we denote by δt n`1 2 " t n`1 ´tn the time steps, and by ∆t " max n"0,...,N δt n`1 2 the maximum time step.

The spatial operators are extended into space-time operators as follows. Let χ represent either p or u. If w " pw n q N n"0 P pX 0 Dχ q N `1, and Ψ Dχ is a spatial GD operator, its space-time extension is defined by Ψ Dχ wp0, ¨q " Ψ Dχ w 0 and, @n P t0, . . . , N ´1u , @t P pt n , t n`1 s, Ψ Dχ wpt, ¨q " Ψ Dχ w n`1 .

For convenience, the same notation is kept for the spatial and space-time operators. Moreover, we define the discrete time derivative as follows: for f : r0, T s Ñ L 1 pΩq piecewise constant on the time discretization, with f n " f |ptn´1,tns and f 0 " f p0q, we set δ t f ptq " fn`1´fn δt n`1 2 for all t P pt n , t n`1 s, n P t0, . . . , N ´1u.

Notice that the space of piecewise constant X 0 Dχ -valued functions f on the time discretization together with the initial value f 0 " f p0q can be identified with pX 0 Dχ q N `1. The same definition of discrete derivative can thus be given for an element w P pX 0 Dχ q N `1. Namely, δ t w P pX 0 Dχ q N is defined by setting, for any n P t0, . . . , N ´1u and t P pt n , t n`1 s, δ t wptq " pδ t wq n`1 -wn`1´wn δt n`1 2

. If Ψ Dχ pt, ¨q is a space-time GD operator, by linearity the following commutativity property holds: Ψ Dχ δ t wpt, ¨q " δ t pΨ Dχ wpt, ¨qq.

The gradient scheme for the system consists in replacing the "continuous" functional space and differential operators in (6)- [START_REF] Antonietti | Mimetic finite difference approximation of flows in fractured porous media[END_REF] by their discrete counterparts. This results in the following discrete problem: find p α P pX 0 Dp q N `1, α P tnw, wu, and u P pX 0 Du q N `1, such that for all ϕ α P pX 0 Dp q N `1, v P pX 0 Du q N `1 and α P tnw, wu,

ż T 0 ż Ω ´δt ´φD Π m Dp s α m ¯Πm Dp ϕ α `ηα m pΠ m Dp s α m qK m ∇ m Dp p α ¨∇m Dp ϕ α ¯dxdt `ż T 0 ż Γ δ t ´df,Du Π f Dp s α f ¯Πf Dp ϕ α dσpxqdt `ż T 0 ż Γ η α f pΠ f Dp s α f q d 3 f,Du 12 ∇ f Dp p α ¨∇f Dp ϕ α dσpxqdt " ż T 0 ż Ω h α m Π m Dp ϕ α dxdt `ż T 0 ż Γ h α f Π f Dp ϕ α dσpxqdt, (8a)
ż T 0 ż Ω ´σDu puq : ε Du pvq ´b Π m Dp p E m div Du pvq ¯dxdt `ż T 0 ż Γ Π f Dp p E f v Du dσpxqdt " ż T 0 ż Ω f ¨ΠDu v dxdt, (8b)
with the closure equations

$ ' ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' ' % p c " p nw ´pw , s α m " S α m pp c q, s α f " S α f pp c q, p E m " ÿ αPtnw,wu p α s α m ´Um pp c q, p E f " ÿ αPtnw,wu p α s α f ´Uf pp c q, φ D ´Πm Dp φ 0 m " b div Du pu ´u0 q `1 M Π m Dp pp E m ´pE,0 m q, d f,Du " ´ u Du , σ Du pvq " 2µε Du pvq `λ div Du pvqI. (8c)
The initial conditions are given by p α 0 " I Dp pα 0 (α P tnw, wu), φ 0 m " I m Dp φ0 m , and the initial displacement u 0 is the solution of (8b) without the time variable and with the equivalent pressures obtained from the initial pressures pp α 0 q αPtnw,wu . Remark 3.2 (Non-homogeneous boundary conditions). The homogeneous Dirichlet boundary conditions are embedded in the discrete spaces X 0 Dp and X 0 Du . Non-homogeneous (or other types of) boundary conditions are equally easy to handle in the GDM setting [30, Section 2.2 and Chapter 3].

Remark 3.3 (GDM framework). As shown above, the GDM framework enables a presentation of the schemes in a way that is almost as compact as the weak formulation itself (compare with Definition 2.1). This presentation is valid for conforming methods, that already have a compact writing but may not be the best suited in practical applications (especially for the flow component), but also for non-conforming methods of practical interest in engineering; explicitly writing, for example, the TPFA formulation for the flow component of the model would lead to much lengthier equations. Additionally, the GDM analysis is also carried out in a compact way, identifying key properties and manipulating discrete equations almost as their continuous counterparts; notwithstanding the fact that this analysis applies to many different methods at once, developing it for a given specific scheme would not lead to any simplification -the complexity in the upcoming analysis comes from the poromechanical model we consider, not from the numerical analysis framework we use.

Properties of gradient discretizations

Let pD l p q lPN and pD l u q lPN be sequences of GDs. We state here the assumptions on these sequences which ensure that the solutions to the corresponding schemes converge. Most of these assumptions are adaptation of classical GD assumptions [START_REF] Droniou | The Gradient Discretisation Method[END_REF], except for the chain-rule, product rule and cut-off properties used in Subsection 4.2 to obtain compactness properties; we note that all these assumptions hold for standard discretizations used in porous media flows.

Following [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF], the spatial GD of the Darcy flow

D p " ´X0 Dp , ∇ m Dp , ∇ f Dp , Π m Dp , Π f
Dp ¯is assumed to satisfy the following coercivity, consistency, limit-conformity and compactness properties.

Coercivity of D p . Let C Dp ą 0 be defined by

C Dp " max 0‰vPX 0 Dp }Π m Dp v} L 2 pΩq `}Π f Dp v} L 2 pΓq }v} Dp . (9)
Then, a sequence of spatial GDs pD l p q lPN is said to be coercive if there exists C p ą 0 such that C D l p ď C p for all l P N. Consistency of D p . Let r ą 8 be given, and for all w P V 0 and v P X 0 Dp let us define

S Dp pw, vq " }∇ m Dp v ´∇w} L 2 pΩq `}∇ f Dp v ´∇τ γw} L r pΓq `}Π m Dp v ´w} L 2 pΩq `}Π f Dp v ´γw} L r pΓq , (10)
and S Dp pwq " min vPX 0 Dp S Dp pw, vq. Then, a sequence of spatial GDs pD l p q lPN is said to be consistent if for all w P V 0 one has lim lÑ`8 S D l p pwq " 0. Moreover, if pD l p q lPN is a sequence of space-time GDs, then it is said to be consistent if the underlying sequence of spatial GDs is consistent as above, and if, for any ϕ P V 0 and ψ P L 2 pΩq, as l Ñ `8,

∆t l Ñ 0 , }Π m D l p I D l p ϕ ´ϕ} L 2 pΩq `}Π f D l p I D l p ϕ ´ϕ} L 2 pΓq Ñ 0 and }Π m D l p I m D l p ψ ´ψ} L 2 pΩq Ñ 0. ( 11
)
Remark 3.4 (Consistency). In [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF], the consistency is only considered for r " 2. As it will appear clear in the analysis, dealing with the coupling and non-linearity of the model requires us to adopt here a slightly stronger consistency assumption. Under standard mesh regularity assumptions, this stronger consistency property is still satisfied for all classical GDs [30, Part III].

Limit-conformity of D p . For all pr m , r f q P C 8 pΩzΓq d ˆC8 pΓq d´1 and v P X 0 Dp , let us define

W Dp pr m , r f , vq " ż Ω ´rm ¨∇m Dp v `Πm Dp v divpr m q ¯dx `żΓ ´rf ¨∇f Dp v `Πf Dp v pdiv τ pr f q ´ r m q ¯dσpxq, (12)
and

W Dp pr m , r f q " max 0‰vPX 0 Dp |W Dp pr m , r f , vq| }v} Dp .
Then, a sequence of spatial GDs pD l p q lPN is said to be limit-conforming if for all pr m , r f q P C 8 pΩzΓq d ˆC8 c pΓq d´1 one has lim lÑ`8 W D l p pr m , r f q " 0. Here C 8 c pΓq d´1 denotes the space of functions whose restriction to each Γ i is in C 8 pΓ i q d´1 tangent to Γ i , compactly supported away from the tips, and satisfying normal flux conservation at fracture intersections not located at the boundary BΩ. Remark 3.5 (Compactly supported fluxes). The role of pr m , r f q is that of test functions (they do not represent the continuous fluxes), to show that the limits of the discrete fluxes are indeed the continuous fluxes, see [START_REF] Brenner | Gradient discretization of hybrid-dimensional Darcy flows in fractured porous media[END_REF]Lemma 5.5].

(Local) compactness of D p . A sequence of spatial GDs pD l p q lPN is said to be locally compact if for all sequences pv l q lPN P pX 0 D l p q lPN such that sup lPN }v l } D l p ă `8 and all compact sets K m Ă Ω and K f Ă Γ, such that K f is disjoint from the intersections pΓ i X Γ j q i "j , the sequences pΠ m D l p v l q lPN and pΠ f D l p v l q lPN are relatively compact in L 2 pK m q and L 2 pK f q, respectively. Remark 3.6 (Local compactness through estimates of space translates). For K m , K f as above, set

T D l p ,Km,K f pξ, ηq " max vPX 0 D l p zt0u }Π m D l p vp¨`ξq ´Πm D l p v} L 2 pKmq `řiPI }Π f D l p vp¨`η i q ´Πf D l p v} L 2 pK f XΓiq }v} D l p ,
where ξ P R d , η " pη i q iPI with η i tangent to Γ i ; for ξ and η small enough, this expression is well defined since K m and K f are compact in Ω and Γ, respectively. Following [30, Lemma 2.21], an equivalent formulation of the local compactness property is: for all K m , K f as above,

lim ξ,ηÑ0 sup lPN T D l p ,Km,K f pξ, ηq " 0.
Remark 3.7 (Usual compactness property for GDs). The standard compactness property for GD is not local but global, that is, on the entire domain and not any of its compact subsets (see, e.g., [START_REF] Droniou | The Gradient Discretisation Method[END_REF]Definition 2.8] and also below for D u ). Two reasons pushed us to consider here the weaker notion of local compactness: firstly, for standard GDs, the global compactness does not seem obvious to establish (or even true) in the fractures, because of the weight d 0 in the norm } ¨}Dp , which prevents us from estimating the translates of the reconstructed function by the gradient near the fracture tips; secondly, we will only prove compactness on saturations, which are uniformly bounded by 1 and for which local and global compactness are therefore equivalent.

In the following, for brevity we refer to the local compactness of pD l p q lPN simply as the compactness of this sequence of GDs.

Chain rule estimate on pD l p q lPN : for any Lipschitz-continuous function F : R Ñ R, there is C F ě 0 such that, for all l P N, v P X 0

D l p , }∇ m D l p F pvq} L 2 pΩq ď C F }∇ m D l p v} L 2 pΩq .
Product rule estimate on pD l p q lPN : there exists C P such that, for any l P N and any u l , v l P X 0

D l p , it holds }∇ m Dp pu l v l q} L 2 pΩq ď C P ´|u l | 8 }∇ m Dp v l } L 2 pΩq `|v l | 8 }∇ m Dp u l } L 2 pΩq ¯,
where |w| 8 -max iPI |w i | whenever w " ř iPI w i e i with pe i q iPI the canonical basis of X 0

D l p .
Cut-off property of pD l p q lPN : for any compact set K Ă ΩzΓ, there exists C K ě 0 and pψ l q lPN P pX 0 D l p q lPN such that p|ψ l | 8 q lPN is bounded and, for l large enough:

Π m D l p ψ l ě 0 on Ω; Π m D l p ψ l " 1 on K; }∇ m D l p ψ l } L 2 pΩq ď C K Π f D l p pv l ψ l q " 0 and ∇ f D l p pv l ψ l q " 0 for all v l P X 0 D l p
Coercivity of pD l u q lPN . Let C Du ą 0 be defined by

C Du " max 0‰vPX 0 D l u }Π D l u v} L 2 pΩq `} v D l u } L 4 pΓq }v} D l u . (13)
Then, the sequence of spatial GDs pD l u q lPN is said to be coercive if there exists C u ą 0 such that C D l u ď C u for all l P N.

Consistency of pD l

u q lPN . For all w P U 0 , it holds lim lÑ`8 S D l u pwq " 0 where

S D l u pwq " min vPX 0 D l u " }ε D l u pvq ´εpwq} L 2 pΩ,S d pRqq `}Π D l u v ´w} L 2 pΩq `› › v D l u ´ w › › L 4 pΓq ı . (14)
Limit-conformity of pD l u q lPN . Let C 8 Γ pΩzΓ, S d pRqq denote the vector space of smooth functions τ : ΩzΓ Ñ S d pRq whose derivatives of any order admit finite limits on each side of Γ, and such that τ `pxqn ``τ ´pxqn ´" 0 and pτ `pxqn `qˆn `" 0 for a.e. x P Γ. For all τ P C 8 Γ pΩzΓ, S d pRqq, it holds lim lÑ`8 W D l u pτq " 0 where

W D l u pτq " max 0‰vPX 0 D l u 1 }v} D l u "ż Ω ´τ : ε D l u pvq `ΠD l u v ¨divpτq ¯dx ´żΓ pτn `q ¨n` v D l u dσpxq  .
Compactness of pD l u q lPN . For any sequence pv l q lPN P pX 0 D l u q lPN such that sup lPN }v l } D l u ă `8, the sequences pΠ D l u v l q lPN and p v l D l u q lPN are relatively compact in L 2 pΩq d and in L s pΓq for all s ă 4, respectively. Remark 3.8 (Compactness through estimates of space translates). Similarly to Remark 3.6 (see also [START_REF] Droniou | The Gradient Discretisation Method[END_REF]Lemma 2.21]), the compactness of pD l u q lPN is equivalent to

lim ξ,ηÑ0 sup lPN T D l u ,s pξ, ηq " 0 @s ă 4,
where

T D l u ,s pξ, ηq " max vPX 0 D l u zt0u }Π D l u vp¨`ξq ´ΠD l u v} L 2 pΩq `řiPI › › v l D l u p¨`η i q ´ v l D l u › › L s pΓiq }v} D l u ,
with ξ P R d , η " pη i q iPI with η i tangent to Γ i , and the functions extended by 0 outside their respective domain Ω or Γ.

Convergence analysis

The main result of this work is the following theorem stating the convergence of the sequence of discrete solutions to a weak solution up to a subsequence.

Theorem 4.1 (Convergence to a weak solution). Let pD l p q lPN , pD l u q lPN , tpt l n q N l n"0 u lPN (where N l is the number of time steps of D l p ), be sequences of space time GDs assumed to satisfy the properties described in Section 3.1. Let φ m,min ą 0 and assume that, for each l P N, the gradient scheme (8a)-(8b) has a solution p α l P pX 0 D l p q N l `1, α P tnw, wu,

u l P pX 0 D l u q N l `1 such that (i) d f,D l u pt
, xq ě d 0 pxq for a.e. pt, xq P p0, T q ˆΓ, (ii) φ D l pt, xq ě φ m,min for a.e. pt, xq P p0, T q ˆΩ. Then, there exist pα P L 2 p0, T ; V 0 q, α P tnw, wu, and ū P L 8 p0, T ; U 0 q satisfying the weak formulation (6)-( 7) such that for α P tnw, wu and up to a subsequence

Π m D l p p α l á pα weakly in L 2 p0, T ; L 2 pΩqq, Π f D l p p α l á γ pα weakly in L 2 p0, T ; L 2 pΓqq, Π D l u u l á ū weakly-‹ in L 8 p0, T ; L 2 pΩq d q, φ D l á φm weakly-‹ in L 8 p0, T ; L 2 pΩqq, d f,D l u Ñ df in L 8 p0, T ; L p pΓqq for 2 ď p ă 4, Π m D l p S α m pp l c q Ñ S α m pp c q in L 2 p0, T ; L 2 pΩqq, Π f D l p S α f pp l c q Ñ S α f pγ pc q in L 2 p0, T ; L 2 pΓqq,
where φm " φ0 m `b divpū ´ū 0 q `1 M pp E m ´p E,0 m q, df " ´ ū , and pc " pnw ´p w .

Remark 4.2 (Discrete porosity and fracture aperture). As mentioned in the introduction, the assumptions that the discrete porosity and fracture aperture remain bounded below is a requirement coming from the model itself (which does not account for possible contact). It is not a fundamental restriction of the numerical framework and analysis.

We first present in Subsections 4.1 and 4.2 a sequence of intermediate results that will be useful for the proof of Theorem 4.1 detailed in Subsection 4.3.

Remark 4.3 (Incompressible limit for the solid matrix). The above convergence result also holds when 1{M " 0, i.e., in the incompressible limit for the grains of the solid matrix (M Ñ `8). Indeed, in this case, Lemma 4.4 below does not ensure L 8 pL 2 q-boundedness of the reconstructed matrix equivalent pressure. Nevertheless, L 2 pL 2 q-boundedness for this quantity (needed in the proof of the above theorem, cf. Subsection 4.3) can be readily inferred, based on the L 2 pL 2 q-boundedness of the reconstructed phase pressures (resulting from Lemma 4.4), the fact that reconstructed saturations are bounded, and the definition (5) of the capillary energy density.

Energy estimates

Using the phase pressures and velocity (time derivative of the displacement field) as test functions, the following a priori estimates can be inferred. Under hypotheses (H1)-(H6), there exists a real number C ą 0 depending on the data, the coercivity constants C Dp , C Du , and φ m,min , such that the following estimates hold:

}∇ m Dp p α } L 2 pp0,T qˆΩq ď C, }d 3 {2 f,Du ∇ f Dp p α } L 2 pp0,T qˆΓq ď C, }U m pΠ m Dp p c q} L 8 p0,T ;L 1 pΩqq ď C, }d 0 U f pΠ f Dp p c q} L 8 p0,T ;L 1 pΓqq ď C, }Π m Dp p E m } L 8 p0,T ;L 2 pΩqq ď C, }ε Du puq} L 8 p0,T ;L 2 pΩ,S d pRqqq ď C, }d f,Du } L 8 p0,T ;L 4 pΓqq ď C. (15)
Proof. For a piecewise constant function v on r0, T s with vptq " v n`1 for all t P pt n , t n`1 s, n P t0, . . . , N ´1u, and the initial value vp0q " v 0 , we define the piecewise constant function v such that vptq " v n for all t P pt n , t n`1 s. We notice the following expression for the discrete derivative of the product of two such functions:

δ t puvqptq " ûptqδ t vptq `vptqδ t uptq. (16)
In (8a), upon choosing ϕ α " p α we obtain T 1 `T2 `T3 `T4 " T 5 `T6 , with

T 1 " ż T 0 ż Ω δ t ´φD Π m Dp s α m ¯Πm Dp p α dxdt, T 2 " ż T 0 ż Ω η α m pΠ m Dp s α m qK m ∇ m Dp p α ¨∇m Dp p α dxdt, T 3 " ż T 0 ż Γ δ t ´df,Du Π f Dp s α f ¯Πf Dp p α dσpxqdt, T 4 " ż T 0 ż Γ η α f pΠ f Dp s α f q d 3 f,Du 12 ∇ f Dp p α ¨∇f Dp p α dσpxqdt, T 5 " ż T 0 ż Ω h α m Π m Dp p α dxdt, T 6 " ż T 0 ż Γ h α f Π f Dp p α dσpxqdt. (17)
First, we focus on the matrix and fracture accumulation terms T 1 and T 3 , respectively. Using ( 16) and the piecewise constant function reconstruction property of Π rt Dp , rt P tm, f u, we can write

δ t pφ D S α m pΠ m Dp p c qq " φD δ t S α m pΠ m Dp p c q `Sα m pΠ m Dp p c qδ t φ D , δ t pd f,Du S α f pΠ f Dp p c qq " df,Du δ t S α f pΠ f Dp p c q `Sα f pΠ f Dp p c qδ t d f,Du .
Summing on α P tw, nwu, we obtain

ÿ α pT 1 `T3 q " ÿ α ´ż T 0 ż Ω φD Π m Dp p α δ t S α m pΠ m Dp p c qdxdt `ż T 0 ż Ω S α m pΠ m Dp p c qΠ m Dp p α δ t φ D dxdt `ż T 0 ż Γ df,Du Π f Dp p α δ t S α f pΠ f Dp p c qdσpxqdt `ż T 0 ż Γ S α f pΠ f Dp p c qΠ f Dp p α δ t d f,Du dσpxqdt ¯. Now, for rt P tm, f u, ÿ α Π rt Dp p α δ t S α rt pΠ rt Dp p c q " Π rt Dp p c δ t S nw rt pΠ rt Dp p c q ě δ t U rt pΠ rt Dp p c q. (18)
Indeed, for n P t0, . . . , N ´1u, by the definition (5) of the capillary energy U rt and letting π rt c,n " Π rt Dp p c,n , we have

π rt c,n`1 pS nw rt pπ rt c,n`1 q ´Snw rt pπ rt c,n qq " U rt pπ rt c,n`1 q ´Urt pπ rt c,n q `ż π rt c,n`1 π rt c,n
pS nw rt pqq ´Snw rt pπ rt c,n qqdq ě U rt pπ rt c,n`1 q ´Urt pπ rt c,n q, where the last inequality holds since S nw rt is a non-decreasing function (see (H2)). Thus, we obtain

ÿ α pT 1 `T3 q ě ż T 0 ż Ω φD δ t U m pΠ m Dp p c qdxdt `ż T 0 ż Γ df,Du δ t U f pΠ f Dp p c qdσpxqdt `ÿ α ´ż T 0 ż Ω S α m pΠ m Dp p c qΠ m Dp p α δ t φ D dxdt `ż T 0 ż Γ S α f pΠ f Dp p c qΠ f Dp p α δ t d f,Du dσpxqdt ¯.
Applying again [START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF], we have

φD δ t U m pΠ m Dp p c q " δ t pφ D U m pΠ m Dp p c qq ´Um pΠ m Dp p c qδ t φ D , df,Du δ t U f pΠ f Dp p c q " δ t pd f,Du U f pΠ f Dp p c qq ´Uf pΠ f Dp p c qδ t d f,Du .
In the light of the closure equations (8c), this allows us to infer that

ÿ α pT 1 `T3 q ě ż T 0 ż Ω δ t pφ D U m pΠ m Dp p c qqdxdt `ż T 0 ż Γ δ t pd f,Du U f pΠ f Dp p c qqdσpxqdt `ż T 0 ż Ω 1 2M δ t ´Πm Dp p E m ¯2dxdt `ż T 0 ż Ω b Π m Dp p E m div Du pδ t uqdxdt ´ż T 0 ż Γ Π f Dp p E f δ t u Du dσpxqdt, (19)
where we have used the fact that

vδ t v ě δ t ˆv2 2
˙( 20
)
for v piecewise constant on r0, T s. Then, taking into account assumptions (H1)-(H6), there exists a real number C ą 0 depending only on the data such that

ÿ α pT 2 `T4 q ě C ´ż T 0 ż Ω ÿ α |∇ m Dp p α | 2 dxdt `ż T 0 ż Γ ÿ α |d 3 {2 f,Du ∇ f Dp p α | 2 dσpxqdt ¯. (21)
On the other hand, upon choosing v " δ t u in (8b), we get T 7 `T8 `T9 " T 10 , with

T 7 " ż T 0 ż Ω σ Du puq : ε Du pδ t uqdxdt, T 8 " ´ż T 0 ż Ω b Π m Dp p E m div Du pδ t uqdxdt T 9 " ż T 0 ż Γ Π f Dp p E f δ t u Du dσpxqdt, T 10 " ż T 0 ż Ω f ¨ΠDu pδ t uqdxdt. (22)
Using [START_REF] Brenner | Gradient Discretization of Hybrid Dimensional Darcy Flows in Fractured Porous Media with discontinuous pressure at matrix fracture interfaces[END_REF] and developing the definition of σ Du puq, we see that

T 7 ě ż T 0 ż Ω δ t ´1 2 σ Du puq : ε Du puq ¯dxdt, (23)
so that, all in all, taking into account that ř α pT 1 `T2 `T3 `T4 q `T7 `T8 `T9 " ř α pT 5 `T6 q `T10 and inequalities ( 19)-( 21)-( 23), we obtain the following estimate for the solutions of ( 8): there is a real number C ą 0 depending on the data such that

ż T 0 ż Ω δ t pφ D U m pΠ m Dp p c qq dxdt `ż T 0 ż Γ δ t pd f,Du U f pΠ f Dp p c qq dσpxqdt `ż T 0 ż Ω δ t ˆ1 2 σ Du puq : ε Du puq `1 2M pΠ m Dp p E m q 2 ˙dxdt `ÿ α ż T 0 ż Ω |∇ m Dp p α | 2 dxdt `ÿ α ż T 0 ż Γ |d 3 {2 f,Du ∇ f Dp p α | 2 dσpxqdt ď C ˜ż T 0 ż Ω f ¨δt Π Du u dxdt `ÿ α ż T 0 ż Ω h α m Π m Dp p α dxdt `ÿ α ż T 0 ż Γ h α f Π f Dp p α dσpxqdt ¸. (24)
Now, we have

ż T 0 ż Ω f ¨δt Π Du u dxdt " ż Ω f ¨pΠ Du upT q ´f ¨ΠDu up0qqdx ď C Du }f } L 2 pΩq p}ε Du puqpT q} L 2 pΩ,S d pRqq `}ε Du puqp0q} L 2 pΩ,S d pRqq q, ÿ α ´ż T 0 ż Ω h α m Π m Dp p α dxdt `ż T 0 ż Γ h α f Π f Dp p α dσpxqdt ď C Dp ÿ α p}h α m } L 2 pp0,T qˆΩq `}h α f } L 2 pp0,T qˆΓq qp}∇ m Dp p α } L 2 p0,T ;L 2 pΩqq `}d 3 {2 f,Du ∇ f Dp p α } L 2 p0,T ;L 2 pΓqq q,
where we have used the coercivity properties of the two gradient discretizations along with the Cauchy-Schwarz inequality and d 0 ď d f,Du . Using Young's inequality in the last two estimates as well as hypotheses (H1)-(H6) and (ii) in the lemma, it is then possible to infer from [START_REF] Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF] the existence of a real number C ą 0 depending on the data and on φ m,min such that

}U m pΠ m Dp p c qpT q} L 1 pΩq `}d 0 U f pΠ f Dp p c qpT q} L 1 pΩq `}pΠ m Dp p E m qpT q} 2 L 2 pΩq `}ε Du puqpT q} 2 L 2 pΩ,S d pRqq `ÿ α ´}∇ m Dp p α } 2 L 2 p0,T ;L 2 pΩqq `}d 3 {2 f,Du ∇ f Dp p α } 2 L 2 p0,T ;L 2 pΓqq ď C ´}f } 2 L 2 pΩq `ÿ α ´}h α m } 2 L 2 pp0,T qˆΩq `}h α f } 2 L 2 pp0,T qˆΓq }U m pΠ m Dp p c qp0q} L 1 pΩq `}d f,Du p0qU f pΠ f Dp p c qp0q} L 1 pΓq `}pΠ m Dp p E m qp0q} 2 L 2 pΩq `}pΠ f Dp p E f qp0q} 2 L 2 pΓq ¯.
The consistency property [START_REF] Bevillon | Stability and convergence analysis of partially coupled schemes for geomechanicalreservoir simulations[END_REF] shows that the terms above involving the discrete initial conditions are bounded and thus, together with the fact that T can be replaced by any t P p0, T s in the left-hand side, this inequality yields the a priori estimates (15) on p α , p c , p E m and u. The estimate on d f,Du follows from its definition and from the definition (13) of C Du .

Compactness properties

Throughout the analysis, we write a À b for a ď Cb with constant C depending only on the coercivity constants C Dp , C Du of the considered GDs, and on the physical parameters.

Estimates on time translates

Proposition 4.5. Let D p , D u , pt n q N n"0 be given space time GDs and φ m,min ą 0. It is assumed that the gradient scheme (8a)-(8b) has a solution p α P pX 0 Dp q N `1, α P tnw, wu, u P pX 0 Du q N `1 such that φ D pt, xq ě φ m,min for a.e. pt, xq P p0, T q ˆΩ and d f,Du pt, xq ě d 0 pxq for a.e. pt, xq P p0, T q ˆΓ. Let τ, τ 1 P p0, T q and, for s P p0, T s, denote by n s the natural number such that s P pt ns , t ns`1 s. For any ϕ P X 0 Dp , it holds

ˇˇxrφDΠ m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q, Π m Dp ϕy L 2 pΩq `xrd f,Du Π f Dp s α f spτ q ´rd f,Du Π f Dp s α f spτ 1 q, Π f Dp ϕy L 2 pΓq ˇÀ n τ 1 ÿ n"nτ `1 δt n`1 2 ´ξp1q,α,n`1 m }∇ m Dp ϕ} L 2 pΩq `ξp1q,α,n`1 f }∇ f Dp ϕ} L 8 pΓq `ξp2q,α,n`1 m }Π m Dp ϕ} L 2 pΩq `ξp2q,α,n`1 f }Π f Dp ϕ} L 2 pΓq ¯, (25)
with

N ´1 ÿ n"0 δt n`1 2 ´ξpjq,α,n`1 rt ¯2 À 1 for rt P tm, f u, j P t1, 2u, and
ξ p1q,α,n`1 m " }∇ m Dp p α n`1 } L 2 pΩq and ξ p1q,α,n`1 f " }pd n`1 f,Du q 3 {2 ∇ f Dp p α n`1 } L 2 pΓq }d n`1 f,Du } 3 {2 L 4 pΓq , ξ p2q,α,n`1 m " › › › 1 δt n`1 2 ż tn`1 tn h α m pt, ¨qdt › › › L 2 pΩq ξ p2q,α,n`1 f " › › › 1 δt n`1 2 ż tn`1 tn h α f pt, ¨qdt › › › L 2 pΓq
.

Proof. For any ϕ P X 0 Dp , writing the difference of piecewise-constant functions at times τ and τ 1 as the sum of their jumps between these two times, one has

ˇˇxrφDΠ m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q, Π m Dp ϕy L 2 pΩq `xrd f,Du Π f Dp s α f spτ q ´rd f,Du Π f Dp s α f spτ 1 q, Π f Dp ϕy L 2 pΓq ˇď n τ 1 ÿ n"nτ `1 δt n`1 2 ˇˇxδtrφDΠ m Dp s α m spt n`1 q, Π m Dp ϕy L 2 pΩq `xδ t rd f,Du Π f Dp s α f spt n`1 q, Π f Dp ϕy L 2 pΓq ˇˇ. (26)
From the gradient scheme discrete variational equation (8a), we deduce that

ˇˇxδtrφDΠ m Dp s α m spt n`1 q, Π m Dp ϕy L 2 pΩq `xδ t rd f,Du Π f Dp s α f spt n`1 q, Π f Dp ϕy L 2 pΓq ˇÀ }∇ m Dp p α n`1 } L 2 pΩq }∇ m Dp ϕ} L 2 pΩq `}pd n`1 f,Du q 3 {2 ∇ f Dp p α n`1 } L 2 pΓq }pd n`1 f,Du q 3 {2 ∇ f Dp ϕ} L 2 pΓq `› › › 1 δt n`1 2 ż tn`1 tn h α m pt, ¨qdt › › › L 2 pΩq }Π m Dp ϕ} L 2 pΩq `› › › 1 δt n`1 2 ż tn`1 tn h α f pt, ¨qdt › › › L 2 pΓq }Π f Dp ϕ} L 2 pΓq À ξ p1q,α,n`1 m }∇ m Dp ϕ} L 2 pΩq `ξp1q,α,n`1 f }∇ f Dp ϕ} L 8 pΓq `ξp2q,α,n`1 m }Π m Dp ϕ} L 2 pΩq `ξp2q,α,n`1 f }Π f Dp ϕ} L 2 pΓq , (27)
where the term }pd n`1 f,Du q 3 {2 ∇ f Dp ϕ} L 2 pΓq has been estimated using the generalized Hölder inequality with exponents p8, 8{3q, which satisfy 1 8 `3 8 " 1 2 . Hence the result follows from ( 26), [START_REF] Antonio | The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications[END_REF], the a priori estimates of Lemma 4.4, and from the assumptions h α m P L 2 pp0, T q ˆΩq, h α f P L 2 pp0, T q ˆΓq. Remark 4.6. Summing the estimate (25) on α P tnw, wu we obtain the following time translate estimates on φ D and d f,Du : ˇˇxφDpτ q ´φD pτ 1 q, Π m Dp ϕy L 2 pΩq `xd f,Du pτ q ´df,Du pτ 1 q, Π f Dp ϕy L 2 pΓq ˇÀ ÿ αPtnw,wu

n τ 1 ÿ n"nτ `1 δt n`1 2 ´ξp1q,α,n`1 m }∇ m Dp ϕ} L 2 pΩq `ξp1q,α,n`1 f }∇ f Dp ϕ} L 8 pΓq `ξp2q,α,n`1 m }Π m Dp ϕ} L 2 pΩq `ξp2q,α,n`1 f }Π f Dp ϕ} L 2 pΓq ¯. (28)

Compactness properties of Π m

Dp s α m Proposition 4.7. Let pD l p q lPN , pD l u q lPN , tpt l n q N l n"0 u lPN be sequences of space time GDs assumed to satisfy the properties described in Section 3.1. Let φ m,min ą 0 and assume that, for each l P N, the gradient scheme (8a)-(8b) has a solution p α l P pX 0 D l p q N l `1, α P tnw, wu, u l P pX 0 D l u q N l `1 such that (i) d f,D l u pt, xq ě d 0 pxq for a.e. pt, xq P p0, T q ˆΓ, (ii) φ D l pt, xq ě φ m,min for a.e. pt, xq P p0, T q ˆΩ. Then, the sequence pΠ m Dp s α,l m q lPN , with s α,l m " S α m pp l c q, is relatively compact in L 2 pp0, T q ˆΩq.

Proof. Let K be a fixed compact set of ΩzΓ and let us consider cut-off functions ψ l as defined in the cut-off property of the sequence of spatial GDs pD l p q lPN . The superscript l P N will be dropped in the proof, and assumed to be large enough. All hidden constants in the following estimates are independent of l. Using that φ D pt, xq ě φ m,min for a.e. pt, xq P p0, T q ˆΩ, the properties of the cut-off functions, and noting that Π m Dp s α,l m " S α m pΠ m Dp p l c q P r0, 1s, we obtain

ż T 0 }Π m Dp s α m p¨`τ, ¨q ´Πm Dp s α m } 2 L 2 pKq dt À τ `ż T ´τ 0 ż Ω pΠ m Dp ψq φ D ´Πm Dp s α m p¨`τ, ¨q ´Πm Dp s α m ¯2 dxdt " τ `T1 `T2 ,
where

T 1 " ż T ´τ 0 ˇˇxrφDΠ m Dp s α m spt `τ q ´rφ D Π m Dp s α m sptq, Π m Dp ζ α m ptqy L 2 pΩq ˇˇdt, T 2 " ż T ´τ 0 ˇˇxφDpt `τ q ´φD ptq, Π m Dp χ α m ptqy L 2 pΩq ˇˇdt,
with ζ α m ptq " ´sα m pt `τ q ´sα m ptq ¯ψ and χ α m ptq " ζ α m ptq s α m pt `τ q. From the cut-off property it results that

Π f Dp ζ α m " 0 and ∇ f Dp ζ α m " 0.
Then, in view of the estimates (25), we have

T 1 À ż T ´τ 0 n pt`τ q ÿ n"nt`1 δt n`1 2 ´ξp1q,α,n`1 m }∇ m Dp ζ α m ptq} L 2 pΩq `ξp2q,α,n`1 m }Π m Dp ζ α m ptq} L 2 pΩq ¯dt À ż T ´τ 0 n pt`τ q ÿ n"nt`1 δt n`1 2 ´pξ p1q,α,n`1 m q 2 `pξ p2q,α,n`1 m q 2 `}∇ m Dp ζ α m ptq} 2 L 2 pΩq `}Π m Dp ζ α m ptq} 2 L 2 pΩq ¯dt.
From Proposition 4.5, we have

N ´1 ÿ n"0 δt n`1 2 ´pξ p1q,α,n`1 m q 2 `pξ p2q,α,n`1 m q 2 ¯À 1.
Using the a priori estimates of Lemma 4.4, h α m P L 2 pp0, T q ˆΩq, the Lipschitz property of S α m , the chain rule and product rule estimates on the sequence of GDs pD l p q lPN , and the cut-off property, we obtain that

ż T ´τ 0 ´}∇ m Dp ζ α m ptq} 2 L 2 pΩq `}Π m Dp ζ α m ptq} 2 L 2 pΩq ¯dt À 1.
We deduce from [5, Lemma 4.1] that T 1 À τ `∆t with a hidden constant depending on K but independent of l.

Similarly, using the time translate estimate [START_REF] Droniou | Uniform-in-time convergence of numerical methods for non-linear degenerate parabolic equations[END_REF], one shows that T 2 À τ `∆t, which provides the time translates estimates on Π m Dp s α m in L 2 p0, T ; L 2 pKqq.

The space translates estimates for Π m Dp s α m in L 2 p0, T ; L 2 pKqq derive from the a priori estimates of Lemma 4.4, the Lipschitz properties of S α m and from the compactness property of the sequence of spatial GDs pD l p q lPN (cf. Remark 3.6). Combined with the time translate estimates, the Fréchet-Kolmogorov theorem implies that Π m Dp s α m is relatively compact in L 2 p0, T ; L 2 pKqq for any compact set K of ΩzΓ. Since Π m Dp s α m P r0, 1s, it results that Π m Dp s α m is relatively compact in L 2 pp0, T q ˆΩq. Under the assumptions of Proposition 4.7, the sequences pφ D l q lPN and pφ D l Π m Dp s α,l m q lPN , with s α,l m " S α m pp l c q, converge up to a subsequence uniformly in time weakly in L 2 pΩq.

Proof. Let K be a fixed compact set of ΩzΓ and let ψ l be cut-off functions for this compact set, as defined in the cut-off property of pD l p q lPN . The superscript l P N will be dropped when not required for the clarity of the proof, and assumed to be large enough.

For w P V 0 we let P Dp w P X 0 Dp be the element that realizes the minimum in S Dp pwq, so that }∇ m Dp P Dp w ´∇w} L 2 pΩq `}∇ f Dp P Dp w ´∇τ γw} L r pΓq `}Π m Dp P Dp w ´w} L 2 pΩq `}Π f Dp P Dp w ´γw} L r pΓq " S Dp pwq. (29)

Let ϕ P C 8 c pΩq and set ϕ " P Dp ϕ. It results from the cut-off property that Π f Dp pψϕq " 0 and ∇ f Dp pψϕq " 0. Using the GD consistency property of pD l p q lPN and ( 29), we see that }∇ m Dp pψϕq} L 2 pΩq and }Π m Dp pψϕq} L 2 pΩq are bounded by constants depending on K and ϕ but independent of l. Then, from Proposition 4.5, we have with hidden constants independent of l but possibly depending on K and ϕ, that

ˇˇxΠ m Dp ψ ´rφ D Π m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q ¯, Π m Dp ϕy L 2 pΩq ˇ" ˇˇxrφDΠ m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q, Π m Dp pψϕqy L 2 pΩq ˇÀ n τ 1 ÿ n"nτ `1 δt n`1 2 ´ξp1q,α,n`1 m }∇ m Dp pψϕq} L 2 pΩq `ξp2q,α,n`1 m }Π m Dp pψϕq} L 2 pΩq À ˜nτ 1 ÿ n"nτ `1 δt n`1 2 ˆ´ξ p1q,α,n`1 m ¯2 `´ξ p2q,α,n`1 m ¯2˙¸1 2 ˜nτ 1 ÿ n"nτ `1 δt n`1 2 ¸1 2 À |τ ´τ 1 | 1 2 `∆t 1 2 .
Since Π m Dp s α m P r0, 1s, φ D is bounded in L 8 p0, T ; L 2 pΩqq (see (8c) and ( 15)), and Π m Dp ψ is uniformly bounded, one has

ˇˇxΠ m Dp ψ ´rφ D Π m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q ¯, ϕy L 2 pΩq ˇˇÀ |τ ´τ 1 | 1 2 `∆t 1 2 `ωDp , (30)
with ω Dp " }ϕ ´Πm Dp ϕ} L 2 pΩq a consistency error term such that lim lÑ`8 ω D l p " 0. It follows from the discontinuous Ascoli-Arzelà theorem [START_REF] Droniou | The Gradient Discretisation Method[END_REF]Theorem C.11] that (up to a subsequence) the sequence pΠ m Dp ψqφ D pΠ m Dp s α m q " φ D Π m Dp ps α m ψq converges uniformly in time weakly in L 2 pΩq.

Let us now take w P C 8 c pΩzΓq and let K be the support of w. For l large enough, by definition of ψ l we have

`φD l Π m D l p s α,l m ˘|K " φ D l Π m D l p pψ l s α,l m q. Hence, xφ D l Π m D l
p s α,l m , wy L 2 pΩq converges uniformly with respect to t P r0, T s.

Since pφ D l Π m D l p s α,l m q lPN is bounded in L 8 p0, T ; L 2 pΩqq, the density of C 8 c pΩzΓq in L 2 pΩq shows that the convergence (31) is valid for any w P L 2 pΩq, which concludes the proof that the sequence φ D l Π m Proposition 4.9. Under the assumptions of Proposition 4.7, the sequences pd f,D l u q lPN and pd f,D l u Π f Dp s α,l f q lPN , with s α,l f " S α f pp l c q, converge up to a subsequence uniformly in time weakly in L 2 pΓq.

Proof. Let K be a fixed compact set of ΩzΓ and let us consider cut-off functions ψ l as defined in the cut-off property of pD l p q lPN . In the following, the superscript l P N (assumed to be large enough) is dropped when not required for the clarity of the proof, and the hidden constants are independent of l. Let ϕ P C 8 c pΩq and set ϕ " P Dp ϕ, with P Dp characterised by [START_REF] Droniou | Uniform-in-time convergence of numerical methods for non-linear degenerate parabolic equations[END_REF]. From Proposition 4.5 we have

ˇˇxrdf,D u Π f Dp s α f spτ q ´rd f,Du Π f Dp s α f spτ 1 q, Π f Dp ϕy L 2 pΓq ˇÀ ˇˇx ´rφ D Π m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q ¯, Π m Dp ϕy L 2 pΩq ˇm ax ´}∇ m Dp ϕ} L 2 pΩq , }∇ f Dp ϕ} L 8 pΓq , }Π m Dp ϕ} L 2 pΩq , }Π f Dp ϕ} L 2 pΓq nτ 1 ÿ n"nτ `1 δt n`1 2 ˆ´ξ p1q,α,n`1 m ¯2 `´ξ p1q,α,n`1 f ¯2 `´ξ p2q,α,n`1 m ¯2 `´ξ p2q,α,n`1 f ¯2˙¸1 2 ˆ˜n τ 1 ÿ n"nτ `1 δt n`1 2 ¸1 2 À ´|τ ´τ 1 | 1 2 `∆t 1 2 ¯`ˇˇˇx´r φ D Π m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q ¯, Π m Dp ϕy L 2 pΩq ˇˇ.
Since φ D Π m Dp s α m is bounded in L 8 p0, T ; L 2 pΩqq (see the proof of Proposition 4.8), we have

ˇˇx ´rφ D Π m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q ¯, Π m Dp ϕy L 2 pΩq ˇÀ } φ ´Πm Dp ϕ} L 2 pΩq `ˇˇx rφ D Π m Dp s α m spτ q ´rφ D Π m Dp s α m spτ 1 q, φy L 2 pΩq ˇǎ nd ˇˇxrdf,D u Π f Dp s α f spτ q ´rd f,Du Π f Dp s α f spτ 1 q, φ ´Πf Dp ϕy L 2 pΓq ˇˇÀ }d f,Du } L 8 p0,T ;L 2 pΓqq } φ ´Πf Dp ϕ} L 2 pΓq .
Using the a priori estimates of Lemma 4.4, and Proposition 4.8 stating the uniform-in-time L 2 pΩq-weak convergence of φ D Π m Dp s α m (which implies the equi-continuity of the functions τ Þ Ñ xrφ D Π m Dp s α m spτ q, φy L 2 pΩq ), we deduce that ˇˇxrdf,D u Π f Dp s α f spτ q ´rd f,Du Π f Dp s α f spτ 1 q, φy L 2 pΓq ˇˇÀ ωp|τ ´τ 1 |q `∆t Under the assumptions of Proposition 4.7, the sequence pd f,D l u q lPN converges up to a subsequence in L 8 p0, T ; L p pΓqq for all 2 ď p ă 4, and the sequences pd f,D l u Π f Dp s α,l f q lPN and pΠ f Dp s α,l f q lPN , with s α,l f " S α f pp l c q, converge up to a subsequence in L 4 p0, T ; L 2 pΓqq.

Proof. By the characterization in Remark 3.8 of the compactness of pD l u q lPN and the estimate on ε Du puq in Lemma 4.4, we have, for all i P I, all η i tangent to Γ i , a.e. t P p0, T q and all s ă 4,

› › d f,D l u pt, ¨`η i q ´df,D l u pt, ¨q› › L s pΓiq ď T D l u ,s p0, ηq}ε Du puqpt, ¨q} L 2 pΩ,S d pRqq À T D l u ,s p0, ηq,
where η " p0, . . . , 0, η i , 0, . . . , 0q and d f,D l u has been extended by 0 in the hyperplane spanned by Γ i . Together with the uniform-in-time L 2 pΓq-weak convergence of d f,D l u from Proposition 4.9, this shows that we can apply Lemma A.2 to d f,D l u with p " `8 and get the convergence of this sequence in L 8 p0, T ; L 2 pΓqq. Since, from the a priori estimates of Lemma 4.4, this sequence d f,D l u is bounded in L 8 p0, T ; L 4 pΓqq, it follows that it converges in L 8 p0, T ; L q pΓqq for all 2 ď q ă 4.

For any compact set K f Ă Γ that is disjoint from the intersections pΓ i X Γ j q i "j , using that Π f Dp s α f P r0, 1s, that }d f,Du pt, ¨q} L 4 pΓq is uniformly bounded in t, and the Lipschitz properties of S α f , it follows that, for all i P I and η i tangent to Γ i small enough,

}rd f,Du Π f Dp s α f spt, ¨`η i q ´rd f,Du Π f Dp s α f spt, ¨q} L 2 pK f XΓiq ď }d f,Du pt, ¨`η i q ´df,Du pt, ¨q} L 2 pK f XΓiq `}Π f Dp s α f pt, ¨`η i q ´Πf Dp s α f pt, ¨q} L 4 pK f XΓiq }d f,Du pt, ¨q} L 4 pK f XΓiq À }d f,Du pt, ¨`η i q ´df,Du pt, ¨q} L 2 pK f XΓiq `}Π f Dp s α f pt, ¨`η i q ´Πf Dp s α f pt, ¨q} 1 2 L 2 pK f XΓiq À }d f,Du pt, ¨`η i q ´df,Du pt, ¨q} L 2 pK f XΓiq `}Π f Dp p c pt, ¨`η i q ´Πf Dp p c pt, ¨q} 1 2 L 2 pK f XΓiq .
From the compactness properties of pD l u q lPN and pD l p q lPN (see Remarks 3.6 and 3.8) it results that

ÿ iPI › › › sup |ηi|ďδ }rd f,Du Π f Dp s α f sp¨, ¨`η i q ´rd f,Du Π f Dp s α f sp¨, ¨q} L 2 pK f XΓiq › › › L 4 p0,T q À T K f pδq ´}ε Du puq} L 8 p0,T ;L 2 pΩqq `ÿ αPtnw,wu p}d 3 {2 0 ∇ f Dp p α } L 2 p0,T ;L 2 pΓqq `}∇ m Dp p α } L 2 p0,T ;L 2 pΩqq q with lim δÑ0 T K f pδq " 0.
From the a priori estimates of Lemma 4.4, and the uniform-in-time L 2 pΓq-weak convergence of d f,Du s α f of Proposition 4.9, it follows from Lemma A.2 that d f,Du Π f Dp s α f converges up to a subsequence in L 4 p0, T ; L 2 pK f qq.

From the assumption d f,Du pt, xq ě d 0 pxq, d f,Du is bounded below by a strictly positive constant on K f . Writing that

Π f Dp s α f " 1 d f,Du pd f,Du Π f Dp s α f q, it follows that Π f Dp s α f converges in L 4 p0, T ; L 2 pK f qq.
Since this is true for any K f compact in Γ that does not touch the fractures intersections, and since Π f Dp s α f P r0, 1s, we deduce that Π f Dp s α f converges in L 4 p0, T ; L 2 pΓqq.

Convergence to a weak solution

Proof of Theorem 4.1. The superscript l will be dropped in the proof, and all convergences are up to appropriate subsequences. From Lemma 4.4 and Proposition 4.10, there exist df P L 8 p0, T ; L 4 pΓqq and sα

f P L 8 pp0, T q ˆΓq such that d f,Du Ñ df in L 8 p0, T ; L p pΓqq, 2 ď p ă 4, Π f Dp S α f pp c q Ñ sα f in L 4 p0, T ; L 2 pΓqq. (32)
From Proposition 4.7, there exists sα m P L 8 pp0, T q ˆΩq such that

Π m Dp S α m pp c q Ñ sα m in L 2 p0, T ; L 2 pΩqq. ( 33
)
The identification of the limit [18, Lemma 5.5], resulting from the limit-conformity property, can easily be adapted to our definition of V 0 , with weight d

3 {2
0 and the use in the definition of limit-conformity of fracture flux functions that are compactly supported away from the tips. Using this lemma and the a priori estimates of Lemma 4.4, we obtain pα P L 2 p0, T ; V 0 q and g α f P L 2 p0, T ; L 2 pΓq d´1 q, such that the following weak limits hold

Π m Dp p α á pα in L 2 p0, T ; L 2 pΩqq weak, Π f Dp p α á γ pα in L 2 p0, T ; L 2 pΓqq weak, ∇ m Dp p α á ∇p α in L 2 p0, T ; L 2 pΩq d q weak, d 3 {2 0 ∇ f Dp p α á d 3 {2 0 ∇ τ γ pα in L 2 p0, T ; L 2 pΓq d´1 q weak, d 3 {2 f,Du ∇ f Dp p α á g α f in L 2 p0, T ; L 2 pΓq d´1 q weak. ( 34
)
Let ϕ P C 0 c pp0, T q ˆΓq d´1 whose support is contained in p0, T q ˆK, with K compact set not containing the tips of Γ. We have

ż T 0 ż Γ d 3 {2 f,Du ∇ f Dp p α ¨ϕ dσpxqdt Ñ ż T 0 ż Γ g α f ¨ϕ dσpxqdt.
On the other hand, it results from [START_REF] Eymard | Finite volume methods[END_REF] and the fact that d 0 is bounded away from 0 on K (because d 0 is continuous and does not vanish outside the tips of Γ) that ∇ f Dp p α á ∇ τ γ pα in L 2 p0, T ; L 2 pKq d´1 q. Combined with the convergence

d 3 {2
f,Du ϕ Ñ p df q 3 {2 ϕ in L 8 p0, T ; L 2 pΓq d´1 q given by ( 32), we infer that

ż T 0 ż Γ d 3 {2 f,Du ∇ f Dp p α ¨ϕ dσpxqdt Ñ ż T 0 ż Γ p df q 3 {2 ∇ τ γ pα ¨ϕ dσpxqdt.
This shows that g α f " p df q 3 {2 ∇ τ γ pα on p0, T q ˆΓ.

Combining the strong convergence of Π m Dp S α m pp c q " S α m pΠ m Dp p c q (resp. of Π f Dp S α f pp c q " S α f pΠ f Dp p c q), the weak convergence of Π m Dp p c (resp. Π f Dp p c ), and the monotonicity of S α m (resp. S α f q, it results from the Minty trick (see e.g. [START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF]Lemma 2.6]) that sα m " S α m pp c q (resp. sα f " S α f pγ pc q) with pc " pnw ´p w . From the a priori estimates of Lemma 4.4 and the limit-conformity property of the sequence of GDs pD l u q lPN (see Lemma A.3), there exists ū P L 8 p0, T ; U 0 q, such that

Π Du u á ū in L 8 p0, T ; L 2 pΩq d q weak ‹, ε Du puq á εpūq in L 8 p0, T ; L 2 pΩ, S d pRqqq weak ‹, div Du u á divpūq in L 8 p0, T ; L 2 pΩqq weak ‹, d f,Du " ´ u Du á ´ ū in L 8 p0, T ; L 2 pΓqq weak ‹, (35)
from which we deduce that df " ´ ū and that σ Du puq converges to σpūq in L 8 p0, T ; L 2 pΩ, S d pRqqq weak ‹.

From the a priori estimates and the closure equations (8c), there exist φm P L 8 p0, T ; L 2 pΩq and pE m P L 8 p0, T ; L 2 pΩq such that

φ D á φm in L 8 p0, T ; L 2 pΩqq weak ‹, Π m Dp p E m á pE m in L 8 p0, T ; L 2 pΩqq weak ‹. ( 36
)
Since 0 ď U rt pzq " ş p 0 zpS nw rt q 1 pzqdz ď 2|p| for rt P tm, f u, it results from the a priori estimates of Lemma 4.4 that there exist pE f P L 2 p0, T ; L 2 pΓqq, Ūf P L 2 p0, T ; L 2 pΓqq and Ūm P L 2 p0, T ; L 2 pΩqq such that

Π f Dp p E f á pE f in L 2 p0, T ; L 2 pΓqq weak, Π f Dp U f pp c q á Ūf in L 2 p0, T ; L 2 pΓqq weak, Π m Dp U m pp c q á Ūm in L 2 p0, T ; L 2 pΩqq weak. ( 37
)
For rt P tnw, wu, it is shown in [START_REF] Droniou | Numerical analysis of a two-phase flow discrete fracture model[END_REF], following ideas from [START_REF] Droniou | Uniform-in-time convergence of numerical methods for non-linear degenerate parabolic equations[END_REF], that U rt ppq " B rt pS nw rt ppqq where B rt : r0, 1s Þ Ñ R is a convex lower semi-continuous function with finite limits at s " 0 and s " 1 (note that B rt is therefore actually continuous on r0, 1s). Since Π m Dp s nw m converges strongly in L 2 pp0, T q ˆΩq to S nw m pp c q, it converges a.e. in p0, T q ˆΩ. It results that B m pΠ m Dp s nw m q converges a.e. in p0, T q ˆΩ to B m pS nw m pp c qq, and hence that Ūm " B m pS nw m pp c qq " U m pp c q. Similarly, Ūf " B f pS nw f pγ pc qq " U f pγ pc q. We deduce that pE m " ÿ αPtnw,wu pα S α m pp c q ´Um pp c q and pE f " ÿ αPtnw,wu γ pα S α f pγ pc q ´Uf pγ pc q.

Using the estimate |U rt pp 2 q ´Urt pp 1 q| " ˇˇˇż p2 p1 zpS nw rt q 1 pzqdz ˇˇˇď |p 2 ´p1 | `|p 2 S nw rt pp 2 q ´p1 S nw rt pp 1 q|, the Lipschitz property of S nw rt , pα 0 P V 0 X L 8 pΩq, γ pα 0 P L 8 pΓq, α P tnw, wu, and the consistency of the sequence of GDs pD l p q lPN , we deduce that

Π m Dp p E,0 m Ñ pE,0 m in L 2 pΩq, Π f Dp p E,0 f Ñ pE,0 f in L 2 pΓq. (38)
Then, from Proposition A.4 it holds that div Du pu 0 q Ñ divpū 0 q in L 2 pΩq,

u 0 Du Ñ ū0 " ´d 0 f in L 2 pΓq. (39)
It results from [START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF], [START_REF] Eymard | Discontinuous Galerkin gradient discretisations for the approximation of secondorder differential operators in divergence form[END_REF], [START_REF] Garipov | Discrete fracture model for coupled flow and geomechanics[END_REF] and [START_REF] Giovanardi | Unfitted FEM for modelling the interaction of multiple fractures in a poroelastic medium[END_REF] and the definition of

φ D that φm " φ0 m `b divpū ´ū 0 q `1 M pp E m ´p E,0 m q.
Let us now prove that the functions pα , α P tnw, wu, and ū satisfy the variational formulation ( 6)-(7) by passing to the limit in the gradient scheme [START_REF] Ayuso De Dios | The nonconforming virtual element method[END_REF].

For θ P C 8 c pr0, T qq and ψ P C 8 c pΩq let us set, with P Dp characterised by [START_REF] Droniou | Uniform-in-time convergence of numerical methods for non-linear degenerate parabolic equations[END_REF], ϕ " pϕ 1 , . . . , ϕ N q P pX 0 Dp q N with ϕ i " θpt i´1 qpP Dp ψq.

From the consistency properties of pD l p q lPN with given r ą 8, we deduce that

Π m Dp P Dp ψ Ñ ψ in L 2 pΩq, Π f Dp P Dp ψ Ñ γψ in L 2 pΓq, Π m Dp ϕ Ñ θψ in L 8 p0, T ; L 2 pΩqq, Π f Dp ϕ Ñ θγψ in L 8 p0, T ; L 2 pΓqq, ∇ m Dp ϕ Ñ θ∇ψ in L 8 p0, T ; L 2 pΩq d q, ∇ f Dp ϕ Ñ θ∇ τ γψ in L 8 p0, T ; L r pΓq d´1 q. (40)
Setting

T 1 " ż T 0 ż Ω δ t ´φD Π m Dp s α m ¯Πm Dp ϕ dxdt T 2 " ż T 0 ż Ω η α m pΠ m Dp s α m qK m ∇ m Dp p α ¨∇m Dp ϕ dxdt T 3 " ż T 0 ż Γ δ t ´df,Du Π f Dp s α f ¯Πf Dp ϕ dσpxqdt T 4 " ż T 0 ż Γ η α f pΠ f Dp s α f q d 3 f,Du 12 ∇ f Dp p α ¨∇f Dp ϕ dσpxqdt T 5 " ż T 0 ż Ω h α m Π m Dp ϕ dxdt `ż T 0 ż Γ h α f Π f Dp ϕ dσpxqdt,
the gradient scheme variational formulation (8a) states that T 1 `T2 `T3 `T4 " T 5 .

For ω P C 8 c pr0, T qq and a smooth function w : ΩzΓ Ñ R d vanishing on BΩ and admitting finite limits on each side of Γ, let us set v " pv 1 , . . . , v N q P pX 0 Du q N with v i " ωpt i´1 qpP Du wq where P Du w realizes the minimum in the definition ( 14) of S Du pwq. From the consistency properties of pD l u q lPN , we deduce that

Π Du v Ñ ωψ in L 8 p0, T ; L 2 pΩq d q, ε Du pvq Ñ ωεpwq in L 8 p0, T ; L 2 pΩ, S d pRqqq, v Du Ñ ω w in L 8 p0, T ; L 2 pΓqq. (41)
Setting

T 6 " ż T 0 ż Ω ´σDu puq : ε Du pvq ´bpΠ m Dp p E m qdiv Du pvq ¯dxdt, T 7 " ż T 0 ż Γ pΠ f Dp p E f q v Du dσpxqdt, T 8 " ż T 0 ż Ω f ¨ΠDu v dxdt.
the gradient scheme variational formulation (8b) states that

T 6 `T7 " T 8 .
Using a discrete integration by part [30, Section D.1.7], we have T 1 " T 11 `T12 with

T 11 " ´ż T 0 ż Ω φ D pΠ m Dp s α m qpΠ m Dp P Dp ψqθ 1 ptq dxdt, T 12 " ´żΩ pΠ m Dp I m Dp φ0 qpΠ m Dp S α m pI Dp pα 0 qqpΠ m Dp P Dp ψqθp0q dx.
Using [START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF] and [START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF], and that Π m Dp s α m P r0, 1s converges to S α m pp c q a.e. in p0, T q ˆΩ (this follows from (33)), it holds that

T 11 Ñ ´ż T 0 ż Ω φm S α m pp c qψθ 1 ptq dxdt.
Using [START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF], that Π m Dp I m Dp φ0 converges in L 2 pΩq to φ0 and that Π m Dp S α m pP Dp pα 0 q P r0, 1s converges a.e. in Ω to S α m pp α 0 q, we deduce that T 12 Ñ ´żΩ φ0 S α m pp α 0 qψθp0q dx.

Writing T 3 " T 31 `T32 with

T 31 " ´ż T 0 ż Γ d f,Du pΠ f Dp s α f qpΠ f Dp P Dp ψqθ 1 ptq dσpxqdt, T 32 " ż Γ u 0 Du pΠ f Dp S α f pI Dp pα 0 qqpΠ f Dp P Dp ψqθp0q dσpxq,
we obtain, using similar arguments and (39), that

T 31 Ñ ´ż T 0 ż Γ df S α f pγ pc qγψθ 1 ptq dσpxqdt,
and

T 32 Ñ ´żΓ d0 f S α f pγ pα 0 qγψθp0q dσpxq.
Using that 0 ď η α m pΠ m Dp s α m q ď η α m,max , the continuity of η α m , the convergence of Π m Dp s α m a.e. in p0, T q ˆΩ to S α m pp c q, (34) and [START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF], it holds that

T 2 Ñ ż T 0 ż Ω η α m pS α m pp c qqK m ∇p α ¨θ∇ψ dxdt.
The convergence

T 4 Ñ ż T 0 ż Γ η α f pS α f pγ pc qq d 3 f 12 ∇ τ γ pα ¨θ∇ τ γψ dσpxqdt is established using 0 ď η α f pΠ f Dp s α f q ď η α f,max
, the continuity of η α f , the convergence of Π f Dp s α f a.e. in p0, T q ˆΓ to S α f pγ pc q, combined with the weak convergence of

d 3 {2 f,Du ∇ f Dp p α to d3 {2 f ∇ τ γ pα in L 2 pp0, T qˆΓq d´1 , the strong convergence of d 3 {2
f,Du to d3 {2 f in L s pp0, T q ˆΓq for all 2 ď s ă 8 3 (resulting from ( 32)), and the strong convergence (40) of ∇ f Dp ϕ to θ∇ τ γψ in L 8 p0, T ; L r pΓqq with r ą 8.

The convergence

T 5 Ñ ż T 0 ż Ω h α m θψ dxdt `ż T 0 ż Γ h α f θpγψq dσpxqdt
is readily obtained from [START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF]. The following convergences of T 6 , T 7 , T 8

T 6 Ñ ż T 0 ż Ω ´σpūq : εpwqω ´bp E m divpwqω ¯dxdt, T 7 Ñ ż T 0 ż Γ pE f w ω dσpxqdt, T 8 Ñ ż T 0 ż Ω f ¨wω dxdt
classically result from the strong convergences (41) combined with the weak convergences [START_REF] Eymard | Discontinuous Galerkin gradient discretisations for the approximation of secondorder differential operators in divergence form[END_REF].

Using the above limits in T 1 `T2 `T3 `T4 " T 5 and T 6 `T7 " T 8 concludes the proof that pα , α P tnw, wu, and ū satisfy the variational formulation ( 6)-( 7).

Two-dimensional numerical example

The objective of this section is to numerically investigate the convergence of the discrete solutions on a simple geometrical configuration based on a cross-shaped fracture network. We refer to [START_REF] Bonaldi | Gradient discretization of two-phase poro-mechanical models with discontinuous pressures at matrix fracture interfaces[END_REF] for the presentation of a more advanced application to the desaturation by suction at the interface between a ventilation tunnel and a Callovo-Oxfordian argilite fractured storage rock.

Setting

Let us consider the system (3)-(4) in the square domain Ω " p0, Lq 2 , with L " 100 m, lying in the xy-plane and containing a cross-shaped fracture network Γ made up of four fractures (cf. Figure 3), each of length L 8 , aligned with the coordinate axes and intersecting at the center of the domain p L 2 , L 2 q. More precisely, the fracture network is defined as follows: Γ "

Ť 4 i"1 Γ i , where Γ 1 " p 3 8 L, L 2 q ˆt L 2 u, Γ 2 " p L 2 , 5 8 Lq ˆt L 2 u, Γ 3 " t L 2 u ˆp 3 8 L, L 2
q, and Γ 4 " t L 2 u ˆp L 2 , 5 8 Lq. The data set employed is inspired by [START_REF] Shu | A dual-porosity model for two-phase flow in deforming porous media[END_REF]Section 5.1.2]. The matrix and fracture network have the following mobility laws: where µ w " 10 ´3 Pa¨s and µ nw " 1.851¨10 ´5 Pa¨s are the dynamic viscosities of the wetting and non-wetting phases, respectively. Notice that η α m and η α f do not satisfy the assumptions of our analysis, as they are not bounded below by a strictly positive number; these choices are however physically relevant, and as the test shows, do not seem to impair the convergence of the numerical scheme. A non-degenerate regularization of these mobilities is also investigated below. The function yielding the saturation in both rock types in terms of the capillary pressure is provided by Corey's law:

η α m ps α q " ps α q 2 µ α , η α f ps α q " s α µ α , α P tw, nwu, (42)
s nw rt " S nw rt pp c q " max ´1 ´exp ´´p c R rt ¯, 0 ¯, rt P tm, f u,
with R m " 10 4 Pa and R f " 10 Pa. The matrix is homogeneous and isotropic, i.e. K m " Λ m I, characterized by a permeability Λ m " 3¨10 ´15 m 2 , an initial porosity φ 0 m " 0.2, effective Lamé parameters λ " 833 MPa, µ " 1250 MPa, effective (drained) bulk modulus1 K dr " λ`µ " 2083 MPa, and solid grain bulk modulus K s " 11244 MPa. From these, one can infer the values of the Biot coefficient b " 1 ´Kdr Ks » 0.81, and of the Biot modulus M " Ks b´φ 0 m » 18.4 GPa. Since we consider a horizontal domain with no gravity effect, we set f " 0 in Ω and no gravity term appears in the Darcy laws as in [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF]. The domain is assumed to be clamped all over its boundary, i.e. u " 0 on p0, T q ˆBΩ; for the flows, we impose a wetting saturation s w m " 1 on the north side of the boundary p0, T q ˆpp0, Lq ˆtLuq, whereas the remaining part of the boundary is considered as impervious (q α m ¨n " 0, α P tnw, wu). The system is subject to the initial conditions p nw 0 " p w 0 " 10 5 Pa, which in turn results in an initial saturation s nw 0,rt " 0, rt P tm, f u. The final time is set to T " 1000 days " 8.64¨10 7 s. The system is excited by the following source term, representing injection of non-wetting fluid at the center of the fracture network:

h nw f pt, xq " gpxq ż Γ gpxq dσpxq V por T {5
, pt, xq P p0, T q ˆΓ, where V por " ş Ω φ 0 m pxq dx is the initial porous volume and gpxq " e ´β|px´x0q{L| 2 , x 0 " p L 2 , L 2 q, with β " 1000 and |¨| the Euclidean norm. The remaining source terms h w f and h α m , α P tw, nwu, are all set to zero. As mentioned in the introduction, the GDM framework covers many possible schemes for both the flow and mechanical components of the model. For the flow, one would typically consider finite volume methods (or mixed finite elements), such as the low-order two-and multi-point flux approximations or hybrid mimetic mixed schemes [START_REF] Droniou | The Gradient Discretisation Method[END_REF]Chapters 12,[START_REF] Bonaldi | Two-phase darcy flows in fractured and deformable porous media, convergence analysis and iterative coupling[END_REF]; even though high-order finite volume methods such as the hybrid high-order scheme [START_REF] Antonio | The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications[END_REF] or non-conforming virtual elements [START_REF] Ayuso De Dios | The nonconforming virtual element method[END_REF] also fit into the GDM [START_REF] Di Pietro | Discontinuous Skeletal Gradient Discretisation methods on polytopal meshes[END_REF], their usage in the current context seem less justified given the expected low regularity of the solution. Given that our simulations are done on triangular meshes, we opted to discretise the flow part using the cheap and robust Two-Point Flux Approximation (TPFA). For the elasticity equation in (3), standard conforming finite element methods in standard displacement formulation as well as other more advanced techniques (such as stabilised nodal strain formulation or Hu-Washizu-based formulations) are known to fit in the GDM [START_REF] Droniou | Gradient schemes for linear and non-linear elasticity equations[END_REF]. Our choice was on the second order P 2 finite element in displacement formulation in the matrix [START_REF] Daïm | A preconditioned conjugate gradient based algorithm for coupling geomechanical-reservoir simulations[END_REF][START_REF] Jeannin | Accelerating the convergence of coupled geomechanicalreservoir simulations[END_REF], adding supplementary unknowns on the fracture faces to account for the discontinuities. It provides a better accuracy than P 1 finite element especially on the normal stresses at fracture tips and intersections.

The adaptation of the TPFA discretization to the hybrid-dimensional two-phase Darcy flow model follows [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF] using mf -linear m-upwind model for matrix-fracture interactions. However, unlike [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF], we consider here a centered approximation of the mobilities, and the scheme used in the test can therefore be written as a gradient scheme (8a)-(8b).

Figure 4: Example of admissible triangular mesh with three fracture edges in bold. The dot lines joining each cell center to the center of each of its edges are assumed orthogonal to the edge. The discrete unknowns are presented for the two-phase flow and the mechanics. Note that the discontinuities of the saturations and of the displacement are captured at matrix fracture interfaces. The matrix and fracture saturations s α m , s α f at matrix-fracture interfaces are computed using a single primary unknown parametrizing the capillary pressure graphs (cf. [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF]). Note also that additional nodal unknowns are defined at intersections of at least three fractures.

The discrete unknowns for the phase pressures, the phase saturations and for the displacement field are shown in Fig- ure 4. The computational domain Ω is decomposed using admissible triangular meshes for the TPFA scheme (cf. [34, Section 3.1.2] and the example Figure 4). Let n P N ‹ denote the time step index. The time stepping is adaptive, defined as δt n`1 2 " mint δt n´1 2 , ∆t max u, where δt 1 2 " 0.025 days is the initial time step, ∆t max " 5 days is the maximal time step (except for the finest mesh for which it is set to ∆t max " 2 days), and " 1.1. At each time step, the flow unknowns are computed by a Newton-Raphson algorithm. At each Newton-Raphson iteration, the Jacobian matrix is computed analytically and the linear system is solved using a GMRes iterative solver. The time step is reduced by a factor 2 whenever the Newton-Raphson algorithm does not converge within 50 iterations, with the stopping criteria defined by the relative residual norm lower than 10 ´5 or a maximum normalized variation of the primary unknowns lower than 10 ´4. On the other hand, given the matrix and fracture equivalent pressures p E m and p E f , the displacement field u is computed using the direct solver MA48 (see [START_REF] Duff | The design of MA48: a code for the direct solution of sparse unsymmetric linear systems of equations[END_REF]). Following [START_REF] Bevillon | Stability and convergence analysis of partially coupled schemes for geomechanicalreservoir simulations[END_REF][START_REF] Daïm | A preconditioned conjugate gradient based algorithm for coupling geomechanical-reservoir simulations[END_REF][START_REF] Jeannin | Accelerating the convergence of coupled geomechanicalreservoir simulations[END_REF][START_REF] Kim | Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits[END_REF][START_REF] Mikelić | Convergence of iterative coupling for coupled flow and geomechanics[END_REF][START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF], the coupling between the two-phase Darcy flow and the mechanical deformation is solved by means of a fixed-point algorithm. This algorithm computes the matrix porosity and the fracture aperture, using discrete versions of the coupling laws (4), at each time step and fixed-point iteration. The algorithm is summarized in the following scheme, where k denotes the current fixed-point iteration and n the current time step.

Iterative coupling algorithm

At each time step n, for k " 1, . . . , until convergence, solve the following Darcy and mechanical subproblems: (i) Compute p α,n,k rt , s α,n,k rt , α P tw, nwu, rt P tm, f u, solving the Darcy flow model using d n,k´1 f in the fracture conductivity and the following porosity and fracture aperture in the accumulation term:

$ & % φ n,k m ´φn´1 m " C r,m pp E,n,k m ´pE,n,k´1 m q `b divpu n,k´1 ´un´1 q `1 M pp E,n,k m ´pE,n´1 m q, d n,k f ´dn´1 f " C r,f pp E,n,k f ´pE,n,k´1 f q ´ u n,k´1 ´un´1 .
(ii) Compute the displacement field u n,k using the equivalent pressures p E,n,k m and p E,n,k f computed at step (i).

Initialization

For given n ą 1, set $ ' ' & ' ' % p E,n,0 rt ´pE,n´1 rt δt n´1 2 " p E,n´1,0 rt ´pE,n´2 rt δt n´3 2 , rt P tm, f u, u n,0 ´un´1 δt n´1 2 " u n´1 ´un´2 δt n´3 2 ; For n " 1, set # p E,´1 rt " p E,0 rt , rt P tm, f u, u ´1 " u 0 .
Here, C r,m and C r,f are positive relaxation parameters mimicking the rock compressibility (see e.g. [START_REF] Bevillon | Stability and convergence analysis of partially coupled schemes for geomechanicalreservoir simulations[END_REF][START_REF] Daïm | A preconditioned conjugate gradient based algorithm for coupling geomechanical-reservoir simulations[END_REF][START_REF] Jeannin | Accelerating the convergence of coupled geomechanicalreservoir simulations[END_REF][START_REF] Kim | Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits[END_REF][START_REF] Mikelić | Convergence of iterative coupling for coupled flow and geomechanics[END_REF][START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF]). For our numerical simulations, we choose C r,m " 16b 2 2µ`2λ (cf. [START_REF] Mikelić | Convergence of iterative coupling for coupled flow and geomechanics[END_REF]), and C r,f " r d f C r,m with r d f " 10 ´3 m. The convergence of this fixed-point algorithm is achieved if the relative norm of the displacement field increment between two successive iterations is lower than 10 ´5.

Numerical convergence

To verify the convergence of the method, we take into account six refined admissible triangular grids with N " N 0 , 4N 0 , 16N 0 , 64N 0 , 256N 0 , 1024N 0 cells, N 0 " 224. All the numerical experiments of this subsection consider the centered approximation of the degenerate mobilities [START_REF] Hanowski | The hydromechanical equilibrium state of poroelastic media with a static fracture: A dimension-reduced model with existence results in weighted Sobolev spaces and simulations with an XFEM discretization[END_REF]. The non-degenerate regularization consisting in replacing the mobilities with µ α η α rt ps α q ` µ α p1 ` q (for rt P tm, f u and α P tw, nwu)

has also been investigated, and found to exhibit significant differences, compared to the degenerate case, mainly on the matrix saturations and only for ě 10 ´3; the differences are small for " 10 ´4 and not observable for ď 10 ´5.

Figure 7 shows the convergence of the displacement field and gas saturation profiles along the line y " 55 m, intersecting the vertical fracture, computed at the final time for the first five grids. In addition, we consider a reference solution (denoted with the subscript ref) computed on the finest (sixth) grid, made up by 1024N 0 " 229376 cells, and used to showcase the time histories of the solution as well as to compute the time histories of the relative errors for each grid. Figure 6 shows the variation with respect to the curvilinear abscissa (x or y, depending on the orientation) of the initial and final apertures for the fractures in the cross-shaped network, based on the reference solution. Note that the non-symmetry of the y plot results from the output boundary condition on the north-side. At time t " 0, the widths of both xand y-oriented fractures coincide. Figure 5 displays the final non-wetting matrix pressure and saturation computed on the fifth grid; as expected, the non-wetting fluid accumulates at the tips, flows through the fracture network and is attracted towards the upper open boundary. Figure 8 showcases the time histories of the average of some relevant physical quantities computed based on the reference solution (the average of a is denoted by a ‹ ). In particular, we notice the increase in width for the fracture network as a result of the gas injection, followed by a decrease after attaining a maximum due to an increasing gas matrix mobility in the neighborhood of the fractures. The same remark holds for the equivalent pressure p E m . The mean saturation in the matrix, as expected, grows linearly with time until the gas front reaches the upper boundary. To illustrate the spatial convergence of the scheme, Figure 7 plots on 4 meshes the cuts at y " 55 m of both components of the displacement field and of the matrix non-wetting saturation. The non-monotone profile of the saturation cut results from the fronts propagating from the different tips of the fractures. Figure 9 shows the convergence of the errors 1: Performance of the method with the centered scheme in terms of the number of mesh elements, the number of successful time steps, the total number of Newton-Raphson iterations, the total number of GMRes iterations, the total number of fixed-point iterations, and the CPU time. as a function of the mesh step for the first five meshes and a " d f , s nw m , p E m . Computations are carried out, again, using averaged quantities (a ‹ N denotes the spatial average of quantity a computed using N triangular elements). A similar convergence rate is observed for all quantities. Finally, we give an insight into the performance of our method in Table 1, where

˜şT 0 pa ‹ N ptq ´a‹ ref ptqq 2 dt ş T 0 a ‹ ref ptq 2 dt ¸1 2 , NbCells N ∆t N Newton N GMRes N FixedPoint CPU (
• NbCells is the number of cells of the mesh,

• N ∆t is the number of successful time steps,

• N Newton is the total number of Newton-Raphson iterations,

• N GMRes is the total number of GMRes iterations,

• N FixedPoint is the total number of fixed point iterations,

• CPU (s) is the CPU time of the simulation in seconds.

The iterative coupling algorithm exhibits good robustness with respect to the mesh size, and has a linear convergence behaviour as proved in [START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF] in the linear case. As expected for incompressible fluids in fractured porous media [START_REF] Girault | Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium[END_REF], the convergence rate is however very sensitive to small initial time steps. This issue is shown in [START_REF] Bonaldi | Two-phase darcy flows in fractured and deformable porous media, convergence analysis and iterative coupling[END_REF] to be efficiently solved by using a Newton Krylov acceleration of the fixed-point algorithm.

x → d f,ref (0, (x, L 2 )) x → d f,ref (T, (x, L 2 )) y → d f,ref (T, ( L 2 , y))
Fracture width vs. x or y at times t = 0 and t = T Table 2: Performance of the method with the upwind scheme in terms of the number of mesh elements, the number of successful time steps, the total number of Newton-Raphson iterations, the total number of GMRes iterations, the total number of fixed-point iterations, and the CPU time.

Comparison between the upwind and centered schemes

In this subsection, the solutions obtained by the upwind approximation of the mobilities [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF] are compared to the solutions obtained in the previous subsection based on the centered approximation. Small differences can be noticed on the upper part of the matrix saturation at final time in Figures 10 and 11 due to the larger numerical diffusion of the upwind scheme. They clearly reduce with the mesh refinement as can be observed from the comparison between the line cut plots in these figures. The cuts at y " 80 m in the right Figure 12 for the fifth and sixth meshes confirm that, compared to the upwind scheme, the centered scheme converges more quickly with the mesh refinement. This is also checked on the left Figure 12 which exhibits the convergence of the error for the average quantities d ‹ f , s nw,‹ m and p E,‹ m with respect to the reference mesh solutions. It shows that the upwind scheme is slightly more accurate on the two coarsest meshes but that the centered scheme has a better convergence rate. It has been checked that the displacement field and fracture aperture for both schemes exhibit very little differences on the finest meshes. Table 2 shows the numerical performances of the upwind scheme. Comparing with the results obtained for the centered scheme (Table 1), the upwind scheme offers a slighly more efficient nonlinear convergence. It was also checked that the upwind scheme can accommodate larger time steps with successful nonlinear convergence than the centered scheme.

We note that a GDM discretisation of ( 3)-( 4) with upwind approximation of the mobilities as in [START_REF] Aghili | Two-phase discrete fracture matrix models with linear and nonlinear transmission conditions[END_REF] for the TPFA scheme requires a scheme-dependent definition of the Darcy fluxes, and hence loses the generality of the GDM. Once such definition is provided, the convergence analysis would entail a specific treatment which needs to be investigated due to the nonlinearity introduced by the upwinding; however, this analysis could benefit from the tools developed here, such as the relative compactness results.

Conclusions

We developed, in the framework of the gradient discretization method, the numerical analysis of a two-phase flow model in deformable and fractured porous media. The model considers a linear elastic mechanical model with open fractures coupled with an hybrid-dimensional two-phase Darcy flow assuming continuity of each phase pressure across the fractures. The model accounts for a general network of planar fractures including immersed, non-immersed fractures and fracture intersections, and considers different rock types in the matrix and fracture network domains.

It is assumed, for the convergence analysis, that the porosity remains bounded below by a strictly positive constant, and that the fracture aperture remains larger than a fixed non-negative continuous function vanishing only at the tips and p E,‹ m with respect to the corresponding reference time histories (cf. Figure 8). of the fracture network. These assumptions stem from the limitations of the continuous model itself. In addition, the mobility functions are assumed to be bounded below by strictly positive constants. However, unlike previous works, the fracture conductivity d 3 f {12 was not frozen and the complete non-linear coupling between the flow and mechanics equations was considered.

Assuming that the gradient discretization meet generic coercivity, consistency, limit-conformity and compactness properties, we proved the weak convergence of the phase pressures and displacement field to a weak continuous solution, as well as the strong convergence of the fracture aperture and of the matrix and fracture saturations. Numerical experiments carried out for a cross-shaped fracture network immersed in a two-dimensional porous medium and using a TPFA finite volume scheme for the flow combined with a P 2 finite element method for the mechanics, confirmed the numerical convergence of the scheme.

A Appendix

A.1 Appendix 1

Proposition A.1. Let X Ă R d be bounded, δ ą 0 and let `Aδ m ˘mPM δ be a covering of X in disjoint cubes of length δ. Let R δ : L 2 pR d q Ñ L 2 pXq be such that, for any v P L 2 pR d q, `Rδ v ˘|A δ m XX "

1 δ d ż A δ m vpxq dx @m P M δ ,
Then, we have }R δ v ´v} L 2 pXq ď 2 d{2 sup |z|ďδ }vp¨`zq ´v} L 2 pXq .

Proof. The proof can be found in [29, p. 756]. Note that the assumption, in this reference, that v is zero outside X is actually not useful.

Lemma A.2. Let X Ă R d be bounded, and U be an open subset of R d such that tx P R d : distpx, Xq ă δ 0 u Ă U for a given δ 0 ą 0, where the distance is considered for the supremum norm in R d . Let pw k q kPN be a bounded sequence in L 8 p0, T ; L 2 pU qq that converges uniformly in time and weakly in L 2 pU q to w P L 8 p0, T ; L 2 pU qq. Let p P r1, `8s and let us define T pδq " sup .

If lim δÑ0 T pδq " 0, then the sequence pw k q kPN converges to w in L p p0, T ; L 2 pXqq.

Proof. For 0 ă δ ă δ 0 , let `Aδ m ˘mPM δ be a covering of X in disjoint cubes of length δ and let R δ be the corresponding L 2 projection operator as defined in Proposition A.1. We write w k ´w " pw k ´Rδ w k q `pR δ w k ´Rδ wq `pR δ w ´wq implying that }w k ´Rδ w k } L p p0,T ;L 2 pXqq À T pδq.

Setting v k " w k ´Rδ w k , k P N, we have, if p " 8, }v k pt, ¨q} L 2 pXq À T pδq for a.e. t P p0, T q. Since Id ´Rδ : L 2 pXq Ñ L 2 pXq is linear, the weak convergence of w k pt, ¨q implies that v k pt, ¨q á vpt, ¨q -wpt, ¨q ´Rδ wpt, ¨q weakly in L 2 pXq, and thus that }wpt, ¨q ´Rδ wpt, ¨q} L 2 pXq ď lim inf kÑ`8 }v k pt, ¨q} L 2 pXq À T pδq.

For p ă 8, we have, using the above weak convergence of pv k pt, ¨qq kPN and Fatou's lemma, Hence, for any p, }w ´Rδ w} L p p0,T ;L 2 pXqq À T pδq.

Proposition A.4. Let ū P U 0 be the solution of (43) and u P X 0 Du the solution of the gradient scheme [START_REF] Jaffré | A discrete fracture model for two-phase flow with matrix-fracture interaction[END_REF]. Then, there exists a hidden constant depending only on the coercivity constant C Du and on the physical data such that the following error estimate holds

As a consequence, if pD l u q lPN is a sequence of coercive, consistent and limit-conforming GDs, if u l is the solution of (44) for D u " D l u , if pD l p q lPN is a sequence of GDs and p α,l P X 0

Proof. We note that even though W Du was considered, in the definition of limit-conformity of a sequence of GDs, only on C 8 Γ pΩzΓ, S d pRqq, it can be defined on H div,Γ pΩzΓ; S d pRqq -σ P L 2 pΩ; S d pRqq : divpσq| Ω β P L 2 pΩ β q d , β P Ξ, σ `n``σ´n´" 0 on Γ, pσ `n`q ˆn`" 0 on Γ ( , where pΩ β q βPΞ are the connected components of ΩzΓ. Setting σ " σpūq ´b pE m I P H div,Γ pΩzΓ; S d pRqq as an argument of W Du and using divσ " ´f , we obtain that for all v P X 0 Combined with the definition of C Du , the estimates above establish [START_REF] Jeannin | Accelerating the convergence of coupled geomechanicalreservoir simulations[END_REF].

Under the assumptions in the second part of the proposition, the hidden constant in (45) is independent of l, the last two terms in the left-hand side of this estimate converge to 0 as l Ñ `8, as well as S D l u puq by definition of the consistency of the sequence of GDs. When its argument σ is in the vector space C 8 Γ pΩzΓ, S d pRqq, W D l u pσq also converges to 0 by definition of limit-conformity; since this space is dense in H div,Γ pΩzΓ; S d pRqq, the arguments in [START_REF] Droniou | The Gradient Discretisation Method[END_REF]Lemma 2.17] show that this convergence also holds for the argument σ " σpūq ´b pE m I. Estimate (45) therefore yields the convergences [START_REF] Jha | Coupled Modeling of Multiphase Flow and Fault Poromechanics during geologic CO 2 storage[END_REF].

Figure 1 :

 1 Figure 1: Example of a 2D domain Ω with three intersecting fractures Γ i , i P t1, 2, 3u.

Lemma 4 . 4 (

 44 A priori estimates). Let p α , u be a solution to problem (8) such that (i) d f,Du pt, xq ě d 0 pxq for a.e. pt, xq P p0, T q ˆΓ, (ii) φ D pt, xq ě φ m,min for a.e. pt, xq P p0, T q ˆΩ, where φ m,min ą 0 is a constant.

4. 2 . 3

 23 Uniform-in-time L 2 -weak convergence of φ D Π m Dp s α m and φ D Proposition 4.8.

  m converges uniformly in time, weakly in L 2 pΩq. We deduce that the sequence φ D l " ř αPtnw,wu φ D l Π m D l p s α,l m also converges uniformly in time, weakly in L 2 pΩq. 4.2.4 Uniform-in-time L 2 -weak convergence of d f,Du Π f Dp s α f and d f,Du

1 2 `

 2 Dp , with lim hÑ0 ωphq " 0 and Dp " }ϕ´Π m Dp ϕ} L 2 pΩq `} φ´Π f Dp ϕ} L 2 pΓq a consistency error term such that lim lÑ`8 D l p " 0. It follows from the discontinuous Ascoli-Arzelà theorem [30, Theorem C.11] that (up to a subsequence) the sequence d f,Du Π f Dp s α f converges uniformly in time weakly in L 2 pΓq. Summing over α P tnw, wu, we also deduce the uniform-intime L 2 pΓq-weak convergence of d f,Du . 4.2.5 Strong convergence of d f,Du , d f,Du Π f Dp s α f , and Π f Dp s α f Proposition 4.10.

Figure 3 :

 3 Figure 3: Computational domain.

  (a) p nw m pT ; x, yq (b) s nw m pT ; x, yq (c) u 1 pT ; x, yq (d) u 2 pT ; x, yq

Figure 5 :

 5 Figure 5: Final non-wetting matrix pressure and saturation (a)-(b), and final (c)-(d) displacement field on the fifth mesh of size 57344 cells.

  x or y (m)

Figure 6 :

 6 Figure6: Initial and final widths the xand y-oriented fractures vs. corresponding curvilinear abscissae, computed using the finest grid (reference solution). The initial width for both the xand y-oriented fractures is the same.NbCells N ∆t N Newton N GMRes N FixedPoint CPU (s) 256N 0 246 4809 82085 4054 2250 1024N 0 537 5486 114763 4136 13600

Figure 7 :Figure 8 :

 78 Figure 7: Convergence of the profiles of the displacement field components (m) u 1 (a), u 2 (b), and gas saturation s nw m (c) at the final time, along the line y " 55 m intersecting the vertical fracture, for four grids, with N triangular elements, and N 0 " 224.

Figure 9 :

 9 Figure 9: Relative L 2 norm of the error as a function of the mesh step computed on the first five meshes for the time histories of the mean quantities d ‹ f , s nw,‹ m

Figure 10 :

 10 Figure 10: On the fifth mesh with 256N 0 cells, final non-wetting phase matrix saturations for the centered scheme (top left) and upwind scheme (top right), and line cuts of both solutions at y " 55 m and y " 80 m (bottom).

Figure 11 :

 11 Figure 11: On the finest mesh with 1024N 0 cells, final non-wetting phase matrix saturations for the centered scheme (top left) and upwind scheme (top right), and line cuts of both solutions at y " 55 m and y " 80 m (bottom).

Figure 12 :

 12 Figure12:(Top): on the fifth (256N 0 cells) and sixth (1024 N 0 cells) meshes, cuts at y " 80 m of the non-wetting phase matrix saturations for both the centered and upwind schemes. (Bottom): for both the centered and upwind schemes, relative L 2 norm of the error as a function of the mesh step computed on the first five meshes for the time histories of the mean quantities d ‹ f , s nw,‹ m and p E,‹ m with respect to the corresponding reference time histories.

  }w k p¨, ¨`zq ´wk p¨, ¨q} L 2 pXq › › › L p p0,T q

and we establish the convergence to 0

 0 of each bracketed term in the right-hand side. First, in view of Proposition A.1 }w k pt, ¨q ´Rδ w k pt, ¨q} L 2 pXq À sup |z|ďδ }w k pt, ¨`zq ´wk pt, ¨q} L 2 pXq

ż T 0 }vpt, ¨q} p L 2

 02 pXq dt ď ż T 0 lim inf kÑ`8 }v k pt, ¨q} p L 2 pXq dt ď lim inf kÑ`8 ż T 0 }v k pt, ¨q} p L 2 pXq dt À T p pδq.

}ε

  Du puq ´εpūq} L 2 pΩ,S d pRqq `}Π Du u ´ū} L 2 pΩq `} u Du ´ ū } L 2 pΓq À S Du pūq `WDu pσpūq ´b pE m Iq `}p E m ´Πm Dp p E m } L 2 pΩq `}p E f ´Πf Dp p E f } L 2 pΓq .

  f in L 2 pΓq, then, as l Ñ `8,ε D l u pu l q Ñ εpūq in L 2 pΩ, S d pRqq, Π D l u u l Ñ ū in L 2 pΩq d , u l D l u Ñ ū in L 2 pΓq.

Ω

  ´pσpūq ´σDu puqq : ε Du pvq ´bpp E m ´Πm Dp p E m qdiv Du pvqq ¯dx `żΓ pp E f ´Πf Dp p E f q v Du dσpxq ˇˇď }v} Du W Du pσpūq ´b pE m Iq. Setting v " P Du ū ´u, where P Du ū realizes the minimum in S Du pūq, we infer }P Du ū ´u} Du À S Du pūq `WDu pσpūq ´b pE m Iq `}p E m ´Πm Dp p E m } L 2 pΩq `}p E f ´Πf Dp p E f } L 2 pΓq .

  Du : X 0 Du Ñ L 2 pΩ, S d pRqq, • a displacement function reconstruction linear operator Π Du : X 0 Du Ñ L 2 pΩq d , • a normal jump function reconstruction linear operator ¨ Du : X 0 Du Ñ L 4 pΓq, where S d pRq is the vector space of real symmetric matrices of size d. Let us define the divergence and stress tensor operators by div Du pvq " Tracepε Du pvqq and σ Du pvq " 2µε Du pvq `λ div Du pvqI, and the fracture width d f,Du " ´ u Du . It is assumed that the following quantity defines a norm on X 0 Du : }v} Du -}ε Du pvq} L 2 pΩ,S d pRqq .

	0 Du
	and
	• a discrete symmetric gradient linear operator ε

Table

						s)
	N 0	246	11902	81037	11163	23
	4N 0	246	4685	47352	4234	35
	16N 0	246	4626	53870	4138	130
	64N 0	246	4713	65293	4063	500
	256N 0	246	4951	88106	4062	2600
	1024N 0	537	5788	125495	4147	14500

In general, K dr " λ `2µ{d, where d P t2, 3u is the space dimension.

Acknowledgements We are grateful to Andra and to the Australian Research Council's Discovery Projects (project DP170100605) funding scheme for partially supporting this work.

Finally, R δ w k ´Rδ w " ÿ mPM δ a δ km ptq1 A δ m XX , with a δ km ptq "

Since the covering `Aδ m ˘mPM δ is finite and since, for all m P M δ , the term a δ km ptq converges uniformly in time to zero, it results that pR δ w k ´Rδ wq converges as k Ñ `8 to zero in L 8 pp0, T q ˆXq.

Gathering the estimates, we have that }w k ´w} L p p0,T ;L 2 pXqq À 2T pδq `}R δ w k ´Rδ w} L p p0,T ;L 2 pXqq .

Passing to the superior limit as k Ñ `8, we deduce that lim sup kÑ`8 }w k ´w} L p p0,T ;L 2 pXqq À 2T pδq which yields, letting δ Ñ 0, lim sup kÑ`8 }w k ´w} L p p0,T ;L 2 pXqq " 0.

A.2 Appendix 2

Lemma A.3. Let pD l u q lPN be a sequence of GDs assumed to satisfy the coercivity and limit-conformity properties. Let pu l q lPN be a sequence of vectors with u l P X 0 D l u such that there exist C independent of l P N with }u l } Du ď C. Then, there exists ū P U 0 such that, up to a subsequence, the following weak limits hold:

Proof. By assumption the sequence p}ε D l u } L 2 pΩ,S d pRqq q lPN is bounded which implies, from the coercivity property, that the sequences p}Π D l u u l } L 2 pΩq q lPN and p} u l D l u } L 2 pΓq q lPN are also bounded. Hence there exist ū P L 2 pΩq d , ε P L 2 pΩ, S d pRqq and ḡ P L 2 pΓq such that, up to a subsequence, one has

Passing to the limit in the definition of the limit-conformity yields, for any σ P C 8 Γ pΩzΓ, S d pRqq, ż Ω ´σ : ε `ū ¨divpσq ¯dx ´żΓ pσn `q ¨n`ḡ dσpxq " 0.

Selecting first σ with a compact support in ΩzΓ, and then a generic σ, it results that ū P U 0 with ε " εpūq and ḡ " ū . Since div D l u pu l q " Tracepε D l u pu l qq, it also holds that div D l u pu l q á divpūq in L 2 pΩq.

Let us fix pα P V 0 , α P tnw, wu, f P L 2 pΩq d , and define

pα S α m pp c q ´Um pp c q and pE f " ÿ αPtnw,wu γ pα S α f pγ pc q ´Uf pγ pc q.

with pc " pnw ´p w . We consider the solution ū P U 0 of the following variational formulation

Let us take p α P D p , α P tnw, wu, p c " p nw ´pw and p E m " ÿ αPtnw,wu p α S α m pp c q ´Um pp c q and p E f " ÿ αPtnw,wu p α S α f pp c q ´Uf pp c q.

We consider the following gradient scheme for (43): Find u P X 0 Du such that, for all v P X 0 Du , ż

The Lax-Milgram theorem ensures that the exact solution ū and approximate solution u exist and are unique. The following proposition provides an error estimate.