Jakub W Both

C L Ément

Canc Ès

A GLOBAL EXISTENCE RESULT FOR WEAKLY COUPLED TWO-PHASE POROMECHANICS

Keywords: 2010 Mathematics Subject Classification. 76S05, 35M33, 35K65 Poromechanics, two-phase porous media flows, global existence . 1

Multiphase poromechanics describes the evolution of multiphase flow in deformable porous media. Mathematical models for such multiphysics system are inheritely nonlinear, potentially degenerate and fully coupled systems of partial differential equations. In this work, we present a thermodynamically consistent multiphase poromechanics model falling into the category of Biot equations and obeying to a generalized gradient flow structure. It involves capillarity effects, degenerate relative permeabilities, and gravity effects. In addition to established models it introduces a Lagrange multiplier associated to a bound constraint on the effective porosity in particular ensuring its positivity. We establish existence of global weak solutions under the assumption of a weak coupling strength, implicitly utilizing the gradient flow structure, as well as regularization, a Faedo-Galerkin approach and compactness arguments. This comprises the first global existence result for multiphase poromechanics accounting for degeneracies that are consistent with the multiphase nature of the flow.

1. Two-phase poromechanics: motivation, formulation and main result 1.1. Introduction. Poromechanics, modeled by the prototypical Biot equations, describes the twoway coupled interaction between flow in porous media and its macro-and microscopic deformation. Since the seminal works in mathematical modelling by von Terzaghi [START_REF] Terzaghi | Erdbaumechanik auf bodenphysikalischer Grundlage[END_REF] and Biot [START_REF] Biot | General theory of three-dimensional consolidation[END_REF], the study of poromechanics has been extended to various applications across engineering [START_REF] Coussy | Poromechanics[END_REF]. With relevance ranging from geotechnical, environmental, over industrial to biomedical engineering, poromechanics plays a paramount role.

In many engineering systems, as e.g. subsurface reservoirs for geological CO 2 storage or geothermal energy, the simultaneous presence of multiple fluids introduces a range of physical phenomena which requires nonlinear and degenerate mathematical models. To capture particular multiphase effects as varying saturations, capillary pressure, drying shrinkage, pore pressure changes and finally related deformation and subsidence, the classical linearized Biot equations are not sufficient and require the introduction of a range of nonlinear constitutive relations. For this, mathematical models for multiphase poromechanics have been introduced combining separate components of multiphase flow and poroelasticity modeling [START_REF] Coussy | Poromechanics[END_REF]. The general model structure has been also studied from a thermodynamical point of view [START_REF] Seguin | Multi-component multiphase flow through a poroelastic medium[END_REF]. Furthermore, there exist many works on the numerical approximation of such systems, e.g., [START_REF] Bonaldi | Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media[END_REF][START_REF] Both | Anderson accelerated fixed-stress splitting schemes for consolidation of unsaturated porous media[END_REF] and they are employed by practitioners [START_REF] Jha | Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering[END_REF].

Despite an increased interest in the modelling and numerics communities, multiphase poromechanics has been studied in much less detail and with less rigor than the much simpler linearized Biot equations for single fluids and several extensions -no well-posedness results exist for twophase poromechanics models. In contrast, in the seminal work [START_REF] Auriault | étude du comportement macroscopique d'un mileu poreux saturé deformable[END_REF], the first well-posedness results for the linear Biot equations have been established, deriving existence of strong solutions. A series of works has followed establishing weak solution [START_REF] Ženíšek | The existence and uniqueness theorem in Biot's consolidation theory[END_REF], well-posedness from the perspectives of semi-group theory [START_REF] Showalter | Diffusion in poro-elastic media[END_REF] as well as generalized gradient flow theory [START_REF] Both | The gradient flow structures of thermo-poro-viscoelastic processes in porous media[END_REF]. The well-posedness for mixed-dimensional extensions to fractured media under contact have been established in [START_REF] Boon | Mixed-dimensional poromechanical models of fractured porous media[END_REF] utilizing monotone operator theory. Additional physics have been considered and rigorously investigated, as linear poroviscoelastic effects [START_REF] Bociu | Mathematical effects of linear visco-elasticity in quasi-static Biot models[END_REF] as well as the coupling with Stokes' equations in a classical fluidstructure-interaction fashion [START_REF] Ambartsumyan | A nonlinear stokes-biot model for the interaction of a non-newtonian fluid with poroelastic media[END_REF][START_REF] Bociu | Multilayered poroelasticity interacting with Stokes flow[END_REF]. Analytical results for nonlinear extensions include single-phase flow accounting for compressibility and positivity-preserving porosity [START_REF] Van Duijn | Mathematical theory of nonlinear single-phase poroelasticity[END_REF], thermal effects [START_REF] Van Duijn | Thermoporoelasticity via homogenization: modeling and formal two-scale expansions[END_REF][START_REF] Brun | Well-posedness of the fully coupled quasi-static thermo-poroelastic equations with nonlinear convective transport[END_REF], as well as displacement-depending permeability laws [START_REF] Bociu | Analysis of nonlinear poro-elastic and poro-visco-elastic models[END_REF]. Finally, closely related to this study, it is worth stressing, that the existence of weak solutions for unsaturated poromechanics models has been established in [START_REF] Showalter | Partially saturated flow in a poroelastic medium[END_REF][START_REF] Both | Global existence of weak solutions to unsaturated poroelasticity[END_REF], where unsaturated media unlike true multiphase systems are partially saturated by one active fluid governing the displacement of a second passive fluid.

