Fully discrete finite element approximation of unsteady flows of implicitly constituted incompressible fluids

Implicit constitutive theory provides a very general framework for fluid flow models, including both Newtonian and generalized Newtonian fluids, where the Cauchy stress tensor and the rate of strain tensor are assumed to be related by an implicit relation associated with a maximal monotone graph. For incompressible unsteady flows of such fluids, subject to a homogeneous Dirichlet boundary condition on a Lipschitz polytopal domain $\Omega \subset \mathbb{R}^d$, $d \in \{2,3\}$, we investigate a fully-discrete approximation scheme, using a spatial mixed finite element approximation combined with backward Euler time-stepping. We show convergence of a subsequence of approximate solutions, when the velocity field belongs to the space of solenoidal functions contained in $L^\infty(0,T;L^2(\Omega)^d)\cap L^q(0,T;W^{1,q}_0(\Omega)^d)$, provided that $q\in \big(\frac{2d}{d+2},\infty\big)$, which is the maximal range for $q$ with respect to existence of weak solutions. This is achieved by a technique based on splitting and regularizing, the use of a solenoidal parabolic Lipschitz truncation method, a local Minty-type monotonicity result, and various weak compactness results.


Introduction
In the mechanics of viscous incompressible fluids typical constitutive relations relate the shear stress tensor to the rate of strain tensor through an explicit functional relationship. In the case of a Newtonian fluid the relationship is linear, and in the case of generalized Newtonian fluids it is usually a powerlaw-like nonlinear, but still explicit, functional relation. Implicit constitutive theory was introduced in order to describe a wide range of non-Newtonian rheology, by admitting implicit and discontinuous constitutive laws, see [33,34]. The existence of weak solutions to mathematical models of this kind was explored, for example, in [10,11] for steady and unsteady flows, respectively. The aim of the present paper is to construct a fully discrete numerical approximation scheme, in the unsteady case, for a class of such implicitly constituted models, where the shear stress and rate of strain tensors are related through a (possibly discontinuous) maximal monotone graph. The scheme is based on a spatial mixed finite element approximation and a backward Euler discretization with respect to the temporal variable. We will show weak convergence (up to subsequences) to a weak solution of the model. The mathematical ideas contained in the paper are motivated by the existence theory formulated, in the unsteady case, in [11], and the convergence theory for finite element approximations of steady implicitly constituted fluid flow models developed in [16].

Implicit Constitutive Law
Statement of the Problem. Let Ω ⊂ R d with d ≥ 2 be a bounded Lipschitz domain and denote by Q = (0, T ) × Ω the parabolic cylinder for a given final time T ∈ (0, ∞). Furthermore, let f : Q → R d be a given external force and let u 0 : Ω → R d be an initial velocity field. We seek a velocity field u : Q → R d , a pressure π : Q → R, and a stress tensor field S S S : Q → R d×d sym satisfying the balance law of linear momentum and the incompressibility condition: ∂ t u + div(u ⊗ u) − div S S S = −∇π + f , on (0, T ) × Ω, div u = 0, on (0, T ) × Ω, (1.1) subject to the following initial condition and no-slip boundary condition: in Ω, (1.2) In order to close the system we need to impose a relation, the so-called constitutive law, G G G(·, D D Du, S S S) = 0 0 0, (1.4) between the stress tensor S S S and the symmetric gradient D D Du = 1 2 (∇u + (∇u) ⊤ ), which represents the shear rate of the fluid. In the following we will refer to the problem consisting of (1.1)-(1.4) as (P), and will use the notation z := (t, x) ∈ Q.
The relation G G G may be fully implicit and we assume that G G G can be identified with a maximal monotone graph A(z) ⊂ R d×d sym × R d×d sym , for z ∈ Q, as G G G(z, D D Du(z), S S S(z)) = 0 0 0 ⇔ (D D Du(z), S S S(z)) ∈ A(z), where A(·) satisfies the following assumption, similarly as in [10, p. 110] and [11,Sec. 1.2].
Assumption 1.1 (Properties of A(·)). We assume that A(·) ⊂ R d×d sym × R d×d sym satisfies the following conditions for a.e. z ∈ Q: (i) In [11] the authors phrase the condition (A4) in the more general context of Orlicz-Sobolev spaces. Here we will restrict ourselves to the usual Sobolev setting. This framework covers explicit relations, including Newtonian fluids, where q = 2, and q-fluids describing shear-thinning and shear-thickening behaviour, for 1 < q < 2 and q > 2, respectively. Also relations, where the the stress is a set-valued or discontinuous function of the symmetric gradient, as for Bingham and Herschel-Bulkley fluids, are included, which is shown in [11,Lem. 1.1]. Furthermore, fully implicitly constituted fluids are covered and the constitutive relation is allowed to depend on (t, x) ∈ Q.

Overview of the Context
For explicit constitutive laws the list of existence results is long, see amongst others [31,32,21,18] and the references therein. For those constitutive laws the main challenge is the lack of admissibility for small q, caused by the presence of the convective term. In the most recent results this difficulty was overcome by the use of the Lipschitz truncation method.
In the case of implicitly constituted fluids, the existence of weak solutions for q > 2d d+2 for steady and unsteady flows was proved in [10,11]. The restriction on q is required to ensure compactness of the embedding W 1,q (Ω) ֒→ L 2 (Ω), and implies that the convective term in the weak formulation is well-defined. In [11], a Navier slip boundary condition and C 1,1 regularity of ∂Ω were assumed to avoid technicalities related to lack of regularity of the pressure in the unsteady case. Due to the weak structural assumptions the existence of short-time strong solutions and uniqueness cannot be expected to hold in general. The proof in [11] is constructive and is based on a three-level approximation using finite-dimensional Galerkin subspaces spanned by eigenfunctions of higher order elliptic operators. These Galerkin spaces are not available for practical computations and therefore we take an alternative route here in the construction of a numerical method for the problem and for its convergence analysis.
Here we shall consider a mixed finite element approximation under minimal regularity hypotheses, so we can only hope for qualitative convergence results rather than qualitative error bounds in terms of the spatial and temporal discretization parameters. The approximation scheme will be constructed for a regularized version of the equations, and after passing to the limit with the discretization parameters we shall pass to the limit with the regularization parameter.
Concerning the numerical analysis of implicitly constituted fluid flow models, to the best of our knowledge the only results available are those contained in [16,30], which deal with the steady case. By means of a discrete Lipschitz truncation method and various weak compactness results the authors of [16] prove the convergence of a large class of mixed finite element methods for q > 2d d+1 for discretely divergence-free finite element functions for the velocity, and for q > 2d d+2 for exactly (i.e., pointwise) divergence-free finite element functions for the velocity field. In the case of discretely divergence-free mixed finite element approximations the more demanding requirement, that q > 2d d+1 , arises from the (numerical) modification of the trilinear form associated with the convective term in the weak formulation of the problem, whose purpose is to reinstate the skew-symmetry of the trilinear form, lost in the course of the spatial approximation. In [30] an a posteriori analysis is performed for implicitly constituted fluid flow models using discretely divergence-free finite element functions. In the unsteady case no convergence result is available for implicitly constituted fluid flow models, and even those contributions that are focussed on explicit constitutive laws, such as [12], for example, are restricted to the case when q > 2(d+1) d+2 .

