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

Jérôme Droniou, Robert Eymard

► To cite this version:

Jérôme Droniou, Robert Eymard. Uniform-in-time convergence result of numerical methods for non-linear parabolic equations. Numerische Mathematik, 2015, 10.1007/s00211-015-0733-6. hal-00949682v3

HAL Id: hal-00949682

https://hal.science/hal-00949682v3

Submitted on 28 Oct 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Uniform-in-time convergence result of numerical methods for non-linear parabolic equations

J. Droniou* and R. Eymard†

October 28, 2014

## Abstract

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.

**AMS Subject Classification:** 46N40

# 1 Introduction

## 1.1 Motivation

The following generic nonlinear parabolic model

$$ 
\begin{aligned}
\partial_t \beta(\bar{u}) - \operatorname{div}(a(x, \nu(\bar{u}), \nabla \zeta(\bar{u}))) &= f \quad \text{in } \Omega \times (0, T), \\
\beta(\bar{u})(x, 0) &= \beta(u_{\text{ini}})(x) \quad \text{in } \Omega, \\
\zeta(\bar{u}) &= 0 \quad \text{on } \partial\Omega \times (0, T)
\end{aligned} 
\quad (1) 
$$

where $\beta, \zeta$ are non-decreasing, $\nu$ is such that $\nu' = \beta'\zeta'$ 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 $\zeta(s) = s$, $\nu = \beta$ and $a(x, \nu(\bar{u}), \nabla \zeta(\bar{u})) = \Lambda(x)K(x, \beta(\bar{u}))\nabla \bar{u}$, which describes the flow of water in a heterogeneous anisotropic underground medium,

2. Stefan's model [4], setting $\beta(s) = s$, $\nu = \zeta$, $a(x, \nu(\bar{u}), \nabla \zeta(\bar{u})) = K(\zeta(\bar{u}))\nabla \zeta(\bar{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 $\beta(s) = \zeta(s) = \nu(s) = s$ and $a(x, \nu(\bar{u}), \nabla \zeta(\bar{u})) = |\nabla \bar{u}|^{p-2}\nabla \bar{u}$, which is involved in the motion of glaciers [29] or flows of incompressible turbulent fluids through porous media [11].