The model we consider here accounts for the motion of two incompressible phases in a porous medium, modeled by the Darcy-Muskat law and some capillary pressure law, see e.g. [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF]. In particular, the system degenerates as one fluid phase vanishes, leading to mathematical difficulties making suitable mathematical reformulations involving the Kirchhoff transform and Chavent global pressure necessary (see for instance [START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution[END_REF]). The porous matrix is supposed to be elastic (small deformation is postulated) and interacts with the fluid by the mean of the equivalent pressure postulated by Coussy [START_REF] Coussy | Poromechanics[END_REF]. As we restrict to linear relative permeabilities, the Coussy equivalent pressure coincides with the Chavent global pressure. Another important difficulty comes from the fact that linear elasticity does not prevent the porosity to become nonpositive. To maintain the model in its regime of validity, we incorporated some Lagrange multiplier χ, so that we do not have to assume the porosity to be positive as for instance in [START_REF] Bonaldi | Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media[END_REF]. Furthermore, we allow for the permeability to depend on the porosity in the line of Kozeny [40] and Carman [START_REF] Carman | Fluid flow through granular beds[END_REF]. Finally, state-dependent gravitational forces are incorporated acting both on the bulk as well as on the fluid phases. To the best of the authors' knowledge, the presented result is the first-ever existence result for multiphase poromechanics in presence of physically relevant degenerate mobilities. Our results comes however with restrictions. Constant Biot coefficients and Biot modulus are assumed for simplicity, as well as specific (but physically meaningful) boundary conditions on the displacement. We also require some weak-coupling condition. Assumptions are detailed and discussed in Section 1.4 later on. The extension of existence of weak solutions for tightly coupled and/or heterogeneous in space systems with discontinuous characteristics remains an open problem.

1.2. The governing equations. Let us first start by stating the equations governing the motions of two immiscible fluids -a wetting one labeled with subscript w and a non-wetting one labeled by n -in a porous medium represented by some bounded open set Ω ⊂ R d . We remain sloppy here concerning regularity, which will be made precise later on. For α ∈ {n, w}, denote by ρ α and µ α the (constant) density and viscosity of the phase α, and by g the gravity vector (pointing downwards). In addition, as in the single-phase Biot model [START_REF] Biot | General theory of three-dimensional consolidation[END_REF], the porous matrix is assumed to be deformable; let u denote the macroscopic displacement of the matrix wrt. Ω. The classical multiphase Darcy law, with the hydrostatic phase pressure extended to the deformable case in a thermodynamically consistent way, then writes (1a)

∂ t φ α -∇ • s α µ α K(φ)∇ (p α -ρ α g • (x + u)) = 0
where, for volume densities φ = (φ n , φ w ) of the two phases, we denote by φ = φ n + φ w the porosity, and by s α = φα φ the saturation (or volume fraction) of the phase α. The intrinsic permeability K of the porous medium may depend on the porosity φ, whereas we restrict our purpose to linear relative permeabilities for technical reasons that will appear in what follows. The pressure of the phase α ∈ {n, w} decomposes into three contributions:

(1b) p α = pα + π + χ.

In the above right-hand side, only the first term pα depends on α. It is related to φ, and more specifically to the saturations, via the formulas (1c) pn = γ(s n ) + s w γ ′ (s n ) and pw = γ(s n )s n γ ′ (s n )

where γ : [0, 1] → R + is convex and increasing, and shall be interpreted as the antiderivative of the capillary pressure function. The last term χ, the introduction of which being somehow artificial, shall be thought as a Lagrange multiplier ensuring that the porosity will not leave the range of validity of the model φ

∈ [φ ♭ , φ ♯ ] for some constants 0 < φ ♭ < φ ♯ < 1. Therefore, χ is related to φ through the maximal monotone graph χ = ∂1 [φ ♭ ,φ ♯ ] by (1d) χ ∈ χ(φ) with χ(φ) =      0 if φ ♭ < φ < φ ♯ , (-∞, 0] if φ = φ ♭ , [0, +∞) if φ = φ ♯ .
Note in particular that since s n + s w = 1, one recovers the capillary pressure relation

p n -p w = γ ′ (s n ).
Moreover, we deduce from (1c) that s n pn + s w pw = γ(s n ), so that (1b) yields

s n p n + s w p w -γ(s n ) = π + χ.
In particular, on the set {φ ♭ < φ < φ ♯ }, where the Lagrange multiplier χ vanishes, π coincides with the equivalent pressure introduced by Coussy [START_REF] Coussy | Poromechanics[END_REF], and is therefore referred to as the Coussy pressure. It encodes the pressure felt by the porous matrix surrounding the fluid.

We assume small displacements u, describing the macroscopic deformation of the bulk, so that we can stick to linear elasticity. Moreover, we assume the system to be quasi-static, i.e. the mechanical response of the matrix is supposed to be instantaneous, so that mechanical forces remain at equilibrium. Denoting by σ the Cauchy stress tensor, decomposing into effective stress σ and pore pressure contribution π [START_REF] Biot | General theory of three-dimensional consolidation[END_REF][START_REF] Coussy | Poromechanics[END_REF] (2)

σ = σ -bπI,
and by f a body force, the balance of linear momentum

-∇ • σ = f writes (3a) ∇ • σ = b ∇π -f ,
where b ∈ (0, 1] is the so-called Biot coefficient. The effective stress tensor σ and the strain tensor ε(u) = (∇u + ∇u T )/2 are related through the isotropic Hooke law

(3b) σ = 2µ ε(u) + λ∇ • u I,
where λ, µ > 0 are the Lamé coefficients. The body force f may depend on the fluid distribution, e.g., the gravitational force typically incorporates an effective volume-averaged bulk density

(4a) f g (φ) =   φ α∈{n,w} s α ρ α + (1 -φ)ρ s (φ)   g =   α∈{n,w} φ α ρ α + (1 -φ r )ρ s,r   g,
where ρ s denotes the rock density set to be a function of the porosity, defined through

(4b) (1 -φ) ρ s (φ) = (1 -φ r ) ρ s,r ,
and thus satisfying satisfying mass conservation d dt E (1φ)ρ s (φ) = 0 for any measurable subset E of Ω, cf. [START_REF] Coussy | Poromechanics[END_REF]. From now on, we consider for f

(4c) f (φ) = f g (φ) + f ext
for some additional external, state-independent body force f ext .

Beside the macroscopic deformation of the matrix encoded by its displacement, the porous structure is assumed to be compressible. This is encoded by the parameter θ representing the microscopic compression of the solid grains of the porous structure. It simply relates to the Coussy pressure by [START_REF] Auriault | étude du comportement macroscopique d'un mileu poreux saturé deformable[END_REF] M θ = π with M > 0 being referred to in the literature as the Biot modulus.

The last equation set on the time-and-space bulk R + × Ω is a constraint on the fact that the fluid and the solid have to share the available space, leading to

(6) φ -b ∇ • u -θ = φ r
for some spatially varying reference porosity with values φ r (x) ∈ (φ ♭ , φ ♯ ), x ∈ Ω, representing the porosity at rest. To close the system, we prescribe some initial conditions ( 7)

φ |t=0 = (φ n| t=0 , φ w |t=0 ) = φ 0 = (φ 0 n , φ 0 w ) on the fluid contents with φ ♭ ≤ φ 0 = φ 0 n + φ 0 w ≤ φ ♯ .
Since the mechanical equilibration is instantaneous, the initial displacement u 0 and microscopic deformations θ 0 are derived from φ 0 as the minimizer of some mechanical energy under the constraint [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF], see (H4) later on. Boundary conditions of mixed type are considered. More precisely, given a partition Γ N , Γ D of ∂Ω, and denoting by n the normal to ∂Ω outward w.r.t. Ω, we prescribe

(8a) u = 0 on R + × Γ D and σ • n = b πn on R + × Γ N
for the solid mechanics equations [START_REF] Kim | Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits[END_REF], i.e., homogeneous boundary conditions for displacement and traction in terms of the effective stress, while the fluid part (1) is complemented by setting

(8b) - s α µ α K(φ)∇ (p α -ρ α g • (x + u)) • n = 0 on R + × Γ N , and (8c)
p n = p D n and p w = p D w on R + × Γ D .
Remark 1.1 (Non-homogeneous Dirichlet data for the displacement). The case of non-homogeneous Dirichlet data u D for the displacement in (8a) can be reduced in a standard way to the case of homogeneous data, considered above. Indeed, in the momentum balance equation (3a), using a lifting, its contribution can be incorporated in the constant, external body force f ext , whereas, in the porosity constraint [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF], the contribution is included in the definition of the reference porosity φ r , consistently with its character representing the medium at rest.

1.3.

A thermodynamic viewpoint. The equations presented in Section 1.2 have strong connections with thermodynamics. In particular, the dynamics prescribed by ( 1)-( 8) can be interpreted as some generalized gradient flow of the Helmholtz free energy augmented with some potential energy related to body forces and in particular to gravity. More precisely, to a set X = (φ, ε, θ) of primary unknowns (here ε = ε(u)), we associate the Helmholtz free energy defined as the sum of three contributions ( 9)

F (X) = F f (φ) + F s (ε, θ) + F c (X).
The energy F f (φ) associated to the fluid is defined as

(10) F f (φ) = Ω φγ(s n ) + 1 [φ ♭ ,φ ♯ ] (φ) ,
where

1 C (v) = 0 if v ∈ C and +∞ otherwise.
The mechanical energy associated to the porous matrix deformation is given by ( 11)

F s (ε, θ) = Ω M 2 θ 2 + µ ε : ε + λ 2 | Tr ε| 2 with Tr ε = ∇ • u.
The last term F c (X) enforces the constraint (6) to hold (almost) everywhere in Ω. We introduce the (linear thus

) convex set K = {X | φ -b Tr ε -θ = 1 -φ ⋆ r }, then we set (12) F c (X) = Ω 1 K (X) = sup w Ω (φ -φ r -b Tr ε -θ) w.
In addition, the potential energy of the bulk associated to mechanical loading, is given by ( 13)

F g (φ, u) = - Ω f • (x + u) = - Ω α φ α ρ α + (1 -φ r )ρ s,r g + f ext • (x + u) .
Note that F g is neither convex nor concave, but smooth as the sum of linear and quadratic terms. The Helmholtz free energy is a convex yet non-smooth function of X because of the constraints φ ♭ ≤ φ ≤ φ ♯ incorporated in the definition of F s and the constraint (6) corresponding to the term F c . As a consequence, its subdifferential is not single valued. In particular, Y = (p, σ, z) belongs to ∂F (X) if p = (p n , p w ) is given by (1b

)(1c)&(1d), if σ = 2µ ε + λ Tr ε I -b πI = σ -b πI and if z = M θ -π.
In the previous expressions, the Coussy pressure π can be interpreted as the Lagrange multiplier for the constraint X ∈ K. For the sum of both the Helmholtz and the gravitational energy, the hydrostatic phase pressures p αρ α g • (x + u), acting as fluid potential in (1a), is derived as as an element of the subdifferential ∂ φα (F + F g ).

Assume now that t → X(t) satisfies the systems (1)-( 8), then

(14a) d dt F (X) = Ω α∈{n,w} p α ∂ t φ α + σ : ∂ t ε + z∂ t θ , and
(14b) d dt F g (φ, u) = - Ω f • ∂ t u - Ω α∈{n,w} (ρ α g • (x + u)) ∂ t φ α .
First, z = 0 as a direct consequence of (5), while -∇ • σ = f owing to [START_REF] Kim | Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits[END_REF]. Therefore, Stokes' theorem provides

Ω σ : ∂ t ε - Ω f • ∂ t u = - ∂Ω ∂ t u • σn (8a) = 0.
Besides, it follows from (1a) together with Stokes' theorem and (8) that Ω α∈{n,w}

(p α -ρ α g • (x + u)) ∂ t φ α = - Ω α∈{n,w} s α µ α K(φ) 1/2 ∇ (p α -ρ a g • (x + u)) 2 + Σ D , with Σ D = Γ D α∈{n,w} p D α -ρ α g • x s α µ α K(φ)∇ (p α -ρ α g • (x + u)) • n
being the work of the force imposed on the fluid at the level of the Dirichlet boundary condition (8c) and hydrostatic phase pressures. Therefore, ( 14) yields

(15) d dt (F (X) + F g (φ, u)) = - Ω α∈{n,w} s α µ α K(φ) 1/2 ∇ (p α -ρ a g • (x + u)) 2 + Σ D ,
The first contribution in the right-hand side of ( 15) is nonpositive since s α ≥ 0 and K(φ) is symmetric definite positive. It encodes the entropy production of the system, in terms of the hydrostatic phase pressures. Note that the evolution of the solid is assumed to be reversible, in the sense that no entropy is produced by the mechanical deformation of the porous matrix.

1.4. Weak formulation and main result. The analysis we propose in the next section for system ( 1)-( 8) strongly builds on the stability estimate [START_REF] Both | Global existence of weak solutions to unsaturated poroelasticity[END_REF], and therefore on some mathematical counterparts of the second principle of thermodynamics. More precisely, system (1)-( 8) can be interpreted as a non-autonomous generalized gradient flow (see for instance [START_REF] Mielke | A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems[END_REF][START_REF] Peletier | Variational modelling: Energies, gradient flows, and large deviations[END_REF]). The particular structure of the model under consideration yields several difficulties. A first one comes from the fact that the equations governing the fluid flow and the solid deformation are coupled through the strong constraint [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF]. A second difficulty comes from the degeneracy of the dissipation term (mainly in the fluid phase pressure contribution)

(16) D = Ω α∈{n,w} s α µ α K(φ)∇p α • ∇p α ≥ 0.
Besides the lack of dissipation for the solid part coming from the reversibility pointed out above, another degeneracy comes from the fact that the prefactor s α in D vanishes in regions where only one fluid phase is present. As a consequence, only a weak control on the phase pressures can be deduced from the control of D, as usual in the two-phase setting. It motivates the introduction of the Kirchhoff transform to carry out the mathematical analysis. Let ξ, ψ : [0, 1] → R be the continuous and increasing functions respectively defined by

(17) ξ(s) = s 0 z(1 -z)γ ′′ (z)dz and ψ(s) = s 0 z(1 -z)γ ′′ (z)dz, s ∈ [0, 1],
then one readily deduces from (1c) that ( 18)

s n ∇p n = ∇ψ(s n ) = √ s n s w ∇ξ(s n ), s w ∇p w = -∇ψ(s n ),
and that

(19) s n |∇p n | 2 + s w |∇p w | 2 = |∇ξ(s n )| 2 .
Therefore, provided K(φ) ≥ K ♭ I in the sense of symmetric matrices for some K ♭ > 0, we have

D ≥ K ♭ µ ♯ Ω α∈{n,w} s α |∇p α | 2 = K ♭ µ ♯ Ω |∇ξ(s n )| 2 + |∇(π + χ)| 2
where we have set µ ♯ = max{µ n , µ w }.

To properly define the notion of weak solution, we still have to introduce some relevant functional spaces. In what follows, we denote by H 1 (Ω) the usual Sobolev space equipped with the norm

v 2 H 1 (Ω) = v 2 L 2 (Ω) + ∇v L 2 (Ω) d . We also denote by V = {v ∈ H 1 (Ω) s.t. v | Γ D = 0},
equipped with v V = ∇v L 2 (Ω) d which defines a norm thanks to the Poincaré inequality, and by V ′ its topological dual. The d-dimensional product space of V is denoted by V d . We also define the closed subspace U of V d as [START_REF] Buzzi | Interface conditions for degenerate two-phase flow equations in one space dimension[END_REF][START_REF] Cancès | Two-phase flows involving capillary barriers in heterogeneous porous media[END_REF][START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF], the inverse capillary function (γ ′ )

U = {u ∈ V d s.t. ∇ • u ∈ H 1 (Ω)}, equipped with the norm u 2 U = u 2 V d + ∇(∇ • u)
-1 is extended into a continuous function on the whole R by ( 21)

S(p) =      0 if p ≤ γ ′ (0) (γ ′ ) -1 (p) if γ ′ (0) ≤ p ≤ γ ′ (1) 1 if p ≥ γ ′ (1).
As ξ • S is 1 2 -Lipschitz continuous, then ξ(s D n ) also belongs to H 1 (Ω). We end up with the following definition of a weak solution.

Definition 1.2. A set of functions (φ, u, θ, χ, π) with

• φ α ∈ L ∞ (R + × Ω; R + ) for α ∈ {n, w} satisfying φ ♭ ≤ φ ≤ φ ♯ a.e. in R + × Ω, such that ξ(s n ) -ξ(s D n ) ∈ L 2 loc (R + ; V ), • u ∈ L 2 loc (R + ; U ) ∩ L ∞ loc (R + ; V d ), • θ, π, χ ∈ L 2 loc (R + ; H 1 (Ω)) with χ ∈ χ(φ)
a.e. in R + × Ω, is said to be a global in time weak solution to the problem (1)-( 8) if (5) and (6) hold almost everywhere in R + × Ω, and if

(22) (0,T )×Ω 2µ ε(u) : ε(v) + λ(∇ • u)(∇ • v) = (0,T )×Ω b π ∇ • v + (0,T )×Ω f (φ) • v for all v ∈ L 2 ((0, T ); V d ), and if (23) R+ Ω φ n ∂ t v + Ω φ 0 n v(0, •) + R+ Ω 1 µ n K(φ) (∇ψ(s n ) + s n ∇ ((π + χ) -ρ n g • (x + u))) • ∇v = 0, and (24)
R+ Ω

φ w ∂ t v+ Ω φ 0 w v(0, •)+ R+ Ω 1 µ n K(φ) (-∇ψ(s n ) + s w ∇ ((π + χ) -ρ w g • (x + u)))•∇v = 0,
for all v ∈ C 1 (R + ; V ), such that there exists some T > 0 such that v(t, •) = 0 for all t ≥ T .

Our main result is the existence of such a global in time weak solution under the assumptions we list below.

(H1) The viscosities µ n , µ w and the densities ρ n , ρ w are positive constants, whereas g ∈ R 3 is constant. The Lamé coefficients λ, µ are positive constants, as well as the Biot modulus M . The Biot coefficient b is also constant and belongs to (0, 1]. The porosity at rest φ r is assumed to have regularity φ r ∈ H 1 (Ω) with values in [φ ♭ , φ ♯ ] for some constants 0

< φ ♭ < φ ♯ < 1. (H2) The function γ ∈ C([0, 1]; R + ) ∩ C 1 ([0, 1)
) is strictly convex and increasing, with s → √ 1s γ ′′ (s) belonging to L 1 (0, 1). It is extended into a lower-semicontinuous convex function γ : R → [0, +∞] by setting γ(s

) = +∞ is s / ∈ [0, 1]. (H3) The intrinsic permeability function K belongs to C([φ ♭ , φ ♯ ]; S ++ d (R)).
In particular, there exist

K ♯ ≥ K ♭ > 0 such that K ♭ I ≤ K(φ) ≤ K ♯ I for all φ ∈ [φ ♭ , φ ♯ ] in the sense of symmetric definite matrices. (H4) The initial fluid content φ 0 = φ 0 n , φ 0 w belongs to L ∞ (Ω; R 2 + ) with φ ♭ ≤ φ 0 n + φ 0 w = φ 0 ≤ φ ♯ almost everywhere in Ω. Besides, the initial displacement u 0 ∈ H 1 (Ω) d , with u 0 ∈ V d ,
and the microscopic compression θ 0 ∈ L 2 (Ω) are the unique solution to the elliptic problem

(u 0 , θ 0 ) = argmin (u,θ) s.t. (6) Ω M 2 θ 2 + µ ε(u) : ε(u) + λ 2 |∇ • u| 2 -f (φ 0 ) • u . It is characterized by φ 0 -b ∇ • u 0 -θ 0 = φ r ,
as well as We assume that Γ D has positive (d -1)-dimensional Hausdorff (or Lebesgue) measure, yielding some Poincaré inequality, as well as Korn's inequality. Moreover, we assume that there exists C 1 depending only on Ω such that the unique solutions

Ω 2µ ε(u 0 ) : ε(v) + λ(∇ • u 0 )(∇ • v) = Ω bπ 0 ∇ • v + Ω f (φ 0 ) • v, ∀v ∈ V d , with π 0 = M θ 0 . ( H5
v 1 , v 2 ∈ V d to ∇ • ε(v 1 ) + λ ∇ • v 1 I = w 1 in Ω, ε(v 1 ) + λ ∇ • v 1 I • n = 0 on Γ N , (25a) ∇ • ε(v 2 ) + λ ∇ • v 2 I = ∇w 2 in Ω, ε(v 2 ) + λ ∇ • v 2 I • n = w 2 n on Γ N , (25b) with λ > 0 and w 1 ∈ L 2 (Ω) d and w 2 ∈ H 1 (Ω) satisfy v 1 , v 2 ∈ U and λ ∇(∇ • v 1 ) L 2 (Ω) d ≤ C 1 w 1 L 2 (Ω) , (25c) λ ∇(∇ • v 2 ) L 2 (Ω) d ≤ C 1 w 2 H 1 (Ω) . (25d)
(H7) We assume that the coefficients of the problem satisfy the following weak coupling condition:

λ > M b 2 C 1 .
where C 1 is the constant appearing in Assumption (H6).

(H8) The body force f is of the form prescribed by [START_REF] Andreianov | A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs[END_REF]. We assume moreover that the external body force satisfies f ext ∈ L 2 (Ω) and the reference density of the rock ρ s,r ∈ L ∞ (Ω) does not depend on time.

The above assumptions deserve some comments. First (H1) requires the domain to be homogeneous in space. The extension to the case of heterogeneous porous media would of course be of great interest.

Rather that prescribing the capillary pressure function, we prescribe its antiderivative γ in (H2), the interpretation of which in terms of energy being a cornerstone of [START_REF] Cancès | The gradient flow structure of immiscible incompressible two-phase flows in porous media[END_REF][START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF], see also Section 1.3. The setting we study does not allow the capillary energy density function γ to depend on the porosity φ, as suggested in the seminal work of Leverett [41]. This choice has been made to stick to the framework of Coussy [START_REF] Coussy | Poromechanics[END_REF]. Extending our result to the case where γ also depends on φ wouldn't lead to major difficulties provided F f in (9) remains convex. However, our framework already encompasses classical models from the literature, as for instance the Brooks-Corey model. In the later, γ(s) ∼ (1s)

1-1 λ BC satisfies (H2) for λ BC > 2, corresponding to a so-called narrow pore size distribution. Note that Assumption (H2) can be easily relaxed by assuming merely that s → (1s) γ ′′ (s) belongs to L 1 (0, 1), at the price of a slightly weaker regularity requirement in Definition

1.2, that is ψ(s n ) -ψ(s D n ) ∈ L 2 loc (R; V ) instead of ξ(s n ) -ξ(s D n ) ∈ L 2 loc (R; V ).
The proof, which can be readily completed by passing to the limit in yet another step of regularization, is left to the reader. Under such a relaxed assumption, the full range λ BC > 1 of Brook-Corey exponents can be recovered.

Assumption (H3) gives a generic framework for the dependance of the permeability with respect to the porosity. This framework encompasses the classical models by Kozeny [START_REF] Kozeny | Uber kapillare Leitung der Wasser in Boden[END_REF] and Carman [START_REF] Carman | Fluid flow through granular beds[END_REF], but also more recent models [START_REF] Costa | Permeability-porosity relationship: A reexamination of the Kozeny-Carman equation based on a fractal pore-space geometry assumption[END_REF][START_REF] Schulz | Beyond Kozeny-Carman: predicting the permeability in porous media[END_REF].

As the mechanical response of the porous matrix is instantaneous, it is natural to require the initial data u 0 and θ 0 to be at equilibrium with the fluid distribution, the later being of finite energy. This is the purpose of Assumption (H4).

Assumption (H6) looks reasonable as it extends to the case of more general boundary conditions a results which is known to hold true for convex, polygonal two-dimensional domains Ω, cf. Brenner and Sung [START_REF] Brenner | Linear finite element methods for planar linear elasticity[END_REF]Section 2] in the pure traction or pure displacement regimes. Here, it is here merely extended to the case of more general boundary conditions of mixed type. Note however that even in the simpler case of the Laplace equation, this may lead to geometrical constraints on the splitting of ∂Ω into Γ D and Γ N , see for instance [START_REF] Grisvard | Singularités en elasticité[END_REF]Section 6.2]. Note also that the full H 1 (Ω) norm appears in the right-hand side of (25d) since a constant w 2 in (25b) possibly yields a non-constant solution v 2 due to the boundary condition on Γ N . Assumption (H5) is there for the sake of simplicity. In the case of time varying boundary conditions on the pressures, suitable regularity assumptions are needed, as for instance in [START_REF] Cancès | Convergence and a posteriori error analysis for energy-stable finite element approximations of degenerate parabolic equations[END_REF].

