Gradient discretization of two-phase poro-mechanical models with discontinuous pressures at matrix fracture interfaces

We consider a two-phase Darcy flow in a fractured and deformable porous medium for which the fractures are described as a network of planar surfaces leading to so-called hybrid-dimensional models. The fractures are assumed open and filled by the fluids and small deformations with a linear elastic constitutive law are considered in the matrix. As opposed to [10], the phase pressures are not assumed continuous at matrix fracture interfaces, which raises new challenges in the convergence analysis related to the additional interfacial equations and unknowns for the flow. As shown in [16, 2], unlike single phase flow, discontinuous pressure models for two-phase flows provide a better accuracy than continuous pressure models even for highly permeable fractures. This is due to the fact that fractures fully filled by one phase can act as barriers for the other phase, resulting in a pressure discontinuity at the matrix fracture interface. The model is discretized using the gradient discretization method [22], which covers a large class of conforming and non conforming schemes. This framework allows for a generic convergence analysis of the coupled model using a combination of discrete functional tools. In this work, the gradient discretization of [10] is extended to the discontinuous pressure model and the convergence to a weak solution is proved. Numerical solutions provided by the continuous and discontinuous pressure models are compared on gas injection and suction test cases using a Two-Point Flux Approximation (TPFA) finite volume scheme for the flows and $P_2$ finite elements for the mechanics.


