Error estimates for the Smagorinsky turbulence model: enhanced stability through scale separation and numerical stabilization

In the present work we show some results on the effect of the Smagorinsky model on the stability of the associated perturbation equation. We show that in the presence of a spectral gap, such that the flow can be decomposed in a large scale with moderate gradient and a small amplitude fine scale with arbitratry gradient, the Smagorinsky model admits stability estimates for perturbations, with exponential growth depending only on the large scale gradient. We then show in the context of stabilized finite element methods that the same result carries over to the approximation and that in this context, for suitably chosen finite element spaces the Smagorinsky model acts as a stabilizer yielding close to optimal error estimates in the $L^2$-norm for smooth flows in the pre-asymptotic high Reynolds number regime.


Introduction
The modelling and accurate approximation of turbulent flows at high Reynolds numbers remain an outstanding challenge. One of the earliest attempts at the design of a turbulence model is the one attributed to Smagorinsky [36,37], where a nonlinear viscosity replaces the Reynolds stresses, thereby providing closure to the filtered Navier-Stokes' equations. The model is believed to be inspired by early work on shock capturing methods for conservation laws due to Von Neumann and Richtmyer [39].
An appealing feature is that it has been shown that the Navier-Stokes' equations regularized system admits a unique solution [31], with enhanced regularity [5,6]. It has also been shown using scaling arguments that away from boundary layers the dissipation rate matches the Kolmogorov decay rate [32] for uniformly turbulent flows. On the other hand the Smagorinsky model is considered to be too dissipative, in particular in the laminar zone. The parameter of the model has to be tuned differently to obtained the best results depending the flow problem [22] and on the numerical method used [34], pointing to a subtle interplay between the flow and the numerical method for the performance of the Smagorinsky model in LES. Nevertheless thanks to its simplicity and its physical consistency, the Smagorinsky model remains an important tool for the modelling of turbulent high Reynolds flows.
The objective of the present work is to show some results on the effect of the Smagorinsky model on the stability of the fluid dynamics and on numerical approximation. We show that in the presence of a spectral gap, such that the flow can be decomposed in a large scale with moderate gradient and a small amplitude fine scale with arbitrary gradient, the Smagorinsky model admits stability estimates with exponential growth depending only on the large scale gradient. We then show in the context of stabilized finite element methods that the same result carries over to the approximation and that in this context, for suitably chosen finite element spaces the Smagorinsky model acts as a stabilizer yielding optimal error estimates for smooth flows in the high Reynolds number regime. This gives some evidence that the Smagorinsky model is efficient in situations where the flow has a clearly defined spectral gap. This is not in general the case for turbulent flows, but is believed to occur at high Reynolds number, near rough surfaces or the atmospheric boundary layer [21], or more generally in atmospheric and geophysical flows [20,40,35]. Interestingly, this is the context for which the Smagoringsky model was originally proposed [36].
The notion of Large Eddy Simulation (LES) is not very well defined. Depending on the context it is considered either as a modelling issue, or a numerical technique. In the first case the partial differential equation is perturbed by adding terms that model the so-called Reynolds stress term expressing the effect of unresolved scales, where the smallest length scale is a model parameter. The regularized model should be well-posed and can then be discretized. The challenge here is on the modelling side: how can one design a model that allows for a more stable representation of quantities of interests? Observe that to distinguish this method from the situation where the turbulence model is used for numerical stability, the model parameters should be fixed before discretization and not coupled to the mesh-size. In case such a coupling is introduced it would need to be justified, in particular in the asymptotic limit. A typical requirement is that the physically relevant solution, or a so-called "suitable" solution is obtained in the limit [24]. For so called numerical LES, or Implicit LES, on the other hand [9], the Navier-Stokes' equations are considered, but a numerically stable discretization method is used ensuring that sufficient dissipation is added for the numerical scheme to remain stable, but not so much that smooth structures of the flow are destroyed. The underpinning idea is that any quantity of interest that has sufficient stability properties will then be computable [27,10], if the numerical method resolves a sufficient range of scales so that the effect of unresolved scales is negligible. To be able to claim that an analysis of a LES method is successful it is necessary to address the question of stability, but there appears to be no stability estimates with moderate exponential growth, for any quantity of interest for the Navier-Stokes equations in the high Reynolds number regime. The best rigorous result of quantitative bound of perturbation growth in the inviscid regime appears to be [1], where a two-dimensional flow is analyzed and the stability constant grows with a double exponential in time.
Finite element analysis typically relies on the velocities having bounded gradients [30,26], with a slight improvement possible in the two-dimensional case using parabolic regularization of the viscous scales [10]. Typically the stability constant is on the form exp(KT ), with K proportional to the maximum gradient of the velocity in the flow.
For the discussion below we will for simplicity consider the p-Laplacian model that acts on the whole gradient instead of the symmetric gradient as in the Smagorinsky case. This regularization was introduced in [33, Chap. 2, sec. 5], where also existence and uniqueness of solutions were studied. The below arguments can however easily be extended to the more common versions of the Smagorinsky model using the deformation tensor. The only essential difference is the need to apply Korn's inequality in L p , see for instance [16,Chapter 7] to get control of the gradient of the solution in L p through the control of its symmetric part. Main results and outline. There appears to be very little known about the stability properties of the Smagorinsky model, or for that matter its qualities as a stabilizing term for numerical methods for high Reynolds flow in the laminar regime. In this paper we will consider the Smagorinsky model first, as a turbulence model and then, in a finite element method, as a stabilizing term. In the former case we prove a stability result (Theorem 4.1) showing that the perturbation growth is independent of high frequency, low amplitude oscillations in the solution and that as the perturbation error becomes larger, its growth is moderated by the nonlinear feedback of the turbulence model.
On the other hand the use of the Smagorinsky model as a stabilizer in the framework of divergence free affine finite elements, allows us to prove the classical L 2 error estimates of order O(h 3 2 ) (where h denotes the mesh-size), that are the best known estimates for stabilized finite element methods using affine approximation for the Navier-Stokes' equations in the high Reynolds laminar regime [30,26,29,12]. Indeed if u denotes the solution to the Navier-Stokes equations andû h denotes the finite element solution stabilized using the Smagorinsky model with characteristic lengthscale O(h) and the flow is laminar in the pre asymptotic regime where the mesh-size is larger than the viscosity, then there (Theorem 6.7) where C(u) is a constant that depends on Sobolev norms of the exact solution to the Navier-Stokes' equations and time. In particular there are no inverse powers of the viscosity in C(u). The exponential growth of the constant inherits the scale separation properties of the continuous equations.
Finally we show the performance of the Smagorinsky finite element method qualitatively of some academic test cases, illustrating that the method returns a solution also in underresloved computations where the standard Galerkin breaks down, but also that this solution is sensitive to the choice of the parameter in the Smagorinsky model. A large parameter gives an overly diffusive approximation on coarse meshes.
The main conclusion of this work is that the Smagorinsky model acts both on the level of stability of perturbations and as a numerical stabilizer. This shows that the nonlinear viscosity can be interpreted both as a turbulence model for LES and as an implicit LES model based on stable numerical simulation, provided the numerical method is chosen carefully to balance and complement the dissipation properties.