The numerical approximation of these models has been extensively studied in the literature (see the fundamental work on the Stefan problem [34], and [35] for a review of some numerical approximations, see [32] for the Richards problem and see [12, 16] 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

<footnote>*School of Mathematical Sciences, Monash University, Victoria 3800, Australia. jerome.droniou@monash.edu.</footnote>

<footnote>†Université Paris-Est, Laboratoire d'Analyse et de Mathématiques Appliquées, UMR 8050, 5 boulevard Descartes, Champs-sur-Marne 77454 Marne-la-Vallée Cedex 2, France. Robert.Eymard@univ-mlv.fr.</footnote>
establishing space-time averaged results (e.g. in $L^2(\Omega \times (0, T)))$, see for example [23, 26]). Yet, the quantity of interest is often not $\bar{u}$ on $\Omega \times (0, T)$ but $\bar{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^\infty(0, T; L^2(\Omega))$ on $u - \bar{u}$, where $u$ is the approximated solution. Establishing such error estimates is however only feasible when uniqueness of the solution $\bar{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 (1) or even $p$-Laplace problems. It moreover requires some regularity assumptions on $\bar{u}$, which clearly fail for (1) (and simpler $p$-Laplace problems) for which, because of the possible plateaux of $\beta$ and $\zeta$, 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^\infty(0, T; L^2(\Omega))$ convergence result can be established for numerical approximations of (1), without having to assume non-physical regularity assumptions on the data. Note that, although Richards' and Stefan's models are formally equivalent when $\beta$ and $\zeta$ are strictly increasing (consider $\beta = \zeta^{-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 $\zeta$ is constant on the range of the solution) and Richards' model can become a non-transient elliptic equation (if $\beta$ is constant on this range). The innovative technique we develop in this paper is nonetheless generic enough to work directly on (1) 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 [7],

* The Crouzeix-Raviart non-conforming Finite Element, with or without mass lumping [9, 20],

* The Raviart-Thomas Mixed Finite Elements [5],

* The Vertex Approximate Gradient scheme [24],

* The Hybrid Mimetic Mixed family [15], which includes Mimetic Finite Differences [6], Mixed Finite Volume [13] and the SUSHI scheme [22],

* The Discrete Duality Finite Volume scheme in dimension 2 [30, 2], and the CeVeFE-Discrete Duality Finite Volume scheme in dimension 3 [8],

* The Multi-Point Flux Approximation O-method [1, 18].

We refer the reader to [14, 16, 21, 27] for more details. Let us finally emphasize that the unified convergence study of numerical schemes for Problem (1), which combines a general Leray-Lions operator and nonlinear functions $\beta$ or $\zeta$, 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 (1) and, in Section 1.3, we give an overview of the ideas involved in the proof
of uniform-in-time convergence. This overview is given not in a numerical analysis context but in a context of pure stability analysis of (1) with very little regularity on the data, for which the uniform-in-time convergence result also seems to be new. Section 2 presents the Gradient Schemes for our generic model (1). We give in Section 3 some preliminaries to the convergence study. Section 4 contains the complete convergence proof of Gradient Schemes for (1), including the uniform-in-time convergence result. An appendix, Section 5, concludes the article with important technical results, and in particular a generalisation of Ascoli-Arzela compactness result to discontinuous functions and a characterisation of uniform convergence of sequences of functions which is critical to establishing our uniform-in-time convergence result.

Note that these results and their proofs have been sketched and illustrated by some numerical examples in [17], for $a(x, \nu(\bar{u}), \nabla\zeta(\bar{u})) = \nabla\zeta(\bar{u})$.

## 1.2 Hypotheses and weak sense for the continuous problem

We consider the evolution problem (1) under the following hypotheses.

$$ \Omega \text{ is an open bounded connected polyhedral subset of } \mathbb{R}^d \ (d \in \mathbb{N}^*) \text{ and } T > 0, \quad (2a) $$

$$ \zeta \in C^0(\mathbb{R}) \text{ is non-decreasing, Lipschitz continuous with Lipschitz constant } L_\zeta > 0 \quad (2b) \\ \text{such that } \zeta(0) = 0 \text{ and, for some } M_0, M_1 > 0, |\zeta(s)| \ge M_0|s| - M_1 \text{ for all } s \in \mathbb{R}. $$

$$ \beta \text{ is a non-decreasing Lipschitz continuous function with Lipschitz constant } L_\beta > 0 \quad (2c) \\ \text{and } \beta(0) = 0. $$

$$ \forall s \in \mathbb{R}, \quad \nu(s) = \int_{0}^{s} \zeta'(q)\beta'(q) dq. \quad (2d) $$

$$ a : \Omega \times \mathbb{R} \times \mathbb{R}^d \to \mathbb{R}^d, \text{ with } p \in (1, +\infty), \text{ is a Caratheodory function}, \quad (2e) $$

(i.e. a function such that, for a.e. $x \in \Omega$, $(s, \xi) \mapsto a(x, s, \xi)$ is continuous and, for any $(s, \xi) \in \mathbb{R} \times \mathbb{R}^d$, $x \mapsto a(x, s, \xi)$ is measurable)

$$ \exists \underline{a} \in (0, +\infty) : a(x, s, \xi) \cdot \xi \ge \underline{a}|\xi|^p, \text{ for a.e. } x \in \Omega, \forall s \in \mathbb{R}, \forall \xi \in \mathbb{R}^d, \quad (2f) $$

$$ (a(x, s, \xi) - a(x, s, \chi)) \cdot (\xi - \chi) \ge 0, \text{ for a.e. } x \in \Omega, \forall s \in \mathbb{R}, \forall \xi, \chi \in \mathbb{R}^d, \quad (2g) $$

$$ \exists \overline{a} \in L^p(\Omega), \exists \mu \in (0, +\infty) : \quad \begin{array}{l} |a(x, s, \xi)| \le \overline{a}(x) + \mu|\xi|^{p-1}, \text{ for a.e. } x \in \Omega, \forall s \in \mathbb{R}, \forall \xi \in \mathbb{R}^d, \quad (2h) \\ \text{ closest } t \text{ to } 0 \text{ such that } \beta(t) = s. \end{array} $$

and

$$ u_{\text{ini}} \in L^2(\Omega), \quad f \in L^{p'}(\Omega \times (0, T)). \quad (2i) $$

We denote by $R_\beta$ the range of $\beta$ and define the pseudo-inverse function $\beta_r : R_\beta \to \mathbb{R}$ of $\beta$ by

$$ \forall s \in R_{\beta}, \beta_r(s) = \begin{cases} \inf\{t \in \mathbb{R} | \beta(t) = s\} & \text{if } s \ge 0, \\ \sup\{t \in \mathbb{R} | \beta(t) = s\} & \text{if } s < 0, \\ \text{closest } t \text{ to } 0 \text{ such that } \beta(t) = s. \end{cases} \quad (3) $$

Since $\beta(0) = 0$, we notice that $\beta_r \ge 0$ on $R_\beta \cap \mathbb{R}^+$ and $\beta_r \le 0$ on $R_\beta \cap \mathbb{R}^-$. We then define $B : R_\beta \to [0, \infty]$ by

$$ B(z) = \int_{0}^{z} \zeta(\beta_{r}(s)) ds. $$

Since $\beta_r$ is non-decreasing, this expression is always well-defined in $[0, \infty)$. The signs of $\beta_r$ also ensure that $B$ is non-decreasing on $R_\beta \cap \mathbb{R}^+$ and non-increasing on $R_\beta \cap \mathbb{R}^-$. We can therefore extend $B$ to $\bar{R}_\beta$ by these limits (possibly $+\infty$) at the potential endpoints of $R_\beta$.
The precise notion of solution to (1) that we consider is then the following:

$$
\left\{
\begin{array}{l}
\bar{u} \in L^p(0, T; L^p(\Omega)), \zeta(\bar{u}) \in L^p(0, T; W_0^{1,p}(\Omega)), \\
B(\beta(\bar{u})) \in L^\infty(0, T; L^1(\Omega)), \beta(\bar{u}) \in C([0, T]; L^2(\Omega)-\text{w}), \partial_t \beta(\bar{u}) \in L^{p'}(0, T; W^{-1,p'}(\Omega)), \\
\beta(\bar{u})(\cdot, 0) = \beta(u_{\text{ini}}) \text{ in } L^2(\Omega), \\
\displaystyle\int_0^T \langle \partial_t \beta(\bar{u})(\cdot, t), \bar{v}(\cdot, t) \rangle_{W^{-1,p'}, W_0^{1,p}} dt & (4) \\
\displaystyle\int_0^T \int_{\Omega} \mathbf{a}(\mathbf{x}, \nu(\bar{u}(\mathbf{x}, t)), \nabla \zeta(\bar{u})(\mathbf{x}, t)) \cdot \nabla \bar{v}(\mathbf{x}, t) d\mathbf{x} dt = \int_0^T \int_{\Omega} f(\mathbf{x}, t) \bar{v}(\mathbf{x}, t) d\mathbf{x} dt, \\
\forall \bar{v} \in L^p(0; T; W_0^{1,p}(\Omega)).
\end{array}
\right.
$$

where $C([0, T]; L^2(\Omega)-\text{w})$ denotes the space of continuous functions $[0, T] \mapsto L^2(\Omega)$ for the weak-* topology of $L^2(\Omega)$. Here and in the following, we denote by $p'$ the dual exponent $\frac{p}{p-1}$ to $p$ and we remove the mention of $\Omega$ in the duality bracket $\langle \cdot, \cdot \rangle_{W^{-1,p'}, W_0^{1,p}} = \langle \cdot, \cdot \rangle_{W^{-1,p'}(\Omega), W_0^{1,p}(\Omega)}$.

**Remark 1.1** The derivative $\partial_t \beta(\bar{u})$ is to be understood in the usual sense of distributions on $\Omega \times (0, T)$. Since the set $\mathcal{T} = \{\sum_{i=1}^q \varphi_i(t)\gamma_i(\mathbf{x}) : q \in \mathbb{N}, \varphi_i \in C_c^\infty(0, T), \gamma_i \in C_c^\infty(\Omega)\}$ of tensorial functions in $C^\infty(\Omega \times (0, T))$ is dense in $L^p(0, T; W_0^{1,p}(\Omega))$, one can ensure that this distribution derivative $\partial_t \beta(\bar{u})$ belongs to $L^{p'}(0, T; W^{-1,p'}(\Omega)) = (L^p(0, T; W_0^{1,p}(\Omega))')'$ by checking that the linear form

$$ \varphi \in \mathcal{T} \mapsto \langle \partial_t \beta(\bar{u}), \varphi \rangle_{\mathcal{D}', \mathcal{D}} = - \int_0^T \int_{\Omega} \beta(\bar{u})(\mathbf{x}, t) \partial_t \varphi(\mathbf{x}, t) d\mathbf{x} dt $$

is continuous for the norm of $L^p(0, T; W_0^{1,p}(\Omega))$.

Note that the continuity property of $\beta(\bar{u})$ natural. Indeed, the PDE in the sense of distributions shows that $T_\varphi : t \mapsto \langle \beta(\bar{u})(t), \varphi \rangle_{L^2}$ belongs to $W^{1,1}(0, T)$, and is therefore continuous, for any $\varphi \in C_c^\infty(\Omega)$. The density in $L^2(\Omega)$ of such $\varphi$, combined with the fact that $\beta(\bar{u}) \in L^\infty(0, T; L^2(\Omega))$ (coming from $B(\beta(\bar{u})) \in L^\infty(0, T; L^1(\Omega))$ and (26)), proves the continuity of $T_\varphi$ for any $\varphi \in L^2(\Omega)$, that is to say the continuity of $\beta(\bar{u}) : [0, T] \to L^2(\Omega)-\text{w}$.

This notion of $\beta(\bar{u})$ as a function continuous in time is nevertheless a subtle one. It is to be understood in the sense that the function $(\mathbf{x}, t) \mapsto \beta(\bar{u}(\mathbf{x}, t))$ has an a.e. representative which is continuous $[0, T] \mapsto L^2(\Omega)-\text{w}$. In other words, there is a function $Z \in C([0, T]; L^2(\Omega)-\text{w})$ such that $Z(t)(\mathbf{x}) = \beta(\bar{u}(\mathbf{x}, t))$ for a.e. $(\mathbf{x}, t) \in \Omega \times (0, T)$. We must however make sure, when dealing with pointwise values in time to separate $Z$ from $\beta(\bar{u}(\cdot, \cdot))$ as $\beta(\bar{u}(\cdot, t_1))$ may not make sense for a particular $t_1 \in [0, T]$.

That being said, in order to adopt a simple notation, we will denote by $\beta(\bar{u})(\cdot, \cdot)$ the function $Z$, and by $\beta(\bar{u}(\cdot, \cdot))$ the a.e.-defined composition of $\beta$ and $\bar{u}$. Hence, it will make sense to talk about $\beta(\bar{u})(\cdot, t)$ for a particular $t_1 \in [0, T]$, and we will only write $\beta(\bar{u})(\mathbf{x}, t) = \beta(\bar{u}(\mathbf{x}, t))$ for a.e. $(\mathbf{x}, t) \in \Omega \times (0, T)$. Note that from this a.e. equality we can ensure that $\beta(\bar{u})(\cdot, \cdot)$ takes its values in the closure $\bar{R}_\beta$ of the range of $\beta$.

## 1.3 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):

$$
\begin{align}
& \partial_t \beta(\bar{u}_\varepsilon) - \operatorname{div}(\boldsymbol{a}_\varepsilon(\boldsymbol{x}, \nu(\bar{u}_\varepsilon), \nabla \zeta(\bar{u}_\varepsilon))) = f \quad \text{in } \Omega \times (0, T), \nonumber \\
& \beta(\bar{u}_\varepsilon)(\boldsymbol{x}, 0) = \beta(u_{\text{ini}})(\boldsymbol{x}) \quad \text{in } \Omega, \nonumber \\
& \zeta(\bar{u}_\varepsilon) = 0 \quad \text{on } \partial\Omega \times (0, T) \tag{5}
\end{align}
$$
where $\mathbf{a}_\varepsilon$ satisfies Assumptions (2e)-(2h) with constants not depending on $\varepsilon$ and $\mathbf{a}_\varepsilon \to \mathbf{a}$ pointwise as $\varepsilon \to 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 $\varepsilon \to 0$), which are compatible with basic compactness results that can be obtained on $(\bar{u}_\varepsilon)_\varepsilon$ and also correspond to the initial convergences (18) that can be obtained on numerical approximations of (1):

$$ 
\begin{aligned}
& \beta(\bar{u}_\varepsilon) \to \beta(\bar{u}) \text{ in } C([0,T]; L^2(\Omega)-\text{w}), \quad \nu(\bar{u}_\varepsilon) \to \nu(\bar{u}) \text{ strongly in } L^1(\Omega \times (0,T)), \\
& \zeta(\bar{u}_\varepsilon) \to \zeta(\bar{u}) \text{ weakly in } L^p(0,T; W_0^{1,p}(\Omega)), \\
& \mathbf{a}_\varepsilon(\cdot, \nu(\bar{u}_\varepsilon), \nabla\zeta(\bar{u}_\varepsilon)) \to \mathbf{a}(\cdot, \nu(\bar{u}), \nabla\zeta(\bar{u})) \text{ weakly in } L^p(\Omega \times (0,T))^d.
\end{aligned}
\quad (6) $$

We will prove from these convergences that, along the same subsequence, $\nu(\bar{u}_\varepsilon) \to \nu(\bar{u})$ strongly in $C([0,T]; L^2(\Omega))$, which is our uniform-in-time convergence result.

We start by noticing that the weak-in-space uniform-in-time convergence of $\beta(\bar{u}_\varepsilon)$ gives, for any $T_0 \in [0, T]$ and any family $(T_\varepsilon)_{\varepsilon>0}$ converging to $T_0$ as $\varepsilon \to 0$, $\beta(\bar{u}_\varepsilon)(T_\varepsilon, \cdot) \to \beta(\bar{u})(T_0, \cdot)$ weakly in $L^2(\Omega)$. 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

$$ \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0)) d\mathbf{x} \leq \liminf_{\varepsilon \to 0} \int_{\Omega} B(\beta(\bar{u}_\varepsilon)(\mathbf{x}, T_\varepsilon)) d\mathbf{x}. \quad (7) $$

The second step is to notice that, by (2g) for $\mathbf{a}_\varepsilon$,

$$ \int_0^{T_\varepsilon} \int_{\Omega} [\mathbf{a}_\varepsilon(\cdot, \nu(\bar{u}_\varepsilon), \nabla\zeta(\bar{u}_\varepsilon)) - \mathbf{a}_\varepsilon(\cdot, \nu(\bar{u}_\varepsilon), \nabla\zeta(\bar{u}))] \cdot [\nabla\zeta(\bar{u}_\varepsilon) - \nabla\zeta(\bar{u})] d\mathbf{x} dt \geq 0. $$

Developing this expression and using the convergences (6), we find that

$$ \liminf_{\varepsilon \to 0} \int_0^{T_\varepsilon} \int_{\Omega} \mathbf{a}_\varepsilon(\cdot, \nu(\bar{u}_\varepsilon), \nabla\zeta(\bar{u}_\varepsilon)) \cdot \nabla\zeta(\bar{u}_\varepsilon)(\mathbf{x}, t) d\mathbf{x} dt \geq \int_0^{T_\varepsilon} \int_{\Omega} \mathbf{a}(\cdot, \nu(\bar{u}), \nabla\zeta(\bar{u})) \cdot \nabla\zeta(\bar{u}) d\mathbf{x} dt. \quad (8) $$

We then establish the following formula:

$$ 
\begin{aligned}
& \int_{\Omega} B(\beta(\bar{u}_{\varepsilon}(\mathbf{x}, T_{\varepsilon}))) d\mathbf{x} + \int_{0}^{T_{\varepsilon}} \int_{\Omega} \mathbf{a}_{\varepsilon}(\mathbf{x}, \nu(\bar{u}_{\varepsilon}(\mathbf{x}, t)), \nabla\zeta(\bar{u}_{\varepsilon})(\mathbf{x}, t)) \cdot \nabla\zeta(\bar{u}_{\varepsilon})(\mathbf{x}, t) d\mathbf{x} dt \\
& = \int_{\Omega} B(\beta(u_{\text{ini}}(\mathbf{x}))) d\mathbf{x} + \int_{0}^{T_{\varepsilon}} \int_{\Omega} f(\mathbf{x}, t) \zeta(\bar{u}_{\varepsilon})(\mathbf{x}, t) d\mathbf{x} dt. \quad (9)
\end{aligned}
 $$

This energy estimate is formally obtained by multiplying (5) by $\zeta(\bar{u}_\varepsilon)$ and integrating by parts, using the fact that $(B \circ \beta)' = \zeta\beta'$ (see Lemma 3.3). The rigorous justification of (9) is however quite technical, see Lemma 3.6 and Corollary 3.8. Thanks to (8), we can pass to the lim sup in (9) and we find, using the same energy estimate with $(\bar{u}, \mathbf{a}, T_0)$ instead of $(\bar{u}_\varepsilon, \mathbf{a}_\varepsilon, T_\varepsilon)$,

$$ \limsup_{\varepsilon \to 0} \int_{\Omega} B(\beta(\bar{u}_{\varepsilon}(\mathbf{x}, T_{\varepsilon}))) d\mathbf{x} \leq \int_{\Omega} B(\beta(\bar{u}(\mathbf{x}, T_0))) d\mathbf{x}. \quad (10) $$

Combined with (7), this shows that $\int_{\Omega} B(\beta(\bar{u}_\varepsilon(\mathbf{x}, T_\varepsilon))) d\mathbf{x} \to \int_{\Omega} B(\beta(\bar{u}(\mathbf{x}, T_0))) d\mathbf{x}$. The uniform convexity of $B$ (see (28)) then allows us to deduce that $\nu(\bar{u}_\varepsilon(\cdot, T_\varepsilon)) \to \nu(\bar{u}(\cdot, T_0))$ strongly in $L^2(\Omega)$ and thus that $\nu(\bar{u}_\varepsilon) \to \nu(\bar{u})$ strongly in $C([0, T]; L^2(\Omega))$ (see Lemma 6.3).

**Remark 1.2** A close examination of this proof indicates that equality in the energy estimate (9) is not required for $\bar{u}_\varepsilon$. An inequality $\le$ 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 $\bar{u}$ satisfies the equivalent of (9) with an equal sign (or $\ge$).

# 2 Gradient discretisations and gradient schemes

## 2.1 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 [14, 24, 16] 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. Throughout this paper, $\Omega$ is an open bounded subset of $\mathbb{R}^d$, $d \in \mathbb{N}^*$, and $p \in (1, +\infty)$.

**Definition 2.1 (Space-Time gradient discretisation for homogeneous Dirichlet boundary conditions)**

We say that $\mathcal{D} = (X_{\mathcal{D},0}, \Pi_{\mathcal{D}}, \nabla_{\mathcal{D}}, \mathcal{I}_{\mathcal{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_{\mathcal{D},0}$ is a finite dimensional real vector space,

2. the linear mapping $\Pi_{\mathcal{D}}: X_{\mathcal{D},0} \to L^{\infty}(\Omega)$ is a piecewise constant reconstruction operator in the sense that there exists a set $I$ of degrees of freedom such that $X_{\mathcal{D},0} = \mathbb{R}^I$ and there exists a family $(\Omega_i)_{i \in I}$ of disjoint subsets of $\Omega$ such that $\bar{\Omega} = \bigcup_{i \in I} \bar{\Omega}_i$ and, for all $u = (u_i)_{i \in I} \in X_{\mathcal{D},0}$ and all $i \in I, \Pi_{\mathcal{D}}u = u_i$ on $\Omega_i$,

3. the linear mapping $\nabla_{\mathcal{D}}: X_{\mathcal{D},0} \to L^p(\Omega)^d$ gives a reconstructed discrete gradient. It must be chosen such that $\|\nabla_{\mathcal{D}} \cdot \|_{L^p(\Omega)^d}$ is a norm on $X_{\mathcal{D},0}$,

4. $\mathcal{I}_{\mathcal{D}}: L^2(\Omega) \to X_{\mathcal{D},0}$ is a linear interpolation operator,

5. $t^{(0)} = 0 < t^{(1)} < t^{(2)} < \dots < t^{(N)} = T$.

We then set $\delta^{(n+\frac{1}{2})} = t^{(n+1)} - t^{(n)}$, for $n = 0, \dots, N-1$, and $\delta t_{\mathcal{D}} = \max_{n=0,...,N-1} \delta t^{(n+\frac{1}{2})}$, and we define the dual semi-norm $|w|_{\star,\mathcal{D}}$ of $w \in X_{\mathcal{D},0}$ by

$$|w|_{\star,\mathcal{D}} = \sup \left\{ \int_{\Omega} \Pi_{\mathcal{D}} w(\mathbf{x}) \Pi_{\mathcal{D}} z(\mathbf{x}) d\mathbf{x} : z \in X_{\mathcal{D},0}, \|\nabla_{\mathcal{D}} z\|_{L^p(\Omega)^d} = 1 \right\}. \quad (11)$$

**Remark 2.2 (Boundary conditions)** Other boundary conditions can be seamlessly handled by Gradient Schemes, see [14].

**Remark 2.3 (Nonlinear function of the elements of $X_{\mathcal{D},0}$)** Let $\mathcal{D}$ be a gradient discretisation in the sense of Definition 2.1. For any $\chi: \mathbb{R} \to \mathbb{R}$ and any $u = (u_i)_{i \in I} \in X_{\mathcal{D},0}$, we define $\chi_I(u) \in X_{\mathcal{D},0}$ by $\chi_I(u) = (\chi_I(u)_i)_{i \in I}$ with $\chi_I(u)_i = \chi(u_i)$. As indicated by the subscript $I$, this definition depends on the choice of the degrees of freedom in $X_{\mathcal{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 $\Pi_D$ is a piecewise constant reconstruction is the following:

$$ \forall \chi : \mathbb{R} \mapsto \mathbb{R}, \forall u \in X_{D,0}, \quad \Pi_D \chi(u) = \chi(\Pi_D u). \qquad (12) $$

It is customary to also use the notations $\Pi_D$ and $\nabla_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} \subset X_{D,0}$, we set

$$ \begin{array}{ll} \text{for a.e. } \boldsymbol{x} \in \Omega, \Pi_D v(\boldsymbol{x}, 0) = \Pi_D v^{(0)}(\boldsymbol{x}) \text{ and } \forall n = 0, \dots, N-1, \forall t \in (t^{(n)}, t^{(n+1)}] : & \\ \Pi_D v(\boldsymbol{x}, t) = \Pi_D v^{(n+1)}(\boldsymbol{x}), \nabla_D v(\boldsymbol{x}, t) = \nabla_D v^{(n+1)}(\boldsymbol{x}) & (13) \\ \text{and } \delta_D v(t) = \delta_D^{(n+\frac{1}{2})} v := \frac{v^{(n+1)} - v^{(n)}}{\delta_D^{(n+\frac{1}{2})}} \in X_{D,0}. & \end{array} $$

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 $\mathcal{D} = (X_{D,0}, \Pi_D, \nabla_D, \mathcal{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} \subset X_{D,0}$ such that

$$ \left\{ \begin{array}{l} u^{(0)} = \mathcal{I}_D u_{\text{ini}} \text{ and, for all } v = (v^{(n)})_{n=1,...,N} \subset X_{D,0}, \\ \displaystyle \int_0^T \int_{\Omega} [\Pi_D \delta_D \beta(u)(\boldsymbol{x},t) \Pi_D v(\boldsymbol{x},t) + \boldsymbol{a}(\boldsymbol{x}, \Pi_D \nu(u)(\boldsymbol{x},t), \nabla_D \zeta(u)(\boldsymbol{x},t)) \cdot \nabla_D v(\boldsymbol{x},t)] d\boldsymbol{x}dt \quad (14) \\ \displaystyle = \int_0^T \int_{\Omega} f(\boldsymbol{x},t) \Pi_D v(\boldsymbol{x},t) d\boldsymbol{x}dt. \end{array} \right. $$

**Remark 2.4** 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 $\alpha \in [\frac{1}{2}, 1]$ and replacing the terms $u^{(n+1)}$ appearing in $\boldsymbol{a}(\boldsymbol{x}, \cdot, \cdot)$ in (14) for $t \in (t^{(n)}, t^{(n+1)}]$ with $u^{(n+\alpha)} = \alpha u^{(n+1)} + (1-\alpha)u^{(n)}$. All the results we establish for (14) would hold for such a scheme (see the treatment in [16]).

## 2.2 Properties of gradient discretisations

In order to establish the convergence of the associated gradient schemes, sequences of space-time gradient discretisations are required to satisfy four properties: coercivity, consistency, limit-conformity and compactness.

**Definition 2.5 (Coercivity)** A sequence $(D_m)_{m \in \mathbb{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 \in \mathbb{N}$ and any $v \in X_{D,m,0}$, $||\Pi_D v||_{L^p(\Omega)} \le C_D ||\nabla_D v||_{L^p(\Omega)^d}$.

**Definition 2.6 (Consistency)** If $\mathcal{D}$ is a space-time gradient discretisation in the sense of Definition 2.1, we define

$$ \forall \varphi \in L^2(\Omega) \cap W_0^{1,p}(\Omega), \quad \hat{S}_D(\varphi) = \min_{w \in X_{D,0}} \left( ||\Pi_D w - \varphi||_{L^{\max(p,2)}(\Omega)} + ||\nabla_D w - \nabla \varphi||_{L^p(\Omega)} \right). \quad (15) $$

A sequence $(D_m)_{m \in \mathbb{N}}$ of space-time gradient discretisations in the sense of Definition 2.1 is said to be consistent if

$$ \bullet \text{ for all } \varphi \in L^2(\Omega) \cap W_0^{1,p}(\Omega), \quad \hat{S}_{D_m}(\varphi) \to 0 \text{ as } m \to \infty, $$
* for all $\varphi \in L^2(\Omega)$, $\Pi_{\mathcal{D}_m} \mathcal{I}_{\mathcal{D}_m} \varphi \to \varphi$ in $L^2(\Omega)$ as $m \to \infty$, and

* $\partial_{\mathcal{D}_m} \to 0$ as $m \to \infty$.

**Definition 2.7 (Compactness)** If $\mathcal{D}$ is a space-time gradient discretisation in the sense of Definition 2.1, we define

$$ \forall \xi \in \mathbb{R}^d, T_{\mathcal{D}}(\xi) = \max_{v \in X_{\mathcal{D},0} \setminus \{0\}} \frac{\|\Pi_{\mathcal{D}} v(\cdot + \xi) - \Pi_{\mathcal{D}} v\|_{L^p(\mathbb{R}^d)}}{\|\nabla_{\mathcal{D}} v\|_{L^p(\Omega)}}, $$

where $\Pi_{\mathcal{D}} v$ has been extended by 0 outside $\Omega$.

A sequence $(\mathcal{D}_m)_{m \in \mathbb{N}}$ of space-time gradient discretisations is said to be compact if

$$ \lim_{\xi \to 0} \sup_{m \in \mathbb{N}} T_{\mathcal{D}}(\xi) = 0. $$

**Definition 2.8 (Limit-conformity)** If $\mathcal{D}$ is a space-time gradient discretisation in the sense of Definition 2.1, we define

$$ \forall \varphi \in W^{\text{div},p'}(\Omega), W_{\mathcal{D}}(\varphi) = \max_{u \in X_{\mathcal{D},0} \setminus \{0\}} \frac{\left| \int_{\Omega} (\nabla_{\mathcal{D}} u(\mathbf{x}) \cdot \varphi(\mathbf{x}) + \Pi_{\mathcal{D}} u(\mathbf{x}) \text{div} \varphi(\mathbf{x})) \, d\mathbf{x} \right|}{\|\nabla_{\mathcal{D}} u\|_{L^p(\Omega)^d}}. \quad (16) $$

A sequence $(\mathcal{D}_m)_{m \in \mathbb{N}}$ of space-time gradient discretisations in the sense of Definition 2.1 is said to be limit-conforming if, for all $\varphi \in W^{\text{div},p'}(\Omega)$, $W_{\mathcal{D}_m}(\varphi)$ tends to 0 as $m \to \infty$.

We refer the reader to [16, 14] for a proof of the following lemma.

**Lemma 2.9 (Consequence of limit-conformity)** Let $(\mathcal{D}_m)_{m \in \mathbb{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 \in \mathbb{N}$, $v_m = (v_m^{(n)})_{n=0,...,N_m} \subset X_{\mathcal{D}_m,0}$ be such that, with the notations in (13), $(\nabla_{\mathcal{D}_m} v_m)_{m \in \mathbb{N}}$ is bounded in $L^p(\Omega \times (0,T))$.

Then there exists $v \in L^p(0,T; W_0^{1,p}(\Omega))$ such that, up to a subsequence as $m \to \infty$, $\Pi_{\mathcal{D}_m} v_m \to v$ weakly in $L^p(\Omega \times (0,T))$ and $\nabla_{\mathcal{D}_m} v_m \to \nabla v$ weakly in $L^p(\Omega \times (0,T))^d$.

## 2.3 Main results

Our first theorem states weak or space-time averaged convergence properties of Gradient Schemes for (1). These results are quite classical and have already been established for Richards' and Stefan's models, see [21, 25]. The convergence proof we provide afterwards however covers more non-linear model, as $\mathbf{a}$ do not need to be linear with respect to $\xi$, and is more compact than the ones available in the literature.

**Theorem 2.10 (Convergence of the Gradient Scheme)** Under Assumptions (2), let $(\mathcal{D}_m)_{mp \in \mathbb{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 \in \mathbb{N}$ there exists a solution $u_m$ to (14) with $\mathcal{D} = \mathcal{D}_m$. Moreover, if we assume that

$$ (\forall s \in \mathbb{R}, \beta(s) = s) \quad \text{or} \quad (\forall s \in \mathbb{R}, \zeta(s) = s), \quad (17) $$
then there exists a solution $\bar{u}$ to (4) such that, up to a subsequence, the following convergences hold
as $m \to \infty$:

$$
\begin{align}
\Pi_{\mathcal{D}_m} \beta(u_m) &\to \beta(\bar{u}) \quad \text{weakly in } L^2(\Omega) \text{ uniformly on } [0,T] \quad \text{(see Definition 6.4)}, \nonumber \\
\Pi_{\mathcal{D}_m} \nu(u_m) &\to \nu(\bar{u}) \quad \text{strongly in } L^1(0,T; L^1(\Omega)), \nonumber \\
\Pi_{\mathcal{D}_m} \zeta(u_m) &\to \zeta(\bar{u}) \quad \text{weakly in } L^p(\Omega \times (0,T)). \tag{18} \\
\nabla_{\mathcal{D}_m} \zeta(u_m) &\to \nabla \zeta(\bar{u}) \quad \text{weakly in } L^p(\Omega \times (0,T))^d. \nonumber
\end{align}
$$

**Remark 2.11** Since $|\nu| \le L_\zeta |\beta|$ and $|\nu| \le L_\beta |\zeta|$, the $L^\infty(0,T; L^2(\Omega))$ bound on $\beta(u_m)$ and the $L^p(\Omega \times (0,T))$ bound on $\zeta(u_m)$ (see Lemma 4.1 and Definition 2.5) allows us to see that the strong convergence of $\Pi_{\mathcal{D}_m} \nu(u_m)$ is also valid in $L^q(0,T; L^r(\Omega))$ for any $(q,r)$ such that $q < \infty$ 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 $\bar{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 non-linear parabolic equations. This assumption is actually not mandatory if $p \ge 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 \in \mathbb{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 \in \mathbb{N}$, $u_m$ is a solution to (14) with $D = D_m$, which converges to a solution $\bar{u}$ of (4) in the sense (18).

Then, as $m \to \infty$, $\Pi_{D_m} \nu(u_m) \to \nu(\bar{u})$ strongly in $L^\infty(0, T; L^2(\Omega))$.

**Remark 2.15** Note that since $(\Pi_{D_m} \nu(u_m))_{m \in \mathbb{N}}$ are piecewise constant in time, their convergence in $L^\infty(0, T; L^2(\Omega))$ is actually a uniform-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 [16].

**Theorem 2.16 (Strong convergence of the gradients)** Under Assumptions (2), let $(D_m)_{mp \in \mathbb{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 \in \mathbb{N}$, $u_m$ is a solution to (14) with $D = D_m$, which converges to a solution $\bar{u}$ of (4) in the sense (18). We also assume that $\mathbf{a}$ is strictly monotone in the sense:

$$ (\mathbf{a}(\mathbf{x}, s, \boldsymbol{\xi}) - \mathbf{a}(\mathbf{x}, s, \boldsymbol{\chi})) \cdot (\boldsymbol{\xi} - \boldsymbol{\chi}) > 0, \text{ for a.e. } \mathbf{x} \in \Omega, \forall s \in \mathbb{R}, \forall \boldsymbol{\xi} \neq \boldsymbol{\chi} \in \mathbb{R}^d. \quad (19) $$

Then, as $m \to \infty$, $\Pi_{D_m} \zeta(u_m) \to \zeta(\bar{u})$ strongly in $L^p(\Omega \times (0, T))$ and $\nabla_{D_m} \zeta(u_m) \to \nabla \zeta(\bar{u})$ strongly in $L^p(\Omega \times (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.
# 3 Preliminaries

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

## 3.1 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(\Omega)$-weak compactness result)** Let $T > 0$, and $(D_m)_{m \in \mathbb{N}} = (X_{D_m,0}, \Pi_{D_m}, \nabla_{D_m}, \mathcal{I}_{D_m}, (t_m^{(n)})_{n=0,...,N_m})_{m \in \mathbb{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 \in \mathbb{N}$, let $v_m = (v_m^{(n)})_{n=0,...,N_m} \subset X_{D_m,0}$. If there exists $q > 1$ and $C > 0$ such that, for any $m \in \mathbb{N}$,

$$ ||\Pi_{D_m} v_m||_{L^\infty(0,T;L^2(\Omega))} \le C \quad \text{and} \quad \int_0^T |\delta_m v_m(t)|_{\star, D_m}^q dt \le C, \quad (20) $$

then the sequence $(\Pi_{D_m} v_m)_{m \in \mathbb{N}}$ is relatively compact uniformly-in-time and weakly in $L^2(\Omega)$, 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] \to L^2(\Omega)$ for the weak topology.

**Remark 3.2** The bound on $|\delta_m v_m|_{\star, D_m}$ is often a consequence on a numerical scheme satisfied by $v_m$ and on bounds on $||\nabla_{D_m} v_m||_{L^p(\Omega \times (0,T))}$, see the proof of Lemma 4.3 for example.

**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(\Omega)$, endowed with the weak topology. We let $(\varphi_l)_{l \in \mathbb{N}} \subset C_c^\infty(\Omega)$ be a dense sequence in $L^2(\Omega)$ and endow $E$ with the metric (78) from these $\varphi_l$, which indeed defines the weak $L^2(\Omega)$ topology (see Proposition 6.5).

The set $E$ is metric compact and therefore complete, and all $\Pi_{D_m} 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\varphi_l \in X_{D,0}$ by

$$ P_D\varphi_l = \underset{w \in X_{D,0}}{\operatorname{argmin}} \left( ||\Pi_D w - \varphi_l||_{L^{\max(p,2)}(\Omega)} + ||\nabla_D w - \nabla\varphi_l||_{L^p(\Omega)} \right). \quad (21) $$

For any $0 \le s \le s' \le T$, writing $\Pi_D v_m(s') - \Pi_D v_m(s)$ as the sum of its jumps $\delta^{(n+\frac{1}{2})}\Pi_D \delta_D^{(n+\frac{1}{2})} v_m$ at the points $(t^{(n)})_{n=n_1,...,n_2}$ lying between $s$ and $s'$, the definition of $|\cdot|_{\star, D}$ and Estimate (20) give

$$ \left| \int_{\Omega} (\Pi_D v_m(\mathbf{x}, s') - \Pi_D v_m(\mathbf{x}, s)) \Pi_D P_D \varphi_l(\mathbf{x}) d\mathbf{x} \right| = \left| \int_{t^{(n_1)}}^{t^{(n_2+1)}} \int_{\Omega} \Pi_D \delta_D v(t)(\mathbf{x}) \Pi_D P_D \varphi_l(\mathbf{x}) d\mathbf{x} dt \right| \le C^{1/q} (t^{(n_2+1)} - t^{(n_1)})^{1/q'} ||\nabla_D P_D \varphi_l||_{L^p(\Omega)^d}. $$

By definition of $P_D$, we have

$$ ||\Pi_D P_D \varphi_l - \varphi_l||_{L^2(\Omega)} \le \hat{S}_D(\varphi_l) $$

and

$$ ||\nabla_D P_D \varphi_l||_{L^p(\Omega)^d} \le \hat{S}_D(\varphi_l) + ||\varphi_l||_{L^p(\Omega)} + ||\nabla \varphi_l||_{L^p(\Omega)^d} \le C_{\varphi_l} $$
with $C_{\varphi_l}$ not depending on $\mathcal{D}$ (and therefore on $m$). Since $t^{(n_2+1)} - t^{(n_1)} \le |s' - s| + \partial t$ and $(\Pi_{\mathcal{D}}v_m)_{m \in \mathbb{N}}$ is bounded in $L^\infty(0, T; L^2(\Omega))$, we deduce that

$$ \left| \int_{\Omega} (\Pi_{\mathcal{D}} v_m(\mathbf{x}, s') - \Pi_{\mathcal{D}} v_m(\mathbf{x}, s)) \varphi_l(\mathbf{x}) d\mathbf{x} \right| \le 2C \hat{S}_{\mathcal{D}}(\varphi_l) + C^{1/q} C_{\varphi_l} |s' - s|^{1/q'} + C^{1/q} C_{\varphi_l} \partial t^{1/q'}. \quad (22) $$

Plugged into the definition (78) of the distance in $E$, this shows that

$$ \begin{aligned} d_E(\Pi_{\mathcal{D}} v_m(s'), \Pi_{\mathcal{D}} v_m(s)) & \le \sum_{l \in \mathbb{N}} \frac{\min(1, C^{1/q'} C_{\varphi_l} |s' - s|^{1/q'})}{2^l} + \sum_{l \in \mathbb{N}} \frac{\min(1, 2C \hat{S}_{\mathcal{D}_m}(\varphi_l) + C^{1/q'} C_{\varphi_l} \partial t_m^{1/q'})}{2^l} \\ & =: \omega(s, s') + \delta_m. \end{aligned} $$

Using the dominated convergence theorem for series, we see that $\omega(s, s') \to 0$ as $s - s' \to 0$ and that $\delta_m \to 0$ as $m \to \infty$ (we invoke the space consistency to establish that $\lim_{m \to \infty} \hat{S}_{\mathcal{D}_m}(\varphi_l) \to 0$ for any $l$). Hence, the assumptions of Theorem 6.2 are satisfied and the proof is complete. ■

## 3.2 Technical results

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

**Lemma 3.3** Under Assumptions (2), there holds

$$ |\nu(a) - \nu(b)| \le L_{\beta} |\zeta(a) - \zeta(b)|, \quad (23) $$

$$ (\nu(a) - \nu(b))^2 \le L_{\beta} L_{\zeta} (\zeta(a) - \zeta(b)) (\beta(a) - \beta(b)). \quad (24) $$

The function $B$ is convex on $\overline{R}_{\beta}$, the function $B \circ \beta : \mathbb{R} \to [0, \infty)$ is continuous,

$$ \forall s \in \mathbb{R}, \quad B(\beta(s)) = \int_{0}^{s} \zeta(q)\beta'(q)dq, \quad (25) $$

$$ \exists K_0, K_1, K_2 > 0 \text{ such that } \forall s \in \mathbb{R}, \quad K_0\beta(s)^2 - K_1 \le B(\beta(s)) \le K_2s^2, \quad (26) $$

$$ \forall a \in \mathbb{R}, \forall r \in \overline{R_{\beta}}, \quad B(r) - B(\beta(a)) \ge \zeta(a)(r - \beta(a)), \quad (27) $$

and

$$ \forall s, s' \in \mathbb{R}, \quad (\nu(s) - \nu(s'))^2 \le 4L_{\beta}L_{\zeta} \left[ B(\beta(s)) + B(\beta(s')) - 2B \left( \frac{\beta(s) + \beta(s')}{2} \right) \right]. \quad (28) $$

**Proof.**

Inequality (23) is a straightforward consequence of the estimate $\nu' = \zeta'\beta' \le L_{\beta}\zeta'$. Note that the same inequality also holds with $\beta$ and $\zeta$ swapped. Inequality (24) is then a direct application of (23) and the similar inequality with $\beta$ and $\zeta$ swapped.

Let us first notice that, since $\beta \ge 0$ on $\mathbb{R}^+$ and $\beta \le 0$ on $\mathbb{R}^-$, $\beta_r(s)$ is always a real number when $s \in R_{\beta}$. Moreover, since $\beta$ is non-decreasing, $\beta_r$ is also non-decreasing on $R_{\beta}$ and therefore locally bounded on $R_{\beta}$. Hence, $B$ is well defined and locally Lipschitz-continuous, with an a.e. derivative $B' = \zeta(\beta_r)$. $B'$ is therefore non-decreasing and $B$ is convex.
To prove (25), we denote by $P \subset R_\beta$ the countable set of plateaux values of $\beta$, i.e. the $y \in \mathbb{R}$ such that $\beta^{-1}(\{y\})$ is not reduced to a singleton. If $s \notin \beta^{-1}(P)$ then $\beta^{-1}(\{\beta(s)\})$ is the singleton $\{s\}$ and therefore $\beta_r(\beta(s)) = s$. Moreover, $\beta_r$ is continuous at $\beta(s)$ and thus $B$ is differentiable at $\beta(s)$. Since $\beta$ is differentiable a.e., we therefore deduce that, for a.e. $s \notin \beta^{-1}(P)$, $(B(\beta))'(s) = B'(\beta(s))\beta'(s) = \zeta(\beta_r(\beta(s)))\beta'(s) = \zeta(s)\beta'(s)$. The set $\beta^{-1}(P)$ is a union of intervals on which $\beta$ and thus $B(\beta)$ are locally constant; hence, for a.e. $s$ in this set, $(B(\beta))'(s) = 0$ and $\zeta(s)\beta'(s) = 0$. Hence, the locally Lipschitz-continuous functions $B(\beta)$ and $s \to \int_0^s \zeta(q)\beta'(q)dq$ have identical derivatives a.e. on $\mathbb{R}$ and take the same value at $s = 0$. They are thus equal on $\mathbb{R}$ and the proof of (25) is complete.

The continuity of $B \circ \beta$ is an obvious consequence of (25). The second inequality in (26) can also be easily deduced from (25) by noticing that $|\zeta(s)\beta'(s)| \le L_\zeta L_\beta |s|$ (we can take $K_2 = \frac{L_\beta L_\zeta}{2}$). For the first inequality in (26), we first infer from (2b) the existence of $S > 0$ such that $|\zeta(q)| \ge \frac{M_0}{2} |q| \ge \frac{M_0}{2L_\beta} |\beta(q)|$ whenever $|q| \ge S$. We then write, for $s \ge S$,

$$B(\beta(s)) = \int_{0}^{s} \zeta(q)\beta'(q)dq + \int_{s}^{s} \zeta(q)\beta'(q)dq \geq \frac{M_0}{2L_\beta} \int_{s}^{s} \beta(q)\beta'(q)dq = \frac{M_0}{4L_\beta} (\beta(s)^2 - \beta(S)^2).$$

A similar inequality holds for $s \le -S$ (with $\beta(-S)$ instead of $\beta(S)$) and the first inequality in (26) therefore holds with $K_0 = \frac{M_0}{4L_\beta}$ and $K_1 = \frac{M_0}{4L_\beta} \max[-S, S] \beta^2$.

We now prove (27), which states that $a$ belongs to the convex sub-differential of $B$ at $\beta(a)$. We first start with the case $r \in R_\beta$, that is $r = \beta(b)$ for some $b \in \mathbb{R}$. If $\beta_r$ is continuous at $\beta(a)$ then this inequality is an obvious consequence of the convexity of $B$ since $B$ is then differentiable at $\beta(a)$ with $B'(\beta(a)) = \zeta(\beta_r(\beta(a))) = \zeta(a)$. Otherwise, a plain reasoning also does the job as

$$B(r) - B(\beta(a)) = B(\beta(b) - B(\beta(a))) \\ \int_{a}^{b} \zeta(q)\beta'(q)dq = \int_{a}^{b} (\zeta(q) - \zeta(a))\beta'(q)dq + \zeta(a)(\beta(b) - \beta(a)) \ge \zeta(a)(r - \beta(a)),$$

the inequality coming from the fact that $\beta' \ge 0$ and $\zeta(q) - \zeta(a)$ has the same sign as $b - a$ when $q$ is between $a$ and $b$. The general case $r \in \bar{R}_\beta$ is obtained by passing to the limit on $b_n$ such that $\beta(b_n) \to r$ and using the fact that $B$ has limits (possibly $+\infty$) at the endpoints of $R_\beta$.
Let us now take $s, s' \in \mathbb{R}$. Let $\bar{s} \in \mathbb{R}$ be such that $\beta(\bar{s}) = \frac{\beta(s)+\beta(s')}{2}$. We notice that

$$B(\beta(s)) + B(\beta(s')) - 2B(\beta(\bar{s})) = \int_{\bar{s}}^{s} (\zeta(q) - \zeta(\bar{s}))\beta'(q)dq + \int_{\bar{s}}^{s'} (\zeta(q) - \zeta(\bar{s}))\beta'(q)dq.$$

We then use $|\zeta(q) - \zeta(\bar{s})| \ge \frac{1}{L_\beta} |\nu(q) - \nu(\bar{s})|$ and $\beta'(q) \ge \beta'(q) \frac{\zeta'(q)}{L_\zeta} = \frac{\nu'(q)}{L_\zeta}$ to write

$$\int_{\bar{s}}^{s} (\zeta(q) - \zeta(\bar{s}))\beta'(q)dq \ge \frac{1}{L_{\beta}L_{\zeta}} \int_{\bar{s}}^{s} \nu'(q)(\nu(q) - \nu(\bar{s}))dq = \frac{1}{2L_{\beta}L_{\zeta}} (\nu(s) - \nu(\bar{s}))^2.$$

Thanks to

$$(\nu(s) - \nu(s'))^2 \le 2((\nu(s) - \nu(\bar{s}))^2 + (\nu(s') - \nu(\bar{s}))^2),$$

we deduce that (28) follows. ■

The next lemma is an easy consequence of Fatou's lemma and the fact that strongly lower semi-continuous convex functions are also weakly lower semi-continuous. We all the same provide its short proof.
**Lemma 3.4** Let $I$ be an interval of $\mathbb{R}$ and $H: I \to \mathbb{R}$ be a convex function. We denote by $L^2(\Omega; I)$ the convex set of functions in $L^2(\Omega)$ with values in $I$. Let $v \in L^2(\Omega; I)$ and $(v_m)_{m \in \mathbb{N}}$ a sequence of functions in $L^2(\Omega; I)$ which converges weakly to $v$ in $L^2(\Omega)$. Then

$$ \int_{\Omega} H(v(\boldsymbol{x})) d\boldsymbol{x} \leq \liminf_{m \to \infty} \int_{\Omega} H(v_m(\boldsymbol{x})) d\boldsymbol{x}. $$

**Proof.**

Let $\Phi: L^2(\Omega; I) \to (-\infty, \infty]$ be defined by $\Phi(w) = \int_{\Omega} H(w(\boldsymbol{x})) d\boldsymbol{x}$. Since $H$ is convex, it is greater than a linear functional and $\Phi(w)$ is thus well defined in $(\infty, \infty]$. Moreover, if $w_k \to w$ strongly in $L^2(\Omega; I)$ then, up to a subsequence, $w_k \to w$ a.e. on $\Omega$ and therefore $H(w_k) \to H(w)$ a.e. on $\Omega$. Combined with the linear lower bound of $H$, we can apply Fatou's lemma to see that $\Phi(w) \leq \liminf_{k \to \infty} \Phi(w_k)$.

Hence, $\Phi$ is lower semi-continuous for the strong topology of $L^2(\Omega; I)$. Since $\Phi$ (as $H$) is convex, we deduce that this lower semi-continuity property is also valid for the weak topology of $L^2(\Omega; I)$, see [19]. The result of the lemma is just the translation of this weak lower semi-continuity of $\Phi$. ■

The last technical result is a consequence of the Minty trick. It has been proved and used in the $L^2$ case in [21, 14], but we need here an extension to the non-Hilbertian case.

**Lemma 3.5 (Minty's trick)** Let $H \in C^0(\mathbb{R})$ be a nondecreasing function. Let $(X, \mu)$ be a measurable set with finite measure and let $(u_n)_{n \in \mathbb{N}} \subset L^p(X)$ with $p > 1$ such that

1. there exists $u \in L^p(X)$ such that $(u_n)_{n \in \mathbb{N}}$ weakly converges to $u$ in $L^p(X)$;

2. $(H(u_n))_{n \in \mathbb{N}} \subset L^1(X)$ and there exists $w \in L^1(X)$ such that $(H(u_n))_{n \in \mathbb{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 \to \infty$ as $k \to \infty$ such that $H(T_k(s)) = T_{h_k}(H(s))$. Thus, $H(T_k(u_n)) \to T_{h_k}(w)$ in $L^1(X)$ as $n \to \infty$. Given that $(H(T_k(u_n)))_{n \in \mathbb{N}}$ remains bounded in $L^\infty(X)$, its convergence to $T_{h_k}(w)$ also holds in $L^{p'}(X)$.

Using fact that $H \circ T_k$ is non-decreasing, we write, for any $g \in L^p(X)$,

$$ \int_X (H(T_k(u_n)) - H(T_k(g)))(u_n - g)d\mu \geq 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 \circ T_k$ is bounded, we can take the limit of this expression as $n \to \infty$ and we find

$$ \int_X (T_{h_k}(w) - H(T_k(g)))(u - g)d\mu \geq 0. \quad (29) $$

We then use Minty's trick, that is pick a generic $\varphi \in L^p(X)$, apply (29) to $g = u - t\varphi$, divide by $t$ and let $t \to \pm 0$ (using the dominated convergence theorem and the fact that $H \circ T_k$ is continuous and bounded) to find

$$ \int_X (T_{h_k}(w) - H(T_k(u)))\varphi d\mu = 0. $$

Selecting $\varphi = \text{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 \to \infty$, we conclude that $w = H(u)$ a.e. on $X$. ■
### 3.3 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 $\partial_t\beta(v) \in L^p(0, T; W_0^{1,p}(\Omega))$, $B(\beta(v)) \in L^\infty(0, T; L^1(\Omega))$, $\beta(v) \in C([0, T]; L^2(\Omega)-w)$ and $\partial_t\beta(v) \in L^{p'}(0, T; W^{-1,p'}(\Omega))$. Then $t \in [0, T] \to \int_\Omega B(\beta(v)(\mathbf{x}, t))d\mathbf{x} \in [0, \infty)$ is continuous and, for all $t_1, t_2 \in [0, T]$,

**Lemma 3.6** Let us assume (2b) and (2c). Let $v: \Omega \times (0, T) \to \mathbb{R}$ be measurable such that $\zeta(v) \in L^p(0, T; W_0^{1,p}(\Omega))$, $B(\beta(v)) \in L^\infty(0, T; L^1(\Omega))$, $\beta(v) \in C([0, T]; L^2(\Omega)-w)$ and $\partial_t\beta(v) \in L^{p'}(0, T; W^{-1,p'}(\Omega))$. Then $t \in [0, T] \to \int_\Omega B(\beta(v)(\mathbf{x}, t))d\mathbf{x} \in [0, \infty)$ is continuous and, for all $t_1, t_2 \in [0, T]$,

$$ \int_{t_1}^{t_2} \langle \partial_t \beta(v)(t), \zeta(v(t)) \rangle_{W^{-1,p'}, W_0^{1,p}} dt = \int_{\Omega} B(\beta(v)(\mathbf{x}, t_2)) d\mathbf{x} - \int_{\Omega} B(\beta(v)(\mathbf{x}, t_1)) d\mathbf{x}. \quad (30) $$

**Remark 3.7** Like at the end of Section 1.2, it is important to in mind the separation of $\beta(v(\cdot, \cdot))$ from its continuous representative $\beta(v)(\cdot, \cdot)$.

**Proof.**

Note that we obviously only need to make the proof when $0 \le t_1 < t_2 \le T$.

**Step 1:** truncation, extensions and approximation of $\beta(v)$.

We define $\overline{\beta(v)}: \mathbb{R} \to L^2(\Omega)$ by setting

$$ \overline{\beta(v)}(t) = \begin{cases} \beta(v)(t) & \text{if } t \in [t_1, t_2], \\ \beta(v)(t_1) & \text{if } t \le t_1, \\ \beta(v)(t_2) & \text{if } t \ge t_2. \end{cases} $$

By the continuity property of $\beta(v)$, this definition and gives $\overline{\beta(v)} \in C(\mathbb{R}; L^2(\Omega)-w)$ such that $\partial_t\overline{\beta(v)} = \mathbf{1}_{(t_1, t_2)}\partial_t\beta(v) \in L^{p'}(\mathbb{R}; W^{-1,p'}(\Omega))$ (no Dirac masses have been introduced at $t = t_1$ or $t = t_2$). This regularity of $\partial_t\overline{\beta(v)}$ ensures that the function

$$ t \in \mathbb{R} \mapsto D_h \overline{\beta(v)} := \frac{1}{h} \int_t^{t+h} \partial_t \overline{\beta(v)}(s) ds = \frac{\overline{\beta(v)}(t+h) - \overline{\beta(v)}(t)}{h} \quad \text{in } W^{-1,p'}(\Omega) \quad (31) $$

tend to $\partial_t\overline{\beta(v)}$ in $L^{p'}(\mathbb{R}; W^{-1,p'}(\Omega))$ as $h \to 0$.

**Step 2:** we prove that $||B(\overline{\beta(v)}(t))||_{L^1(\Omega)} \le ||B(\beta(v))||_{L^\infty(0,T;L^1(\Omega))}$ for all $t \in \mathbb{R}$ (not only for a.e. $t$).

Let $t \in [t_1, t_2]$. Since $\beta(v)(\cdot, \cdot) = \beta(v(\cdot, \cdot))$ a.e. on $\Omega \times (t_1, t_2)$, there exists a sequence $t_n \to t$ such that $\beta(v)(\cdot, t_n) = \beta(v(\cdot, t_n))$ in $L^2$ for all $n$ and $||B(\beta(v)(\cdot, t_n))||_{L^1(\Omega)} \le ||B(\beta(v))||_{L^\infty(0,T;L^1(\Omega))}$. As $\beta(v) \in C([0, T]; L^2(\Omega)-w)$, we have $\beta(v)(\cdot, t_n) \to \beta(v)(\cdot, t)$ weakly in $L^2(\Omega)$. We then use the convexity of $B$ and Lemma 3.4 to write, thanks to our choice of $t_n$,

$$ \int_{\Omega} B(\beta(v)(\mathbf{x}, t)) d\mathbf{x} \le \liminf_{n \to \infty} \int_{\Omega} B(\beta(v)(\mathbf{x}, t_n)) d\mathbf{x} \le ||B(\beta(v))||_{L^{\infty}(0, T; L^1(\Omega))} $$

and the proof is complete for $t \in [t_1, t_2]$. The result for $t \le t_1$ or $t \ge t_2$ is obvious since $\overline{\beta(v)}(t)$ is then either $\beta(v)(t_1)$ or $\beta(v)(t_2)$.

**Step 3:** We prove that for all $\tau \in \mathbb{R}$ and a.e. $t \in (t_1, t_2)$,

$$ \langle \overline{\beta(v)}(\tau) - \beta(v)(t), \zeta(v(\cdot, t)) \rangle_{W^{-1,p'}, W_0^{1,p}} \le \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, \tau)) - B(\beta(v)(\mathbf{x}, t)) d\mathbf{x}. \quad (32) $$
Note that if we could just replace the duality product $W^{-1,p'}-W_0^{1,p}$ with an $L^2$ inner product, then this formula would be a straightforward consequence of (27). The problem is that nothing ensures that $\zeta(v)(t) \in L^2(\Omega)$ for a.e. $t$.

We first notice that $\overline{\beta(v)}(\tau) - \beta(v)(t) = \int_t^\tau \partial_t \overline{\beta(v)}(s)ds$ indeed belongs to $W^{-1,p'}(\Omega)$ so the right-hand side of (32) makes sense provided that $t$ is chosen such that $\zeta(v(\cdot, t)) \in W_0^{1,p}(\Omega)$ (which we do from here on). To deal with the fact that $\zeta(v(\cdot, t))$ does not necessarily belong to $L^2(\Omega)$, 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 $\zeta$ we see that there exists $r_k \to +\infty$ as $k \to +\infty$ such that $\zeta(T_k(v(\cdot, t))) = T_{r_k}(\zeta(v(\cdot, t)))$. Hence, $\zeta(T_k(v(\cdot, t)) \in W_0^{1,p}(\Omega)$ and converges, as $k \to \infty$, to $\zeta(v(\cdot, t))$ in $W_0^{1,p}(\Omega)$.

We can therefore write

$$ \begin{aligned} \langle \overline{\beta(v)}(\tau) - \beta(v)(t), \zeta(v)(t) \rangle_{W^{-1,p'}, W_0^{1,p}} &= \lim_{k \to \infty} \langle \overline{\beta(v)}(\tau) - \beta(v)(t), \zeta(T_k(v(\cdot, t))) \rangle_{W^{-1,p'}, W_0^{1,p}} \\ &= \lim_{k \to \infty} \int_{\Omega} (\overline{\beta(v)}(\mathbf{x}, \tau) - \beta(v(\mathbf{x}, t))) \zeta(T_k(v(\mathbf{x}, t))) d\mathbf{x}, \end{aligned} \quad (33) $$

the replacement of the duality product by an $L^2(\Omega)$ inner product being justified since $\overline{\beta(v)}(\tau) - \beta(v)(t)$ and $\zeta(v(\cdot, t))$ both belong to $L^2(\Omega)$. We also used the fact that, for a.e. $t \in (t_1, t_2)$, $\beta(v)(\cdot, t) = \beta(v(\cdot, t))$ a.e. on $\Omega$, so (33) is valid for a.e. $t \in (t_1, t_2)$.

We then write $\beta(v(\mathbf{x}, t)) = \beta(T_k(v(\mathbf{x}, t))) + (\beta(v(\mathbf{x}, t)) - \beta(T_k(v(\mathbf{x}, t))))$ and apply (27) with $r = \overline{\beta(v)}(\mathbf{x}, \tau)$ and $a = T_k(\beta(v(\mathbf{x}, t)))$ to find

$$ \begin{aligned} & \int_{\Omega} (\overline{\beta(v)}(\mathbf{x}, \tau) - \beta(v(\mathbf{x}, t))) \zeta(T_k(v(\mathbf{x}, t))) d\mathbf{x} \\ &= \int_{\Omega} (\overline{\beta(v)}(\mathbf{x}, \tau) - \beta(T_k(v(\mathbf{x}, t)))) \zeta(T_k(v(\mathbf{x}, t))) d\mathbf{x} \\ & \quad - \int_{\Omega} (\beta(v(\mathbf{x}, t)) - \beta(T_k(v(\mathbf{x}, t)))) \zeta(T_k(v(\mathbf{x}, t))) d\mathbf{x} \\ & \le \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, \tau)) - B(\beta(T_k(v(\mathbf{x}, t)))) d\mathbf{x} - \int_{\Omega} (\beta(v(\mathbf{x}, t)) - \beta(T_k(v(\mathbf{x}, t)))) \zeta(T_k(v(\mathbf{x}, t))) d\mathbf{x}. \end{aligned} $$

Studying the cases $v(\mathbf{x}, t) \ge k$ or $v(\mathbf{x}, t) \le -k$, we notice that the last integrand is always non-negative, so we actually can write

$$ \int_{\Omega} (\overline{\beta(v)}(\mathbf{x}, \tau) - \beta(v(\mathbf{x}, t))) \zeta(T_k(v(\mathbf{x}, t))) d\mathbf{x} \le \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, \tau)) - B(\beta(T_k(v(\mathbf{x}, t)))) d\mathbf{x}. $$

We then use the continuity of $B \circ \beta$ and Fatou's lemma to deduce

$$ \begin{aligned} & \limsup_{k \to \infty} \int_{\Omega} (\overline{\beta(v)}(\mathbf{x}, \tau) - \beta(v(\mathbf{x}, t))) \zeta(T_k(v(\mathbf{x}, t))) d\mathbf{x} \\ & \le \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, \tau)) d\mathbf{x} - \liminf_{k \to \infty} \int_{\Omega} B(\beta(T_k(v(\mathbf{x}, t)))) d\mathbf{x} \\ & \le \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, \tau)) d\mathbf{x} - \int_{\Omega} B(\beta(v(\mathbf{x}, t))) d\mathbf{x} \end{aligned} $$

which, combined with (33), concludes the proof of (32) (recall that $t$ has been chosen such that $\beta(v(\cdot, t)) = \beta(v)(\cdot, t)$ a.e. on $\Omega$).

**Step 4:** proof of the formula
By convergence of $D_h \overline{\beta(v)}$ to $\partial_t \overline{\beta(v)}$ in $L^{p'}(0, T; W^{-1,p'}(\Omega))$ and since $\mathbf{1}_{(t_1,t_2)}\zeta(v) \in L^p(\mathbb{R}; W_0^{1,p}(\Omega))$, we have

$$
\begin{align}
\int_{t_1}^{t_2} \langle \partial_t \beta(v)(t), \zeta(v(t)) \rangle_{W^{-1,p'}, W_0^{1,p}} dt &= \int_{\mathbb{R}} \langle \partial_t \overline{\beta(v)}(t), \mathbf{1}_{(t_1,t_2)}(t) \zeta(v(\cdot, t)) \rangle_{W^{-1,p'}, W_0^{1,p}} dt \nonumber \\
&= \lim_{h \to 0} \int_{\mathbb{R}} \langle D_h \overline{\beta(v)}(t), \mathbf{1}_{(t_1,t_2)}(t) \zeta(v(\cdot, t)) \rangle_{W^{-1,p'}, W_0^{1,p}} dt \nonumber \\
&= \lim_{h \to 0} \frac{1}{h} \int_{t_1}^{t_2} \langle \overline{\beta(v)}(s+h) - \overline{\beta(v)}(t), \zeta(v(\cdot, t)) \rangle_{W^{-1,p'}, W_0^{1,p}} dt. \tag{34}
\end{align}
$$

We then use (32) for a.e. $t \in (t_1, t_2)$ to obtain, for $h$ small enough such that $t_1 + h < t_2$,

$$
\begin{align*}
& \frac{1}{h} \int_{t_1}^{t_2} \langle \overline{\beta(v)}(t+h) - \overline{\beta(v)}(t), \zeta(v(\cdot, t)) \rangle_{W^{-1,p'}, W_0^{1,p}} dt \\
& \leq \frac{1}{h} \int_{t_1}^{t_2} \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, t+h)) - B(\overline{\beta(v)}(\mathbf{x}, t)) d\mathbf{x} dt \\
& = \frac{1}{h} \int_{t_2}^{t_2+h} \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, t)) d\mathbf{x} dt - \frac{1}{h} \int_{t_1}^{t_1+h} \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, t)) d\mathbf{x} dt \\
& = \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, t_2)) d\mathbf{x} - \frac{1}{h} \int_{t_1}^{t_1+h} \int_{\Omega} B(\overline{\beta(v)}(\mathbf{x}, t)) d\mathbf{x} dt
\end{align*}
$$

We now take the lim sup of this inequality, using the fact that $B(\beta(v)(\cdot, t_2))$ is integrable (Step 2) to take its integral out of the lim sup. Coming back to (34) we obtain

$$
\int_{t_1}^{t_2} \langle \partial_t \beta(v)(t), \zeta(v(t)) \rangle_{W^{-1,p'}, W_0^{1,p}} dt \\
\leq \int_{\Omega} B(\beta(v)(\mathbf{x}, t_2)) d\mathbf{x} - \liminf_{h \to 0} \frac{1}{h} \int_{t_1}^{t_1+h} \int_{\Omega} B(\beta(v)(\mathbf{x}, t)) d\mathbf{x} dt. \quad (35)
$$

But since $\beta(v) \in C([0, T]; L^2(\Omega)-w)$, as $h \to 0$ we have $\frac{1}{h} \int_{t_1}^{t_1+h} \beta(v)(t)dt \to \beta(v)(t_1)$ weakly in $L^2(\Omega)$. Hence, the convexity of $B$, Lemma 3.4 and Jensen's inequality give

$$
\begin{align*}
\int_{\Omega} B(\beta(v)(\mathbf{x}, t_1)) d\mathbf{x} &\le \liminf_{h \to 0} \int_{\Omega} B \left( \frac{1}{h} \int_{t_1}^{t_1+h} \beta(v)(\mathbf{x}, t) dt \right) d\mathbf{x} \\
&\le \liminf_{h \to 0} \int_{\Omega} \frac{1}{h} \int_{t_1}^{t_1+h} B(\beta(v)(\mathbf{x}, t)) dt d\mathbf{x}.
\end{align*}
$$

Plugged into (35), this inequality shows that (30) holds with $\le$ instead of $=$. The reverse inequality is obtained by reversing the time. We consider $\tilde{v}(t) = v(t_1 + t_2 - t)$. Then $\zeta(\tilde{v})$, $B(\beta(\tilde{v}))$ and $\beta(\tilde{v})$ have the same properties as $\zeta(v)$, $B(\beta(v))$ and $\beta(v)$, and $\beta(\tilde{v})$ takes values $\beta(v)(t_1)$ at $t = t_2$ and $\beta(v)(t_2)$ at $t = t_1$. Applying (30) with "$\le$" instead of "=" to $\tilde{v}$ and using the fact that $\partial_t \beta(\tilde{v})(t) = -\partial_t \beta(v)(t_1 + t_2 - t)$, we obtain (30) with "$\ge$" instead of "=" and the proof of (30) is complete.

The continuity of $t \in [0, T] \mapsto \int_{\Omega} B(\beta(v)(\mathbf{x}, t)) d\mathbf{x}$ is straightforward from (30) as the left-hand side of this relation is continuous with respect to $t_1$ and $t_2$. $\blacksquare$

The following corollary states continuity properties and an essential formula on the solution to (4).
**Corollary 3.8** Under Assumptions (2a)-(2i), if $\bar{u}$ is a solution of (4) then:

1. the function $t \in [0, T] \mapsto \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, t))d\mathbf{x} \in [0, \infty)$ is continuous and bounded,

2. for any $T_0 \in [0, T]$,

$$ \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0))d\mathbf{x} + \int_{0}^{T_0} \int_{\Omega} \mathbf{a}(\mathbf{x}, \nu(\bar{u}(\mathbf{x}, t)), \nabla\zeta(\bar{u})(\mathbf{x}, t)) \cdot \nabla\zeta(\bar{u})(\mathbf{x}, t)d\mathbf{x}dt \\ = \int_{\Omega} B(\beta(u_{\text{ini}}(\mathbf{x})))d\mathbf{x} + \int_{0}^{T_0} \int_{\Omega} f(\mathbf{x}, t)\zeta(\bar{u})(\mathbf{x}, t)d\mathbf{x}dt, \quad (36) $$

3. $\nu(\bar{u})$ is continuous $[0, T] \to L^2(\Omega)$.

**Remark 3.9** The continuity of $\nu(\bar{u})$ has to be understood in the same sense as the continuity of $\beta(\bar{u})$, that is $\nu(\bar{u})$ is a.e. on $\Omega \times (0, T)$ equal to a continuous function $[0, T] \to L^2(\Omega)$. We use in particular the same notation $\nu(\bar{u}(\cdot, \cdot))$ for the continuous representative $\nu(\bar{u})(\cdot, \cdot)$ as we did for the continuous representative of $\beta(\bar{u})$.

**Proof.**

The continuity of $t \in [0, T] \mapsto \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, t))d\mathbf{x} \in [0, \infty)$ and Formula (36) are straightforward consequences of Lemma 3.6 with $v = \bar{u}$ and using (4) with $v = \zeta(\bar{u})$. Note that the bound on $\int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, t))d\mathbf{x}$ 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 $\nu(\bar{u}) : [0, T] \mapsto L^2(\Omega)$. Let $\mathcal{T}$ be the set of $\tau \in [0, T]$ such that $\beta(\bar{u}(\cdot, \tau)) = \beta(\bar{u})(\cdot, \tau)$ a.e. on $\Omega$, and let $(s_l)_{l \in \mathbb{N}}$ and $(t_k)_{k \in \mathbb{N}}$ be two sequences in $\mathcal{T}$ which converge to the same value $s$. Invoking (28) we can write

$$ \int_{\Omega} (\nu(\bar{u}(\mathbf{x}, s_l)) - \nu(\bar{u}(\mathbf{x}, t_k)))^2 d\mathbf{x} \le 4L_{\beta}L_{\zeta} \left( \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, s_l))d\mathbf{x} + \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, t_k))d\mathbf{x} \right) \\ - 8L_{\beta}L_{\zeta} \int_{\Omega} B \left( \frac{\beta(\bar{u})(\mathbf{x}, s_l) + \beta(\bar{u})(\mathbf{x}, t_k)}{2} \right) d\mathbf{x}. \quad (37) $$

Since $\frac{\beta(\bar{u})(\cdot, s_l) + \beta(\bar{u})(\cdot, t_k)}{2} \to \beta(\bar{u})(\cdot, s)$ weakly in $L^2(\Omega)$ as $l, k \to \infty$, Lemma 3.4 gives

$$ \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, s)) d\mathbf{x} \le \liminf_{l,k \to \infty} \int_{\Omega} B \left( \frac{\beta(\bar{u})(\mathbf{x}, s_l) + \beta(\bar{u})(\mathbf{x}, t_k)}{2} \right) d\mathbf{x}. $$

Taking the limsup as $l, k \to \infty$ of (37) and using the continuity of $t \mapsto \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, t))d\mathbf{x}$ thus shows that

$$ \|\nu(\bar{u}(\cdot, s_l)) - \nu(\bar{u}(\cdot, t_k))\|_{L^2(\Omega)} \to 0 \quad \text{as } l, k \to \infty. \quad (38) $$

The existence of an a.e. representative of $\nu(\bar{u}(\cdot, \cdot))$ which is continuous $[0, T] \mapsto L^2(\Omega)$ is a direct consequence of this convergence. Let $s \in [0, T]$ and $(s_l)_{l \in \mathbb{N}} \subset \mathcal{T}$ which converges to $s$. Applied with $t_k = s_k$, (38) shows that $(\nu(\bar{u}(\cdot, s_l)))_{l \in \mathbb{N}}$ is a Cauchy sequence in $L^2(\Omega)$ and therefore that $\lim_{l \to \infty} \nu(\bar{u}(\cdot, s_l))$ exists in $L^2(\Omega)$. Relation (38) moreover also shows that this limit, that we can call $\nu(\bar{u})(\cdot, s)$, does not depend on the Cauchy sequence in $\mathcal{T}$ which converges to $s$. With $t_k = s$, we also see that whenever $s \in \mathcal{T}$ we have $\nu(\bar{u}(\cdot, s)) = \nu(\bar{u})(\cdot, s)$ a.e. on $\Omega$, and $\nu(\bar{u})(\cdot, \cdot)$ is therefore equal to $\nu(\bar{u}(\cdot, \cdot))$ a.e. on $\Omega \times (0, T)$.

It remains to establish that $\nu(\bar{u})$ thus defined is continuous $[0, T] \mapsto L^2(\Omega)$. For any $(\tau_r)_{r \in \mathbb{N}} \subset [0, T]$ which converges to $\tau \in [0, T]$, we can pick $s_r \in \mathcal{T} \cap (\tau_r - \frac{1}{r}, \tau_r + \frac{1}{r})$ and $t_r \in \mathcal{T} \cap (\tau - \frac{1}{r}, \tau + \frac{1}{r})$ such that

$$ \|\nu(\bar{u})(\cdot, \tau_r) - \nu(\bar{u}(\cdot, s_r))\|_{L^2(\Omega)} \le \frac{1}{r}, \quad \|\nu(\bar{u})(\cdot, \tau) - \nu(\bar{u}(\cdot, t_r))\|_{L^2(\Omega)} \le \frac{1}{r}. $$
We therefore have

$$||\nu(\bar{u})(\cdot, \tau_r) - \nu(\bar{u})(\cdot, \tau)||_{L^2(\Omega)} \le \frac{2}{r} + ||\nu(\bar{u}(\cdot, s_r)) - \nu(\bar{u}(\cdot, t_r))||_{L^2(\Omega)}$$

which proves by (38) with $l = k = r$ that $\nu(\bar{u})(\cdot, \tau_r) \to \nu(\bar{u})(\cdot, \tau)$ in $L^2(\Omega)$ as $r \to \infty$ and completes the proof. ■

# 4 Proof of the convergence theorems

## 4.1 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^\infty(0, T; L^2(\Omega))$ estimate and discrete $L^p(0, T; W_0^{1,p}(\Omega))$ estimate) Under Assumptions (2), let $\mathcal{D}$ be a space-time gradient discretisation in the sense of Definition 2.1. Let $u$ be a solution to Scheme (14).

Then, for any $T_0 \in (0, T]$, denoting by $k = 1, \dots, N$ the index such that $T_0 \in (t^{(k-1)}, t^{(k)}]$ we have

$$\int_{\Omega} B(\Pi_{\mathcal{D}}\beta(u)(\mathbf{x}, T_0))d\mathbf{x} + \int_{0}^{T_0} \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}}\nu(u)(\mathbf{x}, t), \nabla_{\mathcal{D}}\zeta(u)(\mathbf{x}, t)) \cdot \nabla_{\mathcal{D}}\zeta(u)(\mathbf{x}, t)d\mathbf{x}dt \\ \le \int_{\Omega} B(\Pi_{\mathcal{D}}\beta(\mathcal{I}_{\mathcal{D}}u_{\text{ini}})(\mathbf{x}))d\mathbf{x} + \int_{0}^{t^{(k)}} \int_{\Omega} f(\mathbf{x}, t)\Pi_{\mathcal{D}}\zeta(u)(\mathbf{x}, t)d\mathbf{x}dt. \quad (39)$$

Consequently, there exists $C_1 > 0$, only depending on $p$, $L_\beta$, $C_P \ge C_\mathcal{D}$ (see Definition 2.5), $C_{\text{ini}} \ge \|\Pi_{\mathcal{D}}\mathcal{I}_{\mathcal{D}}u_{\text{ini}}\|_{L^2(\Omega)}$, $f$ and $\underline{a}$ such that

$$\begin{align} \|\Pi_{\mathcal{D}} B(\beta(u))\|_{L^{\infty}(0,T;L^1(\Omega))} &\le C_1, \quad \|\nabla_{\mathcal{D}} \zeta(u)\|_{L^p(\Omega \times (0,T))} \le C_1 \tag{40} \\ \text{and } \|\Pi_{\mathcal{D}} \beta(u)\|_{L^{\infty}(0,T;L^2(\Omega))} &\le C_1. \nonumber \end{align}$$

**Proof.** Using (12) and (27), we notice that for any $n = 0, \dots, N-1$, any $t \in (t^{(n)}, t^{(n+1)}]$,

$$\begin{align*} \Pi_{\mathcal{D}} \delta_{\mathcal{D}} \beta(u)(t) \Pi_{\mathcal{D}} \zeta(u^{(n+1)}) &= \frac{1}{\dot{\mathcal{A}}^{(n+\frac{1}{2})}} (\beta(\Pi_{\mathcal{D}} u^{(n+1)}) - \beta(u^{(n)})) \zeta(\Pi_{\mathcal{D}} u^{(n+1)}) \\ &\ge \frac{1}{\dot{\mathcal{A}}^{(n+\frac{1}{2})}} (B(\Pi_{\mathcal{D}} \beta(u^{(n+1)})) - B(\Pi_{\mathcal{D}} \beta(u^{(n)}))). \end{align*}$$

Hence, taking $v = (\zeta(u^{(0)}), \zeta(u^{(1)}), \dots, \zeta(u^{(k)}), 0, \dots, 0) \subset X_{\mathcal{D},0}$ in (14), we find

$$\int_{\Omega} B(\Pi_{\mathcal{D}}\beta(u)(\mathbf{x}, t^{(k)}))d\mathbf{x} + \int_{0}^{t^{(k)}} \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}}\nu(u)(\mathbf{x}, t), \nabla_{\mathcal{D}}\zeta(u)(\mathbf{x}, t)) \cdot \nabla_{\mathcal{D}}\zeta(u)(\mathbf{x}, t)d\mathbf{x}dt \\ \le \int_{\Omega} B(\Pi_{\mathcal{D}}\beta(u^{(0)})(\mathbf{x}))d\mathbf{x} + \int_{0}^{t^{(k)}} f(\mathbf{x}, t)\Pi_{\mathcal{D}}\zeta(u)(\mathbf{x}, t)d\mathbf{x}dt \quad (41)$$