Introduction
Coupled flow and geomechanics play an important role in many subsurface processes such as water management, geothermal energy, CO 2 sequestration, oil and gas production or nuclear waste storage. This is particularly the case in the presence of fractures which have a strong impact both on the flow and on the rock mechanical behavior. This work considers the so called hybrid-dimensional or Discrete Fracture Matrix (DFM) models representing the fractures as a network of co-dimension one surfaces coupled with the surrounding matrix domain. The reduced flow model is then obtained by averaging both the unknowns and the equations in the fracture width and by imposing appropriate transmission conditions at both sides of the matrix fracture interfaces. The mechanical model is set on the matrix domain with appropriate boundary conditions on both sides of the fracture interfaces. 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 considerable impact on the flow and transport of mass and energy in porous media, and on the mechanical behavior of rocks. For the derivation and analysis of such models, let us refer to [4,27,40,44,6,13,15,46] for singlephase Darcy flows, [9,49,45,36,14,23,16,2] for two-phase Darcy flows, and [41,42,38,33,34,29,39,30,52,32] for poroelastic models.
As in [10], this work focuses on a hybrid-dimensional two-phase Darcy flow model coupled with a linear poroelastic deformation of the matrix. The fractures are assumed to remain open and fully filled by the fluids, and their propagation over time is neglected. The Poiseuille law is used for the tangential velocity along the fracture network and extended to two-phase flow based on generalized Darcy laws. As in [17], the concept of equivalent pressure, used to extend the poro-mechanical coupling to two-phase flow, is based on the capillary energy density. This is a crucial choice to obtain the stability of the coupled model.
In [10], the continuity of both phase pressures is used as a transmission condition at matrix fracture interfaces. This is a classical assumption in the case of highly permeable fractures such as open fractures. As shown e.g. in [28], this choice is fully justified for the case of single-phase flows. On the other hand, in the case of two-phase flow, this assumption can lead to inaccurate solutions at the matrix fracture interfaces [16,2]. This is in particular the case when the fractures are fully filled by one phase and act as barriers for the other phase due to its very low relative permeability within the fractures, hence leading to a pressure discontinuity. Let us refer to [2] for striking examples including the desaturation by suction at the interface between the atmosphere and a low permeable and fractured storage rock.
This potential inaccuracy of continuous pressure models motivates us to consider the extension of the analysis carried out in [10] to hybrid-dimensional discontinuous pressure flow models [36,23,16,2]. For such flow models, the Darcy fluxes between the matrix fracture interface and the fracture are modelled using a two-point flux approximation combined with an upwind approximation of the mobilities [23,16,2]. Following [23], the model also includes a layer of damaged rock at matrix fracture interfaces. This additional accumulation term plays a major role in the numerical analysis of the model and also improves the nonlinear convergence at each time step of the simulation [23,12]. It must be kept sufficiently small to maintain the accuracy of the solution (see [23]). Following [10] and [23], this new hybriddimensional poro-mechanical model is discretized using the gradient discretization method [22]. 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 allow a generic convergence analysis based on general properties of the reconstruction operators that hold for a large class of conforming and non conforming discretizations. Let us point out that, with respect to [10], additional trace and jump operators need to be defined in this framework, along with new definitions of coercivity, consistency, limit-conformity, and compactness. The two main ingredients to discretize the coupled model are the discretizations of the hybrid-dimensional discontinuous pressure 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 covers the case of cell-centered finite volume schemes with Two-Point Flux Approximation on strongly admissible meshes [40,6,2], or some symmetric Multi-Point Flux Approximations [51,50,3] 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 [44,15,7,32], and for vertex-based discretizations such as the Vertex Approximate Gradient scheme [15,23,16]. For the discretization of the elastic mechanical model, the gradient discretization framework covers conforming finite element methods such as in [33], the Crouzeix-Raviart discretization [35,20], the Hybrid High Order discretization [19], and the Virtual Element Method [8].
The main objective of this work is to introduce the gradient discretization of the hybrid-dimensional poro-mechanical model with discontinuous pressure at matrix fracture interfaces. Then, we prove the convergence of the discrete solution to a weak solution of the model. Compared with [10], new difficulties arise from the interfacial additional nonlinear flux and accumulation terms including the damaged rock type. Assuming that the fracture normal transmissivity is fixed in the interfacial two-point fluxes, i.e. that its fracture aperture dependence is frozen, we are able to prove the convergence of the gradient scheme solution to a weak solution of the model. This assumption is rather mild since, in practice, the solution depends only weakly on this fracture normal transmissivity as long as it remains much larger than the matrix transmissivity. Concerning compactness estimates, the same techniques as in [10] are used: time translates, uniform-in-time L 2 -weak estimates, and a discrete version of the Ascoli-Arzelà theorem. In [10], where fields defined in matrix and fracture are related, matrix and fracture contributions have to be separated by a cut-off argument (since the fracture width vanishes at tips). On the other hand, in this work, such a separation stems from the model itself, but the damaged rock layer has to be embedded in the time translates of the saturations, by using ad-hoc test functions combining the matrix and damaged layer rock types.
As in [10], the proof additionally assumes that the matrix porosity remains bounded from below by a strictly positive constant, that the fracture aperture remains larger than a fixed aperture vanishing only at the tips, and that the mobility functions are bounded from below by strictly positive constants. The assumptions on the porosity and fracture aperture cannot be avoided since the continuous model does not ensure these properties, which are needed to ensure its well-posedness. The assumption on the mobilities are classical to carry out the stability and convergence analysis of two-phase Darcy flows with heterogeneous rock types (see [26,14,23]).
The second objective of this work is to compare the discontinuous pressure poro-mechanical model investigated in this work to the continuous pressure poro-mechanical model presented in [10]. Two test cases are considered. As in [10], the first test case simulates the gas injection in a cross-shaped fracture network immersed in an initially water saturated porous medium. The second test case models the desaturation of a low permeable medium by suction at the interface with a ventilation tunnel. The data set of this second test case is based on the Callovo-Oxfordian argilite rock properties of the nuclear waste storage prototype facility of Andra. The geometry uses an axisymmetric DFM model based on a simplified version of the fracture network at the interface between the storage rock and the ventilation tunnel. In both cases the discretization is based on the Two-Point Flux Approximation finite volume scheme for the flows and second-order finite elements for the mechanical deformation.
The rest of the article is organized as follows. Section 2 introduces the continuous hybrid-dimensional coupled model with discontinuous pressures at matrix fracture interfaces. 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. This convergence falls short, in general, from identifying the limit matrix-fracture nonlinear fluxes; this issue is discussed in Subsection 4.4, in which an assumption is given on the limit fracture width under which the fluxes can be fully identified. In Section 5, devoted to numerical experiments, the discontinuous pressure model is compared to the continuous pressure model presented in [10].

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 where each fracture Γ i Ă Ω, i P I is a planar polygonal simply connected open domain. Without restriction of generality, we will assume that the fractures may intersect exclusively at their boundaries (see Figure 1), that is, for any i, j P I, i ‰ j one has Γ i X Γ j " H, but not necessarily Γ i X Γ j " H. The two sides of a given fracture of Γ are denoted by˘in the matrix domain, with unit normal vectors n˘oriented outward of the sides˘. We denote by γ a the trace operators on the sides a "˘of Γ for functions in H 1 pΩzΓ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´for 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 fieldv P H 1 pΩzΓq d . Let us fix a continuous function d 0 : Γ Ñ p0,`8q with zero limits at BΓzpBΓ X BΩq (i.e. the tips of Γ) and strictly positive limits at BΓ X BΩ. The fracture aperture, denoted byd f and such thatd f "´ ū for a displacement field u P H 1 pΩzΓq d , will be assumed to satisfy the following open fracture condition d f pxq ě d 0 pxq for a.e. x P Γ.
Let us introduce some relevant function spaces. First, we denote by H 1 d0 pΓq the space made of functions v Γ in L 2 pΓq, such that d 0 ∇ τ v Γ belongs to L 2 pΓq d´1 , and whose traces are continuous at fracture intersections BΓ i XBΓ j , pi, jq P IˆI (i ‰ j) and vanish on the boundary BΓ X BΩ. We then introduce the space (1) Figure 2. Example of a 2D domain Ω with its fracture network Γ, the unit normal vectors n˘to Γ, the phase pressuresp α m in the matrix andp α f 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. (inward with respect to the damaged layer) andQ α f,a (outward with respect to the damaged layer).
for the displacement vector, and where V 0 m " tv P H 1 pΩzΓq | γ BΩv " 0u, for each matrix phase pressure, and V 0 f " H 1 d0 pΓq, for each fracture phase pressure. Forv " pv m ,v f q P V 0 , let us denote by v a " γ avm´vf , the jump operator on the side a "˘of the fractures.
The matrix, fracture and damaged rock types are denoted by the indices rt " m, rt " f , and rt "˘, respectively, and the non-wetting and wetting phases by the superscripts α " nw and α " w, respectively. Finally, for any x P R, we set x`" maxt0, xu and x´"´p´xq`.
The PDEs model reads: find the phase pressuresp α ν , ν P tm, f u, α P tnw, wu, and the displacement vector fieldū, such thatp c,ν "p nw ν´p w ν , and for α P tnw, wu, with the coupling conditions is the capillary energy density function for each rock type rt P tm, f,˘u. As already noticed in [42,38,10], 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´S nw rt . (H3) For a "˘, the widthd a and porosityφ a of the damaged rock are strictly positive constants.
(H4) 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.
(H9) The matrix permeability tensor K m is symmetric and uniformly elliptic on Ω.
Let us denote by C 8 c pr0, T qˆΩzΓq the space of smooth functionsv : r0, T sˆpΩzΓq Ñ R vanishing on BΩ and at t " T , and whose derivatives of any order admit finite limits on each side of Γ. We will also use the boldface notation Definition 2.1 (Weak solution of the model). A weak solution of the model is given byp α " pp α m ,p α f q 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