Smagorinsky model problem
Let Ω ⊂ R d , d = 2, 3, denote an open domain with smooth (or convex polyhedral) boundary ∂Ω. Consider the time interval I = [0, T ] and denote the space time domain Ω × I by Q. We consider the Navier-Stokes-Smagorinsky equations on the form , and ν : R + → R + a monotonically increasing function with ν(0) = 0 and δ a characteristic lengthscale. Below we will consider the choice ν(δ) := δ 2 which corresponds to the well-known Smagorinsky model [37], using the gradient instead of the deformation tensor. Taking δ = 0 on the other hand results in the standard Navier-Stokes' equations. The Smagorinsky model is not dependent on a particular filter and the length scale δ can be allowed to vary in space and time. As a consequence Van Driest wall damping [38] or Germano dynamic models can be included in the arguments below. For simplicity we here consider constant δ and either homogeneous Dirichlet boundary conditions or periodic boundary conditions for (2.1). We introduce the vectorial Bochner space We will assume that the solution to the Navier-Stokes' equations, u, has the additional regularity u ∈ V . Since (2.1) results in an O(δ 2 ) perturbation of the Navier-Stokes' equations, it is reasonable to assume that the convergence from the Smagorinsky-Navier-Stokes' system to the Navier-Stokes' system is at best O(δ 2 ) for a smooth enough solution (to the Navier-Stokes' system). If the solution is not smooth enough the effect of the perturbation will be larger. For instance, as shown below, if u ∈ V the perturbation due to the model is O(δ). However, the error due to the inconsistency is not the whole story, indeed the error can typically be written (see Theorem 4.1 below) where u is the solution to the Navier-Stokes' equations,ũ is the solution to the Navier-Stokes'-Smagorinsky equations, C(u) is a constant only dependent on Sobolev norms of u and α(T ) is a function, depending on the exact solution, quantifying the stability of the nonlinear problem and β is the power of the nonconsistency, typically in the interval (0, 2). The consequence of this for practical computation is that for the situation where δ coincides with the cell-width of the computational mesh, improving the approximation order beyond that of affine finite elements (formally giving an O(h 2 ) L 2 -error bound) will not result in an improved approximation for laminar flows: the consistency error O(δ 2 ) of the Smagorinsky model will dominate.
For laminar flows where the flow field has small spatial variation also the coefficient α(T ), which typically scales as |u|/ν 1 2 or I ∇u(·, t) L ∞ (Ω) , is moderate and one can conclude that the perturbation error is essentially determined by the size of δ in the right hand side. Decreasing δ will lead to a decrease in the error and the natural (and trivial) choice is δ = 0 for which the error is zero, since the two solutions coincide: there is no need for a model in laminar flows.
In the less regular case, u ∈ V , which seems to be the minimum requirement for a bound of the type (2.2) to make sense, the situation is different. Now α can not be assumed small, so the right hand side depends both on the coefficient δ and the exponential growth. Except for very short times the exponential growth will make the bound meaningless.
We will show in the present paper that we can refine our definition of the coefficient α by letting it implicitly depend on δ, so that an increase in δ, although it increases the consistency error, can lead to a decrease in the exponential growth through the nonlinear feedback of the Smagorinsky operator. This leads to a potential moderation of the error growth through an increase in δ, which can explain why this behavior has been observed in computations [22]. Observe that this result depends on the structure of the actual solution u, so it is not a general statement and only of qualitative nature. The analysis is also done in the L 2 -norm which is too strong to be a reasonable target for turbulent flows. Unfortunately the analysis of other target quantities appears to be very difficult due to the complexity of the linearized adjoint to the Smagorinsky perturbation equation.
Nevertheless one can argue that if this stabilizing effect is visible already for the L 2 -norm it is likely to be more accentuated for the approximation of averaged quantities where the effect of fluctuations is expected to be less important due to cancellation.