Aim and Main Result
Our objective here is to establish a convergence result for implicitly constituted fluids in the unsteady case for the whole range q > 2d d+2 . More specifically, the aim is to show weak convergence (up to subsequences) of the sequence of approximate solutions to a weak solution of the problem, the main challenges being to deal with the implicit, possibly discontinuous, relation between stress and shear rate, and with small exponents q arising in the coercivity and boundedness assumption on the implicit relation. To this end, we shall consider a three-level approximation, consisting of approximating the potentially multivalued graph A with a single-valued implicit functional relationship between the shear stress tensor and the rate of strain tensor, performing a regularization of the resulting model, which will enable us to cover the entire range of q > 2d d+2 , and then constructing a fully discrete approximation of the regularized model.
The main contribution of the paper is therefore the following. Let Ω be a Lipschitz polytopal domain, q > 2d d+2 , and assume that we have a pair of inf-sup stable finite element spaces. Also we assume that a suitable approximation of the graph A is available, which will be constructed below. Then, a sequence of approximate solutions to the fully discrete problem exists and the corresponding sequence of finite element approximations converges weakly, up to subsequences, to a weak solution of problem (P), when first taking the graph approximation limit, then the spatial and temporal discretization limits, and finally the regularization limit. The precise formulation of this result is contained in Theorem 4.1 and the notion of weak solution is given in Definition 2.1. The main tools in the proof are a Minty type local monotonicity result proved in [11], and the solenoidal parabolic Lipschitz truncation constructed in [8] to overcome the admissibility problem for small q.
The paper is structured as follows. Section 2 provides the analytical setting, including the graph approximation and the Lipschitz truncation. Section 3 describes the finite element approximation, the approximation of the convective term and the time stepping. In Section 4 we first introduce the approximation levels in detail before giving the convergence proof.

Analytical Preliminaries
By R d×d sym we denote the set of all real-valued symmetric d×d-matrices and we use : for the Frobenius scalar product in R d×d . For ω ⊂ R d we denote by |ω| the d-dimensional Lebesgue measure of ω. By 1 ω we denote the characteristic function of the set ω. For the partial derivatives with respect to time we use the shorthand notation ∂ t f := ∂f ∂t . For ω ⊂ R d open and p ∈ [1, ∞) let (L p (ω), ||·|| L p (ω) ) be the standard Lebesgue space of pintegrable functions and the space of essentially bounded functions when p = ∞. For s ∈ N let (W s,p (ω), ||·|| W s,p (ω) ) be the respective Sobolev spaces. For spaces of vector-valued and tensor-valued functions we use superscripts d and d × d, respectively (except for in norms). By L p 0 (ω) we denote the set of functions in L p (ω) with zero mean integral.
For a general Banach space (X, ||·|| X ), the dual space consisting of all continuous linear functionals on X is denoted by X ′ and the dual pairing is denoted by f, g X ′ ,X , if f ∈ X ′ and g ∈ X. If X is a space of functions defined on ω then we denote the dual pairing by f, g ω := f, g X ′ ,X , in case the space X is known from the context. We also use this notation for the integral of the scalar product f · g of two functions f and g, provided that f · g ∈ L 1 (ω). Furthermore, if ω ⊂ R d is measurable and 0 < |ω| < ∞, then we denote ffl ω f (x) dx := 1 |ω|´ω f (x) dx. Now let Ω ⊂ R d be a bounded domain, let T ∈ (0, ∞) and Q = (0, T ) × Ω. Denote by C ∞ 0 (Ω) the set of all smooth and compactly supported functions on Ω and by C ∞ 0,div (Ω) d the set of all functions in . For a given p ∈ (1, ∞) we let the Hölder exponent p ′ be defined by 1 p + 1 p ′ = 1. Then, if p ∈ (1, ∞), L p ′ (Ω) is the dual space of L p (Ω) and W −1,p ′ (Ω) will denote the dual space of W 1,p 0 (Ω). Further, we define the spaces of divergence-free functions Let C(Ω) be the set of all continuous real-valued functions on Ω and C 0,β (Ω) the set of all Hölder continuous functions on Ω with exponent β ∈ (0, 1]. With C([0, T ]; X) we denote the set of all functions defined on [0, T ] taking values in a Banach space X, which are continuous (with respect to the strong topology in X). Furthermore, we define the space of weakly continuous functions with values in X by We denote by L p (0, T ; X) the standard Bochner space of p-integrable X-valued functions. We use the notation ess lim t→0+ f (t) to indicate that there exists a zero set N (f ) ⊂ [0, T ] such that t ∈ (0, T )\N (f ), when considering the limit of f (t), as t → 0 + .
In the following, c > 0 will denote a generic constant, which can change from line to line and depends only on the given data unless specified otherwise.
for a.e. t ∈ (0, T ), for a.e. z ∈ Q, We choose a pressure-free notion of weak solution because in the unsteady problem subject to homogeneous Dirichlet boundary conditions on Lipschitz domains one can only expect to establish a distributional (in time) pressure, see [39, Ch. III, § 1, pp. 266] and also [38].