The gradient discretization method
The gradient discretization (GD) for the Darcy discontinuous pressure model, introduced in [23], is defined by a finite-dimensional vector space of discrete unknowns • two discrete gradient linear operators on the matrix and fracture domains • two function reconstruction linear operators on the matrix and fracture domains • for a "˘, jump reconstruction linear operators ¨ a Dp : X 0 Dp Ñ L 8 pΓq, and trace reconstruction linear operators T a Dp : X 0

The operators Π m
Dp , Π f Dp , T a Dp are assumed piecewise constant [22,Definition 2.12]. A consequence of the piecewiseconstant property is the following: there is a basis pe i q iPI of X 0 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 D m p by applying g component-wise, then Π m Dp gpvq " gpΠ m Dp vq and T a Dp gpvq " gpT a Dp vq. Note that the basis pe i q iPI is usually canonical and chosen in the design of X 0 D m p . The same property holds for X 0 D f p and Π f Dp . The vector space X 0 Dp is endowed with the following quantity, assumed to define a norm: }v} Dp :" }∇ m Dp v} L 2 pΩq d`}d The gradient discretization for the mechanics is defined by a finite-dimensional vector space of discrete unknowns X 0 Du and • a discrete symmetric gradient linear operator ε Du : X 0 Du Ñ L 2 pΩ, S d pRqq where S d pRq is the vector space of real symmetric matrices of size d, • 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. Let us define the divergence operator div Du p¨q " Tracepε Du p¨qq, the stress tensor operator σ 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 : 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, p 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.
Spatial operators are extended into space-time operators as follows. Let Ψ D be a spatial GDM operator defined in X 0 D with D " D u , D m p or D f p , and let w " pw n q N n"0 P pX 0 D q N`1 . Then, 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 is a space-time GDM operator, by linearity the following commutativity property holds: Ψ D δ t wpt,¨q " δ t pΨ D wpt,¨qq.
The gradient scheme for (3) consists in writing the weak formulation (6a)-(6b) with continuous spaces and operators substituted by their discrete counterparts, after a formal integration by part: find p α " pp α m , p α f q P pX 0 Dp q N`1 , α P tnw, wu, and u P pX 0 Du q N`1 , such that for all ϕ α " pϕ α m , ϕ α f q P pX 0 Dp q N`1 , v P pX 0 Du q N`1 and α P tnw, wu, with the closure equations, for ν P tm, f u and a "˘, The initial conditions are given by Dpφ 0 , and the initial displacement u 0 is the solution in X 0 Du of (8b) without the time variable and with the equivalent pressures obtained from the initial pressures pp α 0 q αPtnw,wu .

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 GDM assumptions [22], except for the chain-rule and product rule 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 [23], the spatial GD of the Darcy flow is assumed to satisfy the following coercivity, consistency, limit-conformity and compactness properties.
Coercivity of D p . Let C Dp ą 0 be defined by 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 pw m , v m q P C 8 and 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 pw ν q " 0, ν P tm, f u. 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ϕ m , ϕ f q P V 0 and ψ P L 2 pΩq, as l Ñ`8, Remark 3.1 (Consistency). In [23], the consistency is only considered for r " 2. We have here to adopt a slightly stronger assumption to deal with the coupling and non-linearity involving the fracture aperture d f . Note that, under standard mesh regularity assumptions, this stronger consistency property is still satisfied for all classical GDs.
Limit-conformity of D p . For all q " pq m , q f q P C 8 pΩzΓq dˆC 8 pΓq d´1 , ϕ a P C 8 pΓq, and v " pv m , v f q P X 0 Dp , let us define and W Dp pq, ϕ a q " max 0‰vPX 0 Dp 1 }v} Dp |W Dp pq, ϕ a , vq|. Then, a sequence of spatial GDs pD l p q lPN is said to be limitconforming if for all q P C 8 pΩzΓq dˆC 8 c pΓq d´1 and ϕ a P C 8 pΓq one has lim lÑ`8 W D l p pq, ϕ a 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Ω.
(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 Remark 3.2 (Local compactness through estimates of space translates). For K m , K f as above, set η " 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 [22,Lemma 2.21], An equivalent formulation of the local compactness property is: for all K m , K f as above, Remark 3.3 (Usual compactness property for GDs). The standard compactness property for GD is not local but global, that is, on the entire domain not any of its compact subsets (see, e.g., [22,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 use 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.
Bounds on reconstruction operators of pD l p q lPN . • Chain rule estimate. For any Lipschitz-continuous function F : R Ñ R, there is C F ě 0 such that, for all l P N, and any v m P X 0 • Product rule estimate. There exists C P ě 0 such that, for any l P N and any u l m , v l m P X 0 • Bound on the jump operator. For any l P N and any v l " pv l m , v l f q P X 0 Coercivity of pD l u q lPN . Let C Du ą 0 be defined by 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.
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.
Compactness of pD l u q lPN . For any sequence pv l q lPN P pX 0  where 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 theoretical result of this work is the following convergence theorem.
Theorem 4.1. 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 " pp α m,l , p α f,l q P pX 0 Then, there existp α " pp α m ,p α f q P L 2 p0, T ; V 0 q, α P tnw, wu,ū P L 8 p0, T ; U 0 q andQ α f,a P L 2 p0, T ; L 2 pΓqq satisfying the weak formulations (6a)-(6b) such that for α P tnw, wu and up to a subsequence 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.2 (Limit interface fluxes). The theorem states that the limit functions satisfy all but the first closure equations in (6c). It does not, however, identify the limit interface fluxesQ α f,a , α P tnw, wu. This identification requires the limit functions to satisfy an energy equality, which is known under some assumption on the limit fracture widthd f . See the discussion in Subsection 4.4 for more details.

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)-(H9), 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: Remark 4.4 (Existence of the discrete solution). Since these a priori estimates are obtained assuming lower bounds on the fracture aperture and porosity, the existence of the discrete solution cannot be deduced from these estimates and will be assumed in the following convergence analysis. As noticed in the introduction, these assumptions on lower bounds of the fracture aperture and porosity are mandatory since the model itself does not account for possible contact of fracture walls nor a nonlinear behavior of pore volume contraction, and thus cannot yield such lower bounds. The analysis of a model with contact is a topic for future work.
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 functionv such thatvptq " 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: In (8a), upon choosing ϕ α " p α we obtain T 1`T2`T3`T4`ř a"˘p T a 1`T a 2 q " T 5`T6 , with where η α a,f " η α a pT a Dp s α a q if p α a Dp ě 0 and η α a,f " η α f pΠ f Dp s α f q otherwise. 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 Summing on α P tw, nwu, we obtain Indeed, for n P t0, . . . , N´1u, by the definition (5) of the capillary energy U rt and letting π n c,rt " Π rt Dp p n c,rt , we have π n`1 c,rt pS nw rt pπ n`1 c,rt q´S nw rt pπ n c,rt qq " U rt pπ n`1 c,rt q´U rt pπ n c,rt q`ż π n`1 c,rt π n c,rt pS nw rt pqq´S nw rt pπ n c,rt qqdq ě U rt pπ n`1 c,rt q´U rt pπ n c,rt q, where the last inequality holds since S nw rt is a non-decreasing function. Thus, we obtain Applying again (16), we havê In the light of the closure equations (8c), this allows us to infer that where we have used the fact that, for v piecewise constant on r0, T s, Using, as in (18), the relation ÿ α T a Dp p α m δ t S α a pT a Dp p c,m q " T a Dp p c,m δ t S nw a pT a Dp p c,m q ě δ t U a pT a Dp p c,m q, Then, taking into account assumptions (H1)-(H9) and (i) in the lemma, there exists a real number C ą 0 depending only on the data such that On the other hand, upon choosing v " δ t u in (8b), we get T 7`T8`T9 " T 10 , with Using (20) and developing the definition of σ Du , we see that so that, all in all, taking into account  (19), (22) and (24), we obtain the following estimate for the solutions of (8): there is a real number C ą 0 depending on the data such that Now, we have 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)-(H9) and (ii) in the lemma, and using telescopic sums on the terms involving δ t , it is then possible to infer from (25) the existence of a real number C ą 0 depending on the data and on φ m,min such that he above inequality,along with the fact that T can be replaced by any t P p0, T s in the left-hand side, and in view of (11)-(36)-(37), 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 α " pp α m , p α f q 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ϕ m , ϕ f q P X 0 Dp , it holdšˇˇx 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 hašˇˇx From the gradient scheme discrete variational equation (8a), we deduce thaťˇˇx 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 . The result follows from (27), (28), the a priori estimates of Lemma 4.3, 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 (26) on α P tnw, wu, and using the fact that the two-phase saturations add up to 1 in each medium, we obtain the following time translate estimates on φ D and d f,Du :ˇˇx

Compactness properties of Π m
Dp s α m and T a Dp s α a 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 coercivity and compactness properties, and such that lim lÑ`8 ∆t l " 0. Let φ m,min ą 0 and assume that, for each l P N, the gradient scheme (8a)-(8b) has a solution p α l " pp α m,l , p α f,l q 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 φ D l pt, xq ě φ m,min for a.e. pt, xq P p0, T qˆΩ and d f,D l u pt, xq ě d 0 pxq for a.e. pt, xq P p0, T qˆΓ. Then, the sequences pΠ m Proof. The superscript l P N will be dropped in the proof and all hidden constants in the following estimates are independent of l. Setting it results from hypothesis (H2) that S α rt ppq´S α rt ppq¯2 ď´S α rt ppq´S α rt ppq¯´F α ppq´F α pqq¯, for rt P tm,˘u. Using that φ D pt, xq ě φ m,min for a.e. pt, xq P p0, T qˆΩ and noting that Π m Dp s α m " S α m pΠ m Dp p c,m q P r0, 1s and T a Dp s α a " S α a pT a Dp p c,m q P r0, 1s, we obtain with ζ α m ptq "´F α pp c,m pt`τ qq´F α pp c,m ptqq¯and χ α m ptq " ζ α m ptq s α m pt`τ q. Let us set ζ α ptq " pζ α m ptq, 0q P X 0 Dp . In view of the estimates (26) for ϕ " ζ α ptq, we have From Proposition 4.5, we have Using the a priori estimates of Lemma 4.3, h α m P L 2 pp0, T qˆΩq, the Lipschitz property and boundedness of S α m and S α a , the chain rule estimate on the sequence of GDs pD l p q lPN , and the bound on the jump operator, we obtain that We deduce from [5, Lemma 4.1] that T 1 À τ`∆t. Similarly, using the time translate estimate (29) and the product rule estimate on the sequence of GDs pD l p q lPN , one shows that T 2 À τ`∆t, which provides the time translates estimates on Π m Dp s α m in L 2 p0, T ; L 2 pΩqq and on T a Dp s α a in L 2 p0, T ; L 2 pΓqq. Let us consider any compact sets K m Ă Ω and K f Ă Γ. The space translates estimates for Π m Dp s α m in L 2 p0, T ; L 2 pK m qq and for T a Dp s α a in L 2 p0, T ; L 2 pK f qq derive from the a priori estimates of Lemma 4.3, the Lipschitz properties of S α m and S α a , and from the local compactness property of the sequence of spatial GDs pD l p q lPN (cf. Remark 3.2). Combined with the time translate estimates above, the Fréchet-Kolmogorov theorem implies that Π m Dp s α m is relatively compact in L 2 p0, T ; L 2 pK m qq and that T a Dp s α a is relatively compact in L 2 p0, T ; L 2 pK f qq. Since Π m Dp s α m P r0, 1s and T a Dp s α a P r0, 1s it results that Π m Dp s α m is relatively compact in L 2 p0, T ; L 2 pΩqq and that T a Dp s α a is relatively compact in L 2 p0, T ; L 2 pΓqq.

Uniform-in-time
Proposition 4.8. 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 coercivity and consistency properties. Let φ m,min ą 0 and assume that, for each l P N, the gradient scheme (8a)-(8b) has a solution p α l " pp α m,l , p α f,l q P pX 0 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 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,f q, converge up to a subsequence uniformly in time and weakly in L 2 pΓq, as per [22,Definition C.14].
Proof. In the following, the superscript l P N is dropped when not required for the clarity of the proof, and the hidden constants are independent of l. Let ϕ f P C 8 c pΓq and let ϕ f P X 0 D f p the element that realizes the minimum of (10). From Proposition 4.5 (with ϕ m " 0) we havěˇˇx

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.3 and Proposition 4.9, there existd f P L 8 p0, T ; L 4 pΓqq ands α f P L 8 pp0, T qˆΓq such that From Proposition 4.7, there exists α m P L 8 pp0, T qˆΩq ands α a P L 8 pp0, T qˆΓq such that Π m Dp S α m pp c,m q Ñs α m in L 2 p0, T ; L 2 pΩqq, T a Dp S α a pp c,m q Ñs α a in L 2 p0, T ; L 2 pΓqq.
The identification of the limit [15, Proposition 3.1], resulting from the limit-conformity property, can easily be adapted to our definition of V 0 , with weight d 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.3, we obtain p α " pp α m ,p α f q 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 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 On the other hand, it results from (32) 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 f,Du ϕ Ñ pd f q 3 {2 ϕ in L 8 p0, T ; L 2 pΓq d´1 q given by (30), we infer that This shows that g α f " pd f q 3 {2 ∇ τp α f on p0, T qˆΓ. Combining the strong convergence of Π m Dp S α m pp c,m q " S α m pΠ m Dp p c,m q, the weak convergence of Π m Dp p c,m , it results from the Minty trick (see, e.g., [26, Lemma 2.6]) thats α m " S α m pp c,m q with pp c,m ,p c,f q "p c " pp nw m´p w m ,p nw f´p w f q "p nw´pw . Using the same arguments, we also haves α f " S α f pp c,f q ands α a " S α a pγ apc,m q. From the a priori estimates of Lemma 4.3 and the limit-conformity property of the sequence of GDs pD l u q lPN (see [10,Lemma A.3]), there existsū P L 8 p0, T ; U 0 q, such that from which we deduce thatd f "´ ū 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 andp E m P L 8 p0, T ; Since 0 ď U rt ppq " ş p 0 qpS nw rt q 1 pqqdq ď 2|p| for rt P tm, f u, it results from the a priori estimates of Lemma 4.3 that there existp E 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 For rt P tnw, wu, it is shown in [23], following ideas from [21], that U rt ppq " B rt pS nw rt ppqq where s P r0, 1s Þ Ñ B rt psq P p´8,`8s is a convex lower semi-continuous function with finite limits at s " 0 and s " 1 (note that B rt is therefore continuous). Since Π m Dp s nw m converges strongly in L 2 pp0, T qˆΩq to S nw m pp c,m 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,m qq, and hence thatŪ m " B m pS nw m pp c,m qq " U m pp c,m q. Similarly,Ū f " B f pS nw f pp c,f qq " U f pp c,f q. We deduce, using the strong convergences of the saturations and the weak convergences of the pressures, that Using the estimate |U rt pp 2 q´U rt pp 1 q| "ˇˇˇˇż p2 p1 qpS nw rt q 1 pqqdqˇˇˇˇď |p 2´p1 |`|p 2 S nw rt pp 2 q´p 1 S nw rt pp 1 q|, the Lipschitz property of S nw rt ,p α 0 " pp α 0,m ,p α 0,f q P V 0 X pL 8 pΩqˆL 8 pΓqq, α P tnw, wu, and the consistency of the sequence of GDs pD l p q lPN , we deduce that Then, from [10,Proposition A.4] it holds that It results from (33), (34), (36), (37) and the definition of φ D that Let us now prove that the functionsp α , α P tnw, wu, andū satisfy the variational formulation (6a)-(6b) by passing to the limit in the gradient scheme (8).
For θ P C 8 c pr0, T qq and ψ " pψ m , ψ f q P C 8 c pΩzΓqˆC 8 c pΓq let us set, with P Dp ψ " pP m Dp ψ m , P f Dp ψ f q P X 0 Dp and P ν Dp ψ ν realising the minimum of S D ν p pψ ν q, ϕ " pϕ 1 , . . . , ϕ N q P pX 0 Dp q N with ϕ i " pϕ i m , ϕ i f q " θpt i´1 qpP Dp ψq.
Let us set ϕ ν " pϕ 1 ν , . . . , ϕ N ν q, ν P tm, f u. From the consistency properties of pD l p q lPN with given r ą 8, we deduce that Setting the gradient scheme variational formulation (8a) states that 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 realises the minimum in the definition (14) of S Du pwq. From the consistency properties of pD l u q lPN , we deduce that Setting Du puq : ε Du pvq´bpΠ m Dp p E m qdiv Du pvq¯dxdt, the gradient scheme variational formulation (8b) states that Using a discrete integration by part [22, Section D.1.7], we have T 1 " T 11`T12 with Using (38) and (34), and that Π m Dp s α m P r0, 1s converges to S α m pp c,m q a.e. in p0, T qˆΩ (this follows from (31)), it holds that Using (38), that Π m Dp J m Dpφ 0 converges in L 2 pΩq toφ 0 and that Π m Dp S α m pI m Dpp α m,0 q P r0, 1s converges a.e. in Ω to S α m pp α m,0 q, we deduce that Writing T 3 " T 31`T32 with we obtain, using similar arguments and (37), that Writing T a 1 " T a 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,m q, (32) and (38), it holds that themselves as test functions in the weak equations (6a)-(6b) they satisfy and using fine integrating-by-parts in time results from [21].
The caveat here is that space of test functions for (6a), which is C 8 c pr0, T qˆΩzΓqˆC 8 c pr0, T qˆΓq, is not obviously dense in the space of trial functions L 2 p0, T ; V 0 q, in which the limit pressures are found. Hence, it is not clear that we can indeed use these limit pressures as test functions in (6a). The density issue comes from the fact that we would need to find smooth functions pϕ α f,k q kě1 such thatd , T qˆΓq; in other words, we would like smooth functions to be dense in the weighted space H 1 d f pΓq. Such a density result has been established in [33], but under an additional assumption on the weight. Specifically: d f is bounded in time, smooth in space away from the fracture tips, and, near the tips,d f pt, x, yq « x 1 2` f pyq with ą 0 (above, x is the distance to the tip, y is the coordinate parallel to the tip, and f is smooth). Under this assumption, the arguments of [23] can be reproduced and the limit fluxesQ α f,a can be shown to satisfy the first equation in (6c).

Numerical experiments
The objective of this numerical section is to compare the discontinuous pressure poro-mechanical model investigated in this work with the continuous pressure poro-mechanical model presented in [10]. Two test cases are considered. The first one already described in [10] considers the injection of gas in a cross-shaped fracture network coupled with the matrix domain initially liquid saturated. The second test case models the desaturation by suction at the interface between a ventilation tunnel and a low-permeability fractured porous medium.
For both test cases, the flow part of system (3) is discretized in space by a Two-Point Flux Approximation (TPFA) cell-centered finite volume scheme with additional face unknowns at matrix fracture interfaces [2]. The mechanical part of (3) is discretized using second-order finite elements (P 2 ) for the displacement field in the matrix [18,37], adding supplementary unknowns on the fracture faces to account for the discontinuities. The computational domain Ω is decomposed using admissible triangular meshes for the TPFA scheme (cf. [25, Section 3.1.2]) as illustrated in Figure 4.
For the TPFA scheme, the GD operators Π m Dp and Π f Dp are respectively cell-wise and fracture face-wise constant. It results that the porosity and the fracture aperture defined by the closure laws (8c) will be projected in the matrix and fracture accumulation terms (cf. first two terms in (8a)) to cell-wise and face-wise constant spaces respectively. For simplicity, this face-wise constant projection of the fracture aperture is also used in the fracture conductivity.
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.001 days is the initial time step, ∆t max " 10 days in the first test case and 10 years in the second one, 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 [24]). The coupled nonlinear system is then solved at each time step using a Newton-Krylov acceleration [48] of the fixed point algorithm which, for a given displacement field, solves the two-phase Darcy flow problem, then computes the new displacement field given the new equivalent pressures (see [11] for more details). The stopping criterion is fixed to 10´5 on the relative displacement field increment. This Newton-Krylov algorithm is compared in [11] to the fixed stress algorithm [43] extended to DFM models in [31]. It is shown to solve the robustness issue of fixed stress algorithms w.r.t. to small initial time steps in the case of incompressible fluids.
For this section, we introduce the following notation for the total stress: where σ 0 is a possible pre-stress state [47, Section 4.2.4].

Gas injection in a cross-shaped fracture network
The data set of the continuous pressure model is the one described in [10]. We recall it briefly here. We consider the square Ω " p0, Lq 2 lying in the xy-plane, with L " 100 m, containing a cross-shaped fracture network Γ made up of four fractures, each one of length L 8 intersecting at p L 2 , L 2 q and aligned with the coordinate axes. The matrix and fracture network have the following mobility laws: η α m ps α q " ps α q 2 µ α , η α f ps α q " s α µ α , α P tw, nwu, 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; nevertheless, the results of the numerical experiments are not affected by this circumstance. The saturation-capillary pressure relation is Corey's law: 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´1 5 m 2 , an initial porosity φ 0 m " 0.2, effective Lamé parameters λ " 833 MPa, µ " 1250 MPa, Biot's coefficient b " 1´K dr Ks » 0.81, and Biot's modulus M " 18.4 GPa. The pre-stress state is assumed null: σ 0 " 0. 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 upper side of the boundary p0, T qˆpp0, LqˆtLuq, whereas the remaining part of the boundary is considered 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: , pt, xq P p0, T qˆΓ, ż Ω φ 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. To define the discontinuous pressure model, we consider additionally the normal fracture transmissivity T f " 10´8 m.
From Figure 5, it is clear that both the continuous and discontinuous pressure models provides roughly the same solutions. Nevertheless, the continuous pressure model provides a rather smoothed non-wetting phase saturation at matrix fracture interfaces while the discontinuous pressure model is more accurate as discussed in [2]. Also, Figure 6 shows the axial total stresses σ T x , σ T y , and the shear total stress σ T xy in the matrix at the final time, for the discontinuous pressure model; from the mechanics viewpoint, in this case the difference between the two models is not so remarkable. As expected, stresses are concentrated in the neighborhood of fracture tips. To conclude this subsection, we give an insight into the performance of our method in Table 1 Table 1: Performance of the method for the plane problem, in terms of the number of mesh elements, the number of successful time steps, the number of time step chops, the total number of Newton-Raphson iterations, the total number of GMRes iterations, the total number of Newton-Krylov iterations, and the total computational time.  Figure 6. Axial normal stresses σ T x and σ T y (Pa) and shear stress σ T xy (Pa) at final time for the discontinuous pressure model.

Desaturation by suction of a low-permeability fractured porous medium
In this test case, we consider a hollow cylinder (Figure 7) made up of a low-permeability porous medium, containing an axisymmetric fracture network, subject to axisymmetric loads -uniform pressures exerted on the internal and external surfaces. Using cylindrical coordinates px, r, θq, the problem can therefore be reduced to a two-dimensional formulation on the diametral section of the medium, shown in Figure 8 along with the fracture network, and the displacement field only consists of its axial and radial components: upx, r, θq " u x px, rqe x`ur px, rqe r pθq, e r pθq " pcos θqe y`p sin θqe z , where we have dropped time dependence for simplicity, and taken into account the system of cylindrical coordinates in Figure 7, denoting by e x the axial unit vector, by e r " e r pθq the radial unit vector, and by e θ the orthoradial unit vector. The final time for this simulation is set to T " 200 years. The geometry is characterized by the following data set: length L " 10 m, internal and external radii R int " 5 m, R ext " 35 m; two consecutive fractures are spaced by 1.25 m. The matrix is characterized by the Lamé parameters λ " 1.5 GPa, µ " 2 GPa, by a permeability Λ m " 5¨10´2 0 m 2 , Biot's coefficient and modulus b " 1 and M " 1 GPa respectively, and by an initial porosity φ m " 0.15. The normal transmissibility of fractures is T f " 10´9 m and the initial fracture aperture is set to 10´2 m. The matrix relative permeabilities of the liquid and gas phases are defined by the following Van Genuchten laws: k nw r,m ps nw q " withs w " s w´S lr 1´S lr´Sgr , and the parameter q " 0.328, the residual liquid and gas saturations S lr " 0.40 and S gr " 0; in the fractures, we take k α r,f psq " s for both phases. The phase mobilities are then η α m ps α q " k α r,m ps α q{µ α and η α f ps α q " k α r,f ps α q{µ α , α P tw, nwu both in the matrix and in the fractures, with the same viscosities as in the previous test case. Again, η α m and η α f are not bounded below by a strictly positive number, but this does not have an influence on the numerical results. The saturation-capillary pressure relation is again Corey's law, as in (41), with R m " 2¨10 8 Pa and R f " 10 2 Pa. Moreover, the medium is supposed to be pre-stressed with the following pre-stress state:   Table 2: Performance of the method for the axisymmetric problem, in terms of the number of mesh elements, the number of successful time steps, the number of time step chops, the total number of Newton-Raphson iterations, the total number of GMRes iterations, the total number of Newton-Krylov iterations, and the total computational time.
Full saturation of the liquid phase is assumed at the initial state, both in the matrix and in the fracture network, with an initial uniform pressure p w 0 " p nw 0 " 4 MPa. Concerning flow boundary conditions, the porous medium is assumed impervious (vanishing fluxes) on the lateral boundaries corresponding to x " 0 and x " L. On the inner surface r " R int , a given gas saturation is imposed: s nw m " 0.35 on the matrix side and s nw f " 1´10´8 at fracture nodes, and atmospheric pressure p atm " 10 5 Pa everywhere. On the outer surface r " R ext , a liquid saturation s w m " 1 and pressure p w " 4 MPa are imposed. As for the mechanical boundary conditions, we impose a vanishing axial displacement u x on the lateral boundaries corresponding to x " 0 and x " L. Moreover, on the same boundaries, the tangential stress is set to zero. On the other hand, external surface loads g (uniform pressures) are applied on the inner and outer surfaces: where n " e r for r " R ext and n "´e r for r " R int . We consider σ T N " 10.95 MPa as the numerical value for the uniform pressure on the outer surface.
As shown in Figures 9 and 10 strong capillary forces induce the desaturation of the matrix in the neighborhood of the inner surface combined with a high negative liquid pressure. As exhibited in Figures 11 and 12 this negative liquid pressure triggers the contraction of the pores as well as the spreading of the fracture sides. Figure 9 also displays a comparison between the matrix non-wetting saturations obtained with the discontinuous and continuous pressure models at final time. It can be clearly seen that, unlike the continuous pressure model, the discontinuous pressure model is able to capture the barrier effect induced on the liquid phase by the fractures almost fully filled by the gas phase. This is particularly remarkable at the intersection of the horizontal and oblique fractures. Figure 10 shows a comparison of matrix equivalent pressures at final time obtained for the continuous and discontinuous pressure models; in the first case, discontinuities at the matrix-fracture interface can be clearly detected. Figure 12 shows the time history of the average fracture aperture for the continuous and discontinuous pressure models, with significant differences induced by the equivalent pressures p E m computed in the two models. Finally, in Figure 13 we display the radial, orthoradial, axial, and shear total stresses σ T r , σ T θ , σ T x , and σ T rx respectively in the matrix at the final time for the discontinuous pressure model. Again, stresses are concentrated in the neighborhood of fracture tips, as expected. The arching effect is clearly visible by comparison of the radial and orthoradial stresses in the neighborhood of the inner surface. As expected, the radial stresses are transmitted across the horizontal fracture as opposed to the orthoradial stresses. The comparison of the results given by the two models is shown in Figure 14, where a different behavior in the orthoradial total stresses σ T θ given by the two models along the vertical line x " 5.5 (intersecting the horizontal fracture) can be detected.
As in the previous subsection, we summarize also here the performance of our method in Table 2.

Conclusions
This work extends the gradient discretization and convergence analysis carried out in [10] to the case of hybriddimensional poro-mechanical models with discontinuous phase pressures at matrix fracture interfaces. The model considers a linear elastic mechanical model with open fractures coupled with a two-phase Darcy flow. The Poiseuille law is used for the tangential fracture conductivity and the dependence of the normal fracture transmissivity on the fracture aperture is frozen. The model accounts for a general network of planar fractures including immersed, nonimmersed fractures and fracture intersections, and considers different rock types in the matrix and fracture network domains as well as at the matrix fracture interfaces. Figure 7. Hollow cylinder of length L and internal and external radii R int and R ext , respectively. The diametral section is highlighted in gray, the fracture network is not shown for simplicity.       Two test cases were considered to compare the continuous pressure hybrid-dimensional poro-mechanical model investigated in [10] to the discontinuous pressure model studied in this work. The first test case simulates the gas injection in a cross-shaped fracture network immersed in a two-dimensional porous medium initially water saturated. The second test case is based on an axisymmetric DFM model and simulates the desaturation by suction at the interface between a ventilation tunnel and a Callovo-Oxfordian argilite fractured storage rock. In both cases, it is shown that the discontinuous pressure model provides a better accuracy at matrix fracture interfaces than the continuous pressure model and allows in particular to account for the barrier effect induced on the liquid phase by the gas filled fractures.