Notation and technical results.
To reduce the number of non-essential constants we will use the notation a b for a ≤ Cb where C is an O(1) constant independent of the viscosity or the mesh parameter h.
We will use standard notations for most Sobolev spaces and norms, but we will not distinguish in the notation between the scalar, vector and tensor valued cases. The space of divergence free functions in V and [L 2 (Ω)] d will be denoted To consider different boundary conditions we add a superscript 0 on spaces of functions with zero trace on the domain boundary, for example V 0 0 . To simplify we define the L 2 -inner product for some subset X ⊂ R s , s = 1, . . . , d + 1, by with associated norm · X := (·, ·) 1 2 X . With some abuse of notation we will not distinguish between the inner product of scalar, vector or tensor valued functions. The vector valued case is defined by where ξ : ζ := d i,j=1 ξ ij ζ ij . The tensor valued L p norm will be defined by The norm on C 0 (Ω) will be denoted, We also recall Hölders inequality and Young's inequality that will be used repeatedly below. For p, q ∈ [1, ∞] with p −1 + q −1 = 1 there holds The following Poincaré-Friedrichs inequality holds for all 1 ≤ p < ∞.
We will use the notation T for a quasi-uniform tesselation of Ω, consisting of simplices T with diameter h T . We also introduce the global mesh parameter h = max T ∈T h T . We will denote the broken L 2 -scalar product and norm over the elements of T by The set of internal faces of the triangulation T will be denoted F and we define the norm of a function defined on the skeleton by The following monotonicity and contintuity results of the p-Laplacian will be useful for the analysis.
Proof. The result is elementary. Adding and subtracting ζ we have The claim then follows by applying the reverse triangle inequality ||ζ| F − |ξ| F | ≤ |ζ − ξ| F so that ) Ω . Lemma 2.3 also implies the continuity

Perturbation equation and scale separation
We are interested in estimating the growth of perturbations in (2.1) quantified as the difference of the solutionû, with δ > 0 and the Navier-Stokes' equation u, with δ = 0. To this end we will study the perturbation equation obtained by taking the difference of the equation with δ > 0 and with δ = 0 on weak form.
Let η = u −û, then η ∈ V 0 and using (2.1) we see that A classical energy estimate can be obtained from this equation by testing with η. Clearly which (as we shall see below) leads to exponential growth with coefficient I ∇u(t) L ∞ (Ω) using Gronwall's estimate. The objective of the present contribution is to show that the large scale velocity that drives the exponential growth can be taken over a much larger set than {u,û}. Indeed we will introduce a scale separation property defining a set of fine and coarse scales. The coarse scales are dependent on the Smagorsinsky solution implicitly through the difference between the Navier-Stokes' solution and the Smagorinsky solution. Under the scale separation assumption, the exponential growth of perturbations will only depend on the coarse scales. Any exact solution may be decomposed in a large scale componentū and a fine scale fluctuation u ′ so that where τ L is a characteristic time scale of the large scales of the flow defined by For u ∈ [L 1 (I; W 1,∞ (Ω))] d the relation always holds for the trivial case u ′ = 0, leading to the classical exponential growth with coefficient proportional to the full gradient of the fluid velocity. We say that the flow is scale separated if the decomposition (3.2) exists with τ L ∼ T .
Remark 3.1. Observe that the decomposition depends onν(η) pointwise and that this quantity can be very large for rough flows, since ∇u, which is not regularized, may be large locally without violating the a priori regularity assumption. This is a favorable property: as the flow gets more rough the bound defining admissible small scales is relaxed and the characteristic time of the resolved scales can increase, leading to a reduction in the exponential growth of perturbations.
The rationale for the decomposition is that the constant in the estimates below will have exponential growth depending only on T /τ L . This shows that small amplitude, high frequency perturbations of scale separated flows will not grow exponentially, even if ∇u L ∞ (Ω) is large. The following Lemma is a key tool to quantify this observation.
Lemma 3.2. Letû ∈ V be the solution of (2.1) with δ > 0 and let u ∈ V be the solution with δ = 0. Assume that u satisfies the scale separation property (3.2). Then, denoting η = u −û, there holds, for all ǫ > 0, Proof. Using the divergence theorem and the scale separation property we see that Using the bound of (3.2) we see that, since by assumption By the definition ofν(η) An application of the Cauchy-Schwarz inequality followed by the inequality (2.5), with p = q = 2 then shows that Collecting the above inequalities we have The result follows by noting that

