François Bouchut
email: francois.bouchut@u-pem.fr

David Doyen
email: david.doyen@u-pem.fr

Robert Eymard
email: robert.eymard@u-pem.fr

Convection and total variation flow

Keywords: Bingham fluid, hyperbolic scalar conservation law, total variation flow, 1-Laplacian, entropy formulation, finite volumes, finite elements

published or not. The documents may come

Introduction

A Bingham fluid, also called rigid viscoplastic fluid, is a material that behaves as a rigid solid below a certain yield stress and as a viscous fluid above this yield; a familiar example of such a material is the tooth paste. For a d-dimensional Bingham fluid, the relation between the stress tensor σ, seen as a d × d matrix, the pressure p and the velocity u is

σ = -pI d + g D(u) |D(u)| F + 2νD(u), (1)
where g and ν are positive constants, I d is the d × d identity matrix, D(u) is the d × d matrix such that D(u) ij := 1 2 (∂ i u j + ∂ j u i ), and | • | F denotes the Frobenius matrix norm. The term g D(u) |D(u)| F enforces the plastic behaviour, g being the plasticity yield stress, while the term νD(u) enforces the viscous behaviour, ν being the viscosity parameter. The mathematical analysis of Bingham fluid flows dates back to the work of [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF], where the problems are formulated as variational inequalities in Sobolev spaces. The numerical approximation of a Bingham fluid flow is usually treated with finite element techniques; we refer to [START_REF] Dean | On the numerical simulation of bingham visco-plastic flow: Old and new results[END_REF] for a recent review. When the viscosity becomes negligible (ν = 0), the analytical and numerical framework described above is no longer suitable -let us mention however an existence result in 2D obtained by [START_REF] Lions | Remarks on some nonlinear evolution problems arising in Bingham flows[END_REF]. Although the study of inviscid Bingham fluids has been initiated in [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF] with the case of an unsteady flow without convection term, the presence of a nonlinear convection term is naturally issued from the inertial term in the momentum conservation equation. Unfortunately, the study of this problem seems to be out of reach in the actual state of the art. Therefore, we consider here a simplified model of unsteady Bingham flow with convection. This simplified model is scalar and consists in seeking u : R d × (0, T ) → R and λ : R d × (0, T ) → R d such that ∂ t u + div F (x, t, u) -div λ = 0 and λ ∈ Sgn(∇u), on

Q T := R d × (0, T ), (2)
u(x, 0) = u ini (x), on R d , (3)
where d ∈ N , T > 0 is given, F : R d × (0, T ) × R → R d is divergence-free with respect to the space variables and u ini is a given function from R d to R. We denote by Sgn the vector sign function, which is the set-valued map from R d to P(R d ) defined by

λ ∈ Sgn(µ) ⇔ |λ| ≤ 1 if µ = 0 λ = µ |µ| if µ = 0,
where | • | denotes the euclidean norm in R d . In equation [START_REF] Bartels | Total variation minimization with finite elements: Convergence and interative solution[END_REF], the term div F describes the convection in the fluid regime, while the term div (Sgn(∇u)) enforces the plastic behaviour (the plasticity yield is taken equal to 1 for simplicity). Problem (2)-( 3) is considered in this work under the following hypotheses, denoted by Hypotheses (HC) in this paper.

(HC1) The initial datum u ini is assumed to belong to L 1 (R d ) ∩ L ∞ (R d ). The essential infimum and supremum of u ini are denoted by a 0 and b 0 , respectively.

(HC2) The flux function

F ∈ C 1 (Q T × R, R d
) is assumed to be divergence-free with respect to the space variables, that is

d i=1 ∂F i ∂x i (x, t, u) = 0, ∀x = (x i ) i=1,...,d ∈ R d , ∀t ∈ [0, T ], ∀u ∈ R. (4)
Furthermore, ∂F ∂u is assumed to be locally Lipschitz continuous and such that, for all compact set K ⊂ R, ∂F ∂u ≤ C K a.e. on Q T × K, where C K is a constant depending on K. Problem ( 2)-( 3) can be viewed as a nonlinear hyperbolic scalar conservation law regularised by the total variation flow operator (or the 1-Laplacian). For nonlinear hyperbolic conservation laws, it is well known that the standard weak formulation fails to ensure the uniqueness of the solution and must be replaced by an entropy formulation; see, e.g., [START_REF] Godlewski | Hyperbolic systems of conservation laws, volume 3/4 of Mathématiques & Applications[END_REF] or [START_REF] Serre | Systems of conservation laws[END_REF]. With some types of regularisation, as for instance the viscous regularisation, the uniqueness is recovered. That is not the case with the total variation flow regularisation, which has no spatial smoothing effect and does not prevent the formation of shocks. Nonlinear hyperbolic problems are usually approximated with finite volumes [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. Unfortunately, finite volumes are not suitable for the approximation of the total variation flow: indeed, if a sequence (u k ) k∈N of piece-wise constant functions converges to u in L 1 , the total variation of u k does not converge in general to the total variation of u (see [START_REF] Bělík | Approximation by piecewise constant functions in a BV metric[END_REF] for an example). The total variation flow must be approximated in W 1,1 -conforming discrete spaces, such as P 1 finite element spaces [START_REF] Bartels | Total variation minimization with finite elements: Convergence and interative solution[END_REF][START_REF] Feng | Analysis of total variation flow and its finite element approximations[END_REF][START_REF] Feng | Rate of convergence of regularization procedures and finite element approximations for the total variation flow[END_REF]. Numerical schemes combining finite volumes and finite element schemes have already been considered for scalar conservation laws with a diffusion term [START_REF] Feistauer | Error estimates for a combined finite volume-finite element method for nonlinear convection-diffusion problems[END_REF] and for degenerate parabolic equations [START_REF] Eymard | A combined finite volume-nonconforming/mixed-hybrid finite element scheme for degenerate parabolic problems[END_REF]. In the present article, we first give an entropy formulation for Problem (2)- [START_REF] Bělík | Approximation by piecewise constant functions in a BV metric[END_REF] and prove the uniqueness of the entropy solution using the doubling variable technique. Note that our entropy formulation of the total variation term is similar to the one developed in [START_REF] Bellettini | The total variation flow in R N[END_REF] to study the total variation flow with L 1 loc initial data (without hyperbolic term). The existence of the entropy solution follows from the convergence of a numerical approximation, based on a splitting scheme. The hyperbolic flow is treated with finite volumes and the total variation flow is treated with P 1 finite elements. The finite volume mesh is built as a dual mesh of the finite element mesh, which makes simple the interpolation step between the two meshes. For the hyperbolic step (or finite volume step), we choose an explicit time discretisation for sake of simplicity. For the total variation flow step (or finite element step), we are led to define an implicit scheme accounting for the nonregularity of the total variation flow operator. To guarantee the maximum principle, which is essential for the stability of the scheme, we use a nonobtuse finite element mesh. A small parabolic regularisation term had to be included in this step; the magnitude of this term is controlled by some function θ(h), where h is the size of the mesh, for proving, using entropies, the lower semi-continuity of the bounded variation norm of the discrete solution (see the use of this term in (81)). The convergence proof of the numerical approximation relies on the Kolmogorov-Riesz compactness theorem, which provides us with the strong convergence in L 1 loc (R d × (0, T )) of the discrete solutions. It requires uniform estimates on the space and time translates of the discrete solutions. To establish these estimates, the total variation term is crucial. For scalar conservation laws without total variation flow regularisation, these estimates are not true and the convergence study of the numerical approximations must be carried out with other tools [START_REF] Cockburn | Convergence of the finite volume method for multidimensional conservation laws[END_REF][START_REF] Chainais-Hillairet | Finite volume schemes for a nonlinear hyperbolic equation. Convergence towards the entropy solution and error estimate[END_REF][START_REF] Eymard | Finite volume methods[END_REF]. The article is organised as follows. In section 2, the concept of entropy solution for Problem (2)-( 3) is defined and its uniqueness is proved. Section 3 describes the numerical approximation and its first properties (wellposedness, maximum principle). A priori estimates on the discrete solutions are provided in Section 4 and a discrete entropy formulation is established in Section 5. The convergence of the numerical approximation (and thus the existence of an entropy solution) is finally proved in Section 6 using the results of the two previous sections. In the last section, some numerical simulations in 1D and 2D are presented.

2 Entropy formulation for nonlinear hyperbolic equation with total variation flow

Functions with bounded variation

Let us first recall basic properties concerning the functions with bounded variation. For a comprehensive presentation, we refer to [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF], [START_REF] Evans | Measure theory and fine properties of functions[END_REF] or [START_REF] Ziemer | Weakly differentiable functions[END_REF].

• Let Ω be an open subset of R d . The total variation of a function u ∈ L 1 (Ω) is defined by

T V Ω (u) := sup Ω u divφ dx; φ ∈ C 1 c (Ω; R d ) with φ L ∞ (Ω) ≤ 1 .
In particular, the total variation of a function u ∈ W 1,1 (Ω) is equal to

T V Ω (u) = Ω |∇u| dx.
The space of functions over Ω with bounded variation, denoted by BV (Ω), is the set of functions u ∈ L 1 (Ω) such that T V Ω (u) < +∞. Equipped with the norm

• BV (Ω) := • L 1 (Ω) + T V Ω (•), the space BV (Ω) is a (nonreflexive) Banach space.
• The distributional derivative of u ∈ BV (Ω), denoted by Du, is the vector Radon measure such that

Ω u divφ dx = - Ω φ Du, ∀φ ∈ C 1 c (Ω, R d ).
• The norm of the vector measure Du is denoted by |Du|. It is a positive Radon measure and there is a measurable function h : R → R d , with |h(x)| = 1 for all x ∈ R, such that Du = h|Du|.

• The norm of Du is linked to the total variation by the identity |Du|(Ω) = T V Ω (u).

• The total variation is lower semi-continuous relatively to the convergence in L 1 loc . In other words, if (u k ) k∈N is a sequence in BV (Ω) converging to u in L 1 loc (Ω), then

lim k→∞ inf |Du k |(Ω) ≥ |Du|(Ω). • If u ∈ BV (Ω) and f ∈ C 1 (R), then f (u) ∈ BV (Ω).
• The space L 1 (0, T ; BV (Ω)) is the set of measurable functions u : (0, T ) → BV (Ω) such that

T 0 u(t) BV (Ω) dt < +∞. Equipped with the norm • L 1 (0,T ;BV (Ω)) := T 0 • BV (Ω) dt, it is a (nonreflexive) Banach space.
For more information on functions valued in Banach spaces, we refer for instance to [27, Chapter III].

Definition of entropy solutions

In the usual entropy formulations of scalar conservation laws, the admissible entropies are the C 1 convex functions or the so-called Kruzhkov entropies. Let us recall that the Kruzhkov entropies are the functions | • -κ| with κ ∈ R, the corresponding entropy fluxes being the functions

F (• ∨ κ) -F (• ∧ κ),
∈ C 0 (Q T × R, R d ) such that Φ(x, t, u) = 1 2 u a0 η (s) ∂F ∂u (x, t, s) ds + 1 2 u b0 η (s) ∂F ∂u (x, t, s) ds + 1 2 (η (a 0 )F (x, t, a 0 ) + η (b 0 )F (x, t, b 0 )) .
We then have ∂Φ ∂u (x, t, u) = η (u) ∂F ∂u (x, t, u), and we remark that, since the flux function F is divergence-free with respect to the space variables, the entropy flux is divergence-free in the same sense as well.

= |•-κ| with a 0 ≤ κ ≤ b 0 , then Φ(x, t, •) = F (x, t, • ∨ κ) -F (x, t, • ∧ κ).
This consistency property is used in the course of the proof of the uniqueness theorem. Moreover, for any convex function η ∈ C ∞ (R), letting Φ be given by Definition 2.1, integrate by parts in [a 0 , u] and [u, b 0 ] shows that the following relations hold for a.e. (x, t) ∈ Q T :

η(u) = 1 2 b0 a0 η (s)|u -s| ds + 1 2 ((η (a 0 ) + η (b 0 ))u + η(a 0 ) -η (a 0 )a 0 + η(b 0 ) -η (b 0 )b 0 ), Φ(x, t, u) = 1 2 b0 a0 η (s) (F (x, t, u ∨ s) -F (x, t, u ∧ s)) ds + η (a 0 ) + η (b 0 ) 2 F (x, t, u), ∀u ∈ [a 0 , b 0 ]. ( 5
)
The following definition results from the computations of Section 2.3, by the vanishing viscosity method.

Definition 2.3 (Entropy solution): Under Hypotheses (HC), a function u

∈ L ∞ (Q T ) ∩ L 1 (0, T ; BV (R d )) is said to be an entropy solution of (2)-(3) if there exists λ ∈ L ∞ (Q T ) d ,
ϕ ∈ C ∞ c (R d × [0, T )), Q T η(u)∂ t ϕ + Φ(x, t, u) -λη (u) • ∇ϕ dxdt - Q T ϕ D[η (u)] dt + R d η(u ini (x))ϕ(x, 0) dx ≥ 0. ( 6
) Since η is in C 1 (R), the function η (u) is in L 1 (0, T ; BV (R d )). Therefore, the term Q T ϕ D[η (u)
] dt is meaningful. The function λ, which is not necessarily unique, is called a multiplier by analogy with a Lagrange multiplier.

Formal derivation of the entropy formulation

In order to enlighten the link between the strong formulation and the entropy formulation, we present below a formal derivation of the entropy formulation by the vanishing viscosity method. Let us consider a viscous regularisation of the equation ( 2). We assume that, for all > 0, there exists (u , λ ) with

∂ t u + div F (x, t, u ) -div λ -∆u = 0 and λ ∈ Sgn(∇u ). ( 7
)
One can justify that this problem is well-posed and, owing to the term div λ , we can expect that u tends to u in L 1 loc (Q T ) when → 0. For a given admissible entropy-entropy flux pairs (η, Φ) in the sense of Definition 2.1, multiplying [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF] by η (u ), we find

∂ t η(u ) + div Φ(x, t, u ) -λ η (u ) + ∇η (u ) • λ -∆η(u ) + η (u )|∇u | 2 = 0. ( 8
)
The entropy η being convex, we have

∇η (u ) • λ = η (u )∇u • λ = η (u )|∇u | = |∇η (u )|.
Still by convexity of η, we have

η (u )|∇u | 2 ≥ 0.
Let us now multiply (8) by a nonnegative test function ϕ ∈ C ∞ c (R d × [0, T )) and integrate over Q T . We thus obtain

Q T η(u )∂ t ϕ + Φ(x, t, u ) -λ η (u ) • ∇ϕ -|∇η (u )|ϕ dxdt + R d η(u ini (x))ϕ(x, 0) dx ≥ 0. ( 9
)
Since u → u in L 1 loc (Q T ), it follows from the semi-continuity of the total variation that lim

→0 inf Q T |∇η (u )|ϕ dxdt ≥ Q T ϕ D[η (u)] dt.
Since the family (λ ) >0 is bounded, there exists λ ∈ L ∞ (Q T ) d such that, up to a subsequence,

λ λ weakly- * in L ∞ (Q T ) d .
Finally, letting → 0 in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], we obtain (6).

