J Droniou
email: jerome.droniou@monash.edu.

R Eymard

Uniform-in-time convergence result of numerical methods for non-linear parabolic equations

Keywords:

We prove that all Gradient Schemes -which include Finite Element, some Mixed Finite Element and Finite Volume methods -converge uniformly in time when applied to a family of nonlinear parabolic equations which contains the Richards, Stefan and Leray-Lions models.

Introduction 1.Motivation

The following generic nonlinear parabolic model ∂ t β(u) -div (a(x, ν(u), ∇ζ(u))) = f in Ω × (0, T ), β(u)(x, 0) = β(u ini )(x) in Ω, ζ(u) = 0 on ∂Ω × (0, T ) [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF] where β, ζ are non-decreasing, ν is such that ν = β ζ and a is a Leray-Lions operator, arises in various frameworks (see next section for precise hypotheses on the data). This model includes 1. Richards' model, setting ζ(s) = s, ν = β and a(x, ν(u), ∇ζ(u)) = Λ(x)K(x, β(u))∇u, which describes the flow of water in a heterogeneous anisotropic underground medium, 2. Stefan's model [START_REF] Bertsch | The Stefan problem with heating: appearance and disappearance of a mushy region[END_REF], setting β(s) = s, ν = ζ, a(x, ν(u), ∇ζ(u)) = K(ζ(u))∇ζ(u), which arises in the study of a simplified heat diffusion in a melting medium, 3. p-Laplace problem (and p-Laplace-like problems), setting β(s) = ζ(s) = ν(s) = s and a(x, ν(u), ∇ζ(u)) = |∇u| p-2 ∇u, which is involved in the motion of glaciers [START_REF] Glowinski | Approximation of a nonlinear elliptic problem arising in a nonnewtonian fluid flow model in glaciology[END_REF] or flows of incompressible turbulent fluids through porous media [START_REF] Diaz | On a nonlinear parabolic problem arising in some models related to turbulent flows[END_REF].

The numerical approximation of these models has been extensively studied in the literature (see the fundamental work on the Stefan problem [START_REF] Nochetto | Approximation of degenerate parabolic problems using numerical integration[END_REF], and [START_REF] Pop | Numerical schemes for degenerate parabolic problems[END_REF] for a review of some numerical approximations, see [START_REF] Maitre | Numerical analysis of nonlinear elliptic-parabolic equations[END_REF] for the Richards problem and see [START_REF] Droniou | Finite volume schemes for fully non-linear elliptic equations in divergence form[END_REF][START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF] and references therein for some studies of convergence of numerical methods for the Leray-Lions problem). However, the convergence analysis of the considered schemes received a much reduced coverage and consists mostly in establishing space-time averaged results (e.g. in L 2 (Ω × (0, T ))), see for example [START_REF] Eymard | Finite volumes and nonlinear diffusion equations[END_REF][START_REF] Eymard | The finite volume method for Richards equation[END_REF]). Yet, the quantity of interest is often not u on Ω × (0, T ) but u at a given time, for example t = T . Existing numerical analysis results therefore do not ensure that this quantity of interest is indeed properly approximated by numerical methods.

The usual way to obtain pointwise-in-time approximation results for numerical schemes is to prove estimates in L ∞ (0, T ; L 2 (Ω)) on u -u, where u is the approximated solution. Establishing such error estimates is however only feasible when uniqueness of the solution u to (1) can be proved, which is the case for Richards' and Stefan's problems but not for more complex non-linear parabolic problems as [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF] or even p-Laplace problems. It moreover requires some regularity assumptions on u, which clearly fail for (1) (and simpler p-Laplace problems) for which, because of the possible plateaux of β and ζ, the solution can develop jumps in its gradient. The purpose of this article is to prove that, using Discrete Functional Analysis techniques (i.e. the translation to numerical analysis of nonlinear analysis techniques), an L ∞ (0, T ; L 2 (Ω)) convergence result can be established for numerical approximations of (1), without having to assume nonphysical regularity assumptions on the data. Note that, although Richards' and Stefan's models are formally equivalent when β and ζ are strictly increasing (consider β = ζ -1 to pass from one model to the other), they change nature when these functions are allowed to have plateaux. Stefan's model can degenerate to an ODE (if ζ is constant on the range of the solution) and Richards' model can become a non-transient elliptic equation (if β is constant on this range). The innovative technique we develop in this paper is nonetheless generic enough to work directly on [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF] and with a vast number of numerical methods. That being said, a particular numerical framework must be selected in order to write precise equations and estimates. The framework we choose is that of Gradient Schemes, which has the double benefit of covering a vast number of numerical methods and of having already been studied for many models -elliptic, parabolic, linear or non-linear, possibly degenerate, etc. -with various boundary conditions. The schemes or family of schemes included in the Gradient Schemes framework, and to which our results therefore directly apply, currently are:

• Galerkin methods, including conforming Finite Element schemes,

• Finite Element with mass lumping [START_REF] Ciarlet | The finite element method[END_REF],

• The Crouzeix-Raviart non-conforming Finite Element, with or without mass lumping [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF],

• The Raviart-Thomas Mixed Finite Elements [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF],

• The Vertex Approximate Gradient scheme [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF],