Perturbation growth on the continuous level for scale separated flows
We will first prove, using a perturbation argument on the continuous equations, that the perturbation induced by the Smagorinsky term allows for an O(ν(δ) 1 2 ) error estimate in the L 2 -norm between the Navier-Stokes' solution and the Navier-Stokes-Smagorinsky solution. In case the solution is smooth this bound can be improved to O(ν(δ)).
Theorem 4.1. Let u ∈ V be the solution to (2.1) with δ = 0 and letû be the solution of (2.1) with δ ≥ 0. Then there holds, with η : Proof. Taking v = η in (3.1) we see that, using the skew symmetry of the convective term and the monotonicity of the p-Laplacian, Lemma 2.6, we get We can bound the right hand side using Hölders inequality (2.4) and Young's inequality (2.5), with p = 3/2 and q = 3, The second term in the right hand side is now absorbed by the fourth term in the left hand side of (4.1). After integration in time we obtain Applying Lemma 3.2 we proceed to bound the second term in the right hand side which is the last term that does not have sign.
Taking ǫ = 1/6 and multiplying through with 2 we obtain the bound Rewriting this as we conclude by applying Gronwall's inequality, leading to where we used that 2 This proves the first inequality. To prove the second observe that proceeding by integration by parts we have instead of (4.2), The conclusion once again follows using Gronwall's Lemma.
Remark 4.2. Writing out the bound on u ′ of (3.2) pointwise, with ν(δ) = δ 2 and µ = 0, leads to From this it follows that increasing δ must have one (or a combination) of the following consequences: (1) the error in ∇(u −û) decreases, the scale separation stays the same; (2) the characteristic time τ L increases and as a consequence the exponential growth decreases; (3) the small scale |u ′ | increases; which implicitly allows for a decrease inū through the definition of the scale separation and a possible increase in τ L .
All these three possibilities point to the fact that either the error is reduced by the increase of the parameter, or the exponential coefficient in the estimate will decrease leading to a decrease in the exponential growth. This effect can offset the effect of the increased consistency introduced by increasing δ. Another salient conclusion is that if δ is coupled to the mesh size, there may be flow configurations where the computational error grows under mesh refinement. Although the consistency error decreases and numerical resolution increases, the set of non-essential fine scale decreases leading to increased exponential growth of perturbations. This hints at a resolution barrier for Smagorinsky LES beyond which a full DNS is required to enhance accuracy further.
The result of Theorem 4.1 shows that the exponential growth of perturbations can be moderated by the Smagorinsky term and hence the turbulence model indeed has a stabilizing effect, which was the first objective of the present work.

The Smagorinsky model as a numerical stabilizer
In this section we will consider the situation where the Smagorinsky term is not a physical model, but a stabilizing term in a numerical method. We will consider a piecewise affine finite element method that fits in a discrete de Rham complex, see for instance [41,25,14]. That means that the space has a divergence free subspace with optimal approximation properties. It turns out that the discrete de Rham complex and piecewise affine approximation are exactly the properties that make the Smagorinsky model a stabilizing term, with well balanced stability versus accuracy. The affine approximation order makes the consistency error of the nonlinear viscosity similar to the accuracy of the approximation and the additional stability obtained through the exact satisfaction of the divergence free condition reduces the need of stabilization so that the second order Smagorinsky term is sufficiently large to control fluctuations inside elements. To counter instabilities due to the lack of C 1 -continuity of the approximation space a penalty term is added on a certain component of the jump of the streamline derivative. Observe that this latter stabilizing term echos the early approach to stabilization of divergence free elements proposed in [13]. In this work however no improved convergence rate was obtained as a consequence of the stabilization. To the best of our knowledge, the analysis below gives the first Reynolds number robust L 2 -error estimate with O(h 3 2 ) convergence for a piecewise affine H 1 -conforming finite element method, satisfying the incompressibility constraint exactly. The analysis draws on recent results using Galerkin-Least Squares stabilization of the vorticity equation [2]. The analysis presented there however does not carry over to the piece affine case considered herein. For the corresponding result using H div -conforming methods we refer to [4].