Assumption (H7) is known in the literature as a weak coupling condition. This regime is realistic in many applications with low Skempton coefficient [START_REF] Coussy | Poromechanics[END_REF][START_REF] Ulm | Is concrete a poromechanics materials?-a multiscale investigation of poroelastic properties[END_REF], and a similar condition appears in papers on numerical methods (see for instance [START_REF] Altmann | A decoupling and linearizing discretization for weakly coupled poroelasticity with nonlinear permeability[END_REF]) in which naive coupling strategies are employed, in opposition to the celebrated fixed stress and undrained split [START_REF] Kim | Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits[END_REF][START_REF] Jakub | Robust fixed stress splitting for biot's equations in heterogeneous media[END_REF] approaches which allow to overpass the weak coupling regime for the simulation of (possibly single phase) poromechanics.

Finally, Assumption (H8) comprises typical practical scenarios. External body forces f ext with weaker regularity, e.g., associated to non-homogeneous traction boundary conditions, may be of interest for practical applications and thus also the analysis. For simpler presentation these are, however, not further discussed here.

The following theorem is the main result of our paper.

Theorem 1.3. Under Assumptions (H4)-(H7), there exists (at least) a global in time weak solution to the problem (1)-(8) in the sense of Definition 1.2.

Our proof for this theorem, to be detailed in what follows, relies on compactness arguments. We apply two successive regularizations of the problem. First, we soften the constraint (1d) into χ = G ǫ (φ), ǫ > 0 where G ǫ is a suitable regularization of the maximal monotone graph χ. A prototypical choice for G ǫ is φ → ǫ log φ-φ ♭ φ ♯ -φ for ǫ > 0. The hard constraint (1d) is recovered when ǫ tends to 0. We also regularize the mobilities in (1a) to remove the degeneracy in pure phase zone {s α = 0}: we replace s α by k ǫ (s α ) = max(ǫ, s α ) in (1a). This modification allows to derive some control on the phase pressure thanks to the control of the entropy production [START_REF] Jakub | Robust fixed stress splitting for biot's equations in heterogeneous media[END_REF]. In order to initiate the process, we establish the well-posedness of the elliptic problem consisting in one step of the backward Euler scheme with time step h > 0 for the regularized problem with ǫ > 0. Yet another regularization is required to justify properly our calculations: we make use of a Faedo-Galerkin (spectral) method to rigorously establish the regularity of the solutions to the discrete problem. Then we recover a global in time weak solution by passing first to the limit ǫ → 0, then h → 0.

The proof strongly builds on the second principle of thermodynamics, in the sense that the main estimate used in our existence proof is the control of the production of the Helmholtz free energy sketched out in Section 1.3, opening the way to possible extensions to more complex (but still thermodynamically consistent) physical settings.

2. The discrete and regularized system

2.1. Regularization. Let (G ǫ ) ǫ>0 ⊂ L 1 (φ ♭ , φ ♯ ) be a family of smooth increasing and onto functions from (φ ♭ , φ ♯ ) to R vanishing at φ ♭ +φ ♯ 2 with (26) G ′ ǫ (y) ≥ ǫ, y ∈ R. and (27) G ǫ -→ ǫ→0 0 in L 1 (φ ♭ , φ ♯ ).
Then one infers deduces from Dini's theorem that G ǫ tends to 0 uniformly on any compact of (φ ♭ , φ ♯ ), and in particular that

(28) G ǫ (y) -→ ǫ→0 0 for all y ∈ (φ ♭ , φ ♯ ).
We also set

(29) k ǫ (s) = max(ǫ, s) for all s ∈ R.
The regularization of the mobility k ǫ makes the problem coercive but yields a difficulty that was originally hidden by the degeneracy in pure phase regions. As suggested by the extension ( 21) of (γ ′ ) -1 outside of [γ ′ (0), γ ′ (1)], and in close connection to what was proposed in [START_REF] Cancès | An existence result for multidimensional immiscible two-phase flows with discontinuous capillary pressure field[END_REF][START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF], it becomes necessary to extend the capillary pressure function γ ′ : [0, 1] → R + into a maximum monotone graph. This amounts to define the monotone and anti-monotone graphs p n , p w by [START_REF] Carman | Fluid flow through granular beds[END_REF]

p n (s) = (-∞, γ(0) + γ ′ (0)] if s = 0, γ(s) + (1 -s)γ ′ (s) if s ∈ (0, 1], p w (s) = γ(s) -sγ ′ (s) if s ∈ [0, 1), (-∞, γ(1) -γ ′ (1)] if s = 1.
Then p np w = S -1 in the sense of the maximal monotone graphs. We also regularize and then extend to the whole R the function γ by defining

γ ǫ (s) = γ(0) + s 0 γ ′ ǫ (z)dz if s ∈ [0, 1], +∞ otherwise,
where