Implicit Constitutive Laws
Approximation of A. The implicit relation encoded by A can be viewed as a set-valued map. In order to perform the analysis we require a single-valued map and thus a measurable selection S S S ⋆ of the graph A is chosen, which may have discontinuities.  and S S S ⋆ is (L(Q) ⊗ B(R d×d sym )) − B(R d×d sym )-measurable. Furthermore, for a.e. z ∈ Q, one has that (a4) Let U be a dense set in R d×d sym and let (D D D, S S S) ∈ R d×d sym × R d×d sym . The following are equivalent: (a5) S S S ⋆ is locally bounded, i.e., for a given r > 0 there exists a constant c = c(r) such that To show the existence of solutions to any of the approximate problems considered below, continuity of the (approximate) stress tensor is required. Hence, we introduce the following assumptions on a sequence of approximations of the selection S S S ⋆ .
(α2) There exists a constant c * > 0 and a nonnegative function g ∈ L 1 (Q) such that, for all k ∈ N, for any A A A ∈ R d×d sym and for a.e. z ∈ Q, one has that In the existence proofs in [10,11] and also in [16] the approximating sequence S S S k is chosen as the convolution of the selection S S S ⋆ in the second argument with a mollification kernel.
Example 2.4 (Approximation by Mollification). Let ρ ∈ C ∞ 0 (R d×d sym ) be a mollification kernel, i.e., a nonnegative, radially symmetric function, the support of which is contained in the unit ball B 1 (0 0 0) ⊂ R d×d sym , and which satisfies´Rd×d and define the mollification of S S S ⋆ with respect to the last argument by Lemma 5.1 in the Appendix shows that S S S k satisfies Assumption 2.3.
A possibly more practicable approximation based on a piecewise affine interpolant can be used in the case of a radially symmetric selection function S S S ⋆ under additional regularity assumptions.
Example 2.5 (Approximation by Affine Interpolation). Assume that S ⋆ : Q × R ≥0 → R ≥0 is a measurable function with S ⋆ (z, 0) = 0, for any z ∈ Q, such that S S S ⋆ : Q × R d×d sym → R d×d sym , defined by is a measurable selection of a graph A satisfying Assumption 1.1. Furthermore, we assume that (a1') S ⋆ (z, ·) : R ≥0 → R ≥0 is monotone for a.e. z ∈ Q; (a2') Let I ∈ N and {a i } i∈{1,...,I} ⊂ R > 0, s.t. a 1 < · · · < a I . Let a 0 = 0 and denote by A := i∈{0,...,I} a i the finite set of possible discontinuities. Then, assume that, for a.e. z ∈ Q, we have We construct the approximation as follows. Since the set A is finite, there exists a k 0 ∈ N such that 2 k0 < min i∈{1,...,I} (a i − a i−1 ). Let k ∈ N, with k ≥ k 0 be arbitrary but fixed. Denote for i ∈ {0, . . . , I} Let z ∈ Q be arbitrary but fixed. First we extend S ⋆ (z, ·) as an odd function to R and we still denote the extension by S ⋆ (z, ·). Since the point evaluations S ⋆ (z, a k i,± ), for i ∈ {0, . . . , I}, are well-defined, we can define (2.12) On A k i the approximation S ⋆ (z, ·) is the affine interpolant between S ⋆ (z, a k i,− ) and S ⋆ (z, a k i,+ ) and otherwise S ⋆ is unchanged.
In Lemma 5.2 and Corollary 5.3 in the Appendix it is proved that the resulting approximating sequence S S S k (z, ·) satisfies Assumption 2.3.
Minty's Trick. The following lemma is one of the crucial tools for the identification of the implicit constitutive law upon passage to the limit. Then, we have that (D D D(z), S S S(z)) ∈ A(z) for a.e. z ∈ Q.