Approximation space and technical results.
For the purposes of the present paper it is sufficient to know that the space allows for a divergence free subspace with optimal approximation properties in H 1 and L 2 . The below analysis also uses that the space is affine to achieve the optimal error estimate. For higher order spaces the weak consistency of the Smagorisky model is insufficient, even if the flow is laminar. We let V h denote the velocity subspace of piecewise affine vector functions in [H 1 (Ω)] d constructed such that it satisfies a discrete inf-sup condition for the divergence free constraint using some pressure space Q h such that ∇ · V h ∈ Q h . Observe that by working in the divergence free space the pressure can be eliminated in the analysis, which we will make use of to reduce the technical detail below and we will therefore not discuss Q h further but refer to [14,25]. Below we will assume that Ω ⊂ R 3 is convex, simply connected, polyhedron. The two-dimensional case is also covered by the analysis, but to reduce notation we omit it from the discussion.
The following classical inverse and trace inequalities are used frequently in the analysis: • Inverse inequalities, Here P 1 (T ) denotes the set of polynomials of degree less than or equal to 1 on the simplex T . For p ≥ q ≥ 1, l = 0, 1, there holds For a proof of (5.1) see [

Approximation error estimates.
To simplify the analysis we introduce the divergence free space with a homogeneous Dirichlet condition on the normal component, we also define the divergence free subspace of V h We introduce the L 2 -orthogonal projections π h : [L 2 (Ω)] d → V 0h and Π h : L 0 → V 0h . By the quasi uniformity of the mesh the following bounds hold using standard finite element approximation arguments We also recall the L p -stability of π h , for p ≥ 1 there holds [15], We To see this, note that using the regularity of v and a trace inequality (5.3) on each face F shared by simplices T 1 and T 2 , we have The following Lemma is an immediate consequence of (5.6).
Proof. By the triangle inequality followed by an inverse inequality, (5.2) and the stability of the L 2 -projection (5.6), we have which completes the proof.
We will now prove some estimates in weaker norms that will be helpful in the analysis. We recall the following surjectivity property of the curl operator on simply connected polyhedral domains ω. For all v ∈ V 0n there exists ϕ ∈ V 0 (ω) ∩ [H 1 (ω)] d such that ∇ × ϕ = v in ω, ϕ × n ω = 0 on ∂ω and satisfying the stability ∇ϕ ω ∇ × ϕ ω . [3, Theorem 3.17] Using the above properties we now introduce the regularized approximation error E ∈ L 0 (Ω) ∩ [H 1 (Ω)] d defined by in Ω, (5.10) n × E = 0 on ∂Ω. E Ω + h Proof. By our assumptions on Ω there holds ∇E Ω ∇ × E Ω and there exists Ψ ∈ L 0 (Ω) ∩ [H 1 (Ω)] d and Ψ · n = 0 on ∂Ω such that ∇ × Ψ = E, with ∇Ψ Ω ∇ × Ψ Ω [3, Theorem 3.12]. Using the E and Ψ we may bound the error in the following way We conclude that On each face we use the trace inequality (5.3) It follows that (5.14) h v − Π h v 2 Ω . Multiplying (5.14) by h and using approximation in the right hand sides of (5.13) and (5.14) leads to the desired inequality.

Vector identity.
For the analysis below the following elementary vector identity [4] will be useful. For two 3 × 3 matrices A and B with rows A i and B i (i = 1, 2, 3) we define the vector quantity c := A × B with the components A simple calculation then gives the following identity.

A linear model problem.
Since the approach to stabilization in the present work is nonstandard we first consider the linear model problem for inviscid flow introduced in [4]. The Smagorinsky bulk term is not strictly necessary for the analysis in the linear case. Nevertheless we keep this term, since in this simplified context the stabilizing mechanisms become clear. The model problem takes the form, find a velocity u and a pressure p satisfying ∇ · (u ⊗ β) + σu + ∇p =f in Ω, (5.15a) in Ω, (5.15b) u · n =0 on ∂Ω. (5.15c) We think of u and β as column vectors and we set u ⊗ β = uβ t . We assume that β ∈ [C 1 (Ω)] d ∩ V 0n and σ ∈ R + . To ensure uniqueness we assume σ sufficiently big compared to ∇β ∞ , for details see [4].
The numerical method we will analyze here reads: Find u h ∈ V 0h such that where γ > 0 is a dimensionless parameter. We define the stabilizing operator as a combination of a face penalty operator and an artificial viscosity term in the bulk. In the fully nonlinear case the bulk viscosity will be replaced by the Smagorinsky model, We define the stabilization semi-norm by |v| s := s(v, v) 1 2 and note that the following approximation estimate holds Proof. First note that for the linearized Smagorinsky term there holds by (5.5) For the term on the faces on the other hand we have ] F . The first claim now follows by summing over all faces, using (5.8). The second claim is immediate from the first claim and the assumptions.
Clearly by taking v h = u h in (5.16), integrating by parts and using the properties of β and a Cauchy-Schwarz inequality we have the a priori estimate Thanks to equation (5.17) the problem (5.16) has a unique solution. Moreover, the method (5.16) is weakly consistent; in fact, for (u, p) ∈ [L 2 (Ω)] d ×L 2 0 (Ω) solving (5.15) we have, using that (p, ∇·v h ) Following the ideas of [4] we now prove an a priori error estimate in the L 2 -norm for the finite element solution to (5.16).
Proposition 5.5. Let u ∈ L 0 ∩ [H 2 (Ω)] d be the solution to (5.15), and u h ∈ V 0h the solution to (5.16). Also assume that ∇β ∞ δ −1 β ∞ , γ = O(1) and δ = O(h). Then there holds Proof. Let η = u − u h and observe that since by definition ∇ · β = 0 and β · n = 0 on ∂Ω σ 1 2 η 2 Ω + |η| 2 s = −(η, β · ∇η) Ω + (ση, η) Ω + s(η, η). We can then use (5.16) and (5.18) and the divergence theorem to obtain For the second term on the right hand side we have (5.20) (ση, u − Πu) Ω ≤ σ 1 2 (u − Πu) 2 Ω , and for the last three terms It follows using the Cauchy-Schwarz inequality and the arithmetic-geometric inequality that For the boundary term we use the Cauchy-Schwarz inequality followed by a global and local trace inequality and approximation to obtain For the second inequality we used that by (5.3), u − Πu ∂Ω h − 1 2 u − Πu Ω + h 1 2 ∇u − Πu Ω and using also (5.1) it is straightforward to prove that It only remains to bound the convective term. To this end we use the regularized error defined by (5.9)-(5.11) and write (β · ∇η, u − Πu) Ω = (β · ∇η, ∇ × E) Ω (5.23) Using now that the finite element functions are affine per element and Lemma 5.3 we have For the first term of the right hand side we see that The second term of the right hand side of (5.24) is bounded in the following way, using the Cauchy-Schwarz inequality, Young's inequality and Lemma 5.2 Finally for the last term of the right hand side of (5.23) we see that . We finish the proof by applying (5.5), Lemma 5.4 and Lemma 5.2 in the right hand side of (5.28) and the assumption on β.
Remark 5.6. In the case of the linear model problem the addition of the linearized Smagorinsky term is not strictly necessary to obtain the order h 3 2 in the estimate of Proposition 5.5. Indeed the bound of equation (5.26) can be modified as follows to obtain the convergence with δ = 0.
Adding and subtracting π h η, using the triangle inequality and the inverse inequality (5.1) gives Then by approximation and stability of the L 2 -projection, Using this in (5.29) followed by an arithmetic-geometric inequality leads to The first term in the right hand side can now be absorbed by the zero'th order term in the left hand side and the remaining terms are of order O(h 4 ), by Lemma 5.2, which shows that they are sufficiently small by a margin of O(h). In the nonlinear case below, β will be replaced by the discrete solution and it appears no longer to be possible to balance the estimate to optimal order without the nonlinear stabilization. So independent of the improved stability through scale separation, the nonlinear stabilization is of interest for the finite element discretization of the Navier-Stokes' equations.