γ ′ ǫ (s) = γ ′ (0) + s 0 γ ′′ ǫ (z)dz and γ ′′ ǫ (s) = min ǫ -1 , max(ǫ, γ ′′ (s) , s ∈ [0, 1].
We infer from the dominated convergence theorem that γ ǫ converges uniformly towards γ on [0, 1] as ǫ tends to 0, and from Dini's theorem that the Lipschitz continuous function S ǫ defined by

(31) S ǫ (p) =      0 if p ≤ γ ′ ǫ (0), (γ ′ ǫ ) -1 (p) if γ ′ ǫ (0) ≤ p ≤ γ ′ ǫ (1), 1 if p ≥ γ ′ ǫ (1)
, converges uniformly towards S. We also incorporate these regularizations into the graphs p n,ǫ and p w,ǫ which are defined by [START_REF] Carman | Fluid flow through granular beds[END_REF] where γ and γ ′ have been replaced by γ ǫ and γ ′ ǫ respectively.

Lemma 2.1. Given p = (p n , p w ) ∈ R 2 and π ∈ R, then there exists a unique φ

= (φ n , φ w ) ∈ R 2 + with φ ♭ ≤ φ = φ n + φ w ≤ φ ♯ such that (32) p α = pα + π + G ǫ (φ) for some pα ∈ p α,ǫ (s n ) with s n = φ n /φ.
Moreover, the mapping Φ ǫ : pπ → φ is Lipschitz continuous with Lipschitz constant possibly blowing up with ǫ -1 .

Proof. It follows from [START_REF] Chen | Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution[END_REF] and from the definition (30) of the graphs p α,ǫ that

s n = S ǫ (p n -p w ) = S ǫ (p n -π -(p w -π)), s w = 1 -s n .
Moreover, the relation

G ǫ (φ) = s n p n + s w p w -γ ǫ (s n ) -π = s n (p n -π) + s w (p w -π) -γ ǫ (s n )
uniquely determines φ since G ǫ is an invertible function. Then we can reconstruct φ α = φs α . The Lipschitz continuity of Φ ǫ follows from the Lipschitz continuity of S ǫ , γ ǫ and G -1 ǫ .

The function Φ ǫ can be interpreted as the differential of a convex function. Before stating our next lemma, we introduce the convex function

G ǫ : [φ ♭ , φ ♯ ] → R + z → z φ ♭ +φ ♯ 2 G ǫ (a)da,
as well as the convex and compact subset of R 2

K φ = {φ = (φ n , φ w ) ∈ R 2 + | φ ♭ ≤ φ n + φ w ≤ φ ♯ }. Lemma 2.2. Define the convex function F ǫ : R 2 → R + by F ǫ (φ) = φγ ǫ (s n ) + G ǫ (φ) if φ ∈ K φ , +∞ otherwise, where φ = φ n + φ w and s n = φ n /φ, then F ǫ is convex. Moreover, for (p, π) ∈ R 2 × R and φ = Φ ǫ (p -π), then (33) (p n -π, p w -π) = DF ǫ (φ).
Proof. The function F ǫ is continuous on K φ and continuously differentiable on the interior

• K φ of K φ . Let φ ∈ • K φ , then DF ǫ (φ) = p n,ǫ (s n ) + G ǫ (φ), p w,ǫ (s n ) + G ǫ (φ) ,
the sets p α,ǫ (s n ) being assimilated to their single value since 0 < s n < 1 as φ ∈

• K φ . In view of formula [START_REF] Chen | Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution[END_REF], we deduce that (33) holds true. Since γ ǫ is increasing, since G ǫ (z) ≥ 0 and since φ ≥ φ ♭ for φ ∈ K φ , one gets the uniform lower bound

F ǫ (φ) ≥ φ ♭ γ ǫ (0) = φ ♭ γ(0) ≥ 0.
Finally, the function φ → φγ ǫ (s n ) is 1-homogeneous and (not strictly) convex on

• K φ , while φ → G ǫ (φ)
is the composition of the convex function G ǫ with the linear one φ → φ, so it is convex too, as well as F ǫ .

Remark 2.3. In the line of Lemma 2.1, one can show that F ǫ is uniformly convex for ǫ > 0. We leave to the reader the proof of this property which will not be used explicitly in what follows.

2.2. Faedo-Galerkin space discretization. Concerning the discretization w.r.t. space, we build on a Faedo-Galerkin approach. Since V is separable, there exists a family (v k ) k≥0 ⊂ V such that, denoting by

V k = span{v ℓ , ℓ ≤ k}, then k≥1 V k V = V.
Besides, we also introduce the complete orthonormal family (w k ) k≥0 of L 2 (Ω) made of eigenvectors of the Laplace equation with Neumann boundary conditions, i.e. ( 34)

Ω ∇w k • ∇v = λ k Ω w k v, ∀v ∈ H 1 (Ω),
with λ 0 = 0, w 0 ≡ √ m Ω (here and in what follows, m Ω stands for the d-dimensional Lebesgue measure of Ω), whereas λ k+1 ≥ λ k > 0 and w k L 2 (Ω) = 1 for k ≥ 1 as well. We denote by W k = span{w ℓ , ℓ ≤ k} ⊂ H 1 (Ω), and classical results from the spectral theory of self-adjoint compact operators (see for instance [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]) show that [START_REF] Droniou | The Gradient Discretisation Method[END_REF]

k≥1 W k L 2 (Ω) = L 2 (Ω). The (topological) dual W ′ k of W k is identified to W k
thanks to the Riesz theorem building on the L 2 (Ω)-scalar product. Eventually, we denote by u k ∈ V d the unique solution to ( 36)

Ω [2µ ε(u k ) + λ∇ • u k ] : ε(v) = b Ω w k ∇ • v, ∀v ∈ V d , k ≥ 0.
where w k is the k th eigenvector of the Laplace operator introduced above in (34). We denote by

U k = span {u ℓ , 1 ≤ ℓ ≤ k} , then U k ⊂ U thanks to Assumption (H6).
2.3. The discrete and regularized problem with frozen mobility and linearized gravity.

Since K φ is convex in R 2 , the orthogonal projection Π : R 2 → K φ φ = (φ n , φ w ) → Π(φ) = (Π n (φ), Π w (φ))
is uniquely defined. For φ = ( φn , φw ) ∈ L 2 (Ω) 2 , we define (37) φ = Π n ( φ) + Π w ( φ) and sα = Π α ( φ)/ φ.

Note that these definitions are consistent with the previous ones in the case where φ ∈ K φ .

Proposition 2.4. 2 , let ũ ∈ U , and given ǫ > 0 and h > 0, then for any k ≥ 1, there exists a unique

Let φ ⋆ = (φ ⋆ n , φ ⋆ w ) ∈ L ∞ (Ω) 2 with φ ⋆ ∈ K φ a.e. in Ω, let φ = ( φn , φw ) ∈ L 2 (Ω)
(φ k , p k , u k , θ k , π k ) such that p k = (p n,k , p w,k ) belongs to p D + V k , such that π k and θ k belong to W k with π k = M θ k , such that u k = û + u o k with û ∈ U and u o k ∈ U k respectively solving Ω σ : ε(v) = Ω f ( φ) • v, ∀v ∈ V d , with σ = 2µ ε( û) + λ∇ • ûI, (38a) Ω σ o k : ε(v) = Ω b π k ∇ • v, ∀v ∈ U k , with σ o k = 2µ ε(u o k ) + λ∇ • u o k I, (38b) such that φ k = (φ n,k , φ w,k ) = Φ ǫ (p k -π k ) satisfies (38c) Ω (φ k -b∇ • u k -θ k )w = Ω φ r w, ∀w ∈ W k ,
and such that (38d)

Ω φ α,k -φ ⋆ α h v + Ω k ǫ (s α ) µ α K( φ)∇ (p α,k -ρ α g • (x + ũ)) • ∇v = 0 for all v ∈ V k , α ∈ {n, w}.
Remark 2.5. As a result of the particular choice [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs-application to a turbulence model[END_REF] for the basis functions (u ℓ ) 1≤ℓ≤k of U k , the relation (38b) holds true for all v ∈ V d and not only for v ∈ U k . Indeed, since π k belongs to W k , it can be written decomposed into π k = k ℓ=0 π k,ℓ w ℓ with w ℓ fulfilling [START_REF] Coussy | Poromechanics[END_REF]. Then we deduce from (36)

that u o k = k ℓ=0 π k,ℓ u ℓ satisfies (39) Ω σ o k : ε(v) = Ω b π k ∇ • v, ∀v ∈ V d .
In other words, u o k is the genuine continuous solution to the linear mechanics system corresponding to the approximate right-hand side b∇π k , and thus so does

u k = û + u o k for the full body force term b∇π k + f ( φ), i.e. ( 40
) Ω σ k : ε(v) = Ω f ( φ) • v + bπ k ∇ • v , ∀v ∈ V d , with σ k = 2µ ε(u k ) + λ ∇ • u k I.
Then owing to Assumption (H6) and (H8), for all k ≥ 0 the triangle inequality yields

(41) ∇(∇ • u k ) L 2 (Ω) d ≤ C 1 λ f ( φ) L 2 (Ω) d + b C 1 λ π k H 1 (Ω) .
Proof of Proposition 2.4. Define

H k : (V k ) 2 × U k × W k × W k → (V ′ k ) 2 × U ′ k × W ′ k × W ′ k (p o k = (p o n,k , p o w,k ), u o k , θ k , π k ) → ((r n,k , r w,k ) , r u,k , r θ,k , r π,k ) by setting φ k = Φ ǫ (p o k + p D -π k ), p α,k = p o α,k + p D α , u k = û + u o k and r α,k , v V ′ k ,V k = Ω (φ α,k -φ ⋆ α )v + h Ω k ǫ (s α ) µ α K( φ)∇(p α,k -ρ α g • (x + ũ)) • ∇v, ∀v ∈ V k , r u,k , v U ′ k ,U k = Ω {2µ ε(u o k ) : ε(v) + λ(∇ • u o k )(∇ • v) -b π k ∇ • v} , ∀v ∈ U k , and
r θ,k = 2(M θ k -π k ), r π,k = -φ k + φ r + θ k + b∇ • u o k + b∇ • û + θ k - π k M . Let Y (i) k = p o,(i) k = (p o,(i) n,k , p o,(i) w,k ), u o,(i) k , θ (i)
k , π

(i) k , i = 1, 2, be two elements of (V k ) 2 × U k × W k ×
W k , and denote by R

(i) k = (r (i) n,k , r (i) w,k ), r (i)
u,k , r

(i) θ,k , r (i) π,k = H Y (i) k
, then one checks that

R (1) k -R (2)
k , Y

-Y (2) k = α∈{n,w} Ω (φ (1) k
α,k -φ (2) α,k )(p (1)
α,kπ

k -p (1)
α,k + π (2) k ) + h α∈{n,w} Ω k ǫ (s α ) µ α K( φ)∇(p (2)
α,k -p (2) α,k ) • ∇(p (1)
α,k -p (1)
α,k ) + Ω 2µ ε(u (2)
k -u (1)
k ) : ε(u (2)
k -u (1)
k ) + λ ∇ • (u (2)
k -u (1)
2

+ Ω M (θ (1)
k -θ (2)
k ) 2 + 1 M M (θ (1)
k -θ (2)
k )π

k + π (1)
.

with gravity contributions canceling, as well as those related to û. Owing to the convexity of Φ ǫ established in Lemma 2.2, one has α∈{n,w} Ω (φ

(1) α,k -φ (2) α,k )(p (1)
α,kπ

(1)

k -p (2)
α,k + π (2) k ) ≥ 0. Since K( φ) ≥ K ♭ I and µ α ≤ µ ♯ , one gets that α∈{n,w} Ω k ǫ (s α ) µ α K( φ)∇(p (1)
α,k -p (2) α,k ) • ∇(p (1)
α,k -p (2) α,k ) ≥ ǫ K ♭ µ ♯ p (1) k -p (2) k 2 V 2 .
Note that V k is equipped with the norm v V = ∇v (L 2 ) d which is a norm since we assumed that Γ D has positive measure thanks to Poincaré inequality. We also deduce from Korn inequality that there exists C 2 depending only on Ω and Γ D (but not on k) such that

Ω 2µ ε(u (1)
k -u (2) k ) : ε(u (1)
k -u (2)
k ) + λ ∇ • (u (1)
k -u (2) k ) 2 ≥ C 2 µ u (1) k -u (2) k 2 V d .
Further, elementary calculations show that

x 2 + (x -y) 2 ≥ 3 - √ 5 2 (x 2 + y 2 ), ∀x, y ∈ R, so that Ω M (θ (1) k -θ (2) k ) 2 + 1 M M (θ (1) k -θ (2) k ) -π (1) k + π (2) k 2 ≥ 3 - √ 5 2 M θ (1) k -θ (2) k L 2 (Ω) + 1 M π (1) k -π (2) k L 2 (Ω) .
Then we infer from previous estimates that there exists C 3 > 0 depending on the data of the continuous problem as well as on ǫ and h (but neither on k nor on ǫ) such that

R (1) k -R (2)
k , Y

-Y (2) k ≥ C 3 p (1) k -p (2) k 2 V 2 + u (1) k -u (2) k 2 V d + θ (1) k -θ (2) k 2 L 2 (Ω) + π (1) k -π (2) k 2 L 2 (Ω) (1) k
. As a consequence, H k is strongly monotone, whence there exists a unique Corollaire 17] (see also [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF][START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]). This also allows to reconstruct

Y k = (p o k , u o k , θ k , π k ) such that H k (Y k ) = 0 owing to [20,
φ k = Φ ǫ (p k -π k ).
The above proof (and thus the statement of Proposition 2.4) is still valid at the limit k → +∞, leading to the well-posedness of the limiting problem in the Hilbert space

V 2 × V d × L 2 (Ω) × L 2 (Ω).
The regularity provided by Proposition 2.4 is however not sufficient to carry out our mathematical study and to pass to the limit to recover the continuous problem. This was the motivation for the space discretization, the regularizing effect of which being needed to establish rigorously the following proposition.

Proposition 2.6. Define the approximate Helmholtz free energy F ǫ by

F ǫ (X) = Ω F ǫ (φ) + µ ε(u) : ε(u) + λ 2 |∇ • u| 2 + M 2 |θ| 2 for X = (φ, u, θ, π). Let φ ⋆ ∈ L ∞ (Ω; R 2 + ) with φ ⋆ ∈ K φ a.e.
in Ω, and let u ⋆ ∈ U and

π ⋆ = M θ ⋆ in L 2 (Ω) be such that (42) φ ⋆ -b∇ • u ⋆ -θ ⋆ = φ r , with φ ⋆ = φ ⋆ n + φ ⋆ w .
Denoting by X ⋆ = (φ ⋆ , u ⋆ , θ ⋆ , π ⋆ ), then for any φ ∈ L 2 (Ω) 2 and ũ ∈ H 1 (Ω) d , the unique solution

X k = (φ k , u k , θ k , π k ) to the discrete regularized problem, cf. Proposition 2.4, satisfies (43) 1 -C 5 hǫ 2 F ǫ (X k ) + hǫC 4   α∈{n,w} ∇p α,k 2
L 2 (Ω) d + ∇(G ǫ (φ k ) + π k ) 2 L 2 (Ω) d + ǫ π k 2 H 1 (Ω)   ≤ F ǫ (X ⋆ ) + Ω f ( φ) • (u k -u ⋆ ) + α∈{n,w} Ω (φ α,k -φ ⋆ α )p D α + C 6 h(1 + ũ 2 V d + φ 2 L 2 (Ω) 2 ),
for positive constants C 4 , C 5 and C 6 which neither depend on k, ǫ and h, nor on φ or ũ. Moreover, for ǫ small enough to ensure that C 5 hǫ 2 ≤ 1, there exists C 7 depending neither on k nor on φ (but possibly on ǫ and h) such that

(44) π k H 1 (Ω) + G ǫ (φ k ) H 1 (Ω) ≤ C 7 .
Proof. Testing (38d) by h (p α,kp D α ) ∈ V k and summing over α ∈ {n, w} provides

A k + B k + D k = R 1,k + R 2,k + R 3 ,
where

A k = α∈{n,w} Ω (φ α,k -φ ⋆ α )(p α,k -π a,k ), B k = α∈{n,w} Ω (φ α,k -φ ⋆ α )π α,k , D k = h α∈{n,w} Ω k ǫ (s α ) µ α K( φ)∇p α,k • ∇p α,k , R 1,k = α∈{n,w} Ω (φ α,k -φ ⋆ α )p D α , R 2,k = h α∈{n,w} Ω k ǫ (s α ) µ α K( φ)∇p α,k • ∇ ρ α g • (x + ũ) + p D α , R 3,k = -h α∈{n,w} Ω k ǫ (s α ) µ α K( φ)∇p D α • ∇ (ρ α g • (x + ũ)) .
Young's inequality implies that

R 2,k ≤ 1 2 D k + 1 2 h α∈{n,w} Ω k ǫ (s α ) µ α K( φ) 1/2 ∇ ρ α g • (x + ũ) + p D α 2
, so that (45)

A k + B k + 1 2 D k ≤ R 1,k + C 6 h(1 + ũ 2 V d )
for some C 6 independent on k,h, and ǫ. As a consequence of Lemma 2.2, we deduce from a convexity inequality that

A k = α∈{n,w} Ω (φ α,k -φ ⋆ α )(p α,k -π k ) ≥ Ω (F ǫ (φ k ) -F ǫ (φ ⋆ )).
On the other hand, we infer from (38c) and ( 42) that

B k = Ω (φ k -φ ⋆ )π k = Ω (b∇ • (u k -u ⋆ ) + θ k -θ ⋆ )π k . Since π k = M θ k , the elementary convexity inequality (x -y)x ≥ 1 2 (x 2 -y 2 ) provides Ω (θ k -θ ⋆ )π k ≥ M 2 θ k 2 L 2 (Ω) -θ * 2 L 2 (Ω) .
Besides, employing the same inequality again, we deduce from [START_REF] Kozeny | Uber kapillare Leitung der Wasser in Boden[END_REF] that

Ω b∇ • (u k -u ⋆ )π k = Ω σ k : ε(u k -u ⋆ ) - Ω f ( φ) • (u k -u ⋆ ) ≥ Ω µ (ε(u k ) : ε(u k ) -ε(u ⋆ ) : ε(u ⋆ )) + λ 2 |∇ • u k | 2 -|∇ • u ⋆ | 2 - Ω f ( φ) • (u k -u ⋆ ).
Collecting the above estimates, we get that (46)

A k + B k ≥ F ǫ (X k ) -F ǫ (X ⋆ ) - Ω f ( φ) • (u k -u ⋆ ).
Due to the assumptions on φ ⋆ and to the regularity of u ⋆ and θ ⋆ , the regularized Helmholtz free energy F ǫ (X ⋆ ) of X ⋆ = (φ ⋆ , u ⋆ , θ ⋆ , M θ ⋆ ) is finite. In the end, the gravity-related contribution is linear in u k which will allow its control.

On the other hand, since k ǫ (s α ) ≥ ǫ, since K( φ) ≥ K ♭ I and since 0 ≤ s α,k ≤ 1, one gets that

D k ≥hǫ K ♭ µ ♯ α∈{n,w} Ω |∇p α,k | 2 ≥hǫ K ♭ 2µ ♯ α∈{n,w} Ω |∇p α,k | 2 + Ω s α,k |∇p α,k | 2 =: D 1,k + D 2,k . (47)
The first term in the right-hand side provides some control on p α,k V . Concerning the second contribution, the relation α∈{n,w} s α,k ∇p α,k = 0, cf. [START_REF] Brenner | Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure[END_REF], allows to reformulate α∈{n,w}

s α,k |∇p α,k | 2 = α∈{n,w} s α,k |∇(p α,k + G ǫ (φ k ) + π k )| 2 = α∈{n,w} s α,k |∇p α,k | 2 + |∇(G ǫ (φ k ) + π k )| 2 ≥ |∇(G ǫ (φ k ) + π k )| 2 . ( 48
)
Introducing the non-decreasing function G ǫ : z → G ǫ (z)ǫ z, the latter term rewrites

|∇(G ǫ (φ k ) + π k )| 2 = |∇( G ǫ (φ k ) + ǫφ k + π k )| 2 = |∇( G ǫ (φ k ) + π k )| 2 + ǫ 2 |∇φ k | 2 + 2ǫ∇φ k • ∇π k + 2∇ G ǫ (φ k ) • ∇φ k . We deduce from the monotonicity of G ǫ that ∇ G ǫ (φ k ) • ∇φ k ≥ 0, hence (49) |∇(G ǫ (φ k ) + π k )| 2 ≥ 2ǫ∇φ k • ∇π k .
The particular choice for the space W k will be used here. Let w ℓ be such that (34) holds true, then since φ k ∈ H 1 (Ω), following from the Lipschitz continuity of Φ ǫ , it holds

Ω ∇φ k • ∇w ℓ = λ ℓ Ω φ k w ℓ , 0 ≤ ℓ ≤ k.
Thus, exploiting the fact that φ r ∈ H 1 (Ω) thanks to (H1) and ∇ • u k ∈ H 1 (Ω) thanks to (H6) and more precisely to (41), we can infer that for 0 ≤ ℓ ≤ k,

Ω ∇φ k • ∇w ℓ (34) = λ ℓ Ω φ k w ℓ (38c) = λ ℓ Ω (φ r + b ∇ • u k + θ k )w ℓ (34) = Ω ∇(φ r + b∇ • u k + θ k ) • ∇w ℓ . Therefore, Ω ∇φ k • ∇w = Ω ∇(φ r + b∇ • u k + θ k ) • ∇w, ∀w ∈ W k .
In particular for w = π k , and bearing in mind that M θ k = π k , this yields

Ω ∇φ k • ∇π k = Ω 1 M |∇π k | 2 + ∇(∇ • u k ) • b∇π k + ∇φ r • ∇π k .
It follows from ( 41) combined with Cauchy-Schwarz that

Ω ∇(∇ • u k ) • ∇π k ≥ - b C 1 λ π k 2 H 1 (Ω) - C 1 λ f ( φ) L 2 (Ω) d π k H 1 (Ω) ,
so that Young's inequality provide

Ω ∇(∇ • u k ) • ∇π k ≥ - b C 1 λ π k 2 H 1 (Ω) - 1 4b 1 M - b 2 C 1 λ π k 2 H 1 (Ω) - 1 M - b 2 C 1 λ -1 b f ( φ) 2 L 2 (Ω) d .
Similarly, we obtain that

Ω ∇φ r • ∇π k ≥ - 1 4 1 M - b 2 C 1 λ π k 2 H 1 (Ω) - 1 M - b 2 C 1 λ -1 ∇φ r 2 L 2 (Ω) d .
Therefore, using again π k = M θ k , we obtain that

Ω ∇φ k • ∇π k ≥ 1 2 1 M - b 2 C 1 λ π k 2 H 1 (Ω) -M θ k 2 L 2 (Ω) (50)
- 1 M - b 2 C 1 λ -1 ∇φ r 2 L 2 (Ω) d + b 2 f ( φ) 2 L 2 (Ω) d ,
the constant between the parentheses being strictly positive owing to Assumption (H7). Combining the above estimate with [START_REF] Terzaghi | Erdbaumechanik auf bodenphysikalischer Grundlage[END_REF], and since

F ǫ (X k ) ≥ M 2 θ k 2 L 2 (Ω) , one gets that (51) ∇(G ǫ (φ k ) + π k ) 2 L 2 (Ω) d ≥ C 8 ǫ π k 2 H 1 (Ω) -4ǫF ǫ (X k )-C -1 8 ǫ ∇φ r 2 L 2 (Ω) d + b 2 f ( φ) 2 L 2 (Ω) d . for C 8 = 1 M - b 2 C 1 λ > 0. As a consequence, (48) yields (52)
Ω α∈{n,w}

s α,k |∇p α,k | 2 ≥ C 8 ǫ π k 2 H 1 (Ω) -4ǫF ǫ (X k )-C -1 8 ǫ ∇φ r 2 L 2 (Ω) d + b 2 f ( φ) 2 L 2 (Ω) d ,
and the last term enters an updated definition of C 6 for a uniformly bounded ǫ. We collect ( 46), ( 47) and ( 52) in [START_REF] Schulz | Beyond Kozeny-Carman: predicting the permeability in porous media[END_REF] to recover [START_REF] Mielke | A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems[END_REF]. For ǫ small enough so that C 5 hǫ 2 ≤ 1, one deduces from ( 43) that ( 53)

∇(G ǫ (φ k ) + π k ) 2 L 2 (Ω) d + ǫ π k 2 H 1 (Ω) ≤ C ǫ 1 + 1 h
for some C depending neither on k, φ, ǫ nor on h (C will vary along the following lines, but will remain independent on the aforementioned parameters). Using the elementary inequality a+b 2 ≤ 2 a 2 + 2 b 2 , one further gets that

(54) ∇G ǫ (φ k ) 2 L 2 (Ω) d ≤ 2 ∇ (G ǫ (φ k ) + π k ) 2 L 2 (Ω) d + 2 ∇π k 2 L 2 (Ω) d ≤ C ǫ 1 + 1 h 1 + 1 ǫ .
Moreover, the relation G ǫ (φ k ) = s n,k p n,k + s w,k p w,kγ ǫ (s n,k )π k holds everywhere in Ω, and in particular also on Γ D , where

|G ǫ (φ k )| ≤ |p D n | + |p D w | + γ ǫ (1) + |π k | ≤ |p D n | + |p D w | + γ(1) + ǫ + |π k |.
We infer from a trace theorem that

π k L 2 (Γ D ) ≤ C π k H 1 (Ω) .
Therefore, we get that

G ǫ (φ k ) L 2 (Γ D ) ≤ C(1 + π k H 1 (Ω)
), and thus that 53) and ( 54) in ( 55) gives [START_REF] Peletier | Variational modelling: Energies, gradient flows, and large deviations[END_REF], concluding the proof of Proposition 2.6.

(55) G ǫ (φ k ) H 1 (Ω) ≤ C(1 + π k H 1 (Ω) + ∇G ǫ (φ k ) L 2 (Ω) d ) since y → ∇y L 2 (Ω) d + y L 2 (Γ D ) is equivalent to the usual H 1 (Ω) norm. Incorporating (
2.4. Passing to the limit k → +∞. This section is devoted to the proof of the following proposition, which is deduced from Propositions 2.4 and 2.6 after letting k tend to +∞.

Proposition 2.7. Let φ ⋆ ∈ L ∞ (Ω; R 2 + ) with φ ⋆ ∈ K φ a.e.
in Ω, and let u ⋆ ∈ V d and π ⋆ = M θ ⋆ in L 2 (Ω) be such that (42) holds. Assume that C 5 hǫ 2 ≤ 1, then for any φ ∈ L 2 (Ω) 2 and ũ ∈ U , there exists a unique solution (φ, p, u, θ, π) with M θ = π and pp D ∈ V 2 to the following problem:

(56a) Ω σ : ε(v) = Ω b π∇ • v + Ω f ( φ) • v, ∀v ∈ V d , with σ = 2µ ε(u) + λ ∇ • u I, such that φ = (φ n , φ w ) = Φ ǫ (p -π) satisfies (56b) φ -b ∇ • u -θ = φ r .
and such that

(56c) Ω φ α -φ ⋆ α h v + Ω k ǫ (s α ) µ α K( φ)∇ (p α -ρ α g • (x + ũ)) • ∇v = 0 for all v ∈ V, α ∈ {n, w}.
Moreover, it satisfies

(57) 1 -C 5 hǫ 2 F ǫ (X) + hǫC 4   α∈{n,w} ∇p α 2 L 2 (Ω) d + ∇(G ǫ (φ) + π) 2 L 2 (Ω) d + ǫ π 2 H 1 (Ω)   ≤ F ǫ (X ⋆ ) + Ω f ( φ) • (u -u ⋆ ) + α∈{n,w} Ω (φ α -φ ⋆ α )p D α + C 6 h 1 + ũ 2 V d + φ 2 L 2 (Ω)
and

(58) π H 1 (Ω) + G ǫ (φ) H 1 (Ω) ≤ C 7 .
Proof. The a priori estimates derived in Proposition 2.6 can be refined. Using Korn's inequality and ( 59)

Ω f ( φ) • (u -u ⋆ ) ≤ 1 2 1 -C 5 hǫ 2 F ǫ (X k ) + C u ⋆ 2 L 2 (Ω) d + f ( φ) 2 L 2 (Ω)
for some constant C > 0 depending Korn's constant, Lamé parameters, as well as 1 -C 5 hǫ 2 -1 , thus being uniformly bounded. Starting from ( 43), the state-dependent contribution can be compensated in the Helmholtz energy F ǫ (X k ) on the left hand side of [START_REF] Mielke | A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems[END_REF]. The remaining terms on the right hand side of (59) are uniformly bounded in k, given φ ∈ L 2 (Ω) with M θ = π. The control of the energy F ǫ (X) provides a uniform control on the H 1 (Ω) d norm of u k thanks to Korn's inequality, and even a uniform control on the H 1 (Ω) d norm of ∇ • u k thanks to Assumption (H6). Therefore, there exists u ∈ U such that

u k -→ k→+∞ u weakly in H 1 (Ω) and ∇ • u k -→ k→+∞ ∇ • u weakly in H 1 (Ω).
Since Φ ǫ is continuous owing to Lemma 2.1, one infers from (60a) and (60b) that

φ k = Φ ǫ (p k -π k ) -→ k→+∞ Φ ǫ (p -π) = φ a.e. in Ω,
with φ ∈ K φ a.e. in Ω since φ k does. The aforementioned convergences are enough to pass to the limit in [START_REF] Kozeny | Uber kapillare Leitung der Wasser in Boden[END_REF], which gives (56a), and in (38d), leading to (56c). Thanks to [START_REF] Droniou | The Gradient Discretisation Method[END_REF], passing to the limit in (38c) provides that

Ω (φ -b∇ • u -θ -φ r )w = 0, ∀w ∈ L 2 (Ω),
which is equivalent to claiming to (56b) holds in L 2 (Ω) and thus almost everywhere and in H 1 (Ω). Inequality (57) is recovered from ( 43) by invoking the weak lower semi-continuity of the left-hand side and linearity of the right hand side.

2.5.

Unfreezing the mobility and recovering nonlinear gravitational energy. In this section, we establish the existence of a solution to the previous problem with the additional constraints that φ = φ and ũ = u.

Proposition 2.8. Let φ ⋆ ∈ L ∞ (Ω) 2 with φ ⋆ ∈ K φ a.e.
in Ω, and let u ⋆ ∈ U and π ⋆ = M θ ⋆ in L 2 (Ω) be such that (42) holds. Let h ∈ 0, 1 C 9 , with C 9 defined below. Then there exists a solution

(φ, p, u, θ, π) with p -p D ∈ V 2 , φ ∈ H 1 (Ω) 2 , u ∈ U , π ∈ H 1 (Ω), G ǫ (φ) ∈ H 1 (Ω)
and M θ = π to the following problem:

(61)

Ω σ : ε(v) = Ω b π∇ • v + Ω f (φ) • v, ∀v ∈ V d , with σ = 2µ ε(u) + λ ∇ • u I, such that φ = (φ n , φ w ) = Φ ǫ (p -π) satisfies (62) φ -b ∇ • u -θ = φ r .
and such that

(63) Ω φ α -φ ⋆ α h v + Ω k ǫ (s α ) µ α K(φ)∇ (p α -ρ α g • (x + u)) • ∇v = 0 for all v ∈ V, α ∈ {n, w}. Proof. Let T : L 2 (Ω) 2 × V d → L 2 (Ω) 2 × V d mapping ( φ, ũ) to (φ, u) as in Proposition 2.7.
Let the product space be equipped by the product norm ||| • ||| defined through

|||(φ, u)||| 2 := φ 2 L 2 (Ω) 2 + u 2 µ,λ (φ, u) ∈ L 2 (Ω) 2 × V d , where we have set u 2 µ,λ = Ω 2µ ε(u) : ε(u) + λ(∇ • u) 2 , ∀u ∈ V d .
The norm • µ,λ is equivalent to the H 1 (Ω) d norm on V d thanks to Poincaré's and Korn's inequalities.

First, since φ takes its values in the bounded set K φ of R 2 , and since Ω is bounded, then φ 2 L 2 (Ω) 2 ≤ R φ for some R φ not depending on φ; here, we recall the use of the orthogonal projection Π onto K φ to freeze the porosities, cf. [START_REF] Grisvard | Singularités en elasticité[END_REF]. Furthermore, from (57) we can infer that

u 2 µ,λ ≤ C 1 + h ũ 2 V d + h φ 2 L 2 (Ω) 2 ≤ C 9 1 + h ũ 2 µ,λ + hR φ .
for suitable constant C, C 9 > 0 independent of u and ũ, where the latter bound follows from Korn's inequality. Thus, assuming C 9 h < 1 and ũ 2 µ,λ ≤ R u := C 9 (1+hR φ ) 1-C 9 h , then

|||(φ, u)||| 2 ≤ R φ + R u .
In particular, T maps the ball of radius (R

φ + R u ) 1/2 of L 2 (Ω) 2 × V d into itself.
Second, let us show that T is compact. For this, let ( φk , ũk

) k≥1 ⊂ L 2 (Ω) 2 × V d be a bounded sequence in L 2 (Ω) 2 × H 1 (Ω) d . Then, up to subsequence, there exist ( φ⋆ , ũ⋆ ) ⊂ L 2 (Ω) 2 × V d such that (64) φk -→ k→∞ φ⋆ weakly in L 2 (Ω) 2 , ũk -→ k→∞ ũ⋆ weakly in H 1 (Ω) d .
The latter implies

(65) ũk -→ k→∞ ũ⋆ strongly in L 2 (Ω) d .
Consider the corresponding sequence (φ k , u k ) k≥1 defined by (φ k , u k ) = T ( φk , ũk ), accompanied with X k solving (56a)-(56c) and the uniform stability bound (57) with constants independent of k and with F (X k ) ≥ 0. From (57), we infer that, up to a subsequence,

u k -→ k→∞ u ⋆ weakly in H 1 (Ω) d and strongly in L 2 (Ω) d (66a) π k -→ k→∞ π ⋆ weakly in H 1 (Ω) and strongly in L 2 (Ω) (66b)
for some u ⋆ ∈ V d and π ⋆ ∈ H 1 (Ω). Due the linearity of (56a), cf. ( 8), we can consider the limit k → ∞, yielding (67)

Ω 2µ ε(u ⋆ ) : ε(v) + λ(∇ • u ⋆ )(∇ • v) = Ω b π ⋆ ∇ • v + Ω f ( φ⋆ ) • v, ∀v ∈ V d .
In addition, testing (56a) with v = u k yields ( 68)

Ω 2µ ε(u k ) : ε(u k ) + λ(∇ • u k ) 2 = Ω b π k ∇ • u k + Ω f ( φk ) • u k .
From the convergences (64) and (66), combined in suitable pairs of strong and weak convergence, we can infer for the right hand side of (68) that

Ω b π k ∇ • u k + Ω f ( φk ) • u k -→ k→∞ Ω b π ⋆ ∇ • u ⋆ + Ω f ( φ⋆ ) • u ⋆ ,
which implies convergence of the left hand side by combining (67) (tested with v = u ⋆ ) and ( 68)

u k 2 µ,λ = Ω 2µ ε(u k ) : ε(u k ) + λ(∇ • u k ) 2 -→ k→∞ Ω 2µ ε(u ⋆ ) : ε(u ⋆ ) + λ(∇ • u ⋆ ) 2 = u ⋆ 2 µ,λ .
Norm convergence in • µ,λ together with weak convergence, implies strong convergence

(69) u k -→ k→∞ u ⋆ strongly in H 1 (Ω) d .
Furthermore, we deduce from the fact that φ k = Φ ǫ (p kπ k ) and from the Lipschitz continuity of Φ ǫ , cf. Lemma 2.1, that

∇φ k L 2 (Ω) d×2 ≤ Λ ǫ ∇p k -∇π k L 2 (Ω) d×2 ≤ Λ ǫ ∇p k L 2 (Ω) d×2 + ∇π k L 2 (Ω) d
for some Λ ǫ depending on ǫ (and the uniform bound on ( φk , ũk ) k≥1 ). Then we infer from (57) that ∇φ k L 2 (Ω) d×2 ≤ C for some C, and we conclude (70) φ k -→ k→∞ φ ⋆ weakly in H 1 (Ω) d and strongly in L 2 (Ω).

Thus, finally, from (69) and (70), it follows that T is compact. Continuity of T follows along the proof of compactness. Indeed, the derived convergences are sufficient to infer respective convergence of the accompanying sequences (θ k ) k≥1 and (s α,k ) k≥1 . Therefore, we can pass to the limit in (56a), (56b) and (56c), showing that (φ ⋆ , u ⋆ ) is the solution to the problem described in Proposition 2.7 corresponding ( φ⋆ , ũ⋆ ). From the uniqueness result stated in Proposition 2.7, we deduce that (φ ⋆ , u ⋆ ) = T ( φ⋆ , ũ⋆ ), and the continuity of T follows.

The operator T then fulfills all the assumptions of the Schauder fixed point theorem, ensuring the existence of (at least) one fixed point for T . This concludes the proof of Proposition 2.8.

The semi-discrete in time system without regularization

Our goal in this section is to get rid of the two regularizations we incorporated in the system in Section 2.1, that are:

(i) the regularization (29) of the mobilities to make them non-degenerate;

(ii) the regularization of the capillary energy density function F ǫ introduced in Lemma 2.2 to make it uniformly convex. To this end, we first derive uniform estimates w.r.t. ǫ, to be used to let it tend to 0.

3.1. Uniform estimates w.r.t. ǫ and h. Let us start from a solution (φ ǫ , p ǫ , u ǫ , θ ǫ , π ǫ ) to the regularized semi-discrete system as constructed in Proposition 2.8. We stress here the dependence of the solution in the regularization parameter ǫ > 0.

The core result of this section is the following energy estimate, which involves the increasing continuous function ξ ǫ : [0, 1] → R + defined by

ξ ǫ (s) = s 0 a(1 -a)γ ′′ ǫ (a)da, ∀s ∈ [0, 1].
Proposition 3.1. Let ǫ ∈ (0, 1 4 ] and h ∈ 0, 1 C 9 , then the solution (φ ǫ , p ǫ , u ǫ , θ ǫ , π ǫ ) to the semi-discrete in time and regularized problem fulfills the uniform w.r.t. ǫ and h estimate

(71) (1 -C 10 h)F ǫ (X ǫ ) + K ♭ 4µ ♯ h 2 ∇ξ ǫ (s n,ǫ ) 2 L 2 (Ω) d + C 11 π ǫ 2 H 1 (Ω) + ∇(G ǫ (φ ǫ ) + π ǫ ) 2 L 2 (Ω) d + ǫ α∈{n,w} ∇p α,ǫ 2 L 2 (Ω) d ≤ F ǫ (X ⋆ ) + Ω f (φ ǫ ) • (u ǫ -u ⋆ ) + α∈{n,w} Ω (φ α,ǫ -φ ⋆ α ) p D α + C 6 h,
with C 10 > 0 and C 11 > 0 depending neither on ǫ nor on X ⋆ and nor on h.

Proof. The proof shares several features with the one of Proposition 2.6. Choosing v = h(p α,ǫp D α ) in (63) and summing over α ∈ {n, w}, then proceeding similarly as in the proof of Proposition 2.6 gives, cf. in particular ( 45), (72)

A ǫ + B ǫ + 1 2 D ǫ ≤ R ǫ + C 6 h(1 + u ǫ 2 V d ),
where we have set

A ǫ = α∈{n,w} Ω (φ α,ǫ -φ ⋆ α )(p α,ǫ -π ǫ ), B ǫ = Ω (φ ǫ -φ ⋆ )π ǫ , D ǫ = h α∈{n,w} Ω k ǫ (s α,ǫ ) µ α K(φ ǫ )∇p α,ǫ • ∇p α,ǫ and R ǫ = α∈{n,w} Ω (φ α,ǫ -φ ⋆ α ) p D α .
As in the proof of Propostion 2.6, cf. also [START_REF] Seguin | Multi-component multiphase flow through a poroelastic medium[END_REF], again utilizing the relation φ ǫ = Φ ǫ (p ǫπ ǫ ) and Lemma 2.2, we end up with (73)

A ǫ + B ǫ ≥ F ǫ (X ǫ ) -F ǫ (X ⋆ ) - Ω f (φ ǫ ) • (u ǫ -u ⋆ ).
Concerning the term D ǫ , we start by noticing that k ǫ (s α,ǫ ) ≥ (s α,ǫ +ǫ)/2. Then taking inspiration on what was done in Section 1.4, we underestimate D ǫ by

D ǫ ≥ D 1,ǫ + D 2,ǫ with D 1,ǫ = K ♭ 2µ ♯ hǫ Ω α∈{n,w} |∇p α,ǫ | 2 and D 2,ǫ = K ♭ 2µ ♯ h Ω α∈{n,w} s α,ǫ |∇p α,ǫ | 2 .
Using the elementary inequality (a+b) 2 ≥ a 2 /2-b 2 and the decomposition p α,ǫ = pα,ǫ +π ǫ +G ǫ (φ ǫ ) of the phase pressures, we get that

(74) D 1,ǫ ≥ K ♭ 4µ ♯ hǫ Ω α∈{n,w} |∇p α,ǫ | 2 - K ♭ µ ♯ hǫ Ω |∇(G ǫ (φ ǫ ) + π ǫ )| 2 .
Besides, from similar arguments to those presented in Section 1.4, we can rewrite (75)

Ω α∈{n,w}

s α,ǫ |∇p α,ǫ | 2 = Ω s n,ǫ s w,ǫ |∇(p n,ǫ -p w,ǫ )| 2 + Ω |∇(G ǫ (φ ǫ ) + π ǫ )| 2 .
So for ǫ ≤ 1/4, we get that

D ǫ ≥ K ♭ 2µ ♯ h Ω s n,ǫ s w,ǫ |∇(p n,ǫ -p w,ǫ )| 2 (76) + K ♭ 4µ ♯ hǫ Ω α∈{n,w} |∇p α,ǫ | 2 + K ♭ 4µ ♯ h Ω |∇(G ǫ (φ ǫ ) + π ǫ )| 2 .
For the first term in the above right-hand side, either p n,ǫ -p w,ǫ = γ ′ (s n,ǫ ) or s n,ǫ s w,ǫ = 0. Therefore,

Ω s n,ǫ s w,ǫ |∇(p n,ǫ -p w,ǫ )| 2 = Ω s n,ǫ (1 -s n,ǫ )|∇γ ′ ǫ (s n,ǫ )| 2 = Ω |∇ξ ǫ (s n,ǫ )| 2 . (77)
Concerning the last term in the right-hand side of (76), the approach adopted in the proof of Proposition 2.6 leads to a control depending on ǫ, so modifications are needed. We can use π ǫ = M θ ǫ , (62), and the fact that φ r ∈ H 1 (Ω) to rewrite

|∇(G ǫ (φ ǫ ) + π ǫ )| 2 = |∇π ǫ | 2 + |∇G ǫ (φ ǫ )| 2 + 2∇G ǫ (φ ǫ ) • ∇π ǫ = |∇π ǫ | 2 + |∇G ǫ (φ ǫ )| 2 + 2M ∇G ǫ (φ ǫ ) • ∇ (φ ǫ -b(∇ • u ǫ ) -φ r ) .
So Young's inequalityfor any ξ 1 > 0 and the monotonicty of G ǫ provide that

|∇(G ǫ (φ ǫ ) + π ǫ )| 2 ≥ |∇π ǫ | 2 -M 2 |∇(b∇ • u ǫ + φ r )| 2 (78) ≥ |∇π ǫ | 2 -(bM ) 2 (1 + ξ 1 ) |∇(∇ • u ǫ )| 2 -M 2 1 + 1 4ξ 1 |∇φ r | 2 ,
On the other hand, we infer from Assumption (H6), recalling [START_REF] Leverett | Capillary behavior in porous solids[END_REF], and applying Young's inequality for any ξ 2 > 0, that

Ω |∇(∇ • u ǫ )| 2 ≤ b 2 C 1 λ 2 (1 + ξ 2 ) π ǫ 2 H 1 (Ω) + C 1 λ 2 1 + 1 4ξ 2 f (φ ǫ ) 2 L 2 (Ω) d .
Then we deduce from (78) that

∇(G ǫ (φ ǫ ) + π ǫ ) 2 L 2 (Ω) d ≥ 1 - C 1 M b 2 λ 2 (1 + ξ 1 ) (1 + ξ 2 ) π ǫ 2 H 1 (Ω) -M 2 θ ǫ 2 L 2 (Ω) -M 2 1 + 1 4ξ 1 |∇φ r | 2 L 2 (Ω) -1 + 1 4ξ 2 f (φ ǫ ) 2 L 2 (Ω) d .
Under the weak coupling assumption (H7), one can choose ξ 1 > 0 and ξ 2 ∈ (0, 1) to satisfy

1 + ξ 1 = 1 1 -ξ 2
, and

ξ 2 = λ 2 -(M C 1 b 2 ) 2 λ 2 + 3(M C 1 b 2 ) 2 ∈ (0, 1),
and together with the uniform bounds on f and φ r following from (H1) and (H8), one obtains (79)

Ω |∇(G ǫ (φ ǫ ) + π ǫ )| 2 ≥ C 11 π ǫ 2 H 1 (Ω) -C 10 M 2 θ ǫ 2 L 2 (Ω) -C 12 .
with C 11 > 0, C 10 > 0, and C 12 independent on ǫ, X ⋆ and h. Combining (73)-( 77) and ( 79) in (72) and noticing again that M θ ǫ 2

L 2 (Ω) ≤ 2F ǫ (X ǫ ), in addition to u ǫ 2 V d ≤ CF ǫ (X ǫ
) due to Korn's inequality (for some C > 0 for simplicity entering the definition of C 10 ), to be used in (72), provides the desired estimate (71). Lemma 3.2. There exists C 13 > 0 such that, for all ǫ ∈ (0, 1 4 ] and all h ∈ (0, 1 C 10 ] (without loss of generality, we assume from now on C 9 ≤ C 10 ), there holds

G ǫ (φ ǫ ) 2 H 1 (Ω) ≤ C 13 1 + 1 h .
Proof. We will adapt here the program to derive [START_REF] Peletier | Variational modelling: Energies, gradient flows, and large deviations[END_REF]. As C 10 h ≤ 1, and since F ǫ (X ǫ ) ≥ 0 while F ǫ (X ⋆ ) is finite, we infer from (71) that

π ǫ 2 H 1 (Ω) + ∇(G ǫ (φ ǫ ) + π ǫ ) 2 L 2 (Ω) d ≤ C 1 + 1 h ,
where C denotes again a generic quantity independent of ǫ and h. We note that the state-dependent terms are either uniformly bounded, or (like the u ǫ -dependent term) can be compensated on the left hand side under (8) also affecting the value of C. Then updating (54) leads to

∇G ǫ (φ ǫ ) 2 L 2 (Ω) d ≤ 2 ∇ (G ǫ (φ ǫ ) + π ǫ ) 2 L 2 (Ω) d + 2 ∇π ǫ 2 L 2 (Ω) d ≤ C 1 + 1 h .
Then reproducing the arguments detailed for obtaining (55), one gets that

G ǫ (φ ǫ ) 2 H 1 (Ω) ≤ C(1 + π ǫ 2 H 1 (Ω) + ∇G ǫ (φ ǫ ) 2 L 2 (Ω) d ) ≤ C 1 + 1 h ,
concluding the proof of Lemma 3.2.

Corollary 3.3. There exist C 14 , C 15 > 0 such that, for all h ∈ (0, 1 2C 10 ] and ǫ ∈ (0, 1 4 ], there holds

F ǫ (X ǫ ) + C 14 h ξ ǫ (s n,ǫ ) 2 H 1 (Ω) + π ǫ 2 H 1 (Ω) + G ǫ (φ ǫ ) 2 H 1 (Ω) ≤ (1 + 2C 10 h) F ǫ (X ⋆ ) + Ω f (φ ǫ ) • (u ǫ -u ⋆ ) + α∈{n,w} Ω (φ α,ǫ -φ ⋆ α ) p D α + C 15 h.
Proof. Dividing (71) by 1-C 10 h, then the elementary inequalities 2 ≥ 1+2C 10 h ≥ (1-C 10 h) -1 ≥ 1 (recall here that we assumed 2C 10 h ≤ 1) provide

F ǫ (X ǫ ) + Ch π ǫ 2 H 1 (Ω) + ∇ξ ǫ (s n,ǫ ) 2 L 2 (Ω) d ≤ (1 + 2C 10 h)   F ǫ (X ⋆ ) + Ω f (φ ǫ ) • (u ǫ -u ⋆ ) + α∈{n,w} Ω (φ α,ǫ -φ ⋆ α ) p D α   + C 6 h. Since φ ǫ , φ ⋆ ∈ K φ and since p D ∈ L ∞ (Ω) 2 too, α∈{n,w} Ω (φ α,ǫ -φ ⋆ α ) p D α ≤ C.
Moreover, it follows from (H2) that the function ξ ǫ is uniformly bounded on [0, 1], hence we get

(80) F ǫ (X ǫ ) + Ch π ǫ 2 H 1 (Ω) + ξ ǫ (s n,ǫ ) 2 H 1 (Ω) ≤ (1 + 2C 10 h) F ǫ (X ⋆ ) + Ω f (φ ǫ ) • (u ǫ -u ⋆ ) α∈{n,w} Ω (φ α,ǫ -φ ⋆ α ) p D α + C ′ h
for some C ′ also uniform in ǫ and h. Combining (80) with Lemma 3.2 provides the desired result.

3.2. Passing to the limit ǫ → 0. Let us get rid of the regularization in the mobilities and in the capillary energy density by passing to the limit ǫ → 0 while keeping h > 0 fixed. The main result of the current section is the following.

Proposition 3.4. Let h ∈ (0, 1], then there exists X = (φ, u, θ, π) with φ ∈ K φ a.e. in Ω, φ ∈ H 1 (Ω) and ξ(s n ) -ξ(s D n ) ∈ V , with u ∈ U , and θ, π ∈ H 1 (Ω), as well as some χ ∈ H 1 (Ω) with χ ∈ χ(φ) such that (81) π = M θ and φ -b ∇ • u -θ = φ r , such that (82)
Ω φ n -φ ⋆ n h v - Ω 1 µ n K(φ) (∇ψ(s n ) + s n ∇ ((π + χ) -ρ n g • (x + u))) • ∇v = 0, ∀v ∈ V, (83)
Ω φ w -φ ⋆ w h v - Ω 1 µ w K(φ) (-∇ψ(s n ) + s w ∇ ((π + χ) -ρ w g • (x + u))) • ∇v = 0, ∀v ∈ V, and (84)
Ω (2µ ε(u) : ε(v) + λ(∇ • u)(∇ • v)) = Ω bπ ∇ • v + Ω f (φ) • v, ∀v ∈ V d .
Moreover, X and χ satisfy

(85) F (X) + C 14 h ξ(s n ) 2 H 1 (Ω) + π 2 H 1 (Ω) + χ 2 H 1 (Ω) ≤ (1 + 2C 10 h) F (X ⋆ ) + Ω f (φ) • (u -u ⋆ ) + α∈{n,w} Ω (φ α -φ ⋆ α ) p D α + C 15 h.
Proof. Let (ǫ k ) k≥0 ⊂ (0, 1 4 ] be a sequence of positive regularization parameters tending to 0 as k tend to +∞, then we investigate, thanks to compactness arguments, the behavior when k goes to ∞ of the sequence (X ǫ k ) k≥0 provided by Proposition 2.8. For the ease of reading, we will omit the subscript k, and we will denote by ǫ → 0 the limit k → +∞.

We deduce from Corollary 3.3 that there exists C not depending on ǫ such that

ξ ǫ (s n,ǫ ) H 1 (Ω) ≤ C, G ǫ (φ ǫ ) H 1 (Ω) ≤ C, π ǫ H 1 (Ω) ≤ C, θ ǫ H 1 (Ω) ≤ C, u ǫ U ≤ C,
while φ ǫ ∈ K φ and 0 ≤ s α,ǫ ≤ 1 holds true for all ǫ > 0. Here, we employ the same arguments as in the proof of Lemma 3.2 and Corollary 3.3, regarding state-dependent terms on the right hand side. Then we infer from (62) that

φ ǫ H 1 (Ω) ≤ C.
As a consequence, there exists Ξ, π, θ, φ and χ in H 1 (Ω) such that, up to a subsequence, (86)

ξ ǫ (s n,ǫ ) -→ ǫ→0 Ξ, π ǫ -→ ǫ→0 π, θ ǫ -→ ǫ→0 θ, φ ǫ -→ ǫ→0 φ and G ǫ (φ ǫ ) -→ ǫ→0 χ,
as well as u ∈ U such that, up to a subsequence, (87)

u ǫ -→ ǫ→0 u, ∇ • u ǫ -→ ǫ→0 ∇ • u,
the above convergences holding almost everywhere in Ω as well as in the weak H 1 (Ω) sense. One further readily deduces from the uniform convergence of ξ ǫ towards ξ that ξ ǫ (s D n ) converges weakly in H 1 (Ω) towards ξ(s D n ), and that Ξξ(s D n ) ∈ V . Moreover, since φ ǫ ∈ K φ and 0 ≤ s α,ǫ ≤ 1 for all ǫ > 0, there exists φ ∈ L ∞ (Ω) 2 with φ ∈ K φ a.e. in Ω and s α ∈ L ∞ (Ω) with 0 ≤ s n ≤ 1 and

s w = 1 -s n such that (88) φ ǫ -→ ǫ→0 φ and s α,ǫ -→ ǫ→0 s α
in the L ∞ (Ω) 2 -and L ∞ (Ω)-weak-⋆ senses respectively. Denote by

Ω ♭ = {x ∈ Ω | φ(x) = φ ♭ }, Ω ♯ = {x ∈ Ω | φ(x) = φ ♯ } and O = Ω \ (Ω ♭ ∪ Ω ♯ ),
and let x ∈ O such that φ ǫ (x) → φ(x) and G ǫ (φ ǫ (x)) → χ(x) as ǫ → 0. Then for all η > 0, for ǫ smaller than some ǫ ⋆ depending on η and on x, φ ǫ (x) ∈ (φ ♭ + η, φ ♯η). As a consequence,

G ǫ (φ ♭ + η) ≤ G ǫ (φ ǫ (x)) ≤ G ǫ (φ ♯ -η),
and G ǫ (φ ǫ (x)) goes to 0 = χ(x) thanks to [START_REF] Cancès | Convergence and a posteriori error analysis for energy-stable finite element approximations of degenerate parabolic equations[END_REF]. On the other hand, for x ∈ Ω ♭ such that φ ǫ (x) → φ ♭ and G ǫ (φ ǫ (x)) → χ(x), then for ǫ small enough, φ ǫ (x) ≤ φ ♭ +φ ♯

2

, so that G ǫ (x) ≤ 0. Therefore, χ(x) ≤ 0. Similarly, χ(x) ≥ 0 for almost all x ∈ Ω ♯ . Finally, we get that χ ∈ χ(φ) a.e. in Ω.

Let us now show that Ξ = ξ(s n ) and that the convergence (88) holds almost everywhere in Ω. Let s ∈ L ∞ (Ω) with 0 ≤ s ≤ 1 be arbitrary, then we infer from the aforementioned convergences that

(ξ ǫ (s n,ǫ ) -ξ(s))(s n,ǫ -s) -→ ǫ→0 (Ξ -ξ(s))(s n -s) weakly in L 2 (Ω).
On the other hand, one has

(ξ ǫ (s n,ǫ ) -ξ(s))(s n,ǫ -s) = (ξ ǫ (s n,ǫ ) -ξ(s n,ǫ ))(s n,ǫ -s) + (ξ(s n,ǫ ) -ξ(s))(s n,ǫ -s).
The first term in the above right-hand side tends uniformly to 0 because of the uniform convergence of ξ ǫ towards ξ which can be deduced from Assumption (H2) on γ and of the uniform boundedness of s n,ǫ , whereas the second term is non-negative because ξ is non-decreasing. Therefore, we obtain that (Ξξ(s))(s ns) ≥ 0 a.e. in Ω for any s ∈ [0, 1]. This implies that Ξ = ξ(s n ) (see for instance [START_REF] Andreianov | A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs[END_REF]). Next,

|ξ(s n,ǫ ) -ξ(s n )| ≤ |ξ(s n,ǫ ) -ξ ǫ (s n,ǫ )| + |ξ ǫ (s n,ǫ ) -ξ(s n )|.
The uniform convergence of ξ ǫ towards ξ implies that the first term in the right-hand side goes to 0 with ǫ uniformly on Ω. The second term converges to 0 almost everywhere owing to (86). As a consequence, ξ(s n,ǫ ) tends almost everywhere to ξ(s n ), and since ξ -1 is continuous, one gets that (89) s n,ǫ -→ ǫ→0 s n a.e. in Ω.

Bearing in mind the almost everywhere convergence of φ ǫ towards φ, then

φ n,ǫ = s n,ǫ φ ǫ -→ ǫ→0 s n φ = φ n a.e. in Ω.
The convergence properties (86), ( 87) and ( 88) are enough to pass to the limit in the linear relations to recover (81) and (84). Moreover, the energy estimate (85) holds true thanks to Corollary 3.3, to the lower semi-continuity of F (see for instance [START_REF] Droniou | The Gradient Discretisation Method[END_REF]Lemma C.6] for the capillary energy part) and the lower-semi continuity of the squared norms for the weak convergences. So our last focus is on the proof of (82), obtaining (83) being similar.

Starting from (63), one has that (90)

Ω φ n,ǫ -φ ⋆ n h v + Ω k ǫ (s n,ǫ ) µ n K(φ ǫ )∇ ((p n,ǫ + π ǫ + G ǫ (φ ǫ )) -ρ n g • (x + u ǫ )) • ∇v = 0 for all v ∈ V.
Fix v ∈ V , then because of (88) we readily get that ( 91)

Ω φ n,ǫ -φ ⋆ n h v -→ ǫ→0 Ω φ n -φ ⋆ n h v.
One also easily shows, thanks to Assumption (H3), to (86), (89) and Lebesgue's dominated convergence theorem that k ǫ (s n,ǫ )K(φ ǫ )∇v converges (strongly) in L 2 (Ω) d towards s n K(φ)∇v. Then we infer from (86) that

Ω k ǫ (s n,ǫ ) µ n K(φ ǫ )∇ ((π ǫ + G ǫ (φ ǫ )) -ρ n g • (x + u ǫ )) • ∇v (92) -→ ǫ→0 Ω s n µ n K(φ)∇ ((π + χ) -ρ n g • (x + u)) • ∇v.
The last term to be studied is

J ǫ = Ω k ǫ (s n,ǫ ) µ n K(φ ǫ )∇p n,ǫ • ∇v.
Recall that pn,ǫ ∈ p n,ǫ (s n,ǫ ), with

p n,ǫ (s) = (-∞, γ ǫ (0) + γ ′ ǫ (0)] if s = 0, f ǫ (s) := γ ǫ (s) + (1 -s)γ ′ ǫ (s) if s ∈ (0, 1]
, the function f ǫ being continuous and increasing. Therefore, s n,ǫ is a continuous function of pn,ǫ , i.e. s n,ǫ = Ŝǫ (p n,ǫ ) with

Ŝǫ (p) = 0 if p ≤ γ ǫ (0) + γ ′ ǫ (0), f -1
ǫ (p) otherwise, hence J ǫ can be rewritten as

J ǫ = Ω 1 µ n K(φ ǫ )∇ ψǫ (p n,ǫ ) • ∇v, with ψǫ (p) = p fǫ(0) k ǫ ( Ŝǫ (a))da.
Define also the nonlinear function

ψ ǫ : [0, 1] → R + by ψ ǫ (s) = ψǫ (f ǫ (s)) = s 0 k ǫ (a)(1 -a)γ ′′ ǫ (a)da, s ∈ [0, 1].
It converges uniformly on [0, 1] towards ψ defined in [START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF] and therefore it is uniformly bounded in [0, 1] thanks to Assumption (H2). The function ψ ǫ (s n,ǫ ) -ψǫ (p n,ǫ ) = ǫ(γ ǫ (0) + γ ′ ǫ (0)pn,ǫ ) + vanishes on {s n,ǫ > 0} and, thus, owing to estimate (71), it satisfies

(93) ∇ ψ ǫ (s n,ǫ ) -ψǫ (p n,ǫ ) 2 L 2 (Ω) d ≤ Cǫ 2 ∇p n,ǫ 2
L 2 (Ω) d ≤ Cǫ for some C possibly depending on h > 0 but not on ǫ. Moreover, since 0 ≤ s n,ǫ ≤ k ǫ (s n,ǫ ) ≤ 1 and k ǫ (s n,ǫ ) ≤ s n,ǫ + ǫ, and recalling [START_REF] Brenner | Linear finite element methods for planar linear elasticity[END_REF], we infer from Proposition 3.1 that

∇ψ ǫ (s n,ǫ ) 2 L 2 (Ω) d ≤ ∇ ψǫ (p n,ǫ ) 2 L 2 (Ω) d ≤ Ω k ǫ (s n,ǫ )|∇p n,ǫ | 2 ≤ ∇ξ ǫ (s n,ǫ ) 2 L 2 (Ω) d + ǫ ∇p n,ǫ 2
L 2 (Ω) d ≤ C.
This implies that there exists some Ψ ∈ H 1 (Ω) such that (up to a subsequence as usual)

ψ ǫ (s n,ǫ ) -→ ǫ→0
Ψ weakly in H 1 (Ω) and almost everywhere.

Then because of (89) and of the uniform convergence of ψ ǫ towards ψ, we can identify Ψ as ψ(s n ) thanks to arguments similar to those used previously for showing that Ξ = ξ(s n ). Finally, (93) shows that

∇ ψǫ (p n,ǫ ) -→ ǫ→0 ∇ψ(s n ) weakly in L 2 (Ω) d .
Since K(φ ǫ )∇v converges strongly in L 2 (Ω) d towards K(φ)∇v thanks to Assumption (H3), we can pass to the limit in J ǫ (94)

J ǫ -→ ǫ→0 Ω 1 µ n K(φ)∇ψ(s n ) • ∇v.
The combination of (91), ( 92) and (94) in (90) gives (82), and ends the proof of Proposition 3.4.

Concluding the proof of Theorem 1.3

Let X 0 = (φ 0 , u 0 , θ 0 , π 0 = M θ 0 ) be as in (H4) and fix h > 0. Then applying Proposition 3.4 recursively, we get sequences (X j ) j≥0 and (χ j ) j≥1 of measurable functions such that, for all j ≥ 1, φ j = (φ j n , φ j w ) ∈ K φ a.e. in Ω, with φ j = φ j n + φ j w ∈ H 1 (Ω) and ξ(s j n )ξ(s D n ) ∈ V setting s j n = φ j n /φ j , such that u j ∈ U , such that θ j , π j , χ j ∈ H 1 (Ω) with χ j ∈ χ(φ j ) a.e. in Ω, fulfilling (95)

π j = M θ j and φ j -b ∇ • u j -θ j = φ r , such that Ω φ j n -φ j-1 n h w - Ω 1 µ n K(φ j ) ∇ψ(s j n ) + s j n ∇ (π j + χ j ) -ρ n g • (x + u j ) • ∇w = 0, ∀w ∈ V, (96)
Ω φ j w -φ j-1 w h w - Ω 1 µ w K(φ j ) -∇ψ(s j n ) + s j w ∇ (π j + χ j ) -ρ w g • (x + u j ) • ∇w = 0, ∀w ∈ V, (97)
and ( 98)

Ω 2µ ε(u j ) : ε(v) + λ(∇ • u j )(∇ • v) = Ω bπ j ∇ • v + Ω f (φ j ) • v, ∀v ∈ V d .
Moreover, we deduce from the energy estimate (85) that (99)

F (X j ) + C 14 h ξ(s j n ) 2 H 1 (Ω) + π j 2 H 1 (Ω) + χ j 2 H 1 (Ω) ≤ (1 + 2C 10 h) F (X j-1 ) + Ω f (φ j ) • (u j -u j-1 ) + α∈{n,w} Ω (φ j α -φ j-1 α ) p D α + C 15 h
with constants C 10 , C 14 and C 15 depending neither on j nor on h, which in particular also provides the uniform bound for the phase porosity increments in the dual norm (100)

φ j α -φ j-1 α V ′ ≤ C 16 h ∀α ∈ {n, w} and j ≥ 1,
for a positive constant C 16 not depending on j and h. We define the piecewise constant in time approximations X h = (φ h = (φ n,h , φ w,h ), u h , θ h , π h ) and χ h by setting

X h (t, x) = X j (x), χ h (t, x) = χ j (x) for (t, x) ∈ (j -1)h, jh × Ω.
We also denote by φ h = φ n,h + φ w,h and s α,h = φ α,h /φ h for α ∈ {n, w}. We extend X h to negative times by setting X h (t) = X 0 for t ≤ 0. Moreover, we define the approximate time derivatives ∂ h t φ α,h , α ∈ {n, w}, and ∂ h t φ h by

∂ h t φ α,h (t, •) = φ j α -φ j-1 α h for t ∈ (j -1)h, jh , and
∂ h t φ h = ∂ h t φ n,h + ∂ h t φ w,h .
We will now let h tend to 0. Once again, our proof relies on compactness arguments. In what follows, the limit h → 0 implicitly refers to the convergence when ℓ goes to +∞ of a decreasing sequence (h ℓ ) ℓ≥0 of positive times steps tending to 0.

Let T > 0 be an arbitrary finite time horizon, then summing (99) for j = 1, . . . , ⌈ T h ⌉ leads to

F (X h (T ))+ C 14 ξ(s n,h ) 2 L 2 ((0,T );H 1 (Ω)) + π h 2 L 2 ((0,T );H 1 (Ω)) + χ h 2 L 2 ((0,T );H 1 (Ω)) (101)
≤ F (X 0 ) + 2C 10

⌈ T h ⌉ j=1 hF (X j-1 ) + (1 + 2hC 10 ) ⌈ T h ⌉ j=1 Ω f (φ j ) • (u j -u j-1 ) + α∈{n,w} Ω (φ α,h (T ) -φ 0 α ) p D α + C 15 (T + 1)
Employing summation by parts, used that the summation can be extended to index j = 0 corresponding to zero due to the above extension to negative times, utilizing that f is uniformly bounded and affine due to (H8), and employing the uniform bound (100), we obtain for a suitable constant

C > 0 that ⌈ T h ⌉ j=1 Ω f (φ j ) • (u j -u j-1 ) (102) = Ω f (φ h (T )) • u h (T ) -f (φ 0 ) • u 0 - ⌈ T h ⌉ j=1 Ω f (φ j ) -f (φ j-1 ) • u j-1 ≤ Ω f (φ h (T )) • u h (T ) -f (φ 0 ) • u 0 + ⌈ T h ⌉ j=1 α∈{n,w} ρ a φ j α -φ j-1 α V ′ g • u j-1 V ≤ C + 1 4 F (X h (T )) + C 10 ⌈ T h ⌉ j=1
hF (X j-1 ).

Furthermore, using the uniform boundedness of φ and of p D results in

(103) α∈{n,w} Ω (φ j α -φ j-1 α ) p D α ≤ C ′
for a suitable constant C ′ > 0. Thus, combining (101)-( 103), and assuming 2hC 10 ≤ 1 results in

1 2 F (X h (T ))+ C 14 ξ(s n,h ) 2 L 2 ((0,T );H 1 (Ω)) + π h 2 L 2 ((0,T );H 1 (Ω)) + χ h 2 L 2 ((0,T );H 1 (Ω)) (104)
≤ F (X 0 ) + 4C 10

⌈ T h ⌉ j=1 hF (X j-1 ) + C,
for C depending on T but not on h. Then we deduce from a discrete Gronwall lemma that, for all T ≥ 0, there holds (105)

F (X h (T ))+ ξ(s n,h ) 2 L 2 ((0,T );H 1 (Ω)) + π h 2 L 2 ((0,T );H 1 (Ω)) + χ h 2 L 2 ((0,T );H 1 (Ω)) ≤ C(1+F (X 0 ))
, where C depends on T but not on h. Assumption (H4) ensures that F (X 0 ) is finite, whereas we deduce from Assumption (H6), from (105) and Korn's inequality that

(106) u h L ∞ ((0,T );H 1 (Ω) d ) + ∇(∇ • u h ) L 2 ((0,T )×Ω) d ≤ C.
Then we deduce from the constraints (95) that (107)

θ h L 2 ((0,T );H 1 (Ω)) ≤ C, φ h L 2 ((0,T );H 1 (Ω)) ≤ C. Since ψ • ξ -1 is 1 2 -Lipschitz continuous, one further deduces from (105) that (108) ψ(s n,h ) L 2 ((0,T );H 1 (Ω)) ≤ C.
Let v : R + × Ω → R be a smooth function such that v(t, x) = 0 for (t, x) ∈ R + × Γ D as well as for t ≥ T . For h > 0 fixed, we define v j j≥0 ⊂ V by

v j = 1 h jh (j-1)h v(t, •)dt for j ≥ 1 and v 0 = v(0, •),
and by v h (t, x) = v j (x) for (t, x) ∈ (j -1)h, jh × Ω.

Choosing w = hv j in (96) and summing over j ≥ 1 gives (109

) R+×Ω ∂ h t φ n,h v = R+×Ω 1 µ n K(φ h ) (∇ψ(s n,h ) + s n,h (∇(π h + χ h ) -ρ n g • (x + u h ))) • ∇v.
Then we infer from (105) and (108) together with 0 ≤ s n,h ≤ 1 and the uniform boundedness of K(φ h ) deduced from Assumption (H3) that (110)

∂ h t φ n,h L 2 ((0,T );V ′ ) ≤ C, ∂ h t φ w,h L 2 ((0,T );V ′ ) ≤ C
, the second estimate being similar, and thus that (111)

∂ h t φ h L 2 ((0,T );V ′ ) ≤ C.
We infer from ( 105) and ( 107) that there exists Ξ, π, χ and θ in L 2 ((0, T ); H 1 (Ω)) such that, up to the extraction of a subsequence, there holds

ξ(s n,h ) -→ h→0 Ξ weakly in L 2 ((0, T ); H 1 (Ω)), (112) π h -→ h→0 π weakly in L 2 ((0, T ); H 1 (Ω)), (113) χ h -→ h→0 χ weakly in L 2 ((0, T ); H 1 (Ω)), (114)
θ h -→ h→0 θ weakly in L 2 ((0, T ); H 1 (Ω)), (115)
whereas (106) ensures the existence of some u ∈ L ∞ ((0, T ); H 1 (Ω) d ) ∩ L 2 ((0, T ); U ) such that, up to a subsequence, (116) u h -→ h→0 u in the L ∞ ((0, T ); H 1 (Ω) d )-weak-⋆ and L 2 ((0, T ); U )-weak sense.

Since φ h belongs to the bounded convex subset K φ a.e. in R + × Ω, we get that (117)

φ h -→ h→0 φ in the L ∞ (R + × Ω) 2 -weak-⋆ sense
with φ(t, x) ∈ K φ for a.e. (t, x) ∈ R + × Ω. We denote by φ = φ n + φ w , so that, owing to (107), (118) φ h -→ h→0 φ in the L ∞ (R + × Ω)-weak-⋆ and L 2 ((0, T ); H 1 (Ω))-weak senses, and φ ♭ ≤ φ ≤ φ ♯ almost everywhere. Note that at this point, we have already enough material to pass to the limit in the linear relations

π h = M θ h , φ h -b∇ • u h -θ h = φ r and (0,T )×Ω 2µ ε(u h ) : ε(v)+λ(∇•u h )(∇•v) = (0,T )×Ω b π h ∇•v+ Ω f (φ h )•v, v ∈ L 2 ((0, T ); V d )
to recover (5), ( 6), as well as [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]. What remains in this section is devoted to the proof of the weak forms ( 23) and ( 24) of the phase volume conservation equations (1). We deduce from (110) the existence of δ n , δ w ∈ L 2 ((0, T ); V ′ ) and of δ = δ n + δ w such that

∂ h t φ α,h -→ h→0 δ α and ∂ h t φ h -→ h→0
δ weakly in L 2 ((0, T ); V ′ ), α ∈ {n, w}.

Let ϕ ∈ C ∞ c ((0, T )×Ω), and assume that ϕ(t) = 0 for t ≤ h and t ≥ T -h, which is not a restriction since h is going to tend to 0. Then for α ∈ {n, w}, there holds We further deduce from ( 118) and (119) that one can apply a discrete (here in time only) Aubin-Simon type result (see for instance [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs-application to a turbulence model[END_REF] or [START_REF] Andreianov | A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs[END_REF]), to get that, up to a subsequence, (120) φ h -→ h→0 φ a.e. in (0, T ) × Ω.

Denote by s n the L ∞ ((0, T )×Ω)-weak-⋆ limit of s n,h , then the convergences (117) and (118) show that s n = φ n /φ. In order to show that s n,h converges pointwise towards s n , we take inspiration on what is done in [START_REF] Bonaldi | Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media[END_REF]. Let f : [0, 1] → R be a continuous increasing function such that f • ψ -1 and f are both Lipschitz continuous, so that, thanks to (108), there holds Bearing in mind the definition of φ n,h , we get that

A τ h ≤ L f φ ♭ (0,T -τ )×Ω ⌈ t+τ h ⌉ j=⌈ t h ⌉+1
φ j nφ j-1 n f (s n,h (t + τ ))f (s n,h (t)) .

Then using (96) leads to

A τ h ≤ h L f φ ♭ (0,T -τ )×Ω ⌈ t+τ h ⌉ j=⌈ t h ⌉+1 F j n ∇ f (s n,h (t + τ )) -f (s n,h (t))
where the flux

F j n = - 1 µ n K(φ j ) ∇ψ(s j n ) + s j n ∇ (π j + χ j ) -ρ n g • (x + u j )
satisfies

Ω ⌈ T τ ⌉ j=1 h F j n 2 ≤ C
in view of the a priori estimates (105) and (108). Applying the Cauchy-Schwarz inequality gives

A τ h ≤ A τ,(1) h A τ, (2) h , where A τ,(1) h
= 2 L f φ ♭   ⌈ t+τ h ⌉ j=⌈ t h ⌉+1 h   f (s n,h ) 2 L 2 ((0,T );H 1 (Ω)) (121)
≤ C(τ + h), and where

A τ,(2) h ≤ L f φ ♭ (0,T -τ )×Ω ⌈ t+τ h ⌉ j=⌈ t h ⌉+1 h F j n 2 ≤ Cτ.
As a consequence, we can use the Riesz-Fréchet-Kolmogorov compactness criterion [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]Theorem 4.26] to claim that, up to a subsequence, f (s n,h ) -→ h→0 f a.e. in (0, T ) × Ω.

As f is increasing, this implies that s n,h converges also almost everywhere, and by uniqueness of the limit we have that (122) s n,h -→ h→0 s n a.e. in (0, T ) × Ω.

In combination with (108), this ensures that, up to a subsequence, there holds (123) ψ(s n,h ) -→ h→0 ψ(s n ) a.e. and weakly in L 2 ((0, T ); H 1 (Ω)).

We can now pass to the limit in (109). Thanks to (119), one has

R+×Ω ∂ h t φ n,h v -→ h→0 R+ ∂ t φ α , v V ′ ,V .
Moreover, the convergences (120) and (122) together with Assumption (H3) on K show that K(φ h )∇v -→ h→0 K(φ)∇v and s n,h K(φ h )∇v -→ In view of (114), ( 118) and ( 124), letting h tend to 0 gives (0,T )×Ω (φa)(χ -A)ϕ ≥ 0, ∀ϕ ∈ L 2 ((0, T ) × Ω), ϕ ≥ 0.

This implies that (φa)(χ -A) ≥ 0 a.e. in (0, T ) × Ω, and since (a, A) is arbitrary in the graph χ, then χ ∈ χ(φ) a.e. in (0, T ) × Ω. This concludes the proof of Theorem 1.3.

Conclusion and prospects

In this paper we showed the first global existence for degenerate multiphase poromechanics under some reasonable assumption of the domain Ω, i.e. Ω is homogeneous in space and is such that the elasticity equation enjoy some regularity (H6). The model includes, among others, capillarity effects, nonlinear pore pressure constitutive law, porosity dependent permeability and gravitational forces, in a thermodynamically consistent manner. For simplicity, we have restricted our analysis to the case of linear relative permeabilities, even though our purpose should transpose to mild nonlinear relative permeabilities. Some weak-coupling condition between the flow and the mechanics is also postulated, cf. (H7). This assumption was necessary in our analysis to get estimates on ∇π and ∇ξ separately. The assumption then becomes needless in the case where φ ♭ < φ < φ ♯ , yielding χ = 0, but we are not able to guaranty this situation a priori. The main a priori estimate is derived from the time evolution of the Helmholtz free energy, which is decreasing up to contributions coming from the boundary and from gravity.

The extension of our global existence result to more complex frameworks is an open problem. Mathematical difficulties have to be bypassed here, in particular, as the proofs of Propositions 2.6 and 3.1 strongly rely on the fact that the problem is essentially homogeneous (or at least smoothly varying) in space and on the weak coupling assumption (H7). Extending the proof of Brenner and Sung [START_REF] Brenner | Linear finite element methods for planar linear elasticity[END_REF] to prove (H6) is also a challenging problem which can be of interest in other contexts.

Finally, numerical experiments should be carried out in future works. The rigorous proof of the convergence of numerical methods could even be achieved thanks to compactness arguments adapting ours.

  ) The boundary data p D = (p D n , p D w ) ∈ W 1,∞ (Ω) 2 do not depend on time. (H6) The domain Ω is a bounded and connected open subset of R d . Its boundary is Lipschitz continuous and splits into ∂Ω = Γ D ∪ Γ N with Γ D ∩ Γ N = ∅.

  (121) f (s n,h ) L 2 ((0,T );H 1 (Ω)) ≤ C for some C not depending on h. Let τ > 0 and y ∈ R d , and denote by Ω y = {x ∈ Ω | [x, y] ⊂ Ω}, then (0,T -τ )×Ωy

h→0s

  n K(φ)∇v strongly in L 2 (R + × Ω).Combining this with the weak convergences (113), (114), (123), and (116), one gets thatR+×Ω 1 µ n K(φ h ) (∇ψ(s n,h ) + s n,h ∇ ((π h + χ h )ρ n g • (x + u h ))) • ∇v ) (∇ψ(s n ) + s n ∇ ((π + χ)ρ n g • (x + u))) • ∇v,so that (109) gives (23) at the limit h → 0. Similar arguments allow to recover[START_REF] Buzzi | Interface conditions for degenerate two-phase flow equations in one space dimension[END_REF]. The last step to conclude the proof of Theorem 1.3 is to show that χ ∈ χ(φ). Combining (111) with (114) and (118), we can apply [4, Proposition 3.8] which gives that (124)(0,T )×Ω φ h χ h ϕ -→ h→0 (0,T )×Ω φ χ ϕfor all ϕ ∈ L 2 ((0, T ) × Ω). This allows to use Minty's trick to show that χ ∈ χ(φ) a.e. in Ω. Indeed, χ h ∈ χ(φ h ) is equivalent to the fact that, for all (a, A)∈ [φ ♭ , φ ♯ ] × R such that A ∈ χ(a), then (φ ha)(χ h -A) ≥ 0.Therefore, for all ϕ ∈ L 2 ((0, T ) × Ω) with ϕ ≥ 0, one has (0,T )×Ω (φ ha)(χ h -A)ϕ ≥ 0.

  2 and ũ ∈ U . Together with φ α,k ∈ K φ and this results in a uniform bound in k. This ensures the existence of some p ∈ p D +V 2 , and π, θ ∈ H 1 (Ω), such that, up to a subsequence,

	(60a)	p k -→
	(60b)	π k -→ k→+∞	π a.e. in Ω and weakly in H 1 (Ω),
		θ k -→

