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

HAL
open science

# Lecture Notes on a One-Equation Turbulence Closure Model: Modeling and Analysis

Roger Lewandowski

► To cite this version:

Roger Lewandowski. Lecture Notes on a One-Equation Turbulence Closure Model: Modeling and Analysis. 2026. hal-05481070

HAL Id: hal-05481070

https://hal.science/hal-05481070v1

Preprint submitted on 28 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Lecture Notes on a One-Equation Turbulence Closure Model: Modeling and Analysis

Roger Lewandowski

Univ Rennes, IRMAR, UMR CNRS 6625, and Odyssey Team, INRIA Rennes, France, E-mail: Roger.Lewandowski@univ-rennes1.fr

## Abstract

These lecture notes are devoted to the modeling and the analysis of turbulent incompressible flows through Reynolds-averaged Navier-Stokes (RANS) models. After recalling the fundamental limitations of direct numerical simulation of turbulence, we focus on a one-equation closure model coupling a Navier-Stokes-type system for the mean velocity and pressure with an evolution equation for the turbulent kinetic energy (TKE). This system, referred to as the NS-TKE model, is widely used in engineering and geophysical applications.

The first part of the notes is dedicated to the derivation of the model, starting from the Navier-Stokes equations and introducing the Reynolds decomposition, the Boussinesq hypothesis, and the concept of eddy viscosity. A formal derivation of the TKE equation is presented in the framework of isotropic and homogeneous turbulence.

The second part is concerned with the mathematical analysis of the NS-TKE model with homogeneous boundary conditions. Due to the low regularity of the source term in the TKE equation and the presence of eddy viscosities, the system exhibits significant analytical difficulties. To address these issues, a truncated and regularized version of the model is introduced. The existence of weak solutions is established via the Galerkin method, and the passage to the limit is carried out using energy methods and compactness arguments, with a distinction between the two- and three-dimensional cases.

These notes are intended for graduate students and researchers in mathematical fluid dynamics and aim to provide a self-contained introduction to both the modeling and analysis of RANS turbulence models.

**Keywords:** Turbulence, Navier-Stokes equations RANS models, turbulent kinetic energy, eddy viscosity, weak solutions, parabolic equation with r.h.s in $L^1$.

**Mathematics Subject Classification:** 35Q30, 76D05, 76F02, 35B45.

## Contents

<table><tr><td>Introduction</td><td>3</td></tr><tr><td>Part I : Modeling</td><td>5</td></tr><tr><td>1 General fluid mechanics</td><td>5</td></tr><tr><td>1.1 Fluid equations</td><td>6</td></tr><tr><td>1.1.1 Framework</td><td>6</td></tr><tr><td>1.1.2 Mass conservation</td><td>7</td></tr></table>
<table><tr><td>1.1.3</td><td>Equation of motion</td><td>8</td></tr><tr><td>1.2</td><td>Stress tensor</td><td>9</td></tr><tr><td>1.2.1</td><td>Definition</td><td>9</td></tr><tr><td>1.2.2</td><td>Decomposition of the stress tensor, dynamic pressure</td><td>10</td></tr><tr><td>1.2.3</td><td>Incompressible Navier-Stokes Equations</td><td>11</td></tr><tr><td>1.3</td><td>Energy balance</td><td>12</td></tr><tr><td>1.3.1</td><td>Boundary conditions</td><td>12</td></tr><tr><td>1.3.2</td><td>Derivation of the energy balance</td><td>13</td></tr></table>

## 2 Turbulence modeling 14

<table><tr><td>2.1</td><td>Reynolds decomposition</td><td>14</td></tr><tr><td>2.2</td><td>Mean field equation and Reynolds Stress</td><td>16</td></tr><tr><td>2.3</td><td>Isotropy</td><td>17</td></tr><tr><td>2.4</td><td>Boussinesq assumption</td><td>17</td></tr><tr><td>2.5</td><td>The turbulent kinetic energy equation</td><td>18</td></tr><tr><td>2.5.1</td><td>Homogeneous flows</td><td>18</td></tr><tr><td>2.5.2</td><td>Derivation of the TKE equation</td><td>19</td></tr><tr><td>2.6</td><td>Closures in the TKE equation</td><td>21</td></tr></table>

# Part II: Mathematical analysis 22

## 3 Mathematical framework and basic estimates 22

<table><tr><td>3.1</td><td>Notation and assumptions</td><td>22</td></tr><tr><td>3.2</td><td>Standard results: reminders</td><td>23</td></tr><tr><td>3.3</td><td>Basic estimates for the velocity</td><td>25</td></tr><tr><td>3.4</td><td>Estimates for the TKE</td><td>26</td></tr></table>

## 4 The truncated NS-TKE system 29

<table><tr><td>4.1</td><td>Truncated NS-TKE model and weak formulation</td><td>29</td></tr><tr><td>4.1.1</td><td>Truncated transport terms</td><td>30</td></tr><tr><td>4.1.2</td><td>Diffusion terms</td><td>32</td></tr><tr><td>4.1.3</td><td>Production and dissipation terms in the TKE equation</td><td>33</td></tr><tr><td>4.1.4</td><td>Weak solutions of the truncated NS-TKE system</td><td>33</td></tr><tr><td>4.1.5</td><td>Existence of a weak solution to the truncated NS-TKE model: statement and plan of the proof</td><td>34</td></tr><tr><td>4.2</td><td>Approximate solutions and estimates</td><td>34</td></tr><tr><td>4.2.1</td><td>Finite dimensional spaces</td><td>35</td></tr><tr><td>4.2.2</td><td>Galerkin method</td><td>36</td></tr><tr><td>4.2.3</td><td>Energy estimates and compactness</td><td>37</td></tr><tr><td>4.3</td><td>Passing to the limit in the equations</td><td>38</td></tr><tr><td>4.3.1</td><td>Extracting subsequences</td><td>38</td></tr><tr><td>4.3.2</td><td>Passing to the limit in the fluid equation</td><td>40</td></tr><tr><td>4.3.3</td><td>The energy method</td><td>43</td></tr><tr><td>4.3.4</td><td>Passing to the limit in the TKE-equation</td><td>46</td></tr><tr><td>4.4</td><td>Additional properties of the limit</td><td>47</td></tr><tr><td>4.4.1</td><td>Entropy calculus</td><td>48</td></tr><tr><td>4.4.2</td><td>Maximum principle</td><td>50</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>5 Existence result for the NS-TKE system</td>
      <td>51</td>
    </tr>
    <tr>
      <td>5.1 Main estimates</td>
      <td>52</td>
    </tr>
    <tr>
      <td>5.1.1 Estimates for {v<sub>α</sub>}<sub>α&gt;0</sub></td>
      <td>52</td>
    </tr>
    <tr>
      <td>5.1.2 Estimates for {k<sub>α</sub>}<sub>α&gt;0</sub> by Boccardo-Gallouët Lemma</td>
      <td>52</td>
    </tr>
    <tr>
      <td>5.2 Estimates and test functions spaces when N = 2</td>
      <td>54</td>
    </tr>
    <tr>
      <td>5.2.1 Variational formulation for the fluid equation</td>
      <td>54</td>
    </tr>
    <tr>
      <td>5.2.2 More estimates for the TKE equation</td>
      <td>55</td>
    </tr>
    <tr>
      <td>5.3 Weak solutions for N = 2</td>
      <td>57</td>
    </tr>
    <tr>
      <td>5.3.1 Compactness and extracting subsequences</td>
      <td>58</td>
    </tr>
    <tr>
      <td>5.3.2 Passing to the limit in the fluid equation.</td>
      <td>59</td>
    </tr>
    <tr>
      <td>5.3.3 Passing to the limit in the TKE equation.</td>
      <td>60</td>
    </tr>
    <tr>
      <td>5.3.4 L<sup>∞</sup>([0, T], L<sup>1</sup>(Ω)) estimate.</td>
      <td>61</td>
    </tr>
    <tr>
      <td>5.4 Case N = 3</td>
      <td>62</td>
    </tr>
    <tr>
      <td>5.4.1 Main result</td>
      <td>62</td>
    </tr>
    <tr>
      <td>5.4.2 Additional estimates for the TKE equation.</td>
      <td>64</td>
    </tr>
    <tr>
      <td>5.4.3 Extracting subsequences and passing to the limit in the equations</td>
      <td>65</td>
    </tr>
    <tr>
      <td>Appendices</td>
      <td>66</td>
    </tr>
    <tr>
      <td>A Tool box</td>
      <td>66</td>
    </tr>
    <tr>
      <td>A.1 Aubin-Lions Lemma</td>
      <td>66</td>
    </tr>
    <tr>
      <td>A.2 Time derivative and initial data</td>
      <td>67</td>
    </tr>
    <tr>
      <td>A.3 Boccardo-Gallouët lemma</td>
      <td>67</td>
    </tr>
    <tr>
      <td>A.4 Miscellaneous</td>
      <td>68</td>
    </tr>
    <tr>
      <td>B Stampacchia Theorem</td>
      <td>70</td>
    </tr>
    <tr>
      <td>B.1 Statement</td>
      <td>70</td>
    </tr>
    <tr>
      <td>B.2 Reduction to the case G(x) = |x|.</td>
      <td>71</td>
    </tr>
    <tr>
      <td>B.3 Proof of the result when G(x) = |x|</td>
      <td>71</td>
    </tr>
    <tr>
      <td>References</td>
      <td>72</td>
    </tr>
  </tbody>
</table>

Introduction

These lecture notes are based on a doctoral course first delivered at the Vietnam Insti-
tute for Advanced Study in Mathematics (VIASM, Hanoï, Vietnam) in May 2025, then
at the Centre International de Rencontres Mathématiques (CIRM, Marseille, France) in
October 2025, and finally at the Department of Mathematics of the University of Rennes
(France) in March 2026. The course is intended for PhD students and researchers in
mathematical fluid dynamics who are interested in numerical models for turbulent incom-
pressible flows.

Turbulence remains one of the most fascinating and challenging problems in both classical and modern science. Its general structure and statistical properties were extensively investigated throughout the twentieth century by pioneering researchers. Among the most influential contributions are those of Prandtl [46], who introduced the concept of the boundary layer; Richardson [48], who proposed the idea of an energy cascade of vortices; Taylor [55], who developed fundamental concepts of isotropic turbulence and experimental techniques; Landau [30], who suggested a phenomenological scenario for the transition to turbulence; Kolmogorov [29], who established the statistical theory of fully developed
turbulence; and Onsager [43], who predicted anomalous energy dissipation in ideal fluids. Later, Batchelor [4] and Frisch [25] provided comprehensive and rigorous analyses of turbulent flows. Finally, the work of Ruelle and Takens [49] introduced the modern viewpoint of turbulence as a manifestation of deterministic chaos in dynamical systems. Together, these foundational studies constitute the backbone of our current understanding of turbulent phenomena.

We do not aim to improve upon all of these works; rather, we focus on the numerical perspective, specifically on how turbulent flows can be simulated. The fundamental equations governing the motion of an incompressible fluid in a domain are the Navier–Stokes equations. Since the seminal work of Leray [33], the existence of solutions that may exhibit turbulence—now commonly referred to as weak solutions—has been established. However, the problem raised by Leray, namely whether a smooth solution can develop a singularity in finite time, remains open [21] and is, at the time of writing, one of the Clay Millennium Prize Problems, endowed with a reward of one million dollars.

Unfortunately, it is well known that directly simulating turbulent flows from the Navier–Stokes equations is not feasible, as it would require an enormous number of degrees of freedom, far beyond the capacity of any computer available today [16]. This limitation motivates the introduction of so-called *turbulence models*. Historically, the first turbulence model is generally attributed to Smagorinsky [51], developed in the context of numerical weather prediction. This model belongs to the class now known as *Large Eddy Simulation* (LES) models, in which only the large-scale eddies of the flow are resolved, while the effects of the smaller scales are parametrized (see also Lesieur, Métais, and Comte [34]). Shortly thereafter, the well-known $k-\epsilon$ model was introduced by Launder and Spalding [31] (see also Mohamadi and Pironneau [40]). This model belongs to the class of *Reynolds-averaged Navier-Stokes* (RANS) models, which aim at capturing the statistical mean behavior of turbulent flows.

This course is devoted to the modeling and analysis of a one-equation RANS closure model, widely used in engineering, oceanography, meteorology, and related applications. The model consists of a Navier–Stokes-like equation with an eddy viscosity, satisfied by the statistical means of the velocity and pressure, denoted by $\bar{v}$ and $\bar{p}$, coupled with an equation for the turbulent kinetic energy $k$ (TKE). The TKE measures the deviation of the instantaneous velocity from its mean value, and the eddy viscosity is a function of $k$. We refer to this system as the *NS–TKE model*.

In the first part of the course, we focus on the modeling of the RANS closure, starting from the Navier–Stokes equations. We then introduce the notion of statistical means, the Reynolds decomposition, the Boussinesq assumption, and the concept of eddy viscosity. This framework allows us to formally derive, from the kinetic energy balance, the equation governing the turbulent kinetic energy in the case of isotropic and homogeneous flows. This modeling part has been designed to be accessible to mathematicians. In this sense, it provides a concise synthesis of the first four chapters of [16], rewritten in a simplified and readable form.

In the second part, we address the mathematical analysis of the NS–TKE model with homogeneous boundary conditions. This line of research began with the seminal paper [36], published in 1997, and was subsequently developed in numerous works, including [6, 11, 37, 44]. It is worth noting that in [16], the NS–TKE model was analyzed with wall laws as boundary conditions, which differs from the homogeneous boundary conditions considered here.

The NS–TKE model exhibits the usual mathematical difficulties associated with the Navier–Stokes equations, together with additional challenges arising from the presence of eddy
viscosities, which are assumed to be bounded throughout the analysis. Moreover, the source term in the equation satisfied by the turbulent kinetic energy belongs only to $L^1$, and not to a higher regularity space, making this equation a nonstandard parabolic problem. To overcome these difficulties, we introduce a regularized system based on truncation functions, referred to as the *truncated NS-TKE model*. After establishing an appropriate variational formulation, we prove the existence of weak solutions using the Galerkin method. We then pass to the limit as the truncation parameter tends to infinity, distinguishing between the two-dimensional and three-dimensional cases, in order to obtain weak solutions to the NS-TKE model in a carefully defined sense.

Throughout the text, we derive a priori estimates step by step. These estimates form the backbone of the analysis and determine the Sobolev and Bochner function spaces used. Estimates for the mean velocity follow the classical energy estimates for the Navier-Stokes equations, provided the eddy viscosity remains strictly nonnegative. Estimates for the TKE equation rely on Stampacchia's theorem and the Boccardo-Gallouët lemma, both of which are carefully presented. Passing to the limit in the equations is based on what we call the *energy method*, relying on the energy balance satisfied by the velocity whenever it can be justified; this energy balance serves as a guiding thread throughout the text.

Finally, in order to make the manuscript as self-contained as possible, we have included an appendix gathering standard technical tools, together with a complete and elementary proof of Stampacchia's theorem, which we were unable to find in the existing literature. A comprehensive bibliography is provided throughout the text.

## Acknowledgements

The author would like to express his sincere gratitude to Dinh-Duong Nguyen and the Vietnam Institute for Advanced Study in Mathematics (VIASM) for their warm hospitality in Hanoi in May 2025, where part of this course was first delivered. He also warmly thanks Francesco Fanelli and the ANR project CRISIS for the invitation to give this course at the CIRM (Marseille) in October 2025.

The author is grateful to the Institut de Recherche Mathématique de Rennes (IRMAR) and to the INRIA project team ODYSSEY for their constant support over the years. Finally, he wishes to acknowledge Claude Bardos and Olivier Pironneau, who introduced him long ago to the theory of turbulence, as well as François Murat, who taught him how to study partial differential equations with measure-valued or $L^1$ source terms.

# Part I : Modeling

## 1 General fluid mechanics

The aim of this section is to provide the reader with a simple and clear derivation of the Navier-Stokes equations from fundamental physical laws: mass conservation and Newton's second law of motion. The challenging part lies in modeling the internal forces, which are generally inferred from experimental observations. We will restrict our focus to incompressible Newtonian fluids, which will be clarified by the end of this section.
# 1.1 Fluid equations

To fix ideas, we work in three-dimensional space, although two-dimensional fluids may also be considered, as we will do in the mathematical analysis part.

## 1.1.1 Framework

Let

* $t \in \mathbb{R}_+^*$ be the time variable,

* $M = (x_1, x_2, x_3)$ be any point in space.

We are interested in the behavior of a fluid moving within a fixed domain $\Omega \subset \mathbb{R}^3$, whose boundary is denoted by $\Gamma = \partial\Omega$. We assume that the domain $\Omega$ is of class $C^1$, which means, that the boundary $\Gamma$ is smooth enough to allow for a proper mathematical treatment. More precisely, this regularity implies that $\Gamma$ can be covered by a finite number of open sets:

$$ \Gamma \subset \bigcup_{i=1}^{n} O_i, $$

with each $O_i$ being an open subset of $\mathbb{R}^3$. For every $i = 1, \dots, n$, there exists a $C^1$ diffeomorphism $H_i$ mapping $O_i$ onto the reference square $C = [-1, 1] \times [-1, 1]$, such that

$$ H_i(\Omega \cap O_i) = [-1, 1] \times [-1, 0] $$

and

$$ H_i(\Gamma \cap O_i) = [-1, 1] \times \{0\}. $$

This formulation allows us to describe the geometry of the boundary in a consistent and controlled way.

At any point $M \in \Gamma$, we define the outward unit normal vector by

$$ \mathbf{n} = \mathbf{n}(M) = \begin{pmatrix} n_1(M) \\ n_2(M) \\ n_3(M) \end{pmatrix}, \quad |\mathbf{n}| = \sqrt{n_1^2 + n_2^2 + n_3^2} = 1. $$

Moreover, if $M$ belongs to the chart intersection $\Gamma \cap O_i$, then this normal vector transforms under $H_i$ into the simple form

$$ H_i(\mathbf{n}) = \begin{pmatrix} 1 \\ 0 \end{pmatrix}. $$

We now recall the Green's formula (also known as the divergence theorem), which connects a volume integral over the domain to a surface integral over its boundary. For any $C^1$ vector field $\mathbf{u}$ defined on $\bar{\Omega}$,

$$ \mathbf{u}(M) = \begin{pmatrix} u_1(M) \\ u_2(M) \\ u_3(M) \end{pmatrix} \in C^1(\bar{\Omega}), $$

the formula reads:

$$ (1.1) \qquad \iiint_{\Omega} \mathrm{div} \mathbf{u}(M) dV(M) = \iint_{\Gamma} \mathbf{u}(M) \cdot \mathbf{n}(M) dS(M), $$
where

$$ \mathrm{div} \mathbf{u} = \frac{\partial u_1}{\partial x_1} + \frac{\partial u_2}{\partial x_2} + \frac{\partial u_3}{\partial x_3} $$

is the divergence of the vector field $\mathbf{u}$.

To simplify notation in what follows, we will write $\int_{\Omega}$ and $\int_{\Gamma}$ instead of $\iiint_{\Omega}$ and $\iiint_{\Gamma}$, whenever the context makes the meaning clear.

### 1.1.2 Mass conservation

Let $\mathbf{v}(t, M)$ denote the Eulerian velocity of the fluid at time $t$ and position $M$:

$$ \mathbf{v}(t, M) = \begin{pmatrix} v_1(t, M) \\ v_2(t, M) \\ v_3(t, M) \end{pmatrix} \in C^2(\mathbb{R}^+ \times \bar{\Omega}, \mathbb{R}^3). $$

The Eulerian description consists in observing the motion of the fluid at a fixed point $M$ in space, measuring the velocity of all fluid particles passing through a small neighborhood of $M$ as time evolves, just like speed cameras on highways, which record the speed of each car at a fixed location.

We also consider the quantity

$$ \rho = \rho(t, M) \in C^1(\mathbb{R}^+ \times \bar{\Omega}, \mathbb{R}), $$

which denotes the fluid density at time $t$ and position $M$, defined as follows. Let $\omega \subset \Omega$. The total mass of the fluid contained in $\omega$ at time $t$ is then given by

$$ (1.2) \qquad m(t, \omega) = \int_{\omega} \rho(t, M) dM. $$

One of the cornerstone equations in fluid mechanics is the mass conservation equation, which expresses the fact that mass cannot be created or destroyed during the fluid's evolution.

**Theorem 1.1.** *The mass conservation equation, valid at every point $(t, M) \in \mathbb{R}_+ \times \Omega$, reads:*

$$ (1.3) \qquad \frac{\partial \rho}{\partial t} + \mathrm{div}(\rho \mathbf{v}) = 0. $$

*Proof.* Let $\omega \subset \subset \Omega$ be a fixed open set strictly contained in $\Omega$, with boundary $\partial\omega$, and let $\mathbf{n} = \mathbf{n}(\mathbf{x})$ denote the outward unit normal vector at any $\mathbf{x} \in \partial\omega$.

Let $t \in \mathbb{R}_+$ and let $\delta t > 0$ be an infinitesimal time increment. We measure the amount of fluid mass leaving $\omega$ through $\partial\omega$ during the interval $[t, t + \delta t]$, denoted by $\Delta m(t, \delta t, \omega)$. Some mass may also enter $\omega$ through $\partial\omega$ during this interval, but it is still counted here as leaving, with a non-positive sign.

We compute $\Delta m(t, \delta t, \omega)$ in two ways. First, we perform a Taylor expansion with respect to $\delta t$. Second, we carry out a local geometric analysis at $\partial\omega$ and express the flux as a boundary integral, which we then convert into a volume integral using Green's formula (1.1). Equating the resulting expressions yields the mass conservation law.

We start with the algebraic mass variation:

$$ \Delta m(t, \delta t, \omega) = m(t, \omega) - m(t + \delta t, \omega). $$
Applying Taylor's formula,

$$ (1.4) \qquad \begin{aligned} \Delta m(t, \delta t, \omega) &= -\delta t \frac{d}{dt} \int_{\omega} \rho(t, \mathbf{x}) \, d\mathbf{x} + O(\delta t^2) \\ &= -\delta t \int_{\omega} \partial_t \rho(t, \mathbf{x}) \, d\mathbf{x} + O(\delta t^2). \end{aligned} $$

Now consider the boundary flux. Fix $\mathbf{x} \in \partial\omega$ and let $\delta S \subset \partial\omega$ be an infinitesimal surface element centered at $\mathbf{x}$, with outward unit normal $\mathbf{n}(\mathbf{x})$. Assume that $(\mathbf{v}, \rho)$ are approximately constant near $(t, \mathbf{x})$ and equal to $(\mathbf{v}(t, \mathbf{x}), \rho(t, \mathbf{x}))$.

A particle located at $y$ near $\mathbf{x}$ at time $t$ moves to $\mathbf{y} + \mathbf{v}(t, \mathbf{x}) \delta t$ at $t + \delta t$. Thus, the mass of fluid crossing $\delta S$ during $[t, t + \delta t]$ is

$$ \delta m = \rho(t, \mathbf{x}) \mathbf{v}(t, \mathbf{x}) \cdot \mathbf{n}(\mathbf{x}) \delta S \delta t. $$

Integrating over $\partial\omega$,

$$ (1.5) \qquad \Delta m(t, \delta t, \omega) = \delta t \int_{\partial\omega} \rho(t, \mathbf{x}) \mathbf{v}(t, \mathbf{x}) \cdot \mathbf{n}(\mathbf{x}) \, dS. $$

Applying Green's formula (1.1) yields,

$$ (1.6) \qquad \Delta m(t, \delta t, \omega) = \delta t \int_{\omega} \mathrm{div}(\rho(t, \mathbf{x}) \mathbf{v}(t, \mathbf{x})) \, d\mathbf{x}. $$

Combining (1.4) and (1.6), dividing by $\delta t$, and letting $\delta t \to 0$, we obtain

$$ (1.7) \qquad \int_{\omega} [\partial_t \rho(t, \mathbf{x}) + \mathrm{div}(\rho(t, \mathbf{x}) \mathbf{v}(t, \mathbf{x}))] \, d\mathbf{x} = 0, $$

for every $t \ge 0$ and every subdomain $\omega \subset \Omega$. By the regularity of $\mathbf{v}$ and $\rho$, this implies the pointwise mass conservation equation (1.3). $\square$

An equivalent form of the mass conservation equation, using the product rule, is

$$ (1.8) \qquad \frac{\partial \rho}{\partial t} + \mathbf{v} \cdot \nabla \rho + \rho \mathrm{div} \mathbf{v} = 0. $$

In the special case where $\rho$ is constant in time and space (i.e. the fluid is homogeneous and incompressible), this equation simplifies to

$$ (1.9) \qquad \mathrm{div} \mathbf{v} = 0. $$

In this situation, we say that the fluid is incompressible.

**Remark 1.1.** On the basis of scale analysis, there are certain fluid flows where the mass density is not actually constant, yet the condition $\mathrm{div} \mathbf{v} = 0$ is still used as an approximation. This is notably the case for large-scale oceanic flows.

### 1.1.3 Equation of motion

The equation of motion is based on Newton's fundamental law. To express it, we first need the acceleration of a fluid particle located at position $M$ at time $t$, whose velocity is given by $\mathbf{v} = \mathbf{v}(t, M)$. At time $t + \delta t$, the particle has approximately moved to the position

$$ M + \delta M \approx M + \mathbf{v} \delta t. $$
Let $i = 1, 2, 3$, and consider the corresponding change in the $i$-th component of velocity:

$$ \delta v_i = \underbrace{v_i(t + \delta t, x_1 + v_1 \delta t, x_2 + v_2 \delta t, x_3 + v_3 \delta t)}_{i^{th} \text{ velocity component at } M + \mathbf{v}\delta t} - \underbrace{v_i(t, x_1, x_2, x_3)}_{i^{th} \text{ velocity component at } M}. $$

Therefore,

$$ \begin{aligned} \delta v_i(t, M) &= \left( \frac{\partial v_i}{\partial t}(t, M) + \mathbf{v}(t, M) \cdot \nabla v_i(t, M) \right) \delta t + o(\delta t) \\ &= \left( \frac{\partial v_i}{\partial t}(t, M) + \sum_{j=1}^{3} v_j(t, M) \frac{\partial v_i}{\partial x_j}(t, M) \right) \delta t + o(\delta t). \end{aligned} $$

At any point $(t, M) \in \mathbb{R}_+ \times \Omega$, the acceleration is given by the limit

$$ \lim_{\delta t \to 0} \frac{\delta v_i}{\delta t} = \frac{\partial v_i}{\partial t} + \mathbf{v} \cdot \nabla v_i. $$

We define the vector field $(\mathbf{v} \cdot \nabla)\mathbf{v}$ by

$$ (\mathbf{v} \cdot \nabla)\mathbf{v} = \begin{pmatrix} \mathbf{v} \cdot \nabla v_1 \\ \mathbf{v} \cdot \nabla v_2 \\ \mathbf{v} \cdot \nabla v_3 \end{pmatrix}, $$

and the vector field $\frac{d\mathbf{v}}{dt}$ by

$$ (1.10) \qquad \frac{d\mathbf{v}}{dt} = \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v}, $$

which is called the total (or material) derivative. This represents the acceleration of the fluid particle. Newton's second law of motion then yields

$$ (1.11) \qquad \rho \frac{d\mathbf{v}}{dt} = \sum \text{Forces} = F_{\text{int}} + F_{\text{ext}}, $$

where the forces acting on the fluid have been decomposed into internal and external components. This equation holds per unit volume, where:

