Invariant measures for the 3D Navier-Stokes-Voigt equations and their Navier-Stokes limit

The Navier-Stokes-Voigt model of viscoelastic incompressible fluid has been recently proposed as a regularization of the three-dimensional Navier-Stokes equations for the purpose of direct numerical simulations. Besides the kinematic viscosity parameter, $\nu>0$, this model possesses a regularizing parameter, $\alpha>0$, a given length scale parameter, so that $\frac{\alpha^2}{\nu}$ is the relaxation time of the viscoelastic fluid. In this work, we derive several statistical properties of the invariant measures associated with the solutions of the three-dimensional Navier-Stokes-Voigt equations. Moreover, we prove that, for fixed viscosity, $\nu>0$, as the regularizing parameter $\alpha$ tends to zero, there exists a subsequence of probability invariant measures converging, in a suitable sense, to a strong stationary statistical solution of the three-dimensional Navier-Stokes equations, which is a regularized version of the notion of stationary statistical solutions - a generalization of the concept of invariant measure introduced and investigated by Foias. This fact supports earlier numerical observations, and provides an additional evidence that, for small values of the regularization parameter $\alpha$, the Navier-Stokes-Voigt model can indeed be considered as a model to study the statistical properties of the three-dimensional Navier-Stokes equations and turbulent flows via direct numerical simulations.


Introduction
In this work, we consider the three-dimensional Navier-Stokes-Voigt equations (NSV), a viscoelastic incompressible fluid model introduced by Oskolkov in [30], and which was proposed by the authors of [2] as a smooth inviscid regularization of the 3D Navier-Stokes equations (NSE) for the purpose of direct numerical simulations (DNS). More specifically, we consider the NSV model subject to periodic or no-slip boundary conditions, and driven by a given force field f. The velocity vector field, v(x, t), and the scalar kinematic pressure, p(x, t), are governed by the system of equations v(x, t) = 0 x ∈ ∂Ω, or v(x, t) is periodic; (1) in the smooth domain Ω ⊂ R 3 , in the case of no-slip Dirichlet boundary condition, or with basic periodic domain Ω = [0, L] 3 ⊂ R 3 , when equipped with periodic boundary conditions. Here, α ≥ 0 is a given length scale parameter, and ν > 0 is a given kinematic viscosity, such that α 2 /ν is the relaxation time of the viscoelastic fluid.
In the periodic case, we assume that the driving force, and the initial velocity u 0 satisfy and it is easy to see that it implies Ω v(x, t) dx = 0, ∀t ≥ 0.
The Navier-Stokes-Voigt model of viscoelastic incompressible fluid, (1), (sometimes written as Navier-Stokes-Voight) was introduced by Oskolkov in [30], and pointed out by O. Ladyzhenskaya as one of the reasonable modifications of the Navier-Stokes equations, see [26]. In [30], A. P. Oskolkov studied and proved its solvability in different functional spaces. We remark that these equations behave like a damped hyperbolic (pseudo-parabolic) system, see [22], and, therefore, the solutions do not experience fast (instantaneous) smoothening of the initial data, as it is for parabolic systems like the Navier-Stokes equations. This fact could prevent the NSV model from being a reasonable modification of the Navier-Stokes equations, however, since we are proposing it as a model for direct numerical simulations of turbulent flows in statistical equilibrium, i.e., after the solutions reach the global attractor, we are mainly interested in its long time behavior. Indeed, it was proved in [21], that solutions in the global attractor are smooth, if the forcing field is smooth enough, even for initial data satisfying only finite kinetic energy and finite enstrophy (i.e. bounded in the Sobolev H 1 -norm). In particular, in [21], it is shown in the periodic case, that if the forcing field is analytic, then the global attractor consists of analytic functions. This result, in conjunction with results proved in this work, proves that if the forcing field is smooth enough, then averaged structure functions, with respect to an invariant measure for the NSV flow, display exponential decaying tail.
As it was observed above, the NSV model presents an extra length scale associated with the viscoelastic properties of the fluid, the parameter α, besides the well known Kolmogorov length scale, η, (see, e.g., [18] for its definition), which is usually associated with the smallest scales of motion in turbulent flows. In [25], based on numerical simulations of the Sabra Shell model, it is observed that for a large range of values of α > 0, i.e. α ≫ η, there are two distinct regions associated with the inertial range of the energy spectrum for the NSV model. The first one obeying the celebrated Kolmogorov k −2/3 power law (with anomalous correction), followed by a second range of wavenumbers, where energy condensates, and it is simply equipartitioned. The range of the second power-law, however, slowly disappears as α/η decreases, restoring the usual Navier-Stokes inertial range regime, when α ≪ η.
The numerical simulations of the Shell model presented in [25] also suggest that for a large range of values of α > 0, when α ≫ η, small scales velocity fluctuations are vigorously damped, due to a slowness of the energy transfer timescales associated with these small length scales. This has an extraordinary effect in reducing the energy dissipation rate intermittency, characterized by violent fluctuations away from the average in the energy dissipation rate. It is therefore suggested in [25] that by tuning the parameter α, one may attenuate the strong velocity fluctuations related to the intermittent events, reducing thus the stiffness of DNS of turbulent flows, with only a small effect on the energy containing scales. This reduction of stiffness in the NSV model has also been observed in numerical experiments in implementing the NSV model for inpainting [8].
The numerical experiments in [25] were performed using the Sabra shell phenomenological model of turbulent flows. However, because these models do not present any spatial structure, it is not clear if such properties persist for full direct numerical simulations of 3D turbulent flows. In fact, a large class of viscoelastic flows present anomalous behavior even for very small viscoelastic parameters, see, e.g., [28], and such phenomenon could be present in flows governed by these equations.
In this work, we rigorously establish several statistical properties of the energy spectrum of the 3D NSV model that were observed in the shell model simulations in [25], using the notion of invariant measures. For example, by considering the results of [21], we can rigorously justify the exponentially decaying dissipation range observed in the simulations presented in [25], despite the fact that these equations behave like a damped hyperbolic (pseudo-parabolic) system [22].
Concerning our main result, we start by defining a notion of strong stationary statistical solution of the Navier-Stokes equations, which is a reg-ularized version of the concept of stationary statistical solution that was introduced and investigated by Foias in [10], [11] (see also [15]). This notion is a generalization of the concept of an invariant measure for the semigroup generated by the three-dimensional Navier-Stokes equations.
Our main result states the following: Fix ν > 0, and let µ αn be a given sequence of invariant measures for the semigroups generated by the NSV model. Then, there exists a subsequence, µ αn j , weakly converging to a strong stationary statistical solution of the Navier-Stokes equations, as α n j → 0.
This fact supports the observations reported in [25], providing a rigorous evidence that, for small values of the regularization parameter α, the NSV model approximates several statistical properties of the 3D Navier-Stokes equations, and, therefore, may be used as a reliable subgrid scale model for direct numerical simulations of turbulent flows.