• The Hybrid Mimetic Mixed family [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF], which includes Mimetic Finite Differences [START_REF] Brezzi | A family of mimetic finite difference methods on polygonal and polyhedral meshes[END_REF], Mixed Finite Volume [START_REF] Droniou | A mixed finite volume scheme for anisotropic diffusion problems on any grid[END_REF] and the SUSHI scheme [START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces[END_REF],

• The Discrete Duality Finite Volume scheme in dimension 2 [START_REF] Hermeline | Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes[END_REF][START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes[END_REF], and the CeVeFE-Discrete Duality Finite Volume scheme in dimension 3 [START_REF] Coudière | A 3d discrete duality finite volume method for nonlinear elliptic equations[END_REF],

• The Multi-Point Flux Approximation O-method [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF][START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF].

We refer the reader to [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF][START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF][START_REF] Eymard | Gradient schemes for the Stefan problem[END_REF][START_REF] Eymard | Gradient scheme approximations for diffusion problems[END_REF][START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF] for more details. Let us finally emphasize that the unified convergence study of numerical schemes for Problem [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF], which combines a general Leray-Lions operator and nonlinear functions β or ζ, seems to be new even without the uniform-in-time convergence result.

The paper is organised as follows. In Section 1.2, we present the assumptions and the notion of weak solution for [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF] and, in Section 1.3, we give an overview of the ideas involved in the proof ζ (q)β (q)dq. (2d)

a : Ω × R × R d → R d , with p ∈ (1, +∞), is a Caratheodory function, (2e)
(i.e. a function such that, for a.e. x ∈ Ω, (s, ξ) → a(x, s, ξ) is continuous and, for any (s, ξ) ∈ R × R d , x → a(x, s, ξ) is measurable) ∃a ∈ (0, +∞) : a(x, s, ξ) • ξ ≥ a|ξ| p , for a.e. x ∈ Ω, ∀s ∈ R, ∀ξ ∈ R d ,

(a(x, s, ξ) -a(x, s, χ)) • (ξ -χ) ≥ 0, for a.e. x ∈ Ω, ∀s ∈ R, ∀ξ, χ ∈ R d , (2g) ∃a ∈ L p (Ω) , ∃µ ∈ (0, +∞) : |a(x, s, ξ)| ≤ a(x) + µ|ξ| p-1 , for a.e. x ∈ Ω, ∀s ∈ R, ∀ξ ∈ R d ,

and

u ini ∈ L 2 (Ω), f ∈ L p (Ω × (0, T )). (2i)
We denote by R β the range of β and define the pseudo-inverse function β r : R β → R of β by

∀s ∈ R β , β r (s) = inf{t ∈ R | β(t) = s} if s ≥ 0, sup{t ∈ R | β(t) = s} if s < 0, = closest t to 0 such that β(t) = s. (3)
Since β(0) = 0, we notice that β r ≥ 0 on R β ∩ R + and β r ≤ 0 on R β ∩ R -. We then define B : R β → [0, ∞] by

B(z) = z 0 ζ(β r (s)) ds.
Since β r is non-decreasing, this expression is always well-defined in [0, ∞). The signs of β r also ensure that that B is non-decreasing on R β ∩ R + and non-increasing on R β ∩ R -. We can therefore extend B to R β by these limits (possibly +∞) at the potential endpoints of R β .

The precise notion of solution to (1) that we consider is then the following:

                       u ∈ L p (0, T ; L p (Ω)) , ζ(u) ∈ L p (0, T ; W 1,p 0 (Ω)) , B(β(u)) ∈ L ∞ (0, T ; L 1 (Ω)), β(u) ∈ C([0, T ]; L 2 (Ω)-w), ∂ t β(u) ∈ L p (0, T ; W -1,p (Ω)), β(u)(•, 0) = β(u ini ) in L 2 (Ω), T 0 ∂ t β(u)(•, t), v(•, t) W -1,p ,W 1,p 0 dt + T 0 Ω a(x, ν(u(x, t)), ∇ζ(u)(x, t)) • ∇v(x, t)dxdt = T 0 Ω f (x, t)v(x, t)dxdt , ∀v ∈ L p (0; T ; W 1,p 0 (Ω)). ( 4
)
where C([0, T ]; L 2 (Ω)-w) denotes the space of continuous functions [0, T ] → L 2 (Ω) for the weak- * topology of L 2 (Ω). Here and in the following, we denote by p the dual exponent p p-1 to p and we remove the mention of Ω in the duality bracket

•, • W -1,p ,W 1,p 0 = •, • W -1,p (Ω),W 1,p 0 (Ω) . Remark 1.1 The derivative ∂ t β(u) is to be understood in the usual sense of distributions on Ω × (0, T ). Since the set T = { q i=1 ϕ i (t)γ i (x) : q ∈ N, ϕ i ∈ C ∞ c (0, T ), γ i ∈ C ∞ c (Ω)} of tensorial functions in C ∞ (Ω × (0, T )) is dense in L p (0, T ; W 1,p 0 (Ω))
, one can ensure that this distribution derivative ∂ t β(u) belongs to L p (0, T ; W -1,p (Ω)) = (L p (0, T ; W 1,p 0 (Ω)) by checking that the linear form

ϕ ∈ T → ∂ t β(u), ϕ D ,D = - T 0 Ω β(u)(x, t)∂ t ϕ(x, t)dxdt
is continuous for the norm of L p (0, T ; W 1,p 0 (Ω)). Note that the continuity property of β(u) natural. Indeed, the PDE in the sense of distributions shows that T ϕ : t → β(u)(t), ϕ L 2 belongs to W 1,1 (0, T ), and is therefore continuous, for any ϕ ∈ C ∞ c (Ω). The density in L 2 (Ω) of such ϕ, combined with the fact that β(u) ∈ L ∞ (0, T ; L 2 (Ω)) (coming from B(β(u)) ∈ L ∞ (0, T ; L 1 (Ω)) and ( 26)), proves the continuity of T ϕ for any ϕ ∈ L 2 (Ω), that is to say the continuity of β(u) : [0, T ] → L 2 (Ω)-w. This notion of β(u) as a function continuous in time is nevertheless a subtle one. It is to be understood in the sense that the function (x, t) → β(u(x, t)) has an a.e. representative which is continuous [0, T ] → L 2 (Ω)-w. In other words, there is a function Z ∈ C([0, T ]; L 2 (Ω)-w) such that Z(t)(x) = β(u(x, t)) for a.e. (x, t) ∈ Ω × (0, T ). We must however make sure, when dealing with pointwise values in time to separate Z from β(u(•, •)) as β(u(•, t 1 )) may not make sense for a particular t 1 ∈ [0, T ]. That being said, in order to adopt a simple notation, we will denote by β(u)(•, •) the function Z, and by β(u(•, •)) the a.e.-defined composition of β and u. Hence, it will make sense to talk about β(u)(•, t) for a particular t 1 ∈ [0, T ], and we will only write β(u)(x, t) = β(u(x, t)) for a.e. (x, t) ∈ Ω × (0, T ). Note that from this a.e. equality we can ensure that β(u)(•, •) takes its values in the closure R β of the range of β.

General principle for the uniform-in-time convergence result

As explained in the introduction, the main innovative result of this article is the uniform-in-time convergence result (Theorem 2.14 below). Although it's stated and proved in the context of numerical approximations of (1), we emphasize that its principle is also applicable to theoretical analysis of PDEs. Let us informally present this principle on the following continuous approximation of (1):

∂ t β(u ε ) -div (a ε (x, ν(u ε ), ∇ζ(u ε ))) = f in Ω × (0, T ), β(u ε )(x, 0) = β(u ini )(x) in Ω, ζ(u ε ) = 0 on ∂Ω × (0, T ) (5)
where a ε satisfies Assumptions (2e)-(2h) with constants not depending on ε and a ε → a pointwise as ε → 0. We want to show here how to deduce from averaged convergences a strong uniform-in-time convergence result. We therefore assume the following convergences (up to a subsequence as ε → 0), which are compatible with basic compactness results that can be obtained on (u ε ) ε and also correspond to the initial convergences [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] that can be obtained on numerical approximations of (1):

β(u ε ) → β(u) in C([0, T ]; L 2 (Ω)-w) , ν(u ε ) → ν(u) strongly in L 1 (Ω × (0, T )), ζ(u ε ) → ζ(u) weakly in L p (0, T ; W 1,p 0 (Ω)) , a ε (•, ν(u ε ), ∇ζ(u ε )) → a(•, ν(u), ∇ζ(u)) weakly in L p (Ω × (0, T )) d . (6)
We will prove from these convergences that, along the same subsequence, ν(u ε ) → ν(u) strongly in C([0, T ]; L 2 (Ω)), which is our uniform-in-time convergence result. We start by noticing that the weak-in-space uniform-in-time convergence of β(u ε ) gives, for any T 0 ∈ [0, T ] and any family (T ε ) ε>0 converging to

T 0 as ε → 0, β(u ε )(T ε , •) → β(u)(T 0 , •) weakly in L 2 (Ω).
Classical strong-weak semi-continuity properties of convex functions (see Lemma 3.4) and the convexity of B (see Lemma 3.3) then ensure that

Ω B(β(u)(x, T 0 ))dx ≤ lim inf ε→0 Ω B(β(u ε )(x, T ε ))dx. ( 7
)
The second step is to notice that, by (2g) for a ε ,

Tε 0 Ω [a ε (•, ν(u ε ), ∇ζ(u ε ) -a ε (•, ν(u ε ), ∇ζ(u))] • [∇ζ(u ε ) -∇ζ(u)] dxdt ≥ 0.
Developing this expression and using the convergences (6), we find that

lim inf ε→0 Tε 0 Ω a ε (•, ν(u ε ), ∇ζ(u ε ) • ∇ζ(u ε )(x, t)dxdt ≥ T0 0 Ω a(•, ν(u), ∇ζ(u)) • ∇ζ(u)dxdt. ( 8
)
We then establish the following formula:

Ω B(β(u ε (x, T ε )))dx + Tε 0 Ω a ε (x, ν(u ε (x, t)), ∇ζ(u ε )(x, t)) • ∇ζ(u ε )(x, t)dxdt = Ω B(β(u ini (x)))dx + Tε 0 Ω f (x, t)ζ(u ε )(x, t)dxdt. (9)
This energy estimate is formally obtained by multiplying (5) by ζ(u ε ) and integrating by parts, using the fact that (B • β) = ζβ (see Lemma 3.3). The rigorous justification of ( 9) is however quite technical, see Lemma 3.6 and Corollary 3.8. Thanks to [START_REF] Coudière | A 3d discrete duality finite volume method for nonlinear elliptic equations[END_REF], we can pass to the lim sup in [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] and we find, using the same energy estimate with (u, a, T 0 ) instead of (u ε , a ε , T ε ),

lim sup ε→0 Ω B(β(u ε (x, T ε )))dx ≤ Ω B(β(u(x, T 0 )))dx. ( 10
)
Combined with [START_REF] Ciarlet | The finite element method[END_REF], this shows that Ω B(β(u ε (x, T ε )))dx → Ω B(β(u(x, T 0 )))dx. The uniform convexity of B (see [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]) then allows us to deduce that ν(u

ε (•, T ε )) → ν(u(•, T 0 )) strongly in L 2 (Ω) and thus that ν(u ε ) → ν(u) strongly in C([0, T ]; L 2 (Ω)) (see Lemma 6.3).
Remark 1.2 A close examination of this proof indicates that equality in the energy estimate [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] is not required for u ε . An inequality ≤ would be sufficient. This is particularly important in the context of numerical methods which may introduce additional numerical diffusion (for example due to an implicit-in-time discretisation) and therefore only provide an upper bound in this energy estimate, see for example Estimate (39).

It is however crucial that the limit solution u satisfies the equivalent of (9) with an equal sign (or ≥).

2 Gradient discretisations and gradient schemes

Definitions

We give here a minimal presentation of gradient discretisations and gradient schemes, limiting ourselves to what is necessary to study the discretisation of (1). We refer the reader to [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF][START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF][START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF] for more details.

A gradient scheme can be viewed as a general formulation of several discretisations of (1) which are based on a nonconforming approximation of the weak formulation of the problem. The approximation of the weak formulation of ( 1) is based on some discrete spaces and mappings, the set of which we call a gradient discretisation. We say that D = (X D,0 , Π D , ∇ D , I D , (t (n) ) n=0,...,N ) is a space-time gradient discretisation for homogeneous Dirichlet boundary conditions if 1. the set of discrete unknowns X D,0 is a finite dimensional real vector space, 2. the linear mapping Π D : X D,0 → L ∞ (Ω) is a piecewise constant reconstruction operator in the sense that there exists a set I of degrees of freedom such that X D,0 = R I and there exists a family (Ω i ) i∈I of disjoint subsets of Ω such that Ω = i∈I Ω i and, for all u = (u i ) i∈I ∈ X D,0 and all i ∈ I,Π D u = u i on Ω i , 3. the linear mapping ∇ D : X D,0 → L p (Ω) d gives a reconstructed discrete gradient. It must be chosen such that

∇ D • L p (Ω) d is a norm on X D,0 , 4. I D : L 2 (Ω) → X D,0 is a linear interpolation operator, 5. t (0) = 0 < t (1) < t (2) < . . . < t (N ) = T .
We then set δt (n+ 1 2 ) = t (n+1) -t (n) , for n = 0, . . . , N -1, and δt D = max n=0,...,N -1 δt (n+ 1 2 ) , and we define the dual semi-norm |w| ,D of w ∈ X D,0 by

|w| ,D = sup Ω Π D w(x)Π D z(x)dx : z ∈ X D,0 , ||∇ D z|| L p (Ω) d = 1 . ( 11
)
Remark 2.2 (Boundary conditions) Other boundary conditions can be seamlessly handled by Gradient Schemes, see [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF].

Remark 2.3 (Nonlinear function of the elements of X D,0 ) Let D be a gradient discretisation in the sense of Definition 2.1. For any χ : R → R and any u = (u i ) i∈I ∈ X D,0 , we define χ I (u) ∈ X D,0 by χ I (u) = (χ I (u) i ) i∈I with χ I (u) i = χ(u i ). As indicated by the subscript I, this definition depends on the choice of the degrees of freedom in X D,0 . That said, these degrees of freedom are usually canonical and we therefore drop the index I. An important consequence of the fact that Π D is a piecewise constant reconstruction is the following:

∀χ : R → R , ∀u ∈ X D,0 , Π D χ(u) = χ(Π D u). (12)
It is customary to also use the notations Π D and ∇ D for space-time dependent functions. We will also need a notation for the jump-in-time of piecewise constant functions in time. Hence, if (v (n) ) n=0,...,N ⊂ X D,0 , we set for a.e.

x ∈ Ω, Π D v(x, 0) = Π D v (0) (x) and ∀n = 0, . . . , N -1 , ∀t ∈ (t (n) , t (n+1) ] : Π D v(x, t) = Π D v (n+1) (x) , ∇ D v(x, t) = ∇ D v (n+1) (x)
and

δ D v(t) = δ (n+ 1 2 ) D v := v (n+1) -v (n) δt (n+ 1 2 ) ∈ X D,0 . (13)
Thanks to Remark 2.3, the related gradient scheme is merely the discretisation of the weak formulation of (1) obtained by using the discrete space and mappings of the gradient discretisation.

If D = (X D,0 , Π D , ∇ D , I D , (t (n) ) n=0,...,N
) is a space-time gradient discretisation in the sense of Definition 2.1, we define the following gradient scheme for Problem (1): we consider a sequence (u (n) ) n=0,...,N ⊂ X D,0 such that

           u (0) = I D u ini and, for all v = (v (n) ) n=1,...,N ⊂ X D,0 , T 0 Ω [Π D δ D β(u)(x, t)Π D v(x, t) + a (x, Π D ν(u)(x, t), ∇ D ζ(u)(x, t)) • ∇ D v(x, t)] dxdt = T 0 Ω f (x, t)Π D v(x, t)dxdt. (14)
Remark 2. [START_REF] Bertsch | The Stefan problem with heating: appearance and disappearance of a mushy region[END_REF] We could as well consider, instead of a fully implicit method, a Crank-Nicolson scheme or any scheme in between those two. Such schemes are defined by taking α ∈ [ 1 2 , 1] and replacing the terms u (n+1) appearing in a(x, •, •) in (14) for t ∈ (t (n) , t (n+1) ] with u (n+α) = αu (n+1) + (1 -α)u (n) . All the results we establish for [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF] would hold for such a scheme (see the treatment in [START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF]).

Properties of gradient discretisations

In order to establish the convergence of the associated gradient schemes, sequences of spacetime gradient discretisations are required to satisfy four properties: coercivity, consistency, limitconformity and compactness.

∀ϕ ∈ L 2 (Ω) ∩ W 1,p 0 (Ω), S D (ϕ) = min w∈X D,0 ||Π D w -ϕ|| L max(p,2) (Ω) + ||∇ D w -∇ϕ|| L p (Ω) . (15)
A sequence (D m ) m∈N of space-time gradient discretisations in the sense of Definition 2.1 is said to be consistent if

• for all ϕ ∈ L 2 (Ω) ∩ W 1,p 0 (Ω), S Dm (ϕ) → 0 as m → ∞,
• for all ϕ ∈ L 2 (Ω), Π Dm I Dm ϕ → ϕ in L 2 (Ω) as m → ∞, and

• δt Dm → 0 as m → ∞.
Definition 2.7 (Compactness) If D is a space-time gradient discretisation in the sense of Definition 2.1, we define

∀ξ ∈ R d , T D (ξ) = max v∈X D,0 \{0} ||Π D v(• + ξ) -Π D v|| L p (R d ) ||∇ D v|| L p (Ω) ,
where Π D v has been extended by 0 outside Ω.

A sequence (D m ) m∈N of space-time gradient discretisations is said to be compact if

lim ξ→0 sup m∈N T D (ξ) = 0.
Definition 2.8 (Limit-conformity) If D is a space-time gradient discretisation in the sense of Definition 2.1, we define

∀ϕ ∈ W div,p (Ω) , W D (ϕ) = max u∈X D,0 \{0} Ω (∇ D u(x) • ϕ(x) + Π D u(x)divϕ(x)) dx ∇ D u L p (Ω) d . (16)
A sequence (D m ) m∈N of space-time gradient discretisations in the sense of Definition 2.1 is said to be limit-conforming if, for all ϕ ∈ W div,p (Ω), W Dm (ϕ) tends to 0 as m → ∞.

We refer the reader to [START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF][START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF] for a proof of the following lemma.

Lemma 2.9 (Consequence of limit-conformity) Let (D m ) m∈N be a sequence of space-time gradient discretisations in the sense of Definition 2.1, which is limit-conforming in the sense of Definition 2.8. Let, for any m ∈ N, v m = (v

(n)
m ) n=0,...,Nm ⊂ X Dm,0 be such that, with the notations in (13), (∇ Dm v m ) m∈N is bounded in L p (Ω × (0, T )). Then there exists v ∈ L p (0, T ; W 1,p 0 (Ω)) such that, up to a subsequence as m → ∞,

Π Dm v m → v weakly in L p (Ω × (0, T )) and ∇ Dm v m → ∇v weakly in L p (Ω × (0, T )) d .

Main results

Our first theorem states weak or space-time averaged convergence properties of Gradient Schemes for [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF]. These results are quite classical and have already been established for Richards' and Stefan's models, see [START_REF] Eymard | Gradient schemes for the Stefan problem[END_REF][START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF]. The convergence proof we provide afterwards however covers more non-linear model, as a do not need to be linear with respect to ξ, and is more compact than the ones available in the literature.

Theorem 2.10 (Convergence of the Gradient Scheme) Under Assumptions (2), let (D m ) mp∈N be a sequence of space-time gradient discretisations in the sense of Definition 2.1, which is coercive, consistent, compact and limit-conforming (see Section 2.2). Then for any m ∈ N there exists a solution u m to (14) with D = D m . Moreover, if we assume that

(∀s ∈ R , β(s) = s) or (∀s ∈ R , ζ(s) = s), (17)
then there exists a solution u to (4) such that, up to a subsequence, the following convergences hold as m → ∞:

Π Dm β(u m ) → β(u) weakly in L 2 (Ω) uniformly on [0, T ] (see Definition 6.4), Π Dm ν(u m ) → ν(u) strongly in L 1 (0, T ; L 1 (Ω)), Π Dm ζ(u m ) → ζ(u) weakly in L p (Ω × (0, T )). ∇ Dm ζ(u m ) → ∇ζ(u) weakly in L p (Ω × (0, T )) d . ( 18
) Remark 2.11 Since |ν| ≤ L ζ |β| and |ν| ≤ L β |ζ|, the L ∞ (0, T ; L 2 (Ω)) bound on β(u m ) and the L p (Ω × (0, T )) bound on ζ(u m ) (see Lemma 4.
1 and Definition 2.5) allows us to see that the strong convergence of Π Dm ν(u m ) is also valid in L q (0, T ; L r (Ω)) for any (q, r) such that q < ∞ and r < 2 or q, r < p (and, of course, any space interpolated between the two cases).

Remark 2.12 Note that we do not assume the existence of a solution u to the continuous problem, our convergence analysis will establish this existence.

Remark 2.13 Assumption (17) covers Richards' and Stefan's models, as well as many other nonlinear parabolic equations. This assumption is actually not mandatory if p ≥ 2, see Section 5.

We decide however to first state and prove Theorem 2.10 with this assumption to simplify the presentation. See also Remark 2.17.

The main innovation of this paper is the following theorem, which states the uniform-in-time convergence of numerical methods for fully non-linear parabolic equations with no regularity assumptions on the data.

Theorem 2.14 (Uniform-in-time convergence) Under Assumptions (2), let (D m ) mp∈N be a sequence of space-time gradient discretisations in the sense of Definition (2.1), which is coercive, consistent, compact and limit-conforming (see Section 2.2). We assume that, for any m ∈ N, u m is a solution to [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF] with D = D m , which converges to a solution u of (4) in the sense [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF].

Then, as m → ∞, Π Dm ν(u m ) → ν(u) strongly in L ∞ (0, T ; L 2 (Ω)).
Remark 2.15 Note that since (Π Dm ν(u m )) m∈N are piecewise constant in time, their convergence in L ∞ (0, T ; L 2 (Ω)) is actually a uniforme-in-time convergence (not "uniform a.e. in time").

The last theorem completes our convergence result by stating the strong space-time averaged convergence of the discrete gradients. Its proof is inspired by the study of Gradient Schemes for Leray-Lions operators made in [START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF].

Theorem 2.16 (Strong convergence of the gradients) Under Assumptions (2), let (D m ) mp∈N be a sequence of space-time gradient discretisations in the sense of Definition (2.1), which is coercive, consistent, compact and limit-conforming (see Section 2.2). We assume that, for any m ∈ N, u m is a solution to [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF] with D = D m , which converges to a solution u of (4) in the sense [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF]. We also assume that a is strictly monotone in the sense:

(a(x, s, ξ) -a(x, s, χ)) • (ξ -χ) > 0, for a.e. x ∈ Ω, ∀s ∈ R, ∀ξ = χ ∈ R d . ( 19
)
Then, as m → ∞, Π Dm ζ(u m ) → ζ(u) strongly in L p (Ω × (0, T )) and ∇ Dm ζ(u m ) → ∇ζ(u) strongly in L p (Ω × (0, T )) d .
Remark 2.17 Theorems 2.14 and 2.16 do not require the structural assumption (17); these theorems only require that the convergences (18) hold.

Preliminaries

We establish here a few results which will be used in the analysis of the Gradient Scheme (14).

Uniform-in-time compactness for space-time Gradient Discretisations

The first result is a consequence, in the framework of gradient discretisations, of the results in Section 6.

Theorem 3.1 (Uniform-in-time L 2 (Ω)-weak compactness result) Let T > 0, and (D m ) m∈N = (X Dm,0 , Π Dm , ∇ Dm , I Dm , (t (n)
m ) n=0,...,Nm ) m∈N be a sequence of space-time discretisation in the sense of Definition 2.1 which is consistent in the sense of Definition 2.6.

For any m ∈ N, let v m = (v (n) m ) n=0,...,Nm ⊂ X Dm,0 . If there exists q > 1 and C > 0 such that, for any m ∈ N, ||Π Dm v m || L ∞ (0,T ;L 2 (Ω)) ≤ C and T 0 |δ m v m (t)| q ,Dm dt ≤ C, (20)
then the sequence (Π Dm v m ) m∈N is relatively compact uniformly-in-time and weakly in L 2 (Ω), i.e. it has a subsequence which converges according to Definition 6.4. Moreover, any limit v of such a subsequence is continuous [0, T ] → L 2 (Ω) for the weak topology. Proof. This result is a consequence of the generalised Ascoli-Arzela theorem (Theorem 6.2) with K = [0, T ] and E the ball of radius C in L 2 (Ω), endowed with the weak topology. We let (ϕ l ) l∈N ⊂ C ∞ c (Ω) be a dense sequence in L 2 (Ω) and endow E with the metric (78) from these ϕ l , which indeed defines the weak L 2 (Ω) topology (see Proposition 6.5).

The set E is metric compact and therefore complete, and all Π Dm v m take their values in E. It remains to estimate d E (v m (s), v m (s )). We drop the index m in the spaces for legibility of notations. Let us define the interpolant P D ϕ l ∈ X D,0 by

P D ϕ l = argmin w∈X D,0 ||Π D w -ϕ l || L max(p,2) (Ω) + ||∇ D w -∇ϕ l || L p (Ω) . (21)
For any 0

≤ s ≤ s ≤ T , writing Π D v m (s ) -Π D v m (s) as the sum of its jumps δt (n+ 1 2 ) Π D δ (n+ 1
2 ) D v m at the points (t (n) ) n=n1,...,n2 lying between s and s , the definition of | • | ,D and Estimate [START_REF] Ern | Theory and practice of finite elements[END_REF] give

Ω (Π D v m (x, s ) -Π D v m (x, s)) Π D P D ϕ l (x)dx = t (n 2 +1) t (n 1 ) Ω Π D δ D v(t)(x)Π D P D ϕ l (x)dxdt ≤ C 1/q (t (n2+1) -t (n1) ) 1/q ||∇ D P D ϕ l || L p (Ω) d .
By definition of P D , we have

||Π D P D ϕ l -ϕ l || L 2 (Ω) ≤ S D (ϕ l )
and

||∇ D P D ϕ l || L p (Ω) d ≤ S D (ϕ l ) + ||ϕ l || L p (Ω) + ||∇ϕ l || L p (Ω) d ≤ C ϕ l
with C ϕ l not depending on D (and therefore on m). Since t (n2+1) -t (n1) ≤ |s -s| + δt and

(Π D v m ) m∈N is bounded in L ∞ (0, T ; L 2 (Ω)), we deduce that Ω (Π D v m (x, s ) -Π D v m (x, s)) ϕ l (x)dx ≤ 2C S D (ϕ l ) + C 1/q C ϕ l |s -s| 1/q + C 1/q C ϕ l δt 1/q . ( 22
)
Plugged into the definition (78) of the distance in E, this shows that

d E Π D v m (s ), Π D v m (s) ≤ l∈N min(1, C 1/q C ϕ l |s -s| 1/q ) 2 l + l∈N min(1, 2C S Dm (ϕ l ) + C 1/q C ϕ l δt 1/q m ) 2 l =: ω(s, s ) + δ m .
Using the dominated convergence theorem for series, we see that ω(s, s ) → 0 as s-s → 0 and that δ m → 0 as m → ∞ (we invoke the space consistency to establish that lim m→∞ S Dm (ϕ l ) → 0 for any l). Hence, the assumptions of Theorem 6.2 are satisfied and the proof is complete.

Technical results

We state here a family of technical lemmas, starting with a few useful properties on ν and B.

(ν(a) -ν(b)) 2 ≤ L β L ζ (ζ(a) -ζ(b))(β(a) -β(b)). ( (23)
) 24
The function B is convex on R β , the function

B • β : R → [0, ∞) is continuous, ∀s ∈ R , B(β(s)) = s 0 ζ(q)β (q)dq , (25)
∃K 0 , K 1 , K 2 > 0 such that ∀s ∈ R , K 0 β(s) 2 -K 1 ≤ B(β(s)) ≤ K 2 s 2 , ( 26
) ∀a ∈ R , ∀r ∈ R β , B(r) -B(β(a)) ≥ ζ(a)(r -β(a)), (27)
and

∀s, s ∈ R , (ν(s) -ν(s )) 2 ≤ 4L β L ζ B(β(s)) + B(β(s )) -2B β(s) + β(s ) 2 . ( 28
)
Proof. Inequality ( 23) is a straightforward consequence of the estimate ν = ζ β ≤ L β ζ . Note that the same inequality also holds with β and ζ swapped. Inequality ( 24) is then a direct application of ( 23) and the similar inequality with β and ζ swapped.

Let us first notice that, since β ≥ 0 on R + and β ≤ 0 on R -, β r (s) is always a real number when s ∈ R β . Moreover, since β is non-decreasing, β r is also non-decreasing on R β and therefore locally bounded on R β . Hence, B is well defined and locally Lipschitz-continuous, with an a.e. derivative B = ζ(β r ). B is therefore non-decreasing and B is convex.

To prove [START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF], we denote by P ⊂ R β the countable set of plateaux values of β, i.e. the y ∈ R such that β -1 ({y}) is not reduced to a singleton. If s ∈ β -1 (P ) then β -1 ({β(s)}) is the singleton {s} and therefore β r (β(s)) = s. Moreover, β r is continuous at β(s) and thus B is differentiable at β(s). Since β is differentiable a.e., we therefore deduce that, for a.e. s ∈ β -1 (P ), (B(β)

) (s) = B (β(s))β (s) = ζ(β r (β(s)))β (s) = ζ(s)β (s)
. The set β -1 (P ) is a union of intervals on which β and thus B(β) are locally constant; hence, for a.e. s in this set, (B(β)) (s) = 0 and ζ(s)β (s) = 0. Hence, the locally Lipschitz-continuous functions B(β) and s → s 0 ζ(q)β (q)dq have identical derivatives a.e. on R and take the same value at s = 0. They are thus equal on R and the proof of ( 25) is complete. The continuity of B • β is an obvious consequence of [START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF]. The second inequality in [START_REF] Eymard | The finite volume method for Richards equation[END_REF] can also be easily deduced from [START_REF] Eymard | Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation[END_REF] by noticing that |ζ(s)

β (s)| ≤ L ζ L β |s| (we can take K 2 = L β L ζ 2 ).
For the first inequality in [START_REF] Eymard | The finite volume method for Richards equation[END_REF], we first infer from (2b) the existence of S > 0 such that |ζ(q

)| ≥ M0 2 |q| ≥ M0 2L β |β(q)| whenever |q| ≥ S.
We then write, for s ≥ S,

B(β(s)) = S 0 ζ(q)β (q)dq + s S ζ(q)β (q)dq ≥ M 0 2L β s S β(q)β (q)dq = M 0 4L β β(s) 2 -β(S) 2 .
A similar inequality holds for s ≤ -S (with β(-S) instead of β(S)) and the first inequality in (26) therefore holds with

K 0 = M0 4L β and K 1 = M0 4L β max [-S,S] β 2 .
We now prove [START_REF] Eymard | Gradient scheme approximations for diffusion problems[END_REF], which states that a belongs to the convex sub-differential of B at β(a). We first start with the case

r ∈ R β , that is r = β(b) for some b ∈ R. If β r is continuous at β(a) then this inequality is an obvious consequence of the convexity of B since B is then differentiable at β(a) with B (β(a)) = ζ(β r (β(a))) = ζ(a).
Otherwise, a plain reasoning also does the job as

B(r) -B(β(a)) = B(β(b) -B(β(a)) b a ζ(q)β (q)dq = b a (ζ(q) -ζ(a))β (q)dq + ζ(a)(β(b) -β(a)) ≥ ζ(a)(r -β(a)),
the inequality coming from the fact that β ≥ 0 and ζ(q) -ζ(a) has the same sign as b -a when q is between a and b. The general case r ∈ R β is obtained by passing to the limit on b n such that β(b n ) → r and using the fact that B has limits (possibly +∞) at the endpoints of R β . Let us now take s, s ∈ R. Let s ∈ R be such that β(s) = β(s)+β(s ) 2 . We notice that

B(β(s)) + B(β(s )) -2B(β(s)) = s s (ζ(q) -ζ(s))β (q)dq + s s (ζ(q) -ζ(s))β (q)dq. We then use |ζ(q) -ζ(s)| ≥ 1 L β |ν(q) -ν(s)| and β (q) ≥ β (q) ζ (q) L ζ = ν (q) L ζ to write s s (ζ(q) -ζ(s))β (q)dq ≥ 1 L β L ζ s s ν (q)(ν(q) -ν(s))dq = 1 2L β L ζ (ν(s) -ν(s)) 2 . Thanks to (ν(s) -ν(s )) 2 ≤ 2((ν(s) -ν(s)) 2 + (ν(s ) -ν(s)) 2 ),
we deduce that (28) follows.

The next lemma is an easy consequence of Fatou's lemma and the fact that strongly lower semicontinuous convex functions are also weakly lower semi-continuous. We all the same provide its short proof. Then Hence, Φ is lower semi-continuous for the strong topology of L 2 (Ω; I). Since Φ (as H) is convex, we deduce that this lower semi-continuity property is also valid for the weak topology of L 2 (Ω; I), see [START_REF] Ekeland | Convex analysis and variational problems[END_REF]. The result of the lemma is just the translation of this weak lower semi-continuity of Φ.

Ω H(v(x))dx ≤ lim inf m→∞ Ω H(v m (x))dx. Proof. Let Φ : L 2 (Ω; I) → (-∞, ∞] be defined by Φ(w) = Ω H(w(x))dx. Since H is convex,
The last technical result is a consequence of the Minty trick. It has been proved and used in the L 2 case in [START_REF] Eymard | Gradient schemes for the Stefan problem[END_REF][START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF], but we need here an extension to the non-Hilbertian case.

Lemma 3.5 (Minty's trick) Let H ∈ C 0 (R) be a nondecreasing function. Let (X, µ) be a measurable set with finite measure and let

(u n ) n∈N ⊂ L p (X) with p > 1 such that 1. there exists u ∈ L p (X) such that (u n ) n∈N weakly converges to u in L p (X); 2. (H(u n )) n∈N ⊂ L 1 (X)
and there exists w ∈ L 1 (X) such that (H(u n )) n∈N strongly converges to w in L 1 (X);

Then w = H(u) a.e. on X.

Proof.

Let k > 0 and T k (s) = max(-k, min(s, k)) be the usual truncation at level k. Since H is non-decreasing, there exists

h k → ∞ as k → ∞ such that H(T k (s)) = T h k (H(s)). Thus, H(T k (u n )) → T h k (w) in L 1 (X) as n → ∞. Given that (H(T k (u n ))) n∈N remains bounded in L ∞ (X), its convergence to T h k (w) also holds in L p (X).
Using fact that H • T k is non-decreasing, we write, for any g ∈ L p (X),

X (H(T k (u n )) -H(T k (g)))(u n -g)dµ ≥ 0.
By strong convergence of H(T k (u n )) in L p (X) and weak convergence of u n in L p (X), as well as the fact that H • T k is bounded, we can take the limit of this expression as n → ∞ and we find

X (T h k (w) -H(T k (g)))(u -g)dµ ≥ 0. ( 29
)
We then use Minty's trick, that is pick a generic ϕ ∈ L p (X), apply (29) to g = u -tϕ, divide by t and let t → ±0 (using the dominated convergence theorem and the fact that H • T k is continuous and bounded) to find

X (T h k (w) -H(T k (u)))ϕdµ = 0. Selecting ϕ = sign(T h k (w) -H(T k (u))), we deduce that T h k (w) = H(T k (u))
a.e. on X. Letting k → ∞, we conclude that w = H(u) a.e. on X.

Integration-by-parts for the continuous solution

The last series of preliminary results are properties on the solution to (4), all based on the following integration-by-part property. This property, used in the proof of Theorem 2.10 and 2.14, allows to compute the value of the linear form ∂ t β(v) ∈ L p (0, T ; W -1,p (Ω)) on the function ζ(u). Because of the lack of regularity on u and the many non-linearities in [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF], justifying this integration-by-parts is however not straightforward at all... Lemma 3.6 Let us assume (2b) and (2c).

Let v : Ω × (0, T ) → R be measurable such that ζ(v) ∈ L p (0, T ; W 1,p 0 (Ω)), B(β(v)) ∈ L ∞ (0, T ; L 1 (Ω)), β(v) ∈ C([0, T ]; L 2 (Ω)-w) and ∂ t β(v) ∈ L p (0, T ; W -1,p (Ω)). Then t ∈ [0, T ] → Ω B(β(v)(x, t))dx ∈ [0, ∞) is continuous and, for all t 1 , t 2 ∈ [0, T ], t2 t1 ∂ t β(v)(t), ζ(v(t)) W -1,p ,W 1,p 0 dt = Ω B(β(v)(x, t 2 ))dx - Ω B(β(v)(x, t 1 ))dx. ( 30
)
Remark 3.7 Like at the end of Section 1.2, it is important to in mind the separation of

β(v(•, •)) from its continuous representative β(v)(•, •).
Proof.

Note that we obviously only need to make the proof when 0 ≤ t 1 < t 2 ≤ T .

Step 1: truncation, extensions and approximation of β(v).

We define β(v) : R → L 2 (Ω) by setting

β(v)(t) =    β(v)(t) if t ∈ [t 1 , t 2 ], β(v)(t 1 ) if t ≤ t 1 , β(v)(t 2 ) if t ≥ t 2 .
By the continuity property of β(v), this definition and gives β(v) ∈ C(R; L 2 (Ω)-w) such that ∂ t β(v) = 1 (t1,t2) ∂ t β(v) ∈ L p (R; W -1,p (Ω)) (no Dirac masses have been introduced at t = t 1 or t = t 2 ). This regularity of ∂ t β(v) ensures that the function

t ∈ R → D h β(v) := 1 h t+h t ∂ t β(v)(s)ds = β(v)(t + h) -β(v)(t) h in W -1,p (Ω) (31)
tend to ∂ t β(v) in L p (R; W -1,p (Ω)) as h → 0.

Step 2: we prove that

||B(β(v)(t)|| L 1 (Ω) ≤ ||B(β(v))|| L ∞ (0,T ;L 1 (Ω)) for all t ∈ R (not only for a.e. t). Let t ∈ [t 1 , t 2 ]. Since β(v)(•, •) = β(v(•, •)) a.e. on Ω × (t 1 , t 2 ), there exists a sequence t n → t such that β(v)(•, t n ) = β(v(•, t n )) in L 2 for all n and ||B(β(v)(•, t n ))|| L 1 (Ω) ≤ ||B(β(v))|| L ∞ (0,T ;L 1 (Ω)) . As β(v) ∈ C([0, T ]; L 2 (Ω)-w), we have β(v)(•, t n ) → β(v)(•, t) weakly in L 2 (Ω).
We then use the convexity of B and Lemma 3.4 to write, thanks to our choice of t n ,

Ω B(β(v)(x, t))dx ≤ lim inf n→∞ Ω B(β(v)(x, t n ))dx ≤ ||B(β(v))|| L ∞ (0,T ;L 1 (Ω))
and the proof is complete for t

∈ [t 1 , t 2 ]. The result for t ≤ t 1 or t ≥ t 2 is obvious since β(v)(t) is then either β(v)(t 1 ) or β(v)(t 2 ).
Step 3: We prove that for all τ ∈ R and a.e. t ∈ (t 1 , t 2 ),

β(v)(τ ) -β(v)(t), ζ(v(•, t)) W -1,p ,W 1,p 0 ≤ Ω B(β(v)(x, τ )) -B(β(v)(x, t))dx. (32)
Note that if we could just replace the duality product W -1,p -W 1,p 0 with an L 2 inner product, then this formula would be a straightforward consequence of [START_REF] Eymard | Gradient scheme approximations for diffusion problems[END_REF]. The problem is that nothing ensures that ζ(v)(t) ∈ L 2 (Ω) for a.e. t. We first notice that β(v)(τ ) -β(v)(t) = τ t ∂ t β(v)(s)ds indeed belongs to W -1,p (Ω) so the righthand side of (32) makes sense provided that t is chosen such that ζ(v(•, t)) ∈ W 1,p 0 (Ω) (which we do from here on). To deal with the fact that ζ(v(•, t)) does not necessarily belong to L 2 (Ω), we replace it with a truncation. Denoting by T k (s) = max(-k, min(s, k)) the classical truncation at level k, by the growth assumption (2b) on ζ we see that there exists [START_REF] Minty | On a monotonicity method for the solution of non-linear equations in Banach spaces[END_REF] the replacement of the duality product by an L 2 (Ω) inner product being justified since β(v)(τ )β(v)(t) and ζ(v(•, t)) both belong to L 2 (Ω). We also used the fact that, for a.e. t ∈ (t 1 , t 2 ), β(v)(•, t) = β(v(•, t)) a.e. on Ω, so [START_REF] Minty | On a monotonicity method for the solution of non-linear equations in Banach spaces[END_REF] is valid for a.e. t ∈ (t 1 , t 2 ). We then write

r k → +∞ as k → +∞ such that ζ(T k (v(•, t))) = T r k (ζ(v(•, t))). Hence, ζ(T k (v(•, t)) ∈ W 1,p 0 (Ω) and converges, as k → ∞, to ζ(v(•, t)) in W 1,p 0 (Ω). We can therefore write β(v)(τ ) -β(v)(t), ζ(v)(t) W -1,p ,W 1,p 0 = lim k→∞ β(v)(τ ) -β(v)(t), ζ(T k (v(•, t)) W -1,p ,W 1,p 0 = lim k→∞ Ω β(v)(x, τ ) -β(v(x, t)) ζ(T k (v(x, t))dx,
β(v(x, t)) = β(T k (v(x, t))) + (β(v(x, t)) -β(T k (v(x, t)))) and apply (27) with r = β(v)(x, τ ) and a = T k (β(v(x, t))) to find Ω β(v)(x, τ ) -β(v(x, t)) ζ(T k (v(x, t))dx = Ω β(v)(x, τ ) -β(T k (v(x, t))) ζ(T k (v(x, t))dx - Ω (β(v(x, t)) -β(T k (v(x, t)))) ζ(T k (v(x, t))dx ≤ Ω B(β(v)(x, τ )) -B(β(T k (v(x, t))))dx - Ω (β(v(x, t)) -β(T k (v(x, t)))) ζ(T k (v(x, t))dx.
Studying the cases v(x, t) ≥ k or v(x, t) ≤ -k, we notice that the last integrand is always nonnegative, so we actually can write

Ω β(v)(x, τ ) -β(v(x, t)) ζ(T k (v(x, t))dx ≤ Ω B(β(v)(x, τ )) -B(β(T k (v(x, t))))dx.
We then use the continuity of B • β and Fatou's lemma to deduce lim sup

k→∞ Ω β(v)(x, τ ) -β(v(x, t)) ζ(T k (v(x, t))dx ≤ Ω B(β(v)(x, τ ))dx -lim inf k→∞ Ω B(β(T k (v(x, t))))dx ≤ Ω B(β(v)(x, τ ))dx - Ω B(β(v(x, t)))dx
which, combined with [START_REF] Minty | On a monotonicity method for the solution of non-linear equations in Banach spaces[END_REF], concludes the proof of (32) (recall that t has been chosen such that β(v(•, t)) = β(v)(•, t) a.e. on Ω).

Step 4: proof of the formula

By convergence of D h β(v) to ∂ t β(v) in L p (0, T ; W -1,p (Ω)) and since 1 (t1,t2) ζ(v) ∈ L p (R; W 1,p 0 (Ω)), we have t2 t1 ∂ t β(v)(t), ζ(v(t)) W -1,p ,W 1,p 0 dt = R ∂ t β(v)(t), 1 (t1,t2) (t)ζ(v(•, t)) W -1,p ,W 1,p 0 dt = lim h→0 R D h β(v)(t), 1 (t1,t2) (t)ζ(v(•, t)) W -1,p ,W 1,p 0 dt = lim h→0 1 h t2 t1 β(v)(s + h) -β(v)(t), ζ(v(•, t) W -1,p ,W 1,p 0 dt. ( 34
)
We then use [START_REF] Maitre | Numerical analysis of nonlinear elliptic-parabolic equations[END_REF] for a.e. t ∈ (t 1 , t 2 ) to obtain, for h small enough such that

t 1 + h < t 2 , 1 h t2 t1 β(v)(t + h) -β(v)(t), ζ(v(•, t)) W -1,p ,W 1,p 0 dt ≤ 1 h t2 t1 Ω B(β(v)(x, t + h)) -B(β(v)(x, t))dxdt = 1 h t2+h t2 Ω B(β(v)(x, t))dxdt - 1 h t1+h t1 Ω B(β(v)(x, t))dxdt = Ω B(β(v)(x, t 2 ))dx - 1 h t1+h t1 Ω B(β(v)(x, t))dxdt
We now take the lim sup of this inequality, using the fact that

B(β(v)(•, t 2 )) is integrable (Step 2)
to take its integral out of the lim sup. Coming back to [START_REF] Nochetto | Approximation of degenerate parabolic problems using numerical integration[END_REF] we obtain

t2 t1 ∂ t β(v)(t), ζ(v(t)) W -1,p ,W 1,p 0 dt ≤ Ω B(β(v)(x, t 2 ))dx -lim inf h→0 1 h t1+h t1 Ω B(β(v)(x, t))dxdt. ( 35
) But since β(v) ∈ C([0, T ]; L 2 (Ω)-w), as h → 0 we have 1 h t1+h t1 β(v)(t)dt → β(v)(t 1
) weakly in L 2 (Ω). Hence, the convexity of B, Lemma 3.4 and Jensen's inequality give

Ω B(β(v)(x, t 1 ))dx ≤ lim inf h→0 Ω B 1 h t1+h t1 β(v)(x, t)dt dx ≤ lim inf h→0 Ω 1 h t1+h t1 B(β(v)(x, t))dtdx.
Plugged into [START_REF] Pop | Numerical schemes for degenerate parabolic problems[END_REF], this inequality shows that (30) holds with ≤ instead of =. The reverse inequality is obtained by reversing the time. We consider

v(t) = v(t 1 + t 2 -t). Then ζ( v), B(β( v)) and β( v)
have the same properties as ζ(v), B(β(v)) and β(v), and β( v) takes values β(v)(t 1 ) at t = t 2 and β(v)(t 2 ) at t = t 1 . Applying [START_REF] Hermeline | Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes[END_REF] with "≤" instead of "=" to v and using the fact that

∂ t β( v)(t) = -∂ t β(v)(t 1 + t 2 -t)
, we obtain [START_REF] Hermeline | Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes[END_REF] with "≥" instead of "=" and the proof of (30) is complete.

The continuity of t ∈ [0, T ] → Ω B(β(v)(x, t))dx is straightforward from (30) as the left-hand side of this relation is continuous with respect to t 1 and t 2 .

The following corollary states continuity properties and an essential formula on the solution to (4).

Corollary 3.8 Under Assumptions (2a)-(2i), if u is a solution of (4) then:

1. the function t ∈ [0, T ] → Ω B(β(u)(x, t))dx ∈ [0, ∞) is continuous and bounded, 2. for any T 0 ∈ [0, T ], Ω B(β(u)(x, T 0 ))dx + T0 0 Ω a(x, ν(u(x, t)), ∇ζ(u)(x, t)) • ∇ζ(u)(x, t)dxdt = Ω B(β(u ini (x)))dx + T0 0 Ω f (x, t)ζ(u)(x, t)dxdt, (36) 3. ν(u) is continuous [0, T ] → L 2 (Ω).
Remark 3.9 The continuity of ν(u) has to be understood in the same sense as the continuity of β(u), that is ν(u) is a.e. on Ω × (0, T ) equal to a continuous function [0, T ] → L 2 (Ω). We use in particular the same notation ν(u(•, •)) for the continuous representative ν(u)(•, •) as we did for the continuous representative of β(u).

Proof.

The continuity of t ∈ [0, T ] → Ω B(β(u)(x, t))dx ∈ [0, ∞) and Formula (36) are straightforward consequences of Lemma 3.6 with v = u and using ( 4) with v = ζ(u). Note that the bound on Ω B(β(u)(x, t))dx can be seen as a consequence of (36), or from Step 2 in the proof of Lemma 3.6. Let us prove the strong continuity of ν(u) : [0, T ] → L 2 (Ω). Let T be the set of τ ∈ [0, T ] such that β(u(•, τ )) = β(u)(•, τ ) a.e. on Ω, and let (s l ) l∈N and (t k ) k∈N be two sequences in T which converge to the same value s. Invoking [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] we can write

Ω (ν(u(x, s l )) -ν(u(x, t k ))) 2 dx ≤ 4L β L ζ Ω B(β(u)(x, s l ))dx + Ω B(β(u)(x, t k ))dx -8L β L ζ Ω B β(u)(x, s l ) + β(u)(x, t k ) 2 dx. (37) Since β(u)(•,s l )+β(u)(•,t k ) 2 → β(u)(•, s) weakly in L 2 (Ω) as l, k → ∞, Lemma 3.4 gives Ω B (β(u)(x, s)) dx ≤ lim inf l,k→∞ Ω B β(u)(x, s l ) + β(u)(x, t k ) 2 dx.
Taking the lim sup as l, k → ∞ of (37) and using the continuity of t → Ω B(β(u)(x, t))dx thus shows that

||ν(u(•, s l )) -ν(u(•, t k ))|| L 2 (Ω) → 0 as l, k → ∞. ( 38
)
The existence of an a.e. representative of ν(u(•, •)) which is continuous [0, T ] → L 2 (Ω) is a direct consequence of this convergence. Let s ∈ [0, T ] and (s l ) l∈N ⊂ T which converges to s. Applied with t k = s k , (38) shows that (ν(u(•, s l ))) l∈N is a Cauchy sequence in L 2 (Ω) and therefore that lim l→∞ ν(u(•, s l )) exists in L 2 (Ω). Relation (38) moreover also shows that this limit, that we can call ν(u)(•, s), does not depend on the Cauchy sequence in T which converges to s. With t k = s, we also see that whenever s ∈ T we have ν(u(•, s)) = ν(u)(•, s) a.e. on Ω, and ν(u)(•, •) is therefore equal to ν(u(•, •)) a.e. on Ω × (0, T ). It remains to establish that ν(u) thus defined is continuous [0, T ] → L 2 (Ω). For any (τ r ) r∈N ⊂ [0, T ] which converges to τ ∈ [0, T ], we can pick s r ∈ T ∩ (τ r -1 r , τ r + 1 r ) and t r ∈ T ∩ (τ

-1 r , τ + 1 r ) such that ||ν(u)(•, τ r ) -ν(u(•, s r ))|| L 2 (Ω) ≤ 1 r , ||ν(u)(•, τ ) -ν(u(•, t r ))|| L 2 (Ω) ≤ 1 r .
We therefore have

||ν(u)(•, τ r ) -ν(u)(•, τ )|| L 2 (Ω) ≤ 2 r + ||ν(u(•, s r )) -ν(u(•, t r ))|| L 2 (Ω) ,
which proves by (38) with l = k = r that ν(u)(•, τ r ) → ν(u)(•, τ ) in L 2 (Ω) as r → ∞ and completes the proof.

4 Proof of the convergence theorems

Estimates on the approximate solution

As it is usual in the study of numerical methods for PDE with strong non-linearities or without regularity assumptions on the data, everything starts with a priori estimates.

Lemma 4.1 (L ∞ (0, T ; L 2 (Ω)) estimate and discrete L p (0, T ; W 1,p 0 (Ω)) estimate) Under Assumptions (2), let D be a space-time gradient discretisation in the sense of Definition 2.1. Let u be a solution to Scheme [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF]. Then, for any T 0 ∈ (0, T ], denoting by k = 1, . . . , N the index such that T 0 ∈ (t (k-1) , t (k) ] we have

Ω B(Π D β(u)(x, T 0 ))dx + T0 0 Ω a(x, Π D ν(u)(x, t), ∇ D ζ(u)(x, t)) • ∇ D ζ(u)(x, t)dxdt ≤ Ω B(Π D β(I D u ini )(x))dx + t (k) 0 Ω f (x, t)Π D ζ(u)(x, t)dxdt. (39)
Consequently, there exists C 1 > 0, only depending on p, L β , C P ≥ C D (see Definition 2.5),

C ini ≥ Π D I D u ini L 2 (Ω)
, f and a such that

Π D B(β(u)) L ∞ (0,T ;L 1 (Ω)) ≤ C 1 , ∇ D ζ(u) L p (Ω×(0,T )) d ≤ C 1 and Π D β(u) L ∞ (0,T ;L 2 (Ω)) ≤ C 1 . ( 40
)
Proof. Using ( 12) and ( 27), we notice that for any n = 0, . . . , N -1, any t ∈ (t (n) , t (n+1) ],

Π D δ D β(u)(t)Π D ζ(u (n+1) ) = 1 δt (n+ 1 2 ) β(Π D u (n+1) ) -β(u (n) ) ζ(Π D u (n+1) ) ≥ 1 δt (n+ 1 2 ) B(Π D β(u (n+1) )) -B(Π D β(u (n) )) . Hence, taking v = (ζ(u (0) ), ζ(u (1) ), . . . , ζ(u (k) ), 0, . . . , 0) ⊂ X D,0 in (14), we find Ω B(Π D β(u)(x, t (k) ))dx + t (k) 0 Ω a(x, Π D ν(u)(x, t), ∇ D ζ(u)(x, t)) • ∇ D ζ(u)(x, t)dxdt ≤ Ω B(Π D β(u (0) )(x))dx + t (k) 0 f (x, t)Π D ζ(u)(x, t)dxdt (41)
Equation ( 39) is a straightforward consequence of this estimate, of the relation [START_REF] Droniou | A mixed finite volume scheme for anisotropic diffusion problems on any grid[END_REF]) and of the fact that the integrand involving a is nonnegative on [T 0 , t (k) ].

β(u)(•, T 0 ) = β(u)(•, t (k) ) (see
Using the Young inequality, we can write

t (k) 0 Ω f (x, t)Π D ζ(u)(x, t)dxdt ≤ 2 1/(p-1) C p D (pa) 1/(p-1) p f p L p (Ω×(0,t (k) )) + a 2C p D Π D ζ(u) p L p (Ω×(0,t (k) ))
and the first two estimates in (40) therefore follow from (41), the coercivity assumption (2f) on a and [START_REF] Eymard | The finite volume method for Richards equation[END_REF]. The estimate on Π

D β(u) = β(Π D u) in L ∞ (0, T ; L 2 (Ω)) is a consequence of the estimate on B(β(Π D u)) in L ∞ (0, T ; L 1 (Ω)
) and of [START_REF] Eymard | The finite volume method for Richards equation[END_REF].

Corollary 4.2 (Existence of a solution to the Gradient Scheme) Under Assumptions (2), if D is a gradient discretisation in the sense of Definition 2.1 then there exists at least a solution to the Gradient Scheme [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF].

Proof. For ρ ∈ [0, 1] we let β ρ (u) = ρu + (1 -ρ)β(u), ζ ρ (u) = ρu + (1 -ρ)ζ(u) and a ρ (x, s, ξ) = ρξ + (1 -ρ)a(x, s, ξ). It is clear that β ρ , ζ ρ
and a ρ satisfy the same assumptions as β, ζ and a for some L β , M 0 , M 1 and a not depending on ρ. We can therefore apply Lemma 4.1 to see that there exists C 2 not depending on s such that any solution u ρ to [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF]

with β = β ρ , ζ = ζ ρ and a = a ρ satisfies ||∇ D ζ ρ (u ρ )|| L p ((0,T )×Ω) d ≤ C 2 . Since ||∇ D • || L p (Ω×(0,T )) d is a norm on X D,0 , this shows that (ζ ρ (u ρ )) ρ∈[0,1]
remains bounded in this finite dimensional space, which implies in particular that for all i ∈ I,

(ζ ρ (u ρ ) i ) ρ∈[0,1]
. Using Assumption 2b for ζ ρ with constants not depending on ρ, we deduce that ((u ρ ) i ) ρ∈[0,1] remains bounded for any i ∈ I, and thus that (u ρ ) ρ∈[0,1] is bounded in X D,0 . But if ρ = 0 then ( 14) is a linear scheme. Any solution to this scheme being bounded in X D,0 , this shows that the underlying linear system is invertible. A topological degree argument [START_REF] Deimling | Nonlinear functional analysis[END_REF] then shows, combined with the uniform bound on (u ρ ) ρ∈[0,1] , that the scheme corresponding to ρ = 1, that is [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF], possesses at least one solution.

Lemma 4.3 (Estimate on the dual semi-norm of the discrete time derivative) Under Assumptions (2), let D be a space-time gradient discretisation in the sense of Definition 2.1. Let u be a solution to Scheme [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF]. Then there exists C 3 , only depending on p, L β ,

C P ≥ C D , C ini ≥ Π D I D u ini L 2 (Ω) , f , a, µ, a and T , such that T 0 |δ D β(u)(t)| p ,D dt ≤ C 3 . (42)
Proof. Let us take a generic v = (v (n) ) n=1,...,N ⊂∈ X D,0 as test function in Scheme [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF]. We have, thanks to Assumption (2h) on a,

T 0 Ω Π D δ D β(u)(x, t)Π D v(x, t)dxdt ≤ T 0 Ω (a(x) + µ|∇ D ζ(u)(x, t)| p-1 )|∇ D v(x, t)|dxdt + T 0 Ω f (x, t)Π D v(x, t)dxdt.
Using Hölder's inequality, Definition 2.5 and Estimates (40), this leads to the existence of C 4 > 0 only depending on p, L β , C P , C ini , f , a, a and µ such that

T 0 Ω Π D δ D β(u)(x, t)Π D v(x, t)dxdt ≤ C 4 ∇ D v L p (0,T ;L p (Ω)) d .
The proof of ( 42) is completed by selecting v = (|δ

(n+ 1 2 ) D β(u)| p -1
,D z (n) ) n=0,...,N with (z (n) ) n=0,...,N ⊂ X D,0 such that, for any n = 0, . . . , N -1, z (n+1) realises the supremum in [START_REF] Diaz | On a nonlinear parabolic problem arising in some models related to turbulent flows[END_REF] with w = δ (n+ 1 2 ) D β(u). Let u be a solution to Scheme [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF]. Then there exists C 5 , only depending on p, L β ,

L ζ , C P ≥ C D , C ini ≥ Π D I D u ini L 2 (Ω) , f , a, µ, a and T , such that Π D ν(u)(•, • + τ ) -Π D ν(u)(•, •) 2 L 2 (Ω×(0,T -τ )) ≤ C 5 (τ + δt), ∀τ ∈ (0, T ). ( 43
)
Proof. Let τ ∈ (0, T ). Thanks to [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF], we get that

Ω×(0,T -τ ) Π D ν(u)(x, t + τ ) -Π D ν(u)(x, t) 2 dxdt ≤ L β L ζ T -τ 0 A(t)dt, (44)
where, for almost every t ∈ (0, T -τ ),

A(t) = Ω Π D ζ(u)(x, t + τ ) -Π D ζ(u)(x, t) Π D β(u)(x, t + τ ) -Π D β(u)(x, t) dx. Let t ∈ (0, T -τ ).
Letting n 0 (t), n 1 (t) = 0, . . . , N -1 be such that t (n0(t)) ≤ t < t (n0(t)+1) and t (n1(t)) ≤ t + τ < t (n1(t)+1) , we may write

A(t) = Ω Π D ζ(u (n1(t)+1) )(x) -Π D ζ(u (n0(t)+1) )(x) n1(t) n=n0(t)+1 δt (n+ 1 2 ) Π D δ (n+ 1 2 ) D β(u)(x) dx,
which also reads [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF]. Using (45), we get

A(t) = Ω Π D ζ(u (n1(t)+1) )(x) -Π D ζ(u (n0(t)+1) )(x) × N -1 n=1 χ n (t, t + τ )δt (n+ 1 2 ) Π D δ (n+ 1 2 ) D β(u)(x) dx, (45) with χ n (t, t + τ ) = 1 if t (n) ∈ (t, t + τ ] and χ n (t, t + τ ) = 0 if t (n) / ∈ (t, t + τ ]. We then let v = (χ n-1 (t, t + τ )(ζ(u (n1(t)+1) ) -ζ(u (n0(t)+1) ))) n=0,...,N in Scheme
A(t) = A 3 (t) - N -1 n=1 χ n (t, t + τ )δt (n+ 1 2 ) Ω a x, Π D ν(u (n+1) ), ∇ D ζ(u (n+1) )(x) • ∇ D ζ(u (n1(t)+1) )(x) -∇ D ζ(u (n0(t)+1) )(x) dx with A 3 (t) = N -1 n=1 χ n (t, t + τ ) Ω t (n+1) t (n) f (x, t)dt Π D ζ(u (n1(t)+1) )(x) -Π D ζ(u (n0(t)+1) )(x) dx.
Using the Young inequality and the inequality a 1 a 2 ≤ (p-1) 1/p p (a p 1 + a p 2 ) for a 1 ≥ 0 and a 2 ≥ 0, this yields:

A(t) ≤ (p -1) 1/p p (A 0 (t) + A 1 (t) + 2A 2 (t)) + A 3 (t), (46)
with

A 0 (t) = N -1 n=1 χ n (t, t + τ )δt (n+ 1 2 ) Ω ∇ D ζ(u (n0(t)+1) )(x) p dx, A 1 (t) = N -1 n=1 χ n (t, t + τ )δt (n+ 1 2 ) Ω ∇ D ζ(u (n1(t)+1) )(x) p dx,
and

A 2 (t) = N -1 n=1 χ n (t, t + τ )δt (n+ 1 2 ) Ω a x, Π D ν(u (n+1) ), ∇ D ζ(u (n+1) )(x) p dx.
Applying (80) of Lemma 6.6 yields

T -τ 0 A 0 (t) + A 1 (t)dt ≤ 2(τ + δt) ∇ D ζ(u) p L p (Ω×(0,T )) . (47)
Applying (79) of Lemma 6.6 yields

T -τ 0 A 2 (t)dt ≤ τ a (•, Π D ν(u), ∇ D ζ(u)) p L p (Ω×(0,T )) , (48)
and, similarly, we obtain using Definition (2.5) and (40), the existence of C 6 such that

T -τ 0 A 3 (t)dt ≤ C 6 ((τ + δt) + τ f p L p (Ω×(0,T )) ). (49)
Using inequalities (44), ( 46) and ( 47)-( 49), (43) is proved.

Proof of Theorem 2.10

Step 1 Application of compactness results. Thanks to Theorem 3.1 and Estimates (40) and (42), we first extract a subsequence such that (Π Dm β(u m )) m∈N converges weakly in L 2 (Ω) uniformly in [0, T ] (in the sense of Definition 6.4) to some function

β ∈ C([0, T ]; L 2 (Ω)-w) which satisfies β(•, 0) = β(u ini ) in L 2 (Ω).
Using again Estimates (40) and applying Lemma 2.9, we again extract a subsequence such that, for some

ζ ∈ L p (0, T ; W 1,p 0 (Ω)), Π Dm ζ(u m ) → ζ weakly in L p (Ω × (0, T )) and ∇ Dm ζ(u m ) → ∇ζ weakly in L p (Ω × (0, T )) d .
From Estimates (40), Definition 2.5 and the growth assumption (2b) on ζ, we also see that (Π Dm u m ) m∈N is bounded in L p (Ω × (0, T )) and we can therefore assume, up to a subsequence, that it converges weakly to some u in this space. We then remark that |ν

||Π Dm ν(u m )(• + ξ) -Π Dm ν(u)|| L p (R d ×(0,T )) ≤ L β C 1 T Dm (ξ) (50)
where Π D ν(u) has been extended by 0 outside Ω and lim ξ→0 sup m∈N T Dm (ξ) = 0. Invoking Lemma 4.4, we also see that the time translates of Π Dm ν(u m ) tend uniformly to 0 in L 1 (Ω × (0, T )). We use the fact that Π Dm β(u m ), and therefore also Π Dm ν(u m ), remains bounded in L ∞ (0, T ; L 2 (Ω)) to control the time translates at both ends of [0, T ]. Hence, applying Kolmogorov's theorem, we deduce that, up to the extraction of another subsequence, Π Dm ν(u m ) → ν in L 1 (Ω × (0, T )).

Under the first case in the structural hypothesis ( 17 [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF]. We notice that this is the only place where we need this structural assumption (17) on β, ζ.

Upon extraction of another subsequence, we can also assume that a (•, Π D ν(u), ∇ D ζ(u)) has a weak limit in L p (Ω × (0, T )) d , which we denote by A.

Finally, for any

T 0 ∈ [0, T ], since Π Dm β(u m (•, T 0 )) → β(u)(•, T 0 ) weakly in L 2 (Ω), Lemma 3.4 gives Ω B(β(u)(x, T 0 ))dx ≤ lim inf m→∞ Ω B(β(Π Dm u m )(x, T 0 ))dx. (51)
This shows in particular that B(β(u)) ∈ L ∞ (0, T ; L 1 (Ω)).

Step 2 Passing to the limit in the scheme. We drop the indices m for legibility reasons. Let ϕ ∈ C 1 c (-∞, T ) and w ∈ W 1,p 0 (Ω) ∩ L 2 (Ω). We introduce v = (ϕ(t (n-1) )P D w) n=0,...,N (where t (-1) = 0 for example, this value being irrelevant) as test function in [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF], with P D defined by [START_REF] Eymard | Gradient schemes for the Stefan problem[END_REF]. We get

T (m) 1 + T (m) 2 = T (m) 3 with T (m) 1 = N -1 n=0 ϕ(t (n) )δt (n+ 1 2 ) Ω Π D δ (n+ 1 2 ) D β(u)(x)Π D P D w(x)dx, T (m) 2 = N -1 n=0 ϕ(t (n) )δt (n+ 1 2 ) Ω a x, Π D ν(u (n+1) ), ∇ D ζ(u (n+1) )(x) • ∇ D P D w(x)dx, and
T (m) 3 = N -1 n=0 ϕ(t (n) ) t (n+1) t (n) Ω f (x, t)Π D P D w(x)dxdt.
Using discrete integrate-by-parts to transform the terms ϕ(t

(n) )(Π D β(u (n+1) ) -Π D β(u (n) )) ap- pearing in T (m) 1 into (ϕ(t (n) ) -ϕ(t (n+1) ))Π D β(u (n+1) ), we have T (m) 1 = - T 0 ϕ (t) Ω Π D β(u)(x, t)Π D P D w(x)dxdt -ϕ(0) Ω Π D β(u (0) )(x)Π D P D w(x)dx. Setting ϕ D (t) = ϕ(t (n) ) for t ∈ (t (n) , t (n+1) ), we have T (m) 2 = T 0 ϕ D (t) Ω a (x, Π D ν(u)(x, t), ∇ D ζ(u)(x, t)) • ∇ D P D w(x)dxdt T (m) 3 = T 0 ϕ D (t) Ω f (x, t)Π D P D w(x)dxdt.
We may thus let m → ∞ in

T (m) 1 + T (m) 2 = T (m) 3 to see that u satisfies                  u ∈ L p (Ω × (0, T )) , ζ(u) ∈ L p (0, T ; W 1,p 0 (Ω)) , B(β(u)) ∈ L ∞ (0, T ; L 1 (Ω)), - T 0 ϕ (t) Ω β(u(x, t))w(x)dxdt -ϕ(0) Ω u ini (x)w(x)dx + T 0 ϕ(t) Ω A(x, t) • ∇w(x)dxdt = T 0 ϕ(t) Ω f (x, t)w(x)dxdt, ∀w ∈ W 1,p 0 (Ω) ∩ L 2 (Ω), ∀ϕ ∈ C ∞ c (-∞, T ). (52)
Linear combinations of this relation show that it also holds with ϕ(t)w(x) replaced by a tensorial functions in C ∞ c (Ω × (0, T )). This proves that ∂ t β(u) ∈ L p (0, T ; W -1,p (Ω)). Standard arguments then show that β(u) can be identified with an element of C 0 ([0, T ]; L 2 (Ω)-w) with the property β(u(•, 0)) = β(u ini ) (cf. discussion at the end of Section 1.2). Using the density of tensorial functions in L p (0, T ; W 1,p 0 (Ω)), we then see that u satisfies

T 0 ∂ t β(u)(•, t), v(•, t) (W -1,p (Ω)) W 1,p 0 (Ω) dt + T 0 Ω A(x, t) • ∇v(x, t)dxdt = T 0 Ω f (x, t)v(x, t)dxdt , ∀v ∈ L p (0, T ; W 1,p 0 (Ω)). ( 53
)
Step 3 Proof that u is solution to (4).

The proof will be completed by showing that

A(x, t) = a(x, ν(u)(x, t), ∇ζ(u)(x, t)), for a.e. (x, t) ∈ Ω × (0, T ). ( 54
)
We take T 0 ∈ [0, T ], write (39) with D = D m and take the supremum limit as m → ∞. We notice that the t (k) = T m from Lemma 4.1 converges to T 0 as m → ∞. Hence, using the quadratic growth of B we obtain

lim sup m→∞ T0 0 Ω a(x, Π Dm ν(u m )(x, t), ∇ Dm ζ(u m )(x, t)) • ∇ Dm ζ(u m )(x, t)dxdt ≤ Ω B(β(u ini )(x))dx + T0 0 Ω f (x, t)ζ(u)(x, t)dxdt -lim inf m→∞ Ω B(β(Π Dm u m )(x, T 0 ))dx. ( 55
)
We then apply Lemma 3.6 and take v = ζ(u) in (53) to get

Ω B(β(u)(x, T 0 ))dx - Ω B(β(u)(x, 0))dx + T0 0 Ω A(x, t) • ∇ζ(u)(x, t)dxdt = T0 0 Ω f (x, t)ζ(u)(x, t)dxdt.
This relation, combined with (55) and using (51), shows that lim sup

m→∞ T0 0 Ω a(x, Π Dm ν(u m )(x, t), ∇ Dm ζ(u m )(x, t)) • ∇ Dm ζ(u m )(x, t)dxdt ≤ T0 0 Ω A(x, t) • ∇ζ(u)x, t)dxdt. ( 56
)
It is now possible to apply Minty's trick. Consider, for G ∈ L p (0, T ;

L p (Ω)) d , T0 0 Ω [a(•, Π D ν(u), ∇ D ζ(u)) -a(•, Π D ν(u), G)] • [∇ D ζ(u) -G] dxdt ≥ 0. ( 57
)
By strong convergence of Π Dm ν(u m ) to ν(u) in L 1 (Ω × (0, T )) and Assumptions (2e), (2h) on a, we see that a(•, Π Dm ν(u m ), G) → a(•, ν(u), G) strongly in L p (Ω × (0, T )). We can therefore develop (57) (with T 0 = T ) and use (56) to pass to the supremum limit on the only "weak-weak" term, and we find, for any G ∈ L p (0, T ;

L p (Ω)) d , T 0 Ω [A(x, t) -a(x, ν(u)(x, t), G(x, t))] • [∇ζ(u)(x, t) -G(x, t)] dxdt ≥ 0.
Application of Minty's method [START_REF] Minty | On a monotonicity method for the solution of non-linear equations in Banach spaces[END_REF] (i.e. taking G = ∇ζ(u) + rϕ for ϕ ∈ L p (0, T ; L p (Ω)) d and letting r → 0) then shows that (54) holds and concludes the proof that u is a weak solution to (4).

Proof of Theorem 2.14

Let T 0 ∈ [0, T ] and (T m ) m≥1 be a sequence in [0, T ] which converges to T 0 . By setting T 0 = T m and G = ∇ζ(u) in the developed form of (57), by taking the infimum limit (thanks to the strong convergence of a(•, Π Dm ν(u m ), ∇ζ(u))) and by using (54), we find

lim inf m→∞ Tm 0 Ω a(x, Π Dm ν(u m )(x, t), ∇ Dm ζ(u m )(x, t)) • ∇ Dm ζ(u m )(x, t)dxdt ≥ T0 0 Ω a(x, ν(u)(x, t), ∇ζ(u)(x, t)) • ∇ζ(u)(x, t)dxdt. ( 58
)
We then write (39) with T m instead of T 0 . We notice that the t (k) = t (k(m)) such that T m ∈ (t (k-1) , t (k) ] converges to T 0 as m → ∞. Using (58) and Corollary 3.8, we therefore obtain lim sup

m→∞ Ω B(β(Π Dm u m (x, T m )))dx ≤ Ω B(β(u)(x, T 0 ))dx. (59)
ν(Π Dm u m (•, T m )) -ν(u)(•, T 0 ) 2 L 2 (Ω) ≤ 2 ν(Π Dm u m (•, T m )) -ν(u(•, s m )) 2 L 2 (Ω) + 2 ν(u(•, s m )) -ν(u)(•, T 0 )) 2 L 2 (Ω) ≤ 8L β L ζ Ω [B(β(Π Dm u m (x, T m ))) + B(β(u(x, s m )))] dx -16L β L ζ Ω B β(Π Dm u m (x, T m )) + β(u(x, s m )) 2 dx + 2 ν(u(•, s m )) -ν(u)(•, T 0 )) 2 L 2 (Ω) .
We then take the lim sup as m → ∞ of this expression. Thanks to (59) and the continuity of Writing (39) for u m with T 0 = T , taking the lim sup as m → ∞, using (61) (with T m = T ) and the continuous integration-by-part formula (36) we find that lim sup

t ∈ [0, T ] → Ω B(β(u)(x, t))dx ∈ [0, ∞) (see
m→∞ T 0 Ω a(x, Π Dm ν(u m )(x, t), ∇ Dm ζ(u m )(x, t)) • ∇ Dm ζ(u m )(x, t)dxdt ≤ T0 0 Ω a(x, ν(u)(x, t), ∇ζ(u)(x, t)) • ∇ζ(u)(x, t)dxdt.
Combined with (58), this shows that

lim m→∞ T 0 Ω a(x, Π Dm ν(u m )(x, t), ∇ Dm ζ(u m )(x, t)) • ∇ Dm ζ(u m )(x, t)dxdt = T0 0 Ω a(x, ν(u)(x, t), ∇ζ(u)(x, t)) • ∇ζ(u)(x, t)dxdt. ( 62
)
Let us define

f m = [a(x, Π Dm ν(u m ), ∇ Dm ζ(u m )) -a(x, Π Dm ν(u m )(•, t), ∇ζ(u))] • [∇ Dm ζ(u m ) -∇ζ(u)] ≥ 0.
Developing this expression and using (62) and ( 18), we see that

T 0 Ω f m (x, t)dxdt → 0 as m → ∞.
This shows that f m → 0 in L 1 (Ω × (0, T )) and therefore a.e. up to a subsequence. We can then reason as in [START_REF] Droniou | Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF], using the strict monotony [START_REF] Ekeland | Convex analysis and variational problems[END_REF] of a, the coercivity assumption (2f) and Vitali's theorem, to deduce that

∇ Dm ζ(u m ) → ∇ζ(u) strongly in L p (Ω × (0, T )) d as m → ∞.
5 About Assumption [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF] We discuss show here that, if p ≥ 2, Theorem 2.10 can be proved without the structural assumption (17) -i.e. without assuming that β = Id or ζ = Id. We first notice that we can always assume that

β + ζ is strictly increasing. ( 63
)
Remark 5.3 Note that it is not proved, in this context, that u is a weak limit of Π Dm u m , but such a limit is not stated in [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] and can actually be considered as irrelevant for the model (1), in which the quantities of interest (physically relevant when this PDE models a natural phenomenon) are β(u), ζ(u) and ν(u).

We now state the two key results that allowed us to replace Assumption [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF] by Assumption [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF]. The first is a discrete version of a compensated compactness result in [START_REF] Kazhikhov | Recent developments in the global theory of two-dimensional compressible Navier-Stokes equations[END_REF]. The second is a Minty-like result, useful to identify weak non-linear limits. We note that Theorem 5.4 is a more general version than the one needed in the proof of Theorem 5.2 (which only requires ϕ ≡ 1), but we state this more general version nevertheless as it is the genuine discrete equivalent of the compensated compactness result in [START_REF] Kazhikhov | Recent developments in the global theory of two-dimensional compressible Navier-Stokes equations[END_REF].

Theorem 5.4 (Discrete compensated compactness) We take T > 0, p ≥ 2 and a sequence

(D m ) m∈N = (X Dm,0 , Π Dm , ∇ Dm , I Dm , (t (n)
m ) n=0,...,Nm ) m∈N of space-time discretisation in the sense of Definition 2.1 which is consistent, limit-conforming and compact in the sense of Definitions 2.6, 2.8 and 2.7. For any m ∈ N, let β m = (β

(n) m ) n=0,...,Nm ⊂ X Dm,0 and ζ m = (ζ (n) m ) n=0,...,Nm ⊂ X Dm,0 such that • The sequences ( T 0 |δ m β m (t)| ,Dm ) m∈N and (||∇ Dm ζ m || L 2 (0,T ;L p (Ω) d ) ) m∈N are bounded, • As m → ∞, Π Dm β m → β and Π Dm ζ m → ζ weakly in L 2 (Ω × (0, T )). Then Π Dm β m Π Dm ζ m → β ζ in the sense of measures on Ω×(0, T ), that is, for all ϕ ∈ C(Ω×[0, T ]), lim m→∞ T 0 Ω Π Dm β m (x, t)Π Dm ζ m (x, t)ϕ(x, t)dxdt = T 0 Ω β(x, t) ζ(x, t)ϕ(x, t)dxdt. (64)
Proof.

The idea of the proof is to reduce to the case where Π Dm ζ m are tensorial functions, in order to separate the space and time variables and make use of the separate compactness of Π Dm ζ m and Π Dm β m with respect to each of these variables. Note that this proof also provides an apparently new proof for the continuous equivalent of this compensated compactness result.

Step 1: reduction of Π Dm ζ m to tensorial functions. Let us take δ > 0 and consider a covering (A δ k ) k=1,...,K of Ω in disjoint cubes of length δ. Let us denote R δ : L 2 (Ω) → L 2 (Ω) the operator defined by: ∀k

= 1, . . . , K , ∀x ∈ A δ k : R δ f (x) = 1 meas(A δ k ) A k δ f (y)dy,
where f has been extended by 0 outside Ω. Let x ∈ A δ k . Using Jensen's inequality, fact that meas(A δ k ) = δ d and the change of variable y

∈ A δ k → ξ = y -x ∈ (-δ, δ) d , we can write |R δ f (x) -f (x)| 2 ≤ δ -d A δ k |f (y) -f (x)| 2 dy ≤ δ -d (-δ,δ) d |f (x + ξ) -f (x)| 2 dξ.
Integrating over x ∈ A δ k and summing over k = 1, . . . , K then gives

||R δ f -f || 2 L 2 (Ω) ≤ δ -d (-δ,δ) d Ω |f (x + ξ) -f (x)| 2 dxdξ ≤ 2 d sup ξ∈(-δ,δ) d Ω |f (x + ξ) -f (x)| 2 dx. ( 65
)
On the other side, the compactness of (D m ) m∈N (Definition 2.7) and the fact that p ≥ 2 give (ξ) such that (ξ) → 0 as ξ → 0 and, for all m ∈ N and all v ∈ X Dm,0 ,

||Π Dm v(• + ξ) -Π Dm v|| 2 L 2 (Ω) ≤ (ξ)||∇ Dm v|| 2 L p (Ω) d .
Combining this with (65) and using the bound on

||∇ Dm ζ m || L 2 (0,T ;L p (Ω) d ) shows that ||R δ Π Dm ζ m -Π Dm ζ m || L 2 (Ω×(0,T )) ≤ C sup |ξ|≤δ (ξ) =: ω(δ) (66
|A m (Π Dm ζ m ) -A(ζ)| ≤ Cω(δ) + |A m (R δ Π Dm ζ m ) -A(R δ ζ)|.
Let us assume that we can prove that, for a fixed δ,

A m (R δ Π Dm ζ m ) → A(R δ ζ) as m → ∞. (67)
Then the previous inequality gives lim sup

m→∞ |A m (Π Dm ζ m ) -A(ζ)| ≤ Cω(δ). Letting δ → 0 in this inequality gives A m (Π Dm ζ m ) → A(ζ) as wanted.
Hence, we only need to prove 67. The definition of R δ shows that R

δ f = K k=1 1 A δ k meas(A δ k ) -1 [f ] A δ k where 1 A δ k is the characteristic function of A δ k and [f ] A = A f (x)dx.
Hence, (67) will follow if we can prove that, for any measurable set A of non-zero measure,

lim m→∞ T 0 Ω Π Dm β m (x, t)[Π Dm ζ m ] A (t)ϕ(t, x)1 A (x)dxdt = T 0 Ω β(x, t)[ ζ ] A (t)ϕ(t, x)1 A (x)dxdt (68)
where, for g ∈ L 2 (Ω × (0, T )), we set [g] A (t) = A g(t, y)dy.

Step 2: further reductions. We now reduce ϕ to a tensorial function and 1 A to a smooth function. It is well-known that there exists tensorial functions ϕ r = Lr l=1 θ l,r (t)γ l,t (x), with θ l,r ∈ C ∞ ([0, T ]) and γ l,r ∈ C ∞ (Ω), such that ϕ r → ϕ uniformly on Ω × (0, T ) as r → ∞. Moreover, there exists

ρ r ∈ C ∞ c (Ω) such that ρ r → 1 A in L 2 (Ω) as r → ∞. Hence, as r → ∞ the function (t, x) → ϕ r (t, x)ρ r (x) converges in L ∞ (0, T ; L 2 (Ω)) to the func- tion (t, x) → ϕ(t, x)1 A (x). Since the sequence of functions (t, x) → Π Dm β(t, x)[Π Dm ζ m ] A (t) is bounded in L 1 (0, T ; L 2 (Ω)) (notice that ([Π Dm ζ m ] A ) m∈N is bounded in L 2 (0, T ) since (Π Dm ζ m ) m∈N is bounded in L 2 (Ω × (0, T ))
), a reasoning similar as the one used in Step 1 shows that we only need to prove (68) with ϕ(t, x)1 A (x) replaced with ϕ r (t, x)ρ r (x) for a fixed r. We have ϕ r (t, x)ρ r (x) = Lr l=1 θ l,r (t)(γ l,t ρ)(x) and γ l,r ρ r ∈ C ∞ c (Ω). Hence, (68) with ϕ(t, x)1 A (x) replaced with ϕ r (t, x)ρ r (x) will follow if we can establish that: for any θ ∈ C ∞ ([0, T ]), any ψ ∈ C ∞ c (Ω) and any measurable set A ⊂ Ω,

lim m→∞ T 0 Ω θ(t)Π Dm β m (x, t)[Π Dm ζ m ] A (t)ψ(x)dxdt = T 0 Ω θ(t)β(x, t)[ ζ ] A (t)ψ(x)dxdt. ( 69
)
Step 3: proof of (69).

We now use the estimates on δ m β m to conclude. We first write

T 0 Ω θ(t)Π Dm β m (x, t)[Π Dm ζ m ] A (t)ψ(x)dxdt = T 0 θ(t)[Π Dm ζ m ] A (t)F m (t). ( 70
) with F m (t) = Ω Π Dm β m (x, t)ψ(x)dx. It is clear from the weak convergence of Π Dm ζ m that [Π Dm ζ m ] A → [ ζ ] A weakly in L 2 (0, T ). Hence, if we can prove that F m → F := Ω β(x, •)ψ(x)dx
strongly in L 2 (0, T ), we can pass to the limit in (70) and obtain (69). Since F m weakly converges to F in L 2 (0, T ) (thanks to the weak convergence in L 2 (Ω × (0, T )) of Π Dm β m ), we only have to prove that (F m ) m∈N is relatively compact in L 2 (0, T ).

To prove the strong convergence of F m , we introduce the interpolant P Dm defined by ( 21) and we define G m as F m with ψ replaced with Π Dm P Dm ψ. We then have

|F m (t) -G m (t)| ≤ ||Π Dm β m (•, t)|| L 2 (Ω) S Dm (ψ).
The consistency of (D m ) then shows that

F m -G m → 0 strongly in L 2 (Ω) as m → ∞. (71)
We now study the strong convergence of G m . This function is, as Π Dm β m , piecewise constant on (0, T ) and, by definition of | • | ,Dm , its discrete derivative satisfies

|δ m G m (t)| ≤ |δ m β m (t)| ,Dm ||P Dm ψ|| X Dm,0 .
Since ||P Dm ψ|| X Dm,0 ≤ S Dm (ψ) + ||∇ψ|| L p (Ω) d is bounded uniformly with respect to m, the assumption on δ m β m proves that (||δ m G m || L 1 (0,T ) ) m∈N is bounded. Hence, G m is bounded in BV (0, T ) ∩ L 1 (0, T ) and therefore relatively compact in L 2 (0, T ) (see [START_REF] Attouch | Variational analysis in Sobolev and BV spaces, volume 6 of MPS/SIAM Series on Optimization[END_REF]Theorem 10.1.4]). Combined with (71), this proves that (F m ) m∈N is relatively compact in L 2 (0, T ) and concludes the proof.

Remark 5.5 If we assume that (Π Dm β m ) m∈N is bounded in L ∞ (0, T ; L 2 (Ω)) and that, for some q > 1, ( T 0 |δ m β m (t)| q ,Dm ) m∈N is bounded, then Step 2 becomes a trivial consequence of Theorem 3.1. Indeed, this theorem then shows that (Π Dm β m ) m∈N is relatively compact uniformly-in-time and weakly in L 2 (Ω), which translates into the relative compactness of (F m ) m∈N in L ∞ (0, T ). Lemma 5.6 Let N ∈ N and V be a non-empty measurable subset of R N . Let β, ζ ∈ C 0 (R) be two nondecreasing functions, such that β(0) = ζ(0) = 0 and β + ζ is strictly increasing. We assume that there exists a sequence (w m ) m∈N of measurable functions on V and two functions β, ζ ∈ L 2 (V ) such that:

• β(w m ) → β and ζ(w m ) → ζ weakly in L 2 (V ),

• there exists ϕ ∈ L ∞ (V ) such that ϕ > 0 a.e. on V and Proof. Let us first notice that the last conclusion of the theorem, i.e. that any adherence value v of (v m ) m∈N in F(K, E) is continuous, is trivially obtained by passing to the limit along this subsequence in (76), showing that the modulus of continuity of v is bounded above by ω.

The proof of the compactness result an easy generalisation of the proof of the classical Ascoli-Arzela compactness result. We start by taking a countable dense subset {s l : l ∈ N} in K (the existence of this set is ensured since K is compact metric). Since each set {v m (s l ) : m ∈ N} is relatively compact in E, by diagonal extraction we can select a subsequence of (v m ) m∈N , denoted the same way, such that for any l ∈ N, (v m (s l )) m∈N converges in E. We then proceed in showing that (v m ) m∈N is a Cauchy sequence in (F(K, E), d F ). Since this space is complete, this will show that this sequence converges in this space and will therefore complete the proof. Let ε > 0 and, using (76), take δ > 0 and M ∈ N such that ω(s, s ) ≤ ε whenever d K (s, s ) ≤ δ and δ m ≤ ε whenever m ≥ M . Select a finite set {s l1 , . . . , s l N } such that any s ∈ K is within distance δ of a s li . Then, for any m, m ≥ M , by (76), The following lemma states an equivalent condition for the uniform convergence of functions, which proves extremely useful to establish uniform-in-time convergence of numerical schemes for parabolic equations when no smoothness is assumed on the data. Let us now prove the converse by contraposition. If (v m ) m∈N does not converge to v for d F then there exists ε > 0 and a subsequence (v m k ) k∈N , such that, for any k ∈ N, sup s∈K d E (v m k (s), v(s)) ≥ ε. We can then find a sequence (r k ) k∈N ⊂ K such that, for any k ∈ N,

d E (v m k (r k ), v(r k )) ≥ ε/2. ( 77
)
K being compact, up to another subsequence denoted the same way, we can assume that r k → s in K as k → ∞. It is then trivial to construct a sequence (s m ) m∈N converging to s and such that s m k = r k (just take s m = s when m is not an m k ). We then have v m (s m ) → v(s) in E and, by continuity of v, v(s m ) → v(s) in E. This shows that d E (v m (s m ), v(s m )) → 0, which contradicts (77) and concludes the proof.

Uniform-in-time convergence of numerical solutions to schemes for parabolic equations starts with a weak convergence with respect to the time variable. This weak convergence is then used to prove a stronger convergence. The following definition and proposition recall standard notions related to the weak topology on L 2 (Ω). t (n) -τ dt = τ , thus (79) is proven. We now turn to the proof of (80). We define the function χ(n, t) by χ(n, t) = 1 if n(t) = n, else χ(n, t) = 0. We have

ζ

  ∈ C 0 (R) is non-decreasing, Lipschitz continuous with Lipschitz constant L ζ > 0 such that ζ(0) = 0 and, for some M 0 , M 1 > 0, |ζ(s)| ≥ M 0 |s| -M 1 for all s ∈ R. (2b) β is a non-decreasing Lipschitz continuous function with Lipschitz constant L β > 0 and β(0) = 0. (2c) ∀s ∈ R , ν(s) = s 0

Definition 2 . 5 (

 25 Coercivity) A sequence (D m ) m∈N of space-time gradient discretisations in the sense of Definition 2.1 is said to be coercive if there exists C D such that, for any m ∈ N and any v ∈ X Dm,0 , ||Π Dm v|| L p (Ω) ≤ C D ||∇ Dm v|| L p (Ω) d . Definition 2.6 (Consistency) If D is a space-time gradient discretisation in the sense of Definition 2.1, we define

Remark 3 . 2

 32 The bound on |δ m v m | ,Dm is often a consequence on a numerical scheme satisfied by v m and on bounds on ||∇ Dm v m || L p (Ω×(0,T )) , see the proof of Lemma 4.3 for example.

Lemma 3 . 3

 33 Under Assumptions (2), there holds |ν(a) -ν(b)| ≤ L β |ζ(a) -ζ(b)|,

Lemma 3 . 4

 34 Let I be an interval of R and H : I → R be a convex function. We denote by L 2 (Ω; I) the convex set of functions in L 2 (Ω) with values in I. Let v ∈ L 2 (Ω; I) and (v m ) m∈N a sequence of functions in L 2 (Ω; I) which converges weakly to v in L 2 (Ω).

Lemma 4 . 4 (

 44 Estimate on the time translates of ν(u)) Under Assumptions (2), let D be a space-time gradient discretisation in the sense of Definition 2.1.

  (a) -ν(b)| ≤ L β |ζ(a) -ζ(b)|, which implies that, using Definition 2.7 with v = ζ(u m ) and (40),

  ), we have β = Id and therefore β = u = β(u), and ν = ζ. The strong convergence of Π Dm ν(u m ) = Π Dm ζ(u m ) to ν = ζ thus allows us to apply Lemma 3.5 to see that ζ = ζ(u). Exchanging the roles of β and ζ, we see that β = β(u) and ζ = ζ(u) still hold in the second case of

By Lemma 6 . 3 ,

 63 the uniform-in-time weak convergence of β(Π Dm u m ) to β(ū) and the weak continuity of β(ū) : [0, T ] → L 2 (Ω), we have β(Π Dm u m )(T m ) → β(ū)(T 0 ) weakly in L 2 (Ω) as m → ∞. Therefore, for any (s m ) m∈N converging to T 0 , 1 2 (β(Π Dm u m (T m )) + β(u)(s m )) → β(u)(T 0 ) weakly in L 2 (Ω) as m → ∞ and Lemma 3.4 gives, by convexity of B, Ω B(β(u)(x, T 0 ))dx ≤ lim inf m→∞ Ω B β(Π Dm u m (x, T m )) + β(u)(x, s m ) 2 dx. (60) Property (28) of B and the two inequalities (59) and (60) allow us to conclude the proof. Let (s m ) m∈N be a sequence in T (see the proof of Corollary 3.8) which converges to T 0 . Then ν(u(•, s m )) → ν(u)(•, T 0 ) in L 2 (Ω) as m → ∞. Using (28), we get

Remark 4 . 5

 45 Corollary 3.8), the first term in the right-hand side has a finite lim sup, bounded above by 16L β L ζ Ω B(β(u)(x, T 0 ))dx. We can therefore split the lim sup of this right-hand side without risking writing ∞ -∞ and we get, thanks to (60),lim sup m→∞ ν(Π Dm u m (•, T m )) -ν(u)(•, T 0 ) 2 L 2 (Ω) ≤ 0. Thus, ν(Π Dm u m (•, T m )) → ν(u)(T 0 ) strongly in L 2 (Ω). By Lemma 6.3 and the continuity of ν(u) stated in Corollary 3.8, this concludes the proof of the convergence of ν(Π Dm u m ) to ν(u) in L ∞ (0, T ; L 2 (Ω)). Since β(Π Dm u m )(T m ) → β(ū)(T 0 ) weakly in L 2 (Ω) as m → ∞, Lemma 3.4 also shows that Ω B(β(u)(x, T 0 ))dx ≤ lim inf m→∞ Ω B(β(Π Dm u m )(x, T m ))dx,and therefore, combined with (59), that lim m→∞ Ω B(β(Π Dm u m (x, T m )))dx = Ω B(β(u)(x, T 0 ))dx. (61) 4.4 Proof of Theorem 2.16

  ) where C does not depend on m and ω(δ) → 0 as δ → 0. Note that a similar estimate holds with Π Dm ζ m replaced with ζ since ζ ∈ L 2 (Ω × (0, T )) and, therefore, its translates tend to 0 in L 2 (Ω × (0, T )). If we respectively denote by A m (Π Dm ζ m ) and A(ζ) the left-hand side and right-hand side (64), then since (Π Dm β m ) m∈N is bounded in L 2 (Ω × (0, T )) we have, by (66),

+ ζ 2 - 1 β + ζ 2 .

 212 )β(w m (z))ζ(w m (z))dz = V ϕ(z)β(z) ζ(z)dz. (72) Then β = β(w) and ζ = ζ(w) a.e. in V , where w = β

d

  E (v m (s), v m (s)) ≤ d E (v m (s), v m (s li )) + d E (v m (s li ), v m (s li )) +d E (v m (s li ), v m (s)) ≤ ω(s, s li ) + δ m + d E (v m (s li ), v m (s li )) + ω(s, s li ) + δ m ≤ 4ε + d E (v m (s li ), v m (s li )).Since {(v m (s li )) m∈N : i = 1, . . . , N } form a finite number of converging sequences in E, we can find M ≥ M such that, whenever m, m ≥ M and i = 1, . . . , N , d E (v m (s li ), v m (s li )) ≤ ε. This shows that, for all m, m ≥ M and all s ∈ K, d E (v m (s), v m (s)) ≤ 5ε and concludes the proof that (v m ) m∈N is a Cauchy sequence in (F(K, E), d F ).

Lemma 6 . 3

 63 Let (K, d K ) be a compact metric space, (E, d E ) be a metric space and (F(K, E), d F ) as in Definition 6.1. Let (v m ) m∈N be a sequence in F(K, E) and v : K → E be continuous. Then v m → v for d F if and only if, for any s ∈ K and any sequence (s m ) m∈N ⊂ K converging to s for d K , we have v m (s m ) → v(s) for d E .Proof. If v m → v for d F then for any sequence (s m ) m∈N converging to s,d E (v m (s m ), v(s)) ≤ d E (v m (s m ), v(s m )) + d E (v(s m ), v(s)) ≤ d F (v m , v) + d(v(s m ), v(s))and the right-hand side tends to 0 by definition of v m → v for d F and by continuity of v.

  ) χ(m, t + ζ)dt, (n) ≤t+τ (t (n+1) -t (n) ) ≤ τ + δt, we can write from (81) n(t+ζ)+1) dt ≤ (τ + δt) m∈Z a (m+1) t (m+1) -ζ t (m) -ζ dt = (τ + δt) m∈Z a (m+1) δt (m+ 1 2 ) ,

  Throughout this paper, Ω is an open bounded subset of R d , d ∈ N , and p ∈ (1, +∞).

	Definition 2.1 (Space-Time gradient discretisation for homogeneous Dirichlet bound-
	ary conditions)

  it is greater than a linear functional and Φ(w) is thus well defined in (∞, ∞]. Moreover, if w k → w strongly in L 2 (Ω; I) then, up to a subsequence, w k → w a.e. on Ω and therefore H(w k ) → H(w) a.e. on Ω. Combined with the linear lower bound of H, we can apply Fatou's lemma to see that Φ(w) ≤ lim inf k→∞ Φ(w k ).

  n(t) the element n ∈ Z such that t ∈ [t(n) , t (n+1) ). Let (a (n) ) n∈Z be a family of non negative real values with a finite number of non zero values. Then (n+ 1 2 ) a (n+1) ), ∀τ ∈ (0, +∞), ∀ζ ∈ R. (80)Proof.Let us define the function χ(t, n, τ ) by χ(t, n, τ ) = 1 if t < t (n) and t + τ ≥ t (n) , else χ(t, n, τ ) = 0.

				n(t+τ )
				(δt (n+ 1 2 ) a (n+1) )dt = τ	(δt (n+ 1 2 ) a (n+1) ), ∀τ ∈ (0, +∞),	(79)
		R	n=n(t)+1	n∈Z
	and
		n(t+τ )		
		δt (n+ 1 2 )	 a n(t+ζ)+1 dt
	R	n=n(t)+1
				≤ (τ + δt)
	We have
				n(t+τ )
				(δt (n+ 1 2 ) a (n+1) )dt =	(δt (n+ 1 2 ) a (n+1) χ(t, n, τ ))dt
		R	n=n(t)+1	R n∈Z
				=	δt (n+ 1 2 ) a (n+1)
				n∈Z

n∈Z (δt R χ(t, n, τ )dt . Since R χ(t, n, τ )dt = t (n)

which is exactly (80).

Indeed, if this is not the case then we have s 1 < s 2 such that (β + ζ)(s 1 ) = (β + ζ)(s 2 ) and, since β and ζ are non-decreasing, that [s 1 , s 2 ] is a common plateau of β and ζ. Denoting by β, ζ and ν the functions obtained from β, ζ and ν by removing this common plateau (by a contraction of the s-ordinate), we see that u is a solution to [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF] if and only if u is a solution of the same problem with β, ζ and ν replaced with β, ζ and ν. Once this reduction has been applied, we can state the following theorem.

Theorem 5.2 Under Assumptions (2), let (D m ) mp∈N be a sequence of space-time gradient discretisations in the sense of Definition 2.1, which is coercive, consistent, compact and limit-conforming (see Section 2.2). Let, for any m ∈ N, u m be the solution to [START_REF] Droniou | Gradient schemes for elliptic and parabolic problems[END_REF] with D = D m , provided by Theorem 2.10. If we assume that Assumption (63) hold and that p ≥ 2, then there exists a solution u to (4) such that, up to a subsequence:

• the convergences stated in (18) hold,

Proof.

We only need to prove the first conclusion of the theorem, i.e. that the convergences (18) hold. Indeed, Theorems 2.14 and 2.16 then provide the last two conclusions. The only difference in the assumptions of Theorem 5.2 and 2.10 is the absence, here, of the structural assumption [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF] and the assumption that β + ζ is strictly increasing. The only place in the proof of Theorem (2.10) where we used Assumption [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF] is in Step 1, to identify the limits β, ζ and ν of Π Dm β(u m ), Π Dm ζ(u m ) and Π Dm ν(u m ). We will show that without assuming [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF] but by assuming (63), we can still identify those limits. We first work with the first two limits. We now turn to ν. Estimates (50) and ( 43) and Kolmogorov's compactness theorem show that the convergence of Π Dm ν(u m ) towards ν is actually strong on L 2 (Ω × (0, T )) (we use p ≥ 2 here).

We can therefore apply Lemma 5.6 with ϕ ≡ 1 and µ, ν instead of β, ζ, and we find a measurable w such that ν = ν(w) and µ = µ(w). The second relation translates into

To summarise, the limits of Π Dm β(u m ), Π Dm ζ(u m ) and Π Dm ν(u m ) have been identified as β(u), ζ(u) and ν(u) for some u. This allows to take over the proof of Theorem 2.10 from after the usage of [START_REF] Droniou | Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models[END_REF] and conclude that u is a solution to (4) and that the convergences (18) hold. The last two conclusions of the theorem follow from Theorems 2.14 and 2.16.

Proof. We first notice that β(w), ζ(w) ∈ L 2 (V ), given that β+ζ 2 ∈ L 2 (V ) and β • ( β+ζ 2 ) -1 and ζ • ( β+ζ 2 ) -1 are sub-linear (the sum of both is 2Id and each one has the same sign as s). Since β and ζ are non-decreasing, we can therefore write

Letting m → ∞ in the above inequality, and using the convergences of β(w m ), ζ(w m ) and (72), we obtain

We then remark that

(w). Hence, (74) leads to

Since ϕ is almost everywhere strictly positive on V , we deduce that

a.e. in V , and (73) follows from (75).

6 Appendix: uniform-in-time compactness results for timedependent problems

We establish in this appendix some generic results, unrelated to the framework of Gradient Schemes, that form the starting point for our uniform-in-time convergence results.

Solutions of numerical schemes for parabolic equations are usually piecewise constant, and therefore not continous, in time. As their jump nevertheless tend to become small with the time step, it is possible to establish some uniform-in-time convergence results using a generalisation to noncontinuous functions of the classical Ascoli-Arzela theorem. Definition 6.1 If (K, d K ) and (E, d E ) are metric spaces, we denote by F(K, E) the space of functions K → E, endowed with the uniform metric d F (v, w) = sup s∈K d E (v(s), w(s)) (note that this metric may take infinite values).

Theorem 6.2 (Generalised Ascoli-Arzela's theorem) Let (K, d K ) be a compact metric space, (E, d E ) be a complete metric space and (F(K, E), d F ) as in Definition 6.1. Let (v m ) m∈N be a sequence in F(K, E) such that there exists a function ω :

We also assume that, for all s ∈ K, {v m (s) :

Then (v m ) m∈N is relatively compact in (F(K, E), d F ) and any adherence value of (v m ) m∈N in this space is continuous K → E.

The next result is classical, but its short proof is recalled for the reader's convenience.

Proposition 6.5 Let E be closed bounded ball in L 2 (Ω) and (ϕ l ) l∈N be a dense sequence in L 2 (Ω).

Then, on E, the weak topology of L 2 (Ω) is the topology given by the metric Proof. The sets

(Ω) and ε > 0, define a basis of neighborhood of 0 for the weak L 2 (Ω) topology on E, and a basis of neighborhood of any other points is obtained by translation of this particular basis. If R is the radius of the ball E then for any ϕ ∈ L 2 (Ω), l ∈ N and v ∈ E we have

By density of (ϕ l ) l∈N we can select l ∈ N such that ||ϕ -ϕ l || L 2 (Ω) < ε/(2R) and we then see that E ϕ l ,ε/2 ⊂ E ϕ,ε . Hence, a basis of neighborhood of 0 in E for the weak L 2 (Ω) is also given by (E ϕ l ,ε ) l∈N, ε>0 .

From the definition of d E we see that, for any l ∈

Hence, any neighborhood of 0 in E for the L 2 (Ω) weak topology is a neighborhood of 0 for d E .

Conversely, for any ε > 0, selecting N ∈ N such that l≥N +1 2 -l < ε/2 gives, from the definition (78

Hence, any ball for d E centered at 0 is a neighborhood of 0 for the L 2 (Ω) weak topology. Since d E and the L 2 (Ω) weak neighborhood are invariant by translation, this concludes the proof that this weak topology is identical to the topology generated by d E .

The conclusion on weak uniform convergence of sequences of functions follows from the preceding result, and more precisely by noticing that all previous inclusions are, when applied to u m (t)-u(t), uniform with respect to t ∈ [0, T ].

The following lemma has been established in [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]Proposition 9.3] but its proof is recalled for the reader's convenience.

Lemma 6.6 Let (t (n) ) n∈Z be a stricly increasing sequence of real values such that δt (n+ 1 2 ) := t (n+1) -t (n) is uniformly bounded by δt > 0, lim n→-∞ t (n) = -∞ and lim n→∞ t (n) = ∞. For all t ∈ R, we denote by