Equation (39) is a straightforward consequence of this estimate, of the relation $\beta(u)(\cdot, T_0) = \beta(u)(\cdot, t^{(k)})$ (see (13)) and of the fact that the integrand involving $\mathbf{a}$ is nonnegative on $[T_0, t^{(k)}]$. Using the Young inequality, we can write
$$ \int_0^{t^{(k)}} \int_{\Omega} f(\mathbf{x}, t) \Pi_{\mathcal{D}} \zeta(u)(\mathbf{x}, t) d\mathbf{x} dt \\ \leq \frac{2^{1/(p-1)} C_{\mathcal{D}}^{p'}}{(p\bar{a})^{1/(p-1)}} \frac{\|f\|_{L^{p'}(\Omega \times (0, t^{(k)}))}^{p'}}{p'} + \frac{a}{2C_{\mathcal{D}}^p} \|\Pi_{\mathcal{D}} \zeta(u)\|_{L^p(\Omega \times (0, t^{(k)}))}^p $$

and the first two estimates in (40) therefore follow from (41), the coercivity assumption (2f) on **a** and (26). The estimate on $\Pi_{\mathcal{D}}\beta(u) = \beta(\Pi_{\mathcal{D}}u)$ in $L^{\infty}(0, T; L^2(\Omega))$ is a consequence of the estimate on $B(\beta(\Pi_{\mathcal{D}}u))$ in $L^{\infty}(0, T; L^1(\Omega))$ and of (26). ■