Mathematical framework
We will follow the standard functional formulation for the Navier-Stokes equations, see, e.g., [3,15,35], and the functional formulation for the NSV equations used in [22]. We denote by V, in the no-slip case, the set In the periodic case, we define V as V = {All the three-dimensional vector valued trigonometric polynomials ϕ with basic periodic domain Ω = [0, L] 3 ; ∇ · ϕ = 0 and Ω ϕ dx = 0 .
(3) The two fundamental functional spaces in this work are defined by The inner products in H and V are denoted, respectively, by and the associated norms by |u| = (u, u) 1/2 , u = ((u, u)) 1/2 (the latter is a norm thanks to the Poincaré inequality (8), below). We denote by P LH the (Leray-Helmholtz) orthogonal projector in L 2 (Ω) 3 onto the subspace H. The Stokes operator is defined by Moreover, because the inverse of the Stokes operator is a positive self-adjoint compact operator in H, we can define its powers A s , s ∈ R, with domain D(A s ). We have V = D(A 1/2 ) and its dual V ′ = D(A −1/2 ), see, e.g., [3], [15], or [35].
The term B(u, v) = P LH ((u · ∇)v) is a bilinear term associated with the inertial term, which satisfies where V ′ denotes the dual space of V , (see, e.g., [3,35]). Thus, taking the duality action in V × V ′ , between the bilinear term, B(u, v) and a vector field w ∈ V yields a trilinear term (see, e.g., [3,35] which is defined for u, v, w ∈ V , where ·, · V,V ′ denotes the duality action between the spaces V , and its dual space V ′ . An important relation for the trilinear term is the orthogonality property, see, e.g., [3,35], for u, v ∈ V . It follows from this relation the anti-symmetry property (see, e.g., [3,35] for all u, v, w ∈ V . Moreover, because the inverse of the Stokes operator is a positive selfadjoint compact operator in H, then there exists a complete orthonormal basis of H formed by eigenvectors, {w j } j , with associated eigenvalues, Aw j = λ j w j , satisfying 0 < λ 1 ≤ λ 2 ≤ . . . ≤ λ j → ∞, as j → ∞, see, e.g., [3,15,35] for details. In this geometry and setting, the Poincaré inequality holds where λ 1 is the first eigenvalue of the Stokes operator. We will need the following inequalities, which are consequences of Hölder inequality, Poincaré inequality, and well-known inequalities of Ladyzhenskaya, see, e.g., [3], [15], or [35] : for every u ∈ H, and v, w ∈ D(A). Also, for every u, v ∈ D(A), and w ∈ H. Throughout this work, if not otherwise stated, we assume a forcing field f ∈ H. If we apply the Leray-Helmholtz operator, P LH , to the first equation of the system (1), we obtain the following equivalent functional differential equation Global existence and uniqueness of (11) was first studied in [30], where it was established that the system (1) generates a continuous semigroup, S α (t) : V → V , t ∈ R + , i.e., a one-parameter family of maps {S α (t)} t≥0 , satisfying the properties of a continuous semigroup, such that the solution of (1) (or equivalently, of (11)), satisfies v(t) = S(t)v 0 , for every v 0 ∈ V . In [2], global regularity of (11) was also proved for the inviscid model, i.e. when ν = 0. Moreover, higher-order regularity of the inviscid case, i.e. the Euler-Voigt model, is established in [24]; furthermore, a new blowup criterion for the three-dimensional Euler equations, by means of this inviscid regularization is introduced.

Energy budget for the NSV model
The global regularity result for solutions of the Navier-Stokes-Voigt equations established in [30] implies that the following energy equality holds for every t ∈ [0, ∞): (in fact the energy equality and global well-posedness holds for all t ∈ R). Similarly, the global regularity results established in [2] imply that the solutions of the NSV equations in the inviscid (i.e. Euler-Voigt model, ν = 0) and unforced setting, f = 0, satisfy for every t ∈ R: Therefore, the conserved quantity in the inviscid and unforced setting of the NSV (i.e. Euler-Voigt) model is which we call the α-energy. The quantity is the usual kinetic energy. We remark that while the kinetic energy is formally conserved for the inviscid and unforced incompressible Navier-Stokes equations, this is not the case for the inviscid and unforced setting of the NSV model. Those differences in the conserved quantity for the inviscid and unforced case are the main source of deviation between the statististical properties of the NSV model and of the Navier-Stokes equations observed in [25]. In Section 6, we establish a result concerning the convergence of the averaged kinetic energy of the NSV model as the parameter α tends to zero. Let w j be the orthonormal basis of H composed by the eigenfunctions of the Stokes operator, A. For a vector field v ∈ H, we define the component v κ , for a wavenumber κ, by We also define the Then, we can write the projected Navier-Stokes-Voigt equations on the shell Let us now obtain the α-energy budget for the shell [κ ′ , κ ′′ ). Multiplying in L 2 by v κ ′ ,κ ′′ , we obtain where is the net rate of α-energy transfer at κ, and represents the net rate of α-energy from the lower modes to the higher modes, and , v κ 1 ,κ ) represents the net rate of α-energy from the higher modes to the lower modes. Remark. Since we are dealing simultaneously with α > 0 and α = 0, then for κ ′′ = ∞, all the calculations above are formal. However, they can all be rigorously recovered for the NSV case, i.e. when α > 0, via the Galerkin approximation procedure. For the Navier-Stokes case, i.e. when α = 0, if κ ′′ = ∞, we can only state equations (18) and (19) in terms of inequalities, due to a possible lack of regularity for the three-dimensional Navier-Stokes equations, see, e.g., [15].

Invariant measures and stationary statistical solutions for the NSV model
A rigorous mathematical framework for investigating the statistical properties of turbulent fluid flows in statistical equilibrium was first considered by Hopf in [19]. In [10,11], Foias established the notion of stationary statistical solutions of the Navier-Stokes equations, which is a generalization of the concept of invariant probability measures for the semigroup associated with the solutions of the equations of motion. This notion is important because the semigroup generated by the Navier-Stokes equations is not known to be well-posed in V , and, therefore, the notion of invariant measure in V is not well-posed either. A detailed discussion of this issue can be found in [15].
Because of the global regularity of solutions of the NSV model in the space V obtained in [30], one can refrain from using the abstract notion of stationary statistical solutions, and work only with the more familiar notion of invariant measure.
However, because the main goal of this work is to show convergence results of the statistical properties of the NSV model to the corresponding ones of the Navier-Stokes equations, as α → 0, we will show first that invariant measures of the NSV model are stationary statistical solutions of the NSV model, and using this fact, we prove a weak convergence theorem in the sense of stationary statistical solutions of the Navier-Stokes equations. This approximation result is the main reason why we are interested in deriving several statistical properties of the invariant measures associated with the NSV model in this section.
We recall that a probability measure µ α on V is called an invariant measure for the semigroup It is easy to see that the Dirac measure concentrated at the steady state solutions of the NSV model, which coincide with steady solutions of the Navier-Stokes equations, are invariant measures for the semigroup {S α (t)} t≥0 generated by the NSV equations. Therefore, because the set of steady states is nonempty (see, e.g., [3], [16] and [35]), we have that the set of invariant measures for {S α (t)} t≥0 is nonempty.
We recall that the support of an invariant measure, µ α , consists only of points that are nonwandering, i.e., for any u ∈ supp µ α , and any Borel set E containing u, there exists a sequence of positive times t n → ∞, such that S α (t n )E ∩ E = ∅, for all n ∈ Z + , see, e.g., [15]. This implies that the support of an invariant measure, µ α , is included in the global attractor of the semigroup S α (t).
In [22], it was proved that S α (t) : V → V has an absorbing ball in V , and it is an asymptotically compact semigroup, implying the existence of a global attractor in V . Moreover, the global compact attractor was shown to be bounded in D(A) and to have finite Hausdorff and fractal dimensions.
Furthermore, despite the fact that the NSV equations behave like a damped hyperbolic (pseudo-parabolic) system, rather than a parabolic equation, it was proved in [21] that the 3D periodic NSV equations possess an asymptotic smoothing property. More specifically, let us define the Gevrey class of functions whereû j are the correspondent Fourier coefficients of the Ω-periodic vector field u, τ > 0, and r ≥ 0. The space is equipped with the corresponding inner product (21) for u, v ∈ G r τ , and with the corresponding norm for u ∈ G r τ . One can prove that the space of real analytic functions C ω (Ω) has the following characterization: for any r ≥ 0, see, e.g., [5], [27] for details. In [21], it was proved that if the driving force, f, belongs to a Gevrey class of functions, G 1 τ 0 , for some τ 0 > 0, then for every solution v(x, t) of the 3D periodic NSV equations, with initial data v 0 ∈ V , there exists t 0 > 0, and a function w(t) ∈ L ∞ (0, ∞, G 2 τ ), for some τ > 0 depending only on |f | 1,τ 0 , ν, λ 1 and α, so that This result proves, in particular, that if f is analytic, then the global attractor consists of analytic functions. We recall that this technique of Gevrey class regularity was first introduced by Foias and Temam [17] for proving the analyticity, in space and time, for strong solutions of the three-dimensional Navier-Stokes equations, for short time (see also [9] for generalization of the technique to parabolic analytic equations).
The following theorem is a straightforward application of the facts discussed above: Theorem 1 Let µ α be an invariant measure for the semigroup {S α (t)} t≥0 generated by the 3D periodic NSV equations, (1). If the forcing field f belongs to H, then the supp µ α is a subset of the global attractor, A, of the semigroup generated by the 3D NSV model, Remark. We remark that the trigonometric polynomials, usually considered as forcing fields in the numerical investigations of turbulent flows, belong to the class of analytic functions. Therefore, our result applies to a wide set of simulations. In particular, it justifies the observed exponential tail in Sabra shell model simulations of the NSV equations in [25] (see similar results concerning the three-dimensional Navier-Stokes equations [6]). Now, we define the notion of stationary statistical solutions of the Navier-Stokes-Voigt equations inspired by Foias [10,11], and [15].
for any test functional Ψ ∈ T α , where T α is defined in Definition 2. (3) Definition 2 We define the class T α of test functions to be the set of all real-valued functionals Ψ = Ψ(u) on H that are bounded on bounded subsets of V and such that the following conditions hold: , and the map u → Ψ ′ (u) is continuous and bounded as a function from H into H.
For example, we can take the cylindrical test functions Ψ : H → R of the form Ψ(u) = ψ ((u, g 1 ), . . . , (u, g m )), where ψ is a C 1 (R m ; R) scalar function on R m , m ∈ N, with compact support, and g 1 , . . . , g m belong to H. For this case we have where ∂ j ψ denotes the derivative of ψ with respect to the j-th variable. In this case, it follows that Ψ ′ (u) ∈ H since it is a linear combination of the g j .
The class of test functions T α are broader than the class considered in the usual definition of stationary statistical solutions of the Navier-Stokes equations, see, e.g., [15]. This is made possible by the fact that solutions of the NSV model (1) are globally regular in V , see [2,30].
Condition 1 in Definition 1 implies that the support of the stationary statistical solutions are included in D(A). Condition 3 is a local energy balance equation, which implies that the stationary statistical solutions have supports which are bounded in V , see Proposition 1 below. Moreover, in Corollary 1, we show that for fixed α 0 > 0, Condition 3 implies that the supports are uniformly bounded H, for all α ∈ (0, α 0 ].
We also remark that because we consider weakly converging subsequences of probability measures in the Section 6, we work with the weak topology of V . However, because V is a separable Hilbert space, the Borel σ-algebra generated by the weakly open sets coincides with that for the open sets in the strong topology of V , see, e.g., [7]. Thus, we identify these two probability spaces in the rest of the work. Now, we prove a result concerning bounds for the supp µ α in the space V . This result will be important because we will derive from it a uniform bound, with respect to α, in the H-norm.
Proposition 1 Let µ α be a stationary statistical solution of the Navier-Stokes-Voigt equations. Then, Proof. We follow the arguments used in [15] for the NSE case. It follows from Definition 1 item (3), that if then, by the Cauchy-Schwarz and Poincaré, (8), inequalities, where λ 1 is the first eigenvalue of the operator A. Thanks to Poncaré inequality, (8), and Condition 1 in Definition 1, the term in the right-hand side of the last inequality is bounded, and therefore, Notice that if we choose, in the set Γ, E 1 = 0, and let E 2 → ∞, and using the fact that µ α is a probability measure, we find that Let 0 < ǫ < 1, we define By Poincaré inequality, (8), we have that |u| 2 + α 2 u 2 ≤ (λ −1 1 + α 2 ) u 2 , and, therefore, We obtain, by applying (26) to the set Γ ǫ , that µ α (Γ ǫ ) = 0. Letting ǫ ↓ 0, we conclude (24). Now, we define a set that will be often used throughout the work The following corollary results immediately from Proposition 1.
Corollary 1 Let µ α be a stationary statistical solution of the Navier-Stokes-Voigt equations, then Next, we show that for every α > 0, invariant measures for the semigroup S α (t) : V → V generated by the NSV model, (1) (or (11)), are stationary statistical solutions for the NSV model. The proof is in the same lines of a related result for the two-dimensional Navier-Stokes presented in [15]. (11)). Then, µ α is a stationary statistical solution of the NSV model, in the sense of Definition 1.
Proof. Let µ α be an invariant measure for the semigroup S α (t) : V → V generated by the NSV model. Then, by Theorem 1, supp µ α is included in the global attractor, A, of S α , which is bounded in D(A). Therefore, the function u → u 2 D(A) is bounded on the supp µ α , which implies that Therefore, Condition 1 in Definition 1 is satisfied. Now, we prove Condition First, we observe that thanks to the Cauchy-Schwarz inequality, the Poincaré inequality, (8), and to (30), for all t ≥ 0. Now, take the average with respect to t over [0, T ], and use the fact that the left-hand side of the above equation is constant in time, to obtain where in the last step we are allowed to use the Fubini theorem, (see, e.g., [15]), because the integrand is bounded in Γ, jointly continuous on [0, T ] × V , and belongs to L 1 ([0, T ] × V ), for all T > 0. By the energy equality (12), we obtain , then by letting T → ∞ we obtain by the Lebesgue dominated convergence theorem that which is Condition 3 in Definition 1. For every m ∈ Z + , let us define the projector P m : H → H by where {w j } is an orthonormal basis of H consisting of eigenvectors of A. We now prove Condition 2 of Definition 1. Let Φ ∈ T α be a test functional.
Now, using the fact that the left-hand side is constant in time, we take the average with respect to t over [0, T ], to obtain where we have used Fubini's theorem in the last step. Now, if u ∈ D(A), then by the global regularity result of the NSV proved in [2,30], see also [22], the function t → (I + α 2 A)S α (t)u belongs to C 1 ((0, ∞), H). Therefore, because the operator (I + α 2 A) −1 : H → H is bounded, see Lemma 2, we obtain, by (11), that S α (t)(u) satisfies Therefore, for all t > 0, and any u ∈ D(A), where F is as in (11). Then, take the average with respect to t over [0, T ], and substituting the result in the last line of (34), we conclude that (37) Because Φ m is bounded, we obtain, as T → ∞, Now, because P m u → u, in H, as m → ∞, for all u ∈ V , we have (40) Therefore, because Φ ′ is bounded, we have by (9) that there exists a constant c > 0, which is independent of α, such that where the right-hand side belongs to L 1 (dµ α ), by (30). Then, by the Lebesgue dominated convergence theorem, we obtain Therefore, Condition 2 holds, completing the proof of the theorem.

Corollary 2
The set of stationary statistical solutions of the NSV equations is nonempty.
Proof. As we have observed earlier, the set of steady state solutions of the NSV, or equivalently of the NSE, is non-empty, see, e.g., [3,35]. Therefore, the Dirac delta measures concentrated on the steady state solutions and all their convex combinations are invariant measures. Consequently, by Proposition 2, they are stationary statistical solutions of the NSV model.
Remark. One could also prove the converse of Proposition 2, i.e., that every stationary statistical solution of the NSV model is an invariant measure, by following the arguments of Foias in [10,11] for the 2D periodic case, see also [15]. Only a slight modification in the treatment of the nonlinear term is necessary, but which can be circumvented by using that supp µ α is included in D(A), and is bounded in V . We do not present it here, since the semigroup generated by the NSV model is well defined in V , and therefore, we can use the more natural notion of invariant measures in our approximation results in Section 6.
In the sequel, we use the symbol . . . α to denote average with respect to an invariant measure of the NSV model µ α , i.e.
We can use the stationary statistical solutions formalism to show that the stationary Reynolds averaged equations hold for the NSV equations in this framework. Indeed, we can define the mean velocity field u α ∈ D(A), thanks to Condition 1 in Definition (1), as follows Moreover, thanks to (10) and Condition 1 of Definition 3, we may define the average B(u, u) α ∈ H by The mean flow u α is a vector field on Ω with u α ∈ D(A), while B(u, u) α ∈ H. The next proposition shows that, assuming statistical equilibrium, the vector field u α also satisfies the Reynolds averaged equations in H. This should be compared to the NSE case, where these equations are known to be valid only in V ′ , (see, e.g., [15]).
Proposition 3 Let µ α be an invariant measure of the NSV equations. Then, the following functional form of the Reynolds averaged equations holds in H: Proof. We will use here the notation of Section 3, and set κ 2 1 = λ 1 . Let ψ be a C 1 real-valued function with compact support on R. For any v ∈ H, any finite wavenumber κ ≥ κ 1 , and every u ∈ D(A), the real-valued function is a cylindrical test function, and it satisfies for every φ ∈ D(A). Thus, considering φ = (I + α 2 A) −1 F (u), which belongs to D(A), we obtain By Proposition 2, µ α is also a stationary statistical solutions of the NSV model, therefore, by substituting the expression above in Condition 2 of Definition 1, we obtain Now, by virtue of Proposition 1, supp µ α is bounded in V , and therefore the map u ∈ supp µ α → ((I + α 2 A) 1/2 u, (I + α 2 A) 1/2 v κ 1 ,κ ) is uniformly bounded in supp µ α , let us say by M > 0. Therefore, we can choose ψ, above, such that ψ ′ ≡ 1 at the interval [−M, M], therefore, ψ ′ (((I + α 2 A) 1/2 u, (I + α 2 A) 1/2 v κ 1 ,κ )) = 1, for every u ∈ supp µ α . This yields For each fixed v ∈ H, we observe that since supp µ α ⊂ D(A), we have by the Cauchy-Schwarz inequality, and by (9), that Thanks to Condition 1 of Definition 1, we can apply the Lebesgue dominated convergence Theorem, and let k → ∞ to obtain which yields the result.

Averaged Energy Budget for the NSV model
In this section, we follow [15] to investigate the energy distribution scale-byscale for the 3D NSV equations. For the Navier-Stokes equations, assuming that there exists an inertial range, i.e. an extensive range of wavenumbers where the viscous dissipation does not play a significant role, one can show that the energy simply cascades through these wavenumbers with a rate equal to the mean energy dissipation rate, ǫ = ν u 2 α . For the NSV equations, a similar scenario holds for the α-energy As it is usual in the studies of homogeneous turbulence, and using the notation of Section 3, we will consider the forcing f with finite eigenmodes, i.e. f = for the orthonormal basis, {w j }, of H consisting of the eigenvectors of A, with κ ≥ κ 1 = √ λ 1 . Now, we argue as in the proof of Proposition 2 to obtain an averaged energy balance equation. We skip the details, but we sketch the derivation: Using the fact that µ α is an invariant measure, using equation (18), averaging with respect to t over [0, T ], and using Fubini's Theorem, we have where Γ is as in the proof of Proposition 2. Now, let κ ′′ ≥ κ ′ >κ, and letting T → ∞ in the above equation, we obtain the following balance equation The expression on the right-hand side of the last equality is the mean net α-energy transfer in the energy shell [κ ′ , κ ′′ ]. Because the supp µ α in included in D(A) and is bounded in V , we can use the Lebesgue dominated convergence Theorem, and let κ ′′ → ∞ to obtain This expression shows that the net α-energy transfer is positive for every κ >κ. We remark that a similar expression holds for the three-dimensional NSE case substituting equality by inequality, see, e.g., [15]. Therefore, we have proved Proposition 4 Let µ α be an invariant measure for the semigroup, {S α (t)} t≥0 , generated by the NSV model, and letκ < κ ′ ≤ κ ′′ . Then Moreover, for all κ >κ, we have The result above shows that the mean net α-energy transfer is positive, which is an important consistency check for turbulence models. Moreover, if we assume that there exists a range of wavenumbers, [κ ′ , κ ′′ ], where the viscous dissipation term satisfies ν u κ ′ ,κ ′′ 2 α ≪ 1, it follows from (48) that the α-energy transfer is nearly constant within this range, that is, there exists an inertial range for the NSV model.

Statistical approximations of the 3D Navier-Stokes equations
In this section we prove that given a sequence of invariant measures for the NSV equations, {µ αn }, with the regularizing parameter α n converging to 0, as n → ∞, there exists a subsequence converging weakly to a strong stationary statistical solutions of the 3D Navier-Stokes equations, a notion that we introduce in Definition 3 below. This is a slight modification of the notion of stationary statistical solutions of the 3D NSE introduced by Foias in [10,11], see also [15]. The main difficulty in proving this result is that the supports of the invariant measures, µ α , of the NSV are not uniformly bounded, with respect to the parameter α, in the norm of V , as α tends to zero. This means that we cannot use in a straightforward manner the standard Prokhorov's Theorem, see, e.g. [34,Thm. 1.12]. This is because the whole space V is not metrizable when endowed with its weak topology. Therefore, we have to consider more general theorems concerning convergence of measures in arbitrary Hausdorff spaces. Fortunately, the family of measures considered in this section satisfy the tightness condition of [36], and the weak convergence holds for our case.
We start by defining the terminology, from [36], that is relevant to the results presented in this section. A paving is a non-empty set consisting of subsets of a given set X. N ⊂ 2 X is said to be a (∅, ∪f, ∩c)−paving of X if N is closed under finite unions and countable intersections, and if in addition ∅ ∈ N . Similarly, N ⊂ 2 X is said to be a (∅, ∪f, ∩f )−paving of X if N is closed under finite unions and finite intersections, and ∅ ∈ N . A paving K is called compact [semicompact] if every family [every countable family] of sets in K, which has the finite intersection property (the intersection of any finite subset of K is non-empty) has a non-empty intersection (the intersection of all elements of K is non-empty). A paving L separates the sets in K if to any pair K 1 , K 2 of disjoint elements in K, we can find a pair G 1 , G 2 of disjoint sets in L such that K 1 ⊂ G 1 and K 2 ⊂ G 2 .
The results stated in [36] consider a Hausdorff topological space X, and a net, {x β }, in X. The notion of net can be found in [36], but it is not important for our present work because the space V is a separable topological space with its weak topology, and thus all statements can be done in terms of sequences. However, we will state a theorem appearing in [36] in its original formulation, in terms of nets. We only recall that a net (x β ) β∈D on a topological space X is compact if every subnet has a further subnet that converges.
Let X be a Hausdorff space, and let L and K be pavings in X. We denote by B = B(X, L), the smallest σ-field containing every set E ⊂ X for which L ∩ E ∈ L, ∀L ∈ L. We denote by M + (X, L) the set of finite, non-negative measures defined on B(X, L). M + (X, L, K) denotes the set of measures in M + (X, L) which are regular with respect to K, i.e., µ ∈ M + (X, L, K) if and only if µ satisfies From now on, we will consider the following axioms 1. X is a Hausdorff topological space.
5. L separates the sets in K.

K is semicompact.
Let µ ∈ M + (X, L), and let µ β be a net in M + (X, L). Then, we say that µ β converges to µ in the weak topology if and only if µ β (X) → µ(X) and lim inf µ β (G) ≥ µ(G), ∀G ∈ L. Now, we are going to state the abstract result that we need to continue with our investigation. It is just the first half of Corollary 1 in [36].
We will now translate all the preceding abstract definitions to our specific context. We consider the topological space X = V , where V is defined in (4), endowed with its weak topology. Thanks to the Hahn-Banach Theorem, the space V is a Hausdorff topological space with the weak topology. Moreover, V is also separable with the weak topology, see, e.g., [15]. This implies that there is no need to use nets in the application of Theorem 2 in this work, and we will consider only sequences in the subsequent results.
We for every h : V → R bounded continuous function. We also recall that the subset of Borel probability measures M 1 is closed in M + (V ) under the weak topology, (see, e.g., [7]), which implies that weakly converging sequences in M 1 converge to elements in M 1 . Now we will prove that every Borel finite measure defined on V endowed with the weak topology satisfies the tightness condition of Theorem 2.
Proof. First, we notice that because V with the strong topology is a separable Hilbert space, the Borel σ-algebra associated with the strong topology is the same as the Borel σ-algebra associated with the weak topology, (see, e.g., [15]).
Moreover, because V with the strong topology is separable, any strongly Borel probability measure, µ, is regular in the sense that for every strongly Borel set E, see, e.g., [7]. Now, because strongly compact sets in V are also weakly compact, we have that if E is a weakly Borel set in V (which also implies that it is a strongly Borel set), then µ satisfies Therefore, which shows that µ ∈ M + (V, L(V ), K(V )). Now, it is easy to see that Theorem 2, Lemma 1, and the preceding discussions imply the following theorem: Because of the convergence result of our main theorem in this section occurs in the space V , endowed with its weak topology, the limit measure thus obtained will be defined on this space as well. Therefore, we need to define the following notion Definition 3 A strong stationary statistical solution of the 3D Navier-Stokes equations is a Borel probability measure µ on V such that for any test functional Ψ ∈ T , where T is defined in Definition 4. (3) As remarked before, we may use the weakly Borel σ-algebra in Definition 3 because V is a separable Hilbert space, and therefore, weakly Borel sets coincide with the strong Borel sets, see, e.g., [7].

Definition 4
The class of test functions T is the set of functions Ψ : V → R of the form Ψ(u) := ψ ((u, g 1 ), . . . , (u, g m )) , where the function ψ is a C 1 scalar valued function defined on R m , with m ∈ Z + , and g 1 , . . . , g m belong to V .

Remark 1:
We remark again that for f ∈ H, there exist steady state solutions of the 3D Navier-Stokes equations which belongs to D(A), see, e.g., [3,35]. Therefore, every Dirac measure concentrated on a steady solution of the 3D Navier-Stokes equations which belongs to D(A) is a strong stationary statistical solution. Moreover, if u(t) is a strong solution of the NSE, bounded in V for every t ≥ 0, then the trajectory of u(t) is included in a closed ball in V , which is compact in the weak topology of V . Therefore, one can also generate a strong stationary statistical solution by the Krylov-Bogolyubov procedure for generating invariant measures; see, e.g., [4], [15] for details on how to use this procedure to generate time-average measures and the Banach limit. Remark 2: We remark that strong stationary statistical solutions are defined over the space V with its corresponding Borel σ-algebra, while that the usual stationary statistical solutions, defined by Foias in [10], [11], are defined over the space H with its Borel σ-algebra. This justifies the denomination strong stationary statistical solutions. Indeed, because Leray-Hopf weak solutions of the 3D NSE are uniformly bounded in the H-norm, with respect to time, one can generate, by the Krylov-Bogolyubov procedure, a Borel probability measure defined over H, which is a stationary statistical solution, (see, e.g., [15]). However, because strong stationary statistical solutions are defined over V , the Krylov-Bogolyubov procedure can generate a strong stationary statistical solution only if applied to strong solutions of the 3D NSE, that are uniformly bounded, with respect to time, in the V -norm, as described in the former remark.
We now proceed with the computation of Ψ ′ for test functions Ψ ∈ T . It is easy to see that for Ψ ∈ T , as defined in (56), we have g 1 ), . . . , (u, g m )) (φ, g j ).
It is easy to see that φ → ∇ u Ψ(u)·φ is a bounded linear continuous functional in H, hence, by the Riesz representation theorem, there exists an element Ψ ′ (u) ∈ H such that This is the identification that is used in Definition 3. More precisely, We now establish a lemma that will play a key role in the rest of this work.
Proof. Indeed, let {w j } be a complete orthonormal basis of H formed by the eigenvectors of the operator A, and let 0 < λ 1 ≤ λ 2 . . . ≤ λ j → ∞, as j → ∞, be the corresponding eigenvalues of the operator A associated with the eigenvectors {w j }. Then, by the Parseval identity, for every α ≥ 0, which proves that (I + α 2 A) −1 H→H ≤ 1. Now, for every φ ∈ H, we have by the Parseval identity that j |A η φ| 2 , for all j ≥ 1, we once again have by the Parseval identity that The asymptotic behavior of the eigenvalues of the Stokes operator defined in Ω ⊂ R n is known to satisfy the Weyl-type formula; see, e.g., [1], [29]: where |Ω| is the n-dimensional Lebesgue measure of Ω, and ω n is the volume of the unit ball in R n , which implies that the infinite series ∞ j=1 λ 1 λ j 2(η−1) converges for η > 7/4, if n = 3. Therefore, by (61), it is easy to see that . The following lemma, which is a trivial consequence of the compact embedding V ⊂⊂ H, will play a key role in the sequel of this work. We prove it here for the sake of completeness.