Stabilized finite element method for the incompressible Navier-Stokes' equations
The finite element space semi-discretization for the approximation of (2.1) reads: findû h ∈ V 0h , withû h (0) = Πu(0) such that Here the term bc is a consistency term added due to the fact that only the normal component of the velocity is set to zero in V 0h defined by where t = I − n⊗ n is the projection onto the tangential plane of the boundary of Ω. The stabilization term s is defined by The parameters γ i , i = 0, 1 are dimensionless, positive, real numbers and U = O(1) is some characteristic velocity of the flow.
The following Lemma is useful for the analysis in the presence of bc. We give the proof of this results in appendix.
The formulation (6.1) corresponds to a dynamical system. It admits a unique solution as shown in the following propositions. Proposition 6.2. For γ 1 large enough, a solutionû h to the system (6.1) satisfies the following stability estimate for all T > 0, where the constant C is independent of T .
Proof. Testing (6.1) with v h =û h yields, for all t ∈ I, Noting that the second term is zero by skew-symmetry we have after integration on (0, r) for all r ∈ I and using Lemma 6.1 Taking now the supremum over r ∈ I in the first term of the left hand side we see that and therefore Finally we see that which finishes the proof.
Proposition 6.3. The system (6.1) admits a unique solutionû h ∈ V 0h on the interval (0, T ], T > 0. Proof. Let N = dim V 0h , U : R → R N , F : R N → R N and let {ϕ i } denotes a basis of V 0 0h . The system (6.1) is equivalent to a dynamical system The matrix M corresponds to the mass matrix of the velocity finite element basis, defined blockwise by M ij := (ϕ j , ϕ i ) Ω . By inspection the function F (U ) is locally Lipschitz. Indeed, for all U , Y ∈ R N there holds To verify this we only need to consider the nonlinear terms. Let u h and y h denote the functions in V 0h associated to the vectors of unknowns U and V and let w h = u h − y h . Then we have The convective term and the face oriented contribution to the stabilization both are [C 1 (R N )] N and the Lipschitz continuity follows from the mean value theorem. For the Smagorinsky term on the other hand it is immediate by the inequality (2.9). It follows that for every h > 0 there exists a T > 0 such that (6.1) admits a unique solution on (0, T ]. To extend this to arbitrary time intervals we need a bound on u h L ∞ (Ω) . We apply an inverse inequality u h (·, T ) L ∞ (Ω) ≤ Ch − d 2 u h (·, T ) Ω . Recalling Lemma 6.2 it follows that for fixed h the solution is globally bounded. This proves the claim.
Remark 6.4. A consequence of the error analysis below is that the L ∞ -bound on the discrete solution actually holds independently of h provided the exact solution is smooth. For details see Corollary 6.8 6.1. Error analysis with exponential growth moderated through scale separation. We will now prove an error estimate for the approximation of the regularized equation. We are interested in underresolved flow so in what follows we assume that µ ≤ U h. Observe that the arguments are valid also for µ = 0, but in this case the no-slip condition must be relaxed on the continuous level. Moreover the same analysis may be used to prove an optimal estimate of O(h) in the H 1 -norm if µ ≥ U h. To moderate the exponential growth we use a scale separation argument similar to that of section 4. 6.1.1. Scale separation for the finite dimensional approximation. When a finite dimensional space is used for the Smagorinsky-Navier-Stokes' model we can introduce a scale separation argument similar to (3.2), but here the nonlinear feedback in the second equation of (3.2) takes place through the computational error. As the computational error grows, the exponential growth of the error is moderated. More precisely defineū and u ′ by (3.2) with η defined by η := u −û h . We also note that then the result of Lemma 3.2 holds with η redefined to be the approximation error. Below we will refer to these results assuming that η is redefined as above.
6.1.2. Error estimates. First we prove a Lemma that is needed to estimate the consistency error of the stabilization term to the right order. Lemma 6.5. For v h ∈ V 0h , and f ∈ L 3 (I) there holds, for all v ∈ V and for all ǫ > 0, Then using an inverse inequality (5.2) followed by Morrey's inequality [19,Theorem B.42] and Lemma 2.1 we see that for all Using also the stability of the L 2 -projection on quasi-uniform meshes we have
One application of the previous lemma is the following consistency bound for the stabilizing term Lemma 6.6. The stabilizing term satisfies the bound Proof. By the definition of the stabilization term we have the bound Then by Lemma 6.5, s admits the bound Using (5.8), and |u| H 2 (Ω) u W 2,3 (Ω) we have The no-slip condition satisfied by u, the trace inequality (5.3) and the approximation (5.5) may be applied to control the second term, ) . The claim follows by collecting the above bounds. Theorem 6.7. Let u ∈ V and that u in addition has sufficient regularity so that C(u) < ∞ if we define (assuming U ≤ u L 3 (I;W 2,3 (Ω)) ), Assume that the parameter δ = O(h) and µ ≤ U h, and that the hypothesis of Theorem 6.2 hold. Then there holds sup The coefficient τ L in the exponential is defined by Proof. The proof follows the ideas of the result for the linear model problem Proposition 5.5, but the nonlinearity of the equation and the stabilization adds a layer of technicality. We proceed in 4 steps.
Step 1. Perturbation equation. Let η = u −û h , η Π = u − Πu and η π = u − πu. Similarly as in the perturbation argument for the Smagorinsky model we have, since ∇ ·û h = 0 andû h · n| ∂Ω = 0, Using the definition of the finite element method (6.1) we have the consistency relation Integrating in time over the interval I yields Step 2. Continuity and approximation properties of the right hand side of the perturbation equation. We now bound the terms in the right hand side. There are two different main difficulties involving nonlinearities. The material derivative and the viscous term on the one hand and the nonlinear stabilization on the other. 2.1. Material derivative and viscous term. The following term will be bounded in this paragraph.
First for the time derivative observe that using the orthogonality of the approximation error η Π , and for the linear viscous term For the convective part we write For the first term I we use similar arguments as those for Lemma 3.2 to get the bound In particular note that By Lemma 5.1 we may write the last term in the right hand side Since δ 2 h − 3 2 = h 1 2 and |u| H 2 (Ω) ≤ |u| W 2,3 (Ω) we have using approximation that, For the second term of the right hand side of (6.10) we first use that there exists E according to (5.9)-(5.11) to write using partial integration For the first term II a we use Cauchy-Schwarz inequality followed by the arithmetic-geometric inequality Applying now Lemma 6.5 to the right hand side we have for ǫ > 0, +h 3 )(U 3 + u 3 L 3 (I;W 2,3 (Ω)) ).
Here we used the bounds ). For the second term II b we use the fact thatû h is piecewise affine, add and subtract u and then integrate by parts H 2 (Ω)) + hτ L ∇ × u 2 L 2 (I;L ∞ (Ω)) |u| 2 L ∞ (I;H 2 (Ω)) ) Where we used that and . For II c we add and subtract u and π h η, Observing that η − π h η = u − π h u it follows, using Hölder's inequality, Sobolev injection and approximation, that . Once again adding and subtracting π h η we have . Using an inverse inequality (5.2) and the stability of the projection π h we see that The term marked ( * ) needs to scale as O(h 3 ). Since E 3 Ω = O(h 9 ) we see that we need h − 1 2 δ −4 ≤ h −6 , which is satisfied for δ = O(h). We arrive at the following bound for II c , where we used the bound h 9/2 |u| 3 W 2,3 (Ω) . Collecting the above bounds for the terms I and II and using that µ ≤ U h, we see that