**Corollary 4.2 (Existence of a solution to the Gradient Scheme) Under Assumptions (2), if $\mathcal{D}$ is a gradient discretisation in the sense of Definition 2.1 then there exists at least a solution to the Gradient Scheme (14).**

**Proof.** For $\rho \in [0, 1]$ we let $\beta_{\rho}(u) = \rho u + (1 - \rho)\beta(u)$, $\zeta_{\rho}(u) = \rho u + (1 - \rho)\zeta(u)$ and $a_{\rho}(\mathbf{x}, s, \xi) = \rho\xi + (1 - \rho)a(\mathbf{x}, s, \xi)$. It is clear that $\beta_{\rho}$, $\zeta_{\rho}$ and $a_{\rho}$ satisfy the same assumptions as $\beta$, $\zeta$ and $a$ for some $L_{\beta}$, $M_0$, $M_1$ and $a$ not depending on $\rho$. We can therefore apply Lemma 4.1 to see that there exists $C_2$ not depending on $s$ such that any solution $u_{\rho}$ to (14) with $\beta = \beta_{\rho}$, $\zeta = \zeta_{\rho}$ and $a = a_{\rho}$ satisfies

$$ \|\nabla_{\mathcal{D}} \zeta_{\rho}(u_{\rho})\|_{L^p((0,T) \times \Omega)^d} \leq C_2. $$