Existence and uniqueness

Theorem 2.4: Under Hypotheses (HC), there exists one and only one entropy solution of (2)-(3) in the sense of Definition 2.3.

Proof The existence of an entropy solution is proved by the convergence of a numerical approximation; see Theorem 6.2 below. To prove the uniqueness, we use the doubling variable technique.

Step 1. Let u be an entropy solution and λ u a corresponding multiplier. We consider the entropy η(• -κ), where η ∈ C ∞ (R) is an even convex function and κ ∈ R. Let Φ κ be the corresponding entropy flux. Then, by definition of an entropy solution, for all nonnegative test functions ϕ

∈ C ∞ c (R d × [0, T )), Q T η(u(x, t) -κ)∂ t ϕ + Φ κ (x, t, u(x, t)) -λ u (x, t)η (u(x, t) -κ) • ∇ x ϕ dxdt - Q T ϕ D x [η (u(x, t) -κ)] dt + R d η(u ini (x) -κ)ϕ(x, 0) dx ≥ 0. ( 10
)
Let v be another entropy solution and λ v a corresponding multiplier. Since η is an even function, η(

• -κ) = η(κ -•) and η (• -κ) = -η (κ -•).
Then, denoting by y and s the space and time variables, for all nonnegative

test functions ϕ ∈ C ∞ c (Q T ), Q T η(κ -v(y, s))∂ s ϕ + Φ κ (y, s, v(y, s)) + λ v (y, s)η (κ -v(y, s)) • ∇ y ϕ dyds - Q T ϕ D y [η (κ -v(y, s))] dt + R d η(u ini (y) -κ)ϕ(y, 0) dy ≥ 0. ( 11
)
Step 2. We now introduce well-chosen test functions. Let {ρ } >0 be a family of mollifiers in R d such that supp ρ ⊂ B(0, ) and {ρ } >0 be a family of mollifiers in R such that supp ρ ⊂ [-, 0]. Let r > 0, τ > 0, and

ψ ∈ C ∞ c (R, R + ). We define φ : Q T × Q T → R + such that φ(x, t, y, s) = ψ(t)ρ r (x -y)ρ τ (t -s), ∀(x, t, y, s) ∈ Q T × Q T . (12)
We take κ = v(y, s) and ϕ = φ(•, •, y, s) in [START_REF] Cockburn | Convergence of the finite volume method for multidimensional conservation laws[END_REF]. Next, integrating with respect to y and s over Q T , and noticing that φ(x, 0, y, s) = 0 for all s > 0, we obtain

Q T ×Q T η(u(x, t) -v(y, s))∂ t φ + Φ v(y,s) (x, t, u(x, t)) -λ u (x, t)η (u(x, t) -v(y, s)) • ∇ x φ dxdtdyds - Q T ×Q T φ D x [η (u(x, t) -v(y, s))] dtdyds ≥ 0. ( 13
)
Similarly, taking κ = u(x, t) and ϕ = φ(x, t, •, •) in [START_REF] Dean | On the numerical simulation of bingham visco-plastic flow: Old and new results[END_REF], then integrating with respect to x and t over Q T , we obtain

Q T ×Q T η(u(x, t) -v(y, s))∂ s φ + Φ u(x,t) (y, s, v(y, s)) + λ v (y, s)η (u(x, t) -v(y, s)) • ∇ y φ dxdtdyds - Q T ×Q T φ D y [η (u(x, t) -v(y, s))] dsdxdt + R d ×Q T η(u ini (y) -u(x, t))φ(x, t, y, 0)dydxdt ≥ 0. ( 14
)
Adding the above relations ( 13) and ( 14), we find

A 1 + A 2 + A 3 + A 4 + A 5 + A 6 + A 7 ≥ 0,
where

A 1 := Q T ×Q T η(u(x, t) -v(y, s))(∂ t φ + ∂ s φ) dxdtdyds, A 2 := Q T ×Q T Φ v(y,s) (x, t, u(x, t)) • ∇ x φdxdtdyds, A 3 := Q T ×Q T Φ u(x,t) (y, s, v(y, s)) • ∇ y φ dxdtdyds, A 4 := Q T ×Q T (λ u (x, t) -λ v (y, s))η (u(x, t) -v(y, s)) • (∇ x φ + ∇ y φ) dxdtdyds, A 5 := Q T ×Q T λ v (y, s)η (u(x, t) -v(y, s)) • ∇ x φ dxdtdyds - Q T ×Q T φ D x [η (u(x, t) -v(y, s))] dtdyds, A 6 := Q T ×Q T -λ u (x, t)η (u(x, t) -v(y, s)) • ∇ y φ dxdtdyds - Q T ×Q T φ D y [η (u(x, t) -v(y, s))] dsdxdt, A 7 := R d ×Q T η(u ini (y) -u(x, t))φ(x, t, y, 0) dydxdt.
Step 3. We now remark that we can get rid of A 4 , A 5 and A 6 , the terms arising from the total variation. The term A 4 vanishes owing to identity ∇ x φ + ∇ y φ = 0. The term A 5 is nonpositive. Indeed, integrating by parts with respect to x the first term of A 5 , we obtain

A 5 = Q T ×Q T -λ v (y, s)φ D x [η (u(x, t) -v(y, s))]dtdyds - Q T ×Q T φ D x [η (u(x, t) -v(y, s))] dtdyds. ( 15
)
Since |λ v (y, s)| ≤ 1 for all (y, s) ∈ Q T , we deduce that A 5 ≤ 0. With the same argument, we prove that A 6 ≤ 0.

Step 4. The resulting equation is now

A 1 + A 2 + A 3 + A 7 ≥ 0. (16)
Let κ ∈ R. Consider the sequence of entropies (η k ) k∈N * such that η k (x -κ) := 1/k + (x -κ) 2 . This sequence of entropies converges uniformly to the Kruzhkov entropy | • -κ|, and the entropy flux converges as well to the Kruzhkov entropy flux (see Remark 2.2). Passing to the limit in [START_REF] Eymard | Approximation by the finite volume method of an elliptic-parabolic equation arising in environmental studies[END_REF] with this sequence of entropies, we obtain,

A 10 + A 20 + A 30 + A 70 ≥ 0, (17)
where

A 10 := Q T ×Q T |u(x, t) -v(y, s)|(∂ t φ + ∂ s φ) dxdtdyds, A 20 := Q T ×Q T [F (x, t, u(x, t) ∨ v(y, s)) -F (x, t, u(x, t) ∧ v(y, s))] • ∇ x φ dxdtdyds, A 30 := Q T ×Q T [F (y, s, v(y, s) ∨ u(x, t)) -F (y, s, v(y, s) ∧ u(x, t))] • ∇ y φ dxdtdyds, A 70 := R d ×Q T |u ini (y) -u(x, t)|φ(x, t, y, 0) dydxdt.
Step 5. The remaining of the proof is identical to the one of Lemma 5 in [START_REF] Chainais-Hillairet | Finite volume schemes for a nonlinear hyperbolic equation. Convergence towards the entropy solution and error estimate[END_REF], since its starting point is precisely [START_REF] Eymard | A combined finite volume-nonconforming/mixed-hybrid finite element scheme for degenerate parabolic problems[END_REF].

3 Numerical approximation

Notation and hypotheses

The finite element mesh, denoted by T h , is a conforming simplicial mesh of R d of size h: T h is a (necessarily countable) set of disjoint open simplices such that K∈T h K = R d , and h is the maximum value of the diameter of all K ∈ T h . The mesh is conforming in the sense that, for two distinct elements K, L of T h , K ∩ L is either empty or a simplex included in an affine subset of R d with dimension strictly lower than d, whose vertices are simultaneously vertices of K and L. Therefore the set of the vertices of the mesh is countable as well, and denoted by {x p , p ∈ N}. In order to ensure the maximum principle, each element of T h is assumed to be nonobtuse; we recall that a simplex is said to be nonobtuse if the angles between any two facets are less than or equal to π/2. For any K ∈ T h , we denote by V K ⊂ N the set of the d + 1 indices of the vertices of K, and by E K the set of the d + 1 faces of K.

The finite volume mesh, denoted by D h , is a polyhedral mesh of R d such that the interface between two cells is a finite union of faces. The mesh D h is a dual mesh of T h in the sense that each cell of D h contains one and only one node of T h . For any p ∈ N, the cell of D h containing the node x p is denoted by Q p . We assume that

∀p ∈ N, Q p ⊂ K∈T h s.t. p∈V K K.
Let us introduce some additional notation about D h : N p is the set containing the indices of the neighbouring cells of Q p , E h is the set of couples (p, q) such that Q p and Q q are neighbours and p < q, σ p,q is the interface between two neighbour cells Q p and Q q , ν p,q is the unit normal vector to σ p,q pointing toward Q q , m p is the measure of Q p , m p,q is the measure of σ p,q .

Remark 3.1: Since a square can be divided into two right triangles and a cube can be divided into six nonobtuse tetrahedra, it is easy to build nonobtuse simplicial meshes of R 2 and R 3 . In fact, it is possible to generate nonobtuse simplicial meshes on any polygonal or polyhedral domain; see [START_REF] Brandts | On nonobtuse simplicial partitions[END_REF] and references therein.

Remark 3.2:

The dual mesh can be defined, for any d ≥ 1, by the following procedure. For any p ∈ N and any K ∈ T h with p ∈ V K , the set Q K p is defined as the set of all points x ∈ K under the form x = q∈V K λ q x q , with q∈V K λ q = 1 and λ p ≥ λ q for all q ∈ V K , and then define Q p as the union of all Q K p , such that p ∈ V K . In R 2 , for every triangle K ∈ T h containing the node x p , a part of the boundary of Q K p is the union of the segments joining the centre of gravity of K with the midpoint of the two edges of K incident to x p . This yields a closed polygonal line which delimits a cell Q p associated to x p .

In our scheme, the unknown function is simultaneously reconstructed from the values v h = (v p ) p∈N at the vertices using a continuous piecewise affine reconstruction denoted by v h , and using a piecewise constant reconstruction denoted by v h :

∀v h ∈ R N , v h ∈ C 0 (R d ) ; v h|K is affine for each K ∈ T h , v h (x p ) = v p , ∀p ∈ N, v h ∈ L 1 loc (R d ) ; v h|Qp = v p , ∀p ∈ N.
Letting K ∈ T h , if we denote by {x p } p∈V K the vertices of K and by {φ p } p∈V K the corresponding Lagrange basis function, we can write, for any

u h , v h ∈ R N , ∇ v h|K • ∇ u h|K = p∈V K q∈V K v p u q ∇φ p|K • ∇φ q|K .
Using the fact that p∈V K φ p|K = 1, and thus p∈V K ∇φ p|K = 0, the preceding equation can be rewritten as

∇ v h|K • ∇ u h|K = p∈V K q∈V K u q (v p -v q )∇φ p|K • ∇φ q|K .
Exchanging the roles of p and q, we get

∇ v h|K • ∇ u h|K = p∈V K q∈V K u p (v q -v p )∇φ p|K • ∇φ q|K .
Adding the two preceding relations provides

2∇ v h|K • ∇ u h|K = - p∈V K q∈V K (u p -u q )(v p -v q )∇φ p|K • ∇φ q|K ,
and therefore

K ∇ v h (x) • ∇ u h (x)dx = p∈V K q∈V K T K pq (u p -u q )(v p -v q ) with T K pq = - 1 2 K ∇φ p (x) • ∇φ q (x)dx. ( 18
)
Since the simplex K is nonobtuse, we have the standard inequality (see for example [START_REF] Brandts | On nonobtuse simplicial partitions[END_REF])

∇φ p|K • ∇φ q|K ≤ 0 and T K pq ≥ 0, ∀p, q ∈ V K , p = q. ( 19
)
Let us observe that, for any u h ∈ R N , we have, denoting for a.e.

x ∈ R d by p(x) ∈ N such that x ∈ Q p , for a.e. x ∈ R d , | u h (x) -u h (x)| = |(x -x p(x) ) • ∇ u h (x)| ≤ h|∇ u h (x)|, ( 20
) which implies, if ∇ u h ∈ L 1 (R d ) d , u h -u h L 1 (R d ) ≤ h ∇ u h L 1 (R d ) d . (21)
For the finite volume step, we need numerical fluxes F n p,q (u n p , u n q ) between two neighbouring cells p and q at time t n , function of u n p and u n q , the respective approximations of the unknown function in Q p and Q q at time t n . We require that the family of numerical fluxes (F n p,q ) p,q,n∈N is admissible and consistent with the flux F in the sense of the two definitions below. The Godunov scheme [24, see, e.g.,] provides, for instance, such numerical fluxes. Definition 3.3: A family of numerical fluxes (F n p,q ) p,q,n∈N is said to be admissible if

• F n p,q ∈ C 0 (R 2
) and there exists L > 0 such that, for all p ∈ N and q ∈ N p and for all n ∈ {0, ..., N -1}, the function F n p,q is Lipschitz continuous with the constant m p,q L with respect to each of its variables. • F n p,q is monotone, in the sense that it is non-decreasing with respect to its first argument and non-increasing with respect to its second argument,

• F n p,q is conservative, i.e. F n p,q (u, v) = -F n q,p (v, u) for all (u, v) ∈ R 2 .
Definition 3.4: Let F be a flux function. A family of numerical fluxes (F n p,q ) p,q,n∈N is said to be consistent with F if

F n p,q (u, u) = 1 δt t n+1 t n σp,q F (x, t, u) • ν p,q dγ(x)dt, ∀u ∈ R 2 . (22)
Let us give two examples of consistent and admissible families of numerical fluxes:

• the Godunov numerical flux, defined by:

F n p,q (u, v) =          min u≤s≤v 1 δt t n+1 t n σp,q F (x, t, s) • ν p,q dγ(x)dt if u ≤ v max v≤s≤u 1 δt t n+1 t n σp,q F (x, t, s) • ν p,q dγ(x)dt if v ≤ u,
• and the Rusanov numerical flux, defined, denoting by L F a Lipschitz constant of F , by:

F n p,q (u, v) = 1 δt t n+1 t n σp,q F (x, t, u + v 2 ) • ν p,q dγ(x)dt + m p,q L F 2 (u -v).
The following lemma, which is a discrete version of the divergence theorem, will be used below.

Lemma 3.5: Let (F n p,q ) p,q,n∈N be a family of numerical fluxes consistent with a flux function F in the sense of Definition

3.4. If F is divergence-free, then ∀u ∈ R, q∈Np F p,q (u, u) = 0. ( 23
)
Proof Owing to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] and to the divergence theorem, we have, for a given u ∈ R,

q∈Np F p,q (u, u) = 1 δt t n+1 t n q∈Np σp,q F (x, t, u) • ν p,q dγ(x)dt = 1 δt t n+1 t n p d i=1 ∂F i ∂x i (x, t, u) dxdt,
which vanishes owing to the assumption that F is divergence-free.

Remark 3.6: The consistency hypothesis for the numerical fluxes could be slightly weakened (see, for instance [8, Section 2.1]).

Description of the numerical scheme and well-posedness

We consider a family of discretisations (F h,δt ) h,δt>0 -by discretisation, we mean a finite element mesh T h , a finite volume mesh D h , a time step δt, a family of numerical fluxes (F n p,q ) p,q,n∈N . We assume that the following hypotheses, denoted in the following by Hypotheses (HD), are satisfied uniformly by any element F h,δt of the family.

(HD1) There exists α > 0 such that, for all p ∈ N, m p ≥ α h d , α |∂Q p | ≤ h d-1 , and |K| ≥ α h d , for all K ∈ T h . (HD2) There exists an admissible family of numerical fluxes (F n p,q ) p,q,n∈N in the sense of Definition 3.3, which is consistent with F in the sense of Definition 3.4. The constant L in Definition 3.3 is assumed to be independent of the discretisation.

(HD3) The time interval [0, T ] is divided into N equal intervals of length δt, such that the following CFL condition holds

δt ≤ α 2 h L . ( 24
)
Note that, thanks to Hypothesis (HD1), the condition [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] implies that

δt ≤ 1 L m p q∈Np m p,q , ∀p ∈ N. ( 25
)
The scheme for approximating (2)-( 3) is given by:

• Initialisation of (u 0 p ) n∈N such that u 0 h ∈ L ∞ (R d ) ∩ L 1 (R d ): u 0 p = 1 m p Qp u ini (x)dx, ∀p ∈ N. (26)
• Finite volume step. Letting (

u n p ) p∈N such that u n h ∈ L ∞ (R d ) ∩ L 1 (R d ), seek (u n+ 1 2 p ) p∈N such that u n+ 1 2 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and m p u n+ 1 2 p -u n p δt + q∈Np F n p,q (u n p , u n q ) = 0, ∀p ∈ N. (27)
• Finite element step. Let θ ∈ C 0 ((0, +∞)) be a positive function such that

lim h→0 θ(h) = lim h→0 h θ(h) = 0, (28)
and let us define

X h = {v h ∈ R N ; ∇ v h ∈ L 1 (R d ) ∩ L 2 (R d ), v h ∈ L 2 (R d )} Λ h := {µ h ∈ L ∞ (R d ) d ; µ h|K is constant for each K ∈ T h }.
Then the finite element step consists in seeking (u

n+1 h , λ n+1 h ) ∈ X h × Λ h such that R d u n+1 h -u n+ 1 2 h δt v h dx + R d (λ n+1 h + θ(h)∇ u n+1 h ) • ∇ v h dx = 0, ∀v h ∈ X h , (29)
λ n+1 h ∈ Sgn(∇ u n+1 h ). ( 30
)
An example of such a function θ is θ(h) = h γ with 0 < γ < 1. Note that, if d = 1, it is possible to let θ(h) = 0.

Remark 3.7: The explicit time discretisation of the hyperbolic step could be replaced by an implicit time discretisation. In this case, the CFL condition [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] on the time step would not be necessary.

For each discretisation F h,δt , we define the approximate solutions u h,δt :

Q T → R, u h,δt : Q T → R, and λ h,δt : Q T → R such that u h,δt (•, t) := u n+1 h if t ∈ (t n , t n+1 ], u h,δt (•, t) := u n+1 h if t ∈ (t n , t n+1 ], λ h,δt (•, t) := λ n+1 h if t ∈ (t n , t n+1 ].
The proposition below proves that the scheme has at least one solution, which is unique with respect to the unknown u h .

Proposition 3.8: Let us assume Hypotheses (HC) and Hypotheses (HD). Then there exists at least one solution to Scheme ( 26)-( 30) such that

u 0 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and, for all n ∈ N, u n+ 1 2 h , u n+1 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and (u n+1 h , λ n+1 h ) ∈ X h × Λ h . Moreover, u 0 h and, for all n ∈ N, u n+ 1 2 h
and u n+1 h are unique and, for all p ∈ N,

a 0 ≤ u n+ 1 2 p ≤ b 0 and a 0 ≤ u n p ≤ b 0 .
Proof Thanks to [START_REF] Serre | Systems of conservation laws[END_REF], we get from Hypothesis (HC1) that

u 0 h L 1 (R d ) ≤ u ini L 1 (R d ) ,
and a 0 ≤ u 0 p ≤ b 0 for all p ∈ N, which completes the proof that

u 0 h ∈ L ∞ (R d ) ∩ L 1 (R d ).
For any n ∈ N, we get from Propositions 3.9 and 3.11, assuming

u n h ∈ L ∞ (R d ) ∩ L 1 (R d ) with a 0 ≤ u n p ≤ b 0 for all p ∈ N, that u n+ 1 2 h ∈ L ∞ (R d ) ∩ L 1 (R d ) with a 0 ≤ u n+ 1 2 p
≤ b 0 for all p ∈ N. Using Proposition 3.12, we get the existence of (u n+1 h , λ n+1 h ) ∈ X h × Λ h such that (30) holds, and we get that u n+1 h is unique. It now suffices to apply Proposition 3.14, for obtaining that

u n+1 h ∈ L ∞ (R d ) ∩ L 1 (R d ) with a 0 ≤ u n+1 p ≤ b 0 for all p ∈ N.
Let us now state and prove the propositions used in the proof of the preceding result. We have first the following result.

Proposition 3.9: Let us assume Hypotheses (HC) and Hypotheses (HD). Let n ∈ N, κ ∈ R and a family (u n p ) p∈N be given such that a 0 ≤ u n p ≤ b 0 for all p ∈ N. Let (u

n+ 1 2 p
) p∈N be given by [START_REF] Showalter | Monotone operators in Banach space and nonlinear partial differential equations[END_REF]. Then there holds

u n+ 1 2 p ∨ κ ≤ u n p ∨ κ - δt m p q∈Np F n p,q (u n p ∨ κ, u n q ∨ κ), ∀p ∈ N, (31)
u n+ 1 2 p ∧ κ ≥ u n p ∧ κ - δt m p q∈Np F n p,q (u n p ∧ κ, u n q ∧ κ), ∀p ∈ N, (32)
m p |u n+ 1 2 p -κ| -|u n p -κ| δt + q∈Np F n p,q (u n p ∨ κ, u n q ∨ κ) -F n p,q (u n p ∧ κ, u n q ∧ κ) ≤ 0, ∀p ∈ N. (33)
Consequently a 0 ≤ u

n+ 1 2 p ≤ b 0 for all p ∈ N.
Proof The proof of this proposition is done in [8, Lemma 3] or [START_REF] Eymard | Finite volume methods[END_REF]Lemma 27.1]. We recall it since it is very brief. We consider the function H n p : R 1+#Np → R, defined by

H n p (a, (b q ) q∈Np ) = a - δt m p q∈Np F n p,q (a, b q ).
We observe that, for any a > a, there holds