Smagorinsky nonlinear viscosity.
We proceed to bound terms due to the Smaginsky nonlinear viscosity, i.e. the fifth and sixth terms on the right hand side of (6.9): Considering first the fourth term of the right hand side of (6.9) we have using (2.4) followed by (2.5) with p = 3, q = 3 2 . Then we use that ab L for all ǫ 0 , ǫ 1 ∈ R + . First observe that using Sobolev injection, and by applying Lemma 5.1, Collecting these bounds and choosing ǫ 0 and ǫ 1 small so that Now we consider the fifth term in the right hand side of (6.9), that quantifies the consistency error. We note that by partial integration and Cauchy-Schwarz inequality followed by Young's inequality, (2.5) and the stability of the L 2 -projection Π, we have where we used, that by the Cauchy-Schwarz inequality and (5.5) there holds and by Hölders inequality with p = 3/2 and q = 3 inequality followed by a global trace inequality and the inverse inequality (5.2) on ∂Ω, This completes the bound of the Smagorinsky terms. Collecting the different contributions we get 2.3. Terms related to boundary conditions and stabilization. These are the last three terms of the right hand side of (6.9).
The first term, related to approximation of initial data is bounded using approximation For the stabilization term finally we proceed using the Cauchy-Schwarz inequality followed by the arithmetic geometric inequality Applying Lemma 6.6 we conclude that Step 3. Application of scale separation argument. Collecting the above bounds (6.13), (6.15) and (6.16) and applying them to (6.9) we get, for γ 1 sufficiently large A naive application of Gronwall's lemma at this stage results in exponential growth with a coefficient proportional to the maximum of the velocity gradient. Instead recall the inequality proven in Lemma 3.2, Applying this inequality, with ǫ = 1/16 in the right hand side of (6.17) we obtain (6.18) + u 3 L 3 (I;W 2,3 (Ω)) + u 4 L 4 (I;W 2,3 (Ω)) + h ∇ · |∇u| F ∇u 2 Q .
Assuming U bounded by u L 3 (I;W 2,3 (Ω)) , we can absorb the first two terms in the right hand side in the u 3 L 3 (I;W 2,3 (Ω)) term.