Since $\|\nabla_{\mathcal{D}} \cdot \|_{L^p(\Omega \times (0, T))^d}$ is a norm on $X_{\mathcal{D},0}$, this shows that $(\zeta_{\rho}(u_{\rho}))_{\rho \in [0,1]}$ remains bounded in this finite dimensional space, which implies in particular that for all $i \in I$, $(\zeta_{\rho}(u_{\rho})_i)_{\rho \in [0,1]}$. Using Assumption 2b for $\zeta_{\rho}$ with constants not depending on $\rho$, we deduce that $((u_{\rho})_i)_{\rho \in [0,1]}$ remains bounded for any $i \in I$, and thus that $(u_{\rho})_{\rho \in [0,1]}$ is bounded in $X_{\mathcal{D},0}$.

But if $\rho = 0$ then (14) is a linear scheme. Any solution to this scheme being bounded in $X_{\mathcal{D},0}$, this shows that the underlying linear system is invertible. A topological degree argument [10] then shows, combined with the uniform bound on $(u_{\rho})_{\rho \in [0,1]}$, that the scheme corresponding to $\rho = 1$, that is (14), possesses at least one solution. ■

**Lemma 4.3 (Estimate on the dual semi-norm of the discrete time derivative)**

Under Assumptions (2), let $\mathcal{D}$ be a space-time gradient discretisation in the sense of Definition 2.1. Let $u$ be a solution to Scheme (14). Then there exists $C_3$, only depending on $p$, $L_{\beta}$, $C_P \geq C_{\mathcal{D}}$, $C_{\text{ini}} \geq \|\Pi_{\mathcal{D}} I_{\mathcal{D}} u_{\text{ini}}\|_{L^2(\Omega)}$, $f$, $a$, $\mu$, $\bar{a}$ and $T$, such that

$$ \int_0^T |\delta_{\mathcal{D}}\beta(u)(t)|_{\star, \mathcal{D}}^{p'} dt \leq C_3. \qquad (42) $$

**Proof.** Let us take a generic $v = (v^{(n)})_{n=1,...,N} \subset X_{\mathcal{D},0}$ as test function in Scheme (14). We have, thanks to Assumption (2h) on **a**,

$$ \int_0^T \int_{\Omega} \Pi_{\mathcal{D}} \delta_{\mathcal{D}} \beta(u)(\mathbf{x}, t) \Pi_{\mathcal{D}} v(\mathbf{x}, t) d\mathbf{x} dt \leq \int_0^T \int_{\Omega} (\bar{a}(\mathbf{x}) + \mu |\nabla_{\mathcal{D}} \zeta(u)(\mathbf{x}, t)|^{p-1}) |\nabla_{\mathcal{D}} v(\mathbf{x}, t)| d\mathbf{x} dt \\ + \int_0^T \int_{\Omega} f(\mathbf{x}, t) \Pi_{\mathcal{D}} v(\mathbf{x}, t) d\mathbf{x} dt. $$

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_{\beta}$, $C_P$, $C_{\text{ini}}$, $f$, $a$, $\bar{a}$ and $\mu$ such that

$$ \int_0^T \int_{\Omega} \Pi_{\mathcal{D}} \delta_{\mathcal{D}} \beta(u)(\mathbf{x}, t) \Pi_{\mathcal{D}} v(\mathbf{x}, t) d\mathbf{x} dt \leq C_4 \|\nabla_{\mathcal{D}} v\|_{L^p(0, T; L^p(\Omega))^d}. $$
The proof of (42) is completed by selecting $v = (|\delta_D^{(n+\frac{1}{2})}\beta(u)|_{\star, D}^{p'-1}z^{(n)})_{n=0,...,N}$ with $(z^{(n)})_{n=0,...,N} \subset X_{D,0}$ such that, for any $n = 0, ..., N-1$, $z^{(n+1)}$ realises the supremum in (11) with $w = \delta_D^{(n+\frac{1}{2})}\beta(u)$. ■

**Lemma 4.4 (Estimate on the time translates of $\nu(u)$)**

Under Assumptions (2), let $\mathcal{D}$ be a space-time gradient discretisation in the sense of Definition 2.1. Let $u$ be a solution to Scheme (14). Then there exists $C_5$, only depending on $p, L_\beta, L_\zeta, C_P \ge C_\mathcal{D}$, $C_{\text{ini}} \ge \|\Pi_\mathcal{D}I_\mathcal{D}u_{\text{ini}}\|_{L^2(\Omega)}$, $f, \underline{a}, \mu, \bar{a}$ and $T$, such that

$$ \|\Pi_{\mathcal{D}}\nu(u)(\cdot, \cdot + \tau) - \Pi_{\mathcal{D}}\nu(u)(\cdot, \cdot)\|_{L^2(\Omega \times (0, T - \tau))}^2 \le C_5(\tau + \|\dot{a}\|), \quad \forall \tau \in (0, T). \quad (43) $$

**Proof.** Let $\tau \in (0, T)$. Thanks to (24), we get that

$$ \int_{\Omega \times (0, T - \tau)} \left( \Pi_{\mathcal{D}} \nu(u)(\mathbf{x}, t + \tau) - \Pi_{\mathcal{D}} \nu(u)(\mathbf{x}, t) \right)^2 d\mathbf{x}dt \le L_{\beta} L_{\zeta} \int_{0}^{T - \tau} A(t)dt, \quad (44) $$

where, for almost every $t \in (0, T - \tau)$,

$$ A(t) = \int_{\Omega} \left( \Pi_{\mathcal{D}} \zeta(u)(\mathbf{x}, t + \tau) - \Pi_{\mathcal{D}} \zeta(u)(\mathbf{x}, t) \right) \left( \Pi_{\mathcal{D}} \beta(u)(\mathbf{x}, t + \tau) - \Pi_{\mathcal{D}} \beta(u)(\mathbf{x}, t) \right) d\mathbf{x}. $$

Let $t \in (0, T - \tau)$. Letting $n_0(t), n_1(t) = 0, ..., N-1$ be such that $t^{(n_0(t))} \le t < t^{(n_0(t)+1)}$ and $t^{(n_1(t))} \le t + \tau < t^{(n_1(t)+1)}$, we may write

$$ A(t) = \int_{\Omega} \left( \Pi_{\mathcal{D}} \zeta(u^{(n_1(t)+1)})(\mathbf{x}) - \Pi_{\mathcal{D}} \zeta(u^{(n_0(t)+1)})(\mathbf{x}) \right) \left( \sum_{n=n_0(t)+1}^{n_1(t)} \dot{\mathcal{A}}^{(n+\frac{1}{2})} \Pi_{\mathcal{D}} \delta_D^{(n+\frac{1}{2})} \beta(u)(\mathbf{x}) \right) d\mathbf{x}, $$

which also reads

$$ A(t) = \int_{\Omega} \left( \Pi_{\mathcal{D}} \zeta(u^{(n_1(t)+1)})(\mathbf{x}) - \Pi_{\mathcal{D}} \zeta(u^{(n_0(t)+1)})(\mathbf{x}) \right) \times \left( \sum_{n=1}^{N-1} \chi_n(t, t+\tau) \dot{\mathcal{A}}^{(n+\frac{1}{2})} \Pi_{\mathcal{D}} \delta_D^{(n+\frac{1}{2})} \beta(u)(\mathbf{x}) \right) d\mathbf{x}, \quad (45) $$

with $\chi_n(t, t + \tau) = 1$ if $t^{(n)} \in (t, t + \tau]$ and $\chi_n(t, t + \tau) = 0$ if $t^{(n)} \notin (t, t + \tau]$. We then let $v = (\chi_{n-1}(t, t + \tau)(\zeta(u^{(n_1(t)+1)})) - \zeta(u^{(n_0(t)+1)})))_{n=0,...,N}$ in Scheme (14). Using (45), we get

$$ A(t) = A_3(t) - \sum_{n=1}^{N-1} \chi_n(t, t + \tau) \dot{\mathcal{A}}^{(n+\frac{1}{2})} \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}}\nu(u^{(n+1)}), \nabla_{\mathcal{D}}\zeta(u^{(n+1)}))(\mathbf{x}) \cdot (\nabla_{\mathcal{D}}\zeta(u^{(n_1(t)+1)})(\mathbf{x}) - \nabla_{\mathcal{D}}\zeta(u^{(n_0(t)+1)})(\mathbf{x})) d\mathbf{x} $$

with

$$ A_3(t) = \sum_{n=1}^{N-1} \chi_n(t, t + \tau) \int_{\Omega} \int_{t^{(n)}}^{t^{(n+1)}} f(\mathbf{x}, t) dt \left( \Pi_{\mathcal{D}}\zeta(u^{(n_1(t)+1)})(\mathbf{x}) - \Pi_{\mathcal{D}}\zeta(u^{(n_0(t)+1)})(\mathbf{x}) \right) d\mathbf{x}dt. $$
Using the Young inequality and the inequality $a_1a_2 \le \frac{(p-1)^{1/p'}}{p}(a_1^p + a_2^{p'})$ for $a_1 \ge 0$ and $a_2 \ge 0$, this yields:

$$A(t) \le \frac{(p-1)^{1/p'}}{p} (A_0(t) + A_1(t) + 2A_2(t)) + A_3(t), \quad (46)$$

with

$$A_0(t) = \sum_{n=1}^{N-1} \chi_n(t, t + \tau) \partial t^{(n+\frac{1}{2})} \int_{\Omega} |\nabla_D \zeta(u^{(n_0(t)+1)})(x)|^p dx,$$

$$A_1(t) = \sum_{n=1}^{N-1} \chi_n(t, t + \tau) \partial t^{(n+\frac{1}{2})} \int_{\Omega} |\nabla_D \zeta(u^{(n_1(t)+1)})(x)|^p dx,$$

and

$$A_2(t) = \sum_{n=1}^{N-1} \chi_n(t, t + \tau) \partial t^{(n+\frac{1}{2})} \int_{\Omega} |\mathbf{a}(x, \Pi_D \nu(u^{(n+1)}), \nabla_D \zeta(u^{(n+1)}))(x)|^{p'} dx.$$

Applying (80) of Lemma 6.6 yields

$$\int_{0}^{T-\tau} A_0(t) + A_1(t) dt \le 2(\tau + \partial t) \|\nabla_D \zeta(u)\|_{L^p(\Omega \times (0,T))}^p. \quad (47)$$

Applying (79) of Lemma 6.6 yields

$$\int_{0}^{T-\tau} A_2(t) dt \le \tau \|\mathbf{a}(\cdot, \Pi_D \nu(u), \nabla_D \zeta(u))\|_{L^{p'}(\Omega \times (0,T))}^{p'} \quad (48)$$

and, similarly, we obtain using Definition (2.5) and (40), the existence of $C_6$ such that

$$\int_{0}^{T-\tau} A_3(t) dt \le C_6 ((\tau + \partial t) + \tau \|f\|_{L^{p'}(\Omega \times (0,T))}^{p'}). \quad (49)$$

Using inequalities (44), (46) and (47)-(49), (43) is proved. ■

## 4.2 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 $(\Pi_{D_m}\beta(u_m))_{m \in \mathbb{N}}$ converges weakly in $L^2(\Omega)$ uniformly in $[0, T]$ (in the sense of Definition 6.4) to some function $\bar{\beta} \in C([0, T]; L^2(\Omega)-w)$ which satisfies $\bar{\beta}(\cdot, 0) = \beta(u_{\text{ini}})$ in $L^2(\Omega)$. Using again Estimates (40) and applying Lemma 2.9, we again extract a subsequence such that, for some $\bar{\zeta} \in L^p(0, T; W_0^{1,p}(\Omega))$, $\Pi_{D_m}\zeta(u_m) \to \bar{\zeta}$ weakly in $L^p(\Omega \times (0, T))$ and $\nabla_{D_m}\zeta(u_m) \to \bar{\nabla}\zeta$ weakly in $L^p(\Omega \times (0, T))^d$. From Estimates (40), Definition 2.5 and the growth assumption (2b) on $\zeta$, we also see that $(\Pi_{D_m}u_m)_{m \in \mathbb{N}}$ is bounded in $L^p(\Omega \times (0, T))$ and we can therefore assume, up to a subsequence, that it converges weakly to some $\bar{u}$ in this space.

We then remark that $|\nu(a) - \nu(b)| \le L_{\beta}|\zeta(a) - \zeta(b)|$, which implies that, using Definition 2.7 with $v = \zeta(u_m)$ and (40),

$$\|\Pi_{D_m} \nu(u_m)(\cdot + \xi) - \Pi_{D_m} \nu(u)\|_{L^p(\mathbb{R}^d \times (0,T))} \le L_{\beta} C_1 T_{D_m}(\xi) \quad (50)$$

where $\Pi_D\nu(u)$ has been extended by 0 outside $\Omega$ and $\lim_{\xi \to 0} \sup_{m \in \mathbb{N}} \sup_{T_{D_m}(\xi)} = 0$. Invoking Lemma 4.4, we also see that the time translates of $\Pi_{D_m}\nu(u_m)$ tend uniformly to 0 in $L^1(\Omega \times (0, T))$. We
use the fact that $\Pi_{\mathcal{D}_m}\beta(u_m)$, and therefore also $\Pi_{\mathcal{D}_m}\nu(u_m)$, remains bounded in $L^\infty(0, T; L^2(\Omega))$ 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, $\Pi_{\mathcal{D}_m}\nu(u_m) \to \bar{\nu}$ in $L^1(\Omega \times (0, T))$. Under the first case in the structural hypothesis (17), we have $\beta = \text{Id}$ and therefore $\bar{\beta} = \bar{u} = \beta(\bar{u})$, and $\nu = \zeta$. The strong convergence of $\Pi_{\mathcal{D}_m}\nu(u_m) = \Pi_{\mathcal{D}_m}\zeta(u_m)$ to $\bar{\nu} = \bar{\zeta}$ thus allows us to apply Lemma 3.5 to see that $\bar{\zeta} = \zeta(\bar{u})$. Exchanging the roles of $\beta$ and $\zeta$, we see that $\bar{\beta} = \beta(\bar{u})$ and $\bar{\zeta} = \zeta(\bar{u})$ still hold in the second case of (17). We notice that this is the only place where we need this structural assumption (17) on $\beta, \zeta$.

Upon extraction of another subsequence, we can also assume that $\mathbf{a}(\cdot, \Pi_{\mathcal{D}}\nu(u), \nabla_{\mathcal{D}}\zeta(u))$ has a weak limit in $L^p'(\Omega \times (0, T))^d$, which we denote by $\mathbf{A}$.

Finally, for any $T_0 \in [0, T]$, since $\Pi_{\mathcal{D}_m}\beta(u_m(\cdot, T_0)) \to \beta(\bar{u})(\cdot, T_0)$ weakly in $L^2(\Omega)$, Lemma 3.4 gives

$$ \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0)) d\mathbf{x} \leq \liminf_{m \to \infty} \int_{\Omega} B(\beta(\Pi_{\mathcal{D}_m} u_m)(\mathbf{x}, T_0)) d\mathbf{x}. \quad (51) $$

This shows in particular that $B(\beta(\bar{u})) \in L^\infty(0, T; L^1(\Omega))$.

**Step 2** Passing to the limit in the scheme. We drop the indices $m$ for legibility reasons.

Let $\varphi \in C_c^1(-\infty, T)$ and $w \in W_0^{1,p}(\Omega) \cap L^2(\Omega)$. We introduce $v = (\varphi(t^{(n-1)}), P_D w)_n$ with $v = (0, \dots, 0)$ (where $t^{(n-1)} = 0$ for example, this value being irrelevant) as test function in (14), with $P_D$ defined by (21).

We get $T_1^{(m)} + T_2^{(m)} = T_3^{(m)}$ with

$$ T_1^{(m)} = \sum_{n=0}^{N-1} \varphi(t^{(n)}) \partial t^{(n+\frac{1}{2})} \int_{\Omega} \Pi_{\mathcal{D}} \delta_{\mathcal{D}}^{(n+\frac{1}{2})} \beta(u)(\mathbf{x}) \Pi_{\mathcal{D}} P_{\mathcal{D}} w(\mathbf{x}) d\mathbf{x}, $$

$$ T_2^{(m)} = \sum_{n=0}^{N-1} \varphi(t^{(n)}) \partial t^{(n+\frac{1}{2})} \int_{\Omega} \mathbf{a} \left( \mathbf{x}, \Pi_{\mathcal{D}}\nu(u^{(n+1)}), \nabla_{\mathcal{D}}\zeta(u^{(n+1)})(\mathbf{x}) \right) \cdot \nabla_{\mathcal{D}} P_{\mathcal{D}} w(\mathbf{x}) d\mathbf{x}, $$

and

$$ T_3^{(m)} = \sum_{n=0}^{N-1} \varphi(t^{(n)}) \int_{t^{(n)}}^{t^{(n+1)}} \int_{\Omega} f(\mathbf{x}, t) \Pi_{\mathcal{D}} P_{\mathcal{D}} w(\mathbf{x}) d\mathbf{x} dt. $$

Using discrete integrate-by-parts to transform the terms $\varphi(t^{(n)})(\Pi_{\mathcal{D}}\beta(u^{(n+1)})) - \Pi_{\mathcal{D}}\beta(u^{(n)})$ appearing in $T_1^{(m)}$ into $(\varphi(t^{(n)}) - \varphi(t^{(n+1)}))\Pi_{\mathcal{D}}\beta(u^{(n+1)})$, we have

$$ T_1^{(m)} = - \int_0^T \varphi'(t) \int_{\Omega} \Pi_{\mathcal{D}} \beta(u)(\mathbf{x}, t) \Pi_{\mathcal{D}} P_{\mathcal{D}} w(\mathbf{x}) d\mathbf{x} dt - \varphi(0) \int_{\Omega} \Pi_{\mathcal{D}} \beta(u^{(0)})(\mathbf{x}) \Pi_{\mathcal{D}} P_{\mathcal{D}} w(\mathbf{x}) d\mathbf{x}. $$

Setting $\varphi_{\mathcal{D}}(t) = \varphi(t^{(n)})$ for $t \in (t^{(n)}, t^{(n+1)})$, we have

$$ T_2^{(m)} = \int_0^T \varphi_{\mathcal{D}}(t) \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}}\nu(u)(\mathbf{x}, t), \nabla_{\mathcal{D}}\zeta(u)(\mathbf{x}, t)) \cdot \nabla_{\mathcal{D}} P_{\mathcal{D}} w(\mathbf{x}) d\mathbf{x} dt $$

$$ T_3^{(m)} = \int_0^T \varphi_{\mathcal{D}}(t) \int_{\Omega} f(\mathbf{x}, t) \Pi_{\mathcal{D}} P_{\mathcal{D}} w(\mathbf{x}) d\mathbf{x} dt. $$

We may thus let $m \to \infty$ in $T_1^{(m)} + T_2^{(m)} = T_3^{(m)}$ to see that $\bar{u}$ satisfies