Lipschitz Approximation
For small q ∈ (1, ∞) a weak solution according to Definition 2.1 is not an admissible test function because of the presence of the convective term. The Lipschitz truncation method helps to identify the implicit relation despite the lack of admissibility. It first appeared in [1] and since then the method was further developed and refined in a series of papers, see, e.g., [29,21,17,18,7,8,16], to mention just a few. For a sequence of solutions to a sequence of divergence-form evolution equations a solenoidal parabolic Lipschitz truncation was developed in [8]. Note that the sets B l,j in the following Lemma satisfy B l,j = O l,j ∩ Q 0 , where O l,j are the "bad sets" in the construction in [8]. Let p ∈ (1, ∞), σ ∈ (1, min(p, p ′ )) and let Q 0 = I 0 × B 0 ⊂ R × R d be a parabolic cylinder, for d = 3, for an open interval I 0 and an open ball B 0 . Let {v l } l∈N be a sequence of weakly divergencefree functions, which is converging to zero weakly in L p (I 0 ; W 1,p (B 0 ) d ), strongly in L σ (Q 0 ) d and is uniformly bounded in L ∞ (I 0 , L σ (B 0 ) d ). Consider a sequence {G G G l 1 } l∈N , converging to zero weakly in L p ′ (Q 0 ) d×d and a second sequence, {G G G l 2 } l∈N , converging to zero strongly in L σ (Q 0 ) d×d . Furthermore, denoting G G G l := G G G l 1 + G G G l 2 , assume that, for any l ∈ N, the equation is satisfied. Then, there exists a j 0 ∈ N, • a double sequence {λ l,j } l,j∈N ⊂ R with λ l,j ∈ 2 2 j , 2 2 j+1 −1 , for any l, j ∈ N, , and supp(v l,j ) ⊂ 1 6 Q 0 , for any j ≥ j 0 and any l ∈ N; (iv) there exists a constant c > 0 such that (vii) there exists a constant c > 0 such that, for any H H H ∈ L p ′ ( 1 6 Q 0 ) d×d , we have that The lemma is stated for d = 3, but according to [8, Rem. 2.1, p. 2692] the result holds for all d ≥ 2, in fact, with minor modifications of the proof.
Proof. For p ≥ d, the reader is referred to [ Then, by integrating over (0, T ), further estimating the second factor and using the continuous embedding of W 1,p (Ω) ֒→ L p * (Ω), we obtain which completes the proof.
Lemma 2.9 (Continuity I). Let Z be a reflexive Banach space. Then, any function v ∈ L 1 (0, T ;

Finite Element Spaces and Assumptions
The setting here is slightly more general than the one in [16]. (i) Each element K ∈ T n is affine-equivalent to the closed standard reference simplex K := conv{0, e 1 , . . . , e d } ⊂ R d , i.e., there exists an affine invertible function F F F K : K → K; (ii) {T n } n∈N is shape-regular, i.e., there exists a constant c r (independent of n ∈ N) such that hK ρK ≤ c r for all K ∈ T n and all n ∈ N, where h K := diam(K) and ρ K := sup{diam(B) : B is a ball contained in K}; (iii) The grid size h n := max{h K : K ∈ T n } vanishes, as n → ∞.
Finite Element Spaces. Let P V ⊂ W 1,∞ ( K) d and let P Q ⊂ L ∞ ( K) be finite-dimensional function spaces on the reference simplex K (with a slight abuse of notation) as in [16]. Further, let V ⊂ C(Ω) d and let Q ⊂ L ∞ (Ω). Then we define the conforming finite element spaces V n and Q n with respect to the partition T n by Let us also introduce the subspace of discretely divergence-free functions of V n and the subspace of zero integral mean functions of Q n by V n div := {V ∈ V n : div V , Q Ω = 0 for all Q ∈ Q n }, Note that the functions in V n div are in general not divergence-free, so in general V n div ⊂ W 1,∞ 0,div (Ω) d .
Assumption 3.2 (Approximability, [16,Assump. 5]). Assume that for all p ∈ [1, ∞), we have that Projectors. For the convergence analysis we require certain projectors to the respective finite element spaces and suitable assumptions on them. Since we do not need local stability of the projector Π n , we assume less than in [16].
Assume that for each n ∈ N there exists a linear projector Π n : W 1,1 0 (Ω) d → V n such that: (ii) (W 1,p -stability) for any p ∈ (1, ∞) there exists a constant c(p) > 0 (independent of n) such that Assumption 3.4 (Projector Π n Q ). Assume that for each n ∈ N there exists a linear projector Π n Q : L 1 (Ω) → Q n such that, for any p ∈ (1, ∞), there exists a constant c(p) > 0 such that for all h ∈ L p (Ω) and all n ∈ N.
The existence of the Bogovskiȋ operator, see [6] and also [17, p. 223] implies that the continuous inf-sup condition holds for any p ∈ (1, ∞). With this and the Assumption 3.3 the corresponding discrete inf-sup condition holds uniformly in n ∈ N; see Fortin's Lemma for Banach spaces in [20,Lem. 4.19]. This means that the framework results in an inf-sup stable pair (V n , Q n ).
Note that the stability in Assumption 3.3 (ii) is only shown for p = 2, if d = 2.
L 2 -Projector to V n div . Let us introduce the projector onto V n div , given by Directly from the definition we have L 2 -stability and optimality of the approximation in L 2 (Ω) d , i.e., By these properties and an approximation argument using the properties of Π n (see Remark 3.5 (i)) one can show that

Convective Term and its Numerical Approximation
Motivated by the form of the convective term in the conservation of momentum equation, we consider the trilinear form b defined by where the second equality follows by integration by parts. Hence for divergencefree functions u the last term vanishes and b(u, This property can be extended to Sobolev functions. As in general V n div ⊂ W 1,∞ 0,div (Ω) d , the second term in (3.12) need not vanish. To preserve the skewsymmetry of the trilinear form associated with the convective term the usual approach in the numerical analysis literature (see, e.g., [39]) is therefore to consider instead the skew-symmetric trilinear form Thus we have that b(u, v, v) = 0 regardless of the solenoidality of u. Note that we have b(u, ·, ·) = b(u, ·, ·) for divergence-free functions u.
In the equations the terms appear in the form b(u, u, v) and b(u, u, v), for the velocity u and a test function v. The natural function space for weak solutions of the problem is L ∞ (0, T ; On the other hand, for the modification (cf. the first term in (3.13)) of the trilinear form b associated with the convective term one obtains Evidently, the source of this more restrictive requirement on q is the modification of the trilinear form b, introduced in order to reinstate the skew symmetry of b, lost in the course of approximating the pointwise divergence-free solution by discretely divergence-free finite element functions. We note in passing that the restriction q ≥ 2(d+1) d+2 in the unsteady case corresponds to the restriction q ≥ 2d d+1 in the steady case in [16].
To motivate the choice of the form of the additional regularization term that we shall add to the weak form to relax the excessive restriction q ≥ 2(d+1) d+2 to the natural restriction on q > 2d d+2 , we note that by Hölder's inequality we have that , without any restrictions on the range of q, other than q ∈ (1, ∞). This justifies the use of a regularizing term guaranteeing additional L 2q ′ -integrability.

Time Discretization
For the purpose of time discretization, let l ∈ N and define the time step by δ l = T /l → 0, as l → ∞. For l ∈ N, we shall use the equidistant temporal grid on [0, T ] defined by {t l i } i∈{0,...,l} , where t l i := iδ l , for i ∈ {0, . . . , l}. In the following we will suppress the superscript l and write t i , i ∈ {0, . . . , l}.

Convergence Proof
Motivated by the approach in [11,Sec. 3.1] we consider the following levels of approximation.
k ∈ N: The selection S S S ⋆ given in Lemma 2.2 is approximated by a family of Carathéodory functions {S S S k } k∈N , which satisfy Assumption 2.3. The approximation of the stress is then explicit and continuous in D D Du. l ∈ N: A time stepping based on the implicit Euler method is introduced similarly as, e.g., in [12,39], see Subsection 3.3. n ∈ N: The velocity u is approximated by a Galerkin approximation in finite element spaces in the spatial variable, see Section 3. m ∈ N: The regularizing term 1 m |u| 2q ′ −2 u is added to the equation to gain admissibility of the approximate solutions if q ≤ 3d+2 d+2 and to enable us to use the bound on b(·, ·, ·) in (3.16), without imposing the restriction q > 3d+2 d+2 . This results in a fully discrete approximation. The limits are taken in the order k → ∞, l, n → ∞, and then m → ∞, and we can take the limits in l, n → ∞ simultaneously. To simplify the notation we shall write to denote the fact that the limits k, (l, n), m are taken successively in the order of indexing (from leftto right) and the space X describes the weakest topology of the three limits. We will use analogous notation for weak and weak* convergence. In each step one has to identify the equation and the implicit relation, which is the most challenging part. The most significant difference compared to [11] lies in the passage to the limits l, n → ∞ and the identification of the implicit law.
As both the external force f and the approximate stress S S S k will be allowed to be time-dependent, and the time-dependence will not be assumed to be continuous, we shall consider an integral-averaged version in the approximate problem. Let us recall the notation in Subsection 3.3 and introduce for f ∈ L q ′ (0, T ; W −1,q ′ (Ω) d ) and S S S k : Q × R d×d sym → R d×d sym as in Assumption 2.3 and l ∈ N the averages with respect to the time grid {t i } i∈{0,...,l} defined, for i ∈ {1, . . . , l}, by and let the corresponding piecewise constant interpolants f and S S S k be defined as in Subsection 3.3 (3.18). Recall that by (3.24) and (3.25) we have that for k, l, n, m ∈ N, i ∈ {1, . . . , l} and b(·, ·, ·) as defined in (3.13).
Approximate Problem: For k, l, n, m ∈ N find a sequence {U k,l,n,m i } i∈{0,...,l} ⊂ V n div such that and for a given U k,l,n,m i−1 ∈ V n div the approximate solution on the next time level, U k,l,n,m . . , l} a fully implicit problem has to be solved, since the numerical solution from the previous time level only appears in the term involving d t U k,l,n,m i , as defined in (3.17).
as k → ∞, (l, n) → ∞ (combined) and m → ∞, when taking the limits successively, without restrictions on the relation between the discretization parameters δ l and h n .

Remark 4.2.
(i) In the proof of Theorem 4.1 it is essential that the limits are taken in the indicated order.
(ii) If S S S ⋆ is a Carathéodory function, then the approximation level corresponding to k ∈ N can be skipped. (iii) For Lipschitz polytopal domains Theorem 4.1 is also a new existence result, since in [11] a Navier slip boundary condition and ∂Ω ∈ C 1,1 are assumed. (iv) The convergence proof is presented for discretely divergence-free velocity functions. If V n div ⊂ W 1,∞ 0,div (Ω) d , then no modification of the convective term is required and the proof that u m is divergence-free is also simpler.
The rest of this section consists of the proof of Theorem 4.1, which relies on Lemmas 4.3-4.5 dealing with the existence of the discrete solution, and the limit k → ∞, Lemmas 4.6 and 4.7 covering the combined limit l, n → ∞, and Lemmas 4.8 and 4.9 the limit m → ∞.
Limit k → ∞. The existence and convergence in Lemmas 4.3 and 4.4 follow by a standard approach presented, e.g., in [39], with minor modifications required to deal with the time-dependence of S S S k . Taking k → ∞ we remain in the finite-dimensional setting and hence strong convergence of the sequence of symmetric gradients follows. Consequently, the identification of the limiting equation is based on the properties of the sequence {S S S k } k∈N according to Assumption 2.3; c.f. [11].
Step 1: A priori estimates. The a priori estimates follow from standard arguments, see [39], in combination with the estimates available for S S S k by Assumption 2.3: testing (4. since the term involving b vanishes by skew-symmetry. By the fact that 2a .17), the first term in (4.10) can be rewritten as (4.11) Using the definition of S S S k i in (4.3) and Assumption 2.3 (α2) one has that where the last inequality follows by Korn's and Poincaré's inequality. On the term on the right-hand side of (4.10) by duality of norms and by Young's inequality with ε > 0 we obtain that where the last inequality follows by (4.4). Applying the estimates (4.11)-(4.13) in (4.10), after rearranging, choosing ε > 0 small enough and multiplying by δ l , we arrive at (4.14) For arbitrary j ∈ {1, . . . , l}, summing over i ∈ {1, . . . , j}, yields because of cancellation in the first term. Applying the estimate and taking the supremum over all j ∈ {1, . . . , l} in (4.15) finishes the proof of (4.9).
Step 2: Existence of {U κ i } i∈{0,...,l} . Let us fix κ ∈ N 4 and since U κ 0 = P n div u 0 by (4.7), we only have to show that for a given this amounts to finding α ∈ R dn such that F (α) = 0. Considering term by term one can see that F is continuous in α. By use of an L 2 -orthonormal basis of V n div one can show that Furthermore, we find that (4.20) The first term vanishes thanks to the skew-symmetry and the other terms can be estimated similarly as in (4.11)-(4.14) to obtain Thus, there exists an R > 0 such that F (α) · α > 0, for any α ∈ R dn with |α| = R. So F is outward normal on ∂B R (0), and hence, as a consequence of Brouwer's fixed point theorem, F has a zero in B R (0), see [25, § 5.7, (G.7), p. 104]. This means that U satisfying (4.17) exists and the claim is proved.
Step 1: Identification of the equation. We have that U κ (0, ·) = U κ 0 = P n div u 0 by definition of U κ and by (4.7), which shows (4.24). The equation (4.23) follows from (4.8) and the fact that for t ∈ (t i−1 , t i ] we have that . . , l}. Step 2: Estimates. Let λ = (l, n, m) ∈ N 3 be arbitrary, but fixed. When taking k → ∞ we stay in the finite-dimensional setting, hence it suffices to focus on estimates for the coefficient functions α κ = α k,λ and α κ = α k,λ , uniformly in k ∈ N. Using the definition of U k,λ in (4.21) with the estimate (4.19), the fact that the interpolants are piecewise constant and also the a priori estimate in (4.9), we obtain for any t ∈ (0, T ) uniformly in k ∈ N. The corresponding estimate holds for | α k,λ (t)| 2 , using also that U k,λ (0, ·) = P n div u 0 by (4.24) and the L 2 -stability of P n div in (3.9). Thus, we have that ≤ c for all k ∈ N.  From the a priori estimate (4.9) it follows directly that ≤ c for all k ∈ N.
for all k ∈ N. This also shows that for any k ∈ N and any i ∈ {1, . . . , l}.
Step 2: Identification of the limiting equation. Let W ∈ V n div and let ϕ ∈ C ∞ 0 ((0, T )) be arbitrary but fixed. With the convergence of ∂ t U k,λ in (4.29) it follows that , as k → ∞. In particular, the strong convergence in (4.25) and in (4.27) allows us to take the limit in the numerical convective term without any restriction. Finally, (4.55) and (4.56) applied in (4.23) imply that (4.50) holds for a.e. t ∈ (0, T ). Recall that S S S λ (t, ·) = S S S λ i , for any t ∈ (t i−1 , t i ] and i ∈ {1, . . . , l} by (4.32). Since now the terms in (4.50) are constant on each interval (t i−1 , t i ], for i ∈ {1, . . . , l}, the equation also holds for all t ∈ (0, T ]. Also, (4.54) follows from (4.50) since ∂ t U λ (t, ·) = d t U λ i , U λ (t, ·) = U λ i and the corresponding holds for S S S λ i and f i , for any t ∈ (t i−1 , t i ] and i ∈ {1, . . . , l}. for all ϕ ∈ C ∞ 0 (Q) such that ϕ ≥ 0 and for all matrices B B B ∈ U . By Lemma 2.2 (a4) this allows us to conclude that (D D DU λ (z), S S S λ (z)) ∈ A(z) for a.e. z ∈ Q, so (4.52) is shown.
Limit l, n → ∞. We are taking the limits l, n → ∞ simultaneously without imposing any condition on δ l and h n . The condition q > 2d d+2 is required to gain compactness. Two additional difficulties, compared to [11], arise from the discretization. The first is that in order to prove a uniform bound on the sequence of approximations to the time derivative one would require the stability of the L 2 -projector onto V n div in Sobolev norms, which would impose stronger requirements on the finite element partition of Ω. To avoid this, instead of the Aubin-Lions lemma we shall employ an alternative compactness result due to Simon (cf. Lemma 2.11), which requires convergence properties of time-increments. The second difficulty is that, in the identification of the implicit relation we have to deal with the discrepancy between S S S λ and S S S λ , since S S S λ appears in the equation (4.50) and S S S λ satisfies the implicit relation in (4.52).

(4.59)
Furthermore, for each m ∈ N there exists a u m ∈ L ∞ (0, T ; L 2 div (Ω) d ) ∩ X q div (Q), S S S m ∈ L q ′ (Q) d×d and subsequences such that, as l, n → ∞, Proof.

(4.71)
The other terms follow immediately and (4.59) is proved.
Step 2: Estimates on the discrete level. Similarly as in the proof of Lemma 4.3, testing (4.54) with The first term on the left-hand side is bounded as in (4.11) and the term on the right-hand side is bounded as in (4.13). The only difference arises in bounding the term involving the stress tensor; cf.
, where again Poincaré's and Korn's inequalities were used. Following the same procedure as in (4.14)-(4.16) we arrive at max j∈{1,...,l} ≤ c for all λ = (l, n, m) ∈ N 3 .  Step 3: Estimates on the continuous level. By the definition of the piecewise constant interpolant according to (3.18) it follows from the discrete estimates that   For the estimates of the continuous, piecewise affine interpolant U λ according to (3.19) one also has to estimate the corresponding norms of U λ 0 ; by (4.53) and the stability of the L 2 -projector in (3.9) we have that Together with the discrete estimate in (4.74) this yields that  For the compactness argument instead of U λ we consider U λ ∈ C([0, T ]; V n div ) defined by ≤ c, (4.81) for all λ = (l, n, m) ∈ N 3 .

By the fact that
. . , l}, and by the discrete estimate (4.74) we obtain   Step 4: Convergence of the time increments (compare [12, pp. 174]). Instead of applying the Aubin-Lions lemma, as in [11], here we apply the compactness result due to Simon, stated in Lemma 2.11. This means that we do not need uniform bounds on the time derivatives but only convergence properties for time increments, which avoids the use of stability results in Sobolev norms for the L 2 -projector onto V n div . We wish to apply Lemma 2.11 to the sequence { U l,n,m } l,n∈N , for fixed m ∈ N, with X = W 1,q (Ω) d , B = L 2 (Ω) d and p = 2. Let us show that ds → 0, as ε → 0, uniformly for l, n ∈ N. Consider the term U λ (s + ε, ·) − U λ (s, ·), W Ω , for W ∈ V n div , s ∈ (0, T ) and ε > 0 such that s + ε < T . If s + ε ≤ δ l , then we have U λ (s + ε) = U λ (s) = U λ 1 , so the term vanishes. Now let s + ε > δ l . By the definition of U λ in (4.80) we have that U λ (s, ·) = U λ (max(s, δ l ), ·). By the continuity of U λ and since ∂ t U λ is integrable, we obtain where in the last line we have used that U λ (t, ·) and U λ (t, ·) coincide on (max(s, δ l ), s + ε) ⊂ [δ l , T ].
Step 5: Convergence as l, n → ∞. Recall that we have λ = (l, n, m) ∈ N 3 and let m ∈ N be fixed.
By the estimates (4.81) and (4.83) we have that { U l,n,m } l,n∈N is bounded in particular in L 2 (Q) d and L 1 (0, T ; W 1,q (Ω) d ). By the condition that q > 2d d+2 , the embedding W 1,q (Ω) ֒→֒→ L 2 (Ω) is compact and with (4.87) all the assumptions in Lemma 2.11 are satisfied for X = W 1,q (Ω) d , B = L 2 (Ω) d and p = 2. Hence, there exists u m ∈ L 2 (Q) d and a subsequence such that U l,n,m → u m strongly in L 2 (Q) d , as l, n → ∞.
≤ cδ l → 0, as l → ∞. for any p ∈ [1, ∞) and any r ∈ [1, η). By the uniform bounds in (4.76) and (4.77) and the (sequential) Banach-Alaoglu theorem, up to subsequences, we have that The argument that u m is divergence-free follow as in [16, p. 1001]: Let h ∈ L q ′ (Ω) and note that by the Assumption 3.4 on the projector Π n Q we have that Π n Q h → h in particular in L q ′ (Ω), as n → ∞, compare Remark 3.5 (ii). Also, let ϕ ∈ C ∞ 0 (0, T ). By (4.98) we have that div U l,n,m ⇀ div u m weakly in L q (Ω), and hence div U l,n,m , ϕΠ n Q h Q → div u m , ϕh Q , as l, n → ∞. Since U l,n,m ∈ P l 0 (0, T ; V n div ) the left-hand side vanishes for all l, n ∈ N, and hence we have that div u m , hϕ Q = 0 for all h ∈ L q ′ (Ω) and all ϕ ∈ C ∞ (0, T ), so by density u m is (weakly) divergencefree.
By (4.76) it follows that {|U l,n,m | 2q ′ −2 U l,n,m } l,n∈N is bounded in L (2q ′ ) ′ (Q) d and thus, by the (sequential) Banach-Alaoglu theorem there exists a subsequence and By the strong convergence in (4.96), there exists a subsequence, which converges a.e. in Q, and hence we can identify ψ m = |u m | 2q ′ −2 u m , which shows (4.67).
Step 2: Bound on the time-derivative. The distributional derivative of u m satisfies, by definition and using (4.122), that for all w ∈ C ∞ 0,div (Ω d ) and all ϕ ∈ C ∞ 0 ((0, T )), since supp ϕ ⊂ (0, T ). Using this equation we wish to show that ∂ t u m ∈ L τ (0, T ; (X q div (Ω)) ′ ) (not uniformly in m ∈ N), for τ as in (2.3) and X q div (Ω) as in (2.4). For L m as defined in (4.106), using the fact that u m ∈ L 2q ′ (Q) d and S S S m ∈ L q ′ (Q) d×d let us first estimate for all ϕ ∈ C ∞ 0 ((0, T )) and all w ∈ C ∞ 0,div (Ω) d , since τ ′ = max(2q ′ , q). By the density of the respective test function spaces, L m [u m , ·], · (0,T ) represents a bounded linear functional on L τ ′ (0, T ; X q div (Ω)), and thus we have that ∂ t u m ∈ L τ (0, T ; (X q div (Ω)) ′ ) by (4.123) and by reflexivity of the function space. Consequently, ∂ t u m , w Ω is integrable for w ∈ C ∞ 0,div (Ω) d , and thus, we can rephrase (4.123) by the fundamental lemma of calculus of variations in the pointwise sense in time, so (4.107) is proved.
Step 4: Energy identity. Recall that u m ∈ X q div (Q) ֒→ L min(q,2q ′ ) (0, T ; X q div (Ω) d ) and ∂ t u m ∈ L τ (0, T ; (X q div (Ω)) ′ ), where τ = min(q ′ , (2q ′ ) ′ ), and equation (4.107) is satisfied. Because of the lack of integrability in time an approximation procedure by means of mollification in time can be applied to show the energy identity for a.e. s ∈ (0, T ), where also the attainment of the initial datum in the sense of (4.109) is used. The proof follows by a standard procedure and we therefore omit the details; see, e.g., [32, Ch. 2.5].
where the last inequality is based on the following arguments. By Limit m → ∞. In this step we loose the admissibility of the solution as a test function, and we have to use Lipschitz truncation to identify the implicit relation. The availability of the solenoidal Lipschitz truncation allows to simplify the arguments in [11], since no pressure has to be reconstructed. Let us denote µ := min q(d + 2) 2d , q ′ , (2q ′ ) ′ = min (q ′ , τ ) , (4.145) where τ defined in (2.3). Note that since q > 2d d+2 , we have that µ > 1.
Step 1: Estimates. Recall that by (4.134) we have the following energy identity for a.e. t ∈ (0, T ): where we have used Poincaré's and Korn's inequalities in the last line. Similarly as before, we use duality of norms and Young's inequality with ε > 0 to bound for a.e. t ∈ (0, T ) and all m ∈ N. Taking the essential supremum over t ∈ (0, T ) and also applying the parabolic interpolation from Lemma 2.8 shows (4.146).
Step 2: Bound on the time derivative. In order to derive a uniform bound on the time derivative let us estimate L m [u m ; v]. Since no uniform bounds on ||u m || L 2q ′ (Q) are available at this point, we use the bound (3.14) on the convective term, withq as defined in (2.1), to deduce that | u(t, ·) ⊗ u(t, ·), ∇v Ω | ≤ c ||u(t, ·)|| 2 L q(d+2) d
where S S S ∈ L q ′ (Q) d×d is the limiting function introduced in Lemma 4.8. i.e., (u, S S S) is a weak solution according to Definition 2.1.

Proof.
Step 1: Identification of the limiting equation. Let w ∈ C ∞ 0,div (Ω) d and ϕ ∈ C ∞ 0 ((0, T )) and let us consider each of the terms in (4.107) and (4.162). By the weak convergence results in (4.151) and (4.152) we have that as m → ∞. Since by (4.147) we have that u m → u in L r (Q) d×d for all r ∈ 1, q(d+2) d it follows that d+2 , this set is nonempty and the convergence holds in particular in L 1 (Q), so as m → ∞. Taking the results in (4.165)-(4.167) and (4.153) shows that (4.107) implies (4.162).
Step 2: Identification of the initial condition. With similar arguments as in the proof of Lemma 4.7, Step 3, it follows that u ∈ C w ([0, T ]; L 2 div (Ω) d ), that u 0 = u(0, ·) ∈ L 2 div (Ω) d and that the initial datum is attained in the sense of (4.164).
Step 3: Higher integrability of the time derivative. As in Step 2 in the proof of Lemma 4.6 we can improve the integrability of ∂ t u using the fact that (4.162) is satisfied. This yields that ∂ t u ∈ Lq ′ (0, T ; (W 1,q 0,div (Ω) d ) ′ ), forq as defined in (2.1).
Step 4: Identification of the implicit relation (compare [11] and [8,Sec. 3]). Recall that D D Du m ⇀ D D Du weakly in L q (Q) d×d by (4.149), that S S S m ⇀ S S S weakly in L q ′ (Q) d×d by (4.152) and that we have that (D D Du m (z), S S S m (z)) ∈ A(z) for a.e. z ∈ Q by (4.108). Hence, by Lemma 2.6, it suffices to show that Since there is no energy identity available for u, in order to identify the implicit relation one has to truncate the elements of the approximating sequence of velocity fields suitably so as to be able to use them as test functions. In contrast with [11] we will not use a parabolic Lipschitz truncation after locally reconstructing the approximations to the pressure, but work with the solenoidal Lipschitz truncation introduced subsequently in [8] and stated in Lemma 2.7, as the argument is then more direct.
We wish to truncate v m := u m − u, which satisfies, for all ξ ∈ C ∞ 0,div (Q) d , the equality by (4.107) and (4.162) and by the density of In order to rewrite the equation in divergence form as required for Lemma 2.7, we adapt the last term locally.
With the aid of the parabolic solenoidal Lipschitz truncation we show that      [2]) there exists a further subsequence, a function H ∈ L 1 ( 1 8 Q 0 ) and a nonincreasing sequence of measurable subsets , as m → ∞, for each fixed i ∈ N. Now let i ∈ N be arbitrary but fixed. With the Dunford-Pettis compactness criterion (see [19,Ch. 8,Thm. 1.3]) it follows that the sequence {H m } m∈N is equi-integrable on 1 8 Q 0 \E i . By the a.e. convergence of H m to zero in particular on 1 8   This shows (4.168) for Q = 1 8 Q 0 \E i , and thus we find that (D D Du(z), S S S(z)) ∈ A(z) for a.e. z ∈ 1 8 Q 0 \E i . Since |E i | → 0, as i → ∞, we have that (D D Du(z), S S S(z)) ∈ A(z) for a.e. z ∈ 1 8 Q 0 . Finally let us consider a cover of Q consisting of (open) cylinders Q j = I j × B j , j ∈ J, for an index set J such that Q = j∈J 1 8 Q j . This can be, for example, chosen as a Whitney type cover, compare, e.g., [18]. Then we can identify the implicit relation a.e. on 1 8 Q j for all j ∈ J by the above and thus, have that (D D Du(z), S S S(z)) ∈ A(z) for a.e. z ∈ Q, which proves (4.163).
Remark 4.10 (Steady Problem). The same regularization and splitting approach can be applied to show convergence in the steady case to cover the range q ∈ ( 2d d+2 , 2d d+1 ] for discretely divergence-free finite element functions, which is missing in [16]. In [16] no regularization term was used, and hence the restriction q > 2d d+1 was required in the case of discretely divergence-free finite element functions. Considering the two results together, for q ∈ ( 2d d+2 , 2d d+1 ] one should either use exactly divergence-free finite element functions as was done in [16], or introduce a regularization term and pass to the limit m → ∞ with the regularization parameter, as we have done here.