$$ (1.12) \qquad \left\{ \begin{array}{ll} \rho \frac{d\mathbf{v}}{dt} : & \text{mass} \times \text{acceleration (per unit volume),} \\ F_{\text{int}} : & \text{internal forces (per unit volume)} \rightarrow \text{to be modeled,} \\ F_{\text{ext}} : & \text{external forces (per unit volume)} \rightarrow \text{given.} \end{array} \right. $$

## 1.2 Stress tensor

### 1.2.1 Definition

If the external forces are assumed to be known, as is the case with gravity and/or magnetic forces in geophysical flows, the challenge lies in modeling the internal forces, i.e. the forces exerted on a fluid particle by the surrounding fluid.

It can be shown that there exists a 3 × 3 symmetric tensor [5, 27, 57]

$$ \boldsymbol{\sigma} = (\sigma_{ij})_{1 \le i,j \le 3} = \boldsymbol{\sigma}(t, M) \in C^1(\mathbb{R}^+ \times \bar{\Omega}, M_{3 \times 3}(\mathbb{R})), $$
such that, for any volume $\omega \subset \Omega$, the force exerted by the fluid on the boundary $\partial\omega$ at time $t$ is given by

$$ (1.13) \qquad \underbrace{F_{int}(t, \omega)}_{\text{physical unit: Newton.Meter}^{-3}} = \int_{\partial\omega} \boldsymbol{\sigma} \cdot \mathbf{n} = \underbrace{\div \boldsymbol{\sigma}}_{\text{by Green's formula}} , $$

where $\mathbf{n}$ denotes the outward unit normal vector on $\partial\omega$, and

$$ \boldsymbol{\sigma} \cdot \mathbf{n}|_{\omega} \in C^1(\mathbb{R}^+ \times \partial\omega, \mathbb{R}^3), \quad \div \boldsymbol{\sigma} \in C^0(\mathbb{R}^+ \times \bar{\Omega}, \mathbb{R}^3) $$

are vector fields expressed by (1.14) below. We refer to $\boldsymbol{\sigma}$ as the stress tensor, and let

$$ \underbrace{\Sigma_i = (\sigma_{i1}, \sigma_{i2}, \sigma_{i3})}_{\text{$i$\textsuperscript{th} row of $\boldsymbol{\sigma}$}} \in C^1(\mathbb{R}^+ \times \bar{\Omega}, M_{1,3}(\mathbb{R})). $$

Then the vector fields $\boldsymbol{\sigma} \cdot \mathbf{n}$ and $\div \boldsymbol{\sigma}$ are given by

$$ (1.14) \qquad \boldsymbol{\sigma} \cdot \mathbf{n} = \begin{pmatrix} \Sigma_1 \cdot \mathbf{n} \\ \Sigma_2 \cdot \mathbf{n} \\ \Sigma_3 \cdot \mathbf{n} \end{pmatrix}, \quad \div \boldsymbol{\sigma} = \begin{pmatrix} \operatorname{div} \Sigma_1 \\ \operatorname{div} \Sigma_2 \\ \operatorname{div} \Sigma_3 \end{pmatrix}. $$

Now take $\omega = B(M, r)$, for a sufficiently small $r > 0$. Then,

$$ \forall M' \in B(M, r), \quad \operatorname{div} \boldsymbol{\sigma}(t, M') = \operatorname{div} \boldsymbol{\sigma}(t, M) + o(r), $$

and therefore,

$$ \int_{\omega} \operatorname{div} \boldsymbol{\sigma}(t, M') dM' \approx \frac{4}{3} \pi r^3 \operatorname{div} \boldsymbol{\sigma}(t, M) + o(r^4). $$

This justifies the definition:

$$ \underbrace{F_{int}}_{\text{physical unit: Newton.Meter}^{-3}} = \text{force per unit volume} = \operatorname{div} \boldsymbol{\sigma}. $$

In summary, Equation (1.11) combined with (1.13) becomes:

$$ (1.15) \qquad \rho \frac{d\mathbf{v}}{dt} = \sum \text{Forces} = \operatorname{div} \boldsymbol{\sigma} + F_{\text{ext}}, $$

### 1.2.2 Decomposition of the stress tensor, dynamic pressure

We now have to model $\boldsymbol{\sigma}$. It is convenient to decompose it as

$$ (1.16) \qquad \boldsymbol{\sigma} = \mathbb{D} - p\mathbf{I}, $$

where $p$ denotes the dynamic pressure, defined by

$$ (1.17) \qquad p = -\frac{1}{3}\operatorname{tr} \boldsymbol{\sigma}. $$

Here, $\mathbf{I}$ represents the identity tensor, while $\mathbb{D}$ is the deviatoric part of $\boldsymbol{\sigma}$.

In the case of Newtonian fluids such as water, air, or oil, the constitutive relation is given by [5, 20]:

$$ (1.18) \qquad \mathbb{D} = 2\mu D\mathbf{v}, $$

where

$$ D\mathbf{v} = \frac{1}{2} (\nabla\mathbf{v} + (\nabla\mathbf{v})^T), \quad \nabla\mathbf{v} = \left(\frac{\partial v_i}{\partial x_j}\right)_{1\le i,j\le 3} \in C^1(\mathbb{R}^+ \times \bar{\Omega}, M_{3\times 3}(\mathbb{R})). $$
• $D\mathbf{v}$ is the rate-of-deformation tensor, represented as a matrix field,

• $\mu$ is the dynamic viscosity, a scalar quantity that typically depends on temperature.

Relation (1.18) is primarily derived from experimental observations. Therefore, (1.15) becomes

$$ (1.19) \qquad \rho \frac{d\mathbf{v}}{dt} = \text{div}(2\mu D\mathbf{v} - p\mathbf{I}) + F_{\text{ext}} $$

**Remark 1.2.** More explicitly, the tensor $D\mathbf{v}$ takes the form

$$ D\mathbf{v} = \begin{pmatrix} D_1 \\ D_2 \\ D_3 \end{pmatrix} = \begin{pmatrix} \frac{\partial v_1}{\partial x_1} & \frac{1}{2} \left( \frac{\partial v_1}{\partial x_2} + \frac{\partial v_2}{\partial x_1} \right) & \frac{1}{2} \left( \frac{\partial v_1}{\partial x_3} + \frac{\partial v_3}{\partial x_1} \right) \\ \frac{1}{2} \left( \frac{\partial v_2}{\partial x_1} + \frac{\partial v_1}{\partial x_2} \right) & \frac{\partial v_2}{\partial x_2} & \frac{1}{2} \left( \frac{\partial v_2}{\partial x_3} + \frac{\partial v_3}{\partial x_2} \right) \\ \frac{1}{2} \left( \frac{\partial v_3}{\partial x_1} + \frac{\partial v_1}{\partial x_3} \right) & \frac{1}{2} \left( \frac{\partial v_3}{\partial x_2} + \frac{\partial v_2}{\partial x_3} \right) & \frac{\partial v_3}{\partial x_3} \end{pmatrix}. $$

therefore,

$$ \sigma = \begin{pmatrix} \Sigma_1 \\ \Sigma_2 \\ \Sigma_3 \end{pmatrix} = \begin{pmatrix} 2\mu \frac{\partial v_1}{\partial x_1} - p & \mu \left( \frac{\partial v_1}{\partial x_2} + \frac{\partial v_2}{\partial x_1} \right) & \mu \left( \frac{\partial v_1}{\partial x_3} + \frac{\partial v_3}{\partial x_1} \right) \\ \mu \left( \frac{\partial v_2}{\partial x_1} + \frac{\partial v_1}{\partial x_2} \right) & 2\mu \frac{\partial v_2}{\partial x_2} - p & \mu \left( \frac{\partial v_2}{\partial x_3} + \frac{\partial v_3}{\partial x_2} \right) \\ \mu \left( \frac{\partial v_3}{\partial x_1} + \frac{\partial v_1}{\partial x_3} \right) & \mu \left( \frac{\partial v_3}{\partial x_2} + \frac{\partial v_2}{\partial x_3} \right) & 2\mu \frac{\partial v_3}{\partial x_3} - p \end{pmatrix}. $$

### 1.2.3 Incompressible Navier-Stokes Equations

We now turn to incompressible flows, where $\mu$ is constant and

$$ \text{div}\mathbf{v} = 0. $$

In this setting, we observe that

$$ (1.20) \qquad \text{div}(2\mu D\mathbf{v}) = 2\mu \begin{pmatrix} \text{div}D_1 \\ \text{div}D_2 \\ \text{div}D_3 \end{pmatrix} = \mu(\Delta\mathbf{v} + \nabla(\text{div}\mathbf{v})) = \mu\Delta\mathbf{v}. $$

That is,

$$ \text{div}\mathbf{v} = 0 \Rightarrow \Delta\mathbf{v} = \text{div}(2D\mathbf{v}). $$

Moreover, one easily verifies that

$$ (1.21) \qquad \text{div}(p\mathbf{I}) = \nabla p. $$

Assuming a constant fluid density $\rho$, we introduce the kinematic viscosity

$$ \nu = \frac{\mu}{\rho}, $$

and denote by $\mathbf{f}$ the external force per unit mass:

$$ \mathbf{f} = \frac{F_{\text{ext}}}{\rho}. $$

Combining (1.19), (1.20) and (1.21), we obtain the incompressible Navier-Stokes equations:

$$ (1.22) \qquad \begin{cases} \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} - \nu \Delta \mathbf{v} + \nabla p = \mathbf{f} & \text{in } \Omega, \\ \text{div} \mathbf{v} = 0 & \text{in } \Omega. \end{cases} $$
Componentwise, System (1.22) is written as

$$ (1.23) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial v_i}{\partial t} + \sum_{j=1}^{3} v_j \frac{\partial v_i}{\partial x_j} - \nu \sum_{j=1}^{3} \frac{\partial^2 v_i}{\partial x_j^2} + \frac{\partial p}{\partial x_i} = f_i, \quad i = 1, 2, 3, \\[1.5em] \displaystyle \operatorname{div} \mathbf{v} = \sum_{k=1}^{3} \frac{\partial v_k}{\partial x_k} = 0, \end{array} \right. $$

where

$$ \mathbf{f} = \begin{pmatrix} f_1 \\ f_2 \\ f_3 \end{pmatrix}. $$

**Remark 1.3.** Since $\operatorname{div} \mathbf{v} = 0$, for each $i \in \{1, 2, 3\}$ we have

$$ [(\mathbf{v} \cdot \nabla)\mathbf{v}]_i = \sum_{j=1}^{3} v_j \frac{\partial v_i}{\partial x_j} = \sum_{j=1}^{3} \frac{\partial}{\partial x_j} (v_i v_j) = [\operatorname{div}(\mathbf{v} \otimes \mathbf{v})]_i, $$

so the momentum equation can be rewritten as

$$ (1.24) \qquad \frac{\partial \mathbf{v}}{\partial t} + \operatorname{div}(\mathbf{v} \otimes \mathbf{v}) - \nu \Delta \mathbf{v} + \nabla p = \mathbf{f}, $$

where

$$ (1.25) \qquad \mathbf{v} \otimes \mathbf{v} = (v_i v_j)_{1 \le i,j \le 3}. $$

## 1.3 Energy balance

One of the central tools in mathematical fluid dynamics is the energy balance, which provides the main estimates on the velocity field. We derive it formally in what follows. Before doing it, we discuss the boundary conditions.

Throughout the rest of the paper, given two matrices $A = (a_{ij})$ and $B = (b_{ij})$, we set

$$ (1.26) \qquad A : B = \sum_{i,j} a_{ij} b_{ij}, \quad |A| = \left( \sum_{i,j} a_{ij}^2 \right)^{1/2} . $$

Moreover, as above, $\mathbf{n} = \begin{pmatrix} n_1 \\ n_2 \\ n_3 \end{pmatrix}$ denotes the outward normal vector at the boundary $\Gamma$.

### 1.3.1 Boundary conditions

Boundary conditions encompass:

* The initial data that describes the state of the system at time $t = 0$: $\forall M \in \Omega$, $\mathbf{v}_0 = \mathbf{v}(0, M) = \mathbf{v}_{t=0}$,

* The behavior of the velocity at the boundary $\Gamma$ of the flow domain $\Omega$.
In practice, the initial data typically comes from observational measurements. The behavior of the fluid at the boundary has been the subject of much discussion since the establishment of the Navier-Stokes equations, originally formulated by Navier in 1832. Two main types of boundary conditions are considered: friction-based conditions, such as Navier or Manning laws [16], and Dirichlet-type conditions.

When the boundary is fixed, it is common to assume that the fluid adheres to the boundary $\Gamma$, meaning the velocity vanishes there. This is referred to as the no-slip condition, which we will adopt throughout this paper.

Accordingly, the Navier-Stokes boundary value problem is formulated as follows:

$$ (1.27) \qquad \left\{ \begin{array}{ll} \displaystyle \frac{\partial \mathbf{v}}{\partial t} + \operatorname{div}(\mathbf{v} \otimes \mathbf{v}) - \nu \Delta \mathbf{v} + \nabla p = \mathbf{f} & \text{in } \Omega, \\ \operatorname{div} \mathbf{v} = 0 & \text{in } \Omega, \\ \mathbf{v} = 0 & \text{on } \Gamma, \\ \mathbf{v}_{t=0} = \mathbf{v}_0 & \text{in } \Omega. \end{array} \right. $$

**Remark 1.4.** There is no need to prescribe a boundary condition for the pressure, as it is not a prognostic variable. From a mathematical standpoint, it is treated as a Lagrange multiplier associated with the incompressibility constraint, particularly when the no-slip boundary condition is applied to the velocity field.

### 1.3.2 Derivation of the energy balance

Assume that $\mathbf{v} \in C^2(\mathbb{R}_+ \times \bar{\Omega}, \mathbb{R}^3)$ and $p \in C^1(\mathbb{R}_+ \times \bar{\Omega}, \mathbb{R})$. We take the scalar product of the momentum equation in (1.27) with $\mathbf{v}$, and integrate over $\Omega$ at a fixed time by using the Green formula. Using Einstein's summation convention and recalling that $\operatorname{div}\mathbf{v} = 0$ and $\mathbf{v}|_{\Gamma} = 0$, we obtain:

$$ (i) \quad \int_{\Omega} \frac{\partial \mathbf{v}}{\partial t} \cdot \mathbf{v} = \frac{1}{2} \frac{d}{dt} \int_{\Omega} |\mathbf{v}|^2; $$

$$ \begin{align*} (ii) \quad \int_{\Omega} (\mathbf{v} \cdot \nabla) \mathbf{v} \cdot \mathbf{v} &= \int_{\Omega} v_j \frac{\partial v_i}{\partial x_j} \cdot v_i = \frac{1}{2} \int_{\Omega} v_j \frac{\partial |\mathbf{v}|^2}{\partial x_j} \\ &= \int_{\Gamma} |\mathbf{v}|^2 \mathbf{v} \cdot \mathbf{n} - \frac{1}{2} \int_{\Omega} |\mathbf{v}|^2 \frac{\partial v_j}{\partial x_j} = -\frac{1}{2} \int_{\Omega} |\mathbf{v}|^2 \operatorname{div} \mathbf{v} = 0; \end{align*} $$

$$ (iii) \quad \int_{\Omega} \nabla p \cdot \mathbf{v} = \int_{\Omega} \frac{\partial p}{\partial x_i} v_i = \int_{\Gamma} p \mathbf{v} \cdot \mathbf{n} - \int_{\Omega} \frac{\partial v_i}{\partial x_i} p = - \int_{\Omega} p \operatorname{div} \mathbf{v} = 0. $$

Finally, using the boundary condition $\mathbf{v}|_{\Gamma} = 0$, we obtain:

$$ -\nu \int_{\Omega} \Delta \mathbf{v} \cdot \mathbf{v} = -\nu \int_{\Omega} \Delta v_i v_i = \underbrace{-\nu \int_{\Gamma} v_i \frac{\partial v_i}{\partial x_j} n_j}_{=0} + \nu \int_{\Omega} \frac{\partial v_i}{\partial x_j} \frac{\partial v_i}{\partial x_j} = \nu \int_{\Omega} |\nabla \mathbf{v}|^2 $$

We deduce from the previous formal calculations that at any time $t > 0$, for all $s \in [0, t]$,

$$ \frac{1}{2} \frac{d}{dt} \int_{\Omega} |\mathbf{v}(s, M)|^2 dM + \nu \int_{\Omega} |\nabla \mathbf{v}(s, M)|^2 dM = \int_{\Omega} \mathbf{f}(s, M) \cdot \mathbf{v}(s, M) dM. $$
Integrating over time on $[0, t]$, we obtain the energy balance:

$$ (1.28) \qquad \underbrace{\frac{1}{2} \int_{\Omega} |\mathbf{v}(t, M)|^2 dM}_{\text{kinetic energy}} + \nu \underbrace{\int_{0}^{t} \int_{\Omega} |\nabla \mathbf{v}(s, M)|^2 dM ds}_{\text{molecular dissipation}} = \frac{1}{2} \int_{\Omega} |\mathbf{v}_0(M)|^2 dM + \underbrace{\int_{0}^{t} \int_{\Omega} \mathbf{f}(s, M) \cdot \mathbf{v}(s, M) dM ds}_{\text{work of the external force}} $$

## 2 Turbulence modeling

The behavior of a fluid is governed by the Reynolds number:

$$ R = \frac{UL}{\nu}, $$

where $U$ is a characteristic velocity scale of the flow, and $L$ is a characteristic length scale. For example, in the case of a flow around a boat moving at a steady velocity, $U$ can be taken as the boat's speed, and $L$ as its length.

According to the similarity law, two flows with the same Reynolds number exhibit similar behavior [5, 16]. When the Reynolds number is small, the flow is smooth and orderly: this is known as laminar flow. As the Reynolds number increases, instabilities develop, leading to the formation of waves and eddies at various scales. The flow becomes increasingly chaotic, resembling the turbulent motion of a mountain stream. This regime is referred to as turbulent flow.

Experimentally, turbulence typically arises when the Reynolds number exceeds approximately $R = 2000–3000$. In geophysical contexts, Reynolds numbers often exceed $10^3$.

According to Kolmogorov's theory, simulating a turbulent flow requires approximately $R^{9/4}$ degrees of freedom, which far exceeds the capabilities of current computational resources [45]. As a result, the Navier-Stokes equations cannot be directly solved for turbulent flows, and turbulence models must be developed for practical simulations.

There are two main approaches to modeling turbulence:

* **Large Eddy Simulation (LES):** Only the large-scale structures of the flow are resolved, while smaller scales are modeled.

* **Reynolds-Averaged Navier-Stokes (RANS):** The equations are averaged to describe the mean flow, with turbulence effects incorporated through additional modeling.

In the following sections, we will focus on constructing RANS models.

### 2.1 Reynolds decomposition

Assume that we measure the velocity of a turbulent flow $N$ times at a given point $M$ under identical conditions, $\mathbf{v}_1 = \mathbf{v}(t_1, M), \dots, \mathbf{v}_N = \mathbf{v}(t_N, M)$. Even when using a high-precision technique for fluid-velocity measurement, such as Particle Image Velocimetry [47], the results still vary because of uncontrollable infinitesimal perturbations, which means that the $\mathbf{v}_k$'s are not identical. However, the statistical properties remain stable. Namely, the ensemble average

$$ \bar{\mathbf{v}}_N = \frac{1}{N} \sum_{k=1}^{N} \mathbf{v}_k $$
converges, for large $N$, to a vector $\bar{v}$ that remains approximately constant for each series of measurements. This is as if

$$ \mathbf{v} = \mathbf{v}(t, M, \omega), \quad \omega \in (E, P), $$

where $E$ is a probabilistic space, with $P$ as probability measure, and

$$ (2.1) \qquad \bar{\mathbf{v}} = \mathbb{E}_{\omega}(\mathbf{v}) \quad \text{i.e.} \quad \bar{\mathbf{v}}(t, M) = \int_{E} \mathbf{v}(t, M, \omega)dP(\omega). $$

We decompose the velocity and the pressure as

$$ \mathbf{v} = \bar{\mathbf{v}} + \mathbf{v}' \quad \text{and} \quad \mathbf{p} = \bar{\mathbf{p}} + \mathbf{p}', $$

where $\mathbf{v}'$ and $\mathbf{p}'$ are the fluctuations. This decomposition is known as the Reynolds decomposition.

In the following, we will deal with $k$ and $\varepsilon$, which are essential in turbulence modeling and are given by:

* The turbulent kinetic energy (TKE):

$$ (2.2) \qquad k = \frac{1}{2} \overline{|\mathbf{v}'|^2}. $$

* The turbulent dissipation $\varepsilon$:

$$ (2.3) \qquad \varepsilon = 2\nu \overline{|D\mathbf{v}'|^2}, $$

recalling that

$$ D\mathbf{v}' = \frac{1}{2}(\nabla\mathbf{v}' + (\nabla\mathbf{v}')^t). $$

From dimensional analysis, we can obtain a length scale $\ell$ from $k$ and $\varepsilon$, given by

$$ (2.4) \qquad \ell = \frac{k\sqrt{k}}{\varepsilon}, $$

which is known as the Prandtl mixing length. It represents the typical distance an eddy can travel before being mixed into the surrounding flow by turbulence. One of the main challenges in turbulence modeling is the determination of the equations satisfied by $k$ and $\varepsilon$. Here we will consider a one order closure model, which means that we will just derive an equation for $k$, $\varepsilon$ being calculated by the relation (2.4). In numerical simulations, $\ell$ is often taken to be equal to the mesh size [14]. In [44], we derived an analytical formula for $\ell$ in a boundary layer, based on direct numerical simulations.

We conclude this subsection with an important property of the probabilistic filter given in (2.1), namely that it is linear and commutes with any differential operator. This is known as the Reynolds rules. In particular, it satisfies

$$ \overline{\partial_t \mathbf{v}} = \partial_t \overline{\mathbf{v}}, \quad \overline{\nabla \mathbf{v}} = \nabla \overline{\mathbf{v}} \quad \text{and} \quad \overline{\nabla p} = \nabla \overline{\mathbf{p}}. $$

We will need the following simple result.

**Lemma 2.1.** Assume the filter is given by (2.1). Then the fluctuation mean vanishes, i.e., $\overline{\mathbf{v}}' = 0$.
Proof. We have

$$ (2.5) \qquad \bar{\mathbf{v}} = \bar{\overline{\mathbf{v}}} + \overline{\mathbf{v}'} = \bar{\overline{\mathbf{v}}} + \overline{\mathbf{v}'} $$

using the linearity of the expectation. We observe that $\overline{\mathbf{v}} = \overline{\mathbf{v}}$ because

$$ \overline{\mathbf{v}}(t, M) = \int_E \overline{\mathbf{v}}(t, M)dP(\omega) = \overline{\mathbf{v}}(t, M) \int_E dP(\omega) = \overline{\mathbf{v}}(t, M), $$

since $P$ is a probability measure. This yields

$$ \overline{\mathbf{v}'} = \int_E \mathbf{v}' dP(\omega) = 0. $$

□

## 2.2 Mean field equation and Reynolds Stress

We aim to derive a PDE satisfied by the mean field $\overline{\mathbf{v}}$ and $\overline{\mathbf{p}}$, thanks to Lemma 2.1. To do so, we take the expectation of the NSE by using the Reynolds rules. The issue comes from the nonlinear term $\overline{\mathbf{v}} \otimes \overline{\mathbf{v}}$. Since $\mathbf{v} = \overline{\mathbf{v}} + \mathbf{v}'$, we have

$$ \overline{\mathbf{v}} \otimes \overline{\mathbf{v}} = \overline{\mathbf{v}} \otimes \overline{\mathbf{v}} + \overline{\mathbf{v}'} \otimes \overline{\mathbf{v}'} + \overline{\mathbf{v}} \otimes \overline{\mathbf{v}'} + \overline{\mathbf{v}'} \otimes \overline{\mathbf{v}}. $$

Notice that for any scalar fields $F$ and $G$, it is easily checked that

$$ \overline{FG} = \overline{FG} = \overline{F} \overline{G}. $$

Therefore, we obtain by Lemma 2.1, $\overline{v_i v_j'} = \overline{v_i} \overline{v_j'} = 0$ and $\overline{v_i' v_j} = \overline{v_i'} \overline{v_j} = 0$, which yields

$$ \overline{\mathbf{v}} \otimes \overline{\mathbf{v}'} = \overline{\mathbf{v}} \otimes \overline{\mathbf{v}'} = 0, \quad \overline{\mathbf{v}'} \otimes \overline{\mathbf{v}} = \overline{\mathbf{v}'} \otimes \overline{\mathbf{v}} = 0. $$

By consequence, we get

$$ \overline{\mathbf{v}} \otimes \overline{\mathbf{v}} = \overline{\mathbf{v}} \otimes \overline{\mathbf{v}} + \overline{\mathbf{v}'} \otimes \overline{\mathbf{v}'} $$

which yields by setting

$$ \partial_t = \frac{\partial}{\partial t}, $$

the averaged NSE as follows

$$ (2.6) \qquad \left\{ \begin{array}{ll} \partial_t \overline{\mathbf{v}} + \operatorname{div}(\overline{\mathbf{v}} \otimes \overline{\mathbf{v}}) + \operatorname{div}(\overline{\mathbf{v}'} \otimes \overline{\mathbf{v}'}) - \nu \Delta \overline{\mathbf{v}} + \nabla \overline{\mathbf{p}} & = \overline{\mathbf{f}} \quad \text{in } \Omega, \\ \operatorname{div} \overline{\mathbf{v}} & = 0 \quad \text{in } \Omega, \\ \overline{\mathbf{v}} & = 0 \quad \text{on } \Gamma, \\ \overline{\mathbf{v}} & = \overline{\mathbf{v}}_0 \quad \text{at } t = 0. \end{array} \right. $$

This addresses the issue of determining the nonlinear term

$$ (2.7) \qquad \mathbf{R} = \overline{\mathbf{v}'} \otimes \overline{\mathbf{v}'}, $$

which is called as the Reynolds stress tensor. Since $\mathbf{v}' = (v'_1, v'_2, v'_3)$, the Reynolds stress is the following matrix

$$ \mathbf{R} = \begin{pmatrix} \overline{(v'_1)^2} & \overline{v'_1 v'_2} & \overline{v'_1 v'_3} \\ \overline{v'_2 v'_1} & \overline{(v'_2)^2} & \overline{v'_2 v'_3} \\ \overline{v'_3 v'_1} & \overline{v'_3 v'_2} & \overline{(v'_3)^2} \end{pmatrix}. $$

**Remark 2.1.** The TKE is related to the trace of the Reynolds stress by the following relation

$$ (2.8) \qquad k = \frac{1}{2} (\overline{(v'_1)^2} + \overline{(v'_2)^2} + \overline{(v'_3)^2}) = \frac{1}{2} \operatorname{tr}(\mathbf{R}). $$
## 2.3 Isotropy

Most turbulence laws apply to isotropic flows, meaning that the statistical properties of the flow are identical in all spatial directions, whether one considers correlations, velocity fluctuations, or turbulent structures. To be more specific, let us consider the second-order velocity correlation tensor defined by

$$ \mathbb{B} = \mathbb{B}(x, r) = (B_{i,j}(x, r))_{1\le i,j\le 3}, $$

where

$$ (2.9) \quad \forall x \in \Omega, \quad \forall r \in \mathbb{R}^3 \text{ s.t. } x+r \in \Omega, \quad B_{i,j}(x,r) = \overline{v_i(x+r)v_j(x)}. $$

Then isotropy means

$$ (2.10) \quad \begin{array}{l} \forall u, v \in \mathbb{R}^3, \quad \forall x \in \Omega, \quad \forall r \in \mathbb{R}^3 \text{ s.t. } x+r \in \Omega, \\ \forall Q \in O_3(\mathbb{R}), \quad (\mathbb{B}(x, Qr)Qu, Qv) = (\mathbb{B}(x, r)u, v), \end{array} $$

where $O_3(\mathbb{R})$ denotes the orthogonal group in $\mathbb{R}^3$. Of course, most turbulent flows are not isotropic. However, it is conjectured that this property holds at sufficiently small scales.

## 2.4 Boussinesq assumption

We aim to express the Reynolds stress tensor (2.7) in terms of mean quantities that can be calculated, based on the Boussinesq assumption. The key idea of Boussinesq is the following. Instead of trying to compute the Reynolds stress tensor directly (which is unknown and messy), assume turbulence behaves somewhat like an ordinary viscous fluid.

* In laminar flow, viscous stresses are proportional to the rate of deformation $D\mathbf{v}$.

* In turbulence, we make a conceptual copy-paste: we introduce an “effective viscosity”, also the turbulent or eddy viscosity $\nu_t$, that plays the same role as the kinematic viscosity, but for the averaged fields.

Physically:

* Turbulent eddies stir and mix momentum like little spoons homogenizing soup.

* The stronger the velocity gradient, the more intense this “stirring,” so the correlations $v'_i v'_j$ must also grow.

Therefore, the Boussinesq assumption consists of postulating the following relation:

$$ (2.11) \quad \mathbf{R} = -\nu_{\text{turb}} D \bar{\mathbf{v}} + \frac{2}{3} k \mathbf{I}, $$

where $\nu_{\text{turb}}$ is an eddy viscosity, the $\frac{2}{3}k\mathbf{I}$ term restores the isotropic turbulent pressure, and also because of (2.8) and $\text{tr}D\bar{\mathbf{v}} = 0$ due to the incompressibility constraint. Notice that in [7] we were able to prove that the Reynolds stress is dissipative in the case of the long-time average filter, which in some sense validates Formula (2.11).

This assumption is not exact. Indeed, it presumes turbulence is “nearly isotropic” and that its momentum-diffusing effect resembles an ordinary viscosity. But it makes the equations “closed” and solvable, so far we get a formula for $\nu_{\text{turb}}$, which is one of turbulence modeling aim.
According to Formula (2.11) combined with Equations (2.6), the momentum equation for the means becomes

$$ (2.12) \qquad \partial_t \bar{\nu} + \mathrm{div}(\bar{\nu} \otimes \bar{\nu}) - \mathrm{div}((2\nu + \nu_{\mathrm{turb}})D\mathbf{v}) + \nabla \bar{p} = \bar{\mathbf{f}}, $$

where we still denote by $\bar{p}$ the modified pressure $\bar{p} + \frac{2}{3}k$, so far without causing confusion. In this presentation, we assume that

$$ (2.13) \qquad \nu_{\mathrm{turb}} = \nu_{\mathrm{turb}}(k), $$

which is a classical assumption (among others) and which yields what is called a first-order closure model. By definition of the physical dimensions,

$$ \underbrace{\mathrm{kinematic viscosity}}_{\mathrm{unit: } m^2 s^{-1}} = \underbrace{\mathrm{length}}_{\mathrm{unit: } m} \times \underbrace{\mathrm{velocity}}_{\mathrm{unit: } ms^{-1}}. $$

Hence, by dimensional analysis,

$$ (2.14) \qquad \nu_{\mathrm{turb}} = C \ell \sqrt{k}, $$

where $C$ is a dimensionless constant determined by experimental data, and the Prandtl mixing length $\ell$ is given by Formula (2.4).

## 2.5 The turbulent kinetic energy equation

From now, we may denote the velocity by $\mathbf{v}(t, x_1, x_2, x_3)$ or $\mathbf{v}(t, x)$ instead of $\mathbf{v}(t, M)$ when needed, as well as for any field related to the flow.

### 2.5.1 Homogeneous flows

Generally speaking, turbulence models are established for isotropic and homogeneous flows. We have already mentioned the notion of isotropy in section 2.3. The aim of this subsection is to discuss homogeneity (see also [4]).

**Definition 2.1.** The flow is said to be homogeneous if, for any fields $\psi, \varphi$ related to the flow (velocity component, fluctuation, energy, pressure...), for any $r \in \mathbb{R}$, $x \in \Omega$ such that $x + r \in \Omega$, and for any $t \in \mathbb{R}_+$, the correlation function

$$ (2.15) \qquad \overline{\psi(t, x + r)\varphi(t, x)} $$

depends only on $t$ and $r$, and not on $x$.

To derive an equation for $k$, we need the following result on homogeneous flows, which is analogous to an integration by parts formula for functions with compact support.

**Lemma 2.2.** Assume the flow homogeneous. Then for any fields $\psi, \varphi$ related to it, we have

$$ (2.16) \qquad \forall i = 1, \dots, N, \quad \forall (t, x) \in \mathbb{R}_+ \times \Omega, \quad \overline{\psi(t, x) \partial_i \varphi(t, x)} = - \overline{\partial_i \psi(t, x) \varphi(t, x)}. $$

*Proof.* Let $i = 1, \dots, N$, $h_i \in \mathbb{R}^*$, and for $i = 2, \dots, N-1$, let $\mathbf{r}_i = (0, \dots, h_i, \dots, 0)$, $\mathbf{r}_1 = (h_1, 0, \dots, 0)$, $\mathbf{r}_N = (0, \dots, 0, h_N)$. Therefore,

$$ (2.17) \qquad \overline{\psi(t, x) \partial_i \varphi(t, x)} = \lim_{h_i \to 0} \left( \frac{1}{h_i} \overline{\psi(t, x) [\varphi(t, x + \mathbf{r}_i) - \varphi(t, x)]} \right). $$
We have

$$ \overline{\psi(t, x)[\varphi(t, x + \mathbf{r}_i) - \varphi(t, x)]} = \overline{\psi(t, x)\varphi(t, x + \mathbf{r}_i)} - \overline{\psi(t, x)\varphi(t, x)}. $$

By the homogeneity of the flow,

$$ \overline{\psi(t, x)\varphi(t, x + \mathbf{r}_i)} = \overline{\psi(t, x - \mathbf{r}_i)\varphi(t, x)}, $$

which yields

$$ (2.18) \qquad \frac{1}{h_i} \overline{\psi(t, x)[\varphi(t, x + \mathbf{r}_i) - \varphi(t, x)]} = -\left(-\frac{1}{h_i}\right) \overline{[\psi(t, x - \mathbf{r}_i) - \psi(t, x)]\varphi(t, x)}. $$

As

$$ \lim_{h_i \to 0} \left( -\frac{1}{h_i} \overline{[\psi(t, x - \mathbf{r}_i) - \psi(t, x)]\varphi(t, x)} \right) = \overline{\partial_i \psi(t, x)\varphi(t, x)}, $$

Formula (2.16) follows from (2.17) combined with (2.18). □

### 2.5.2 Derivation of the TKE equation

We are now able to derive an equation for $k$ for an homogeneous flow, following what was initially done in [40] (see also in [16]), by applying the integration by parts formula (2.16), starting from the kinetic energy of the fluctuation $e$, defined by

$$ (2.19) \qquad e = \frac{1}{2} |\mathbf{v}'|^2 = e(t, M, \omega), $$

satisfying

$$ (2.20) \qquad k = \bar{e}, $$

so that

$$ (2.21) \qquad e = k + e'. $$

We will establish the general equation for $k$ in three steps:

1. write the equation satisfied by the fluctuating field $(\mathbf{v}', p')$,

2. derive from it an equation for $e$,

3. take the expectation of the equation for $e$.

To carry out this program, we need to consider the total dissipation $\mathbb{D}$ given by

$$ (2.22) \qquad \mathbb{D} = 2\nu |D\bar{\mathbf{v}}|^2. $$

Following the Reynolds rules, we decompose $\mathbb{D}$ as

$$ (2.23) \qquad \overline{\mathbb{D}} = 2\nu |D\bar{\mathbf{v}}|^2 + 2\nu |\overline{D\mathbf{v}'}|^2 = 2\nu |D\bar{\mathbf{v}}|^2 + \varepsilon. $$

**Theorem 2.1.** Assume that the flow is homogeneous. Then the TKE $k$ satisfies the following equation

$$ (2.24) \qquad \partial_t k + \bar{\mathbf{v}} \cdot \nabla k + \operatorname{div}(e' \mathbf{v}') = -\mathbf{R} : \nabla \bar{\mathbf{v}} - \varepsilon. $$
Proof. We start from the Reynolds decomposition $\mathbf{v} = \bar{\mathbf{v}} + \mathbf{v}'$ and $p = \bar{p} + p'$. We subtract the system satisfied by $(\bar{\mathbf{v}}, \bar{p})$ from the NSE, and by noting first that

$$ (\mathbf{v} \cdot \nabla)\mathbf{v} - (\bar{\mathbf{v}} \cdot \nabla)\bar{\mathbf{v}} = (\mathbf{v} \cdot \nabla)(\mathbf{v} - \bar{\mathbf{v}}) + (\mathbf{v} \cdot \nabla)\bar{\mathbf{v}} - (\bar{\mathbf{v}} \cdot \nabla)\bar{\mathbf{v}} \\ = \mathbf{v} \cdot \nabla \mathbf{v}' + \mathbf{v}' \cdot \nabla \bar{\mathbf{v}}, $$

which leads to

$$ \partial_t \mathbf{v}' + \mathbf{v} \cdot \nabla \mathbf{v}' + \mathbf{v}' \cdot \nabla \bar{\mathbf{v}} - \operatorname{div}(2\nu D \mathbf{v}') + \nabla p' - \operatorname{div} \mathbf{R} = 0, \\ \operatorname{div} \mathbf{v}' = 0. $$

Taking the dot product of the above system with $\mathbf{v}'$ yields the following equation for $e$

$$ (2.25) \qquad \partial_t e + \underbrace{\mathbf{v} \cdot \nabla e}_{I} + \underbrace{(\mathbf{v}' \otimes \mathbf{v}') : \nabla \bar{\mathbf{v}}}_{II} - \underbrace{(\operatorname{div}(2\nu D \mathbf{v}'))}_{III} \cdot \mathbf{v}' + \underbrace{\nabla p' \cdot \mathbf{v}'}_{IV} - \underbrace{(\operatorname{div} \mathbf{R})}_{V} \cdot \mathbf{v}' = 0. $$

We take the expectation of this equation, and study each term one after the other.

I. We have the identity

$$ (2.26) \qquad \mathbf{v} \cdot \nabla e = \operatorname{div}(e \mathbf{v}) $$

since $\operatorname{div} \mathbf{v} = 0$. Using the Reynolds decomposition, we find that

$$ (2.27) \qquad \overline{e \mathbf{v}} = \overline{\mathbf{v}} \overline{e} + \overline{e' \mathbf{v}'} = \overline{\mathbf{v}} k + \overline{e' \mathbf{v}'} $$

and because also $\operatorname{div} \mathbf{v}' = 0$, combining the relation $\mathbf{v} \cdot \nabla k = \operatorname{div}(k \mathbf{v})$ with (2.26) and (2.27), we obtain

$$ (2.28) \qquad \overline{\mathbf{v} \cdot \nabla e} = \overline{\mathbf{v}} \cdot \nabla k + \operatorname{div}(\overline{e' \mathbf{v}'}) $$

II. We have

$$ (2.29) \qquad \overline{(\mathbf{v}' \otimes \mathbf{v}') : \nabla \bar{\mathbf{v}}} = \overline{(\mathbf{v}' \otimes \mathbf{v}') : \nabla \bar{\mathbf{v}}} = \mathbf{R} : \nabla \bar{\mathbf{v}}. $$

III. Similarly,

$$ -\overline{(\operatorname{div}(2\nu D \mathbf{v}')) \cdot \mathbf{v}'} = 2\nu \overline{D \mathbf{v}'} : \nabla \mathbf{v}', $$

which yields

$$ (2.30) \qquad -\overline{(\operatorname{div}(2\nu D \mathbf{v}')) \cdot \mathbf{v}' 2\nu \overline{D \mathbf{v}' : \nabla \mathbf{v}'}} = 2\nu \overline{D \mathbf{v}'} : \overline{D \mathbf{v}'} = 2\nu \overline{|D \mathbf{v}'|^2} = \varepsilon. $$

(because $|\nabla \mathbf{v}'|^2 = |\nabla(\mathbf{v}')^t|^2$, then $D\mathbf{v}' : \nabla \mathbf{v}' = D\mathbf{v}' : \nabla(\mathbf{v}')^t$, hence the formula)

IV. Because of the homogeneity assumption, we have

$$ (2.31) \qquad \overline{\nabla p' \cdot \mathbf{v}'} = -\overline{p' \operatorname{div} \mathbf{v}'} = 0, $$

V. From $\overline{\mathbf{v}'} = 0$ and since $\operatorname{div} \mathbf{R}$ does not depend on $\omega$, we deduce

$$ (2.32) \qquad \overline{(\operatorname{div} \mathbf{R}) \cdot \mathbf{v}'} = (\operatorname{div} \mathbf{R}) \cdot \overline{\mathbf{v}'} = 0. $$

To conclude, Equation (2.24) for $k$ results from the identity (2.20) combined with the Reynolds rules, that gives,

$$ \overline{\partial_t e} = \partial_t k, $$

applied in the equation (2.25) combined with (2.28), (2.29), (2.30), (2.31) and (2.32) $\square$
## 2.6 Closures in the TKE equation

Our work is not yet complete because of the terms

$$ \mathrm{div}(\overline{e'\mathbf{v}'}), \quad -\mathbf{R}: \nabla \overline{\mathbf{v}} $$

in the TKE equation (2.24), which must be expressed in terms of mean quantities. This procedure is referred to as closing the equation, which is the aim of this section. Recall that we have set, according to the Boussinesq assumption,

$$ \mathbf{R} = -\nu_{\mathrm{turb}} D\overline{\mathbf{v}} + \frac{2}{3} k \mathbf{I}, $$

which yields

$$ -\mathbf{R}: \nabla \overline{\mathbf{v}} = \nu_{\mathrm{turb}} D\overline{\mathbf{v}}: \nabla \overline{\mathbf{v}} - \frac{2}{3} k \mathbf{I}: \nabla \overline{\mathbf{v}}. $$

On the one hand,

$$ \frac{2}{3} k \mathbf{I}: \nabla \overline{\mathbf{v}} = \frac{2}{3} k \mathrm{tr} \nabla \overline{\mathbf{v}} = \frac{2}{3} k \mathrm{div} \overline{\mathbf{v}} = 0. $$

On the other hand,

$$ \nu_{\mathrm{turb}} D\overline{\mathbf{v}}: \nabla \overline{\mathbf{v}} = \nu_{\mathrm{turb}} D\overline{\mathbf{v}}: D\overline{\mathbf{v}} = \nu_{\mathrm{turb}} |D\overline{\mathbf{v}}|^2. $$

Therefore, the TKE equation (2.25) becomes

$$ \partial_t k + \overline{\mathbf{v}} \cdot \nabla k + \mathrm{div}(\overline{e'\mathbf{v}'}) = \nu_{\mathrm{turb}} |D\overline{\mathbf{v}}|^2 - \varepsilon. $$

The remaining task is to treat the term

$$ \mathrm{div}(\overline{e'\mathbf{v}'}). $$

Following the Boussinesq assumption, generally speaking we assume that transport by fluctuations in the mean flow is modeled as diffusion. This question was discussed in many works, in particular for periodic structures by means of homogenization [15, 39] We thus introduce an eddy diffusivity $\mu_{\mathrm{turb}} = \mu_{\mathrm{turb}}(k)$ such that

$$ \overline{e'\mathbf{v}'} = -\mu_{\mathrm{turb}} \nabla \overline{e} = -\mu_{\mathrm{turb}} \nabla k, $$

which yields

$$ (2.33) \qquad \partial_t k + \overline{\mathbf{v}} \cdot \nabla k - \mathrm{div}(\mu_{\mathrm{turb}} \nabla k) = \nu_{\mathrm{turb}} |D\overline{\mathbf{v}}|^2 - \varepsilon. $$

One option is to derive an additional equation for $\varepsilon$, yielding the full $k - \varepsilon$ model. Since in this paper we restrict ourselves to a single-equation closure, we do not establish here the equation for $\varepsilon$ (see [40]). Instead, we use the relation (2.4):

$$ \varepsilon = \frac{k^{\frac{3}{2}}}{\ell}, $$

where $\ell$ denotes the Prandtl mixing length, and as $k$ is an energy, we assume here that $k \ge 0$. The eddy diffusivity $\mu_{\mathrm{turb}}$ remains to be specified. By the same dimensional argument as in formula (2.14), we write

$$ (2.34) \qquad \mu_{\mathrm{turb}} = C_{\mu} \ell \sqrt{k}. $$
In conclusion, combining the equation (2.12) satisfied by $\bar{\mathbf{v}}$ and $\bar{p}$ with the incompressibility constraint and the TKE equation (2.33), we obtain what we call the NS-TKE model:

$$ (2.35) \quad \left\{ \begin{array}{ll} \partial_t \bar{\mathbf{v}} + \operatorname{div}(\bar{\mathbf{v}} \otimes \bar{\mathbf{v}}) - \operatorname{div}((2\nu + \nu_{\text{turb}}(k))D\bar{\mathbf{v}}) + \nabla \bar{p} = \bar{\mathbf{f}} & \text{in } \Omega, \\ \operatorname{div}\bar{\mathbf{v}} = 0 & \text{in } \Omega, \\ \partial_t k + \bar{\mathbf{v}} \cdot \nabla k - \operatorname{div}(\mu_{\text{turb}}(k)\nabla k) = \nu_{\text{turb}}(k)|D\bar{\mathbf{v}}|^2 - \frac{k^3}{\ell} & \text{in } \Omega. \end{array} \right. $$

We assume that $\bar{\mathbf{v}}$ also satisfies homogeneous boundary conditions at $\Gamma$, as $\mathbf{v}$. Normally, $k = |\bar{\mathbf{v}}|^2$ at $\Gamma$, therefore $k|_{\Gamma} = 0$. To summerize, the boundary conditions are the following:

$$ (2.36) \quad \left\{ \begin{array}{ll} \bar{\mathbf{v}} = 0 & \text{at } \Gamma, \\ k = 0 & \text{at } \Gamma, \\ \bar{\mathbf{v}}_{t=0} = \bar{\mathbf{v}}_0 & \text{in } \Omega, \\ k_{t=0} = k_0 & \text{in } \Omega. \end{array} \right. $$

This concludes the modeling part. In the next part we will prove the existence of weak solutions to the NS-TKE system (2.35)-(2.36).

**Remark 2.2.** Many mathematical works aim to provide rigorous justifications for laws of turbulence and related equations; see for example [3, 7, 17, 18, 19, 23, 24, 32]. In particular, the work of Bardos and Titi [3] provides important insights on the connection between weak solutions of the Navier-Stokes equations and Reynolds-averaged quantities, showing under which conditions such averages satisfy a form of the energy balance. Despite these advances, deriving the turbulent kinetic energy (TKE) equation rigorously from a Leray solution of the Navier-Stokes equations remains an open problem, highlighting the ongoing challenges in giving a fully rigorous mathematical foundation to turbulence modeling.

# Part II: Mathematical analysis

## 3 Mathematical framework and basic estimates

### 3.1 Notation and assumptions

From now on, we write $(\mathbf{v}, p, \mathbf{f})$ instead of $(\bar{\mathbf{v}}, \bar{p}, \bar{\mathbf{f}})$ for simplicity, and $\nu_{\text{turb}}$ instead of $2\nu + \nu_{\text{turb}}$. We assume that $\Omega \subset \mathbb{R}^N$, $N = 2, 3$, is bounded and of class $C^1$, and we still set $\Gamma = \partial\Omega$. For simplicity, and without loss of generality, we replace the deformation operator $Du$ by the gradient operator $\nabla u$ in the fluid equation, since this does not significantly affect the results thanks to Korn's inequality (see [16]). Therefore, the NS-TKE model (2.35), together with the boundary conditions (2.36), yields the following boundary value problem,

$$ (3.1) \quad \left\{ \begin{array}{ll} \partial_t \mathbf{v} + \operatorname{div}(\mathbf{v} \otimes \mathbf{v}) - \operatorname{div}(\nu_{\text{turb}}(k)\nabla \mathbf{v}) + \nabla p = \mathbf{f} & \text{in } \Omega, \\ \operatorname{div}\mathbf{v} = 0 & \text{in } \Omega, \\ \partial_t k + \mathbf{v} \cdot \nabla k - \operatorname{div}(\mu_{\text{turb}}(k)\nabla k) = \nu_{\text{turb}}(k)|D\mathbf{v}|^2 - \frac{k|k|^{\frac{1}{2}}}{\ell} & \text{in } \Omega, \\ \mathbf{v} = 0 & \text{on } \Gamma, \\ k = 0 & \text{on } \Gamma, \\ \mathbf{v}_{t=0} = \mathbf{v}_0 & \text{in } \Omega, \\ k_{t=0} = k_0 & \text{in } \Omega. \end{array} \right. $$
The main aim is to establish the existence of a solution — at least in the distributional sense — to this system, if any. We denote by $x = (x_1, \dots, x_N) \in \Omega$ the standard point in $\Omega$ (instead of $M$), let $T$ be the final time with $t \in [0, T]$, and define

$$ (3.2) \qquad Q_T = [0, T] \times \Omega. $$

Assume that the following conditions on $\nu_{\text{turb}}$, $\mu_{\text{turb}}$, $\ell$ and $f$ are satisfied:

1. The functions $\nu_{\text{turb}}$ and $\mu_{\text{turb}}$ are of class $C^1$ over $\mathbb{R}$, nonnegative, and bounded above and below: $\exists \nu > 0, \beta > 0$ such that

$$ (3.3) \qquad \forall k \in \mathbb{R}, \quad 0 < \nu \le \nu_{\text{turb}}(k), \mu_{\text{turb}}(k) \le \beta < \infty. $$

2. The function $\ell = \ell(t, x)$ is continuous and bounded below: $\exists \ell_0 > 0$ such that

$$ (3.4) \qquad \forall (t, x) \in Q_T, \quad 0 < \ell_0 \le \ell(t, x). $$

3. The source term $f$ satisfies

$$ (3.5) \qquad f \in C([0, T], L^2(\Omega)^N). $$

## 3.2 Standard results: reminders

This section is devoted to recalling some basic results that we will use in the remainder, beginning with the standard framework for the Navier–Stokes equations. Let

$$ (3.6) \qquad \mathcal{D}(\Omega) = \{u : \Omega \to \mathbb{R} \in C^\infty(\Omega) \text{ with compact support}\}, $$

and $W^{m,p}(\Omega)$ denotes the standard Sobolev spaces. We denote the $W^{m,p}$ norm (see in [10, 54]) by

$$ (3.7) \qquad \|u\|_{W^{m,p}(\Omega)} = \|u\|_{\Omega; m, p}. $$

Functions in $W^{1,p}(\Omega)$ have a trace at $\Gamma$ in the space $W^{1-\frac{1}{p},p}(\Gamma)$. The Green’s formula (1.1) yields the Stokes formula, also known the integration by parts formula,

$$ (3.8) \qquad \forall u \in W^{1,p}(\Omega), v \in W^{1,p'}(\Omega), \forall i = 1, \dots, N, \quad \int_{\Omega} u \partial_i v = \int_{\Gamma} u v n_i - \int_{\Omega} v \partial_i u, $$

where $\frac{1}{p} + \frac{1}{p'} = 1$, and $n = (n_1, \dots, n_N)$ is the outward normal vector at $\Gamma$.

We denote by $W_0^{1,p}(\Omega)$ the closure of $\mathcal{D}(\Omega)$ in $W^{1,p}(\Omega)$, and we have

$$ W_0^{1,p}(\Omega) = \{u \in W^{1,p}(\Omega) \text{ s.t. } u|_{\Gamma} = 0\}. $$

The Poincaré’s inequality reads

$$ \forall u \in W_0^{1,p}(\Omega), \quad \|u\|_{\Omega;0,p} \le C_P \|\nabla u\|_{\Omega;0,p}, $$

therefore, $u \to \|\nabla u\|_{\Omega;0,p}$ is a norm on $W_0^{1,p}(\Omega)$ equivalent to the $W^{m,p}$ norm given by $u \to \|u\|_{\Omega;0,p} + \|\nabla u\|_{\Omega;0,p}$

Let us consider the following functional spaces:

$$ \begin{align} \mathcal{V} &= \{v \in \mathcal{D}(\Omega)^N : \operatorname{div} v = 0\}, \tag{3.9} \\ \mathcal{V} &= \text{the closure of } \mathcal{V} \text{ in } H_0^1(\Omega)^N, \nonumber \\ \mathcal{H} &= \text{the closure of } \mathcal{V} \text{ in } L^2(\Omega)^N. \nonumber \end{align} $$
**Remark 3.1.** It can be proved that

$$ (3.10) \qquad H = \{ \mathbf{u} \in L^2(\Omega)^N : \underbrace{\mathrm{div} \mathbf{u} = 0}_{\text{in } \mathcal{D}'(\Omega)} , \underbrace{\mathbf{u} \cdot \mathbf{n}|_{\Gamma} = 0}_{\text{in } H^{-1/2}(\Gamma)} \}. $$

See, for instance, Girault–Raviart [28].

**Remark 3.2.** We have

$$ (3.11) \qquad V \subset H \subset V', $$

each embedding being continuous and dense. Since $H_0^1(\Omega) \subset L^2(\Omega)$ with dense, continuous and compact embedding, we naturally have $V \subset H$, and by construction of these spaces the embedding is continuous, compact and dense. Therefore, we also have $H' \subset V'$, because both $V$ and $H$ are Hilbert spaces. Thus, (3.11) holds, identifying $H$ with $H'$ via the Riesz theorem (which we do not do for $V$ and $V'$).

**Theorem 3.1** (De Rham). Let $F \in H^{-1}(\Omega)^N$ such that $\forall \mathbf{v} \in V$, $\langle F, \mathbf{v} \rangle = 0$. Then there exists a function $p \in L^2(\Omega)$ such that

$$ (3.12) \qquad F = \nabla p. $$

**Theorem 3.2.** We have

$$ (3.13) \qquad V = \{ \mathbf{v} \in H_0^1(\Omega)^N, \mathrm{div} \mathbf{v} = 0 \}. $$

The proofs of Theorems 3.1 and 3.2 can be found in [28]. We now state the Stampacchia theorem, which is essential in the analysis of the equation for $k$, which is a parabolic equation with r.h.s in $L^1$.

**Theorem 3.3** (Stampacchia). Let $G : \mathbb{R} \to \mathbb{R}$ be a Lipschitz function such that $G'$ has a finite number of discontinuities, with $G(0) = 0$. Let $u \in H_0^1(\Omega)$. Then $G(u) \in H_0^1(\Omega)$ and the following holds:

$$ (3.14) \qquad \nabla(G(u)) = G'(u)\nabla u \text{ in } \mathcal{D}'(\Omega). $$

**Remark 3.3.** A Lipschitz function has a derivative defined almost everywhere, so there is no ambiguity in the statement of the theorem.

The original proof of Stampacchia's theorem can be found in [52]. We provide a very easy proof in Appendix B when the discontinuities of $G'$ are of the first kind, which is sufficient for our purposes.

Let $\mathcal{S}$ denote the set of all Lipschitz functions $G : \mathbb{R} \to \mathbb{R}$ such that $G'$ has a finite number of discontinuities.

To conclude this section, we recall the definition of the Bochner spaces. Let $E$ be a Banach space, $T > 0$, $1 \le p \le \infty$, and let $L^p([0, T], E)$ denotes the set of measurable functions $u : [0, T] \to E$ such that the function $t \to \|u(t)\|_E$ is in $L^p([0, T])$. This is a Banch space, equipped with the norm, when $1 \le p < \infty$,

$$ \|u\|_{L^p([0,T],E)} = \left( \int_0^T \|u(t)\|_E dt \right)^{\frac{1}{p}}, $$

and when $p = \infty$,

$$ \|u\|_{L^\infty([0,T],E)} = \operatorname{esssup}_{t \in [0,T]} \|u(t)\|_E, $$

and when $E = W^{m,q}(\Omega)$, we put

$$ \|u\|_{L^p([0,T], W^{m,q}(\Omega))} = \|u\|_{p; \Omega; m, q}. $$

**Remark 3.4.** Let $1 < p < \infty$, $\frac{1}{p} + \frac{1}{p'} = 1$, then $(L^p([0, T], E))' = L^{p'}([0, T], E')$.
### 3.3 Basic estimates for the velocity

Throughout the following, we assume that assumptions (3.3), (3.4), and (3.5) hold. Recall that $Q_T = [0, T] \times \Omega$, and for a given $t \in [0, T]$, we set $Q_t = [0, t] \times \Omega$.

The first task is to derive a priori estimates from formal energy balances for strong solutions to the NS-TKE model (3.1), in order to identify the appropriate functional spaces in which to seek possible solutions later on, starting with the velocity. Strong solutions means solutions that are smooth enough to verify the system in the classical sense. Recall that the equations for the velocity are:

$$ (3.15) \qquad \left\{ \begin{array}{ll} \partial_t \mathbf{v} + \operatorname{div}(\mathbf{v} \otimes \mathbf{v}) - \operatorname{div}(\nu_{\text{turb}}(k)\nabla \mathbf{v}) + \nabla p = \mathbf{f} & \text{in } Q_T, \\ \operatorname{div} \mathbf{v} = 0 & \text{in } Q_T. \end{array} \right. $$

**Proposition 3.1.** Let $\mathbf{v} \in C^2(Q_T)$, $p \in C^1(Q_T)$ and $k \in C^0(Q_T)$, such that $\mathbf{v}|_{\Gamma} = 0$ and $(\mathbf{v}, p, k)$ satisfies (3.15). Then, $\forall t \in [0, T]$,

$$ (3.16) \qquad \int_{\Omega} |\mathbf{v}(t, x)|^2 dx + \nu \iint_{Q_t} |\nabla \mathbf{v}(s, x)|^2 dxds \le \frac{C_p}{\nu} \|\mathbf{f}\|_{Q_{T;0,2}}^2 + \int_{\Omega} |\mathbf{v}_0(x)|^2 dx, $$

where $\mathbf{v}_0(x) = \mathbf{v}(0, x)$, and $C_p$ is Poincaré's constant.

*Proof.* We take the scalar product of the fluid equation with $\mathbf{v}$ and integrate by parts, as we already did in section 1.3.2 for the Navier-Stokes equations. The difference here lies in the term $-\operatorname{div}(\nu_{\text{turb}}(k)\nabla \mathbf{v})$, in which by Stokes Formula (3.8),

$$ \int_{\Omega} -\operatorname{div}(\nu_{\text{turb}}(k)\nabla \mathbf{v}) \cdot \mathbf{v} = - \underbrace{\int_{\Gamma} \nu_{\text{turb}}(k) \frac{\partial \mathbf{v}}{\partial \mathbf{n}} \cdot \mathbf{v}}_{=0 \text{ because } \mathbf{v}|_{\Gamma}=0} + \int_{\Omega} \nu_{\text{turb}}(k) \nabla \mathbf{v} : \nabla \mathbf{v} = \int_{\Omega} \nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2. $$

Since $0 < \nu \le \nu_{\text{turb}}(k)$ for all $t \in [0, T]$, after integrating in time over $[0, t]$, we obtain

$$ \nu \iint_{Q_t} |\nabla \mathbf{v}|^2 \le \iint_{Q_t} \nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2. $$

Thus,

$$ (3.17) \qquad \frac{1}{2} \int_{\Omega} |\mathbf{v}(t, x)|^2 + \nu \iint_{Q_t} |\nabla \mathbf{v}|^2 \le \iint_{Q_t} \mathbf{f} \cdot \mathbf{v} + \frac{1}{2} \int_{\Omega} |\mathbf{v}_0|^2. $$

Combining the Cauchy-Schwarz and Young inequalities, for a given $\eta > 0$ that will be fixed later,

$$ \left| \iint_{Q_t} \mathbf{f} \cdot \mathbf{v} \right| \le \|\mathbf{f}\|_{Q_{T;0,2}} \|\mathbf{v}\|_{Q_{T;0,2}} \le \frac{1}{2\eta} \|\mathbf{f}\|_{Q_{T;0,2}}^2 + \frac{\eta}{2} \|\mathbf{v}\|_{Q_{T;0,2}}^2, $$

and by Poincaré's inequality, since $\mathbf{v}|_{\Gamma} = 0$,

$$ \|\mathbf{v}\|_{Q_{T;0,2}}^2 = \int_0^t \int_{\Omega} |\mathbf{v}(s, x)|^2 dxds \le C_p \int_0^t \int_{\Omega} |\nabla \mathbf{v}(s, x)|^2 dxds. $$

Choosing

$$ \eta = \frac{\nu}{C_p}, $$

we finally obtain (3.16).

□
**Remark 3.5.** This result implies:

$$ (3.18) \qquad \text{The function } t \to \int_{\Omega} |\mathbf{v}(t, x)|^2 dx \text{ is bounded on } [0, T], $$

in other words, $\mathbf{v} \in L^\infty([0, T], L^2(\Omega))$, and

$$ (3.19) \qquad \text{The function } t \to \iint_{Q_t} |\nabla \mathbf{v}(s, x)|^2 dx ds \text{ is bounded on } [0, T], $$

which also means that $\mathbf{v} \in L^2([0, T], V)$. Moreover, we obtain $|\nabla \mathbf{v}|^2 \in L^1(Q_T)$. Thus, the equation for the TKE in system (3.1) is a transport-diffusion equation with a right-hand side in $L^1$.

**Remark 3.6.** Since

$$ H_0^1(\Omega) \hookrightarrow L^p(\Omega) \quad \text{for} \quad \begin{cases} p < \infty & \text{if } N = 2, \\ p \le 6 & \text{if } N = 3, \end{cases} $$

it follows from the bound in $L^2([0, T], V) \cap L^\infty([0, T], L^2(\Omega))$ that

$$ (3.20) \qquad \mathbf{v} \in L^2([0, T], L^p(\Omega)^N) \quad \text{for} \quad \begin{cases} p < \infty & \text{if } N = 2, \\ p \le 6 & \text{if } N = 3. \end{cases} $$

By combining (3.20) with the $L^\infty([0, T], L^2(\Omega)^N)$ estimate and applying Lemma A.4 and formula (A.21) with $\zeta(2, 6, 2) = 10/3$, we get

$$ (3.21) \qquad \mathbf{v} \in L^{\frac{10}{3}}(Q_T) \quad \text{if } N = 3. $$

The case $N = 2$ is a little bit better, thanks to the Ladyzenskaya inequality (see a proof in [56], Lemma 3.3):

$$ \forall v \in H_0^1(\Omega), \quad \|v\|_{\Omega;0,4} \le 2^{\frac{1}{4}} \|v\|_{\Omega;0,2}^{\frac{1}{2}} \|\nabla v\|_{\Omega;0,2}^{\frac{1}{2}}, $$

which yields by the bound in $L^2([0, T], V) \cap L^\infty([0, T], L^2(\Omega))$,

$$ (3.22) \qquad \mathbf{v} \in L^4(Q_T) \quad \text{if } N = 2. $$

**Remark 3.7.** We also obtain estimates for the pressure, following the classical analysis for the Navier-Stokes equations (see [1, 13, 22, 28]), since the eddy viscosity terms do not modify the usual arguments. Therefore,

$$ (3.23) \qquad p \in \begin{cases} L^2(Q_T), & \text{if } N = 2, \\ L^{5/3}(Q_T), & \text{if } N = 3. \end{cases} $$

Estimates for the pressure can also be deduced from the Nečas inequalities (see in [9, 42]).

## 3.4 Estimates for the TKE

We will use the Stokes formula (3.8) in the remainder, without mention it systematically. Recall that the equation for the TKE in the NS-TKE model (3.1) is as follows:

$$ (3.24) \qquad \partial_t k + \mathbf{v} \cdot \nabla k - \operatorname{div}(\mu_{\text{turb}}(k) \nabla k) = \nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2 - \frac{k|k|^{\frac{1}{2}}}{\ell}. $$

We prove the following result.
**Proposition 3.2.** Let $\mathbf{v} \in C^1(Q_T)$ such that $\mathrm{div}\mathbf{v} = 0$ and $\mathbf{v}|_{\Gamma} = 0$. Let $k \in C^2(\overline{Q_T})$ such that $k|_{\Gamma} = 0$. Let $k_0(x) = k(0, t)$. Assume that $k$ satisfies (3.24). Then there exists a constant

$$C = C(\beta, \|k_0\|_{\Omega;0,1}, \|\nabla \mathbf{v}\|_{Q_r;0,2})$$

such that

$$(3.25) \qquad \forall t \in [0, T], \quad \int_{\Omega} |k(t, x)|dx \le C,$$

in other words $k \in L^{\infty}([0, T], L^1(\Omega))$.

**Remark 3.8.** The idea of the proof is to use test functions of the form $H(k)$ in equation (3.24), for suitably chosen functions $H \in \mathbb{S}$, in order to apply Stampacchia's theorem. Note that we have chosen $k \in C^2(\bar{\Omega})$ with $k|_{\Omega} = 0$ in the statement of Proposition 3.2 for convenience, to ensure that there is no ambiguity regarding the differential terms such as $\partial_t k$, $\nabla k$, etc., and to work with strong solutions. However, everything carried out below still holds for weaker solutions.

*Proof.* Let $\delta > 0$, and let $G_{\delta} \in \mathbb{S}$ be defined by

$$(3.26) \qquad G_{\delta} : u \mapsto \begin{cases} -1 & \text{if } u < -\delta, \\ \frac{u}{\delta} & \text{if } -\delta \le u \le \delta, \\ 1 & \text{if } u > \delta. \end{cases}$$

<figure><img src="image_3.png" /><figcaption>Figure 1: Graph of the function $G_{\delta}$.</figcaption></figure>

We set

$$H_{\delta}(u) = \int_{0}^{u} G_{\delta}(v) \, dv \in \mathbb{S}.$$

We take $G_{\delta}(k)$ as a test function in the Equation (3.24) and integrate by parts. We treat each term separately. Note that $H_{\delta} \in C^1(\mathbb{R})$, and since $k \in C^2(\Omega)$, we have by the standard result,

$$\nabla k G_{\delta}(k) = \nabla H_{\delta}(k).$$

Because $\mathbf{v}|_{\Gamma} = 0$ and $\mathrm{div}\mathbf{v} = 0$, we obtain

$$(3.27) \qquad \int_{\Omega} \mathbf{v} \cdot \nabla k G_{\delta}(k) = \int_{\Omega} \mathbf{v} \cdot \nabla H_{\delta}(k) = \int_{\Gamma} H_{\delta}(k) \mathbf{v} \cdot \mathbf{n} - \int_{\Omega} H_{\delta}(k) \mathrm{div}\mathbf{v} = 0.$$

We now consider the diffusion term. Since $G_{\delta} \in \mathbb{S}$ and $G(0) = 0$, and because our assumptions on $k$ imply in particular that $k \in H_0^1(\Omega)$, we have by Stampacchia's theorem,

$$\nabla G_{\delta}(k) = G'_{\delta}(k) \nabla k,$$
which yields

$$ (3.28) \qquad \begin{aligned} \int_{\Omega} \mathrm{div}(\mu_{\mathrm{turb}}(k) \nabla k) G_{\delta}(k) &= \underbrace{\int_{\Gamma} \mu_{\mathrm{turb}}(k) \frac{\partial k}{\partial \mathbf{n}} G_{\delta}(k)}_{=0 \text{ because } k|_{\Gamma}=0} - \int_{\Omega} \mu_{\mathrm{turb}}(k) \nabla k \cdot \nabla G_{\delta}(k) \\ &= - \int_{\Omega} \mu_{\mathrm{turb}}(k) G'_{\delta}(k) |\nabla k|^2. \end{aligned} $$

Moreover,

$$ (3.29) \qquad \int_{\Omega} \partial_t k G_{\delta}(k) = \int_{\Omega} \partial_t (H_{\delta}(k)) = \frac{d}{dt} \int_{\Omega} H_{\delta}(k). $$

Therefore, combining (3.27), (3.28) and (3.32), we obtain

$$ (3.30) \qquad \frac{d}{dt} \int_{\Omega} H_{\delta}(k) + \int_{\Omega} \mu_{\mathrm{turb}}(k) G'_{\delta}(k) |\nabla k|^2 = \int_{\Omega} \nu_{\mathrm{turb}}(k) |\nabla \mathbf{v}|^2 G_{\delta}(k) - \int_{\Omega} G_{\delta}(k) \frac{k|k|^{\frac{1}{2}}}{\ell}. $$

We then integrate in time over $[0, t]$, for a fixed $t$, which gives:

$$ (3.31) \qquad \int_{\Omega} H_{\delta}(k(t, x)) + \iint_{Q_t} \mu_{\mathrm{turb}}(k) G'_{\delta}(k) |\nabla k|^2 = \int_{\Omega} H_{\delta}(k_0) + \iint_{Q_t} \nu_{\mathrm{turb}}(k) |\nabla \mathbf{v}|^2 G_{\delta}(k) - \iint_{Q_t} G_{\delta}(k) \frac{k\sqrt{|k|}}{\ell}. $$

We observe that:

- $G_{\delta}$ is non-decreasing, and $0 < \nu \le \mu_{\mathrm{turb}}(k) \le \beta$, hence

$$ 0 \le \iint_{Q_t} \mu_{\mathrm{turb}}(k) G'_{\delta}(k) |\nabla k|^2; $$

- $G_{\delta}$ is an odd function, and $0 \le \ell$, hence

$$ 0 \le \iint_{Q_t} G_{\delta}(k) \frac{k\sqrt{|k|}}{\ell}; $$

- $|G_{\delta}| \le 1$ and $0 < \nu \le \nu_{\mathrm{turb}}(k) \le \beta$, hence

$$ \iint_{Q_t} \nu_{\mathrm{turb}}(k) |\nabla \mathbf{v}|^2 G_{\delta}(k) \le \beta \|\nabla \mathbf{v}\|_{Q_{T;0,2}}^2. $$

Therefore, (3.31) yields

$$ (3.32) \qquad \int_{\Omega} H_{\delta}(k(t, x)) \le \int_{\Omega} H_{\delta}(k_0) + \beta \|\nabla \mathbf{v}\|_{Q_{T;0,2}}^2. $$

It remains to let $\delta \to 0$ in (3.32). To do this, we argue as follows.

- Let $u = u(x) : \Omega \to \mathbb{R}$ be any given function. Since $0 \le H_{\delta}$, we have

$$ H_{\delta}(u(x)) \underset{\delta \to 0}{\longrightarrow} |u(x)| \quad \text{a.e. in } \Omega. $$

Then by Fatou's lemma, $\forall t \in [0, T]$,

$$ (3.33) \qquad \int_{\Omega} |k(t, x)| \le \liminf_{\delta \to 0} \int_{\Omega} H_{\delta}(k(t, x)). $$
- It is easily checked that for all $0 < \delta \le 1/2$,

$$0 \le H_{\delta}(u) \le K(u) \quad \text{where} \quad K(u) = \begin{cases} |u| & \text{if } 1 \le |u|, \\ 1 & \text{if } |u| \le 1. \end{cases}$$

Therefore, since $k \in C^2(\overline{Q}_T)$, in particular $k_0 \in L^1(\Omega)$ (because $\Omega$ is bounded), we have

$$0 \le H_{\delta}(k_0) \le K(k_0) \in L^1(\Omega).$$

Hence, by Lebesgue's theorem,

$$\lim_{\delta \to 0} \int_{\Omega} H_{\delta}(k_0) = \int_{\Omega} |k_0|.$$

Finally, taking $\delta \to 0$ in (3.32), we find that $\forall t \in [0, T]$,

$$(3.34) \qquad \int_{\Omega} |k(t, x)| dx \le \int_{\Omega} |k_0| + \beta \|\nabla \mathbf{v}\|_{Q_T; 0,2} = C(\beta, \|k_0\|_{\Omega; 0,1}, \|\nabla \mathbf{v}\|_{Q_T; 0,2}),$$

hence (3.25).

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

The aim of the remainder of this paper is to prove the existence of a "weak solution" when $N = 2$, and of a "weak subsolution" to the full NS-TKE model (3.1) when $N = 3$. Our approach proceeds in three main steps:

1. Introduce a regularized system using truncations, and define what is meant by a "weak solution".

2. Prove the existence of a solution to the truncated system by means of the Galerkin method.

3. Pass to the limit as the truncation parameter tends to $\infty$.

# 4 The truncated NS-TKE system

The main issue in the NS-TKE model (3.1) comes from the transport terms in both equations, and from the production term in the right-hand side of the TKE equation, which can only be estimated in $L^1(Q_T)$ and not much more. Therefore, it is natural to truncate this production term.

Regarding the transport term, we will use the truncated transport principle, initially introduced in [35].

In this section, we introduce the Truncated NS-TKE model and carefully establish its weak formulation. We then prove the existence of a weak solution for this truncated model by means of the Galerkin method, under assumptions (3.3), (3.4), and (3.5).

## 4.1 Truncated NS-TKE model and weak formulation

Let $0 < \alpha$ and the truncation function $T_{\alpha} \in \mathbb{S}$ be given by

$$(4.1) \qquad T_{\alpha}: u \mapsto \begin{cases} -\alpha & \text{if } u < -\alpha, \\ u & \text{if } -\alpha \le u \le \alpha, \\ \alpha & \text{if } u > \alpha. \end{cases}$$
<figure><img src="image_5.png" /><figcaption>Figure 2: Graph of the function $T_\alpha$.</figcaption></figure>

For a given vector field $\mathbf{v} = \mathbf{v}(t, x) = (v_1, v_2, v_3)$, let

$$ (4.2) \qquad T_{\alpha} \mathbf{v} = (T_{\alpha} v_1, T_{\alpha} v_2, T_{\alpha} v_3). $$

The truncated NS-TKE system (T-NS-TKE in the remainder) of the NS-TKE is the following:

$$ (4.3) \qquad \left\{ \begin{array}{ll} \partial_t \mathbf{v} + (\mathbf{v} \cdot \nabla) T_\alpha(\mathbf{v}) - \operatorname{div}(\nu_{\text{turb}}(k) \nabla \mathbf{v}) + \nabla p = \mathbf{f} & \text{in } \Omega, \\ \operatorname{div} \mathbf{v} = 0 & \text{in } \Omega, \\ \partial_t k + \mathbf{v} \cdot \nabla T_\alpha(k) - \operatorname{div}(\mu_{\text{turb}}(k) \nabla k) = T_\alpha(\nu_{\text{turb}}(k) |D\mathbf{v}|^2) - \frac{k|k|^{\frac{1}{2}}}{\ell} & \text{in } \Omega, \\ \mathbf{v} = 0 & \text{on } \Gamma, \\ k = 0 & \text{on } \Gamma, \\ \mathbf{v}_{t=0} = \mathbf{v}_0 & \text{in } \Omega, \\ k_{t=0} = T_\alpha(k_0) = k_{\alpha,0} & \text{in } \Omega. \end{array} \right. $$

As we have already said, we aim to prove the existence of a weak solution to the T-NS-TKE system (4.3) for a given $\alpha > 0$. We start by giving a weak sense to each term in the equations, one after another, defined as continuous forms on appropriate spaces. We set, for simplicity,

$$ (4.4) \qquad E_v = L^2([0, T], V), \quad E_k = L^2([0, T], H_0^1(\Omega)), $$

where we recall that

$$ V = \{ \mathbf{v} \in H_0^1(\Omega)^N, \operatorname{div} \mathbf{v} = 0 \}. $$

Throughout this section, $\alpha > 0$ is fixed, and assumptions (3.3), (3.4), and (3.5) hold.

### 4.1.1 Truncated transport terms

Given $\mathbf{v}, \mathbf{w} \in E_v, k \in E_k$, let $B_\alpha(\mathbf{v}, \mathbf{w}) \in E'_v$ and $\tilde{B}_\alpha(\mathbf{v}, k) \in E'_k$ be given by

$$ (4.5) \qquad \left\{ \begin{array}{ll} \forall \mathbf{u} \in E_v, & \langle B_\alpha(\mathbf{v}, \mathbf{w}), \mathbf{u} \rangle = \iint_{Q_T} (\mathbf{v} \cdot \nabla) T_\alpha \mathbf{w} \cdot \mathbf{u}, \\ \\ \forall q \in E_k, & \langle \tilde{B}_\alpha(\mathbf{v}, k), q \rangle = \iint_{Q_T} q \mathbf{v} \cdot \nabla T_\alpha k. \end{array} \right. $$

**Lemma 4.1.** Let $\mathbf{v}, \mathbf{w} \in E_v, k \in E_k$, then

$$ (4.6) \qquad \left\{ \begin{array}{ll} B_{\alpha}(\mathbf{v}, \mathbf{w}) \in E'_{v} = L^2([0, T], V'), \\ \tilde{B}_{\alpha}(\mathbf{v}, k) \in E'_{k} = L^2([0, T], H^{-1}(\Omega)), \end{array} \right. $$
and

$$
(4.7) \qquad \left\{
\begin{array}{l}
\|B_{\alpha}(\mathbf{v}, \mathbf{w})\|_{2;\Omega;-1,2} \le \alpha C_P \|\mathbf{v}\|_{2;\Omega;1,2}, \\
\|\tilde{B}_{\alpha}(\mathbf{v}, k)\|_{2;\Omega;-1,2} \le \alpha C_P \|\mathbf{v}\|_{2;\Omega;1,2},
\end{array}
\right.
$$

where $C_P$ denotes the Poincaré's constant.

Proof. Let $\mathbf{u} \in E_v$. By the Stokes formula,

$$
\int_{\Omega} (\mathbf{v} \cdot \nabla) T_{\alpha} \mathbf{w} \cdot \mathbf{u} = \int_{\Omega} v_j \partial_j T_{\alpha} w_i \cdot u_i = \\
\int_{\Gamma} \mathbf{v} \cdot \mathbf{n} T_{\alpha} w_i \cdot u_i - \int_{\Omega} T_{\alpha} w_i \partial_j (v_j u_i) = - \int_{\Omega} (\mathbf{v} \cdot \nabla) \mathbf{u} \cdot T_{\alpha} \mathbf{w},
$$

by using $\mathbf{v}|_{\Gamma} = 0$ and $\mathrm{div}\mathbf{v} = 0$. Then

$$
\langle B_{\alpha}(\mathbf{v}, \mathbf{w}), \mathbf{u} \rangle = - \iint_{Q_T} (\mathbf{v} \cdot \nabla) \mathbf{u} \cdot T_{\alpha} \mathbf{w},
$$

hence, by Cauchy-Schwarz inequality,

$$
|\langle B_{\alpha}(\mathbf{v}, \mathbf{w}), \mathbf{u} \rangle| \le \underbrace{\|T_{\alpha}\mathbf{w}\|_{\infty;\Omega;0,\infty}}_{\le \alpha} \|\mathbf{v}\|_{2;\Omega;0,2} \underbrace{\|\nabla\mathbf{u}\|_{2;\Omega;0,2}}_{=\|\mathbf{u}\|_{2;\Omega;1,2}}
$$

hence the result by the Poincaré's inequality. The second inequality is proved in the same
way. $\square$

**Remark 4.1.** Notice that because $\mathrm{div}\mathbf{v}=0$,

$$
(4.8) \quad (\mathbf{v} \cdot \nabla) T_{\alpha} \mathbf{v} = \mathrm{div}(T_{\alpha} \mathbf{v} \otimes \mathbf{v}) \quad \text{and} \quad \mathbf{v} \cdot \nabla T_{\alpha} k = \mathrm{div}(\mathbf{v} T_{\alpha} k).
$$

We now prove that the truncated transport ”does not work”, in the sense that it does not
contribute to the energy balance, which ensures the consistency of the model.

**Lemma 4.2.** Let $\alpha > 0$ be fixed. Then

$$
(4.9) \qquad \forall \mathbf{v}, \mathbf{w} \in E_v, \langle B_\alpha(\mathbf{v}, \mathbf{w}), \mathbf{w} \rangle = 0,
$$

and for all $G \in \mathbb{S}$,

$$
(4.10) \qquad \forall \mathbf{v} \in E_v, \forall k \in E_k, \quad \langle \tilde{B}_{\alpha}(\mathbf{v}, k), G(k) \rangle = 0.
$$

Proof. It can be shown that $C^{\infty}([0, T], \mathcal{V})$ is dense in $E_v = L^2([0, T], \mathcal{V})$. We first prove the result for $\mathbf{w} \in C^{\infty}([0, T], \mathcal{V})$; the general case $\mathbf{w} \in E_v$ then follows by density and the continuity of $B_{\alpha}$.

Therefore, let $\mathbf{v} \in E_v$, $\mathbf{w} \in C^\infty([0, T], \mathcal{V})$. Observe that for all $t \in [0, T]$, as for $i = 1, 2, 3$, the function $x \mapsto w_i(t, x)$ is continuous, the sets

$$
A_i = \{x \in \Omega, |w_i(t, x)| > \alpha\}
$$

are open sets for any given $t \in [0, T]$. Therefore, as $T_{\alpha}w_i = \alpha$ over $A_i$, and we have

$$
\forall j = 1, 2, 3, \quad \partial_j T_\alpha w_i = 0 \text{ over } A_i.
$$

By writing,

$$
\Omega = \{x \in \Omega, |w_i(t, x)| \le \alpha\} \cup A_i = (\{|w_i| \le \alpha\} \cup A_i,
$$
we get

$$ (4.11) \qquad \begin{aligned} \langle B_{\alpha}(\mathbf{v}, \mathbf{w}), \mathbf{w} \rangle &= \int_{0}^{T} \int_{\Omega} v_{j} \partial_{j} T_{\alpha} w_{i} \cdot w_{i} = \int_{0}^{T} \underbrace{\int_{\{|w_{i}| \le \alpha\}}}_{\text{since } \partial_{j} T_{\alpha} w_{i} = 0} v_{j} \partial_{j} T_{\alpha} w_{i} \cdot w_{i} \\ &= \int_{0}^{T} \underbrace{\int_{\Omega} v_{j} \partial_{j} T_{\alpha} w_{i} \cdot T_{\alpha} w_{i}}_{\text{on } A_{i}} \cdot \underbrace{\text{since } T_{\alpha} w_{i} = w_{i} \text{ on } \{|w_{i}| \le \alpha\}} \\ &= \int_{0}^{T} \int_{\Omega} v_{j} \partial_{j} T_{\alpha} w_{i} \cdot w_{i} \end{aligned} $$

Therefore,

$$ (4.12) \qquad \begin{aligned} \langle B_{\alpha}(\mathbf{v}, \mathbf{w}), \mathbf{w} \rangle &= \frac{1}{2} \int_{0}^{T} \int_{\Omega} \mathbf{v} \cdot \nabla |T_{\alpha} \mathbf{w}|^{2} \\ &= \int_{0}^{T} \int_{\Gamma} |T_{\alpha} \mathbf{w}|^{2} \mathbf{v} \cdot \mathbf{n} - \int_{0}^{T} \int_{\Omega} |T_{\alpha} \mathbf{w}|^{2} \mathrm{div} \mathbf{v} \\ &= 0 \end{aligned} $$

because $\mathbf{v}|_{\Gamma} = 0$ and $\mathrm{div} \mathbf{v} = 0$, hence (4.9).

Let us prove (4.10). Let $G \in \mathbb{S}$,

$$ H(u) = \int_{0}^{u} G(v) dv \in C^{1}(\mathbb{R}), $$

and $H(0) = 0$. By the same arguments as above, since $H(0) = 0$,

$$ \begin{aligned} \langle \tilde{B}_{\alpha}(\mathbf{v}, k), G(k) \rangle &= \int_{0}^{T} \int_{\Omega} v_{j} \partial_{j} T_{\alpha} k \cdot G(k) \\ &= \int_{0}^{T} \int_{\Omega} \mathbf{v} \cdot \nabla H(T_{\alpha} k) \\ &= \int_{0}^{T} \int_{\Gamma} H(T_{\alpha} k) \mathbf{v} \cdot \mathbf{n} - \int_{0}^{T} \int_{\Omega} H(T_{\alpha} k) \mathrm{div} \mathbf{v} \\ &= 0. \end{aligned} $$

□

### 4.1.2 Diffusion terms

We now turn to the diffusion terms. Their variational expressions follow from Stokes formula together with the homogeneous boundary conditions, which formally yields

$$ \begin{aligned} & - \int_{\Omega} \mathrm{div}(\nu_{\mathrm{turb}}(k) \nabla \mathbf{v}) \cdot \mathbf{u} = \int_{\Omega} \nu_{\mathrm{turb}}(k) \nabla \mathbf{v} : \nabla \mathbf{u}, \\ & - \int_{\Omega} \mathrm{div}(\mu_{\mathrm{turb}}(k) \nabla r) q = \int_{\Omega} \mu_{\mathrm{turb}}(k) \nabla r \cdot \nabla q. \end{aligned} $$

Therefore, given $\mathbf{v} \in E_v$ and $k, r \in E_k$, let $A(k, \mathbf{v}) \in E'_v$ and $\tilde{A}(k, r) \in E'_k$ be defined by

$$ (4.13) \qquad \left\{ \begin{array}{ll} \forall \mathbf{u} \in E_v, & \langle A(k, \mathbf{v}), \mathbf{u} \rangle = \iint_{Q_T} \nu_{\mathrm{turb}}(k) \nabla \mathbf{v} : \nabla \mathbf{u} = \langle -\mathrm{div}(\nu_{\mathrm{turb}}(k) \nabla \mathbf{v}), \mathbf{u} \rangle, \\ \\ \forall q \in E_k, & \langle \tilde{A}(k, r), q \rangle = \iint_{Q_T} \mu_{\mathrm{turb}}(k) \nabla r \cdot \nabla q = \langle -\mathrm{div}(\mu_{\mathrm{turb}}(k) \nabla r), q \rangle. \end{array} \right. $$
The proof of the following is straightforward, because of assumption (3.3).

**Lemma 4.3.** We have $A(k, \mathbf{v}) \in L^2([0, T], V')$, $\tilde{A}(k, k) \in L^2([0, T], H^{-1}(\Omega))$ and

$$ (4.14) \qquad \left\{ \begin{array}{l} \|A(k, \mathbf{v})\|_{2;\Omega;-1,2} \le \beta \|\mathbf{v}\|_{2;\Omega;1,2}, \\ \|\tilde{A}(k, k)\|_{2;\Omega;-1,2} \le \beta \|k\|_{2;\Omega;1,2}. \end{array} \right. $$

### 4.1.3 Production and dissipation terms in the TKE equation

Next, we focus on the production and dissipation terms. Given $(v, k) \in E_v \times E_k$, let $D_\alpha(\mathbf{v}, k) \in E'_k$ and $\varepsilon(k, r) \in E'_k$ be given by for $\forall q \in E_k$

$$ (4.15) \qquad \left\{ \begin{array}{l} \langle D_\alpha(\mathbf{v}, k), q \rangle = \iint_{Q_T} T_\alpha(\nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2) q, \\ \langle \varepsilon(k, r), q \rangle = \iint_{Q_T} \frac{k|r|^{\frac{1}{2}}q}{\ell}. \end{array} \right. $$

The proof of the following is also straightforward, by using the Sobolev inequalities.

**Lemma 4.4.** We have $D_\alpha(k, \mathbf{v})$, $\varepsilon(k, r) \in L^2([0, T], H^{-1}(\Omega))$, and

$$ (4.16) \qquad \left\{ \begin{array}{l} \|D_\alpha(k, \mathbf{v})\|_{2;\Omega;-1,2} \le C_D \alpha, \\ \|\varepsilon(k, r)\|_{2;\Omega;-1,2} \le \frac{C_\varepsilon}{\ell_0} \|k\|_{2;\Omega;1,2} \|r\|_{2;\Omega;1,2}^{\frac{1}{2}}, \end{array} \right. $$

where $C_D = C_D(Q_T) > 0$, $C_\varepsilon = C_\varepsilon(Q_T) > 0$.

### 4.1.4 Weak solutions of the truncated NS-TKE system

Throughout the rest of this section, we assume that (3.3), (3.4), and (3.5) hold, in particular $\mathbf{f} \in C([0, T], L^2(\Omega)^N)$. The reason for this time-continuity regularity of the source term will become clear shortly.

Recall that

$$ (4.17) \qquad \forall \mathbf{u} \in V, \forall p \in L^2(\Omega), \quad \underbrace{\langle \nabla p, \mathbf{u} \rangle}_{\text{in } D'(\Omega)} = - \int_{\Omega} p \operatorname{div} \mathbf{u} = 0 $$

Therefore, when we take the scalar product of the fluid equation with a vector field in $L^2([0, T], V)$, the resulting pressure term vanishes. The pressure is recovered by the De Rham theorem, Theorem 3.1. The initial data for the fluid equation will be taken in the space $H$ defined in (3.10). Recall that in this section, $\alpha > 0$ is fixed. The notion of weak solution of the T-NS-TKE is given as follows.

**Definition 4.1.** We say that $(\mathbf{v}, k) \in E_v \times E_k$ is a weak solution to the T-NS-TKE system (4.3) if

1. Continuity of trajectories:

$$ (4.18) \qquad \mathbf{v} \in C([0, T], H), \quad k \in C([0, T], L^2(\Omega)). $$

2. Time derivative regularity:

$$ (4.19) \qquad \partial_t \mathbf{v} \in L^2([0, T], V'), \quad \partial_t k \in L^2([0, T], H^{-1}(\Omega)); $$
3. Variational formulation: $\forall (\mathbf{w}, q) \in E_v \times E_k$,

$$ (4.20) \quad \langle \partial_t \mathbf{v}, \mathbf{w} \rangle + \langle B_\alpha(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle + \langle A(k, \mathbf{v}), \mathbf{w} \rangle = \iint_{Q_T} \mathbf{f}(t, x) \cdot \mathbf{w}(t, x) dx dt, $$

$$ (4.21) \quad \langle \partial_t k, q \rangle + \langle \tilde{B}_\alpha(\mathbf{v}, k), q \rangle + \langle \tilde{A}(k, k), q \rangle = \langle \mathbb{D}_\alpha(k, \mathbf{v}) - \varepsilon(k, k), q \rangle. $$

4. Initial data:

$$ (4.22) \quad \lim_{t \to 0} \|\mathbf{v}(t, \cdot) - \mathbf{v}_0\|_{\Omega; 0, 2} = 0, $$

$$ (4.23) \quad \lim_{t \to 0} \|k(t, \cdot) - k_{\alpha, 0}\|_{\Omega; 0, 2} = 0. $$

In (4.22), we mean by $\|\mathbf{v}(t, \cdot) - \mathbf{v}_0\|_{\Omega; 0, 2}$,

$$ \|\mathbf{v}(t, \cdot) - \mathbf{v}_0\|_{\Omega; 0, 2} = \left( \int_{\Omega} (\mathbf{v}(t, x) - \mathbf{v}_0(x))^2 dx \right)^{\frac{1}{2}}, $$

and the same thing in (4.23). Notice that there is no ambiguity with (4.22) and (4.23) since we are looking for $\mathbf{v}$ and $k$ in $C([0, T], H)$ and $C([0, T], L^2(\Omega))$.

### 4.1.5 Existence of a weak solution to the truncated NS-TKE model: statement and plan of the proof

We will prove the existence of a weak solution to the T-NS-TKE system (4.3) in the following.

**Theorem 4.1.** Assume that (3.3), (3.4), and (3.5) hold. Let $\mathbf{v}_0 \in H$, $k_0 \in L^1(\Omega)$, and $\alpha > 0$. Then the T-NS-TKE system (4.3) admits a weak solution.

**Remark 4.2.** We assumed that $\mathbf{f} \in C([0, T], L^2(\Omega)^N)$ in order to apply the Cauchy-Lipschitz theorem in the Galerkin method. However, by a density argument, the result still holds when $\mathbf{f} \in L^2([0, T], V')$. We will not elaborate further on this point.

The proof of Theorem 4.1 is based on the Galerkin method, proceeding in six main steps:

1. Project the problem onto finite-dimensional spaces and obtain approximate solutions by the Cauchy-Lipschitz theorem;

2. Energy estimates and compactness;

3. Extract suitable subsequences;

4. Pass to the limit as the dimension tends to infinity in the fluid equation;

5. Energy method to handle the production term in the TKE equation;

6. Pass to the limit in the TKE equation.

## 4.2 Approximate solutions and estimates

In what follows, we set

$$ (u, v) = \int_{\Omega} u(x)v(x) dx. $$
### 4.2.1 Finite dimensional spaces

We now describe the finite-dimensional spaces onto which we will project the equations, starting with the one for the velocity. By the compact operator theory applied to the Stokes operator, there exists a Hilbert basis of $V$, $\{\mathbf{u}_1, \dots, \mathbf{u}_n, \dots\}$, consisting of solutions to the eigenvalue problem on $V$, called a “special basis” (see [38]):

$$ (4.24) \qquad \left\{ \begin{array}{l} -\Delta \mathbf{u}_n + \nabla p_n = \lambda_n \mathbf{u}_n, \\ \mathrm{div} \mathbf{u}_n = 0, \\ \mathbf{u}_n|_{\Gamma} = 0, \end{array} \right. $$

where $0 < \lambda_1 \le \lambda_2 \le \dots$, and in addition $\lim_{n \to \infty} \lambda_n = +\infty$. Moreover, let

$$ \tilde{\mathbf{u}}_j = \sqrt{\lambda_j} \mathbf{u}_j. $$

Then $\{\tilde{\mathbf{u}}_1, \dots, \tilde{\mathbf{u}}_n, \dots\}$ is a Hilbert basis of $H$. We set

$$ (4.25) \qquad V_n = \mathrm{vect}\{\mathbf{u}_1, \dots, \mathbf{u}_n\} \subset V, $$

$$ (4.26) \qquad H_n = \mathrm{vect}\{\tilde{\mathbf{u}}_1, \dots, \tilde{\mathbf{u}}_n\} \subset H, $$

and we project the initial data $\mathbf{v}_0$ for the velocity onto $H_n$ as follows:

$$ (4.27) \qquad \mathbf{v}_{0,n} = \sum_{j=1}^{n} (\tilde{\mathbf{u}}_j, \mathbf{v}_0) \tilde{\mathbf{u}}_j = \sum_{j=1}^{n} \lambda_j (\mathbf{u}_j, \mathbf{v}_0) \mathbf{u}_j. $$

We now consider the spaces associated with the TKE. Let $\{\kappa_1, \dots, \kappa_n, \dots\}$ be a Hilbert basis of $H_0^1(\Omega)$ made of eigenfunctions of the $-\Delta$ operator, namely

$$ (4.28) \qquad \forall n \in \mathbb{N}^*, \kappa_n \in H_0^1(\Omega) \quad \text{and} \quad -\Delta \kappa_n = \mu_n \kappa_n, $$

where $0 < \mu_1 \le \dots \le \mu_n \le \dots$, $\lim_{n \to \infty} \mu_n = +\infty$. We also consider

$$ \{\tilde{\tau}_n\}_{n \in \mathbb{N}^*} = \{\sqrt{\mu_n} \tau_n\}_{n \in \mathbb{N}^*}, $$

which form a basis of $L^2(\Omega)$. We set

$$ (4.29) \qquad K_n = \mathrm{vect}\{\tau_1, \dots, \tau_n\} \subset H_0^1(\Omega), $$

$$ (4.30) \qquad L_n = \mathrm{vect}\{\tilde{\tau}_1, \dots, \tilde{\tau}_n\} \subset L^2(\Omega) $$

As we did for the velocity, we project the initial data $k_{\alpha,0}$ for the TKE as follows:

$$ (4.31) \qquad k_{\alpha,0,n} = \sum_{j=1}^{n} (\tilde{\tau}_j, k_{\alpha,0}) \tilde{\tau}_j = \sum_{j=1}^{n} \mu_j(\tau_j, k_{\alpha,0}) \tau_j. $$

**Remark 4.3.** We have

$$ (4.32) \qquad \lim_{n \to \infty} \mathbf{v}_{0,n} = \mathbf{v}_0 \quad \text{in } H, \quad \text{and} \quad \lim_{n \to \infty} k_{\alpha,0,n} = k_{\alpha,0} \quad \text{in } L^2(\Omega). $$

**Remark 4.4.** By a standard bootstrap argument combined with the elliptic regularity theory, we see that the $\mathbf{u}_j$'s and the $\tau_j$'s are of class $C^\infty$.
### 4.2.2 Galerkin method

We project the variational equations (4.20) and (4.21) onto $V_n$ and $K_n$, seeking a solution in $V_n$ and $K_n$, with test functions in the same finite-dimensional spaces, noting that it is sufficient to test against the elements of the basis of each space. Therefore, we consider the following variational problem: find

$$ (4.33) \qquad \begin{aligned} \mathbf{v}_n &= \mathbf{v}_n(t, x) = \sum_{j=1}^{n} g_{j,n}(t) \mathbf{u}_j(x) \in V_n, \\ k_n &= k_n(t, x) = \sum_{j=1}^{n} h_{j,n}(t) \tau_j(x) \in K_n, \end{aligned} $$

where $g_{j,n}, h_{j,n} \in C^1(\mathbb{R}_+)$ are such that, according to (4.27) and (4.31),

$$ (4.34) \qquad \forall j = 1, \dots, n, \quad g_{j,n}(0) = \lambda_j(\mathbf{u}_j, \mathbf{v}_0), \quad h_{j,n}(0) = \mu_j(\tau_j, k_{\alpha,0}), $$

and which satisfy, for all $j = 1, \dots, n$, and $\forall t \in [0, T]$,

$$ (4.35) \qquad \left\{ \begin{aligned} & g'_{j,n}(t) + \int_{\Omega} (\mathbf{v}_n(t,x) \cdot \nabla) T_{\alpha} \mathbf{v}_n(t,x) \cdot \mathbf{u}_j(x) dx \\ & \qquad + \int_{\Omega} \nu_{\text{turb}}(k_n(t,x)) \nabla \mathbf{v}_n(t,x) : \nabla \mathbf{u}_j(x) dx = \int_{\Omega} \mathbf{f}(t,x) \cdot \mathbf{u}_j(x) dx, \\ & h'_{j,n}(t) + \int_{\Omega} \mathbf{v}_n(t,x) \cdot \nabla T_{\alpha} k_n(t,x) \tau_j(x) dx \\ & \qquad + \int_{\Omega} \mu_{\text{turb}}(k_n(t,x)) \nabla k_n(t,x) \cdot \nabla \tau_j(x) dx = \\ & \qquad \int_{\Omega} T_{\alpha}(\nu_{\text{turb}}(k_n(t,x))) |\nabla \mathbf{v}_n(t,x)|^2 \tau_j(x) dx - \int_{\Omega} \frac{k_n(t,x) |k_n(t,x)|^{\frac{1}{2}} \tau_j(x)}{\ell(x) dx}. \end{aligned} \right. $$

By noting $\mathbf{v}_n(t,x) = \mathbf{v}_n(t) : [0, T] \to V_n$, and $k_n(t,x) = k_n(t) : [0, T] \to K_n$, considered as trajectories in $V_n$ and $K_n$, therefore the system (4.34)–(4.35) is a nonlinear differential system that can be written in the form

$$ (4.36) \qquad (P_n) : \left\{ \begin{aligned} \begin{pmatrix} \mathbf{v}_n(t) \\ k_n(t) \end{pmatrix}' &= F_n \begin{pmatrix} t, \begin{pmatrix} \mathbf{v}_n(t) \\ k_n(t) \end{pmatrix} \end{pmatrix}, \\ \begin{pmatrix} \mathbf{v}_n(0) \\ k_n(0) \end{pmatrix} &= \begin{pmatrix} \mathbf{v}_{0,n} \\ k_{\alpha,0,n} \end{pmatrix}. \end{aligned} \right. $$

As we have assumed that $\nu_{\text{turb}}, \mu_{\text{turb}}$ are of class $C^1$, that $\mathbf{f} \in C([0, T], L^2(\Omega)^N)$, and that the truncation is a Lipschitz function, the map

$$ F_n : \left\{ \begin{array}{rcl} [0, T] \times (V_n \times K_n) & \to & \mathbb{R}, \\ (t, \mathbf{U}) & \mapsto & F_n(t, \mathbf{U}), \end{array} \right. $$

is continuous and locally Lipschitz with respect to $\mathbf{U}$. Then, by the Cauchy-Lipschitz theorem, problem $(P_n)$ has a unique maximal solution

$$ \begin{pmatrix} \mathbf{u}_n(t) \\ k_n(t) \end{pmatrix} $$

of class $C^1$ on an interval $[0, T_n^*] \subset [0, T]$. It is at this point that one understands why we have assumed that $\mathbf{f}$ is continuous in time. Moreover, for any time $t < T_n^*$, we have $\partial_t \mathbf{u}_n, \partial_t k_n \in C(Q_t)$.
### 4.2.3 Energy estimates and compactness

We now have to prove that, for every $n, T_n^* = T$. To do so, we use the no-blow-up principle, deduced from the energy balance. Let

$$ \mathbf{w} = \mathbf{w}(t, x) = \sum_{j=1}^{n} w_j(t) \mathbf{u}_j(x), \quad q = \sum_{j=1}^{n} q_j(t) \tau_j(x), $$

such that $w_j, q_j \in L^2([0, T])$, for all $j = 1, \dots, n$, so that $\mathbf{w} \in E_v$ and $q \in E_k$. Multiplying the first equation of (4.34) by $w_j(t)$ and the second one by $q_j(t)$, summing over $j$ from 1 to $n$, and integrating with respect to time over $[0, t]$, yields

$$ (4.37) \quad \iint_{Q_t} \partial_t \mathbf{v}_n \cdot \mathbf{w} + \langle B_\alpha(\mathbf{v}_n, \mathbf{v}_n), \mathbf{w} \rangle + \langle A(k_n, \mathbf{v}_n), \mathbf{w} \rangle = \iint_{Q_t} \mathbf{f} \cdot \mathbf{w}, $$

$$ (4.38) \quad \iint_{Q_t} \partial_t k_n q + \langle \tilde{B}_\alpha(\mathbf{v}_n, k_n), q \rangle + \langle \tilde{A}(k_n, k_n), q \rangle = \langle \mathbb{D}_\alpha(k_n, \mathbf{v}_n), q \rangle - \langle \varepsilon(k_n, k_n), q \rangle, $$

where, so far no risk of confusion occurs, the duality brakets hold in $L^2([0, t], V')$ and $L^2([0, t], H^{-1}(\Omega))$. In particular, we can take $\mathbf{w} = \mathbf{v}_n$ in (4.37) and $r = k_n$ in (4.38). As $\mathbf{v}_n$ is $C^1$ in time, $C^\infty$ in space (c.f. Remark 4.4), we have

$$ \iint_{Q_t} \partial_t \mathbf{v}_n(s, x) \cdot \mathbf{v}_n(s, x) dx ds = \frac{1}{2} \int_{\Omega} |\mathbf{v}_n(t, x)|^2 dx - \frac{1}{2} \int_{\Omega} |\mathbf{v}_{0,n}(x)|^2 dx, $$

$$ \iint_{Q_t} \partial_t k_n \cdot k_n = \frac{1}{2} \int_{\Omega} |\mathbf{v}_n(t, x)|^2 dx - \frac{1}{2} \int_{\Omega} k_{0, \alpha, n}(x)^2 dx. $$

Moreover, we also have by Parseval-Bessel inequality

$$ \int_{\Omega} |\mathbf{v}_{0,n}(x)|^2 dx \leq \int_{\Omega} |\mathbf{v}_0(x)|^2 dx, \quad \int_{\Omega} k_{0,\alpha,n}(x)^2 dx \leq \int_{\Omega} k_{\alpha,0}(x)^2 dx. $$

Then, in the same way as in Section 3.3, using the results of Lemma 4.2 concerning the operators $B_\alpha$ and $\tilde{B}_\alpha$, combined with $\langle \varepsilon(k_n, k_n), k_n \rangle \geq 0$, we obtain, after applying the Cauchy-Schwarz, Young, and Poincaré inequalities together with assumption (3.3), $\forall t \in [0, T]$,

$$ (4.39) \quad \int_{\Omega} |\mathbf{v}_n(t, x)|^2 dx + \nu \iint_{Q_t} |\nabla \mathbf{v}_n(s, x)|^2 dx ds \leq \frac{C_p}{\nu} \|\mathbf{f}\|_{Q_{T;0,2}}^2 + \int_{\Omega} |\mathbf{v}_0(x)|^2 dx, $$

$$ \int_{\Omega} |k_n(t, x)|^2 dx + \nu \iint_{Q_t} |\nabla k_n(s, x)|^2 dx ds \leq \frac{\alpha^2}{\nu} + \int_{\Omega} k_{\alpha,0}(x)^2 dx, $$

by recalling that $Q_t = [0, t] \times \Omega$. Therefore,

$$ t \mapsto \left\| \begin{pmatrix} \mathbf{u}_n(t) \\ k_n(t) \end{pmatrix} \right\| $$

remains bounded on $[0, T_n^*[$, where the norm considered is the Euclidean one. By the no-blow-up criterion, we obtain $T_n^* = T$. Moreover,

$$ (4.40) \quad \left\{ \begin{array}{ll} \text{the sequence } (\mathbf{v}_n)_{n \in \mathbb{N}} \text{ is bounded in } L^\infty([0, T], L^2(\Omega)^N) & \\ \text{and in } E_v = L^2([0, T], V), \\ \text{the sequence } (k_n)_{n \in \mathbb{N}} \text{ is bounded in } L^\infty([0, T], L^2(\Omega)) & \\ \text{and in } E_k = L^2([0, T], H_0^1(\Omega)). & \end{array} \right. $$
In order to get more compactness, we aim to apply the Aubin-Lions lemma (Lemma A.1), which requires that we can obtain an estimate for $\partial_t \mathbf{v}_n$ and $\partial_t k_n$. Let $\mathbb{P}_n, \tilde{\mathbb{P}}_n$ denotes the orthogonal projectors onto $V_n$ and $K_n$. We note that

$$ \partial_t \mathbf{v}_n = \sum_{j=1}^{n} g'_{j,n}(t) \mathbf{u}_j, \quad \partial_t k_n = \sum_{j=1}^{n} h'_{j,n}(t) \tau_j. $$

Therefore, given $\mathbf{w} \in E_v, q \in E_k$,

$$ \begin{aligned} \langle \partial_t \mathbf{v}_n, \mathbf{w} \rangle &= \langle \partial_t \mathbf{v}_n, \mathbb{P}_n \mathbf{w} \rangle = \langle -B_\alpha(\mathbf{v}_n, \mathbf{v}_n) - A(k_n, \mathbf{v}_n) + \mathbf{f}, \mathbb{P}_n \mathbf{w} \rangle \\ \langle \partial_t k_n, q \rangle &= \langle \partial_t k_n, \tilde{\mathbb{P}}_n q \rangle = \langle -\tilde{B}_\alpha(\mathbf{v}_n, k_n) - \tilde{A}(k_n, k_n) + \mathbb{D}_\alpha(k_n, \mathbf{v}_n) - \varepsilon(k_n, k_n), \tilde{\mathbb{P}}_n q \rangle \end{aligned} $$

We deduce from Lemmas 4.1, 4.3, 4.4, combined with (4.40), that

$$ |\langle \partial_t \mathbf{v}_n, \mathbf{w} \rangle| \le \|\mathbb{P}_n \mathbf{w}\|_{2;\Omega;1,2} \|\mathbf{v}_n\|_{2;\Omega;1,2} (C_P(\alpha + \|\mathbf{f}\|_{Q_T;0,2}) + \beta) $$

As a consequence, we find that

$$ (4.41) \qquad \left\{ \begin{array}{l} \text{the sequence } (\partial_t \mathbf{v}_n) \text{ is bounded in } L^2([0, T], V'), \\ \text{the sequence } (\partial_t k_n) \text{ is bounded in } L^2([0, T], H^{-1}(\Omega)). \end{array} \right. $$

Therefore, by combining (4.40) and (4.41), we can apply the Aubin-Lions lemma (Lemma A.1) with the following triplets of spaces: (see Remark 3.2)

$$ (4.42) \qquad V \subset H \subset V' \quad \text{and} \quad H_0^1(\Omega) \subset L^2(\Omega) \subset H^{-1}(\Omega), $$

where the embeddings are dense, and continuous, and the first ones are compact (see Remark 3.2). Since $L^2([0, T], L^2(\Omega)) \approx L^2(Q_T)$, we obtain that

$$ (4.43) \qquad (\mathbf{v}_n, k_n)_{n \in \mathbb{N}} \text{ is} \left\{ \begin{array}{l} \text{weakly compact in } E_v \times E_k, \\ \text{strongly compact in } L^2(Q_T)^N \times L^2(Q_T). \end{array} \right. $$

Moreover, we also deduce from (4.40) that

$$ (4.44) \qquad \begin{array}{l} (\mathbf{v}_n)_{n \in \mathbb{N}} \text{ is weakly-star compact in } L^\infty([0, T], L^2(\Omega)^N) = [L^1([0, T], L^2(\Omega)^N)]', \\ (k_n)_{n \in \mathbb{N}} \text{ is weakly-star compact in } L^\infty([0, T], L^2(\Omega)) = [L^1([0, T], L^2(\Omega))]'. \end{array} $$

# 4.3 Passing to the limit in the equations

## 4.3.1 Extracting subsequences

Thanks to (4.43), from the sequence $(\mathbf{v}_n, k_n)_{n \in \mathbb{N}}$, we can extract a subsequence, still denoted by $(\mathbf{v}_n, k_n)_{n \in \mathbb{N}}$, such that

$$ (4.45) \qquad \begin{array}{r@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c} \mathbf{v}_n & \xrightarrow{n \to \infty} & \mathbf{v} & \begin{cases} \text{weakly in } E_v, \\ \text{weakly-star in } L^\infty([0, T], L^2(\Omega)^N), \\ \text{strongly in } L^2(Q_T)^N, \\ \text{almost everywhere in } Q_T, \end{cases} \\ & & & \\ k_n & \xrightarrow{n \to \infty} & k & \begin{cases} \text{weakly in } E_k, \\ \text{weakly-star in } L^\infty([0, T], L^2(\Omega)), \\ \text{strongly in } L^2(Q_T), \\ \text{almost everywhere in } Q_T. \end{cases} \end{array} $$
The almost everywhere convergences in (4.45) are obtained by means of the inverse Lebesgue theorem<fnref n="1" />.

**Remark 4.5.** Throughout the rest of the paper, we will use the following fact: since $Q_T$ is bounded, any sequence that is compact in $L^q(Q_T)$ and bounded in $L^p(Q_T)$ with $1 \le q < p$ is actually compact in every intermediate space $L^r(Q_T)$ for all $q \le r < p$.

By Remark 4.5 and the interpolation inequalities stated in Remark 3.6, we have:

$$ (4.46) \quad (\mathbf{v}_n, k_n) \to (\mathbf{v}, k) \quad \text{strongly in } L^p(Q_T)^N \times L^p(Q_T), \quad \begin{cases} \forall p < 4, & \text{if } N = 2, \\ \forall p < \frac{10}{3}, & \text{if } N = 3. \end{cases} $$

Furthermore, by the inverse Lebesgue theorem, there exists $g \ge 0$ in $L^2(Q_T)$ such that

$$ \forall n \in \mathbb{N}, \quad |\mathbf{v}_n| \le g \quad \text{a.e. in } Q_T. $$

Finally, from (4.41), the subsequence can be chosen so that there exist $G \in L^2([0, T], V')$ and $H \in L^2([0, T], H^{-1}(\Omega))$ satisfying

$$ (4.47) \quad \begin{aligned} \partial_t \mathbf{v}_n &\to G & \text{weakly in } L^2([0, T], V'), \\ \partial_t k_n &\to H & \text{weakly in } L^2([0, T], H^{-1}(\Omega)). \end{aligned} $$

It remains to check that $G = \partial_t v$ and $H = \partial_t k$, in the sense of the distribution, therefore in $L^2([0, T], V')$ and $L^2([0, T], H^{-1}(\Omega))$ respectively. To do so, let $\mathbf{w} \in \mathcal{D}(Q_T)^N$. Then we have

$$ \langle \partial_t \mathbf{v}_n, \mathbf{w} \rangle = -\langle \mathbf{v}_n, \partial_t \mathbf{w} \rangle \underset{n \to \infty}{\longrightarrow} -\langle \mathbf{v}, \partial_t \mathbf{w} \rangle = \langle \partial_t \mathbf{v}, \mathbf{w} \rangle, $$

where we have used the convergences (4.45). Therefore,

$$ (4.48) \quad G = \partial_t \mathbf{v} \text{ in } \mathcal{D}'(Q_T)^N \text{ hence in } L^2([0, T], V'). $$

With the same reasoning,

$$ (4.49) \quad H = \partial_t k \text{ in } L^2([0, T], H^{-1}(\Omega)). $$

We now need to prove that $(\mathbf{v}, k)$ is a weak solution of the T-NS-TKE system (4.3), according to Definition 4.1.

First, we observe that by combining $\mathbf{v} \in E_v = L^2([0, T], V)$ and $k \in E_k = L^2([0, T], H_0^1(\Omega))$ with (4.48) and (4.49), it follows from Theorem A.2 that

$$ (4.50) \quad \mathbf{v} \in C([0, T], H), \quad k \in C([0, T], L^2(\Omega)), $$

which corresponds to condition (4.18) in Definition 4.1 of a weak solution. Moreover, the time regularity condition (4.19) appearing in the definition is also satisfied.

According to our plan, it remains to:

* Pass to the limit as the dimension tends to infinity in the fluid equation;

* Detail the energy method to handle the production term in the TKE equation;

* Pass to the limit in the TKE equation.

<footnote n="1">¹The Lebesgue inverse Theorem says that given a sequence $(f_n)_{n \in \mathbb{N}}$ that converge to $f$ in $L^1(\Omega)$, then from this sequence one can extract a subsequence that converges a.e to $f$ in $\Omega$, and for this subsequence, there exists $F \in L^1(\Omega)$ such that $|f_n| \le F$, see in [10]. This also works in $L^p(\Omega)$.</footnote>
### 4.3.2 Passing to the limit in the fluid equation

Let $(\mathbf{w}, q) \in E_v \times E_k$ be test functions, projected onto $V_n$ and $K_n$ as

$$ (4.51) \qquad \left\{ \begin{array}{l} \mathbf{w}_n = \sum_{j=1}^{n} (\nabla \mathbf{w}, \nabla \mathbf{u}_j) \mathbf{u}_j \in L^2([0, T], V_n), \\ q_n = \sum_{j=1}^{n} (\nabla q, \nabla \tau_j) \tau_j \in L^2([0, T], K_n). \end{array} \right. $$

We obviously have

$$ (\mathbf{w}_n, q_n) \underset{n \to \infty}{\longrightarrow} (\mathbf{w}, q) \quad \text{in } E_v \times E_k. $$

We first focus on the fluid equation (4.37), passing to the limit in each term one by one. However, when the corresponding terms are handled in the same way in the TKE equation (4.38), we proceed directly in order to avoid unnecessary repetitions.

**Item 4.3.i) Evolution term.** Weak convergence combined with strong convergence yields convergence in the duality brackets, and hence in $L^2([0, T], E_v \times E_k)$:

$$ (4.52) \qquad \left( \frac{\langle \partial_t \mathbf{v}_n, \mathbf{w}_n \rangle}{\langle \partial_t k_n, q_n \rangle} \right) \underset{n \to \infty}{\longrightarrow} \left( \frac{\langle \partial_t \mathbf{v}, \mathbf{w} \rangle}{\langle \partial_t k, q \rangle} \right). $$

**Item 4.3.ii) Transport term.** We know that

$$ \langle B_{\alpha}(\mathbf{v}_n, \mathbf{v}_n), \mathbf{w}_n \rangle = - \iint_{Q_T} T_{\alpha} \mathbf{v}_n \otimes \mathbf{v}_n : \nabla \mathbf{w}_n. $$

Since $T_{\alpha}$ is continuous and $|\mathbf{v}_n| \le g \in L^2(Q_T)$ a.e. in $Q_T$, we have

$$ T_{\alpha} \mathbf{v}_n \otimes \mathbf{v}_n \to T_{\alpha} \mathbf{v} \otimes \mathbf{v} \quad \text{a.e. in } Q_T, \quad |T_{\alpha} \mathbf{v}_n \otimes \mathbf{v}_n| \le \alpha g \in L^2(Q_T). $$

Therefore, by Lebesgue's theorem,

$$ T_{\alpha} \mathbf{v}_n \otimes \mathbf{v}_n \underset{n \to \infty}{\longrightarrow} T_{\alpha} \mathbf{v} \otimes \mathbf{v} \quad \text{in } L^2(Q_T)^N. $$

As $\nabla \mathbf{w}_n \underset{n \to \infty}{\longrightarrow} \nabla \mathbf{w}$ in $L^2(Q_T)^N$, it follows that

$$ (4.53) \qquad \lim_{n \to \infty} \langle B_{\alpha}(\mathbf{v}_n, \mathbf{v}_n), \mathbf{w}_n \rangle = \iint_{Q_T} T_{\alpha} \mathbf{v} \otimes \mathbf{v} : \nabla \mathbf{w} = \langle B_{\alpha}(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle. $$

By the same reasoning, we also obtain

$$ (4.54) \qquad \lim_{n \to \infty} \langle \tilde{B}_{\alpha}(\mathbf{v}_n, k_n), q_n \rangle = \iint_{Q_T} \mathbf{v} \cdot \nabla q T_{\alpha}(k) = \langle \tilde{B}_{\alpha}(\mathbf{v}, k), q \rangle. $$

**Item 4.3.iii) Diffusion term.**

$$ \langle A(\mathbf{v}_n, k_n), \mathbf{w}_n \rangle = \iint_{Q_T} \nu_{\text{turb}}(k_n) \nabla \mathbf{v}_n : \nabla \mathbf{w}_n. $$

We start by showing the following lemma.
**Lemma 4.5.** Let

$$
(4.55) \qquad Z_n = \nu_{\text{turb}}(k_n) \nabla \mathbf{v}_n, \quad Z = \nu_{\text{turb}}(k) \nabla \mathbf{v}.
$$

Then

$$
(4.56) \qquad Z_n \xrightarrow{n \to \infty} Z \quad \text{weakly in } L^2(Q_T)^{N \times N}.
$$

Proof. Let $U \in L^2(Q_T)^{N \times N}$. We must prove that

$$
(4.57) \qquad \iint_{Q_T} Z_n(t, x) : U(t, x) dxdt \xrightarrow{n \to \infty} \iint_{Q_T} Z(t, x) : U(t, x) dxdt.
$$

We deduce from (4.45) that

$$
(4.58) \qquad \nabla \mathbf{v}_n \xrightarrow{n \to \infty} \nabla \mathbf{v} \quad \text{weakly in} \quad L^2(Q_T)^N.
$$

Let us define

$$
(4.59) \qquad W_n = \nu_{\text{turb}}(k_n)U.
$$

Since $\nu_{\text{turb}}$ is bounded and continuous and $k_n \to k$ a.e. in $Q_T$, we have

$$
W_n \xrightarrow{n \to \infty} \nu_{\text{turb}}(k)U \quad \text{a.e. in } Q_T, \quad |W_n| \le \beta|U| \in L^2(Q_T).
$$

Hence, by Lebesgue’s theorem,

$$
W_n \xrightarrow{n \to \infty} \nu_{\text{turb}}(k)U \quad \text{in } L^2(Q_T)^{N \times N}.
$$

which combined with (4.58) yields

$$
\lim_{n \to \infty} \iint_{Q_T} \nu_{\text{turb}}(k_n) U : \nabla \mathbf{v}_n = \iint_{Q_T} \nu_{\text{turb}}(k) U : \nabla \mathbf{v},
$$

proving (4.57) thanks to the definition (4.55).

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

We are now able to handle the diffusion term. Indeed, since $\nabla \mathbf{w}_n \to \nabla \mathbf{w}$ strongly in $L^2(Q_T)^{N \times N}$, it follows that

$$
(4.60) \quad \lim_{n \to \infty} \langle A(\mathbf{v}_n, k_n), \mathbf{w}_n \rangle = \iint_{Q_T} \nu_{\text{turb}}(k) \nabla \mathbf{v} : \nabla \mathbf{w} = \langle A(\mathbf{v}, k), \mathbf{w} \rangle.
$$

Exactly in the same way,

$$
(4.61) \quad \lim_{n \to \infty} \langle \tilde{A}(k_n, k_n), q_n \rangle = \iint_{Q_T} \nu_{\text{turb}}(k) \nabla k \cdot \nabla r = \langle \tilde{A}(k, k), q \rangle.
$$

**Item 4.3.iv) Source term and partial conclusion.** Since

$$
\lim_{n \to \infty} \int_{Q_T} \mathbf{f} \cdot \mathbf{w}_n = \int_{Q_T} \mathbf{f} \cdot \mathbf{w},
$$

by combining (4.47), (4.48), (4.52), and (4.60), we obtain

$$
(4.62) \quad \forall \mathbf{w} \in E_v, \quad \langle \partial_t \mathbf{v}, \mathbf{w} \rangle + \langle B_\alpha(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle + \langle A(\mathbf{v}, k), \mathbf{w} \rangle = \int_{Q_T} \mathbf{f}(t, x) \cdot \mathbf{w}(t, x) dxdt,
$$
which proves that $(\mathbf{v}, k)$ satisfies (4.20) in Definition 4.1.

**Item 4.3.v)** *Initial data.* It remains to verify the initial data, namely (4.22). Since $\mathbf{v} \in C([0, T], H)$, $\mathbf{v}(0, x) = \mathbf{V}_0(x) \in H$ and

$$ \lim_{t \to 0} \|\mathbf{v}(t, \cdot) - \mathbf{V}_0\|_{\Omega;0,2} = 0. $$

We must prove that $\mathbf{v}_0 = \mathbf{V}_0$. We return back to (4.27), (4.33) and (4.34). Recall that

$$ \mathbf{v}_n = \sum_{j=1}^{n} g_{j,n}(t) \mathbf{u}_j, \quad \mathbf{v}'_n = \partial_t \mathbf{v}_n = \sum_{j=1}^{n} g_{j,n}(t) \mathbf{u}_j, $$

so that the fluid part of (4.35), for a given $j \le n$, yields

$$ (4.63) \quad \forall \mathbf{w} \in V_j, \quad (\partial_t \mathbf{v}_n, \mathbf{w}) + ((\mathbf{v}_n \cdot \nabla)T_\alpha(\mathbf{v}_n), \mathbf{w}) + (\nu_{\text{turb}}(k_n)\nabla\mathbf{v}_n, \mathbf{w}) = (\mathbf{f}, \mathbf{w}). $$

However, as $\mathbf{v}_n \in C^1(Q_T)$, then for a given $\mathbf{w} \in V_j \subset C^0(\Omega)$, we have

$$ (\partial_t \mathbf{v}_n, \mathbf{w}) = \frac{d}{dt}(\mathbf{v}_n, \mathbf{w}). $$

Therefore, (4.63) becomes

$$ (4.64) \quad \frac{d}{dt}(\mathbf{v}_n, \mathbf{w}) = (\mathbf{f}, \mathbf{w}) - ((\mathbf{v}_n \cdot \nabla)T_\alpha(\mathbf{v}_n), \mathbf{w}) - (\nu_{\text{turb}}(k_n)\nabla\mathbf{v}_n, \mathbf{w}). $$

Let $\Phi \in C^\infty([0, T])$ be such that $\Phi(T) = 0$. We multiply (4.64) by $\Phi$ and integrate by parts in time over $[0, T]$:

$$ (4.65) \quad \int_{\Omega} \mathbf{v}_{n,0}(x) \cdot \mathbf{w}(x) \Phi(0) dx + \iint_{Q_T} \mathbf{v}_n(t,x) \cdot \mathbf{w}(x) \Phi'(t) dx dt = \\ \langle B_\alpha(\mathbf{v}_n, \mathbf{v}_n), \mathbf{w}\Phi \rangle + \langle A(\mathbf{v}_n, k_n), \mathbf{w}\Phi \rangle - \iint_{Q_T} \mathbf{f}(t,x) \cdot \mathbf{w}(x) \Phi(t) dx dt. $$

Since $\Phi \mathbf{w} \in E_v$, letting $n \to \infty$ and using (4.32), (4.45), (4.60) and (4.53), we obtain

$$ (4.66) \quad \int_{\Omega} \mathbf{v}_0(x) \cdot \mathbf{w}(x) \Phi(0) dx + \iint_{Q_T} \mathbf{v}(t,x) \cdot \mathbf{w}(x) \Phi'(t) dx dt = \\ \langle B_\alpha(\mathbf{v}, \mathbf{v}), \mathbf{w}\Phi \rangle + \langle A(\mathbf{v}, k), \mathbf{w}\Phi \rangle - \iint_{Q_T} \mathbf{f}(t,x) \cdot \mathbf{w}(x) \Phi(t) dx dt. $$

This holds for all $\mathbf{w} \in V_j$ and all $j \in \mathbb{N}$, hence for all $\mathbf{w} \in V$. Taking $\Phi = \psi \in D([0, T])$ in (4.66), we infer by Theorem A.1 that

$$ (4.67) \quad \forall \mathbf{w} \in V, \quad \frac{d}{dt}(\mathbf{v}, \mathbf{w}) = (\mathbf{f}, \mathbf{w}) - ((\mathbf{v} \cdot \nabla)T_\alpha(\mathbf{v}), \mathbf{w}) - (\nu_{\text{turb}}(k_n)\nabla\mathbf{v}, \mathbf{w}) \quad \text{in } \mathcal{D}'([0, T]). $$

As $\mathbf{v} \in C([0, T], H)$, for a given $\mathbf{w} \in V$, the map $t \mapsto (\mathbf{v}(t, \cdot), \mathbf{w}(\cdot))$ is continuous. Let $C^\infty([0, T])$ be equipped with the topology given by the seminorms $\|\psi^{(n)}\|_\infty$, $n \in \mathbb{N}$. Because of the continuity of $t \mapsto (\mathbf{v}(t, \cdot), \mathbf{w}(\cdot))$, we can consider $\frac{d}{dt}(\mathbf{v}, \mathbf{w})$ as an element of $C^\infty([0, T])'$, defined by

$$ \langle \frac{d}{dt}(\mathbf{v}, \mathbf{w}), \psi \rangle = (\mathbf{v}(T, \cdot), \mathbf{w}(\cdot))\psi(T) - \underbrace{(\mathbf{v}(0, \cdot), \mathbf{w}(\cdot))}_{(\mathbf{v}_0, \mathbf{w})} \psi(0) - \int_0^T (\mathbf{v}(t, \cdot), \mathbf{w}(\cdot))\psi'(t) dt. $$
Therefore, taking $\psi = \Phi$ in (4.67) with $\Phi(T) = 0$, we obtain

$$ (4.68) \qquad \int_{\Omega} \mathbf{V}_0(x) \cdot \mathbf{w}(x)\Phi(0) \, dx + \iint_{Q_T} \mathbf{v}(t,x) \cdot \mathbf{w}(x)\Phi'(t) \, dx \, dt = \\ \qquad \langle B_{\alpha}(\mathbf{v}, \mathbf{v}), \mathbf{w}\Phi \rangle + \langle A(\mathbf{v}, k), \mathbf{w}\Phi \rangle - \iint_{Q_T} \mathbf{f}(t,x) \cdot \mathbf{w}(x)\Phi(t) \, dx \, dt. $$

Comparing (4.66) and (4.68), and choosing $\Phi$ with $\Phi(0) = 1$, we obtain

$$ \forall \mathbf{w} \in V, \quad (\mathbf{V}_0 - \mathbf{v}_0, \mathbf{w}) = 0. $$

Since $V$ is dense in $H$, we deduce that $\mathbf{V}_0 = \mathbf{v}_0$, which completes the proof.

**Remark 4.6.** Let $\mathbb{P}_n$ denotes the orthogonal projector onto $H_n$, and let $\mathbf{w} \in C^1([0, T], V)$ such that $w(T, \cdot) = 0$. It is easily checked that $\mathbb{P}_n \partial_t \mathbf{w} = \partial_t \mathbb{P}_n \mathbf{w}$. Since,

$$ \langle \partial_t \mathbf{v}_n, \mathbf{w} \rangle = \langle \partial_t \mathbf{v}_n, \mathbb{P}_n \mathbf{w} \rangle, $$

we have

$$ \begin{aligned} \langle \partial_t \mathbf{v}_n, \mathbf{w} \rangle &= - \int_0^T \int_{\Omega} \partial_t \mathbb{P}_n \mathbf{w}(t,x) \cdot \mathbf{v}_n(t,x) dx dt - \int_{\Omega} \mathbb{P}_n \mathbf{w}(0,x) \mathbf{v}_{n,0}(x) dx \\ &= - \int_0^T \int_{\Omega} \mathbb{P}_n \partial_t \mathbf{w}(t,x) \cdot \mathbf{v}_n(t,x) dx dt - \int_{\Omega} \mathbb{P}_n \mathbf{w}(0,x) \mathbf{v}_{n,0}(x) dx \\ &= - \int_0^T \int_{\Omega} \partial_t \mathbf{w}(t,x) \cdot \mathbf{v}_n(t,x) dx dt - \int_{\Omega} \mathbf{w}(0,x) \mathbf{v}_{n,0}(x) dx. \end{aligned} $$

Therefore, we obtain by passing to the limit

$$ \langle \partial_t \mathbf{v}_n, \mathbf{w} \rangle \xrightarrow[n \to \infty]{} - \int_0^T \int_{\Omega} \partial_t \mathbf{w}(t,x) \cdot \mathbf{v}(t,x) dx dt - \int_{\Omega} \mathbf{w}(0,x) \mathbf{v}_0(x) dx. $$

Let

$$ \mathcal{E}_v = \{ \mathbf{w} \in C^1([0, T], V) \text{ s.t. } \forall x \in \Omega, \mathbf{w}(T, x) = 0 \}. $$

Then $(\mathbf{v}, k)$ also satisfies: $\forall \mathbf{w} \in \mathcal{E}_v$,

$$ (4.69) \qquad - \iint_{Q_T} \partial_t \mathbf{w}(t,x) \cdot \mathbf{v}(t,x) dx dt - \int_{\Omega} \mathbf{w}(0,x) \mathbf{v}_0(x) dx + \\ \qquad \langle B_{\alpha}(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle + \langle A(\mathbf{v}, k), \mathbf{w} \rangle = \iint_{Q_T} \mathbf{f}(t,x) \cdot \mathbf{w}(t,x) dx dt. $$

### 4.3.3 The energy method

The main difficulty of our problem lies in passing to the limit in the production term $D_{\alpha}(k_n, \mathbf{v}_n)$, because the sequence $(\nu_{\text{turb}}(k_n) |\nabla \mathbf{v}_n|^2)_{n \in \mathbb{N}}$ is only bounded in $L^1(Q_T)$. To handle this term, we use what we call the energy method, based on energy balances. Since this method will be used throughout the rest of the paper, we present it carefully in this section, so that we will not need to repeat these details later. Recall that

$$ (4.70) \qquad Z_n = \nu_{\text{turb}}(k_n) \nabla \mathbf{v}_n, \quad Z = \nu_{\text{turb}}(k) \nabla \mathbf{v}. $$

We already know by Lemma 4.5 that the sequence $(Z_n)_{n \in \mathbb{N}}$ converges weakly to $Z$ in $L^2(Q_T)^N$. The main result of this section is the following.
**Proposition 4.1.** For all $0 < T' < T$, we have

$$Z_n \xrightarrow{n \to \infty} Z \quad \text{strongly in } L^2(Q_T)^{N \times N}.$$

**Remark 4.7.** The strong convergence holds only on a time interval $[0, T']$ with $T' < T$, for reasons that will be made clear later in the proof. However, it is sufficient to assume that

$$\mathbf{f} \in C([0, T + \delta], L^2(\Omega))$$

for some $\delta > 0$. Then, the approximation can be constructed over $[0, T + \delta]$ with similar estimates, which allows us to obtain the strong convergence result over $[0, T]$. We will therefore take this extension for granted once the proof is complete.

*Proof.* As $Z_n \xrightarrow{n \to \infty} Z$ weakly in $L^2(Q_T)^{N \times N}$, we would like to prove that

$$(4.71) \qquad \forall 0 < T' < T, \quad \|Z_n\|_{Q_{T'}}; 0,2} \xrightarrow{n \to \infty} \|Z\|_{Q_{T'}}; 0,2,$$

which is not exactly what we will do, but almost. The key point of the proof is that we can take $\mathbf{w} = \mathbf{v} \in E_v$ as a test vector field in Equation (4.62). Recall that $B_\alpha(\mathbf{v}, \mathbf{v}, \mathbf{v}) = 0$, and since $\mathbf{v} \in E_v$, it can indeed be taken as a test vector field. Moreover, we also know that $\partial_t \mathbf{v} \in L^2([0, T], V')$. We are therefore in a position to apply Theorem A.2 in Section A.2, and according to Formula (A.6), we have for any $t \in [0, T]$,

$$(4.72) \qquad \langle \partial_t \mathbf{v}, \mathbf{v} \rangle = \frac{1}{2} \int_{\Omega} |\mathbf{v}(t, x)|^2 dx - \frac{1}{2} \int_{\Omega} |\mathbf{v}_0|^2 dx,$$

where the duality brackets hold between $L^2([0, t], V')$ and $L^2([0, t], V)$. Therefore, the following energy balance holds for any time $t$:

$$(4.73) \qquad \frac{1}{2} \int_{\Omega} |\mathbf{v}(t, x)|^2 dx + \iint_{Q_t} \nu_{\text{turb}}(k(s, x)) |\nabla \mathbf{v}(s, x)|^2 dx ds = \\ \frac{1}{2} \int_{\Omega} |\mathbf{v}_0(x)|^2 dx + \iint_{Q_t} \mathbf{f}(s, x) \cdot \mathbf{v}(s, x) dx ds.$$

Similarly, for any $n \in \mathbb{N}^*$, the energy balance satisfied by $(\mathbf{v}_n, k_n)$ reads

$$(4.74) \qquad \frac{1}{2} \int_{\Omega} |\mathbf{v}_n(t, x)|^2 dx + \iint_{Q_t} \nu_{\text{turb}}(k_n(s, x)) |\nabla \mathbf{v}_n(s, x)|^2 dx ds = \\ \frac{1}{2} \int_{\Omega} |\mathbf{v}_{0,n}(x)|^2 dx + \iint_{Q_t} \mathbf{f}(s, x) \cdot \mathbf{v}_n(s, x) dx ds.$$

We would like to pass to the limit in the energy balance (4.74) to prove that

$$(4.75) \qquad \iint \nu_{\text{turb}}(k_n(s, x)) |\nabla \mathbf{v}_n(s, x)|^2 dx ds \xrightarrow{n \to \infty} \iint \nu_{\text{turb}}(k(s, x)) |\nabla \mathbf{v}(s, x)|^2 dx ds,$$

which is unfortunately not possible because we cannot directly pass to the limit in the kinetic energy term $\frac{1}{2} \int_{\Omega} |\mathbf{v}_n(t, x)|^2 dx$.

To overcome this difficulty, we integrate the energy balance once more with respect to time. By Fubini's theorem,

$$\int_{0}^{T} \iint_{Q_t} \nu_{\text{turb}}(k_n(s, x)) |\nabla \mathbf{v}_n(s, x)|^2 dx ds dt = \int_{\Omega} \int_{0}^{T} \int_{0}^{t} \nu_{\text{turb}}(k_n(s, x)) |\nabla \mathbf{v}_n(s, x)|^2 ds dt dx.$$
Integrating by parts yields

$$ \int_0^T \left( \int_0^t \nu_{\text{turb}}(k_n(s,x)) |\nabla \mathbf{v}_n(s,x)|^2 ds \right) dt = \int_0^T (T-t) \nu_{\text{turb}}(k_n(t,x)) |\nabla \mathbf{v}_n(t,x)|^2 dt, $$

so that the energy balance (4.74) satisfied by $(\mathbf{v}_n, k_n)$ becomes

$$ (4.76) \qquad \frac{1}{2} \iint_{Q_T} |\mathbf{v}_n(t,x)|^2 dxdt + \iint_{Q_t} (T-t)\nu_{\text{turb}}(k_n(t,x)) |\nabla \mathbf{v}_n(t,x)|^2 dxdt = \frac{T}{2} \int_{\Omega} |\mathbf{v}_{0,n}(x)|^2 dx + \iint_{Q_t} (T-t)\mathbf{f}(t,x) \cdot \mathbf{v}_n(t,x) dxdt. $$

Similarly, we also find

$$ (4.77) \qquad \frac{1}{2} \iint_{Q_T} |\mathbf{v}(t,x)|^2 dxdt + \iint_{Q_t} (T-t)\nu_{\text{turb}}(k(t,x)) |\nabla \mathbf{v}(t,x)|^2 dxdt = \frac{T}{2} \int_{\Omega} |\mathbf{v}_0(x)|^2 dx + \iint_{Q_t} (T-t)\mathbf{f}(t,x) \cdot \mathbf{v}(t,x) dxdt. $$

From (4.45) we deduce that

$$ (4.78) \qquad \left\{ \begin{array}{l} \lim_{n \to \infty} \frac{1}{2} \iint_{Q_T} |\mathbf{v}_n(t,x)|^2 dxdt = \frac{1}{2} \iint_{Q_T} |\mathbf{v}(t,x)|^2 dxdt, \\ \\ \lim_{n \to \infty} \iint_{Q_t} (T-t)\mathbf{f}(t,x) \cdot \mathbf{v}_n(t,x) dxdt = \iint_{Q_t} (T-t)\mathbf{f}(t,x) \cdot \mathbf{v}(t,x) dxdt. \end{array} \right. $$

Therefore, letting $n \to \infty$ in (4.76) and comparing the result with (4.77), we obtain

$$ (4.79) \qquad \lim_{n \to \infty} \iint_{Q_t} (T-t)\nu_{\text{turb}}(k_n) |\nabla \mathbf{v}_n|^2 = \iint_{Q_t} (T-t)\nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2. $$

Recall that $Z_n$ and $Z$ are defined in (4.70), and that our goal is to prove the convergence of the sequence of norms of $Z_n$ to the norm of $Z$, namely (4.71), which is almost established in (4.79). In fact, it works for the sequence $(G_n)_{n \in \mathbb{N}}$ and $G$ given by

$$ (4.80) \qquad G_n = \sqrt{T-t} Z_n \quad \text{and} \quad G = \sqrt{T-t} Z. $$

As $Z_n \rightharpoonup Z$ weakly in $L^2(Q_T)^{N \times N}$, it is easily checked that

$$ (4.81) \qquad G_n \xrightarrow{n \to \infty} G \quad \text{weakly in} \quad L^2(Q_T)^{N \times N}. $$

Moreover, we deduce from (4.79) that,

$$ \|G_n\|_{Q_T;0,2} \xrightarrow{n \to \infty} \|G\|_{Q_T;0,2}, $$

hence the strong convergence of the sequence $(G_n)_{n \in \mathbb{N}}$ to $G$ in $L^2(Q_T)$. To conclude the $L^2$ convergence of the sequence $(Z_n)_{n \in \mathbb{N}}$ to $Z$, we use the Lebesgue inverse theorem. Hence, from the sequence $(G_n)_{n \in \mathbb{N}}$ we can extract a subsequence, still denoted by $(G_n)_{n \in \mathbb{N}}$, and there exists $0 \le J \in L^2(Q_T)$ such that

$$ (4.82) \qquad G_n \xrightarrow{n \to \infty} G \quad \text{a.e. in} \quad Q_T, \quad |G_n| \le J. $$
Let $T' < T$. As

$$ (4.83) \qquad Z_n = \frac{G_n}{\sqrt{(T-t)}}, $$

we have

$$ (4.84) \qquad Z_n \xrightarrow{n \to \infty} Z \quad \text{a.e. in } Q_{T'}, \quad |Z_n| \le \frac{J}{\sqrt{(T-T')}} \in L^2(Q_{T'}), $$

hence the strong convergence of $(Z_n)_{n \in \mathbb{N}}$ to $Z$ in $L^2(Q_{T'})^N$ by the Lebesgue theorem. As the limit is uniquely determined, the whole sequence converges. $\square$

### 4.3.4 Passing to the limit in the TKE-equation

We are now able to pass to the limit in the approximated TKE-equation (4.38). Passing to the limit in the transport term, the diffusion term, and the initial data follows exactly the same procedure as for the fluid equation, as already discussed in the previous subsection, see items Item 4.3.i, Item 4.3.ii, Item 4.3.iv and Item 4.3.v. It remains to deal with the dissipation and the production terms in the equation.

**Item 4.3.vi) Dissipation term $\varepsilon(k_n, k_n)$**. Let $q \in E_k$,

$$ q_n = \sum_{j=1}^{n} (\nabla q, \nabla \tau_j) \tau_j. $$

Recall that

$$ q_n \xrightarrow{n \to \infty} q \quad \text{in } E_k, $$

therefore this convergence holds also in $L^2(Q_T)$.

For both $N=2$ and $N=3$, we have $k_n \xrightarrow{n \to \infty} k$ strongly in $L^3(Q_T)$. Then

$$ (4.85) \qquad k_n |k_n|^{\frac{1}{2}} \xrightarrow{n \to \infty} k |k|^{\frac{1}{2}} \quad \text{strongly in } L^2(Q_T). $$

Since, by assumption, $\frac{1}{\ell} \in L^\infty(Q_T)$, we obtain

$$ (4.86) \qquad \langle \varepsilon(k_n, k_n), q_n \rangle = \iint_{Q_T} \frac{k_n |k_n|^{\frac{1}{2}} q_n}{\ell} \xrightarrow{n \to \infty} \iint_{Q_T} \frac{k |k|^{\frac{1}{2}} q}{\ell} = \langle \varepsilon(k, k), q \rangle. $$

**Item 4.3.vii) Production term $D_\alpha(k_n, v_n)$**. From Proposition 4.1, combined with Remark 4.7, we know that

$$ (4.87) \qquad \nu_{\text{turb}}(k_n) |\nabla v_n|^2 \xrightarrow{n \to \infty} \nu_{\text{turb}}(k) |\nabla v|^2 \quad \text{strongly in } L^1(Q_T). $$

It is easily deduced from (4.87), that for any given $\alpha > 0$,

$$ (4.88) \qquad T_\alpha(\nu_{\text{turb}}(k_n) |\nabla v_n|^2) \xrightarrow{n \to \infty} T_\alpha(\nu_{\text{turb}}(k) |\nabla v|^2) \quad \text{strongly in } L^2(Q_T). $$

Therefore,

$$ (4.89) \qquad \langle D_\alpha(k_n, v_n), q_n \rangle = \iint_{Q_T} T_\alpha(\nu_{\text{turb}}(k_n) |\nabla v_n|^2) q_n \xrightarrow{n \to \infty} \iint_{Q_T} T_\alpha(\nu_{\text{turb}}(k) |\nabla v|^2) q = \langle D_\alpha(k, v), q \rangle. $$
In conclusion, combining (4.54), (4.61), (4.86) and (4.89), we obtain, for all $q \in E_k$,

$$ (4.90) \qquad \begin{cases} \langle \partial_t k, q \rangle + \langle \tilde{B}_\alpha(\mathbf{v}, k), q \rangle + \langle \tilde{A}(k, k), q \rangle = \langle \mathbb{D}_\alpha(k, \mathbf{v}), q \rangle - \langle \varepsilon(k, k), q \rangle, \\ k|_{t=0} = k_{\alpha,0}. \end{cases} $$

Thus, we have shown that $(\mathbf{v}, k)$ also satisfies (4.21) and (4.23) in Definition 4.1, so that the proof is complete.

**Remark 4.8.** What we did in Remark 4.6 also works in this case. Then, let

$$ (4.91) \qquad \mathcal{E}_k = \{ q \in C^1(Q_T) \mid \forall (t, x) \in [0, T] \times \Gamma, q(t, x) = 0 \text{ and } \forall x \in \Omega, q(T, x) = 0 \}. $$

Therefore, for all $q \in \mathcal{E}_k$,

$$ (4.92) \qquad - \iint_{Q_T} \partial_t q(t, x) k(t, x) dxdt + \int_{\Omega} k_{\alpha,0}(x)q(0, x)dx + \langle \tilde{B}(\mathbf{v}, k), q \rangle + \\ \qquad \qquad \qquad \qquad \langle \tilde{A}(k, k), q \rangle = \langle \mathbb{D}(k, \mathbf{v}) - \varepsilon(k, k), q \rangle. $$

**Remark 4.9.** Assumption (3.3) on the boundedness of the eddy viscosities $v_{\text{turb}}$ and $\mu_{\text{turb}}$ plays a role at several stages of the proof. It is needed to obtain a proper weak formulation that yields an energy balance in the fluid equation, which is required in order to apply the energy method. It is also needed to pass to the limit in the diffusion terms, although, for this latter point, the passage to the limit could be handled under a suitable growth assumption. However, with the present method, we are not able to deal with physical eddy viscosities given by formulas (2.14) and (2.34). In the case of steady-state scalar toy models, there exist some results based on renormalization (see, for instance, [26], [35, Chapter 5]), but these techniques apparently do not extend to the full evolutionary model.

**Remark 4.10.** Assumption (3.5), namely $f \in C([0, T], L^2(\Omega))$, is essential in order to carry out the Galerkin method and to apply the Cauchy-Lipschitz theorem. However, once the result has been established for such source terms $f$, a result analogous to Theorem 4.1 can be obtained for $f \in L^2([0, T], V')$ by approximation, since $C([0, T], L^2(\Omega))$ is dense in $L^2([0, T], V')$.

## 4.4 Additional properties of the limit

From now, let $(\mathbf{v}_\alpha, k_\alpha)$ denote the weak solution to the T-NS-TKE system (4.3) corresponding to a given $\alpha > 0$, whose existence is now well established. We denote $(\mathbf{v}_{n,\alpha}, k_{n,\alpha})$ the corresponding smooth approximation constructed by the Galerkin method above. For all $(\mathbf{w}, q) \in E_v \times E_k$, one has

$$ (4.93) \qquad \langle \partial_t \mathbf{v}_\alpha, \mathbf{w} \rangle + \langle B_\alpha(\mathbf{v}_\alpha, \mathbf{v}_\alpha), \mathbf{w} \rangle + \langle A(k_\alpha, \mathbf{v}_\alpha), \mathbf{w} \rangle = \iint_{Q_T} \mathbf{f}(t, x) \cdot \mathbf{w}(t, x) dxdt, $$

$$ (4.94) \qquad \langle \partial_t k_\alpha, q \rangle + \langle \tilde{B}_\alpha(\mathbf{v}_\alpha, k_\alpha), q \rangle + \langle \tilde{A}(k_\alpha, k_\alpha), q \rangle = \langle \mathbb{D}_\alpha(k_\alpha, \mathbf{v}_\alpha), q \rangle - \langle \varepsilon(k_\alpha, k_\alpha), q \rangle, $$
with the initial data $\mathbf{v}_\alpha|_{t=0} = \mathbf{v}_0$ and $k_\alpha|_{t=0} = T_\alpha(k_0) = k_{\alpha,0}$.

**Remark 4.11.** Once the strong convergence of $(Z_n)_{n \in \mathbb{N}}$ to $Z$ in $L^2(Q_T)$ is established in Section 4.3.3, which yields the $L^1$ convergence of the source terms in the TKE equation, namely (4.88), we also obtain the strong convergence of the sequence $(k_{n,\alpha})_{n \in \mathbb{N}}$ to $k_\alpha$ in the space $E_k = L^2([0, T], H_0^1(\Omega))$ by the energy method. We omit the details.

This section is devoted to prove some additional properties satisfied by $(\mathbf{v}_\alpha, k_\alpha)$ that will be useful in the remainder, such as entropy calculus and the maximum principle.
### 4.4.1 Entropy calculus

In the next section, we let $\alpha$ tend to infinity in order to pass to the limit in Equations (4.93). The delicate point concerns the passage to the limit in the TKE equation. To this end, we need to apply the Boccardo-Gallouët lemma, Lemma A.2, which requires entropy balances. The purpose of this section is to provide the necessary tools to achieve this. Recall that $\mathcal{S}$ denotes the set of all Lipschitz functions $S: \mathbb{R} \to \mathbb{R}$ such that $S'$ has a finite number of discontinuities. The Stampacchia theorem is given by Theorem 3.3. We first need the following entropy properties.

In the following, $S \in \mathcal{S}$ such that $S(0) = 0$, and let

$$ (4.95) \qquad \forall u \in \mathbb{R}, \quad V(u) = \int_0^u S(v)dv. $$

**Lemma 4.6.** Let $\mathbb{P}_n$ be the orthogonal projector onto $K_n$ (see (4.29)). Therefore, $\forall n \in \mathbb{N}$, $S(k_{\alpha,n}) \in E_k$, $S(k_{\alpha}) \in E_k$, and

$$ (4.96) \qquad \mathbb{P}_n S(k_{\alpha,n}) \xrightarrow{n \to \infty} S(k_{\alpha}) \quad \text{weakly in } E_k. $$

*Proof.* As $k_{\alpha,n}$ and $k_{\alpha}$ are in $E_k$, for almost all $t \in [0, T]$, $k_{\alpha,n}(t, \cdot)$ and $k_{\alpha}(t, \cdot)$ are in $H_0^1(\Omega)$. Since $S \in \mathcal{S} \cap W^{1,\infty}(\mathbb{R})$ such that $S(0) = 0$, by the Stampacchia theorem, Theorem 3.3, we know that for such times $t$, $S(k_{n,\alpha})(t, \cdot)$ and $S(k_{\alpha}(t, \cdot))$ are in $H_0^1(\Omega)$ and we have a.e. in $Q_T$,

$$ \nabla S(k_{n,\alpha}) = S'(k_{n,\alpha})\nabla k_{n,\alpha}, \quad \nabla S(k_{\alpha}) = S'(k_{\alpha})\nabla k_{n,\alpha}. $$

As $S \in \mathcal{S}$ is a Lipchitz function, $S' \in L^\infty(\mathbb{R})$, therefore,

$$ (4.97) \qquad |\nabla S(k_{n,\alpha})| \le \|S'\|_{\infty} |\nabla k_{n,\alpha}|, \quad |\nabla S(k_{\alpha})| \le \|S'\|_{\infty} |\nabla k_{\alpha}|. $$

In particular, $S(k_{n,\alpha})$, $S(k_{\alpha}) \in E_k$ and the sequence $(S(k_{n,\alpha}))_{n \in \mathbb{N}}$ is bounded in $E_k$. Moreover, since $S(0) = 0$, $|S(x)| \le C|x|$ for some constant $C$, and by (4.45) we deduce

$$ (4.98) \qquad S(k_{n,\alpha}) \xrightarrow{n \to \infty} S(k_{\alpha}) \quad \begin{cases} \text{strongly in } L^2(Q_T), \\ \text{almost everywhere in } Q_T. \end{cases} $$

We also deduce from (4.97) combined with (4.45) that the sequence $(\nabla S(k_{n,\alpha}))_{n \in \mathbb{N}}$ is bounded in $L^2(Q_T)^N$. Then we can extract a subsequence, still denoted by $(\nabla S(k_{n,\alpha}))_{n \in \mathbb{N}}$, which converges weakly to some $G \in L^2(Q_T)$. Let $\psi \in \mathcal{D}(Q_T)^N$. Then, by (4.98),

$$ \iint_{Q_T} \nabla S(k_{n,\alpha}) \cdot \psi = - \iint_{Q_T} S(k_{n,\alpha}) \mathrm{div} \psi \xrightarrow{n \to \infty} \\ - \iint_{Q_T} S(k_{\alpha}) \mathrm{div} \psi = \iint_{Q_T} \nabla S(k_{\alpha}) \cdot \psi. $$

By consequence, $G = \nabla S(k_{\alpha})$ in $\mathcal{D}'(Q_T)^N$ and therefore in $L^2(Q_T)^N$, so that we have proved that

$$ (4.99) \qquad S(k_{\alpha,n}) \xrightarrow{n \to \infty} S(k_{\alpha}) \quad \text{weakly in } E_k. $$

Let now $q \in L^2([0, T], K_j)$ for a given $j \in \mathbb{N}^*$, where the finite dimensional space $K_j$ is defined by (4.29). For a given $n \ge j$, we have thanks to (4.99),

$$ \int_0^T (\mathbb{P}_n S(k_{\alpha,n}), q) dt = \int_0^T (S(k_{\alpha,n}), q) dt = \iint_{Q_T} \nabla q \cdot \nabla S(k_{\alpha,n}) \xrightarrow{n \to \infty} \iint_{Q_T} \nabla q \cdot \nabla S(k_{\alpha}), $$

which is true for any $q \in L^2([0, T], K_j)$ and any $j \in \mathbb{N}$, hence (4.96). $\square$
**Remark 4.12.** We have not been able to prove the strong convergence of the sequence $(S(k_{\alpha,n}))_{n \in \mathbb{N}}$ to $S(k_{\alpha})$ in $E_k$. This property holds true in the steady-state case, when the time variable $t$ is not involved (see for instance [41]), but things do not work anymore as soon as $t$ is involved, and only weak convergence can be proved, which is enough for our purpose.

**Remark 4.13.** With the same argument, we also have

$$ (4.100) \qquad \mathbb{P}_n S(k_{\alpha,n}) \xrightarrow{n \to \infty} S(k_{\alpha}) \quad \text{weakly in } L^2(Q_T). $$

We conjecture that in this case the convergence is strong; however, weak convergence is sufficient for the sequel.

**Lemma 4.7.** We have $\forall t \in ]0, T]$, (4.101)

$$ \begin{aligned} \int_{\Omega} V(k_{\alpha}(t, x))dx + \iint_{Q_t} S'(k_{\alpha}(s, x))\mu_{\text{turb}}(k_{\alpha}(s, x))|\nabla k_{\alpha}(s, x)|^2 dxds = \\ \int_{\Omega} V(k_{\alpha,0}(x))dx + \iint_{Q_t} T_{\alpha}(\nu_{\text{turb}}(k_{\alpha}(s, x))|\nabla \mathbf{v}_{\alpha}(s, x)|^2)S(k_{\alpha}(s, x))dxds - \\ \iint_{Q_t} \frac{k_{\alpha}(s, x)|k_{\alpha}(s, x)|^{\frac{1}{2}}S(k_{\alpha}(s, x))}{\ell} dxds. \end{aligned} $$

*Proof.* We argue by approximations. Let $t \in [0, T]$. We take $\mathbb{P}_n S(k_{\alpha,n}) \in L^2([0, T], K_n)$ as test in the TKE equation (4.38), by noting that as $\partial_t k_{\alpha,n} \in C^0([0, T], K_n)$,

$$ \int_0^t (\partial_t k_{\alpha,n}, \mathbb{P}_n S(k_{\alpha,n})) = \iint_{Q_t} \partial_t k_{\alpha,n}(s,x) S(k_{\alpha,n}(s,x)) dxds, $$

and since $k_{\alpha,n} \in C^1([0, T], K_n)$ and $K_n \subset C^{\infty}(\Omega)$, and $S \in C^0(\mathbb{R})$, we have

$$ \begin{aligned} \iint_{Q_t} \partial_t k_{\alpha,n}(s,x) S(k_{\alpha,n}(s,x)) dxds &= \int_0^t \left( \frac{d}{dt} \int_{\Omega} V(k_{\alpha,n}(s,x)) dx \right) ds \\ &= \int_{\Omega} V(k_{\alpha,n}(t,x)) dx - \int_{\Omega} V(k_{\alpha,0,n}(x)) dx, \end{aligned} $$

where $k_{\alpha,0,n}$ is defined by (4.31). Therefore we get (4.102)

$$ \begin{aligned} & \int_{\Omega} V(k_{\alpha,n}(t,x))dx - \iint_{Q_t} \mathbf{v}_{n,\alpha}(s,x)T_{\alpha}(k_{\alpha,n}(s,x)) \cdot \nabla \mathbb{P}_n S(k_{\alpha,n}(s,x))dxds + \\ & \iint_{Q_t} \mu_{\text{turb}}(k_{\alpha,n}(s,x))\nabla k_{\alpha,n}(s,x) \cdot \nabla \mathbb{P}_n S(k_{\alpha,n}(s,x))dxds = \\ & \int_{\Omega} V(k_{\alpha,0,n}(x))dx + \iint_{Q_t} T_{\alpha}(\nu_{\text{turb}}(k_{\alpha,n}(s,x))|\nabla \mathbf{v}_{\alpha,n}(s,x)|^2)\mathbb{P}_n S(k_{\alpha,n}(s,x))dxds - \\ & \iint_{Q_t} \frac{k_{\alpha,n}(s,x)|k_{\alpha,n}(s,x)|^{\frac{1}{2}}\mathbb{P}_n S(k_{\alpha}(s,x))}{\ell} dxds. \end{aligned} $$

In order to obtain Formula (4.101), we have to pass to the limit in Formula (5.14) when $n \to \infty$. Each term in (5.14) is handled by combining the convergence results above (4.96) and (4.100), Lemma 4.2, Remark 4.11, Lemma A.5, (4.32), (4.45), (4.46), (4.85), (4.88),
and standard arguments. The only difficulty arises from the term $\int_{\Omega} V(k_{\alpha,n}(t,x)) dx$, which we detail below.

We observe that

$$ \forall x \in \mathbb{R}, \quad |V(x)| \leq C|x|^2. $$

Therefore, by (4.45),

$$ (4.103) \qquad V(k_{\alpha,n}) \xrightarrow{n \to \infty} V(k_{\alpha}) \quad \text{in } L^1(Q_T). $$

Let

$$ u_n(t) = \int_{\Omega} V(k_{\alpha,n}(t,x)) dx. $$

We deduce from (4.103) that the sequence $(u_n)_{n \in \mathbb{N}}$ converges to $u(t) = \int_{\Omega} V(k_{\alpha}(t,x)) dx$ in $L^1([0, T])$. By the Lebesgue inverse theorem, we can extract from the sequence $(u_n)_{n \in \mathbb{N}}$ a subsequence, still denoted by $(u_n)_{n \in \mathbb{N}}$, that converges to $u$ almost everywhere in $[0, T]$. In other words, for this subsequence, and for a.e. $t \in [0, T]$,

$$ \int_{\Omega} V(k_{\alpha,n}(t,x)) dx \xrightarrow{n \to \infty} \int_{\Omega} V(k_{\alpha}(t,x)) dx, $$

which makes it possible to pass to the limit in Formula (5.14) as $n \to \infty$. $\square$

### 4.4.2 Maximum principle

The TKE $k$ is, by definition, an energy. Therefore, it must be a non-negative function. This is why we establish the following maximum principle.

**Lemma 4.8.** Assume that $k_0 \geq 0$ a.e. in $\Omega$, then $k_{\alpha} \geq 0$ a.e. in $Q_T$.

*Proof.* Let

$$ k_{\alpha}^{+} = \frac{|k_{\alpha}| + k_{\alpha}}{2}, \quad k_{\alpha}^{-} = \frac{|k_{\alpha}| - k_{\alpha}}{2}, $$

so that

$$ k_{\alpha} = k_{\alpha}^{+} - k_{\alpha}^{-}, \quad |k_{\alpha}| = k_{\alpha}^{+} + k_{\alpha}^{-}. $$

As $k_0 \geq 0$, $k_{0,\alpha} = T_{\alpha}(k_0) \geq 0$ and $T_{\alpha}(k_0)^- = 0$ a.e. in $\Omega$. The aim of the proof is to prove that $k_{\alpha}^- = 0$ a.e. in $Q_T$, by using the result of Lemma 4.7.

According to Stampacchia Theorem (see also in appendice B), we see that $k_{\alpha}^{+}, k_{\alpha}^{-} \in E_k$, and since $\nabla k_{\alpha} = 0$ a.e. on $\{k_{\alpha} = 0\}$ (see also Remark B.2), which combined with formula (B.1) and $k_{\alpha}^{+} = 0$ on $\{k_{\alpha} < 0\}$, $k_{\alpha}^{-} = 0$ on $\{k_{\alpha} > 0\}$, yields the following formula

$$ (4.104) \qquad \iint_{Q_T} \mu_{\text{turb}}(k_{\alpha}) \nabla k_{\alpha}^{+} \cdot \nabla k_{\alpha}^{-} = \iint_{\{k_{\alpha}=0\}} \mu_{\text{turb}}(k_{\alpha}) \nabla k_{\alpha}^{+} \cdot \nabla k_{\alpha}^{-} + \\ \iint_{\{k_{\alpha}>0\}} \mu_{\text{turb}}(k_{\alpha}) \nabla k_{\alpha}^{+} \cdot \nabla k_{\alpha}^{-} + \iint_{\{k_{\alpha}<0\}} \mu_{\text{turb}}(k_{\alpha}) \nabla k_{\alpha}^{+} \cdot \nabla k_{\alpha}^{-} = 0. $$

Let

$$ S(x) = \frac{x - |x|}{2} \in \mathbb{S}, \quad V(x) = \begin{cases} 0 & \text{if } 0 \leq x, \\ \frac{1}{2}x^2 & \text{if } x \leq 0. \end{cases} $$

In particular

$$ S(k_{\alpha}) = -k_{\alpha}^{-}, \quad V(k_{\alpha}) = \frac{1}{2}|k_{\alpha}|^2 \chi_{\{k_{\alpha} < 0\}} = \frac{1}{2}|k_{\alpha}^{-}|^2. $$
Therefore, as $k_{\alpha}^{-}(0, x) = 0$ a.e. in $\Omega$ and $V(0) = 0$, applying the formula (4.101) combined with the formula (4.104) yields

$$ (4.105) \qquad \frac{1}{2} \int_{\Omega} |k_{\alpha}^{-}(t, x)|^2 dx + \iint_{Q_t} \mu_{\text{turb}}(k_{\alpha}(s, x)) |\nabla k_{\alpha}^{-}(s, x)|^2 dx ds = \\ \iint_{Q_t} T_{\alpha}(\nu_{\text{turb}}(k_{\alpha}(s, x))) |\nabla \mathbf{v}_{\alpha}(s, x)|^2 S(k_{\alpha}(s, x)) dx ds - \\ \iint_{Q_t} \frac{k_{\alpha}(s, x) |k_{\alpha}(s, x)|^{\frac{1}{2}} S(k_{\alpha}(s, x))}{\ell} dx ds. $$

As $S \le 0$ and $0 \le xS(x)$, we deduce from (4.105) that at any time $t \in [0, T]$,

$$ \int_{\Omega} |k_{\alpha}^{-}(t, x)|^2 dx \le 0, $$

hence $k_{\alpha}^{-} = 0$ a.e. in $Q_T$ and the positivity of $k_{\alpha}$ follows. $\square$

# 5 Existence result for the NS-TKE system

We aim to prove that, as $\alpha \to \infty$, the family $\{(\mathbf{v}_{\alpha}, k_{\alpha})\}_{\alpha>0}$ (which we will also refer to as a sequence when no confusion arises) converges to a solution (case $N = 2$), or a subsolution (case $N = 3$), $(\mathbf{v}, k)$ of the NS-TKE system (3.1) in a sense to be defined precisely. We still assume that Assumptions (3.3), (3.4) and (3.5) hold. Moreover, we also assume that

$$ (5.1) \qquad \mathbf{v}_0 \in H, \quad k_0 \in L^1(\Omega), \quad k_0 \ge 0 \text{ a.e. in } \Omega. $$

We begin with some elementary remarks that motivate what will be done in this section. Since

$$ \forall x \in \mathbb{R}, \quad \lim_{\alpha \to \infty} T_{\alpha}(x) = x, $$

we have, at least formally,

$$ (5.2) \qquad B(\mathbf{v}, \mathbf{v}) = \lim_{\alpha \to \infty} B_{\alpha}(\mathbf{v}, \mathbf{v}), \quad \tilde{B}(k, \mathbf{v}) = \lim_{\alpha \to \infty} \tilde{B}_{\alpha}(k, \mathbf{v}), $$

where the truncated transport terms $B_{\alpha}$ and $\tilde{B}_{\alpha}$ are given by (4.5), and the transport terms $B$ and $\tilde{B}$ are given by

$$ (5.3) \qquad \begin{aligned} \langle B(\mathbf{v}, \mathbf{w}), \mathbf{u} \rangle &= \iint_{Q_T} (\mathbf{v} \cdot \nabla) \mathbf{w} \cdot \mathbf{u} = - \iint_{Q_T} \mathbf{w} \otimes \mathbf{v} : \nabla \mathbf{u}, \\ \langle \tilde{B}(\mathbf{v}, k), q \rangle &= \iint_{Q_T} (\mathbf{v} \cdot \nabla k) q = - \iint_{Q_T} k \mathbf{v} \cdot \nabla q. \end{aligned} $$

Moreover,

$$ (5.4) \qquad \mathbb{D}(k, \mathbf{v}) = \lim_{\alpha \to \infty} \mathbb{D}_{\alpha}(k, \mathbf{v}), $$

where

$$ (5.5) \qquad \langle \mathbb{D}(k, \mathbf{v}), r \rangle = \iint_{Q_T} \nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2 r. $$
In the remainder of this section, we start with the main estimates for {$v_α$}$_{α>0}$ derived from the energy balance (4.73). Then we obtain estimates for {$k_α$}$_{α>0}$ by means of the Boccardo-Gallouët lemma (Lemma A.2). From this point, we distinguish between the cases $N = 2$ and $N = 3$. The difference lies in the fact that for $N = 2$, the energy method developed in Section 4.3.3 can still be applied, while for $N = 3$ this is no longer possible: the limit velocity cannot be used as a test vector field because the transport term lacks sufficient regularity. Consequently, it satisfies only an energy inequality, not a full energy balance.

## 5.1 Main estimates

### 5.1.1 Estimates for {$v_α$}$_{α>0}$

We know from Section 4.3.3 that, for any $α > 0$ and any $t \in [0, T]$, $v_α$ satisfies the energy balance

$$ (5.6) \qquad \frac{1}{2} \int_{\Omega} |v_{\alpha}(t, x)|^2 dx + \iint_{Q_t} \nu_{\text{turb}}(k_{\alpha}(s, x)) |\nabla v_{\alpha}(s, x)|^2 dx ds = \\ \frac{1}{2} \int_{\Omega} |v_0(x)|^2 dx + \iint_{Q_t} \mathbf{f}(s, x) \cdot v_{\alpha}(s, x) dx ds, $$

which, combined with the interpolation inequalities recalled in Remark 3.6, estimates (3.21) and (3.22), allows us to deduce that {$v_α$}$_{α>0}$ satisfies:

$$ (5.7) \qquad \{v_\alpha\}_{\alpha>0} \text{ is bounded in} \quad \left\{ \begin{array}{l} E_v = L^2([0, T], V), \\ L^\infty([0, T], L^2(\Omega)^N), \\ L^4(Q_T) \quad \text{if } N = 2, \\ L^{10/3}(Q_T) \quad \text{if } N = 3. \end{array} \right. $$

### 5.1.2 Estimates for {$k_α$}$_{α>0}$ by Boccardo-Gallouët Lemma

The process we carried out in Section 3.4 holds for each $k_α$, combined with the result of Lemma 4.7. In particular, for any $α > 0$, $k_α$ satisfies the estimate (3.34), which, combined with (5.7), yields for any $t \in [0, T]$,

$$ (5.8) \qquad \int_{\Omega} |k_{\alpha}(t, x)| dx \le \int_{\Omega} |k_{\alpha,0}(x)| dx + \beta \|\nabla v_{\alpha}\|^2_{Q_{T;0,1}} \le \\ \|k_0\|_{\Omega;0,1} + C(\|v_0\|_{\Omega;0,2}, \|f\|_{Q_T;0,2}, \nu, \beta). $$

We need a bound for {$k_α$}$_{α>0}$ in some $L^q([0, T], W^{1,q}(\Omega))$. To do so, we will use the Boccardo-Gallouët lemma (Lemma A.2 in Appendix A). We first notice that, since $k_α \in L^2(Q_T)$ for any fixed $α > 0$, the assumptions (A.7), (A.9) and (A.10) of Section A.3 about the Boccardo-Gallouët Lemma are satisfied by each $k_α$, by setting

$$ (5.9) \qquad \mathbb{B}_n(k_\alpha) = \{(t, x) \in Q_T; n \le |k_\alpha(t, x)| < n+1\}. $$

Let $M_n(k_α)$ denotes the integral

$$ (5.10) \qquad M_n(k_\alpha) = \iint_{\mathbb{B}_n(k_\alpha)} |\nabla k_\alpha(t, x)|^2 dx dt. $$

According to assumption (A.13) in the Boccardo-Gallouët lemma, we must obtain a uniform bound for $M_n(k_α)$ in $n$ and $α$. More precisely,
**Lemma 5.1.** There exists a constant $C = C(\|\mathbf{v}_0\|_{\Omega,0,2}, \|k_0\|_{\Omega,0,1}, \|\mathbf{f}\|_{Q_{T;0,2}}, \nu, \beta)$ which does not depend on $\alpha$ nor on $n$ such that

$$ (5.11) \qquad \forall n \in \mathbb{N}^*, \forall \alpha > 0, \quad 0 \le M_n(k_\alpha) \le C. $$

Therefore, for all $1 \le q < \frac{N+2}{N+1}$, there exists a constant $C_q = C(q, \Omega, \|\mathbf{v}_0\|_{\Omega,0,2}, \|\mathbf{f}\|_{Q_{T;0,2}}, \nu)$ such that

$$ (5.12) \qquad \|k_\alpha\|_{q;\Omega;1,q} \le C_q. $$

*Proof.* Let $n \in \mathbb{N}^*$, and define $S_n \in \mathbb{S}$ as the function whose graph is given below.

<figure><img src="image_7.png" /><figcaption>Figure 3: Graph of the odd function $S_n$.</figcaption></figure>

It is easily checked that the following relation holds:

$$ S_n = T_{n+1} - T_n. $$

By the Stampacchia theorem (Theorem 3.3), we have

$$ (5.13) \qquad \forall x \in \mathbb{R}, \quad S'_n(x) = \chi_{[n,n+1]}(x) + \chi_{[-n-1,-n]}(x), \quad |S_n(x)| \le 1, \\ \forall f \in H_0^1(\Omega), \quad \left\{ \begin{array}{l} \nabla S_n(f) = (\chi_{[n,n+1]}(f) + \chi_{[-n-1,-n]}(f))\nabla f \\ = \chi_{\mathbb{B}_n(f)} \nabla f. \end{array} \right. $$

Let

$$ V_n(x) = \int_0^x S_n(x)dx. $$

Combining Formula (4.101) with (5.13) yields

$$ (5.14) \qquad \int_{\Omega} V_n(k_{\alpha}(t, x))dx + \iint_{\mathbb{B}_n(k_{\alpha})} |\nabla k_{\alpha}(s, x)|^2 dxds = \\ \qquad \int_{\Omega} V_n(k_{\alpha,0}(x))dx + \iint_{Q_t} T_{\alpha}(\nu_{\text{turb}}(k_{\alpha}(s, x))|\nabla \mathbf{v}_{\alpha}(s, x)|^2)S_n(k_{\alpha}(s, x))dxds - \\ \qquad \qquad \iint_{Q_t} \frac{k_{\alpha}(s, x)|k_{\alpha}(s, x)|^{\frac{1}{2}}S_n(k_{\alpha}(s, x))}{\ell}dxds. $$

Observe that

$$ (5.15) \qquad \forall x \in \mathbb{R}, \quad 0 \le V_n(x) \le |x|. $$

Therefore,

$$ (5.16) \qquad 0 \le \int_{\Omega} V_n(k_{\alpha}(T, x)) dx, \quad 0 \le \int_{\Omega} V_n(k_{\alpha,0}(x)) dx \le \|k_{\alpha,0}\|_{\Omega;0,1}. $$
Concerning the production term, we deduce from the energy balance (5.6), together with the Cauchy–Schwarz and Poincaré inequalities, that

$$ (5.17) \qquad \left| \iint_{Q_t} T_\alpha(\nu_{\text{turb}}(k_\alpha(s, x)) |\nabla \mathbf{v}_\alpha(s, x)|^2) S_n(k_\alpha(s, x)) dx ds \right| \\ \le \|\nu_t(k_\alpha) |\nabla \mathbf{v}_\alpha|^2\|_{Q_T; 0,1} \le C (\|\mathbf{v}_0\|_{\Omega; 0,2}, \|\mathbf{f}\|_{Q_T; 0,2}, \nu). $$

Moreover, since $S_n$ is an odd function,

$$ (5.18) \qquad 0 \le \iint_{Q_t} \frac{k_\alpha(s, x) |k_\alpha(s, x)|^{\frac{1}{2}} S_n(k_\alpha(s, x))}{\ell} dx ds. $$

Finally, we also have

$$ (5.19) \qquad \forall \alpha > 0, \quad \|k_{\alpha,0}\|_{\Omega;0,1} \le \|k_0\|_{\Omega;0,1}. $$

Combining (5.14), (5.16), (5.17), (5.18) and (5.19), we obtain for all $n \in \mathbb{N}^*$,

$$ 0 \le M_n(k_\alpha) = \iint_{\mathbb{B}_n(k_\alpha)} |\nabla k_\alpha|^2 \le \frac{1}{\nu} \left[ C(\|\mathbf{v}_0\|_{\Omega;0,2}, \|\mathbf{f}\|_{Q_T;0,2}, \nu) + \|k_0\|_{\Omega;0,1} \right]. $$

Hence (5.11).

By the Boccardo–Gallouët lemma, we finally obtain

$$ \forall 1 \le p < \frac{N+2}{N+1}, \quad \{k_\alpha\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{1,p}(\Omega)). $$

□

From there, we distinguish the cases $N = 2$ and $N = 3$.

## 5.2 Estimates and test functions spaces when $N = 2$

We begin this section by deriving additional estimates from the previous ones, with particular emphasis on the transport, production, and dissipation terms appearing on the right-hand side of the TKE equation. These estimates allow us to give later a precise meaning to the notion of a weak solution to the NS–TKE system (3.1).

### 5.2.1 Variational formulation for the fluid equation

We specify first additional bounds for the transport and diffusion terms in the fluid equation, then the time derivative. We give the expected variational formulation in the limit as $\alpha \to \infty$.

**Item 5.2.i) Transport term.** Let $(B_\alpha(\mathbf{v}_\alpha, \mathbf{v}_\alpha))_{\alpha>0}$ be the sequence of truncated transport terms (see (4.5) and (4.93)), which we estimate below. Since $N = 2$, we already know that the sequence $\{\mathbf{v}_\alpha\}_{\alpha>0}$ is bounded in $L^4(Q_T)$ (see (5.7)); hence

$$ (5.20) \qquad \mathbf{v}_\alpha \otimes \mathbf{v}_\alpha \text{ is bounded in } L^2(Q_T). $$

In particular,

$$ (5.21) \qquad B_\alpha(\mathbf{v}_\alpha, \mathbf{v}_\alpha) \text{ is bounded in } L^2([0, T], H^{-1}(\Omega)^N), $$
and therefore in $L^2([0, T], V')$, where we recall that

$$ \langle B_{\alpha}(\mathbf{v}_{\alpha}, \mathbf{v}_{\alpha}), \mathbf{w} \rangle = - \iint_{Q_T} T_{\alpha}(\mathbf{v}_{\alpha}) \otimes \mathbf{v}_{\alpha} : \nabla \mathbf{w}. $$

Moreover, since $|T_{\alpha}(v_i)| \le |v_i|$ for $i = 1, 2, 3$, we have

$$ (5.22) \qquad \|B_{\alpha}(\mathbf{v}_{\alpha}, \mathbf{v}_{\alpha})\|_{2;\Omega;-1,2} \le \|\mathbf{v}_{\alpha}\|_{Q_{T;0,4}}^2 \le C(\|\mathbf{v}_0\|_{\Omega;0,2}, \|\mathbf{f}\|_{Q_{T;0,2}}, \nu). $$

**Item 5.2.ii) Diffusion term.** Let

$$ (5.23) \qquad Z_{\alpha} = \sqrt{\nu_{\text{turb}}(k_{\alpha})} \nabla \mathbf{v}_{\alpha}. $$

We deduce from the energy balance (5.6) that

$$ (5.24) \qquad \{Z_{\alpha}\}_{\alpha>0} \text{ is bounded in } L^2(Q_T)^{N \times N}. $$

As

$$ \langle A(k_{\alpha}, \mathbf{v}_{\alpha}), \mathbf{w} \rangle = \iint_{Q_T} Z_{\alpha} : \nabla \mathbf{w}, $$

we deduce that

$$ (5.25) \qquad \{A(k_{\alpha}, \mathbf{v}_{\alpha})\}_{\alpha>0} \text{ is bounded in } L^2([0, T], V'). $$

**Item 5.2.iii) Time derivative.** We write the fluid equation as

$$ (5.26) \qquad \partial_t \mathbf{v}_\alpha = -B_\alpha(\mathbf{v}_\alpha, \mathbf{v}_\alpha) - A(k_\alpha, \mathbf{v}_\alpha) + \mathbf{f}. $$

We deduce from (5.21) and (5.25), that

$$ (5.27) \qquad \{\partial_t \mathbf{v}_\alpha\}_{\alpha>0} \text{ is bounded in } L^2([0, T], V'). $$

**Item 5.2.iv) Variational formulation.** Therefore, the weak formulation of the fluid equation in the limit holds in the same way as its $\alpha$-approximation. Denoting by $(\mathbf{v}, k)$ the limit of the sequence $\{\mathbf{v}_\alpha, k_\alpha\}_{\alpha>0}$ in a sense that remains to be specified, we have

$$ \mathbf{v} \in E_v = L^2([0, T], V), \quad \partial_t \mathbf{v} \in L^2([0, T], V'), \quad \mathbf{v} \in C([0, T], H), \quad \mathbf{v}|_{t=0} = \mathbf{v}_0, $$

and

$$ (5.28) \qquad \forall \mathbf{w} \in E_v, \quad \langle \partial_t \mathbf{w}, \mathbf{w} \rangle + \langle B(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle + \langle A(k, \mathbf{v}), \mathbf{w} \rangle = (\mathbf{f}, \mathbf{w}). $$

We will discuss the regularity of $k$ in the next section. Note that in this formulation, one can take $\mathbf{w} = \mathbf{v}$ as a test function, so that the energy balance is satisfied. This is essential for passing to the limit in the production term in the TKE equation by the energy method (see Section 4.3.3).

## 5.2.2 More estimates for the TKE equation

Things are different for the transport term in the TKE equation (4.94). We must estimate the product $k_\alpha \mathbf{v}_\alpha$ in some $L^q(Q_T)$ space, in order to apply Lemma A.3. We also estimate the diffusion term, and the production and the dissipation terms in the r.h.s of the TKE equation in suitable dual spaces.
**Item 5.2.v)** $L^p(W^{1,p})$ estimate. When $N=2$, we deduce from Lemma 5.1 that

(5.29) $\forall 1 \le p < \frac{4}{3}$, $\{k_\alpha\}_{\alpha>0}$ is bounded in $L^p([0, T], W^{1,p}(\Omega)) \cap L^\infty([0, T], L^1(\Omega))$.

Therefore, it is natural to seek $k = \lim_{\alpha \to \infty} k_\alpha$ in the set

$$ (5.30) \qquad \mathcal{K}_{k,2} = \bigcap_{1 \le p < \frac{4}{3}} L^p([0, T], W^{1,p}(\Omega)). $$

**Item 5.2.vi)** $L^\zeta(Q_T)$ estimate. When $N=2$, by the Sobolev theorem we know that

$$ \forall 1 \le p \le 4, \quad W^{1,\frac{4}{3}}(\Omega) \hookrightarrow L^p(\Omega). $$

Therefore, since we know that $\{k_\alpha\}_{\alpha>0}$ is bounded in $L^\infty([0, T], L^1(\Omega))$, we deduce from Lemma A.4 and formula (A.21), with $\zeta(4/3, 4, 1) = 2$ that

$$ (5.31) \qquad \forall 1 \le p < 2, \quad \{k_\alpha\}_{\alpha>0} \text{ is bounded in } L^p(Q_T). $$

**Item 5.2.vii)** Transport term. since $\{\mathbf{v}_\alpha\}_{\alpha>0}$ is bounded in $L^4(Q_T)$, by the Hölder inequality combined with (5.31), we find that

$$ (5.32) \qquad \forall 1 \le p < \frac{4}{3}, \quad \{k_\alpha \mathbf{v}_\alpha\}_{\alpha>0} \text{ is bounded in } L^p(Q_T)^2, $$

since $|T_\alpha(k_\alpha)| \le |k_\alpha|$, we deduce from Lemma A.3 that

$$ (5.33) \qquad \forall 1 \le p < \frac{4}{3}, \quad \{\tilde{B}_\alpha(k_\alpha, \mathbf{v}_\alpha)\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{-1,p}(\Omega)), $$

where we recall that

$$ \langle B_\alpha(k_\alpha, \mathbf{v}_\alpha), q \rangle = - \iint_{Q_T} T_\alpha k_\alpha \mathbf{v}_\alpha \cdot \nabla q. $$

**Item 5.2.viii)** Diffusion term. Let

$$ \tilde{Z}_\alpha = \sqrt{\mu_{\text{turb}}(k_\alpha)} \nabla k_\alpha, $$

so that

$$ \langle \tilde{A}(k_\alpha, \mathbf{v}_\alpha), q \rangle = \iint_{Q_T} \tilde{Z}_\alpha \cdot \nabla q. $$

As $\mu_{\text{turb}}$ is bounded, we have by (5.29) that $\{Z_\alpha\}_{\alpha>0}$ is bounded in $L^p(Q_T)$ for any $1 \le p < 4/3$. Therefore,

$$ (5.34) \qquad \forall 1 \le p < \frac{4}{3}, \quad \{\tilde{A}(k_\alpha, \mathbf{v}_\alpha)\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{-1,p}(\Omega)). $$

**Item 5.2.ix)** Production term. From the energy balance (5.6), together with the Cauchy–Schwarz and Poincaré inequalities, we deduce that

$$ (5.35) \qquad \{\mathbb{D}_\alpha(k_\alpha, \mathbf{v}_\alpha)\}_{\alpha>0} \text{ is bounded in } L^1(Q_T). $$

Let $\alpha$ be given, and let $\mathbb{D}_\alpha(k_\alpha, \mathbf{v}_\alpha) = \mathbb{D}_\alpha$ for the simplicity. As $L^1(Q_T) \approx L^1([0, T], L^1(\Omega))$, for almost all $t \in [0, T]$, $D_\alpha(t, \cdot) \in L^1(\Omega)$. Let $r > 4$, $q = q(x) \in W_0^{1,r}(\Omega)$. In this case $W_0^{1,r}(\Omega) \hookrightarrow L^\infty(\Omega)$. Therefore, for a.e. $t \in [0, T]$,

$$ \left| \int_{\Omega} \mathbb{D}_{\alpha}(t, x) q(x) dx \right| \le \|\mathbb{D}_{\alpha}(t, \cdot)\|_{\Omega;0,1} \|q\|_{\Omega;0,\infty} \le C \|\mathbb{D}_{\alpha}\|_{\Omega;0,1} \|q\|_{\Omega;1,r}, $$
which yields

$$ \mathbb{D}_{\alpha}(t, \cdot) \in W^{-1,r'}(\Omega), \quad \text{and} \quad \|\mathbb{D}_{\alpha}(t, \cdot)\|_{\Omega;-1,r'} \le C \|\mathbb{D}_{\alpha}(t, \cdot)\|_{\Omega;0,1}, $$

which we integrate on time over $[0, T]$ to get

$$ \int_{0}^{T} \|\mathbb{D}_{\alpha}(t, \cdot)\|_{\Omega;-1,r'} dt \le C \|D_{\alpha}\|_{Q_T;0,1}, $$

which is bounded uniformly in $\alpha$ by (5.35). In other words,

$$ (5.36) \quad \forall 1 \le p < \frac{4}{3}, \quad \{\mathbb{D}_{\alpha}(k_{\alpha}, \mathbf{v}_{\alpha})\}_{\alpha>0} \text{ is bounded in } L^1([0, T], W^{-1,p}(\Omega)). $$

**Item 5.2.x)** *Dissipation term.* Recall that as we assume $k_0 \ge 0$, we know that $k_\alpha \ge 0$ a.e. in $Q_T$ (see Lemma 4.8 and the assumption (5.1)). By (5.31), we have that $k_\alpha \sqrt{k_\alpha}$ is bounded in $L^p(Q_T)$ for any $1 \le p < \frac{4}{3}$, hence

$$ (5.37) \quad \forall 1 \le p < \frac{4}{3}, \quad \{\varepsilon(k_{\alpha}, k_{\alpha})\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{-1,p}(\Omega)). $$

**Item 5.2.xi)** *Time derivative.* By writing

$$ (5.38) \quad \partial_t k_{\alpha} = -\tilde{B}_{\alpha}(\mathbf{v}_{\alpha}, k_{\alpha}) - \tilde{A}(k_{\alpha}, k_{\alpha}) + \mathbb{D}_{\alpha}(k_{\alpha}, \mathbf{v}_{\alpha}) - \varepsilon(k_{\alpha}, k_{\alpha}). $$

and combining (5.33), (5.34), (5.36) and (5.37), we obtain

$$ (5.39) \quad \forall 1 \le p < \frac{4}{3}, \quad \{\partial_t k_{\alpha}\}_{\alpha>0} \text{ is bounded in } L^1([0, T], W^{-1,p}(\Omega)). $$

**Item 5.2.xii)** *Conclusion : suitable test function space.* Considering the transport diffusion and dissipation terms, we are inclined to choose test functions in $L^q([0, T], W_0^{1,q}(\Omega))$ for $q > 4$. However, we also have to deal with the production term. Although

$$ W_0^{1,q}(\Omega) \hookrightarrow L^{\infty}(\Omega), $$

the time regularity of $k$ does not match, and we should take functions in the space $C^1([0, T], W_0^{1,q}(\Omega))$, $q > 4$. Therefore, for simplicity, and in order to directly take the initial data into consideration in the variational formulation of the TKE equation, we will take test functions in the space

$$ (5.40) \quad \mathcal{E}_k = \{q \in C^1(Q_T) \mid \forall (t, x) \in [0, T] \times \Gamma, q(t, x) = 0 \text{ and } \forall x \in \Omega, q(T, x) = 0\}. $$

## 5.3 Weak solutions for $N = 2$

We define in this section the notion of weak solutions to the NS-TKE system (3.1) when $N = 2$. We then prove that the sequence $\{\mathbf{v}_\alpha, k_\alpha\}_{\alpha>0}$ (see (4.93) and (4.94)) converges, as $\alpha \to \infty$, to a weak solution of the NS-TKE system, which satisfies the energy balance and thereby establishes the existence of a weak solution to the NS-TKE system obtained by approximation.

**Definition 5.1.** A pair $(\mathbf{v}, k) \in E_v \times K_{k,2}$ is a weak solution to the NS-TKE system (3.1)

1. $\partial_t \mathbf{v} \in L^2([0, T], V')$, $\mathbf{v} \in C([0, T], H)$ and

$$ (5.41) \quad \lim_{t \to 0} \|\mathbf{v}(t, \cdot) - \mathbf{v}_0\|_{\Omega;0,2} = 0. $$

$$ (5.42) \quad \forall \mathbf{w} \in E_v, \quad \langle \partial_t \mathbf{v}, \mathbf{w} \rangle + \langle B(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle + \langle A(k, \mathbf{v}), \mathbf{w} \rangle = \iint_{Q_T} \mathbf{f} \cdot \mathbf{w}. $$
2. $k \in L^\infty([0, T], L^1(\Omega))$ and for all $q \in \mathcal{E}_k$,

$$ (5.43) \qquad - \iint_{Q_T} \partial_t q(t, x) k(t, x) dxdt + \int_{\Omega} k_0(x)q(0, x)dx + \langle \tilde{B}(\mathbf{v}, k), q \rangle + \langle \tilde{A}(k, k), q \rangle = \langle \mathbb{D}(k, \mathbf{v}) - \varepsilon(k, k), q \rangle. $$

**Remark 5.1.** Let $(\mathbf{v}, k) \in E_v \times \mathcal{K}_{k,2}$ be a weak solution of the NS-TKE system. As already mentioned in Section 5.2.1, since $\mathbf{v} \in E_v$, it can be taken as a test vector field. As soon as we will have proved $\partial_t \mathbf{v} \in L^2([0, T], V')$, according to Theorem A.2 in Section A.2, we will deduce that $\mathbf{v} \in C([0, T], H)$ and by Formula (A.6), the energy balance is satisfied: for all $t \in [0, T]$,

$$ (5.44) \qquad \frac{1}{2} \int_{\Omega} |\mathbf{v}(t, x)|^2 dx + \iint_{Q_t} \nu_{\text{turb}}(k(s, x)) |\nabla \mathbf{v}(s, x)|^2 dxds = \frac{1}{2} \int_{\Omega} |\mathbf{v}_0(x)|^2 dx + \iint_{Q_t} \mathbf{f}(s, x) \cdot \mathbf{v}(s, x) dsdx. $$

**Remark 5.2.** Uniqueness remains an open problem. According to [6, 11] the conjecture is that if $\|\nu'_{\text{turb}}\|_\infty$ and $\|\mu'_{\text{turb}}\|_\infty$ are sufficiently small, then uniqueness can be proved.

We are now able to prove the following existence result in the two dimensional case.

**Theorem 5.1.** The NS-TKE system (3.1) admits a weak solution $(\mathbf{v}, k) \in E_v \times \mathcal{K}_{k,2}$ in the sense of Definition 3.2, which is constructed as a limit of the solution $\{\mathbf{v}_\alpha, k_\alpha\}_{\alpha>0}$ to the T-NS-TKE (4.3) as $\alpha \to \infty$.

We have to pass to the limit in Equations (4.93) and (4.94) as $\alpha \to \infty$. The proof follows the usual procedure and is divided into four steps:

1. Compactness of the family $\{\mathbf{v}_\alpha, k_\alpha\}_{\alpha>0}$ in suitable $L^p$ spaces, by means of the Aubin-Lions lemma and extraction of subsequences,

2. Passing to the limit, term by term, in the fluid equation (4.93),

3. Passing to the limit, term by term, in the TKE equation (4.94),

4. $L^\infty(L^1)$ bound for the TKE.

The proof of (5.41) concerning the initial data is the same as that of Item 4.3.v in Section 4.3.2. Therefore, we will not repeat it in what follows.

### 5.3.1 Compactness and extracting subsequences

We begin with the fluid equation (4.93). Recall that

$$ (5.45) \qquad V \subset H \subset V', $$

where the injections are dense and the first one is compact. As $\{\mathbf{v}_\alpha\}_{\alpha>0}$ is bounded in $E_v$, $L^2([0, T], V)$, and by (5.27), $\{\partial_t \mathbf{v}_\alpha\}_{\alpha>0}$ is bounded in $L^2([0, T], V)$, we deduce from the Aubin-Lions lemma, Lemma A.1, that $\{\mathbf{v}_\alpha\}_{\alpha>0}$ is compact in $L^2(Q_T)$. To summarize we have:

* $\{\mathbf{v}_\alpha\}_{\alpha>0}$ is weakly compact in $E_v$
• {$v_\alpha$}$_{\alpha>0}$ is strongly compact in $L^2([0, T], L^2(\Omega)) = L^2(Q_T)$,

• {$\partial_t v_\alpha$}$_{\alpha>0}$ is weakly compact in $L^2([0, T], V')$.

Let us consider the equation (4.94) for the TKE. We already know from the previous sections that

• The sequence {$k_\alpha$}$_{\alpha>0}$ is bounded in the space $L^p([0, T], W^{1,p}(\Omega))$, $1 \le p < 4/3$, as well as in $L^\infty([0, T], L^1(\Omega))$.

• The sequence {$\partial_t k_\alpha$}$_{\alpha>0}$ is bounded in $L^1([0, T], W^{-1,p}(\Omega))$, $1 \le p < 4/3$.

Let us consider a given $p < 4/3$, "close from 4/3". As

$$ (5.46) \qquad W_0^{1,p}(\Omega) \subset L^p(\Omega) \subset W^{-1,p}(\Omega), $$

where the injections being dense and the first one being compact. the Aubin-Lions lemma gives us:

• {$k_\alpha$}$_{\alpha>0}$ is weakly compact in $L^p([0, T], W_0^{1,p}(\Omega))$,

• {$k_\alpha$}$_{\alpha>0}$ is strongly compact in $L^p(Q_T)$,

Then there exists $(v, k) \in E_v \times K_2$ such that from {$v_\alpha, k_\alpha$}$_{\alpha>0}$ we can extract a subfamily still denoted {$v_\alpha, k_\alpha$}$_{\alpha>0}$ with

$$ v_\alpha \to v \quad \left\{ \begin{array}{ll} \text{weak in } E_v, & \\ \text{strong in } L^p(Q_T)^N, & 1 \le p < 4, \\ \text{almost everywhere (inverse Lebesgue Theorem) in } \Omega, & \end{array} \right. $$

$$ (5.47) \qquad \partial_t v_\alpha \underset{\alpha \to \infty}{\longrightarrow} \partial_t v \quad \text{weakly in } L^2([0, T], V'), $$

$$ k_\alpha \to k \quad \left\{ \begin{array}{ll} \text{weak in } L^p([0, T], W^{1,p}(\Omega)), & 1 \le p < \frac{4}{3}, \\ \text{strong in } L^q(Q_T), & 1 \le q < 2, \\ \text{almost everywhere in } \Omega. & \end{array} \right. $$

**Remark 5.3.** Let $p < q < 4/3$. Everything that has been done above can be repeated for the extracted subsequence {$k_\alpha$}$_{\alpha>0}$ in $L^q([0, T], W^{1,q}(\Omega))$, and a further subsequence can be extracted satisfying the convergence above. However, in this case, the limit is uniquely determined, so by a standard compactness argument the whole sequence converges. Therefore, everything written in (5.47) about {$k_\alpha$}$_{\alpha>0}$ still holds for all $1 \le p < 4/3$, and by (5.31) we have

$$ (5.48) \qquad k_\alpha \underset{\alpha \to \infty}{\longrightarrow} k \quad \text{strongly in } L^p(Q_T), \quad \forall p < 2. $$

### 5.3.2 Passing to the limit in the fluid equation.

The process of passing to the limit in the fluid equation (4.93) follows in the same way as in the case of the Galerkin method, except for the transport term. Indeed, arguing word by word as in Section 4.3.2, Item 4.3.iv, we see that

$$ Z_\alpha = \sqrt{\nu_{\text{turb}}(k_\alpha)} \nabla v_\alpha \underset{\alpha \to \infty}{\longrightarrow} Z = \sqrt{\nu_{\text{turb}}(k)} \nabla v \quad \text{weakly in } L^2(Q_T)^{N \times N}. $$

Therefore, $\forall \mathbf{w} \in E_v$,

$$ (5.49) \qquad \langle A(k_\alpha, v_\alpha), \mathbf{w} \rangle \underset{\alpha \to \infty}{\longrightarrow} \langle A(k, v), \mathbf{w} \rangle. $$

Moreover, we deduce from (5.20) that {$B_\alpha(v_\alpha, v_\alpha)$}$_{\alpha>0}$ is bounded in $L^2([0, T], V')$, since $\forall i = 1, 2, 3$, $|T_\alpha(v_{\alpha,i})| \le |v_{\alpha,i}|$. We have:
1. By Lemma A.5 in Section A.4, we have $T_\alpha(\mathbf{v}_\alpha) \xrightarrow{\alpha \to \infty} \mathbf{v}$ strongly in $L^2(Q_T)$.

2. Up to a subsequence, there exists $\mathcal{B} \in L^2([0, T], V')$ such that

$$B_{\alpha}(\mathbf{v}_{\alpha}, \mathbf{v}_{\alpha}) \xrightarrow{\alpha \to \infty} \mathcal{B}$$

weakly in $L^2([0, T], V')$.

It remains to check that $\mathcal{B} = B(\mathbf{v}, \mathbf{v})$, where $B(\mathbf{v}, \mathbf{v}) = \mathbf{v} \otimes \mathbf{v}$ was already defined in (5.3). We deduce from the above, combined with (5.47), that

$$T_{\alpha} \mathbf{v}_{\alpha} \otimes \mathbf{v}_{\alpha} \xrightarrow{\alpha \to \infty} \mathbf{v} \otimes \mathbf{v} \text{ strongly in } L^1(Q_T)^{N \times N},$$

which gives

$$(5.50) \qquad \mathcal{B} = B(\mathbf{v}, \mathbf{v}) \text{ in } \mathcal{D}'(Q_T) \text{ and therefore also in } L^2([0, T], V').$$

Consequently, (5.42) holds, which means

$$(5.51) \qquad \forall \mathbf{w} \in E_v, \quad \langle \partial_t \mathbf{v}, \mathbf{w} \rangle + \langle B(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle + \langle A(k, \mathbf{v}), \mathbf{w} \rangle = \iint_{Q_T} \mathbf{f} \cdot \mathbf{w}.$$

Regarding the initial data, arguing word by word as in Item 4.3.v of Section 4.3.2, we find that (5.41) holds. Moreover, according to Remark 5.1, the energy balance holds $\forall t \in [0, T]$,

$$(5.52) \qquad \frac{1}{2} \int_{\Omega} |\mathbf{v}(t, x)|^2 dx + \iint_{Q_t} \nu_{\text{turb}}(k(s, x)) |\nabla \mathbf{v}(s, x)|^2 dx ds = \frac{1}{2} \int_{\Omega} |\mathbf{v}_0(x)|^2 dx + \iint_{Q_t} \mathbf{f}(s, x) \cdot \mathbf{v}(s, x) dx ds.$$

From there, combining the energy balances (5.6) and (5.52), we can apply the energy method of Section 4.3.3, and we obtain:

$$(5.53) \qquad \nu_{\text{turb}}(k_\alpha) |\nabla \mathbf{v}_\alpha|^2 \xrightarrow{\alpha \to \infty} \nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2 \text{ strongly in } L^1(Q_T).$$

We deduce from Lemma A.5 that

$$(5.54) \qquad T_{\alpha}(\nu_{\text{turb}}(k_{\alpha}) |\nabla \mathbf{v}_{\alpha}|^2) \xrightarrow{\alpha \to \infty} \nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2 \text{ strongly in } L^1(Q_T).$$

### 5.3.3 Passing to the limit in the TKE equation.

Let $q \in \mathcal{E}_k$.

**Item 5.3.i)** Evolutionary term. The initial data is handled according to Remark 4.8, and passing to the limit in the equation (4.92). From the strong convergence of $\{k_\alpha\}_{\alpha>0}$ in $L^p(Q_T)$, $1 \le p < 4/3$, and that of $\{k_{\alpha,0}\}_{\alpha>0}$ to $k_0$ in $L^1(\Omega)$, we find

$$(5.55) \qquad \left\{ \begin{array}{l} \lim_{\alpha \to \infty} \iint_{Q_T} \partial_t q(t, x) k_\alpha(t, x) dx dt = \iint_{Q_T} \partial_t q(t, x) k(t, x) dx dt, \\ \\ \lim_{\alpha \to \infty} \int_{\Omega} k_{\alpha,0}(x) q(0, x) dx = \int_{\Omega} k_0(x) q(0, x) dx. \end{array} \right.$$

**Item 5.3.ii)** The transport term. We have
1. {$v_\alpha$}$_{\alpha>0}$ converges strongly to $v$ in $L^3(Q_T)^N$,

2. $\forall 1 \le p < 2$, {$k_\alpha$}$_{\alpha>0}$ converges strongly to $k$ in $L^p(Q_T)$ when $\alpha \to \infty$ (see Remark 5.3 and (5.48)), as well as {$T_\alpha(k_\alpha)$}$_{\alpha>0}$ by Lemma A.5.

Therefore, as $q \in C^\infty(Q_T)$, $\nabla q \in L^\infty(Q_T)^N$, then

$$
(5.56) \qquad T_{\alpha}(k_{\alpha})\mathbf{v}_{\alpha} \cdot \nabla q \xrightarrow[\alpha \to \infty]{} k\mathbf{v} \cdot \nabla q \quad \text{strongly in } L^1(Q_T),
$$

in particular,

$$
(5.57) \qquad \langle \tilde{B}_{\alpha}(k_{\alpha}, \mathbf{v}_{\alpha}), q \rangle \xrightarrow{\alpha \to \infty} \langle \tilde{B}(k, \mathbf{v}), q \rangle.
$$

**Item 5.3.iii) The dissipation term.** Again by (5.48), {$\sqrt{k_\alpha}$}$_{\alpha>0}$ converges to $\sqrt{k}$ in $L^p(Q_T)$ for all $1 \le p < 4$. Therefore,

$$
(5.58) \qquad \frac{k_{\alpha} \sqrt{k_{\alpha}}}{\ell} \xrightarrow{\alpha \to \infty} \frac{k \sqrt{k}}{\ell} \quad \text{strongly in } L^1(Q_T),
$$

which yields

$$
(5.59) \qquad \langle \varepsilon(k_{\alpha}, k_{\alpha}), q \rangle \xrightarrow{\alpha \to \infty} \langle \varepsilon(k, k), q \rangle.
$$

We also have from the energy method, in particular (5.54),

$$
(5.60) \qquad \langle \mathbb{D}_{\alpha}(k_{\alpha}, \mathbf{v}_{\alpha}), q \rangle \xrightarrow{\alpha \to \infty} \langle \mathbb{D}(k, \mathbf{v}), q \rangle.
$$

The diffusion term is handled in the same way as what we did in the Galerkin method,
using the continuity of $\mu_{\text{turb}}$, its boundness, the a.e. convergence of $\{k_\alpha\}_{\alpha>0}$ to $k$, since

$$
(5.61) \qquad \forall 1 \le p < 4/3, \quad \mathcal{E}_k \subset L^{p'}([0, T], W_0^{1,p'}(\Omega)).
$$

Therefore, for all $q \in \mathcal{E}_k$

$$
(5.62) \quad - \iint_{Q_T} \partial_t q(t, x) k(t, x) dxdt + \int_{\Omega} k_0(x)q(0, x)dx + \langle \tilde{B}(\mathbf{v}, k), q \rangle + \\
\qquad \qquad \qquad \langle \tilde{A}(k, k), r \rangle = \langle \mathbb{D}(k, \mathbf{v}), q \rangle - \langle \varepsilon(k, k), q \rangle.
$$

which was what we aimed to prove.

**5.3.4 L$^\infty$([0, T], L$^1(\Omega)$) estimate.**

It remains to prove that $k \in L^\infty([0, T], L^1(\Omega))$. We cannot directly apply the method of Section 3.4. We start from the estimate (5.8), which shows that the sequence {$k_\alpha$}$_{\alpha>0}$ is bounded in $L^\infty([0, T], L^1(\Omega))$. We already know from (5.47) that {$k_\alpha$}$_{\alpha>0}$ converges to $k$ in $L^1(Q_T)$. Let

$$
K_{\alpha}(t) = \int_{\Omega} |k_{\alpha}(t, x)| dx, \quad K(t) = \int_{\Omega} |k(t, x)| dx.
$$

Therefore, we deduce from Fubini's theorem that {$K_\alpha$}$_{\alpha>0}$ converges to $K$ in $L^1([0, T])$.
Hence, by the inverse Lebesgue theorem, we can extract from {$K_\alpha$}$_{\alpha>0}$ (still denoted by
{$K_\alpha$}$_{\alpha>0}$) a subsequence that converges a.e. to $K$. Then, by Fatou's lemma combined
with (5.8),

$$
\|k\|_{\infty; \Omega; 0,1} = \int_0^T K(t)dt \le \liminf_{\alpha \to \infty} \int_0^T K_\alpha(t)dt = \liminf_{\alpha \to \infty} \|k_\alpha\|_{\infty; \Omega; 0,1} \le \\
\|k_0\|_{\Omega; 0,1} + C(\|\mathbf{v}_0\|_{\Omega; 0,2}, \|\mathbf{f}\|_{Q_T; 0,2}, \nu, \beta)
$$

hence the result.

**Remark 5.4.** We conjecture that $k \in C([0, T], L^1(\Omega))$, which remains an open question.
## 5.4 Case $N = 3$

When $N = 3$, most of the arguments for passing to the limit in the sequence $\{\mathbf{v}_\alpha, k_\alpha\}_{\alpha>0}$ work as in the case $N = 2$, except that:

1. The critical exponent in the Boccardo-Galouët inequality is $\frac{5}{4}$.

2. From the estimates in $L^2([0, T], V)$ and $L^\infty([0, T], H)$ for $\{\mathbf{v}_\alpha\}_{\alpha>0}$, we also have by Lemma A.4, in particular (A.19) as $\rho(4) = 8/3$,

$$ (5.63) \qquad \{\mathbf{v}_\alpha\}_{\alpha>0} \text{ is bounded in } L^{8/3}([0, T], L^4(\Omega)), $$

which yields

$$ (5.64) \qquad \{\mathbf{v}_\alpha \otimes \mathbf{v}_\alpha\}_{\alpha>0} \text{ is bounded in } L^{4/3}([0, T], L^2(\Omega)). $$

As a consequence, we obtain

$$ (5.65) \qquad \{\partial_t \mathbf{v}_\alpha\}_{\alpha>0} \text{ is bounded in } L^{4/3}([0, T], V'). $$

In particular, we are not able to use Theorem A.2 for the limit $(\mathbf{v}, k)$ because of (5.65). Indeed, the dual product $\langle \partial_t \mathbf{v}, \mathbf{v} \rangle$ is not well defined since, a priori, $\mathbf{v} \notin L^4([0, T], V)$, and $\mathbf{v}$ cannot be taken as a test vector field in the limit fluid equation. Therefore, we cannot ensure that $\mathbf{v}$ satisfies the energy balance. As a consequence, we are not able to prove the convergence of the sequence $\{\nu_{\text{turb}}(k_\alpha)|\nabla \mathbf{v}_\alpha|^2\}_{\alpha>0}$ to $\nu_{\text{turb}}(k)|\nabla \mathbf{v}|^2$ in $L^1(Q_T)$, and we cannot apply the energy method of Section 4.3.3 anymore. This is why we cannot pass to the limit in the equation (4.94) for the TKE in the three-dimensional case. Hence, only very weak formulations are possible in this case, and the price to pay is the introduction of an additional non-negative measure in the right-hand side of the limit equation for the TKE.

### 5.4.1 Main result

We state the main result of this section about the NS-TKE system (3.1) when $N = 3$. Recall that the space $\mathcal{E}_k$ was defined in (5.40), and remains the same when $N = 3$. We will also use the following spaces:

$$ (5.66) \qquad \left\{ \begin{array}{l} \mathcal{K}_3 = \bigcap_{1 \le p < \frac{5}{4}} L^p([0, T], W^{1,p}(\Omega)), \\ \mathcal{E}_v = \{\mathbf{w} \in C^1([0, T], V) \text{ s.t. } \forall x \in \Omega, \mathbf{w}(T, x) = 0\}. \end{array} \right. $$

When $N = 3$, by passing to the limit when $\alpha \to \infty$, we obtain the following result.

**Theorem 5.2.** *There exists*

1. A non negative measure $m \in \mathcal{M}(Q_T)$,

2. A couple $(\mathbf{v}, k) \in E_v \times \mathcal{K}_3$, obtained as a limit of a subsequence of the sequence $\{\mathbf{v}_\alpha, k_\alpha\}_{\alpha>0}$, such that $\mathbf{v} \in C_w([0, T], H)$, $\partial_t \mathbf{v} \in L^{\frac{4}{3}}([0, T], V')$,

and such that $\forall (\mathbf{w}, q) \in \mathcal{E}_v \times \mathcal{E}_k$,

$$ (5.67) \qquad - \int_{Q_T} \partial_t \mathbf{w}(t, x) \cdot \mathbf{v}(t, x) dxdt - \int_{\Omega} \mathbf{w}(0, x) \mathbf{v}_0(x) dx + \\ \qquad \langle B_\alpha(\mathbf{v}, \mathbf{v}), \mathbf{w} \rangle + \langle A(\mathbf{v}, k), \mathbf{w} \rangle = \iint_{Q_T} \mathbf{f}(t, x) \cdot \mathbf{w}(t, x) dxdt, $$
and

$$
(5.68) \quad - \iint_{Q_T} \partial_t q(t, x) k(t, x) dxdt + \int_{\Omega} k_0(x)q(0, x)dx + \langle \tilde{B}(\mathbf{v}, k), q \rangle + \\
\qquad \qquad \qquad \qquad \langle \tilde{A}(k, k), r \rangle = \langle \mathbb{D}(k, \mathbf{v}), q \rangle - \langle \varepsilon(k, k), q \rangle + \langle m, r \rangle.
$$

Finally, the energy inequality holds, for almost all $t \in [0, T]$:

$$
(5.69) \quad \frac{1}{2} \int_{\Omega} |\mathbf{v}(t, x)|^2 dx + \iint_{Q_t} \nu_{\text{turb}}(k(s, x)) |\nabla \mathbf{v}(s, x)|^2 dxds \le \\
\qquad \qquad \qquad \qquad \frac{1}{2} \int_{\Omega} |\mathbf{v}_0(x)|^2 dx + \iint_{Q_t} \mathbf{f}(s, x) \cdot \mathbf{v}(s, x), dxds.
$$

**Remark 5.5.** As the measure *m* is non negative, we say that $(\mathbf{v}, k)$ is a weak subsolution of the NS-TKE system.

**Remark 5.6.** Recall that $\mathbf{v} \in C_w([0, T], H)$ means that $\mathbf{v}$ is weakly continuous from $[0, T]$ into $H$, namely for all $\mathbf{w} \in H$, the application $t \to (\mathbf{v}(t, \cdot), \mathbf{w}(\cdot))$ is continuous. The weak continuity in Theorem 5.2 is obtained by applying Theorem A.3, since $H \hookrightarrow V'$ with a continuous and dense injection, $\mathbf{v} \in L^\infty([0, T], H)$ and because $\partial_t \mathbf{v} \in L^{4/3}([0, T], V')$, then $\mathbf{v} \in C[0, T], V'$.

**Remark 5.7.** In this case, we cannot define the initial data for the fluid equation by $\lim_{t \to 0} \|\mathbf{v}(t, \cdot) - \mathbf{v}_0\|_{\Omega;0,2} = 0$, since we cannot repeat what was done in Item 4.3.v in Section 4.3.2. This is why we have chosen a distributional formulation for the evolutionary term in the fluid equation, following J. Leray [33]. However, we could also have used the weak continuity with values in $H$, and defined the initial data as:

$$
\forall \mathbf{w} \in H, \quad (\mathbf{v}(t, \cdot), \mathbf{w}(\cdot)) \xrightarrow{t \to 0} (\mathbf{v}_0, \mathbf{w}(\cdot)).
$$

The proof is the same as in the case $N=2$, carried out in Section 5.3, although the exponents are not the same (see, for example, (5.65)). The major difference between the cases $N=2$ and $N=3$ lies in the production term $\nu_{\text{turb}}(k_\alpha) |\nabla \mathbf{v}_\alpha|^2$. Therefore, to avoid unnecessary repetitions, in this proof we will only present how to derive the required estimates and how to perform the extraction of subsequences. Passing to the limit in the remaining terms of the equations follows the same pattern as in Section 5.3, except for the production term, which we detail below. We also explain how to obtain the weak continuity for $\mathbf{v}$ and the energy inequality (5.69).

To start with, we consider the velocity. Estimates (5.7) still hold. Moreover, we deduce
from (5.64)

$$
(5.70) \quad \partial_t \mathbf{v}_\alpha = \underbrace{-B_\alpha(\mathbf{v}_\alpha, \mathbf{v}_\alpha)}_{\text{bounded in } L^{4/3}([0,T], V')} - \underbrace{A(k_\alpha, \mathbf{v}_\alpha)}_{\text{bounded in } L^2([0,T], V')} + \underbrace{\mathbf{f}}_{\text{belongs to } L^2([0,T], V'(\Omega)^N)}
$$

Therefore, as already mentioned,

$$
(5.71) \qquad \{\partial_t \mathbf{v}_\alpha\}_{\alpha>0} \text{ is bounded in } L^{4/3}([0,T], V').
$$

Recall that

$V \subset H \subset V'$,

where the injections are dense, and the first one is compact. By the Aubin–Lions lemma, the sequence ${\mathbf{v}_\alpha}_{\alpha>0}$ is strongly compact in $L^2(Q_T)$, and therefore in $L^q(Q_T)$ for any $1 \le q < \frac{10}{3}$.
### 5.4.2 Additional estimates for the TKE equation.

Combining Lemma 5.1 with (5.8) and interpolation inequalities, we obtain

$$ (5.72) \qquad \left\{ \begin{array}{l} \forall 1 \le p < \frac{5}{4}, \quad \{k_\alpha\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{1,p}(\Omega)), \\ \{k_\alpha\}_{\alpha>0} \text{ is bounded in } L^\infty([0, T], L^1(\Omega)), \end{array} \right. $$

We have to estimate $\{\partial_t k_\alpha\}_{\alpha>0}$. To do so, each term of the TKE equation (4.94) is treated separately, in the same way as in the 2D case.

**Item 5.4.i) $L^\zeta(Q_T)$ estimate.** When $N=3$, by the Sobolev theorem we know that

$$ \forall 1 \le p \le \frac{15}{7}, \quad W^{1, \frac{5}{4}}(\Omega) \hookrightarrow L^p(\Omega). $$

Therefore, since we know that $\{k_\alpha\}_{\alpha>0}$ is bounded in $L^\infty([0, T], L^1(\Omega))$, we deduce from Lemma A.4 and formula (A.21), with $\zeta(5/4, 15/7, 1) = 5/3$ that

$$ (5.73) \qquad \forall 1 \le p < \frac{5}{3} = p_c, \quad \{k_\alpha\}_{\alpha>0} \text{ is bounded in } L^p(Q_T). $$

**Item 5.4.ii) Transport term.** Recall that

$$ \langle B_\alpha(\mathbf{v}_\alpha, k_\alpha), q \rangle = \iint_{Q_T} T_\alpha(k_\alpha)\mathbf{v}_\alpha \cdot \nabla q. $$

Since, for all $i=1,2,3$, $|T_\alpha(v_{i,\alpha})| \le |v_{i,\alpha}|$, we only need to determine in which $L^p$ space the product $k_\alpha \mathbf{v}_\alpha$ is bounded, using (5.7) and (5.72). Applying Hölder's inequality and noting that

$$ \underbrace{\frac{3}{10}}_{\text{critical exponent for } \mathbf{v}_\alpha} + \underbrace{\frac{3}{5}}_{\text{critical exponent for } k_\alpha} = \frac{9}{10} = \frac{1}{p_0} < 1, $$

we deduce that $\forall 1 \le p < p_0$, $\{k_\alpha \mathbf{v}_\alpha\}_{\alpha>0}$ is bounded in $L^p(Q_T)$. Hence, by Lemma A.3,

$$ (5.74) \qquad \forall 1 \le p < p_0, \quad \{\tilde{B}(k_\alpha, \mathbf{v}_\alpha)\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{-1,p}(\Omega)). $$

**Item 5.4.iii) Diffusion term.** Recall that

$$ \langle \tilde{A}(k_\alpha, k_\alpha), q \rangle = \iint_{Q_T} \mu_{\text{turb}}(k_\alpha) \nabla k_\alpha \cdot \nabla q. $$

Since $\mu_{\text{turb}}$ is bounded and (5.72) gives

$$ (5.75) \qquad \forall 1 \le p < \frac{5}{4}, \quad \{\nabla k_\alpha\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{1,p}(\Omega)^N), $$

we deduce

$$ (5.76) \qquad \forall 1 \le p < \frac{5}{4}, \quad \{\tilde{A}(k_\alpha, k_\alpha)\}_{\alpha>0} \text{ is bounded in } L^p([0, T], W^{-1,p}(\Omega)). $$

**Item 5.4.iv) Production term.** From the energy balance (5.6), $\{\nu_{\text{turb}}(k_\alpha)|\nabla \mathbf{v}_\alpha|^2\}_{\alpha>0}$ is bounded in $L^1(Q_T)$. Arguing as in Item 5.2.ix about the 2D case, but with the right exponent and noting that $\forall r > 5$, $W_0^{1,r}(\Omega) \hookrightarrow L^\infty(\Omega)$, we have

$$ (5.77) \qquad \forall 1 \le p < \frac{5}{4}, \quad \{\mathbb{D}_\alpha(k_\alpha, \mathbf{v}_\alpha)\}_{\alpha>0} \text{ is bounded in } L^1([0, T], W^{-1,p}(\Omega)). $$
**Item 5.4.v) Dissipation term.** We deduce from (5.73)

$$ \left\{ \frac{k_{\alpha} \sqrt{k_{\alpha}}}{\ell} \right\}_{\alpha>0} \text{ is bounded in } L^p(Q_T), \quad \forall 1 \le p < \frac{10}{9} = p_0, $$

hence

(5.78) $\forall 1 \le p < p_1$, $\{\varepsilon(k_\alpha, k_\alpha)\}_{\alpha>0}$ is bounded in $L^1([0, T], W^{-1,p}(\Omega))$.

**Item 5.4.vi) Time derivative.** We have

(5.79) $\partial_t k_\alpha = -\tilde{B}(k_\alpha, \mathbf{v}_\alpha) - \tilde{A}(k_\alpha, k_\alpha) + \mathcal{D}(k_\alpha, \mathbf{v}_\alpha) - \varepsilon(k_\alpha, k_\alpha)$.

Combining (5.74), (5.76), (5.77) and (5.78) we

(5.80) $\forall 1 \le p < \inf\{p_0, 5/4\} = p_0$, $\{\partial_t k_\alpha\}_{\alpha>0}$ is bounded in $L^1([0, T], W^{-1,p}(\Omega))$.

Let $1 \le p < p_0$, Since

$$ W_0^{1,p}(\Omega) \subset L^p(\Omega) \subset W^{-1,p}(\Omega), $$

with continuous dense injections and the first one compact, the Aubin-Lions lemma yields

(5.81) $\{k_\alpha\}_{\alpha>0}$ is compact in $L^p(Q_T)$.

### 5.4.3 Extracting subsequences and passing to the limit in the equations

We combine (5.7), (5.71), (5.72), (5.81) Remark 4.5 and Remark 5.3). Then there exists $(\mathbf{v}, k) \in E_v \times K_3$ and a subsequence, again denoted $\{\mathbf{v}_\alpha, k_\alpha\}_{\alpha>0}$, such that

$$ \mathbf{v}_\alpha \to \mathbf{v} \begin{cases} \text{weakly in } E_v, \\ \text{strongly in } L^p(Q_T)^N, \forall p < \frac{10}{3}, \\ \text{a.e. in } \Omega, \end{cases} $$

(5.82) $\partial_t \mathbf{v}_\alpha \xrightarrow[\alpha \to \infty]{} \partial_t \mathbf{v}$ weakly in $L^\frac{4}{3}([0, T], V')$,

$$ k_\alpha \to k \begin{cases} \text{weakly in } L^p([0, T], W^{1,p}(\Omega)), \forall 1 \le p < \frac{5}{4}, \\ \text{strongly in } L^p(Q_T), \forall p < p_c = \frac{5}{3}, \\ \text{a.e. in } \Omega, \end{cases} $$

The limiting process proceeds exactly as in the case $N=2$ (Section 5.3), with only the exponents changing, which determines the admissible test function spaces. However, special attention must be paid to the initial data in order to obtain (5.67). This issue is addressed in Remark 4.6, and it is then easy to pass to the limit in relation (4.69).

As already mentioned, **v** cannot be used as a test function in the limit equation. Therefore, the energy method no longer applies.

Since

$$ \{\sqrt{\nu_{\text{turb}}(k_\alpha)} \nabla \mathbf{v}_\alpha\}_{\alpha>0} \text{ converges weakly to } \sqrt{\nu_{\text{turb}}(k)} \nabla \mathbf{v} \text{ in } L^2(Q_T), $$

we obtain

(5.83) $\|\nu_{\text{turb}}(k)|\nabla\mathbf{v}|^2\|_{Q_T;0,1} \le \liminf_{\alpha \to \infty} \|\nu_{\text{turb}}(k_\alpha)|\nabla\mathbf{v}_\alpha|^2\|_{Q_T;0,1}$.
Thus, up to a subsequence, there exists a non-negative measure $m$ on $Q_T$ such that, in the sense of measures,

$$ (5.84) \qquad \lim_{\alpha \to \infty} \nu_{\text{turb}}(k_\alpha) |\nabla \mathbf{v}_\alpha|^2 = \nu_{\text{turb}}(k) |\nabla \mathbf{v}|^2 + m, $$

hence (5.68).

It remains to check the energy inequality (5.69). We start from the energy balance (5.6) satisfied by $(\mathbf{v}_\alpha, k_\alpha)$. It is easy to pass to the limit in the right-hand side of (5.6). We already have the inequality (5.83). It remains to deal with the kinetic energy term. Let $t \in [0, T]$, and let us define

$$ u_{\alpha}(t) = \frac{1}{2} \int_{\Omega} |\mathbf{v}_{\alpha}(t, x)|^2 dx, \quad u(t) = \frac{1}{2} \int_{\Omega} |\mathbf{v}(t, x)|^2 dx. $$

We have

$$ \begin{aligned} \int_{0}^{T} |u_{\alpha}(t) - u(t)| dt & \le \int_{0}^{T} \int_{\Omega} |(\mathbf{v}_{\alpha}(t, x) - \mathbf{v}(t, x)) \cdot (\mathbf{v}_{\alpha}(t, x) + \mathbf{v}(t, x))| dx dt \\ & \le \|\mathbf{v}_{\alpha} + \mathbf{v}\|_{Q_T;0,2} \|\mathbf{v}_{\alpha} - \mathbf{v}\|_{Q_T;0,2}. \end{aligned} $$

As $\mathbf{v}_\alpha \xrightarrow{\alpha \to \infty} \mathbf{v}$ in $L^2(Q_T)$, the term $\|\mathbf{v}_\alpha + \mathbf{v}\|_{Q_T;0,2}$ is bounded and $\|\mathbf{v}_\alpha - \mathbf{v}\|_{Q_T;0,2} \xrightarrow{\alpha \to \infty} 0$. Therefore, $u_\alpha \xrightarrow{\alpha \to \infty} u$ in $L^1([0, T])$. We deduce from the inverse Lebesgue theorem that, from the sequence $\{u_\alpha\}_{\alpha>0}$, we can extract a subsequence which converges to $u$ almost everywhere in $[0, T]$. Hence, by taking the limit in (5.6) for such times $t$, and combining this with (5.83), we conclude the proof of Theorem 5.2.

**Remark 5.8.** Building on Remark 4.9 concerning the boundedness of the eddy viscosities $\nu_{\text{turb}}$ and $\mu_{\text{turb}}$ (Assumption (3.3)), we point out that, in the three-dimensional case, we were able in [12] to remove the boundedness assumption on $\nu_{\text{turb}}$ and $\mu_{\text{turb}}$ when the velocity $\mathbf{v}$ satisfies Navier boundary conditions. In that setting, we obtain subsolutions for which an equation for the total energy can be derived. The result relies on a careful estimate of the pressure (see also [16, Chapter 8]), which no longer applies in the case of Dirichlet boundary conditions considered in the present work.

# Appendices

## A Tool box

This appendix collects some standard results from functional analysis and PDEs, used in this paper, in order to make it self-contained.

### A.1 Aubin-Lions Lemma

The Aubin-Lions lemma is a compactness result that is essential for proving the existence of solutions to system (3.1). Detailed proofs can be found in [2, 38, 50, 53]. Let $E$, $F$, and $G$ be Banach spaces such that

$$ E \hookrightarrow F \hookrightarrow G, $$

with continuous embeddings, and where $E \hookrightarrow F$ is also compact.

**Lemma A.1 (Aubin-Lions).** Let $1 \le r, p < \infty$, and let the sequence $(u_n)_{n \in \mathbb{N}}$ be bounded in $L^p([0, T], E)$, such that $(\partial_t u_n)_{n \in \mathbb{N}}$ is bounded in $L^r([0, T], G)$. Then one can extract from $(u_n)_{n \in \mathbb{N}}$ a subsequence that converges in $L^p([0, T], F)$.
## A.2 Time derivative and initial data

The proofs of the results in this section can be found in Temam [56].

**Theorem A.1.** Let $X$ be a Banach space, and let $u, g \in L^1([0, T], X)$. Then the following three conditions are equivalent.

i) The function $u$ is equal a.e. to a primitive of $g$,

$$ (A.1) \qquad \exists \xi \in X \text{ s.t. } \forall t \in [0, T], \quad u(t) = \xi + \int_0^t g(s)ds. $$

ii) $\forall \varphi \in D([0, T])$,

$$ (A.2) \qquad \int_0^T u(t)\varphi'(t)dt = -\int_0^T g(t)\varphi(t)dt. $$

iii) $\forall \eta \in X'$,

$$ (A.3) \qquad \frac{d}{dt}\langle\eta, u\rangle = \langle\eta, g\rangle \quad \text{in } D'([0, T]). $$

Moreover, if one of these three conditions is satisfied, $u$ is equal a.e. to a continuous function with values in $X$.

**Theorem A.2.** Let $V$ and $H$ two hilbert spaces such that

$$ (A.4) \qquad V \subset H \subset V', $$

each space is dense in the following one and the injections are continuous. Let $u \in L^2([0, T], V)$ such that $u' \in L^2([0, T], V')$. Then $u$ is a.e. equal to a function in $C([0, T], H)$ (we say $u \in C([0, T], H)$ so far there is no risk of confusion) and we have

$$ (A.5) \qquad \frac{d}{dt}\|u\|_H^2 = 2\langle u', u \rangle \quad \text{in } D'([0, T]). $$

**Remark A.1.** According to Theorem A.1, (A.5) is equivalent to

$$ (A.6) \qquad \forall t \in [0, T], \quad \|u(t)\|_H^2 = \|u(0)\|_H^2 + \int_0^T 2\langle u'(s), u(s) \rangle ds. $$

**Theorem A.3.** Let $X$ and $Y$ two Banach spaces such that $X \subset Y$ with continuous injections. Let $\phi \in L^\infty([0, T], X)$ continuous with values in $Y$. Then $\phi$ is weakly continuous with values in $X$.

## A.3 Boccardo–Gallouët lemma

The Boccardo–Gallouët lemma [8], is one of the main tools for dealing with elliptic and parabolic equations with a right-hand side in $L^1$ (or a measure), as is the case for the TKE equation.

Let $u : Q_T \to \mathbb{R}$ be a measurable function such that

$$ (A.7) \qquad |u| < \infty \quad \text{a.e. in } Q_T. $$
Let

$$
(A.8) \qquad \mathbb{B}_n(u) = \{(t, x) \in Q_T \text{ s.t. } n \le |u(t, x)| < n+1\}.
$$

Moreover, we assume that

$$
(A.9) \qquad \text{for almost all } t \in [0, T], \quad u(t, \cdot) \in \mathcal{D}'(\Omega),
$$

and

$$
(A.10) \qquad \forall n \in \mathbb{N}, \quad \nabla u \in L^2(B_n(u)).
$$

Let

$$
(A.11) \quad \forall n \in \mathbb{N}, \quad M_n(u) = \iint_{\mathbb{B}_n(u)} |\nabla u(t, x)|^2 dx dt.
$$

Finally, set

$$
(A.12) \qquad \mathbf{M}(u) = \sup_{n \in \mathbb{N}} M_n(u).
$$

**Lemma A.2** (Boccardo-Gallouët). Let $u : Q_T \to \mathbb{R}$ be a measurable function that satisfies (A.7), (A.8) and (A.10). Let $\mathbf{M}(u)$ be given by (A.12). Assume that

$$
(A.13) \qquad \mathbf{M}(u) < \infty.
$$

Then for all

$$
(A.14) \qquad 1 \le p < \frac{N+2}{N+1} = \tilde{p},
$$

$u \in L^p([0, T], W^{1,p}(\Omega))$ and there exists a constant $C_p = C(p, M_n(u), \Omega)$, such that

$$
(A.15) \qquad \|u\|_{p;\Omega;1,p} \le C_p.
$$

**Remark A.2.** We have

$$
(A.16) \qquad \lim_{\substack{p \to \tilde{p} \\ p < \tilde{p}}} C_p = +\infty.
$$

**A.4 Miscellaneous**

We start by the recall of a simple principle that will be useful in the sequel. The proof is
straightforward.

**Lemma A.3.** Let $F \in L^p([0, T], L^q(\Omega)^N)$. Then $F \in L^p([0, T], W^{-1,q}(\Omega))$, where

$$
(A.17) \quad \forall u \in L^{p'}([0, T], W_0^{1,q'}(\Omega)), \quad \langle F, u \rangle = \iint_{Q_T} F(t,x) \cdot \nabla u(t,x) \,dx\,dt,
$$

and moreover,

$$
(A.18) \qquad \|F\|_{p;\Omega; -1,q} \le \|F\|_{p;\Omega; 1,q}.
$$

We now give the interpolation principle used throughout the text.
**Lemma A.4.** Let $1 \le r < p < \infty$, $1 \le p$, and consider

$$u \in L^{\infty}([0, T], L^r(\Omega)) \cap L^p([0, T], L^q(\Omega)).$$

then $\forall r < s < q$,

$$(A.19) \qquad u \in L^{\rho(s)}([0, T], L^s(\Omega)), \quad \rho(s) = \frac{ps}{q} \left(\frac{q-r}{s-r}\right),$$

and we have

$$(A.20) \qquad \|u\|_{\rho(s); \Omega; 0, s} \le \|u\|_{\infty; \Omega; 0, r}^{\frac{\alpha}{\rho(s)}} \|u\|_{p; \Omega; 0, q}^{\frac{p}{\rho(s)}},$$

where

$$\alpha = \frac{rp(q-s)(s-r)}{(q-r)^2}.$$

Moreover, if $1 < p \le q$ then $u \in L^{\zeta}(Q_T)$ where

$$(A.21) \qquad \zeta = \zeta(p, q, r) = p + r \left(1 - \frac{p}{q}\right) > 1.$$

*Proof.* Let us write

$$s = \theta q + (1 - \theta)r, \quad \theta \in [0, 1].$$

By Hölder inequality, we have

$$(A.22) \qquad \int_{\Omega} |u|^s \le \left( \int_{\Omega} |u|^q \right)^{\theta} \left( \int_{\Omega} |u|^r \right)^{1-\theta} \le \|u\|_{\infty; \Omega; 0, r}^{(1-\theta)r} \|u\|_{\Omega; 0, q}^{q\theta}.$$

As

$$\theta = \frac{s-r}{q-r},$$

we deduce from (A.22)

$$\| u \|_{\Omega; 0, s}^{\rho(s)} \le \| u \|_{\infty; \Omega; 0, r}^{\alpha} \| u \|_{\Omega; 0, q}^{p},$$

hence (A.20) by integrating with respect to time. Formula (A.21) results from solving the equation $\rho(s) = s$. $\square$

We conclude this section with the following result. Recall that $T_\alpha$ is the truncation function at eight $\alpha$, Let $0 < \alpha$ and the truncation function $T_\alpha \in \mathbb{S}$ be given by

$$(A.23) \qquad T_{\alpha}: u \mapsto \begin{cases} -\alpha & \text{if } u < -\alpha, \\ u & \text{if } -\alpha \le u \le \alpha, \\ \alpha & \text{if } u > \alpha. \end{cases}$$

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