Lemma 3
The function u → |u| is weakly continuous with respect to the V -topology, and bounded over B α ⊂ V , α ∈ R + , where B α is as defined in (28).
Proof. The bound follows immediatelly from the definition of B α in (28). Now, we prove that u → |u| is weakly continuous with the V -topology. Indeed, suppose that u j ⇀ u weakly in V , as j → ∞, but u j does not converge strongly to u in H. Then, there exists ǫ > 0, and a subsequence for all k ∈ Z + . Now, because u j k ⇀ u weakly in V , as k → ∞, we have, by the compact embedding V ⊂⊂ H, that there exists a further subsequence u j k ℓ of {u j k } converging strongly (of course, also weakly) in H to a vector field v ∈ V , as ℓ → ∞. Because weak convergence in V implies weak convergence in H, u j k ℓ ⇀ u weakly in H, as ℓ → ∞. By the uniqueness of the weak limit, u = v. Therefore, |u j k ℓ − u| → 0, as ℓ → ∞, which contradicts (62). This proves the lemma.
In order to prove the convergence results in the sequel of this work, we need first to consider a subclass of smoother test functions, T 1 ⊂ T , consisting of functionals satisfying (56), so that the vector fields g i appearing in their definitions possess a higher regularity, namely, g i ∈ V, where V is given in (2) or (3).
Proof. We first check the bound (63). Since ψ is of class C 1 , there exists C > 0 such that for all u ∈ B α 0 , we have Moreover, because for every vector field g ∈ V, we have that and |A η ∂ ℓ ∂x ℓ j g| < C, for some C > 0, for every ℓ ∈ Z + , and j = 1, 2, 3.
Throughout this proof, let {u j } ⊂ B α 0 denote a sequence converging weakly in V to a vector field u in V . Of course, this also implies that u j converges weakly in H to u. By now, it is easy to see by inspection of (57), and by using (58) and (67), that the weak convergence in H of u j to u implies that strongly in H, as j → ∞, for every k ∈ Z + , and η ≥ 0. Now, because the sequence of real numbers given by the inner product in H of a strongly convergent sequence in H, and a weakly convergent sequence in H is convergent, we have as j → ∞. This shows that F α 1 and F α 2 are weakly continuous in V . The uniform bounds follow by (63), and by inspection of (64). Now, concerning the map F α 3 , it is easy to see by the strong convergence in H of u j to u established in Lemma 3, and by the strong convergence in H stated in (69), that strongly in H, as j → ∞. Hence, by the anti-symmetry of the trilinear term, , by the strong convergence of u j to u, and by (72), it is easy to see that as j → ∞. This shows that F α 3 is weakly continuous in V . The uniform bounds follow from the following inequality, which is an easy consequence of (7), for every u ∈ V , and v ∈ W 1,∞ (Ω), and from the following Sobolev's imbedding inequality in 3D: Indeed, because ∂ x j Ψ ′ (u) ∈ D(A), for all u ∈ B α 0 , and j ∈ {1, 2, 3} we have by (63), and by (58), that where C α 0 is uniform for α ∈ [0, α 0 ]. This implies the uniform boundedness for F α 3 for u ∈ B α 0 .