Appendix: Results about the Constitutive Laws
Lemma 5.1. (Properties of S S S k ) For each k ∈ N the function S S S k : Q × R d×d sym → R d×d sym defined in (2.11) is measurable with respect to its first argument and smooth with respect to the second argument. Furthermore, the sequence {S S S k } k∈N satisfies the Assumption 2.3 with g = g, c * = c * and q ∈ (1, ∞), as in Lemma 2.2.
Proof. Smoothness and measurability follow by the definition of the convolution and Fubini's theorem. To show the properties (α1) and (α2) is straightforward. For (α3) we will follow [11,Sec. 3.2]: Let z ∈ Q be arbitrary but fixed and not in any of the zero-sets for which the properties of S S S ⋆ do not hold. Then we first use the definition of S S S k in (2.11) and the fact that ρ k integrates to 1 on R d×d sym and then the monotonicity of S S S ⋆ in Lemma 2.2 (a2), which gives that where the constant c is independent of z ∈ Ω and of k ∈ N.
Then, since S S S ⋆ (z, ·) is locally bounded (compare Lemma 2.2 (a5)), there exists a constant c > 0 independent of k and z such that where we have again used that ρ k integrates to 1. Now we take (5.2) and (5.4) together, multiply by ϕ ∈ C ∞ 0 (Q) such that ϕ ≥ 0, and then integrate in Q. Recalling that the constant is independent of z ∈ Q this givesˆQ  . The family of functions S k : Q × R → R, k ≥ k 0 , defined in (2.12) has the following properties: each S k is a Carathéodory function and (α1') for a.e. z ∈ Q the function S k (z, ·) is monotone; (α2') there exists a constant c * > 0 and a nonnegative g ∈ L 1 (Q) such that S k (z, B)B ≥ − g(z) + c * |B| q + S k (z, B) q ′ , for any B ∈ R ≥0 and for a.e. z ∈ Q and all k ≥ k 0 ; (α3') for any sequence {D k } k∈N bounded in L ∞ (Q), for any B ∈ R ≥0 and all ϕ ∈ C ∞ 0 (Q) such that ϕ ≥ 0, we have lim inf k→∞ˆQ S k (·, D k ) − S ⋆ (·, B) : D k − B ϕ dz ≥ 0.
Proof. Since S ⋆ (·, ·) is measurable in the first argument by Lemma 2.2, so is S k (·, ·) and by construction S k (·, ·) is continuous in the second argument. Hence, S k (·, ·) is a Carathéodory function.
(α1') By Lemma 2.2 (a2), S ⋆ (z, ·) is monotone for a.e. z ∈ Q. Then, by construction S k (z, ·) is piecewise monotone for a.e. z ∈ Q, and consequently monotone for a.e. z ∈ Q. (α2') Let z ∈ Q be arbitrary but fixed such that all properties hold for S ⋆ (z, ·), excluding zero sets where necessary. Also let B ∈ R ≥0 and k ≥ k 0 be arbitrary but fixed. If B / ∈ A k , then we have that S k (z, B) = S ⋆ (z, B) and hence the claim holds by the estimate for S S S ⋆ (z, ·) in Lemma 2.2 (a3), which is equivalent to the one for S ⋆ (z, ·), with the same g and c * . If B ∈ A k , we make use of the fact that A and hence also A k is bounded. In particular we know that B ≤ a k I,+ ≤ a I + 1. Hence, again with the estimate for S ⋆ (z, ·) corresponding to Lemma 2.2 (a3) we obtain . Now using that B < a I + 1, and monotonicity of both S ⋆ (z, ·) and S k (z, ·) we find that |S ⋆ (z, B)| ≤ |S ⋆ (z, a I + 1)| , S k (z, B) ≤ S k (z, a k I,+ ) = S ⋆ (z, a k I,+ ) ≤ |S ⋆ (z, a I + 1)| . In order to estimate |S ⋆ (z, B)| q ′ further, we use the inequality (a + b) s ≤ 2 s−1 (a s + b s ), for a, b ∈ R ≥0 , s ∈ [1, ∞).
By the triangle inequality and again applying (5.6), this yields Rearranging and applying this in (5.7) gives S k (z, B)B ≥ − 2 |S ⋆ (z, a I + 1)| (a I + 1) − g(z) + c * |B| q + |S ⋆ (z, B)| q ′ ≥ − 2 |S ⋆ (z, a I + 1)| (a I + 1) − g(z) − c * 2 q ′ |S ⋆ (z, a I + 1)| q ′ + c * |B| q + 2 −(q ′ −1) c * S S S k (z, B) We set g(z) := 2 |S ⋆ (z, a I + 1)| (a I + 1) + c * 2 q ′ |S ⋆ (z, a I + 1)| q ′ + g(z) ≥ g(z) ≥ 0, which is in L 1 (Q) thanks to the local boundedness of S ⋆ and we choose c * := 2 −(q ′ −1) c * . (α3') Let {D k } k∈N be bounded in L ∞ (Q) by C, let B ∈ R ≥0 and let z ∈ Q be arbitrary but fixed such that all properties hold for S ⋆ (z, ·), (α1') and (α2') hold for S k (z, ·) and such that D k (z) ≤ C, excluding possibly zero sets. We wish to show that there exists a constant c > 0 depending on S S S ⋆ , a I , B and C such that To this end, we distinguish the following three cases. B / ∈ A k : In this case we have that S k (z, B) = S ⋆ (z, B) and by monotonicity of S k (z, ·), it follows that (⋆) ≥ 0. D k (z) / ∈ A k : Similarly, in this case we have that S k (z, D k (z)) = S ⋆ (z, D k (z)), so by monotonicity of S ⋆ (z, ·) it follows that (⋆) ≥ 0. B, D k (z) ∈ A k : Assume that B ∈ A k i and D k (z) ∈ A k j and further distinguish the cases i = j and i = j. i = j: In this case we make use of the fact that D k (z) − B ≤ 2 k and that A k is bounded, i.e., B, D k (z) ≤ a I + 1. Then, using again the estimates in (5.6) we have ≥ −2S ⋆ (z, a I + 1) 2 k ≥ − c k , since S ⋆ is locally bounded. i = j: In this case we make use of the fact that S k (z, ·) and S ⋆ (z, ·) agree on a k i,± , a k j,± and are monotone. Let i < j, i.e., B ≤ D k (z); it then follows that S ⋆ (z, B) ≤ S ⋆ (z, a k i,+ ) = S k (z, a k i,+ ) ≤ S k (z, a k j,− ) ≤ S k (z, D k (z)), which implies that S k (z, D k (z)) − S ⋆ (z, B) : D k (z) − B ≥ 0.
The case i > j is shown analogously. Altogether, these imply (5.8), where the right-hand side is independent of z. Multiplying with ϕ ∈ C ∞ 0 (Q) such that ϕ ≥ 0, and integrating over Q giveŝ Q S k (·, D k ) − S ⋆ (·, B) : D k − B ϕ dz ≥ − c k .
Then applying lim inf k→∞ yields the assertion. where S k is defined in (2.12), satisfies Assumption 2.3.
Proof. It is straightforward to show that monotonicity of S k implies monotonicity of S S S k and that the growth and coercivity bounds are equivalent for S k and S S S k .