k→+∞ p a.e. in Ω and weakly in H 1 (Ω) 2 , k→+∞ θ a.e. in Ω and weakly in H 1 (Ω),

  ∈ (j -1)h, jh × Ω.Since ϕ is assumed to be regular, one readily checks that ∂h t ϕ converges uniformly towards ∂ t ϕ, hence

		(0,T )×Ω	∂ h t φ α,h ϕ = -	(0,T )×Ω	φ α,h	∂h t ϕ
	where we have set
		∂h t ϕ(t, x) = dτ if (t, x) (0,T )×Ω 1 h jh (j-1)h ϕ(τ + h, x) -ϕ(τ, x) h T ∂ h t φ α,h ϕ -→ h→0 0 δ α , ϕ V ′ ,V = -	(0,T )×Ω	φ α ∂ t ϕ,
	so that δ α = ∂ t φ α , thus
	(119)	∂ h t φ α,h -→

h→0 ∂ t φ α and ∂ h t φ h -→ h→0 ∂ t φ weakly in L 2 ((0, T ); V ′ ), α ∈ {n, w}.

  Because of (121) and of the classical characterization[START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF] Proposition 9.3] of the space H 1 (Ω), one can estimate B y h by B y h ≤ C|y| 2 for some C not depending on h. Let us now adapt the ideas of[START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF] to the term A τ . As f is Lipschitz continuous with Lipschitz constant denoted by L f , there holds (t + τ )s n,h (t) f (s n,h (t + τ ))f (s n,h (t)) (t + τ )φ n,h (t) f (s n,h (t + τ ))f (s n,h (t)) .

		2 ≤ 2(A τ h + B y h )
	with
	A τ h =	2 ,
	B y h =	(0,T )×Ωy
	A τ h ≤ L f s n,h ≤ (0,T -τ )×Ω L f (0,T -τ )×Ω φ ♭ φ n,h

f (s n,h (t + τ, x + y))f (s n,h (t, x)) (0,T -τ )×Ω f (s n,h (t + τ, x))f (s n,h (t, x)) f (s n,h (t, x + y))f (s n,h (t, x))

2 .

Acknowledgements. JWB and CC warmly thank Konstantin Brenner for stimulating exchanges on multiphase poromechanics. CC also thanks Jérôme Droniou and Roland Masson for enlightening discussions. This work was done in the framework of the GradFlowPoro project supported by the Campus France (project N°48356QE) and Research council of Norway (RCN project #331960). CC also acknowledges partial support from Labex CEMPI (ANR-11-LABX-0007-01). JWB also acknowledges partial support from the UoB Akademia-project FracFlow.

Statements. The authors have no conflict of interest related to this work. This manuscript has no associated data.