H n p (a , (b q ) q∈Np ) -H n p (a, (b q ) q∈Np ) ≥ (a -a)(1 - δtL q∈Np m p,q m p ) ≥ 0,
thanks to Definition 3.3 of admissible fluxes and to condition [START_REF] Lions | Remarks on some nonlinear evolution problems arising in Bingham flows[END_REF] implied by [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. Therefore the function H n p is non-decreasing with respect to all its arguments. Noticing that κ = H n p (κ, (κ) q∈Np ) and u

n+ 1 2 p = H n p (u n p , (u n q ) q∈Np ), we get that κ ≤ H n p (u n p ∨ κ, (u n q ∨ κ) q∈Np ) and u n+ 1 2 p ≤ H n p (u n p ∨ κ, (u n q ∨ κ) q∈Np ), which
implies (31). The proof of (32) is similar, and (33) is obtained by the difference between (31) and (32). Letting κ = b 0 in (31) and using [START_REF] Hintermüller | Total bounded variation regularization as a bilaterally constrained optimization problem[END_REF] on one hand, letting κ = a 0 in (32) on the other hand complete the proof that a 0 ≤ u

n+ 1 2 p
≤ b 0 for all p ∈ N. We then deduce the following result. Proposition 3.10: Let us assume Hypotheses (HC) and Hypotheses (HD). Let (η, Φ) be an entropy-entropy flux pair in the sense of Definition 2.1, and let n ∈ N be given. Then, the family (Φ n p,q ) p,q,n∈N of admissible numerical fluxes defined by

Φ n p,q (x, y) := 1 2 b0 a0 η (κ) F n p,q (x ∨ κ, y ∨ κ) -F n p,q (x ∧ κ, y ∧ κ) dκ + η (a 0 ) + η (b 0 ) 2 F n p,q (x, y), (34)
is consistent with Φ in the sense of Definition 3.4, and is such that, if u

n+ 1 2 h ∈ R N is obtained from u n h ∈ R N by (27) with a 0 ≤ u n p ≤ b 0 , for all p ∈ N, then m p η(u n+ 1 2 p ) -η(u n p ) δt + q∈Np Φ n p,q (u n p , u n q ) ≤ 0, ∀n ∈ {0, ..., N -1}, ∀p ∈ N. (35)
Furthermore, there is a constant L , depending only on L, η, a 0 and b 0 , such that, for all (p, q) ∈ E h and for all n ∈ {0, ..., N -1}, the function Φ n p,q is Lipschitz continuous with respect of each of its variables with the constant m p,q L . Proof Thanks to Proposition 3.9, we have that all values u n p and u

n+ 1 2 p
belong to [a 0 , b 0 ], which enables to use relations [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF]. The consistency of Φ n p,q with Φ is a consequence of ( 5) and of the consistency of F n p,q with F . Multiplying (33) by η (κ) and integrating on κ ∈ [a 0 , b 0 ] implies (35), owing to ( 5) and (34). From the above result, we get the following one. Proposition 3.11: Let us assume Hypotheses (HC) and Hypotheses (HD). Assume that, for a given n ∈ N,

a 0 ≤ u n p ≤ b 0 for all p ∈ N and u n h ∈ L 1 (R d ).
Then for all entropy η in the sense of Definition 2.1 such that η(0) = 0 and η (0) = 0,

R d η(u n+ 1 2 h (x))dx ≤ R d η(u n h (x))dx, (36)
and consequently u

n+ 1 2 h L 1 (R d ) ≤ u n h L 1 (R d ) and u n+ 1 2 h L 2 (R d ) ≤ u n h L 2 (R d ) . (37)
Proof We get, from (35) and applying Lemma 3.5,

m p η(u n+ 1 2 p ) -η(u n p ) δt + q∈Np (Φ n p,q (u n p , u n q ) -Φ n p,q (0, 0)) ≤ 0, ∀n ∈ {0, ..., N -1}, ∀p ∈ N,
For ε > 0, we multiply the preceding inequality by δt exp(-ε|x p |) and we sum the result on p ∈ N. We get p∈N m p (η(u

n+ 1 2 p ) -η(u n p )) exp(-ε|x p |) + p∈N q∈Np Φ n p,q (u n p , u n q ) -Φ n p,q (0, 0) exp(-ε|x p |) ≤ 0.
Defining, for any (p, q) ∈ E h , the point x pq = 1 2 (x p + x q ) and using the property Φ n p,q (u, v) = -Φ n q,p (v, u) for all (u, v) ∈ [a 0 , b 0 ] 2 (which therefore implies Φ n p,q (u n p , u n q ) -Φ n p,q (0, 0) + Φ n q,p (u n q , u n p ) -Φ n q,p (0, 0) = 0), we get

Φ n p,q (u n p , u n q ) -Φ n p,q (0, 0) exp(-ε|x p |) + Φ n q,p (u n q , u n p ) -Φ n q,p (0, 0) exp(-ε|x q |) = A n pq + A n qp ,

with

A n pq = Φ n p,q (u n p , u n q ) -Φ n p,q (0, 0) (exp(-ε|x p |) -exp(-ε|x pq |)). Hence the two preceding relations imply p∈N m p (η(u

n+ 1 2 p ) -η(u n p )) exp(-ε|x p |) + p∈N q∈Np A n pq ≤ 0,
Observing that Proposition 3.10 implies

|Φ n p,q (u n p , u n q ) -Φ n p,q (0, 0)| ≤ L m pq (|u n p | + |u n q |), we get, using | exp(-ε|x p |) -exp(-ε|x pq |)| ≤ ε||x p | -|x pq || ≤ ε 1 2 |x p -x q |, that A n pq ≥ - L 2 m pq (|u n p | + |u n q |)εh.
This leads to p∈N m p η(u

n+ 1 2 p ) exp(-ε|x p |) ≤ p∈N m p η(u n p ) exp(-ε|x p |) + ε L 2 p∈N q∈Np m pq h(|u n p | + |u n q |).
Since we have 1 2

p∈N q∈Np

m pq h(|u n p | + |u n q |) = p∈N q∈Np m pq h|u n p | = p∈N |u n p | |∂Q p |h,
we can write, using Hypothesis (HD1), p∈N m p η(u

n+ 1 2 p ) exp(-ε|x p |) ≤ p∈N m p η(u n p ) + ε L α 2 p∈N m p |u n p |.
Letting ε → 0, we get by monotonous convergence (recall that thanks to the hypothesis

u n h ∈ L 1 (R d ), we also have η(u n h ) ∈ L 1 (R d ) and η(s) ≥ η(0) = 0), p∈N m p η(u n+ 1 2 p ) ≤ p∈N m p η(u n p ),
which concludes the proof of (36).

Letting η tend to η(s) = |s| and letting η(s) = s 2 allow for concluding (37).

The proposition below proves that the finite element step is well-posed, provided that u

n+ 1 2 h ∈ L 2 (R d ), which holds if u n+ 1 2 h ∈ L ∞ (R d ) ∩ L 1 (R d )
, and gives a variational characterisation of u n+1 h . Proposition 3.12: Let us assume Hypotheses (HC) and Hypotheses (HD). Let n ∈ N be given and let us assume that u

n+ 1 2 h ∈ L 2 (R d ). Then, equations (29)-(30) admit a solution (u n+1 h , λ n+1 h ) ∈ X h × Λ h . Furthermore, u n+1
h is unique and is the minimiser of the functional J n+1 h : X h → R defined by

J n+1 h (v h ) := 1 2δt R d v h -u n+ 1 2 h 2 dx + R d (|∇ v h | + 1 2 θ(h)|∇ v h | 2 ) dx. ( 38
)
Proof We remark that, defining the following norm on the Banach space X h ,

v h X h = v h L 2 (R d ) + ∇ v h L 1 (R d ) + ∇ v h L 2 (R d ) ,
we obtain that the stricly convex function

J n+1 h is such that lim v h X h →∞ J n+1 h (v h ) = +∞.
Hence J n+1 h reaches its unique minimum value at some point u n+1 h ∈ X h . Let us prove that this point is characterised by ( 29)- [START_REF] Ziemer | Weakly differentiable functions[END_REF].

Let us first assume that there exists (u n+1 h , λ n+1 h ) ∈ X h × Λ h , solution to ( 29)- [START_REF] Ziemer | Weakly differentiable functions[END_REF]. Writing, for any

w h ∈ X h , J n+1 h (w h ) = J n+1 h (u n+1 h + v h ) with v h = w h -u n+1 h
, we have using (29)

J n+1 h (w h ) -J n+1 h (u n+1 h ) = 1 2δt R d v h (x) 2 dx + A + R d 1 2 θ(h)|∇ v h (x)| 2 dx, with A = R d (|∇( u n+1 h + v h )(x)| -|∇ u h (x)| -λ n+1 h (x) • ∇ v h (x))dx = K∈T h |K|(|∇( u n+1 h + v h ) |K | -|(∇ u n+1 h ) |K | -λ n+1 K • (∇ v h ) |K ).
Recall that, if (∇ u n+1 h

) |K = 0, then λ n+1 K = (∇ u n+1 h ) |K |(∇ u n+1
: X h → R by L n+1 h (v h ) := 1 δt R d u n+1 h -u n+ 1 2 h v h dx + {∇ u n+1 h =0} ∇ u n+1 h |∇ u n+1 h | • ∇ v h dx + θ(h) R d ∇ u n+1 h • ∇ v h dx.
For a given ε ∈ R with ε = 0 and for any v h ∈ X h , we write that

J n+1 h (u n+1 h + εv h ) -J n+1 h (u n+1 h ) ≥ 0.
Assuming ε > 0, dividing the above equation by ε and letting ε → 0, we get

L n+1 h (v h ) + {∇ u n+1 h =0} |∇ v h |dx ≥ 0.
Assuming ε < 0, dividing the above equation by ε and letting ε → 0, we get

L n+1 h (v h ) - {∇ u n+1 h =0} |∇ v h |dx ≤ 0.
Hence we get that

∀v h ∈ X h , |L n+1 h (v h )| ≤ {∇ u n+1 h =0} |∇ v h |dx = K∈T n+1 h,0 |K| |∇ v h|K |. (39)
We consider the set E of all functions f from T n+1 h,0 → R d which are bounded for the norm

f E = K∈T n+1 h,0 |K| |f K |.
We observe that, defining

F = {f ∈ E, ∃v h ∈ X h , f = ∇ v h }, the linear form B : F → R, such that B(f ) = L n+1 h (v h ) is well defined (since if ∇ v h = ∇ w h on all K ∈ T n+1 h,0 , then we get from (39) that L n+1 h (v h ) = L n+1 h (w h ))
and continuous (again from (39), which proves that

B F := sup f ∈F \{0} B(f ) f E ≤ 1)
. Applying the Hahn-Banach theorem, this linear form can be extended on E by a linear form, again denoted by B, with the same norm (hence lower or equal to 1). Hence there exists (

λ n+1 K ) K∈T n+1 h,0 with ∀f ∈ E, B(f ) = - K∈T n+1 h,0 |K| λ n+1 K • f K .
and 29)- [START_REF] Ziemer | Weakly differentiable functions[END_REF] holds.

B E = sup K∈T n+1 h,0 |λ n+1 K | ≤ 1. Therefore we have ∀v h ∈ X h , L n+1 h (v h ) + K∈T n+1 h,0 |K| λ n+1 K • (∇ v h ) |K = 0, which concludes, denoting λ n+1 K = (∇ u n+1 h ) |K |(∇ u n+1 h ) |K | if (∇ u n+1 h ) |K = 0, the proof that λ n+1 h ∈ Λ h is such that (
Remark 3.13: Solving the saddle-point problem ( 29)- [START_REF] Ziemer | Weakly differentiable functions[END_REF] or, which is equivalent, minimising the functional J n+1 h is not an easy task. For instance, if one tries to minimise J n+1 h (or more precisely a regularisation of J n+1 h , since J n+1 h is nondifferentiable) with conventional algorithms, such as the gradient or Newton algorithms, one obtains very slow convergences. Motivated by the applications in image denoising where this kind of problem arises, specific algorithms have been developed since the 1990s. Let us mention in particular the lagged diffusivity algorithm [START_REF] Vogel | Iterative methods for total variation denoising[END_REF], saddle-point algorithms [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Bartels | Total variation minimization with finite elements: Convergence and interative solution[END_REF] or semi-smooth Newton methods [START_REF] Hintermüller | Total bounded variation regularization as a bilaterally constrained optimization problem[END_REF]. Proposition 3.14: Let us assume Hypotheses (HC) and Hypotheses (HD). Let n ∈ N. Let η ∈ C 2 (R) such that η(0) = 0, η (0) = 0 and there exists M ∈ R + with η (s) ∈ [0, M ] for all s ∈ R. Assume that u

n+ 1 2 h ∈ R N is such that p∈N m p η(u n+ 1 2 p ) < ∞. Then p∈N m p η(u n+1 p ) ≤ p∈N m p η(u n+ 1 2 p ). ( 40
)
As a consequence, if a 0 ≤ u

n+ 1 2 p ≤ b 0 for all p ∈ N, then a 0 ≤ u n+1 p ≤ b 0 for all p ∈ N, and if u n+ 1 2 h ∈ L 1 (R d ) then u n+1 h ∈ L 1 (R d ).
Proof We remark that v h , defined by v p = η (u n+1 p ), is such that v h ∈ X h . Indeed, there holds |v p | ≤ M |u n+1 p |, and, using ( 18)-( 19), we have

∇ v h 2 L 2 (R) d = K∈T h p∈V K q∈V K T K pq (η (u n+1 p ) -η (u n+1 q )) 2 ≤ M 2 K∈T h p∈V K q∈V K T K pq (u n+1 p -u n+1 q ) 2 = M 2 ∇ u n+1 h 2 L 2 (R) d ,
which implies v h ∈ X h . We now remark that, for any K ∈ T h , if ∇ u n+1 h|K = 0, then u n+1 p = u n+1 q for any p, q ∈ V K , and therefore ∇ v h|K = 0. So one can write

λ n+1 h|K • ∇ v h|K = 0 = ∇ u n+1 h|K • ∇ v h|K .
If ∇ u n+1 h|K = 0, we then have

λ n+1 h|K • ∇ v h|K = 1 |∇ u n+1 h|K | ∇ u n+1 h|K • ∇ v h|K .
Hence, defining α n+1 K by α n+1 K = 1 if ∇ u n+1 h|K = 0, and

α n+1 K = 1 |∇ u n+1 h|K | otherwise, we get R d (λ n+1 h + θ(h)∇ u n+1 h ) • ∇ v h dx = K∈T h (α n+1 K + θ(h)) p∈V K q∈V K T K pq (u n+1 p -u n+1 q )(η (u n+1 p ) -η (u n+1 q )) ≥ 0.
Hence we can write from ( 29)

p∈N m p (u n+1 p -u n+ 1 2 p )η (u n+1 p ) ≤ 0. Applying η(b) -η(a) = η (a)(b -a) + η (c) (b-a) 2 2
for c between a and b, we get

p∈N m p (η(u n+1 p ) -η(u n+ 1 2 p )) ≤ p∈N m p (u n+1 p -u n+ 1 2 p )η (u n+1 p ) ≤ 0,
which concludes the proof of (40).

Then, assuming u

A priori estimates on the approximate solutions

A L ∞ (Q T ) estimate on the family of approximate velocities (u h,δt ) h,δt>0 has already been proved in Proposition 3.8. The aim of this section is to establish additional estimates on ( u h,δt ) h,δt>0 , namely a L 1 (0, T ; BV (R d )) estimate, and L 1 loc (Q T ) estimates on the space and time translates. The estimates on the space and time translates are deduced from the L 1 (0, T ; BV (R d )) estimate.

Remark 4.1: Hypothesis (HD1) implies that each cell of D h has a finite number of neighbours (and this number is bounded independently of F h,δt ). 26)- [START_REF] Ziemer | Weakly differentiable functions[END_REF]. Then there holds

u 0 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and, for all n ∈ N, u n+ 1 2 h , u n+1 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and (u n+1 h , λ n+1 h ) ∈ X h × Λ h be a solution to Scheme (
1 2 u h,δt 2
L ∞ (0,T ;L 2 (R d )) + N n=1 δt R d |∇ u n h | dx ≤ 1 2 u ini 2 L 2 (R d ) , (41)
and

N -1 n=0 δt R d θ(h)|∇ u n+1 h | 2 dx = θ(h) N -1 n=0 δt K∈T h p∈V K q∈V K T K pq (u n+1 p -u n+1 q ) 2 ≤ 1 2 u ini 2 L 2 (R d ) . (42)
Proof Owing to (37) in Proposition 3.11, we have

R d u n+ 1 2 h (x) 2 dx ≤ R d u n h (x) 2 dx. (43)
We test [START_REF] Vogel | Iterative methods for total variation denoising[END_REF] with v h = δt u n+1 This leads to

p∈N m p 1 2 (u n+1 p ) 2 + 1 2 (u n+1 p -u n+ 1 2 p ) 2 - 1 2 (u n+ 1 2 p ) 2 + δt R d (|∇ u n+1 h | + θ(h)|∇ u n+1 h | 2 ) dx = 0, hence giving, thanks to (43) p∈N m p 1 2 (u n+1 p ) 2 - 1 2 (u n p ) 2 + δt R d (|∇ u n+1 h | + θ(h)|∇ u n+1 h | 2 ) dx ≤ 0. ( 44
)
Summing (44) over n ∈ {0, ..., m} for any m = 1, . . . , N -1 we obtain the inequality

1 2 p∈N m p (u m+1 p ) 2 + m n=0 δt R d (|∇ u n+1 h | + θ(h)|∇ u n+1 h | 2 ) dx ≤ p∈N m p 1 2 (u 0 p ) 2 . ( 45
)
We conclude the proof of the lemma since the above inequality holds for any m.

Time translate estimate

Proposition 4.4: Let us assume Hypotheses (HC) and Hypotheses (HD). Let 26)- [START_REF] Ziemer | Weakly differentiable functions[END_REF]. For all R > 0, there exists a constant C, independent of the family (F h,δt ) h,δt>0 , such that,

u 0 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and, for all n ∈ N, u n+ 1 2 h , u n+1 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and (u n+1 h , λ n+1 h ) ∈ X h × Λ h be a solution to Scheme (
T -s 0 B(0,R) |u h,δt (x, t + s) -u h,δt (x, t)| dxdt ≤ C √ s, ∀s ∈ [0, T ]. (46) Proof 1. Let s ∈ [0, T ].
Using the Cauchy-Schwarz inequality, we have

T -s 0 B(0,R) |u h,δt (x, t + s) -u h,δt (x, t)| dxdt ≤ CR d T T -s 0 R d (u h,δt (x, t + s) -u h,δt (x, t)) 2 dxdt 1 2
.

(47) We define the function ν : R → Z such that ν(t) = n + 1 if t ∈ (nδt, (n + 1)δt] and the function

χ n (t, s) : R × R → {0, 1} such that χ n (t, s) = 1 if ν(t) ≤ n < ν(t + s) and χ n (t, s) = 0 otherwise; so that T -s 0 R d (u h,δt (x, t + s) -u h,δt (x, t)) 2 dxdt = T -s 0 R d u ν(t+s) h -u ν(t) h 2 dxdt = T -s 0 R d u ν(t+s) h -u ν(t) h N -1 n=1 χ n (t, s) u n+1 h -u n h dxdt. ( 48
)
2. We denote τ = 0 or τ = s, and we multiply [START_REF] Showalter | Monotone operators in Banach space and nonlinear partial differential equations[END_REF], for n ≥ 1, by δt u ν(t+τ ) p

and sum over all p ∈ N:

p∈N m p u ν(t+τ ) p u n+ 1 2 p -u n p + δt p∈N q∈Np u ν(t+τ ) p F n p,q (u n p , u n q ) -F n p,q (u n p , u n p ) = 0. ( 49
)
Testing [START_REF] Vogel | Iterative methods for total variation denoising[END_REF] with v h = δtu ν(t+τ ) h

, we get:

p∈N m p u ν(t+τ ) p u n+1 p -u n+ 1 2 p + δt R d (λ n+1 h • ∇ u ν(t+τ ) h + θ(h)∇ u n+1 h • ∇ u ν(t+τ ) h ) dx = 0. ( 50
)
Adding the above equalities (49) and (50), we find

p∈N m p u ν(t+τ ) p u n+1 p -u n p + δt p∈N u ν(t+τ ) p q∈Np F n p,q (u n p , u n q ) -F n p,q (u n p , u n p ) + δt R d (λ n+1 h • ∇ u ν(t+τ ) h + θ(h)∇ u n+1 h • ∇ u ν(t+τ ) h ) dx = 0. ( 51
)
3. Let p ∈ N, q ∈ N p and K an element of T h for which x p and x q are vertices. In view of (HD1) and (HD2), for n ≥ 1,

F n p,q (u n p , u n q ) -F n p,q (u n p , u n p ) ≤ 2m p,q L|u n p -u n q | ≤ 2m p,q Lh|∇ u n h|K | ≤ Ch d |∇ u n h|K |. ( 52
)
Therefore, owing to (HD1), to the finite number of neighbours in D h and to the L ∞ estimate stated in Proposition 3.8,

p∈N u ν(t+τ ) p q∈Np F n p,q (u n p , u n q ) -F n p,q (u n p , u n p ) ≤ C R d |∇ u n h | dx. ( 53
)
The property

|λ n+1 h | ≤ 1 gives immediately R d λ n+1 h • ∇ u ν(t+τ ) h dx ≤ R d |∇ u ν(t+τ ) h | dx. ( 54
)
From ( 51), ( 53) and ( 54), it follows

R d u ν(t+τ ) h u n+1 h -u n h dx ≤ Cδt R d (|∇ u n h | + |∇ u ν(t+τ ) h | + θ(h) 2 (|∇ u n+1 h | 2 + |∇ u ν(t+τ ) h | 2 )) dx . (55)
4. Using the above estimate (55), we have

T -s 0 R d u ν(t+τ ) h N -1 n=1 χ n (t, s) u n+1 h -u n h dxdt ≤ C N n=1 T -s 0 χ n (t, s)δt(A n + B ν(t+τ ) ) dt,
with, for all n = 1, . . . , N ,

A n = R d (|∇ u n h | + θ(h) 2 |∇ u n+1 h | 2 ) dx
and

B n = R d (|∇ u n h | + θ(h) 2 |∇ u n h | 2 ) dx.
We then apply Lemma 4.5 below, which leads to

T -s 0 R d u ν(t+τ ) h N -1 n=1 χ n (t, s) u n+1 h -u n h dxdt ≤ Cs N n=1 δt(A n + B n ). ( 56
)
5. Recalling (41) and ( 42) which provide bounds on N n=1 δt(A n + B n ), and collecting ( 48) and ( 56), we obtain the desired estimate (46). Let us now state a lemma proved in [START_REF] Eymard | Approximation by the finite volume method of an elliptic-parabolic equation arising in environmental studies[END_REF]Lemma 4.6] and used in the preceding proof. 26)- [START_REF] Ziemer | Weakly differentiable functions[END_REF]. There exists a constant C, independent of the family (F h,δt ) h,δt>0 , such that

Space translate estimate

0 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and, for all n ∈ N, u n+ 1 2 h , u n+1 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and (u n+1 h , λ n+1 h ) ∈ X h × Λ h be a solution to Scheme (
T 0 R d | u h,δt (x + y, t) -u h,δt (x, t)| dxdt ≤ C|y|, ∀y ∈ R d . ( 57
)
and, for any R > 0, there exists C , independent of the family (F h,δt ) h,δt>0 , we find the inequality

T 0 B(0,R) |u h,δt (x + y, t) -u h,δt (x, t)| dxdt ≤ C(|y| + 2CT R d h), ∀y ∈ R d . ( 58
| u n h (x + y) -u n h (x)| ≤ K∈T h χ K (x, x + y)|∇ u n h|K |. ( 59
)
Next, an integration with respect to x yields

R d | u n h (x + y) -u n h (x)| dx ≤ K∈T h |∇ u n h|K | R d χ K (x, x + y) dx. ( 60
)
For any K ∈ T h , the function x → χ K (x, x + y) is bounded by min(h, |y|) and is zero outside a domain of measure lower than h d-1 (h + |y|). Therefore,

R d χ K (x, x + y) dx ≤ h d-1 (h + |y|) min(h, |y|) ≤ 2h d |y|, (61)
and

R d | u n h (x + y) -u n h (x)| dx ≤ C|y| K∈T h h d |∇ u n h|K |. ( 62
)
Using the hypothesis (HD1), we obtain

K∈T h h d |∇ u n h|K | ≤ C R d |∇ u n h | dx. ( 63
)
Summing the above inequality (63) over {1, ..., N }, and using (41), we find the desired estimate (57). We then deduce (58) using [START_REF] Feng | Rate of convergence of regularization procedures and finite element approximations for the total variation flow[END_REF].

Entropy formulation for the approximate solutions

The aim of this section is to establish an entropy formulation, similar to ( 6), for the approximate solutions. We first prove a discrete entropy inequality for the finite volume step (Proposition 3.10). Then, we take into account the finite element step to obtain the complete discrete entropy formulation (Proposition 5.1). Error terms occur in this formulation. Proposition 5.1 ensures that they tend to zero when the meshsize and the time step tend to zero.

Proposition 5.1: Let us assume Hypotheses (HC) and Hypotheses (HD). Let

u 0 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and, for all n ∈ N, u n+ 1 2 h , u n+1 h ∈ L ∞ (R d ) ∩ L 1 (R d ) and (u n+1 h , λ n+1 h
) ∈ X h × Λ h be a solution to Scheme (26)- [START_REF] Ziemer | Weakly differentiable functions[END_REF]. Let (η, Φ) be an entropy-entropy flux pair. Then, for all nonnegative test functions ϕ ∈ C ∞ c (R d × [0, T )), there holds

Q T η(u h,δt )∂ t ϕ dxdt + Q T Φ(x, t, u h,δt ) • ∇ϕ dxdt - Q T η ( u h,δt )λ h,δt • ∇ϕ dxdt - Q T ϕ|∇η ( u h,δt )| dxdt + R d η(u 0 h )ϕ(x, 0) dx + e h,δt ≥ 0, (64)
where e h,δt satisfies lim

h→0,δt→0 e h,δt = 0. ( 65
)
Proof First step: proof of (64)

1. Let ϕ n+ 1 2 ∈ C ∞ c (R d , R + ) be defined by ϕ n+ 1 2 (x) := 1 δt t n+1 t n
ϕ(x, s) ds. We test [START_REF] Vogel | Iterative methods for total variation denoising[END_REF] with the function

v n+1 h ∈ X h such that v n+1 p = η (u n+1 p )ϕ n+ 1 2 (x p ), for all p ∈ N. We get: p∈N m p u n+1 p -u n+ 1 2 p δt η (u n+1 p )ϕ n+ 1 2 (x p ) + R d (λ n+1 h • ∇ v n+1 h + θ(h)∇ u n+1 h • ∇ v n+1 h ) dx = 0. ( 66
) Since η is convex, u n+1 p -u n+ 1 2 p η (u n+1 p ) ≥ η(u n+1 p ) -η(u n+ 1 2 p ) ∀p ∈ N,
and thus (66) leads to

p∈N m p η(u n+1 p ) -η(u n+ 1 2 p ) δt ϕ n+ 1 2 (x p ) + R d (λ n+1 h • ∇ v n+1 h + θ(h)∇ u n+1 h • ∇ v n+1 h ) dx ≤ 0. ( 67
)
Defining

v n+1 := η ( u n+1 h )ϕ n+ 1 2
, the second term of (67) can be rewritten as

R d λ n+1 h • ∇ v n+1 h dx = R d λ n+1 h • ∇ v n+1 dx + R d λ n+1 h • ∇ v n+1 h -v n+1 dx = R d λ n+1 h • ∇ η ( u n+1 h ) ϕ n+ 1 2 dx + R d λ n+1 h • ∇ϕ n+ 1 2 η ( u n+1 h ) dx + R d λ n+1 h • ∇ v n+1 h -v n+1 dx. (68)
We have

R d λ n+1 h • ∇ η ( u n+1 h ) ϕ n+ 1 2 dx = R d η ( u n+1 h )λ n+1 h • ∇( u n+1 h )ϕ n+ 1 2 dx = R d |∇η ( u n+1 h )|ϕ n+ 1 2 dx. ( 69
)
The third term of (67) can be rewritten as

R d θ(h)∇ u n+1 h •∇ v n+1 h dx = θ(h) K∈T h p∈V K q∈V K T K pq (u n+1 p -u n+1 q )(η (u n+1 p )ϕ n+ 1 2 (x p )-η (u n+1 q )ϕ n+ 1 2 (x q )) = T 1 +T 2 ,
with

T 1 = θ(h) K∈T h p∈V K q∈V K T K pq (u n+1 p -u n+1 q )(η (u n+1 p ) -η (u n+1 q )) ϕ n+ 1 2 (x p ) + ϕ n+ 1 2 (x q ) 2 ≥ 0 T 2 = θ(h) K∈T h p∈V K q∈V K T K pq (u n+1 p -u n+1 q )(ϕ n+ 1 2 (x p ) -ϕ n+ 1 2 (x q )) η (u n+1 p ) + η (u n+1 q ) 2 . (70)
Collecting ( 67), ( 68), ( 70) and (69), we obtain

p∈N m p η(u n+1 p ) -η(u n+ 1 2 p ) δt ϕ n+ 1 2 (x p ) + R d λ n+1 h • ∇ϕ n+ 1 2 η ( u n+1 h ) dx + R d |∇η ( u n+1 h )|ϕ n+ 1 2 dx ≤ R d λ n+1 h • ∇ v n+1 -v n+1 h dx + θ(h) p∈N q∈N T pq |u n+1 p -u n+1 q | |ϕ n+ 1 2 (x p ) -ϕ n+ 1 2 (x q )| η (u n+1 p ) + η (u n+1 q ) 2 . (71)
Multiplying by δt and summing over n ∈ {0, ..., N -1}, we eventually find

N -1 n=0 p∈N η(u n+1 p ) -η(u n+ 1 2 p ) δt t n+1 t n m p ϕ(x p , t) dt + Q T η ( u h,δt )λ h,δt • ∇ϕ dxdt + Q T ϕ|∇η ( u h,δt )| dxdt ≤ e (1) h,δt + e (2) h,δt , (72) where e (1)
h,δt := N -1 n=0 δt R d λ n+1 h (x) • ∇( v n+1 (x) -v n+1 h (x)) dx, (73) and e (2)
h,δt :=

N -1 n=0 δtθ(h) K∈T h p∈V K q∈V K T K pq |u n+1 p -u n+1 q ||ϕ n+ 1 2 (x p ) -ϕ n+ 1 2 (x q )| η (u n+1 p ) + η (u n+1 q ) 2 . ( 74
)
2. Multiplying (35) by ϕ n p := t n+1 t n ϕ(x p , s) ds, then adding to (72), we obtain

N -1 n=0 p∈N η(u n+1 p ) -η(u n p ) δt m p ϕ n p + N -1 n=0 p∈N q∈Np Φ n p,q (u n p , u n q )ϕ n p + Q T η ( u h,δt )λ h,δt • ∇ϕ dxdt + Q T ϕ|∇η ( u h,δt )| dxdt -e (1) h,δt -e (2)
h,δt ≤ 0. (75)

Now, observing that, since ϕ(•, t) = 0 for t ≥ T , we have ϕ N p = 0 and therefore

N -1 n=0 p∈N η(u n+1 p ) -η(u n p ) δt m p ϕ n p = - N -1 n=0 p∈N η(u n+1 p )m p ϕ n+1 p -ϕ n p δt - p∈N η(u 0 p )m p 1 δt ϕ 0 p
we can rewrite (75) as

Q T η(u h,δt )∂ t ϕ dxdt + Q T Φ(x, t, u h,δt ) • ∇ϕ dxdt - Q T η ( u h,δt )λ h,δt • ∇ϕ dxdt - Q T ϕ|∇η ( u h,δt )| dxdt + R d η(u 0 h )ϕ(x, 0) dx + e (0) h,δt + e (1)
h,δt + e

h,δt + e

h,δt ≥ 0, (

h,δt := p∈N 1 δt δt 0 Qp η(u 0 p )(ϕ(x p , t) -ϕ(x, 0))dxdt, 76) where e (0)
h,δt :=

N -1 n=0 p∈N η(u n+1 p ) t n+1 t n Qp ϕ(x p , t + δt) -ϕ(x p , t) δt -∂ t ϕ(x, t) dxdt, (78) e (4)
h,δt := -

Q T Φ(x, t, u h,δt ) • ∇ϕ dxdt - N -1 n=0 p∈N q∈Np Φ n p,q (u n p , u n q )ϕ n p . (79) Hence, setting e h,δt = e (0) h,δt + e (1) h,δt + e (2) h,δt + e (3) h,δt + e (4) h,δt , (80)
we obtain (64).

Second step: proof of (65)

Study of e

h,δt : use of the vanishing viscous term We define χ K = 0 if ϕ(x, t) = 0 on K × [0, T ) and 1 otherwise. We observe that, for a given n = 0, . . . , N -1,

R d λ n+1 h • ∇ v n+1 -v n+1 h dx = K∈T h χ K λ n+1 h • K ∇ v n+1 -v n+1 h dx,
and we have

K ∇ v n+1 -v n+1 h dx = σ∈E K σ v n+1 (x) -v n+1 h (x) dγ(x).
We use the fact that, for all f ∈ C 2 (σ) and for all x ∈ σ, we have |f

(x) -f h (x)| ≤ max x∈σ |D 2 f (x)|h 2
, denoting by f h (x) the affine function equal to f at the vertices of the simplex σ. We apply this inequality to the function

f (x) = v n+1 (x) = η ( u n+1 h (x))ϕ n+ 1 2 (x). Using that u n+1 h is affine on σ with tangential gradient bounded by |∇ u n+1 h|K |, we get, for all x ∈ σ, letting C η 3 be a bound of |η |, |η | and |η | and C ϕ 2 be a bound of ϕ n+ 1 2 , |∇ϕ n+ 1 2 | and |D 2 ϕ n+ 1 2 |, | v n+1 (x) -v n+1 h (x)| ≤ h 2 C η 3 C ϕ 2 (|∇ u n+1 h|K | 2 + 2|∇ u n+1 h|K | + 1
). Since the above expression is integrated over σ,we get

R d λ n+1 h • ∇ v n+1 -v n+1 h dx ≤ K∈T h σ∈E K χ K |σ|h 2 C η 3 C ϕ 2 (|∇ u n+1 h|K | 2 + 2|∇ u n+1 h|K | + 1).
Multiplying by δt and summing on n = 0, . . . , N -1, we obtain e

h,δt ≤ N -1 n=0 δt K∈T h σ∈E K χ K |σ|h 2 C η 3 C ϕ 2 (|∇ u n+1 h|K | 2 + 2|∇ u n+1 h|K | + 1). (1)
Owing to the geometrical hypotheses and to (42), we get

N -1 n=0 δt K∈T h σ∈E K χ K |σ|h 2 |∇ u n+1 h|K | 2 ≤ Ch N -1 n=0 δt K∈T h |K||∇ u n+1 h|K | 2 ≤ h θ(h) C. (81)
Besides, we have, thanks to (41),

N -1 n=0 δt K∈T h σ∈E K χ K |σ|h 2 |∇ u n+1 h|K | ≤ Ch N -1 n=0 δt K∈T h |K||∇ u n+1 h|K | ≤ hC. Finally N -1 n=0 δt K∈T h σ∈E K χ K |σ|h 2 ≤ ChT K∈T h χ K |K| ≤ ChT (|supp(ϕ)| + h),
where supp(ϕ) = {x ∈ R d , ∃t ∈ [0, T ), ϕ(x, t) = 0}. Hence each of the above terms tends to 0 with h thanks to Hypotheses (HD), thus completing the proof that lim

h,δt→0 e (1)
h,δt = 0.

Study of e

h,δt : proof that the viscous term is vanishing

Letting C ϕ be a Lipschitz constant for ϕ n+ 1 2 , and C a bound for η (u n+1 p ) (this term remains bounded since η is continuous and a 0 ≤ u n+1 p ≤ b 0 , see Proposition 3.8), we get, thanks to the Cauchy-Schwarz inequality

N -1 n=0 δt K∈T h p∈V K q∈V K T K pq |u n+1 p -u n+1 q | |ϕ n+ 1 2 (x p ) -ϕ n+ 1 2 (x q )| η (u n+1 p ) + η (u n+1 q ) 2 ≤ CC ϕ   N -1 n=0 δt K∈T h p∈V K q∈V K T K pq (u n+1 p -u n+1 q ) 2   1/2   N -1 n=0 δt K∈T h p∈V K q∈V K T K pq χ K |x p -x q | 2   1/2 .
Thanks to the geometrical hypotheses and to (42), we get, multiplying by θ(h), e

h,δt ≤ CC ϕ θ(h) C θ(h) C T (|supp(ϕ)| + h) 1/2 (2)
.

We get lim h,δt→0 e

h,δt = 0 thanks to the hypotheses (HD).

Study of e (0) h,δt

We observe that, for all x ∈ Q p and t ∈ [0, δt],

|ϕ(x p , t) -ϕ(x, 0)| ≤ h max (x,
h,δt = 0.

Study of e (4) h,δt

Since we are using below the BV estimate (41), which only involves values n ≥ 1, we define

I 0 := p∈N q∈Np Φ 0 p,q (u 0 p , u 0 q )ϕ 0 p , I := N -1 n=1 p∈N q∈Np
Φ n p,q (u n p , u n q )ϕ n p , and

I * 0 := - δt 0 R d Φ(x, t, u h,δt ) • ∇ϕ dxdt, I * := - T δt R d Φ(x, t, u h,δt ) • ∇ϕ dxdt.
We then have e

h,δt = -I 0 -I + I * 0 + I * . We write, as in the proof of Proposition 3.11, again using the property Φ n p,q (u, v) = -Φ n q,p (v, u) for all (u, v) ∈ [a 0 , b 0 ] 2 ,

I 0 = p∈N q∈Np (Φ 0 p,q (u 0 p , u 0 q ) -Φ 0 p,q (0, 0))ϕ 0 p = p∈N q∈Np (Φ 0 p,q (u 0 p , u 0 q ) -Φ 0 p,q (0, 0))(ϕ 0 p - 1 2 (ϕ 0 p + ϕ 0 q )).
We observe that Proposition 3.10 implies

|Φ 0 p,q (u 0 p , u 0 q ) -Φ 0 p,q (0, 0)| ≤ L m pq (|u 0 p | + |u 0 q |) ≤ 2(b 0 -a 0 )L m pq ,
and that

|ϕ 0 p -ϕ 0 q | ≤ δt h max (x,t)∈Q T |∇ϕ(x, t)|χ pq ,
where we denote by χ pq = 0 if ϕ(x, t) = 0 on (Q p ∪ Q q ) × [0, T ) and 1 otherwise. This leads to

|I 0 | ≤ (b 0 -a 0 )L δt max (x,t)∈Q T |∇ϕ(x, t)| p∈N q∈Np
χ pq m pq h, which leads, thanks to geometrical hypotheses (HD), to

|I 0 | ≤ Cδt.
We find as well that |I * 0 | ≤ Cδt. Let us now turn to the study of I -I * . From (3.5), it follows q∈Np Φ n p,q (u, u) = 0 and thus

I = N -1 n=1 p∈N q∈Np Φ n p,q (u n p , u n q ) -Φ n p,q (u n p , u n p ) ϕ n p = I 1 -I 2 ,
where

I 1 := N -1 n=1 (p,q)∈E h Φ n p,q (u n p , u n q ) -Φ n p,q (u n p , u n p ) ϕ n p , I 2 := N -1 n=1 (p,q)∈E h Φ n p,q (u n p , u n q ) -Φ n p,q (u n q , u n q ) ϕ n q .
Applying the divergence theorem on each cell in the expression of I * , we find

I * = - N -1 n=1 p∈N t n+1 t n ∂Qp Φ(x, t, u n p ) • ν p,q ϕ(x, t) dγ(x)dt = - N -1 n=1 (p,q)∈E h t n+1
t n σp,q Φ(x, t, u n p ) • ν p,q -Φ(x, t, u n q ) • ν p,q ϕ(x, t) dγ(x)dt

= I * 1 -I * 2 ,
where

I * 1 := N -1 n=1 (p,q)∈E h t n+1
t n σp,q 1 m p,q Φ n p,q (u n p , u n q ) -Φ(x, t, u n p ) • ν p,q ϕ(x, t) dγ(x)dt,

I * 2 := N -1 n=1 (p,q)∈E h t n+1
t n σp,q 1 m p,q Φ n p,q (u n p , u n q ) -Φ(x, t, u n q ) • ν p,q ϕ(x, t) dγ(x)dt.

We first rewrite I * 1 as

I * 1 = N -1 n=1 (p,q)∈E h Φ n p,q (u n p , u n q ) -Φ n p,q (u n p , u n p ) 1 m p,q t n+1 t n σp,q ϕ(x, t) dγ(x)dt + N -1 n=1 (p,q)∈E h t n+1 t n σp,q 1 m p,q Φ n p,q (u n p , u n p ) -Φ(x, t, u n p ) • ν p,q ϕ(x, t) dγ(x)dt.
Then, the consistency of the family {Φ n p,q } with Φ allows us to turn the above equation into

I * 1 = N -1 n=1 (p,q)∈E h Φ n p,q (u n p , u n q ) -Φ n p,q (u n p , u n p ) 1 m p,q t n+1 t n σp,q ϕ(x, t) dγ(x)dt + N -1 n=1 (p,q)∈E h t n+1 t n σp,q Φ(x, t, u n p ) • ν p,q 1 δtm p,q t n+1 t n σp,q
(ϕ(y, s) -ϕ(x, t)) dγ(y)dsdγ(x)dt.

Hence we can now write I 1 -

I * 1 = A 1 + B 1 with A 1 := N -1 n=1 (p,q)∈E h Φ n p,q (u n p , u n q ) -Φ n p,q (u n p , u n p ) 1 m p,q t n+1
t n σp,q (ϕ(x p , t) -ϕ(x, t)) dγ(x)dt,

B 1 := N -1 n=1 (p,q)∈E h t n+1
t n σp,q Φ(x, t, u n p ) • ν p,q 1 δtm p,q t n+1 t n σp,q (ϕ(x, t) -ϕ(y, s)) dγ(y)dsdγ(x)dt, and similarly

I 2 -I * 2 = A 2 + B 2 with A 2 := N -1 n=1 (p,q)∈E h Φ n p,q (u n p , u n q ) -Φ n p,q (u n q , u n q ) 1 m p,q t n+1
t n σp,q (ϕ(x q , t) -ϕ(x, t)) dγ(x)dt,

B 2 := N -1 n=1 (p,q)∈E h t n+1 t n σp,q Φ(x, t, u n q ) • ν p,q 1 δtm p,q t n+1 t n
σp,q (ϕ(x, t) -ϕ(y, s)) dγ(y)dsdγ(x)dt.

From the mean value theorem and the uniform Lipschitz continuity of {Φ n p,q }, we thus derive the estimates

|A 1 | ≤ D, |A 2 | ≤ D and |B 1 -B 2 | ≤ D, with D = C(h + δt) N -1 n=1 (p,q)∈E h δtm p,q |u n p -u n q |.
Using the finite number of neighbours in D h , hypothesis (HD3), and estimate (41), we obtain the bound

N -1 n=1 (p,q)∈E h δtm p,q |u n p -u n q | ≤ C N -1 n=1 K∈T h δth d |∇ u n h|K | ≤ C N -1 n=1 δt R d |∇ u n h | dx ≤ C. (82)
Finally we deduce the estimate |I -I * | ≤ C(h + δt), which yields |e

h,δt | ≤ C(h + δt).

Convergence of the approximate solutions

The following lemma is used in the course of the proof of the convergence theorem.

Lemma 6.1: For all n ∈ N, let u n ∈ L 1 (0, T ; BV (R d )) be such that:

1. there exists C ≥ 0 such that

T 0 T V (u n (•, t))dt ≤ C for all n ∈ N, 2. there exists u ∈ L 1 loc (Q T ) such that u n → u in L 1 loc (Q T ). Then u ∈ L 1 (0, T ; BV (R d )) and ∀ϕ ∈ C c (R d × [0, T ), R + ), lim inf n→∞ T 0 R d ϕ(x, t)|Du n (•, t)|(dx)dt ≥ T 0 R d ϕ(x, t)|Du(•, t)|(dx)dt. (83)
Proof Since u n → u in L 1 loc (Q T ), up to a subsequence, for a.e. t ∈ (0, T ), we have that

u n (•, t) → u(•, t) in L 1 loc (R d ). Therefore, for a.e. t ∈ (0, T ), T V (u(•, t)) ≤ lim inf n→∞ T V (u n (•, t)),
where the quantities in the above inequality may be equal to +∞. Integrating the above inequality and applying Fatou's lemma, we get

T 0 T V (u(•, t))dt ≤ T 0 lim inf n→∞ T V (u n (•, t))dt ≤ lim inf n→∞ T 0 T V (u n (•, t))dt ≤ C. Therefore u ∈ L 1 (0, T ; BV (R d ))
. Including the weight ϕ in the definition of the total variation, the proof works the same for getting (83).

We may now state the convergence theorem.

Theorem 6.2: Let us assume Hypotheses (HC). Let (F h k ,δt k ) k∈N be a sequence of discretisations which uniformly satisfies the conditions (HD1)-(HD3) and whose meshsize and time step tend to zero. For any k ∈ N, let (u h k ,δt k ) k∈N , ( u h k ,δt k ) k∈N and (λ h k ,δt k ) k∈N be a solution to Scheme ( 26)- [START_REF] Ziemer | Weakly differentiable functions[END_REF]. Then there exists u ∈ L ∞ (Q T ) ∩ L 1 (0, T ; BV (R d )), and λ ∈ L ∞ (Q T ) d , with |λ| ≤ 1 almost everywhere on Q T , such that, up to a subsequence,

u h k ,δt k → u in L 1 loc (Q T ), u h k ,δt k → u in L 1 loc (Q T ), λ h k ,δt k λ weakly- * in L ∞ (Q T ). (84)
Moreover, u is the unique entropy solution of (2)-(3), and the whole sequences (u h k ,δt k ) k∈N , ( u h k ,δt k ) k∈N converge to u in the above sense.

Proof 1. The estimates (46) and (58) allow us to apply Kolmogorov theorem to the sequence (u h k ,δt k ) k∈N . Thus, there exist u ∈ L 1 loc (Q T ) and a subsequence of (F h k ,δt k ) k∈N , again denoted by (F h k ,δt k ) k∈N , such that (u h k ,δt k ) k∈N converges to u in L 1 loc (Q T ). Applying ( 21) and (41), we get that

u h k ,δt k → u in L 1 loc (Q T ) as well. Applying Lemma 6.1, we get that u ∈ L ∞ (Q T ) ∩ L 1 (0, T ; BV (R d )), since the sequence ( u h k ,δt k ) k∈N is uniformly bounded in L ∞ (Q T ) and L 1 (0, T ; BV (R d )). 2. As the sequence (λ h k ,δt k ) k∈N is bounded in L ∞ (Q T ), there exists λ ∈ L ∞ (Q T ) d such that, up to a subse- quence, (λ h k ,δt k ) k∈N converges weakly- * to λ in L ∞ (Q T ). Furthermore, |λ| ≤ 1 almost everywhere on Q T , since |λ h k ,δt k | < 1 on Q T . 3.
Let us now consider an entropy η and a test function ϕ. Since the sequence (u h k ,δt k ) k∈N is bounded and converges to u in L 1 loc (Q T ), we have

Q T η(u h k ,δt k )∂ t ϕ dxdt + Q T Φ(x, t, u h k ,δt k ) • ∇ϕ dxdt → Q T η(u)∂ t ϕ + Φ(x, t, u) • ∇ϕ dxdt. Since u h k ,δt k → u in L 1 loc (Q T ) and λ h k ,δt k λ weakly- * in L ∞ (Q T ), Q T η ( u h k ,δt k )λ h k ,δt k • ∇ϕ dtdx → Q T η (u)λ • ∇ϕ dtdx.
Again applying Lemma 6.1 to the family η ( u h k ,δt k ), we obtain lim inf

k→+∞ Q T ϕ|∇η ( u h k ,δt k )| dtdx ≥ Q T ϕ D[η (u)] dt. (85)
From ( 26), we get that u 0 h k → u ini in L 1 loc (R d ), and thus

R d η(u 0 h )ϕ(x, 0) dx → R d η(u ini )ϕ(x, 0) dx.
Finally, using the above limits and Proposition 5.1, we can pass to the limit in (64) and find

Q T η(u)∂ t ϕ + Φ(x, t, u) -λη (u) • ∇ϕ dxdt - Q T ϕ D[η (u)] dt + R d η(u ini )ϕ(x, 0) dx ≥ 0,
which proves that u is the entropy solution. Owing to the uniqueness of the entropy solution (proved in Section 1.4), we conclude that, in fact, the whole sequence ( u h k ,δt k ) k∈N converges to u.

Remark 6.3:

In the present article, we have proposed a finite volume-finite element approximation of a simplified scalar model of inviscid Bingham flow and proved its convergence. The convergence analysis relies on L ∞ (Q T ) and L 1 (0, T ; BV (R d )) estimates of the discrete solution. The L ∞ (Q T ) estimate is obtained by a maximum principle in both the hyperbolic flow step and the total variation flow step. The extension of the finite volume-finite element approximation to the vector-valued inviscid Bingham flow is quite straightforward. We refer to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for the finite volume approximation of hyperbolic systems and [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF][START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF] for the approximation of the vectorial total variation flow. However, the convergence analysis in this case seems out of reach. The difficulties do not come from the vectorial total variation flow: the BV estimate is still valid and, although there is no longer a maximum principle, a L ∞ estimate can still be proven [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF]. On the contrary, most of the tools used to prove the convergence of a finite volume approximation of a nonlinear hyperbolic scalar equation (maximum principle, L ∞ estimate, weak BV estimates) are no longer valid in the vectorial case.

Numerical examples

To illustrate the behaviour of the numerical scheme, simulations are performed on two examples: a onedimensional problem where the flux function F describes a Burgers flow (we use Scilab for these simulations) and a two-dimensional problem where the flux function F describes a mere advection (we use FreeFem++ for these simulations). We numerically solve the two-dimensional problem posed on R 2 by considering periodic boundary conditions. Since in this case, an analytical solution is known, a quantitative study of the convergence is provided. In both cases, we do not implement the additional term in θ(h), which allows for keeping a compact support for the approximation of u in the numerical tests.

Burgers flow and total variation flow in 1D

We consider the equation

∂ t u + 1 2 ∂ x (u 2 ) -g ∂ x Sgn(∂ x u) = 0, (x, t) ∈ R × (0, T ), (86)
with g = 5 • 10 -4 . The initial datum is a bump (see Figure 1, left). The computational domain is taken large enough for the solution not to touch the boundary in the course of the simulation. The finite element and finite volume meshes, of size h, are staggered. The finite volume step is advanced using the Godunov scheme [24, see, e.g.,], which is stable under the CFL condition sup p |u n p |δt ≤ 1 2 h, for all n ∈ {0, ..., N -1}. As in the previous test case, the finite element step is solved using the lagged diffusivity algorithm The numerical simulation is performed with h = 0.002, δt = 0.001, = 10 -8 and 20 iterations at each time step for the lagged diffusivity algorithm. In Figure 1, the computed solution is represented at different times. As in a Burgers flow (without total variation flow), a shock wave forms at the right end, while a rarefaction wave forms at the left end. The main difference introduced by the total variation flow is the progressive clipping of the solution.

Advection and total variation flow in 2D

We consider the problem ∂ t u + c x ∂ x u + c y ∂ y u -g div Sgn(∇u) = 0, (x, y, t) ∈ R 2 × (0, T ), (87) where c x , c y , and g are constants. The initial condition is u(x, y, 0) = 1 D0 (x, y),

(x, y) ∈ R 2 , ( 88
)
where 1 D0 denotes the indicator function of D 0 , the disk of radius r 0 centred at point (x 0 , y 0 ). It is possible to determine the exact solution to this problem. Indeed, for some bounded sets S of finite perimeter in R 2 , including disks, it is proved in [START_REF] Bellettini | The total variation flow in R N[END_REF] that the solution of ∂ t u -div Sgn(∇u) = 0, (x, y, t) ∈ R 2 × (0, T ), u(x, y, 0) = 1 S (x, y), (x, y) ∈ R 2 , is the function u : (x, y, t) → (1 -|∂S|/|S|t) + 1 S (x, y). The proof can be easily adapted to deal with the advection and the exact solution of (87)-( 88) is found to be u : (x, y, t) → 1 -2gt r 0 + 1 D0 (x -c x t, y -c y t).

The numerical simulations use the following parameters: (x 0 , y 0 ) = (0.3, 0.3), r 0 = 0.2, c x = 0.8, c y = 0.2, g = 0.03. To mimick the space R 2 , square domains with periodic boundary conditions are used for both the finite volume step and the finite element step. The finite element domain, denoted by Ω, is meshed with squares of size h divided into two triangles. The finite volume domain, slightly larger than the finite element domain, is meshed with square cells of size h, centred at each node of the finite element mesh. Such a pair of meshes is represented in Figure 2. The advection step is solved using the so-called corner-transport upwind scheme [24, Fig. 2: Example of meshes used for the two-dimensional simulations. The finite element mesh is in solid line and the finite volume mesh is in dashed line.

see, e.g.,], which is slightly more accurate than the standard upwind scheme. This scheme is stable under the CFL condition δt ≤ min h cx , h cy . The finite element step is solved using the lagged diffusivity algorithm, whose principle is recalled in Appendix A.

Remark 7.1: The corner-transport upwind scheme does not fit exactly into the theoretical framework defined in Section 3.1 since there is a numerical flux between the cells that touch at corner. However, the well-posedness and convergence results remain valid (the proofs can be readily adapted). The first simulation is performed with Ω = (0, 1) × (0, 1), h = 0.01, CFL= 1, = 10 -6 and 20 iterations at each time step for the lagged diffusivity algorithm. The numerical solution, represented at different times in Figure 3, is in good agreement with the exact solution: the contours are sharp and the support is little deformed. Interestingly, the total variation minimisation in the finite element step limits the diffusion due to the finite volume scheme; see Figure 4. Table 1 collects the relative error between the exact solution and the discrete solution in L 1 -norm for various meshsizes (keeping a CFL equal to 1). Precisely, the relative error is given by the expression u per -u h,δt L 1 (Ω×(0,T )) u per L 1 (Ω×(0,T ))

where u per is the periodic version of u and the norms are computed with a second-order accurate quadrature formula. The results show, as expected, a decrease of the error when the mesh becomes finer. The computed convergence rate is clearly sublinear (0.28 between the first mesh and the second mesh, 0.21 between the second mesh and the third mesh). It is probably due to the discontinuity in the solution. Indeed, the optimal L 1 -norm convergence rate for a monotone scheme applied to a linear advection problem with discontinuous data is only

O(h 1
2 ) [START_REF] Tang | The sharpness of Kuznetsov's O( √ ∆x) L 1 -error estimate for monotone difference schemes[END_REF]. The convergence rate for the total variation minimisation with finite elements seems also to be sublinear in the case of a discontinuous source term (the best L 2 -norm convergence rate proved in the literature is O(h Remark 7.2: The regularisation parameter is taken very small and the number of iterations is relatively large in order to make negligible the error due to the lagged diffusivity algorithm. It is actually possible to choose a larger regularisation parameter (and thus a smaller number of iterations) without deteriorating the convergence rate of the finite element step. We refer to [START_REF] Feng | Rate of convergence of regularization procedures and finite element approximations for the total variation flow[END_REF] for the scaling law relating the meshsize to the regularisation parameter.

Definition 2 . 1 :

 21 where a ∨ b denotes the maximum of a and b and a ∧ b denotes the minimum of a and b. The entropy formulation of the problem (2)-(3), owing to the term div Sgn(∇u), requires more regular entropies. Under Hypotheses (HC), an admissible entropy is a convex Lipschitz continuous function η ∈ C ∞ (R). The corresponding entropy flux is the locally Lipschitz continuous function Φ

Remark 2 . 2 (

 22 Consistency with Kruzhkov entropy pairs): Definition 2.1 is such that, for η

h)

  |K | and otherwise that |λ n+1 K | ≤ 1. Using that for any a, b ∈ R d with a = 0, we have |a + b| -|a| -a•b |a| = |a||a+b|-a•(a+b) |a| ≥ 0, and |b| -λ n+1 K • b ≥ 0, we get A ≥ 0, which proves that J n+1 h (w h ) -J n+1 h (u n+1 h ) ≥ 0, and therefore that J n+1 h reaches its minimum value at u n+1 h . Reciprocally, let us assume that J n+1 h reaches its minimum value at some point u n+1 h ∈ X h . Let us prove that there exists λ n+1 h ∈ Λ h such that (29)-(30) holds. We denote by T n+1 h,0 = {K ∈ T h , (∇ u n+1 h ) |K = 0}, and we define the linear form L n+1 h

n+ 1 2 p

 2 ≤ b 0 for all p ∈ N and letting η(s) tend to s ∨ b 0 -b 0 (this is possible since u ini ∈ L 1 (R d ) implies a 0 ≤ 0 ≤ b 0 ), we get that η(u n+1 p ) = 0, which shows that u n+1 p ≤ b 0 . The same reasoning holds, letting η(s) tend to a 0 -s ∧ a 0 , which shows that u n+1 p ≥ a 0 . Finally, letting η(s) tend to |s|, we conclude that u n+1 h ∈ L 1 (R d ).

Remark 4 . 2 : 4 . 3 :

 4243 Throughout this section and the next one, C denotes a generic constant independent of the discretisation F h,δt .4.1 L 1 (0, T ; BV (R d )) estimate Proposition Let us assume Hypotheses (HC) and Hypotheses (HD). Let

Lemma 4 . 5 :

 45 Let T > 0, δt ∈ (0, T ) and (a n ) n∈N be a family of non negative real values. Then, defining the function ν : R → Z such that ν(t) = n + 1 if t ∈ (nδt, (n + 1)δt]

Proposition 4 . 6 :

 46 Let us assume Hypotheses (HC) and Hypotheses (HD). Let u

)

  Proof For a given element K ∈ T h and a given couple of points (a, b) ∈ R d × R d , we denote by χ K (a, b) the length of the segment [a, b]∩K. In particular, if [a, b] does not intersect K, then χ K (a, b) = 0. Let n ∈ {1, ..., N } and (x, y) ∈ R d ×R d . By applying the mean value theorem on each element intersected by the segment [x, x+y],

Fig. 1 :

 1 Fig. 1: Burgers flow and total variation flow. Numerical solution at times t = 0, t = 0.05, and t = 0.15.

Tab. 1 :

 1 Relative error between the discrete and exact solutions in L 1 -norm.

Fig. 3 :

 3 Fig.3: Advection and total variation flow. Numerical solution at times t = 0, t = 0.5, t = 1, t = 1.5. The mesh is warped according to the value of the solution and the color scale goes from blue (0) to red (1).

Fig. 4 :

 4 Fig. 4: Numerical solution of the advection and total variation flow (left) and numerical solution of the advection without total variation flow (right) at time t = 1.5. The color scale goes from blue (0) to red (1). The black circle indicates the support of the exact solution.

  with |λ| ≤ 1 almost everywhere on Q T , such that, for all admissible entropy-entropy flux pairs (η, Φ) in the sense of Definition 2.1 and all nonnegative test functions

A Lagged diffusivity algorithm

The lagged diffusivity algorithm is a simple and efficient algorithm for solving the discrete total variation minimisation problem [START_REF] Vogel | Iterative methods for total variation denoising[END_REF][START_REF] Dobson | Convergence of an iterative method for total variation denoising[END_REF]. Under the notation of Section 3.1, consider the minimisation over Y h of the functional J h : Y h → R, defined by

, δt > 0 are given and where

The lagged diffusivity algorithm replaces the nondifferentiable functional J h by the regularised functional J h :

) is a small positive function on R d . Then the respective minimisers u h and u h are the respective solutions of the following variational inequalities:

Letting v h = u h in the first one and v h = u h in the second one, adding the resulting equations, using

We now remark that the minimiser u h of J h is also solution of the variational problem

The preceding problem is then computed with a fixed-point algorithm. The algorithm is initialised with u

) dx = 0, ∀v h ∈ Y h .

The diffusivity field in the above diffusion problem depends on u (k) h , whence the name lagged diffusivity.