Numerical examples
In the numerical examples below we use the minimal compatible element of [14] with piecewise affine velocity and piecewise constant pressure on macro elements (see [11] for implementation details and application to linear incompressible problems), together with the backwards differentiation formula BDF2 for time stepping. On each time step we use a linearized formula, using the solution from the previous time step, and solve only once. We use only Smagorinsky stabilization, i.e., we set γ 0 = γ 1 = 0 in (6.1). The Smagorinsky term is set by ν := γ |T | |∇u| F where |T | is the element area and γ controls the amount of dissipation in the model. 7.1. Shear layers. We shall first consider an example using the Euler equations, i.e., µ = 0 in (2.1), the double shear layer problem [7]. The computational domain is Ω) = (0, 2π) × (0, 2π) and initial conditions (7.1) u y (x, 0) = δ sin (x), u x (x, 0) = tanh ((y − π/2)/ρ), y ≤ π tanh ((3π/2 − y)/ρ), y > π which gives two horizontal shear layers perturbed by a small vertical velocity. We take ρ = π/15 and δ = 0.05 and apply periodic boundary conditions. The problem is solved on a Union Jack (macro element) mesh of 100 × 100 (boundary) nodes using a timestep size k = 1/100. In Figs. 1-3 we show the effect of dissipation on vorticity by increasing γ from zero to 10 −1 . Th unstabilized solution shows oscillations in the solution which worsen with time, whereas the choice γ = 10 −1 compares well with [7]. 7.2. Vortex shedding. The second example shows the effect of dissipation on von Karman vortex shedding around a cylider. We used µ = 3 × 10 −4 (which is close to the limit for stable solutions with γ = 0 on this mesh). The outer domain is (−1/2, 2) × (−1/, 1/2), and the cylinder has center at the origin and radius r = 1/10. The boundary conditions are u = 0 at y = ±1/2 and at the cylinder, homogeneous Neumann conditions at x = 2, and u = (3/2 − 6y 2 , 0) at x = −1/2. The initial conditions correspond to the stationary Stokes solution, and the timestep size k = 1/100. The computational (macro) mesh is shown in Fig. 4, and the solution at time t = 10 is shown in Fig. 5 for increasing γ. We note that too much dissipation severely affects the shape of the vortex street, whereas smaller amounts of dissipation only serve to stabilize the solution.
In Fig. 6 we show that the method can be unstable for high Reynolds numbers. We used µ = 10 −6 and show the instability evolving for t = 0.15, t = 0.2 and t = 0.3 with γ = 0. In Fig. 7 we show the corresponding stabilized solution with γ = 10 −1 . The velocities are shown in the nodes of the macro mesh. Finally, in Fig. 8 we show the relative streamlines at time t = 10 for the stabilized model.

Conclusion
We have considered a p-Laplacian Smagorinsky model for high Reynolds flow problems. We showed using a scale separation argument that the Smagorinsky model has stability properties that only depend on the large scales of the flow, provided a spectral gap exists for the particular flow configuration. The set of non-essential fine scales grows as the perturbation due to the Smagorinsky model increases, hence moderating the exponential growth. In a second part we considered the Smagorinsky model as a stabilizing term in a low order finite element method and we showed that the resulting stabilized method has optimal properties for smooth solutions in the laminar high Reynolds number regime. The stabilized finite element method also inherits the reduced exponential growth through scale separation from the continuous case.
To the best of our knowledge these results are the first that give quantitative evidence that the Smagorinsky model enhances both the accuracy and the stability for computations of high Reynolds number flows.