$$ \left\{ \begin{array}{l} \bar{u} \in L^p(\Omega \times (0, T)), \zeta(\bar{u}) \in L^p(0, T; W_0^{1,p}(\Omega)), B(\beta(\bar{u})) \in L^\infty(0, T; L^1(\Omega)), \\ - \int_0^T \varphi'(t) \int_{\Omega} \beta(\bar{u}(\mathbf{x}, t)) w(\mathbf{x}) d\mathbf{x} dt - \varphi(0) \int_{\Omega} u_{\text{ini}}(\mathbf{x}) w(\mathbf{x}) d\mathbf{x} \\ + \int_0^T \varphi(t) \int_{\Omega} \mathbf{A}(\mathbf{x}, t) \cdot \nabla w(\mathbf{x}) d\mathbf{x} dt = \int_0^T \varphi(t) \int_{\Omega} f(\mathbf{x}, t) w(\mathbf{x}) d\mathbf{x} dt, \\ \forall w \in W_0^{1,p}(\Omega) \cap L^2(\Omega), \forall \varphi \in C_c^\infty(-\infty, T). \end{array} \right. \quad (52) $$
Linear combinations of this relation show that it also holds with $\varphi(t)w(x)$ replaced by a tensorial functions in $C_c^\infty(\Omega \times (0, T))$. This proves that $\partial_t\beta(\bar{u}) \in L^{p'}(0, T; W^{-1,p'}(\Omega))$. Standard arguments then show that $\beta(\bar{u})$ can be identified with an element of $C^0([0, T]; L^2(\Omega)-w)$ with the property $\beta(\bar{u}(\cdot, 0)) = \beta(u_{\text{ini}})$ (cf. discussion at the end of Section 1.2). Using the density of tensorial functions in $L^p(0, T; W_0^{1,p}(\Omega))$, we then see that $\bar{u}$ satisfies

$$ \int_0^T \langle \partial_t \beta(\bar{u})(\cdot, t), \bar{v}(\cdot, t) \rangle_{(W^{-1,p'}(\Omega))' W_0^{1,p}(\Omega)} dt + \int_0^T \int_{\Omega} \mathbf{A}(\mathbf{x}, t) \cdot \nabla \bar{v}(\mathbf{x}, t) d\mathbf{x} dt = \int_0^T \int_{\Omega} f(\mathbf{x}, t) \bar{v}(\mathbf{x}, t) d\mathbf{x} dt, \quad \forall \bar{v} \in L^p(0, T; W_0^{1,p}(\Omega)). \quad (53) $$

**Step 3** Proof that $\bar{u}$ is solution to (4).

The proof will be completed by showing that

$$ \mathbf{A}(\mathbf{x}, t) = \mathbf{a}(\mathbf{x}, \nu(\bar{u})(\mathbf{x}, t), \nabla\zeta(\bar{u})(\mathbf{x}, t)), \text{ for a.e. } (\mathbf{x}, t) \in \Omega \times (0, T). \quad (54) $$

We take $T_0 \in [0, T]$, write (39) with $\mathcal{D} = \mathcal{D}_m$ and take the supremum limit as $m \to \infty$. We notice that the $t^{(k)} = T_m$ from Lemma 4.1 converges to $T_0$ as $m \to \infty$. Hence, using the quadratic growth of $B$ we obtain

$$ \limsup_{m \to \infty} \int_0^{T_0} \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}_m} \nu(u_m)(\mathbf{x}, t), \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t)) \cdot \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t) d\mathbf{x} dt \\ \le \int_{\Omega} B(\beta(u_{\text{ini}})(\mathbf{x})) d\mathbf{x} + \int_0^{T_0} \int_{\Omega} f(\mathbf{x}, t) \zeta(\bar{u})(\mathbf{x}, t) d\mathbf{x} dt \\ - \liminf_{m \to \infty} \int_{\Omega} B(\beta(\Pi_{\mathcal{D}_m} u_m)(\mathbf{x}, T_0)) d\mathbf{x}. \quad (55) $$

We then apply Lemma 3.6 and take $\bar{v} = \zeta(\bar{u})$ in (53) to get

$$ \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0)) d\mathbf{x} - \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, 0)) d\mathbf{x} \\ + \int_0^{T_0} \int_{\Omega} \mathbf{A}(\mathbf{x}, t) \cdot \nabla \zeta(\bar{u})(\mathbf{x}, t) d\mathbf{x} dt = \int_0^{T_0} \int_{\Omega} f(\mathbf{x}, t) \zeta(\bar{u})(\mathbf{x}, t) d\mathbf{x} dt. $$

This relation, combined with (55) and using (51), shows that

$$ \limsup_{m \to \infty} \int_0^{T_0} \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}_m} \nu(u_m)(\mathbf{x}, t), \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t)) \cdot \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t) d\mathbf{x} dt \\ \le \int_0^{T_0} \int_{\Omega} \mathbf{A}(\mathbf{x}, t) \cdot \nabla \zeta(\bar{u})(\mathbf{x}, t) d\mathbf{x} dt. \quad (56) $$

It is now possible to apply Minty's trick. Consider, for $\mathbf{G} \in L^p(0, T; L^p(\Omega))^d$,

$$ \int_0^{T_0} \int_{\Omega} [\mathbf{a}(\cdot, \Pi_{\mathcal{D}}\nu(u), \nabla_{\mathcal{D}}\zeta(u)) - \mathbf{a}(\cdot, \Pi_{\mathcal{D}}\nu(u), \mathbf{G})] \cdot [\nabla_{\mathcal{D}}\zeta(u) - \mathbf{G}] d\mathbf{x} dt \ge 0. \quad (57) $$

By strong convergence of $\Pi_{\mathcal{D}_m} \nu(u_m)$ to $\nu(\bar{u})$ in $L^1(\Omega \times (0, T))$ and Assumptions (2e), (2h) on $\mathbf{a}$, we see that $\mathbf{a}(\cdot, \Pi_{\mathcal{D}_m} \nu(u_m), \mathbf{G}) \to \mathbf{a}(\cdot, \nu(\bar{u}), \mathbf{G})$ strongly in $L^p'(\Omega \times (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 \in L^p(0, T; L^p(\Omega))^d$,

$$ \int_0^T \int_{\Omega} [A(\mathbf{x}, t) - a(\mathbf{x}, \nu(\bar{u})(\mathbf{x}, t), G(\mathbf{x}, t))] \cdot [\nabla \zeta(\bar{u})(\mathbf{x}, t) - G(\mathbf{x}, t)] \, d\mathbf{x} dt \geq 0. $$

Application of Minty's method [33] (i.e. taking $G = \nabla\zeta(\bar{u}) + r\varphi$ for $\varphi \in L^p(0, T; L^p(\Omega))^d$ and letting $r \to 0$) then shows that (54) holds and concludes the proof that $\bar{u}$ is a weak solution to (4).

## 4.3 Proof of Theorem 2.14

Let $T_0 \in [0, T]$ and $(T_m)_{m \ge 1}$ be a sequence in $[0, T]$ which converges to $T_0$. By setting $T_0 = T_m$ and $G = \nabla\zeta(\bar{u})$ in the developed form of (57), by taking the infimum limit (thanks to the strong convergence of $a(\cdot, \Pi_{D_m}\nu(u_m), \nabla\zeta(\bar{u}))$) and by using (54), we find

$$ \begin{aligned} & \liminf_{m \to \infty} \int_0^{T_m} \int_{\Omega} a(\mathbf{x}, \Pi_{D_m} \nu(u_m)(\mathbf{x}, t), \nabla_{D_m} \zeta(u_m)(\mathbf{x}, t)) \cdot \nabla_{D_m} \zeta(u_m)(\mathbf{x}, t) d\mathbf{x} dt \\ & \ge \int_0^{T_0} \int_{\Omega} a(\mathbf{x}, \nu(\bar{u})(\mathbf{x}, t), \nabla\zeta(\bar{u})(\mathbf{x}, t)) \cdot \nabla\zeta(\bar{u})(\mathbf{x}, t) d\mathbf{x} dt. \end{aligned} \quad (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 \in (t^{(k-1)}, t^{(k)}]$ converges to $T_0$ as $m \to \infty$. Using (58) and Corollary 3.8, we therefore obtain

$$ \limsup_{m \to \infty} \int_{\Omega} B(\beta(\Pi_{D_m} u_m(\mathbf{x}, T_m))) d\mathbf{x} \le \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0)) d\mathbf{x}. \quad (59) $$

By Lemma 6.3, the uniform-in-time weak convergence of $\beta(\Pi_{D_m} u_m)$ to $\beta(\bar{u})$ and the weak continuity of $\beta(\bar{u}) : [0, T] \to L^2(\Omega)$, we have $\beta(\Pi_{D_m} u_m)(T_m) \to \beta(\bar{u})(T_0)$ weakly in $L^2(\Omega)$ as $m \to \infty$. Therefore, for any $(s_m)_{m \in \mathbb{N}}$ converging to $T_0$, $\frac{1}{2}(\beta(\Pi_{D_m} u_m(T_m)) + \beta(\bar{u})(s_m)) \to \beta(\bar{u})(T_0)$ weakly in $L^2(\Omega)$ as $m \to \infty$ and Lemma 3.4 gives, by convexity of $B$,

$$ \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0)) d\mathbf{x} \le \liminf_{m \to \infty} \int_{\Omega} B \left( \frac{\beta(\Pi_{D_m} u_m(\mathbf{x}, T_m)) + \beta(\bar{u})(\mathbf{x}, s_m)}{2} \right) d\mathbf{x}. \quad (60) $$

Property (28) of $B$ and the two inequalities (59) and (60) allow us to conclude the proof. Let $(s_m)_{m \in \mathbb{N}}$ be a sequence in $\mathcal{T}$ (see the proof of Corollary 3.8) which converges to $T_0$. Then $\nu(\bar{u}(\cdot, s_m)) \to \nu(\bar{u})(\cdot, T_0)$ in $L^2(\Omega)$ as $m \to \infty$. Using (28), we get

$$ \begin{aligned} & \|\nu(\Pi_{D_m} u_m(\cdot, T_m)) - \nu(\bar{u})(\cdot, T_0)\|_{L^2(\Omega)}^2 \\ & \le 2\|\nu(\Pi_{D_m} u_m(\cdot, T_m)) - \nu(\bar{u}(\cdot, s_m))\|_{L^2(\Omega)}^2 + 2\|\nu(\bar{u}(\cdot, s_m)) - \nu(\bar{u})(\cdot, T_0)\|_{L^2(\Omega)}^2 \\ & \le 8L_{\beta}L_{\zeta} \int_{\Omega} [B(\beta(\Pi_{D_m} u_m(\mathbf{x}, T_m))) + B(\beta(\bar{u}(\mathbf{x}, s_m)))] d\mathbf{x} \\ & \quad - 16L_{\beta}L_{\zeta} \int_{\Omega} B \left( \frac{\beta(\Pi_{D_m} u_m(\mathbf{x}, T_m)) + \beta(\bar{u}(\mathbf{x}, s_m))}{2} \right) d\mathbf{x} \\ & \quad + 2\|\nu(\bar{u}(\cdot, s_m)) - \nu(\bar{u})(\cdot, T_0)\|_{L^2(\Omega)}^2. \end{aligned} $$

We then take the $\limsup$ as $m \to \infty$ of this expression. Thanks to (59) and the continuity of $t \in [0, T] \mapsto \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, t)) d\mathbf{x} \in [0, \infty)$ (see Corollary 3.8), the first term in the right-hand side
has a finite lim sup, bounded above by $16L_{\beta}L_{\zeta} \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0))d\mathbf{x}$. We can therefore split the lim sup of this right-hand side without risking writing $\infty - \infty$ and we get, thanks to (60),

$$ \limsup_{m \to \infty} \|\nu(\Pi_{\mathcal{D}_m} u_m(\cdot, T_m)) - \nu(\bar{u})(\cdot, T_0)\|_{L^2(\Omega)}^2 \le 0. $$

Thus, $\nu(\Pi_{\mathcal{D}_m} u_m(\cdot, T_m)) \to \nu(\bar{u})(T_0)$ strongly in $L^2(\Omega)$. By Lemma 6.3 and the continuity of $\nu(\bar{u})$ stated in Corollary 3.8, this concludes the proof of the convergence of $\nu(\Pi_{\mathcal{D}_m} u_m)$ to $\nu(\bar{u})$ in $L^{\infty}(0, T; L^2(\Omega))$.

**Remark 4.5** Since $\beta(\Pi_{\mathcal{D}_m} u_m)(T_m) \to \beta(\bar{u})(T_0)$ weakly in $L^2(\Omega)$ as $m \to \infty$, Lemma 3.4 also shows that $\int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0))d\mathbf{x} \le \liminf_{m \to \infty} \int_{\Omega} B(\beta(\Pi_{\mathcal{D}_m} u_m)(\mathbf{x}, T_m))d\mathbf{x}$, and therefore, combined with (59), that

$$ \lim_{m \to \infty} \int_{\Omega} B(\beta(\Pi_{\mathcal{D}_m} u_m(\mathbf{x}, T_m)))d\mathbf{x} = \int_{\Omega} B(\beta(\bar{u})(\mathbf{x}, T_0))d\mathbf{x}. \quad (61) $$

## 4.4 Proof of Theorem 2.16

Writing (39) for $u_m$ with $T_0 = T$, taking the lim sup as $m \to \infty$, using (61) (with $T_m = T$) and the continuous integration-by-part formula (36) we find that

$$ \begin{aligned} & \limsup_{m \to \infty} \int_0^T \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}_m} \nu(u_m)(\mathbf{x}, t), \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t)) \cdot \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t) d\mathbf{x} dt \\ & \le \int_0^{T_0} \int_{\Omega} \mathbf{a}(\mathbf{x}, \nu(\bar{u})(\mathbf{x}, t), \nabla \zeta(\bar{u})(\mathbf{x}, t)) \cdot \nabla \zeta(\bar{u})(\mathbf{x}, t) d\mathbf{x} dt. \end{aligned} $$

Combined with (58), this shows that

$$ \begin{aligned} & \lim_{m \to \infty} \int_0^T \int_{\Omega} \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}_m} \nu(u_m)(\mathbf{x}, t), \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t)) \cdot \nabla_{\mathcal{D}_m} \zeta(u_m)(\mathbf{x}, t) d\mathbf{x} dt \\ & = \int_0^{T_0} \int_{\Omega} \mathbf{a}(\mathbf{x}, \nu(\bar{u})(\mathbf{x}, t), \nabla \zeta(\bar{u})(\mathbf{x}, t)) \cdot \nabla \zeta(\bar{u})(\mathbf{x}, t) d\mathbf{x} dt. \end{aligned} \quad (62) $$

Let us define

$$ f_m = [\mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}_m} \nu(u_m), \nabla_{\mathcal{D}_m} \zeta(u_m)) - \mathbf{a}(\mathbf{x}, \Pi_{\mathcal{D}_m} \nu(u_m)(\cdot, t), \nabla \zeta(\bar{u}))] \cdot [\nabla_{\mathcal{D}_m} \zeta(u_m) - \nabla \zeta(\bar{u})] \ge 0. $$

Developing this expression and using (62) and (18), we see that $\int_0^T \int_{\Omega} f_m(\mathbf{x}, t)d\mathbf{x}dt \to 0$ as $m \to \infty$. This shows that $f_m \to 0$ in $L^1(\Omega \times (0, T))$ and therefore a.e. up to a subsequence. We can then reason as in [16], using the strict monotony (19) of $\mathbf{a}$, the coercivity assumption (2f) and Vitali's theorem, to deduce that $\nabla_{\mathcal{D}_m} \zeta(u_m) \to \nabla \zeta(\bar{u})$ strongly in $L^p(\Omega \times (0, T))^d$ as $m \to \infty$.

# 5 About Assumption (17)

We discuss show here that, if $p \ge 2$, Theorem 2.10 can be proved without the structural assumption (17) – i.e. without assuming that $\beta = \text{Id}$ or $\zeta = \text{Id}$.

We first notice that we can always assume that

$$ \beta + \zeta \text{ is strictly increasing.} \quad (63) $$
Indeed, if this is not the case then we have $s_1 < s_2$ such that $(\beta + \zeta)(s_1) = (\beta + \zeta)(s_2)$ and, since $\beta$ and $\zeta$ are non-decreasing, that $[s_1, s_2]$ is a common plateau of $\beta$ and $\zeta$. Denoting by $\tilde{\beta}$, $\tilde{\zeta}$ and $\tilde{\nu}$ the functions obtained from $\beta$, $\zeta$ and $\nu$ by removing this common plateau (by a contraction of the $s$-ordinate), we see that $u$ is a solution to (1) if and only if $u$ is a solution of the same problem with $\beta$, $\zeta$ and $\nu$ replaced with $\tilde{\beta}$, $\tilde{\zeta}$ and $\tilde{\nu}$.

**Remark 5.1** Note that a rigorous and global way to remove all common plateaux of $\beta$ and $\zeta$ at once is to actually consider, for $f = \beta$, $\zeta$ or $\nu$, $\tilde{f} = f \circ (\beta + \zeta) \circ (\beta + \zeta)_r$, where $(\beta + \zeta)_r$ is the pseudo-inverse of $\alpha + \beta$ constructed as in (3).

Once this reduction has been applied, we can state the following theorem.