Theorem 4 Given a sequence of invariant probability measures of the 3D
Navier-Stokes-Voigt model, {µ αn }, with α n → 0, as n → ∞, there exists a subsequence, denoted also by {µ αn }, and a Borel probability measure µ on V , such that lim for all Φ : V → R weakly continuous bounded real-valued functionals. Furthermore, the weak limit measure µ is a strong stationary statistical solution of the 3D Navier-Stokes equations.
Proof. Because {µ αn } are invariant measures, by Proposition 1, we know that they are also stationary statistical solutions of the NSV model, and, therefore, they satisfy Conditions 1, 2 and 3 in Definition 1. We prove first that the sequence {µ αn } satisfies the tightness condition (55). Rewriting more explicitly inequality (26), we obtain for every E 1 , E 2 > 0. Now, fixed R > 0, choose E (n) 1 = α 2 n R 2 , and let E 2 → ∞, which is possible because supp µ αn is bounded in V . This yields The first inequality above holds because of the relation and the fact that u 2 is positive. Therefore, denoting by B R 2 the closed ball with radius R 2 in V , we have by (84) that Because the closed balls B R 2 are weakly compact in V , this yields the tightness condition (55), necessary for Theorem 3. Therefore, it is easy to see that Theorem 3 applies to the sequence µ αn , yielding the existence of a subsequence, denoted also by {µ αn }, and a Borel probability measure µ on V , such that (82) holds. Now, we prove that µ is indeed a strong stationary statistical solution of the NSE. Condition 1 of Definition 3 follows directly from (83), and by Fatou's Lemma The proof that Condition 2 of Definition 3 is valid follows directly from Lemma 4 and Lemma 5. Indeed, let us assume first that Ψ ∈ T 1 ⊂ T . Then, by Lemma 4, the functions F 0 i are weakly continuous in V , and uniformly bounded for u ∈ supp µ α , α ∈ [0, α 0 ]. Therefore, Now, because the functions F αn i satisfy the conditions of Lemma 5, we have (86) Now, we are ready to remove the extra hypothesis that assumed Ψ ∈ T 1 ⊂ T . LetΨ ∈ T be defined as in (56) with the functionψ, and the vector fieldsg i ∈ V , instead of ψ and g i . Because V is dense in V , we have that for eachg i , we can find a sequence g j →g i , strongly in V, as j → ∞. Now, let us defineΨ j ∈ T 1 ⊂ T by simply replacing the vector fieldsg i in the definition ofΨ by g We also define (F 0 i ) j as in (64), (65) and (66), withΨ replaced byΨ j . Notice that (86) holds forΨ j . It is easy to check that for every u ∈ V fixed, the following convergence holds: which implies Moreover, it is easy to see that the sequence g can be chosen so that there exists C > 0 such that uniformly in j ∈ Z + . Now, because the right-hand side of the above expression is integrable, by (3), we can use (89) and the Lebesgue Dominated Convergence Theorem to take the limit j → ∞ in (86), yielding This proves that Condition 2 of Definition 3 holds for everyΨ ∈ T . Now, we prove that Condition 3 holds. Let us define and ψ (m) ∈ C ∞ (R m ) as Let us also consider a real-valued function χ ∈ C ∞ (R), and let {w j } ⊂ (C ∞ (Ω)) 3 ∩ V be the orthonormal basis of H composed of the eigenvectors of the Stokes operator, A. Thus, for fixed m, we define for u ∈ B α 0 Ψ χ m (u) := 1 2 χ(ψ (m) ((u, w 1 ), . . . , (u, w m ))).
It is easy to see that Ψ χ m ∈ T , and that Now, because χ ∈ C ∞ (R), we have by the Parseval identity that for every u, h ∈ H, the following limit holds (u, w j )(h, w j ) = χ ′ (|u| 2 )(u, h).
(92) Notice that it follows from the definition of weak convergence that supp µ ⊂ B α 0 , which is bounded in H. Therefore, because f ∈ H, we have that ((Ψ χ m ) ′ (u), f) is uniformly bounded for all u ∈ supp µ, and m ∈ Z + . Thus, we may apply the Lebesgue dominated convergence theorem to obtain Now, let us denote by 1 [E 1 ,E 2 ] : R → R the characteristic step function on the interval of real numbers [E 1 , E 2 ] ⊂ R, i.e.
Consider a sequence of smooth real-valued functions χ i ∈ C ∞ (R), so that χ ′ i (x) → 1 [E 1 ,E 2 ] (x) pointwise, as i → ∞, with |χ i (y)| + |χ ′ i (y)| ≤ C, for every i ∈ Z + . Thus, we may again apply the Lebesgue dominated convergence theorem to deduce that Now, we investigate the limits for the F 0 2 term. We obtain for every u ∈ supp µ ⊂ V , (95) Therefore, because supp µ ⊂ V , we have Thus, we may again apply the Lebesgue dominated convergence theorem to deduce that Now, we calculate the limits for the F 0 3 term. Now, from (5), for u ∈ supp µ ⊂ V , we have B(u, u) ∈ V ′ , and A −1/2 B(u, u) ∈ H. It is easy to prove that the following identity holds B(u, u), w j V ′ ,V = (A −1/2 B(u, u), A 1/2 w j ) = λ 1/2 j (A −1/2 B(u, u), w j ), (99) for every j = 1, 2, . . .. We refer the reader to [3], [15], or [35], for a discussion on the properties of the negative powers of the Stokes operator. Therefore, (u, w j )λ for every u ∈ supp µ. Therefore, by Fatou's lemma: Therefore, using (90), (94), (98), and (102), we deduce that This proves that Condition 3 of Definition 3 is valid, concluding the proof that µ is indeed a strong stationary statistical solution of the Navier-Stokes equations. Now, it follows immediately from Lemma 3 that the function u → |u| 2 is weakly continuous in V , and bounded for u ∈ supp µ α , for α ∈ [0, α 0 ]. Thus, we can prove strong convergence for the kinetic energy.
Now, we prove our result concerning the convergence of the net energy transfer.