**Theorem 5.2** Under Assumptions (2), let $(\mathcal{D}_m)_{mp \in \mathbb{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 \in \mathbb{N}$, $u_m$ be the solution to (14) with $\mathcal{D} = \mathcal{D}_m$, provided by Theorem 2.10.

If we assume that Assumption (63) hold and that $p \ge 2$, then there exists a solution $\bar{u}$ to (4) such that, up to a subsequence:

* the convergences stated in (18) hold,

* $\Pi_{\mathcal{D}_m} \nu(u_m) \to \nu(\bar{u})$ strongly in $L^\infty(0, T; L^2(\Omega))$ as $m \to \infty$,

* under the strict monotony of $\mathbf{a}$ (i.e. Assumption (19)), as $m \to \infty$, $\Pi_{\mathcal{D}_m} \zeta(u_m) \to \zeta(\bar{u})$ strongly in $L^p(\Omega \times (0, T))$ and $\nabla_{\mathcal{D}_m} \zeta(u_m) \to \nabla \zeta(\bar{u})$ strongly in $L^p(\Omega \times (0, T))^d$.

**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 (17) and the assumption that $\beta + \zeta$ is strictly increasing. The only place in the proof of Theorem (2.10) where we used Assumption (17) is in Step 1, to identify the limits $\bar{\beta}$, $\bar{\zeta}$ and $\bar{\nu}$ of $\Pi_{\mathcal{D}_m} \beta(u_m)$, $\Pi_{\mathcal{D}_m} \zeta(u_m)$ and $\Pi_{\mathcal{D}_m} \nu(u_m)$. We will show that without assuming (17) but by assuming (63), we can still identify those limits.

We first work with the first two limits. Lemmas 4.1 and 4.3 show that $\beta_m = \beta(u_m)$ and $\zeta_m = \zeta(u_m)$ satisfy the assumptions of the discrete compensated compactness theorem 5.4 below (we also use that $p \ge 2$ here). Hence, $\Pi_{\mathcal{D}_m} \beta(u_m) \Pi_{\mathcal{D}_m} \zeta(u_m) \to \bar{\beta} \bar{\zeta}$ in the sense of measures on $\Omega \times (0, T)$. We can then apply Lemma 5.6 with $\varphi \equiv 1$ and $w_m = \Pi_{\mathcal{D}_m} u_m$. This gives a measurable $\bar{u}$ such that $\bar{\beta} = \beta(\bar{u})$ and $\bar{\zeta} = \zeta(\bar{u})$ a.e. on $\Omega \times (0, T)$.

We now turn to $\bar{\nu}$. Estimates (50) and (43) and Kolmogorov's compactness theorem show that the convergence of $\Pi_{\mathcal{D}_m} \nu(u_m)$ towards $\bar{\nu}$ is actually strong on $L^2(\Omega \times (0, T))$ (we use $p \ge 2$ here). Setting $\mu = \beta + \zeta$, we saw that $\mu(\Pi_{\mathcal{D}_m} u_m) \to \bar{\mu} := \beta(\bar{u}) + \zeta(\bar{u})$ weakly in $L^2(\Omega \times (0, T))$. We can therefore apply Lemma 5.6 with $\varphi \equiv 1$ and $\mu, \nu$ instead of $\beta, \zeta$, and we find a measurable $w$ such that $\bar{\nu} = \nu(w)$ and $\bar{\mu} = \mu(w)$. The second relation translates into $(\beta + \zeta)(\bar{u}) = (\beta + \zeta)(w)$, that is $w = \bar{u}$ since $\beta + \zeta$ is strictly increasing. Hence, $\bar{\nu} = \nu(\bar{u})$.

To summarise, the limits of $\Pi_{\mathcal{D}_m} \beta(u_m)$, $\Pi_{\mathcal{D}_m} \zeta(u_m)$ and $\Pi_{\mathcal{D}_m} \nu(u_m)$ have been identified as $\beta(\bar{u})$, $\zeta(\bar{u})$ and $\nu(\bar{u})$ for some $\bar{u}$. This allows to take over the proof of Theorem 2.10 from after the usage of (17) and conclude that $\bar{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. ■
**Remark 5.3** Note that it is not proved, in this context, that $\bar{u}$ is a weak limit of $\Pi_{\mathcal{D}_m} u_m$, but such a limit is not stated in (18) 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 $\beta(\bar{u})$, $\zeta(\bar{u})$ and $\nu(\bar{u})$.

We now state the two key results that allowed us to replace Assumption (17) by Assumption (17). The first is a discrete version of a compensated compactness result in [31]. 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 $\varphi \equiv 1$), but we state this more general version nevertheless as it is the genuine discrete equivalent of the compensated compactness result in [31].

**Theorem 5.4 (Discrete compensated compactness)** We take $T > 0$, $p \ge 2$ and a sequence $(\mathcal{D}_m)_{m \in \mathbb{N}} = (X_{\mathcal{D}_m,0}, \Pi_{\mathcal{D}_m}, \nabla_{\mathcal{D}_m}, \mathcal{I}_{\mathcal{D}_m}, (t_m^{(n)})_{n=0,...,N_m})_{m \in \mathbb{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 \in \mathbb{N}$, let $\beta_m = (\beta_m^{(n)})_{n=0,...,N_m} \subset X_{\mathcal{D}_m,0}$ and $\zeta_m = (\zeta_m^{(n)})_{n=0,...,N_m} \subset X_{\mathcal{D}_m,0}$ such that

* The sequences $(\int_0^T |\delta_m \beta_m(t)|_{\star, \mathcal{D}_m})_{m \in \mathbb{N}}$ and $((|\nabla_{\mathcal{D}_m} \zeta_m|)_{L^2(0,T;L^p(\Omega)^d)})_{m \in \mathbb{N}}$ are bounded,

* As $m \to \infty$, $\Pi_{\mathcal{D}_m} \beta_m \to \bar{\beta}$ and $\Pi_{\mathcal{D}_m} \zeta_m \to \bar{\zeta}$ weakly in $L^2(\Omega \times (0,T))$.

Then $\Pi_{\mathcal{D}_m} \beta_m \Pi_{\mathcal{D}_m} \zeta_m \to \bar{\beta} \bar{\zeta}$ in the sense of measures on $\Omega \times (0, T)$, that is, for all $\varphi \in C(\bar{\Omega} \times [0, T])$,

$$ \lim_{m \to \infty} \int_0^T \int_{\Omega} \Pi_{\mathcal{D}_m} \beta_m(\mathbf{x}, t) \Pi_{\mathcal{D}_m} \zeta_m(\mathbf{x}, t) \varphi(\mathbf{x}, t) d\mathbf{x} dt = \int_0^T \int_{\Omega} \bar{\beta}(\mathbf{x}, t) \bar{\zeta}(\mathbf{x}, t) \varphi(\mathbf{x}, t) d\mathbf{x} dt. \quad (64) $$

**Proof.**

The idea of the proof is to reduce to the case where $\Pi_{\mathcal{D}_m} \zeta_m$ are tensorial functions, in order to separate the space and time variables and make use of the separate compactness of $\Pi_{\mathcal{D}_m} \zeta_m$ and $\Pi_{\mathcal{D}_m} \beta_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 $\Pi_{\mathcal{D}_m} \zeta_m$ to tensorial functions.

Let us take $\delta > 0$ and consider a covering $(A_k^\delta)_{k=1,...,K}$ of $\Omega$ in disjoint cubes of length $\delta$. Let us denote $R_\delta : L^2(\Omega) \to L^2(\Omega)$ the operator defined by:

$$ \forall k = 1, \dots, K, \quad \forall \mathbf{x} \in A_k^\delta : R_\delta f(\mathbf{x}) = \frac{1}{\text{meas}(A_k^\delta)} \int_{A_k^\delta} f(\mathbf{y}) d\mathbf{y}, $$

where $f$ has been extended by 0 outside $\Omega$. Let $\mathbf{x} \in A_k^\delta$. Using Jensen's inequality, fact that $\text{meas}(A_k^\delta) = \delta^d$ and the change of variable $\mathbf{y} \in A_k^\delta \mapsto \mathbf{\xi} = \mathbf{y} - \mathbf{x} \in (-\delta, \delta)^d$, we can write

$$ |R_\delta f(\mathbf{x}) - f(\mathbf{x})|^2 \le \delta^{-d} \int_{A_k^\delta} |f(\mathbf{y}) - f(\mathbf{x})|^2 d\mathbf{y} \le \delta^{-d} \int_{(-\delta, \delta)^d} |f(\mathbf{x} + \mathbf{\xi}) - f(\mathbf{x})|^2 d\mathbf{\xi}. $$

Integrating over $\mathbf{x} \in A_k^\delta$ and summing over $k = 1, \dots, K$ then gives

$$ \begin{aligned} ||R_\delta f - f||_{L^2(\Omega)}^2 & \le \delta^{-d} \int_{(-\delta, \delta)^d} \int_{\Omega} |f(\mathbf{x} + \mathbf{\xi}) - f(\mathbf{x})|^2 d\mathbf{x} d\mathbf{\xi} \\ & \le 2^d \sup_{\mathbf{\xi} \in (-\delta, \delta)^d} \int_{\Omega} |f(\mathbf{x} + \mathbf{\xi}) - f(\mathbf{x})|^2 d\mathbf{x}. \end{aligned} \quad (65) $$
On the other side, the compactness of $(\mathcal{D}_m)_{m \in \mathbb{N}}$ (Definition 2.7) and the fact that $p \ge 2$ give $\epsilon(\xi)$ such that $\epsilon(\xi) \to 0$ as $\xi \to 0$ and, for all $m \in \mathbb{N}$ and all $v \in X_{\mathcal{D}_m, 0}$,

$$||\Pi_{\mathcal{D}_m} v(\cdot + \xi) - \Pi_{\mathcal{D}_m} v||_{L^2(\Omega)}^2 \le \epsilon(\xi)||\nabla_{\mathcal{D}_m} v||_{L^p(\Omega)^d}^2.$$

Combining this with (65) and using the bound on $||\nabla_{\mathcal{D}_m} \zeta_m||_{L^2(0,T;L^p(\Omega)^d)}$ shows that

$$||R_{\delta}\Pi_{\mathcal{D}_m}\zeta_m - \Pi_{\mathcal{D}_m}\zeta_m||_{L^2(\Omega \times (0,T))} \le C \sup_{|\xi| \le \delta} \sqrt{\epsilon(\xi)} =: \omega(\delta) \quad (66)$$

where $C$ does not depend on $m$ and $\omega(\delta) \to 0$ as $\delta \to 0$. Note that a similar estimate holds with $\Pi_{\mathcal{D}_m} \zeta_m$ replaced with $\bar{\zeta}$ since $\bar{\zeta} \in L^2(\Omega \times (0, T))$ and, therefore, its translates tend to 0 in $L^2(\Omega \times (0, T))$.

If we respectively denote by $A_m(\Pi_{\mathcal{D}_m} \zeta_m)$ and $A(\bar{\zeta})$ the left-hand side and right-hand side (64), then since $(\Pi_{\mathcal{D}_m} \beta_m)_{m \in \mathbb{N}}$ is bounded in $L^2(\Omega \times (0, T))$ we have, by (66),

$$|A_m(\Pi_{\mathcal{D}_m} \zeta_m) - A(\bar{\zeta})| \le C\omega(\delta) + |A_m(R_{\delta}\Pi_{\mathcal{D}_m} \zeta_m) - A(R_{\delta}\bar{\zeta})|.$$

Let us assume that we can prove that, for a fixed $\delta$,

$$A_m(R_{\delta}\Pi_{\mathcal{D}_m}\zeta_m) \to A(R_{\delta}\bar{\zeta}) \text{ as } m \to \infty. \quad (67)$$

Then the previous inequality gives $\limsup_{m \to \infty} |A_m(\Pi_{\mathcal{D}_m} \zeta_m) - A(\bar{\zeta})| \le C\omega(\delta)$. Letting $\delta \to 0$ in this inequality gives $A_m(\Pi_{\mathcal{D}_m} \zeta_m) \to A(\bar{\zeta})$ as wanted.

Hence, we only need to prove 67. The definition of $R_{\delta}$ shows that $R_{\delta}f = \sum_{k=1}^K \mathbf{1}_{A_k^{\delta}} \text{meas}(A_k^{\delta})^{-1}[f]_{A_k^{\delta}}$ where $\mathbf{1}_{A_k^{\delta}}$ is the characteristic function of $A_k^{\delta}$ and $[f]_A = \int_A f(\mathbf{x})d\mathbf{x}$. Hence, (67) will follow if we can prove that, for any measurable set $A$ of non-zero measure,

$$\lim_{m \to \infty} \int_0^T \int_{\Omega} \Pi_{\mathcal{D}_m} \beta_m(\mathbf{x}, t) [\Pi_{\mathcal{D}_m} \zeta_m]_A(t) \varphi(t, \mathbf{x}) \mathbf{1}_A(\mathbf{x}) d\mathbf{x} dt \\ = \int_0^T \int_{\Omega} \bar{\beta}(\mathbf{x}, t) [\bar{\zeta}]_A(t) \varphi(t, \mathbf{x}) \mathbf{1}_A(\mathbf{x}) d\mathbf{x} dt \quad (68)$$

where, for $g \in L^2(\Omega \times (0, T))$, we set $[g]_A(t) = \int_A g(t, \mathbf{y})d\mathbf{y}$.

**Step 2:** further reductions.

We now reduce $\varphi$ to a tensorial function and $\mathbf{1}_A$ to a smooth function. It is well-known that there exists tensorial functions $\varphi_r = \sum_{l=1}^{L_r} \theta_{l,r}(t)\gamma_{l,t}(\mathbf{x})$, with $\theta_{l,r} \in C^{\infty}([0, T])$ and $\gamma_{l,r} \in C^{\infty}(\bar{\Omega})$, such that $\varphi_r \to \varphi$ uniformly on $\Omega \times (0, T)$ as $r \to \infty$. Moreover, there exists $\rho_r \in C_c^{\infty}(\Omega)$ such that $\rho_r \to \mathbf{1}_A$ in $L^2(\Omega)$ as $r \to \infty$.

Hence, as $r \to \infty$ the function $(t, \mathbf{x}) \mapsto \varphi_r(t, \mathbf{x})\rho_r(\mathbf{x})$ converges in $L^\infty(0, T; L^2(\Omega))$ to the function $(t, \mathbf{x}) \mapsto \varphi(t, \mathbf{x})\mathbf{1}_A(\mathbf{x})$. Since the sequence of functions $(t, \mathbf{x}) \mapsto \Pi_{\mathcal{D}_m} \beta(t, \mathbf{x})[\Pi_{\mathcal{D}_m} \zeta_m]_A(t)$ is bounded in $L^1(0, T; L^2(\Omega))$ (notice that $[\Pi_{\mathcal{D}_m} \zeta_m]_A)_{m \in \mathbb{N}}$ is bounded in $L^2(0, T)$ since $(\Pi_{\mathcal{D}_m} \zeta_m)_{m \in \mathbb{N}}$ is bounded in $L^2(\Omega \times (0, T))$), a reasoning similar as the one used in Step 1 shows that we only need to prove (68) with $\varphi(t, \mathbf{x})\mathbf{1}_A(\mathbf{x})$ replaced with $\varphi_r(t, \mathbf{x})\rho_r(\mathbf{x})$ for a fixed $r$.

We have $\varphi_r(t, \mathbf{x})\rho_r(\mathbf{x}) = \sum_{l=1}^{L_r} \theta_{l,r}(t)(\gamma_{l,t}\rho)(\mathbf{x})$ and $\gamma_{l,r}\rho_r \in C_c^{\infty}(\Omega)$. Hence, (68) with $\varphi(t, \mathbf{x})\mathbf{1}_A(\mathbf{x})$ replaced with $\varphi_r(t, \mathbf{x})\rho_r(\mathbf{x})$ will follow if we can establish that: for any $\theta \in C^{\infty}([0, T])$, any $\psi \in C_c^{\infty}(\Omega)$ and any measurable set $A \subset \Omega$,

$$\lim_{m \to \infty} \int_0^T \int_{\Omega} \theta(t) \Pi_{\mathcal{D}_m} \beta_m(\mathbf{x}, t) [\Pi_{\mathcal{D}_m} \zeta_m]_A(t) \psi(\mathbf{x}) d\mathbf{x} dt = \int_0^T \int_{\Omega} \theta(t) \bar{\beta}(\mathbf{x}, t) [\bar{\zeta}]_A(t) \psi(\mathbf{x}) d\mathbf{x} dt. \quad (69)$$
**Step 3:** proof of (69).

We now use the estimates on $\delta_m\beta_m$ to conclude. We first write

$$ \int_0^T \int_{\Omega} \theta(t) \Pi_{\mathcal{D}_m} \beta_m(x, t) [\Pi_{\mathcal{D}_m} \zeta_m]_A(t) \psi(x) \mathrm{d}x \mathrm{d}t = \int_0^T \theta(t) [\Pi_{\mathcal{D}_m} \zeta_m]_A(t) F_m(t). \quad (70) $$

with $F_m(t) = \int_{\Omega} \Pi_{\mathcal{D}_m} \beta_m(x, t) \psi(x) \mathrm{d}x$. It is clear from the weak convergence of $\Pi_{\mathcal{D}_m} \zeta_m$ that $[\Pi_{\mathcal{D}_m} \zeta_m]_A \to [\zeta]_A$ weakly in $L^2(0, T)$. Hence, if we can prove that $F_m \to F := \int_{\Omega} \bar{\beta}(x, \cdot) \psi(x) \mathrm{d}x$ 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(\Omega \times (0, T))$ of $\Pi_{\mathcal{D}_m} \beta_m$), we only have to prove that $(F_m)_{m \in \mathbb{N}}$ is relatively compact in $L^2(0, T)$.

To prove the strong convergence of $F_m$, we introduce the interpolant $P_{\mathcal{D}_m}$ defined by (21) and we define $G_m$ as $F_m$ with $\psi$ replaced with $\Pi_{\mathcal{D}_m} P_{\mathcal{D}_m} \psi$. We then have

$$ |F_m(t) - G_m(t)| \leq ||\Pi_{\mathcal{D}_m} \beta_m(\cdot, t)||_{L^2(\Omega)} S_{\mathcal{D}_m}(\psi). $$

The consistency of $(\mathcal{D}_m)$ then shows that

$$ F_m - G_m \to 0 \text{ strongly in } L^2(\Omega) \text{ as } m \to \infty. \quad (71) $$

We now study the strong convergence of $G_m$. This function is, as $\Pi_{\mathcal{D}_m} \beta_m$, piecewise constant on $(0, T)$ and, by definition of $|\cdot|_{\star, \mathcal{D}_m}$, its discrete derivative satisfies

$$ |\delta_m G_m(t)| \leq |\delta_m \beta_m(t)|_{\star, \mathcal{D}_m} ||P_{\mathcal{D}_m} \psi||_{X_{\mathcal{D}_m, 0}}. $$

Since $||P_{\mathcal{D}_m} \psi||_{X_{\mathcal{D}_m, 0}} \leq S_{\mathcal{D}_m}(\psi) + ||\nabla \psi||_{L^p(\Omega)^d}$ is bounded uniformly with respect to $m$, the assumption on $\delta_m \beta_m$ proves that $(||\delta_m G_m||_{L^1(0,T)})_{m \in \mathbb{N}}$ is bounded. Hence, $G_m$ is bounded in $BV(0,T) \cap L^1(0,T)$ and therefore relatively compact in $L^2(0,T)$ (see [3, Theorem 10.1.4]). Combined with (71), this proves that $(F_m)_{m \in \mathbb{N}}$ is relatively compact in $L^2(0,T)$ and concludes the proof. ■

**Remark 5.5** If we assume that $(\Pi_{\mathcal{D}_m} \beta_m)_{m \in \mathbb{N}}$ is bounded in $L^\infty(0, T; L^2(\Omega))$ and that, for some $q > 1$, $(\int_0^T |\delta_m \beta_m(t)|_{\star, \mathcal{D}_m}^q)_{m \in \mathbb{N}}$ is bounded, then Step 2 becomes a trivial consequence of Theorem 3.1. Indeed, this theorem then shows that $(\Pi_{\mathcal{D}_m} \beta_m)_{m \in \mathbb{N}}$ is relatively compact uniformly-in-time and weakly in $L^2(\Omega)$, which translates into the relative compactness of $(F_m)_{m \in \mathbb{N}}$ in $L^\infty(0, T)$.

**Lemma 5.6** Let $N \in \mathbb{N}^*$ and $V$ be a non-empty measurable subset of $\mathbb{R}^N$. Let $\beta, \zeta \in C^0(\mathbb{R})$ be two nondecreasing functions, such that $\beta(0) = \zeta(0) = 0$ and $\beta + \zeta$ is strictly increasing. We assume that there exists a sequence $(w_m)_{m \in \mathbb{N}}$ of measurable functions on $V$ and two functions $\bar{\beta}, \bar{\zeta} \in L^2(V)$ such that:

* $\beta(w_m) \to \bar{\beta}$ and $\zeta(w_m) \to \bar{\zeta}$ weakly in $L^2(V)$,

* there exists $\varphi \in L^\infty(V)$ such that $\varphi > 0$ a.e. on $V$ and

$$ \lim_{m \to \infty} \int_V \varphi(z) \beta(w_m(z)) \zeta(w_m(z)) \mathrm{d}z = \int_V \varphi(z) \bar{\beta}(z) \bar{\zeta}(z) \mathrm{d}z. \quad (72) $$

Then

$$ \bar{\beta} = \beta(w) \text{ and } \bar{\zeta} = \zeta(w) \text{ a.e. in } V, \text{ where } w = \left( \frac{\beta + \zeta}{2} \right)^{-1} \left( \frac{\bar{\beta} + \bar{\zeta}}{2} \right). \quad (73) $$
**Proof.** We first notice that $\beta(w), \zeta(w) \in L^2(V)$, given that $\bar{\beta}+\bar{\zeta} \in L^2(V)$ and $\beta \circ (\frac{\beta+\zeta}{2})^{-1}$ and $\zeta \circ (\frac{\beta+\zeta}{2})^{-1}$ are sub-linear (the sum of both is 2Id and each one has the same sign as $s$). Since $\beta$ and $\zeta$ are non-decreasing, we can therefore write

$$ \int_V \varphi(z) [\beta(w_m(z)) - \beta(w(z))] [\zeta(w_m(z)) - \zeta(w(z))] dz \geq 0. $$

Letting $m \to \infty$ in the above inequality, and using the convergences of $\beta(w_m)$, $\zeta(w_m)$ and (72), we obtain

$$ \int_V \varphi(z) [\bar{\beta}(z) - \beta(w(z))] [\bar{\zeta}(z) - \zeta(w(z))] dz \geq 0. \quad (74) $$

We then remark that

$$ \frac{\bar{\beta} + \bar{\zeta}}{2} = \frac{\beta(w) + \zeta(w)}{2}, \quad (75) $$

which gives $\beta(w) = \frac{\bar{\beta}+\bar{\zeta}}{2} + (\frac{\beta-\zeta}{2}) (w)$ and $\zeta(w) = \frac{\bar{\beta}+\bar{\zeta}}{2} - (\frac{\beta-\zeta}{2}) (w)$. Hence, (74) leads to

$$ -\int_V \varphi(z) \left( \frac{\bar{\beta} - \bar{\zeta}}{2}(z) - \left( \frac{\beta - \zeta}{2} \right) (w(z)) \right)^2 dz \geq 0. $$

Since $\varphi$ is almost everywhere strictly positive on $V$, we deduce that $\frac{\bar{\beta}-\bar{\zeta}}{2} = \frac{\beta(w)-\zeta(w)}{2}$ a.e. in $V$, and (73) follows from (75). ■

# 6 Appendix: uniform-in-time compactness results for time-dependent 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 continuous, 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 non-continuous functions of the classical Ascoli-Arzela theorem.

**Definition 6.1** If $(K, d_K)$ and $(E, d_E)$ are metric spaces, we denote by $\mathcal{F}(K, E)$ the space of functions $K \to E$, endowed with the uniform metric $d_{\mathcal{F}}(v, w) = \sup_{s \in 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 $(\mathcal{F}(K, E), d_{\mathcal{F}})$ as in Definition 6.1.

Let $(v_m)_{m \in \mathbb{N}}$ be a sequence in $\mathcal{F}(K, E)$ such that there exists a function $\omega: K \times K \to [0, \infty]$ and a sequence $(\delta_m)_{m \in \mathbb{N}} \subset [0, \infty)$ satisfying

$$ \lim_{d_K(s,s') \to 0} \omega(s,s') = 0, \quad \lim_{m \to \infty} \delta_m = 0, \\ \forall (s,s') \in K^2, \forall m \in \mathbb{N}, d_E(v_m(s), v_m(s')) \le \omega(s,s') + \delta_m. \quad (76) $$

We also assume that, for all $s \in K$, $\{v_m(s) : m \in \mathbb{N}\}$ is relatively compact in $(E, d_E)$.

Then $(v_m)_{m \in \mathbb{N}}$ is relatively compact in $(\mathcal{F}(K, E), d_{\mathcal{F}})$ and any adherence value of $(v_m)_{m \in \mathbb{N}}$ in this space is continuous $K \to E$.
**Proof.** Let us first notice that the last conclusion of the theorem, i.e. that any adherence value $v$ of $(v_m)_{m \in \mathbb{N}}$ in $\mathcal{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 $\omega$.

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 \in \mathbb{N}\}$ in $K$ (the existence of this set is ensured since $K$ is compact metric). Since each set $\{v_m(s_l) : m \in \mathbb{N}\}$ is relatively compact in $E$, by diagonal extraction we can select a subsequence of $(v_m)_{m \in \mathbb{N}}$, denoted the same way, such that for any $l \in \mathbb{N}$, $(v_m(s_l))_{m \in \mathbb{N}}$ converges in $E$. We then proceed in showing that $(v_m)_{m \in \mathbb{N}}$ is a Cauchy sequence in $(\mathcal{F}(K, E), d_{\mathcal{F}})$. Since this space is complete, this will show that this sequence converges in this space and will therefore complete the proof.

Let $\varepsilon > 0$ and, using (76), take $\delta > 0$ and $M \in \mathbb{N}$ such that $\omega(s, s') \le \varepsilon$ whenever $d_K(s, s') \le \delta$ and $\delta_m \le \varepsilon$ whenever $m \ge M$. Select a finite set $\{s_{l_1}, \dots, s_{l_N}\}$ such that any $s \in K$ is within distance $\delta$ of a $s_{l_i}$. Then, for any $m, m' \ge M$, by (76),

$$ 
\begin{aligned}
d_E(v_m(s), v_{m'}(s)) &\le d_E(v_m(s), v_m(s_{l_i})) + d_E(v_m(s_{l_i}), v_{m'}(s_{l_i})) \\
&\quad + d_E(v_{m'}(s_{l_i}), v_{m'}(s)) \\
&\le \omega(s, s_{l_i}) + \delta_m + d_E(v_m(s_{l_i}), v_{m'}(s_{l_i})) + \omega(s, s_{l_i}) + \delta_{m'} \\
&\le 4\varepsilon + d_E(v_m(s_{l_i}), v_{m'}(s_{l_i})).
\end{aligned}
$$

Since $\{(v_m(s_{l_i}))_{m \in \mathbb{N}} : i = 1, \dots, N\}$ form a finite number of converging sequences in $E$, we can find $M' \ge M$ such that, whenever $m, m' \ge M'$ and $i = 1, \dots, N$, $d_E(v_m(s_{l_i}), v_{m'}(s_{l_i})) \le \varepsilon$. This shows that, for all $m, m' \ge M'$ and all $s \in K$, $d_E(v_m(s), v_{m'}(s)) \le 5\varepsilon$ and concludes the proof that $(v_m)_{m \in \mathbb{N}}$ is a Cauchy sequence in $(\mathcal{F}(K, E), d_{\mathcal{F}})$. ■

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.

**Lemma 6.3** Let $(K, d_K)$ be a compact metric space, $(E, d_E)$ be a metric space and $(\mathcal{F}(K, E), d_{\mathcal{F}})$ as in Definition 6.1. Let $(v_m)_{m \in \mathbb{N}}$ be a sequence in $\mathcal{F}(K, E)$ and $v : K \to E$ be continuous. Then $v_m \to v$ for $d_{\mathcal{F}}$ if and only if, for any $s \in K$ and any sequence $(s_m)_{m \in \mathbb{N}} \subset K$ converging to $s$ for $d_K$, we have $v_m(s_m) \to v(s)$ for $d_E$.

**Proof.** If $v_m \to v$ for $d_{\mathcal{F}}$ then for any sequence $(s_m)_{m \in \mathbb{N}}$ converging to $s$,

$$ d_E(v_m(s_m), v(s)) \le d_E(v_m(s_m), v(s_m)) + d_E(v(s_m), v(s)) \le d_{\mathcal{F}}(v_m, v) + d(v(s_m), v(s)) $$

and the right-hand side tends to 0 by definition of $v_m \to v$ for $d_{\mathcal{F}}$ and by continuity of $v$. Let us now prove the converse by contraposition. If $(v_m)_{m \in \mathbb{N}}$ does not converge to $v$ for $d_{\mathcal{F}}$ then there exists $\varepsilon > 0$ and a subsequence $(v_{m_k})_{k \in \mathbb{N}}$, such that, for any $k \in \mathbb{N}$, $\sup_{s \in K} d_E(v_{m_k}(s), v(s)) \ge \varepsilon$. We can then find a sequence $(r_k)_{k \in \mathbb{N}} \subset K$ such that, for any $k \in \mathbb{N}$,

$$ d_E(v_{m_k}(r_k), v(r_k)) \ge \varepsilon/2. \qquad (77) $$

$K$ being compact, up to another subsequence denoted the same way, we can assume that $r_k \to s$ in $K$ as $k \to \infty$. It is then trivial to construct a sequence $(s_m)_{m \in \mathbb{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) \to v(s)$ in $E$ and, by continuity of $v$, $v(s_m) \to v(s)$ in $E$. This shows that $d_E(v_m(s_m), v(s_m)) \to 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(\Omega)$.
**Definition 6.4 (Uniform-in-time $L^2(\Omega)$-weak convergence)** A sequence of functions $u_m : [0, T] \to L^2(\Omega)$ converges weakly in $L^2(\Omega)$ uniformly on $[0, T]$ to a function $u : [0, T] \to L^2(\Omega)$ if, for all $\varphi \in L^2(\Omega)$, as $m \to \infty$ the sequence of functions $t \in [0, T] \to \langle u_m(t), \varphi \rangle_{L^2(\Omega)}$ converges uniformly on $[0, T]$ to $t \in [0, T] \to \langle u(t), \varphi \rangle_{L^2(\Omega)}$, where $\langle \cdot, \cdot \rangle_{L^2(\Omega)}$ is the inner product in $L^2(\Omega)$.

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(\Omega)$ and $(\varphi_l)_{l \in \mathbb{N}}$ be a dense sequence in $L^2(\Omega)$. Then, on E, the weak topology of $L^2(\Omega)$ is the topology given by the metric

$$d_E(v, w) = \sum_{l \in \mathbb{N}} \frac{\min(1, |\langle v - w, \varphi_l \rangle_{L^2(\Omega)}|)}{2^l}. \qquad (78)$$

Moreover, a sequence of functions $u_m : [0, T] \to E$ converges uniformly to $u : [0, T] \to E$ for the weak topology of $L^2(\Omega)$ if and only if, as $m \to \infty$, $d_E(u_m, u) : [0, T] \to [0, \infty)$ converges uniformly to 0.

**Proof.** The sets $E_{\varphi, \varepsilon} = \{v \in E : |\langle v, \varphi \rangle_{L^2(\Omega)}| < \varepsilon\}$, for $\varphi \in L^2(\Omega)$ and $\varepsilon > 0$, define a basis of neighborhood of 0 for the weak $L^2(\Omega)$ 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 $\varphi \in L^2(\Omega)$, $l \in \mathbb{N}$ and $v \in E$ we have

$$|\langle v, \varphi \rangle_{L^2(\Omega)}| \le R ||\varphi - \varphi_l||_{L^2(\Omega)} + |\langle v, \varphi_l \rangle_{L^2(\Omega)}|.$$

By density of $(\varphi_l)_{l \in \mathbb{N}}$ we can select $l \in \mathbb{N}$ such that $||\varphi - \varphi_l||_{L^2(\Omega)} < \varepsilon/(2R)$ and we then see that $E_{\varphi_l, \varepsilon/2} \subset E_{\varphi, \varepsilon}$. Hence, a basis of neighborhood of 0 in E for the weak $L^2(\Omega)$ is also given by $(E_{\varphi_l, \varepsilon})_{l \in \mathbb{N}, \varepsilon > 0}$.

From the definition of $d_E$ we see that, for any $l \in \mathbb{N}$, $\min(1, |\langle v, \varphi_l \rangle_{L^2(\Omega)}|) \le 2^l d_E(0, v)$. If $d_E(0, v) < 2^{-l}$ this shows that $|\langle v, \varphi_l \rangle_{L^2(\Omega)}| \le 2^l d_E(0, v)$ and therefore that

$$B_{d_E}(0, \min(2^{-l+1}, \varepsilon 2^{-l})) \subset E_{\varphi_l, \varepsilon}.$$

Hence, any neighborhood of 0 in E for the $L^2(\Omega)$ weak topology is a neighborhood of 0 for $d_E$. Conversely, for any $\varepsilon > 0$, selecting $N \in \mathbb{N}$ such that $\sum_{l \ge N+1} 2^{-l} < \varepsilon/2$ gives, from the definition (78) of $d_E$,

$$\bigcap_{l=1}^{N} E_{\varphi_l, \varepsilon/4} \subset B_{d_E}(0, \varepsilon).$$

Hence, any ball for $d_E$ centered at 0 is a neighborhood of 0 for the $L^2(\Omega)$ weak topology. Since $d_E$ and the $L^2(\Omega)$ 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 \in [0, T]$. ■

The following lemma has been established in [28, Proposition 9.3] but its proof is recalled for the reader's convenience.

**Lemma 6.6**

Let $(t^{(n)})_{n \in \mathbb{Z}}$ be a strictly increasing sequence of real values such that $\delta t^{(n+1/2)} := t^{(n+1)} - t^{(n)}$ is uniformly bounded by $\delta t > 0$, $\lim_{n \to -\infty} t^{(n)} = -\infty$ and $\lim_{n \to \infty} t^{(n)} = \infty$. For all $t \in \mathbb{R}$, we denote by
$n(t)$ the element $n \in \mathbb{Z}$ such that $t \in [t^{(n)}, t^{(n+1)})$. Let $(a^{(n)})_{n \in \mathbb{Z}}$ be a family of non negative real values with a finite number of non zero values. Then

$$ \int_{\mathbb{R}} \sum_{n=n(t)+1}^{n(t+\tau)} (\delta t^{(n+\frac{1}{2})} a^{(n+1)}) dt = \tau \sum_{n \in \mathbb{Z}} (\delta t^{(n+\frac{1}{2})} a^{(n+1)}), \quad \forall \tau \in (0, +\infty), \quad (79) $$

and

$$ \int_{\mathbb{R}} \left( \sum_{n=n(t)+1}^{n(t+\tau)} \delta t^{(n+\frac{1}{2})} \right) a^{n(t+\zeta)+1} dt \\ \le (\tau + \delta t) \sum_{n \in \mathbb{Z}} (\delta t^{(n+\frac{1}{2})} a^{(n+1)}), \quad \forall \tau \in (0, +\infty), \quad \forall \zeta \in \mathbb{R}. \quad (80) $$

**Proof.**

Let us define the function $\chi(t, n, \tau)$ by $\chi(t, n, \tau) = 1$ if $t < t^{(n)}$ and $t + \tau \ge t^{(n)}$, else $\chi(t, n, \tau) = 0$. We have

$$ \int_{\mathbb{R}} \sum_{n=n(t)+1}^{n(t+\tau)} (\delta t^{(n+\frac{1}{2})} a^{(n+1)}) dt = \int_{\mathbb{R}} \sum_{n \in \mathbb{Z}} (\delta t^{(n+\frac{1}{2})} a^{(n+1)} \chi(t, n, \tau)) dt \\ = \sum_{n \in \mathbb{Z}} (\delta t^{(n+\frac{1}{2})} a^{(n+1)} \int_{\mathbb{R}} \chi(t, n, \tau) dt). $$

Since $\int_{\mathbb{R}} \chi(t, n, \tau) dt = \int_{t^{(n)}-\tau}^{t^{(n)}} dt = \tau$, thus (79) is proven.

We now turn to the proof of (80). We define the function $\tilde{\chi}(n, t)$ by $\tilde{\chi}(n, t) = 1$ if $n(t) = n$, else $\tilde{\chi}(n, t) = 0$. We have

$$ \int_{\mathbb{R}} \left( \sum_{n=n(t)+1}^{n(t+\tau)} \delta t^{(n+\frac{1}{2})} \right) a^{(n(t+\zeta)+1)} dt \\ = \int_{\mathbb{R}} \left( \sum_{n=n(t)+1}^{n(t+\tau)} \delta t^{(n+\frac{1}{2})} \right) \sum_{m \in \mathbb{Z}} a^{(m+1)} \tilde{\chi}(m, t+\zeta) dt, $$

which yields

$$ \int_{\mathbb{R}} \left( \sum_{n=n(t)+1}^{n(t+\tau)} \delta t^{(n+\frac{1}{2})} \right) a^{(n(t+\zeta)+1)} dt \\ = \sum_{m \in \mathbb{Z}} a^{(m+1)} \int_{t^m - \zeta}^{t^{m+1} - \zeta} \left( \sum_{n=n(t)+1}^{n(t+\tau)} \delta t^{(n+\frac{1}{2})} \right) dt. \quad (81) $$

Since we have

$$ \sum_{n=n(t)+1}^{n(t+\tau)} \delta t^{(n+\frac{1}{2})} = \sum_{n \in \mathbb{Z}, t < t^{(n)} \le t+\tau} (t^{(n+1)} - t^{(n)}) \le \tau + \delta t, $$

we can write from (81)

$$ \int_{\mathbb{R}} \left( \sum_{n=n(t)+1}^{n(t+\tau)} \delta t^{(n+\frac{1}{2})} \right) a^{(n(t+\zeta)+1)} dt \\ \le (\tau + \delta t) \sum_{m \in \mathbb{Z}} a^{(m+1)} \int_{t^{(m)} - \zeta}^{t^{(m+1)} - \zeta} dt \\ = (\tau + \delta t) \sum_{m \in \mathbb{Z}} a^{(m+1)} \delta t^{(m+\frac{1}{2})}, $$
which is exactly (80).

References

[1] I. Aavatsmark, T. Barkve, O. Boe, and T. Mannseth. Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media. *J. Comput. Phys.*, 127(1):2–14, 1996.

[2] B. Andreianov, F. Boyer, and F. Hubert. Discrete duality finite volume schemes for Leray-Lions-type elliptic problems on general 2D meshes. *Numer. Methods Partial Differential Equations*, 23(1):145–195, 2007.

[3] H. Attouch, G. Buttazzo, and G. Michaille. *Variational analysis in Sobolev and BV spaces*, volume 6 of *MPS/SIAM Series on Optimization*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Programming Society (MPS), Philadelphia, PA, 2006.

[4] M. Bertsch, P. De Mottoni, and L. Peletier. The Stefan problem with heating: appearance and disappearance of a mushy region. *Trans. Amer. Math. Soc*, 293:677–691, 1986.

[5] F. Brezzi and M. Fortin. *Mixed and hybrid finite element methods*, volume 15 of Springer Series in Computational Mathematics. Springer-Verlag, New York, 1991.

[6] F. Brezzi, K. Lipnikov, and V. Simoncini. A family of mimetic finite difference methods on polygonal and polyhedral meshes. *Math. Models Methods Appl. Sci.*, 15(10):1533–1551, 2005.

[7] P. Ciarlet. The finite element method. In P. G. Ciarlet and J.-L. Lions, editors, *Part I*, Handbook of Numerical Analysis, III. North-Holland, Amsterdam, 1991.

[8] Y. Coudière and F. Hubert. A 3d discrete duality finite volume method for nonlinear elliptic equations. *SIAM Journal on Scientific Computing*, 33(4):1739–1764, 2011.

[9] M. Crouzeix and P.-A. Raviart. Conforming and nonconforming finite element methods for solving the stationary Stokes equations. I. *Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge*, 7(R-3):33–75, 1973.

[10] K. Deimling. *Nonlinear functional analysis*. Springer-Verlag, Berlin, 1985.

[11] J. Diaz and F. de Thelin. On a nonlinear parabolic problem arising in some models related to turbulent flows. *SIAM J. Math. Anal.*, 25(4):1085–1111, 1994.

[12] J. Droniou. Finite volume schemes for fully non-linear elliptic equations in divergence form. *ESAIM: Mathematical Modelling and Numerical Analysis*, 40(6):1069, 2006.

[13] J. Droniou and R. Eymard. A mixed finite volume scheme for anisotropic diffusion problems on any grid. *Numer. Math.*, 105(1):35–71, 2006.

[14] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. Gradient schemes for elliptic and parabolic problems. 2014. To appear.

[15] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin. A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods. *Math. Models Methods Appl. Sci.*, 20(2):265–295, 2010.
[16] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin. Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations. *Math. Models Methods Appl. Sci.* (M3AS), 23(13):2395-2432, 2013.

[17] J. Droniou, R. Eymard, and C. Guichard. Uniform-in-time convergence of numerical schemes for Richards' and Stefan's models. In *Finite Volumes for Complex Applications VII*. Springer, 2014.

[18] M. G. Edwards and C. F. Rogers. Finite volume discretization with imposed flux continuity for the general tensor pressure equation. *Comput. Geosci.*, 2(4):259-290, 1998.

[19] I. Ekeland and R. Témam. *Convex analysis and variational problems*, volume 28 of *Classics in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, english edition, 1999. Translated from the French.

[20] A. Ern and J.-L. Guermond. *Theory and practice of finite elements*, volume 159. Springer, 2004.

[21] R. Eymard, P. Féron, T. Gallouet, R. Herbin, and C. Guichard. Gradient schemes for the Stefan problem.

[22] R. Eymard, T. Gallouët, and R. Herbin. Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces. *IMA J. Numer. Anal.*, 30(4):1009-1043, 2010.

[23] R. Eymard, T. Gallouët, D. Hilhorst, and Y. Naït Slimane. Finite volumes and nonlinear diffusion equations. *RAIRO Modél. Math. Anal. Numér.*, 32(6):747-761, 1998.

[24] R. Eymard, C. Guichard, and R. Herbin. Small-stencil 3d schemes for diffusive flows in porous media. *M2AN*, 46:265-290, 2012.

[25] R. Eymard, C. Guichard, R. Herbin, and R. Masson. Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation. *ZAMM*, page accepted for publication, 2013.

[26] R. Eymard, M. Gutnic, and D. Hilhorst. The finite volume method for Richards equation. *Comput. Geosci.*, 3(3-4):259-294, 1999.

[27] R. Eymard and R. Herbin. Gradient scheme approximations for diffusion problems. *Finite Volumes for Complex Applications VI Problems & Perspectives*, pages 439-447, 2011.

[28] R. Eymard, R. Herbin, and A. Michel. Mathematical study of a petroleum-engineering scheme. *M2AN Math. Model. Numer. Anal.*, 37(6):937-972, 2003.

[29] R. Glowinski and J. Rappaz. Approximation of a nonlinear elliptic problem arising in a non-newtonian fluid flow model in glaciology. *M2AN Math. Model. Numer. Anal.*, 37(1):175-186, 2003.

[30] F. Hermeline. Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes. *Computer methods in applied mechanics and engineering*, 192(16):1939-1959, 2003.

[31] A. V. Kazhikhov. Recent developments in the global theory of two-dimensional compressible Navier-Stokes equations. 25:ii+52, 1998.
[32] E. Maitre. Numerical analysis of nonlinear elliptic-parabolic equations. *M2AN Math. Model. Numer. Anal.*, 36(1):143-153, 2002.

[33] G. Minty. On a monotonicity method for the solution of non-linear equations in Banach spaces. *Proceedings of the National Academy of Sciences of the United States of America*, 50(6):1038, 1963.

[34] R. H. Nochetto and C. Verdi. Approximation of degenerate parabolic problems using numerical integration. *SIAM J. Numer. Anal.*, 25(4):784-814, 1988.

[35] I. S. Pop. Numerical schemes for degenerate parabolic problems. In *Progress in industrial mathematics at ECMI 2004*, volume 8 of *Math. Ind.*, pages 513-517. Springer, Berlin, 2006.