# Derivation and Implementation of a Non-Local Term to Improve the Oceanic Convection Representation Within the k-ε Parameterization

Alexandre Legay, Bruno Deremble, Hans Burchard

## ► To cite this version:

Alexandre Legay, Bruno Deremble, Hans Burchard. Derivation and Implementation of a Non-Local Term to Improve the Oceanic Convection Representation Within the k-ε Parameterization. Journal of Advances in Modeling Earth Systems, 2025, 17 (1), pp.e2024MS004243. 10.1029/2024MS004243. hal-05034814

HAL Id: hal-05034814

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

Submitted on 15 Apr 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2024MS004243

**Key Points:**

* Analytical derivation of a non-local term within the *k*–*ε* model to capture the effect of convective plumes

* In convective cases (free convection, diurnal cycle), the impact of the non-local term results in a better representation of the mixed layer

* For wind-driven mixing, the influence of the non-local term is negligible

**Correspondence to:**

A. Legay,
alexandre.legay@univ-grenoble-alpes.fr

**Citation:**

Legay, A., Deremble, B., & Burchard, H. (2025). Derivation and implementation of a non-local term to improve the oceanic convection representation within the *k*–*ε* parameterization. *Journal of Advances in Modeling Earth Systems*, 17, e2024MS004243. https://doi.org/10.1029/2024MS004243

Received 31 JAN 2024

Accepted 10 NOV 2024

© 2025 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# Derivation and Implementation of a Non-Local Term to Improve the Oceanic Convection Representation Within the *k*–*ε* Parameterization

Alexandre Legay¹ ID, Bruno Deremble¹ ID, and Hans Burchard² ID

¹CNRS, INRAE, IRD, Grenoble INP, IGE, Université Grenoble Alpes, Grenoble, France, ²Leibniz Institute for Baltic Sea Research Warnemünde, Rostock, Germany

**Abstract** The representation of turbulent fluxes during oceanic convective events is important to capture the evolution of the oceanic mixed layer. To improve the accuracy of turbulent fluxes, we examine the possibility of adding a non-local component in their expression in addition to the usual downgradient part. To do so, we extend the *k*–*ε* algebraic second-moment closure by relaxing the assumption on the equilibrium of the temperature variance θ̂². With this additional evolution equation for the temperature variance, we obtain a *k*–*ε*–θ̂² model (the “*k*et” model) which includes a non-local term for the temperature flux. We validate this new model against Large Eddy Simulations (LES) in three test cases: free convection (FC), wind-driven mixing, and diurnal cycle (DC). For wind-driven mixing, *k*et is equivalent to *k*–*ε*. However, in the presence of a buoyancy flux (FC and DC), we find that the vertical profile of temperature of the LES is better captured by *k*et than *k*–*ε*. Particularly, the non-local term increases the fraction of the mixed layer that is stably stratified. For FC, this fraction is near 50% for both *k*et and the LES, whereas the *k*–*ε* value is 20%. We show that this improvement is due to a better representation of the temperature variance in the inner part of the mixed layer. This better representation is mainly caused by the diffusion of temperature variance, which is described by *k*et and not by *k*–*ε*.

**Plain Language Summary** Cooling of the ocean surface creates dense cold water that tends to sink and vertically mix the water column. In numerical models, this vertical mixing is often represented by a diffusion of heat. This choice results in a “step-by-step” mixing from the ocean surface. However, in the case of strong surface cooling, the dense water parcels created at the top of the ocean sink rapidly, keeping their “coherence” during the descent. They only finally mix with the surrounding waters at depths of hundreds of meters. This cannot therefore be represented by a “step-by-step” mixing from the ocean surface. This kind of phenomenon is often referred to as coherent eddies or non-local turbulence. In this article, we perform an analytical derivation to give a mathematical expression of the impact of non-local turbulence. We then compare our new model with more realistic three-dimensional models of the ocean turbulence and conclude that the new term derived here is important to reproduce the vertical temperature profile in the surface ocean.

## 1. Introduction

In the realm of climate modeling, the oceanic mixed layer plays a critical role because it is responsible for regulating oceanic heat uptake and carbon storage. This storage depends on the maximum value of mixed layer depth reached during the year (see for example Luyten et al., 1983; Williams et al., 1995). These deepest mixed layers are often obtained in late winter, when ocean surface cooling has been frequent. Surface cooling triggers convective events in which convective plumes can mix the ocean over hundreds of meters and therefore put the mixed layer in direct contact with the deep ocean. Heat and carbon are then stored when the mixed layer becomes shallower and the direct contact with the atmosphere is lost, until same depths will be reached again at another time. Our representation of carbon and heat storage affects, for instance, our ability to make accurate predictions about future climate patterns (Treguier et al., 2023). This is just one example showing why accurately representing the mixed layer and the convective events is crucial in ocean models.

The depth of the mixed layer changes in response to various factors: it deepens when turbulent mixing is triggered by the mechanical effect of the wind or waves; or triggered by air-sea buoyancy fluxes: heat flux (cooling) or freshwater flux (evaporation, sea-ice formation). Conversely, the mixed layer becomes shallower typically when calm weather allows mixed layer instabilities to develop; or when there is a stabilizing buoyancy flux due to warming (e.g., sunny conditions) or freshwater input (e.g., precipitation, sea-ice melt, or river discharge). This
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Normalized profile from Large Eddy Simulations data of Mironov et al. (2000) of (a) the temperature and (b) the vertical turbulent temperature flux. The depth is normalized by the mixed layer depth $h_m$, defined here as the minimum of the temperature flux. The temperature flux is normalized by its surface value $w'\theta'|_0$. The temperature is normalized in $(\Theta - \Theta_{max})/\Theta^*$ with $\Theta_{max}$ the maximum of the temperature over the vertical and $\Theta^* = w'\theta'|_0/w^*$ a scaling of the temperature, with $w^* = (w'\theta'|_0 h_m)^{1/3}$ a scaling of the velocity of the convective plumes (Marshall & Schott, 1999; Willis & Deardorff, 1974). Red dashed lines highlight the location $h_t$ of the zero of the gradient $\partial_z\Theta$ and the location $h_f$ of the zero of the temperature flux.</figcaption></figure>

restratification allows the surface layer to separate from the denser, deeper water (Stull, 1988). In climate models, it is of course impossible to explicitly capture the small-scale turbulence responsible for the vertical mixing near the air-sea interface because the horizontal grid spacing is often on the order of tens of kilometers and the vertical grid spacing is on the order of meters. Instead, the ocean modeling community has developed parameterizations whose goal is to represent the mean effect of the turbulent fluctuations (Burchard & Bolding, 2001; Cheng et al., 2002; Fox-Kemper et al., 2008; Gaspar et al., 1990; Large et al., 1994; Reichl & Hallberg, 2018). The main purpose of a mixed layer parameterization is to propose a closure for the turbulent vertical fluxes $w'x'$, where $w'$ is the turbulent vertical velocity, $x'$ the turbulent fluctuation of a property $x$ (momentum, temperature, salinity, phytoplankton, etc...) and the overline denotes the ensemble averaging over small-scale fluctuations (see Stull, 1988). These turbulent fluxes, and all the other covariances $x'y'$ (with $y$ another property), are called the second-order moments. The traditional approach to close this problem consists of expressing these turbulent fluxes as a function of the vertical gradient of the mean property $X = \bar{x}$ (i.e., a downgradient parameterization), as shown here for the temperature

$$w'\theta' = -K_t\partial_z\Theta, \qquad (1)$$

with $K_t$ an eddy diffusivity coefficient. Among all the possibilities to compute $K_t$ we would like to emphasize the Generic Length Scale (GLS) approach (Umlauf & Burchard, 2003) and more precisely the $k-\epsilon$ closure (Burchard & Bolding, 2001; Umlauf & Burchard, 2005). This closure consists of deriving two equations: one for the evolution of turbulent kinetic energy (TKE) $k$, and one for dissipation $\epsilon$. In this formalism, the downgradient formulation (Equation 1) results from more complex algebraic second-moment closures even if it is not assumed a priori (Burchard & Baumert, 1995). The eddy diffusivity is obtained analytically and is a function of TKE, dissipation, buoyancy frequency, and shear frequency. While this eddy diffusivity approach has been successfully applied in the oceanic and atmospheric modeling communities, it has also been quickly recognized that the shape of the temperature profile during a convective event is not well captured by this closure. In fact, Deardorff (1972) was among the first to realize that after a convective event, the stratification profile in the mixed layer is not neutral as one would expect for a perfectly well-mixed layer but is instead slightly stable. To illustrate this observation, we plot in Figure 1a the typical shape of a normalized temperature profile in the mixed layer from a numerical model that explicitly resolves convection (see Mironov et al. (2000); details about the normalization are provided henceforth; we only wish to focus here on the shape of the temperature profile). This profile can be decomposed into two well-defined zones. Just below the air-sea interface, there is an unstable zone with cold
water above warmer water ($\partial_z\Theta < 0$). Such layer is sometimes called the thermal layer (Lazier, 2001) and we define it here as the layer between the surface and the depth $h_t$ at which $\partial_z\Theta = 0$. Below that depth $h_t$, we find the convective layer; a slightly stable layer that extends until the base of the mixed layer $h_m$. Both layers form the mixed layer. Different criteria for finding the mixed layer depth $h_m$ exist and, in Figure 1, we defined it as the position of the minimum of the temperature flux. The position of $h_t$ has been documented to be near $z = -0.4h_m$ (see Zhou et al., 2018) such that more than half of the mixed layer is stably stratified. The presence of such stable stratification in the convective layer has been attributed to downward propagating plumes which remain coherent during their descent and deposit their negative buoyancy anomaly at their neutral level, thus creating a stable stratification (see Arakawa and Schubert (1974) or Emanuel (1991) for the atmospheric scenario).

Several options have emerged in the literature to reproduce this vertical temperature profile with a stable stratification. The atmospheric community has favored the use of a mass flux parameterization which simulates the vertical movement of air parcels within convective clouds. The so-called “mass-flux” represents the ascent and descent of parcels, which transport heat, moisture, and other properties (Siebesma et al., 2007). These mass flux parameterizations have recently been introduced in ocean models (Garanaik et al., 2024; Giordani et al., 2020). Another perhaps older approach taken by Large et al. (1994) is to add a positive “non-local” term $\Gamma$ in the parameterization of the flux in Equation 1 (see also Troen and Mahrt (1986); or Burchard and Petersen (1999) where the problem of missing non-local fluxes in downgradient parameterization is stated):

$$ \overline{w'\theta'} = -K_t \partial_z\Theta + \Gamma. \qquad (2) $$

$\Gamma$ being positive, it represents a positive turbulent temperature flux, that is, a flux that follows the buoyancy effect (cold going down and hot going up). $\Gamma$ can thus be viewed as representing coherent structures (“non-local eddies,” “coherent plumes”) that are subjected to the buoyancy force. Particularly, we see in Equation 2 that $\Gamma$ allows to keep a positive turbulent temperature flux in situations of neutral ($\partial_z\Theta = 0$) or slightly stable ($\partial_z\Theta > 0$) temperature profiles. In other words, this means that, in stably stratified conditions, coherent structures can be strong enough to counter the downgradient flux that acts in a counter-buoyancy direction. Note that this term is sometimes written $\overline{w'\theta'} = -K_t(\partial_z\Theta - \gamma)$ with $\gamma = \Gamma/K_t$ (e.g., Deardorff, 1972; Large et al., 1994). In this formulation, $\gamma$ corresponds to the maximal stable stratification where a positive turbulent temperature flux can be maintained even if the downgradient flux generates a counter-buoyancy effect. In Large et al. (1994), $\Gamma$ is defined with some constraints: to be zero at the surface and at the base of the mixed layer such that it is merely a redistribution of heat. The magnitude and the exact shape of this term are however chosen in a relatively ad hoc way to respect some empirical rules of convection.

The term $\Gamma$ is often referred to in the literature as a “non-local” term (Ghannam et al., 2017; Large et al., 1994). As mentioned earlier, the denomination “non-local” refers to the fact that it is supposed to represent non-local eddies (coherent plumes). However, Zhou et al. (2018) argue that this often-implied association of the downgradient term with local turbulence and the $\Gamma$ term with non-local eddies is partially wrong. Another possibility is to call $\Gamma$ the “countergradient” term (Deardorff, 1972; Gibbs et al., 2011; Troen & Mahrt, 1986). This refers to the fact that, in the lower part of the mixed layer which is stable, this term acts with an opposite sign compared to the mean gradient. However, in the upper part of the mixed layer which is unstable, the denomination “countergradient” is incorrect since this term acts in a similar way as a downgradient term. We have opted here to keep the “non-local” denomination which is commonly used in the literature.

A key aspect of the addition of the non-local term is to relax the downgradient dependence and particularly the constraint that the depth $h_f$ at which $\overline{w'\theta'}$ vanishes is equal to the depth $h_t$ at which the gradient of the temperature profile vanishes (see Equation 1). To better understand why this matters, we plot in Figure 1b the vertical turbulent heat flux $\overline{w'\theta'}$ obtained in the same numerical model as presented before (Mironov et al., 2000). In this figure, we recover the traditional form of a linear decrease from the surface value (which corresponds to the magnitude of the surface flux) to a cancellation near the bottom of the mixed layer, which has been observed and described elsewhere (e.g., Burchard & Bolding, 2001; Large et al., 1994; Van Roekel et al., 2018). The exact depth at which the heat flux vanishes depends on the surface boundary conditions (wind and heat fluxes) but it has been documented to be close to $h_f = -0.8h_m$ in the free convection (FC) scenario (Garcia & Mellado, 2014). There is thus an obvious discrepancy between $h_t = -0.4h_m$ and $h_f = -0.8h_m$ such that Equation 1 cannot hold
in most of the mixed layer and the addition of an extra term in the definition of the flux is physically relevant. Even if there is a consensus on the need to add a non-local component in the definition of the flux, the exact formulation of this flux remains a matter of debate. To develop a framework that is accurate, robust, and consistent with existing parameterizations, we have opted to focus on extending the k-ε parameterization. There are at least two advantages of this approach compared to existing formulations of non-local terms. First, unlike mass-flux schemes where the interaction between the plumes and the environment depends on the entrainment/detainment parameters that are known to be difficult to parameterize, our analytical derivation ensures that the environment flux (downgradient part) aligns consistently with the representation of the plumes (non-local term). Second, the development in the k-ε framework allows us not to adopt the assumption of instantaneous equilibrium of the water column to surface forcings that is strongly embedded in the KPP framework, so that our formulation allows for memory effects in the development of turbulence in the boundary layer.

We first perform an analytical derivation of the non-local term. Since Deardorff (1972) and Cheng et al. (2020), we know that the non-local term is related to the small-scale temperature variance $\overline{\theta'^2}$. We will therefore derive a second-moment closure that uses a full transport equation for the temperature variance $\overline{\theta'^2}$, in addition to the second-moment transport equations for k and ε, thus extending the k-ε model to a k-ε-θ'^2 model (henceforth called the “ket” model). In this model, we get an analytical expression of a non-local term that shares several properties with the KPP non-local term: it is positive, and vanishes at the surface and at the bottom of the mixed layer. Last, we test the numerical implementation of ket against Large Eddy Simulations (LES) and further compare its results to the predictions of a standard k-ε model and KPP simulations.

## 2. Derivation and Implementation of the kεt Parameterization

This section introduces the second-order moments equations. We recall the hypotheses made in the GLS model to solve this system of equations. Then, we explain how we derive the ket parameterization in the same formalism.

### 2.1. Formalism and Second-Order Moments Equations

The Reynolds Averaged Navier Stokes (RANS) equations used in ocean models are written for the mean velocities $U = (U, V, W)$ and the mean temperature $\Theta$. As in the original derivation of the k-ε model, we consider here only one active tracer (temperature) that enters the equation of state. The RANS equations include the effect of turbulent fluctuations through the second-order moments $\overline{u'_i u'_j}$ and $\overline{u'_i \theta'}$; with $\overline{u'_i} = u'$, $\overline{u'_2} = v'$, and $\overline{u'_3} = w'$. To close the system, we need to provide equations for these moments. We focus here on the procedure derived in Burchard and Bolding (2001). After adopting their closure assumptions for non-closed terms, and neglecting the rotational and viscous effects, the equations of second-order moments are

$$ 
\begin{aligned}
\partial_t \overline{u'_i u'_j} + \partial_l (U_i \overline{u'_i u'_j} + \overline{u'_i u'_j} \overline{u'_j}) &= -c_1 \frac{\epsilon}{k} (\overline{u'_i u'_j} - \frac{2}{3} \delta_{ij} k) \\
&\quad + P_{ij} - c_2 \left( P_{ij} - \frac{2}{3} \delta_{ij} P \right) \\
&\quad + B_{ij} - c_3 \left( B_{ij} - \frac{2}{3} \delta_{ij} B \right) \\
&\quad - c_4 k S_{ij} \\
&\quad - c_5 Z_{ij} \\
&\quad - \frac{2}{3} \delta_{ij} \epsilon,
\end{aligned}
\quad (3) $$

$$ 
\begin{aligned}
\partial_t \overline{u'_i \theta'} + \partial_j (U_j \overline{u'_i \theta'} + \overline{u'_i u'_j \theta'}) &= -c_{1T} \frac{\epsilon}{k} \overline{u'_i \theta'} \\
&\quad - (1 - c_{2T}) \overline{u'_i \theta'} \partial_j U_i - \overline{u'_i u'_j} \partial_j \Theta \\
&\quad + (1 - c_{3T}) \beta_i \overline{\theta'^2} \\
&\quad + c_{4T} \overline{u'_j \theta'} V_{ij},
\end{aligned}
\quad (4) $$
<table><caption>**Table 1**
Values of the Coefficients Appearing in the Second-Order Moment Equations</caption><thead><tr><th>c<sub>1</sub></th><th>c<sub>2</sub></th><th>c<sub>3</sub></th><th>c<sub>4</sub></th><th>c<sub>5</sub></th><th>c<sub>1T</sub></th><th>c<sub>2T</sub></th><th>c<sub>3T</sub></th><th>c<sub>4T</sub></th><th>c<sub>T</sub></th></tr></thead><tbody><tr><td>2.5</td><td>0.984</td><td>0.5</td><td>0.512</td><td>0.416</td><td>5.95</td><td>0.6</td><td>0.33</td><td>0.4</td><td>1.44</td></tr></tbody></table>

$$ \partial_t \overline{\theta^2} + \partial_j (U_j \overline{\theta^2} + \overline{u'_j \theta^2}) = -2 \overline{u'_j \theta'} \partial_j \Theta - 2 \frac{\epsilon}{c_T} \frac{\epsilon}{k} \overline{\theta^2}, \quad (5) $$

with

* $P_{ij} = -\partial_t U_i \overline{u'_j u'_i} - \partial_t U_j \overline{u'_i u'_j}$: Production of $\overline{u'_i u'_j}$ by the shear

* $B_{ij} = \beta_i \overline{u'_j \theta'} + \beta_j \overline{u'_i \theta'}$: Production/destruction of $\overline{u'_i u'_j}$ by the buoyancy

* $S_{ij} = \frac{1}{2}(\partial_t U_j + \partial_j U_i)$: Shear tensor

* $V_{ij} = \frac{1}{2}(\partial_t U_j - \partial_j U_i)$: Vorticity tensor

* $Z_{ij} = V_{il}(\overline{u'_i u'_j} - \frac{2}{3}\delta_{ij}k) + V_{jl}(\overline{u'_j u'_i} - \frac{2}{3}\delta_{il}k)$: Symmetric tensor associated to the vorticity

* $k = \frac{1}{2}(\overline{u^2} + \overline{v^2} + \overline{w^2})$: TKE

* $P = \frac{1}{2}P_{ii}$: Production of TKE by the shear

* $B = \frac{1}{2}B_{ii}$: Production/destruction of TKE by the buoyancy

* $\epsilon$: Dissipation of TKE

Further definitions are $\delta_{ij}$ the Kronecker delta, $\beta = (0,0,\alpha g)$, $\alpha$ the thermal expansion coefficient and $g$ the gravitational acceleration. In all the previous equations and definitions, the Einstein summation convention is adopted. Note that buoyancy depends only on temperature. The inclusion of salinity and a non-linear equation of state add a significant level of complexity (Canuto et al., 2002) and is beyond the scope of this study (see also the mention of salinity in conclusion).

Closure assumptions for non-closed terms result in empirical coefficients $c_1, c_2, c_3, c_4, c_5$ for the parameterization of the pressure-velocity correlation tensor $\Pi_{ij} = \overline{u'_i \partial_j p} + \overline{u'_j \partial_i p}$, coefficients $c_{1T}, c_{2T}, c_{3T}, c_{4T}$ for the parameterization of the pressure-temperature correlations $\Pi_i^\theta = \overline{\theta' \partial_i p}$, and $c_T$ for the parameterization of the temperature variance dissipation. Further details about these parameterizations can be found in Canuto et al. (2001). We report the values of these coefficients in Table 1. These values are the ones of Canuto et al. (2001) model A, converted into the notations used here (it is the same as the values reported in Table 1 of Burchard and Bolding (2001) except for minor typos on $c_3$ and $c_4$ that have been identified. Exact formulations of these coefficients are given in Appendix D).

We are now going to explain the classic procedure used in the GLS models for closing the system, where the new model differs and what are the consequences.

## 2.2. GLS Procedure (k-ε Model)

The GLS procedure is as follows. First, we consider the boundary layer approximation where the vertical scale is much less than the horizontal scale. Horizontal gradients are then neglected in comparison to the vertical gradients. A direct consequence is the simplification of the continuity equation in $\partial_z W = 0$. The resulting expressions of the tensors $P_{ij}, B_{ij}, S_{ij}, V_{ij}$ and $Z_{ij}$ are given in Appendix E.

Second, we consider that the moments $\overline{u'_i \theta'}$ and $\overline{\theta'^2}$ are in local equilibrium, meaning that the sum of the time variations, the advective transports and the turbulent transports of these moments is zero (i.e., the left-hand sides of Equations 4 and 5 are zero). Concerning the moments $\overline{u'_i u'_j}$, the trick is to not make this assumption directly for $\overline{u'_i u'_j}$ but rather to the anisotropic part of these moments $\overline{u'_i u'_j} - 2/3 \delta_{ij}k$ to keep the time variation and the transports of the TKE to be non-zero. These assumptions correspond to the level $2 \frac{1}{2}$ in the hierarchy of models proposed by Mellor and Yamada (1982). This hierarchy has been derived with scaling arguments based on the level of anisotropy of every term. The scaling at level 3 results naturally in neglecting transports and time variations for $\overline{u'_i u'_j} - 2/3 \delta_{ij}k$ and $\overline{u'_i \theta'}$. However, neglecting these terms for the $\overline{\theta'^2}$ equation is not justified by the scaling process and is much more an ad hoc practical hypothesis that results in obtaining this so-called level $2 \frac{1}{2}$ in which the system of equations is now algebraic. Indeed, making all these assumptions on Equations 3–5 leads to the following set of equations
$$0 = -c_1 \frac{\epsilon}{k} \left( \overline{u'_i u'_j} - \frac{2}{3} \delta_{ij} k \right) + (1 - c_2) \left( P_{ij} - \frac{2}{3} \delta_{ij} P \right) + (1 - c_3) \left( B_{ij} - \frac{2}{3} \delta_{ij} B \right) - c_4 k S_{ij} - c_5 Z_{ij}, \quad (6)$$

$$0 = -c_{1T} \frac{\epsilon}{k} \overline{u'_i \theta'} - (1 - c_{2T}) \overline{u'_j \theta'} \partial_j U_i - \overline{u'_i u'_j} \partial_j \Theta + (1 - c_{3T}) \beta_i \overline{\theta'^2} + c_{4T} \overline{u'_j \theta'} V_{ij}, \quad (7)$$

$$0 = -2 \overline{u'_j \theta'} \partial_j \Theta - \frac{2}{c_T} \frac{\epsilon}{k} \overline{\theta'^2}, \quad (8)$$

which, if we assume that $k$ and $\epsilon$ are known, is a linear system of 10 equations with 10 unknowns: $(\overline{u'^2}, \overline{v'^2}, \overline{w'^2}, \overline{u'^v}, \overline{u'^w}, \overline{v'^w}, \overline{u'^\theta}, \overline{v'^\theta}, \overline{w'^\theta}, \overline{\theta'^2})$. For clarity, these 10 equations are written explicitly in Appendix F. We solved this system thanks to the symbolic calculus software Mathematica and we confirmed the expressions obtained by Burchard and Bolding (2001):

$$\overline{u'w'} = -f_m \frac{k^2}{\epsilon} \partial_z U, \quad (9)$$

$$\overline{v'w'} = -f_m \frac{k^2}{\epsilon} \partial_z V, \quad (10)$$

$$\overline{w'\theta'} = -f_h \frac{k^2}{\epsilon} \partial_z \Theta, \quad (11)$$

which reflect downgradient fluxes with an eddy viscosity $K_m = f_m \frac{k^2}{\epsilon}$ and an eddy diffusivity $K_t = f_h \frac{k^2}{\epsilon}$. The dimensionless functions $f_m$ and $f_h$ are the so-called “stability functions” and can be expressed in the following forms

$$f_m = \frac{n_0 + n_1 \alpha_N + n_2 \alpha_M}{d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2}, \quad (12)$$

$$f_h = \frac{n_{0T} + n_{1T} \alpha_N + n_{2T} \alpha_M}{d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2}, \quad (13)$$

with $\alpha_N = \frac{k_m^2}{\epsilon^2} N^2$, $\alpha_M = \frac{k_m^2}{\epsilon^2} M^2$, $N^2 = -g/\rho_0 \partial_z \rho = g\alpha \partial_z \Theta$ the (squared) buoyancy frequency, $M^2 = (\partial_z U)^2 + (\partial_z V)^2$ the (squared) shear frequency, and $\rho_0$ the reference density. Coefficients $n_i$, $n_{iT}$, and $d_i$ depend on the coefficients $c_i$ and $c_{iT}$. Their analytical expressions and numerical values are given in Appendix G.

To compute the fluxes in Equations 9–11, we still need to know $k$ and $\epsilon$. In a GLS model, we solve two prognostic equations, one for $k$ and one for another variable that can be linked to $\epsilon$. The choice of this second equation is the main difference between the different GLS models ($k-\epsilon$: Hanjalić and Launder (1972); Rodi (1987), $k-kl$: Mellor and Yamada (1982), $k-\omega$: Wilcox (1988), $k-\tau$: Zeierman and Wolfshtein (1986); Thangam et al. (1992)). In this paper, we focus on the $k-\epsilon$ model which solves directly the equation for $\epsilon$. The TKE equation and the $\epsilon$ equation are as follows

$$\partial_t k = P + G - \epsilon + D_k, \quad (14)$$

$$\partial_t \epsilon = \frac{\epsilon}{k} (c_{\epsilon 1} P + c_{\epsilon 3} G - c_{\epsilon 2} \epsilon) + D_{\epsilon}, \quad (15)$$

with

* $D_k = \partial_z (\frac{K_m}{\sigma_k} \partial_z k)$ and $D_\epsilon = \partial_z (\frac{K_m}{\sigma_\epsilon} \partial_z \epsilon)$: Diffusion terms

* $\sigma_k$ and $\sigma_\epsilon$: Schmidt numbers for TKE and dissipation

* $P \equiv (-\overline{u'w'} \partial_z U - \overline{v'w'} \partial_z V) = f_m \alpha_M \epsilon$: Production of TKE by the shear

* $G \equiv \beta_3 w' \theta' = -f_h \alpha_N \epsilon$: Production/destruction of TKE by the buoyancy

* $c_{\epsilon 1}$, $c_{\epsilon 2}$ and $c_{\epsilon 3}$: Empirical coefficients
The TKE Equation 14 was obtained by taking the trace of the Reynolds stress Equation 3. With the boundary layer approximation which neglects the horizontal gradient in comparison to the vertical ones, taking this trace gives $\partial_t k + \frac{1}{2}(\partial_z w' u'_i u'_i) = P + G - \epsilon$. We then consider downgradient formulations for the third-order moments $w'u'_i u'_i$ and we finally obtain Equation 14. We want to highlight that the diffusion term thus comes from the divergence of the third-order moments.

An exact equation for $\epsilon$ can be derived but, in practice, this equation needs additional assumptions to be closed. The classic assumptions discussed in the literature (well detailed in equations 3.1–3.5 of Hanjalić and Launder (1972)) are equivalent to scaling the sources and sinks of $\epsilon$ with the ones of the TKE through empirical coefficients $c_{\epsilon1}$, $c_{\epsilon2}$, and $c_{\epsilon3}$ (see Burchard & Bolding, 2001). This gives Equation 15.

Values $\sigma_k = 1$, $c_{\epsilon1} = 1.44$ and $c_{\epsilon2} = 1.92$ are frequently used in the literature (Rodi, 1987). Value $\sigma_\epsilon = 1.20$ is found according to Equation 14 of Umlauf and Burchard (2003). Finally, for $c_{\epsilon3}$, it is common practice to consider two different values in order to keep $c_{\epsilon3}G$ always as a source term of $\epsilon$ (Burchard & Bolding, 2001; Reffray et al., 2015; Rodi, 1987; Umlauf & Burchard, 2003; Warner et al., 2005). A positive value $c_{\epsilon3}^+ = 1$ is used when $G$ is positive (unstable stratification) and a negative value $c_{\epsilon3}^- = -0.65$ is used when $G$ is negative (stable stratification). This value $c_{\epsilon3}^- = -0.65$ is obtained according to Equation 26 of Umlauf et al. (2003) (by considering a steady-state Richardson number equal to 0.25).

## 2.3. Procedure for the k\epsilon Model

The new procedure differs from the GLS one by considering that the temperature variance $\overline{\theta^2}$ is not at equilibrium anymore. Relaxing this assumption takes us from the level 2 1/2 to the level 3 in the hierarchy of Mellor and Yamada (1982). Beyond this mathematical justification, the idea of keeping the non-equilibrium $\overline{\theta^2}$ equation originates from the fact that the $\overline{\theta^2}$ dependence appears only in the $\overline{w'\theta'}$ equation (see Equations 3 and 4). Thus, a physical change in the shape of the $\overline{\theta^2}$ profile will directly impact $\overline{w'\theta'}$. Because we now have an equation for the temperature variance, we are left with Equations 6 and 7 that form a system of nine equations with nine unknowns: $(\overline{u'^2}, \overline{v'^2}, \overline{w'^2}, \overline{u'v'}, \overline{u'w'}, \overline{v'w'}, \overline{u'\theta'}, \overline{v'\theta'}, \overline{w'\theta'})$. For clarity, these nine equations are written explicitly in Appendix H. We solve this system with Mathematica and we obtain

$$ \overline{u'w'} = -f_m \frac{k^2}{\epsilon} \partial_z U, \qquad (16) $$

$$ \overline{v'w'} = -f_m \frac{k^2}{\epsilon} \partial_z V, \qquad (17) $$

$$ \overline{w'\theta'} = -f_h \frac{k^2}{\epsilon} \partial_z \Theta + f_h^* \frac{k}{\epsilon} \beta_3 \overline{\theta^2}. \qquad (18) $$

The momentum fluxes are still downgradient with an eddy viscosity $K_m = f_m \frac{k^2}{\epsilon}$ whereas the temperature flux now has a “non-local” contribution $\Gamma_{ket} = f_h^* \frac{k}{\epsilon} \beta_3 \overline{\theta^2}$ related to the temperature variance in addition to the downgradient part with eddy diffusivity $K_t = f_h \frac{k^2}{\epsilon}$. The stability functions $f_m$, $f_h$ and $f_h^*$ can be expressed in the following forms

$$ f_m = \frac{n_0 + n_1 \alpha_N + n_2 \alpha_M + n_3 \alpha_T}{d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2}, \qquad (19) $$

$$ f_h = \frac{n_{0T} + n_{1T} \alpha_N + n_{2T} \alpha_M}{d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2}, \qquad (20) $$

$$ f_h^* = \frac{n_{0T}^* + n_{1T}^* \alpha_N + n_{2T}^* \alpha_M}{d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2}, \qquad (21) $$

with $\alpha_N = \frac{k^2}{\epsilon^2} N^2$, $\alpha_M = \frac{k^2}{\epsilon^2} M^2$, and $\alpha_T = \frac{k}{\epsilon^2} \beta_3^2 \overline{\theta^2}$. Coefficients $n_i$, $n_{iT}$, $n_{iT}^*$, and $d_i$ depends on the coefficients $c_i$ and $c_{iT}$. Their analytical expressions and numerical values are given in Appendix I.
As in the GLS procedure, the TKE and $\epsilon$ equations (Equations 14 and 15) are solved prognostically. The only differences in these equations are about the terms G and P that differ through the new expressions of the stability functions $f_m$ and $f_h$, and by the additional influence of $f_h^*$ in $G = -f_h \alpha_N \epsilon + f_h^* \alpha_T \epsilon$.

The key change of the $k\epsilon t$ model is that the temperature variance is now also solved prognostically through:

$$ \partial_t \overline{\theta'^2} = -2 \overline{w'\theta'} \partial_z \Theta - \frac{2}{c_T} \frac{\epsilon}{k} \overline{\theta'^2} + D_{\overline{\theta'^2}} \quad (22) $$

with $D_{\overline{\theta'^2}} = \partial_z \left( \frac{K_m}{\sigma_{\overline{\theta'^2}}} \partial_z \overline{\theta'^2} \right)$ the diffusion and $\sigma_{\overline{\theta'^2}}$ the Schmidt number for the temperature variance. As for the TKE equation, the diffusion term $D_{\overline{\theta'^2}}$ results from the closure of the third-order moment $\overline{w'\theta'\theta'}$ by a downgradient formulation (Moeng & Wyngaard, 1989). Using LES data, we have compared $\overline{w'\theta'\theta'}$ and $\frac{K_m}{\sigma_{\overline{\theta'^2}}} \partial_z \overline{\theta'^2}$ and found that $\sigma_{\overline{\theta'^2}}$ is order unity (not shown). We thus opted for taking $\sigma_{\overline{\theta'^2}} = 1$.

We can add three general remarks before closing this section. The expression of $\Gamma_{ket} = f_h^* \frac{k}{\epsilon} \beta_3 \overline{\theta'^2}$ obtained for the temperature flux can be compared to the one found by Deardoff (1972). By reasoning with the $\overline{w'\theta'}$ equation, Deardoff (1972) finds a non-local term $\Gamma_{Deardoff} \propto l/k^{1/2} \overline{\theta'^2}$ with l a mixing length introduced for the parame- terization of the pressure-temperature correlation. If we consider the classic scaling $l \propto k^{3/2}/\epsilon$ (see for example Rodi, 1987; Umlauf & Burchard, 2003, 2005), we obtain $\Gamma_{Deardoff} \propto k/\epsilon \overline{\theta'^2}$. The non-local expressions of $\Gamma_{ket}$ and $\Gamma_{Deardoff}$ thus both exhibit the same dependence on the turbulence time scale $k/\epsilon$ and on the temperature variance $\overline{\theta'^2}$.

In a $k-\epsilon$ model, $\overline{\theta'^2}$ is at equilibrium (i.e., it follows Equation 22 with $\partial_t \overline{\theta'^2} = 0$ and $D_{\overline{\theta'^2}} = 0$). It is possible to diagnose the part of $\Gamma_{ket}$ that is associated with this equilibrium value and therefore implicitly present in the $k-\epsilon$ model. This point is detailed in Appendix A.

Finally, we point out that, just as we retained the non-equilibrium equation of $\overline{\theta'^2}$ to obtain a non-local term for $\overline{w'\theta'}$, it would be tempting to retain the non-equilibrium equation of $\overline{w'^2}$ to obtain non-local terms for the velocity fluxes $\overline{u'w'}$ and $\overline{v'w'}$. We solved this problem and found that the velocity fluxes $\overline{u'w'}$ and $\overline{v'w'}$ in this context are still downgradient. Results of this $k-\epsilon-\overline{\theta'^2}-\overline{w'^2}$ model are detailed in Appendix C.

## 2.4. 1D Models Simulations

We implemented the $k\epsilon t$ parameterization, with the formalism described in Section 2.3, in the 1D code presented in Fearon et al. (2020). This code is a standalone 1D vertical version of the Coastal and Regional Ocean Community model (CROCO, https://www.croco-ocean.org/) and allows to run simulations with KPP, TKE, and several GLS schemes (note that we also re-implemented the $k-\epsilon$ model with the formalism presented in Section 2.2, that is equivalent to using the Canuto et al. (2001) stability functions).

In this code, the marching steps are as follows. For a given time step, we first solve the equations for the TKE and $\epsilon$. We then calculate the associated diffusion coefficients $K_m$ and $K_t$. Finally, the equations for the mean variables $(U, V, \Theta)$ are solved. For the implementation of the $k\epsilon t$ model, we added the resolution of the $\overline{\theta'^2}$ equation at the same time as the TKE and the $\epsilon$ equations.

For the time-discretization of the temperature variance Equation 22, the $k, \epsilon, K_t, K_m, f_h^*$, and $\partial_z \Theta$ values used are the ones of the previous time step (explicit expression). For the terms involving $\overline{\theta'^2}$, the following choices are made:

$$ \frac{1}{\Delta t} \left( \underbrace{\overline{\theta'^2}}_{n+1} - \underbrace{\overline{\theta'^2}}_{n} \right) = 2 K_t (\partial_z \Theta)^2 - \beta_3 f_h^* k/\epsilon \partial_z \Theta \underbrace{\overline{\theta'^2}}_{\text{see text}} - \frac{2}{c_T} \frac{\epsilon}{k} \underbrace{\overline{\theta'^2}}_{n+1} + \partial_z \left( \frac{K_m}{\sigma_{\overline{\theta'^2}}} \underbrace{\partial_z \overline{\theta'^2}}_{n+1} \right) \quad (23) $$
<table><caption>**Table 2**
Summary of the Parameters Used for the Three Test Cases</caption><thead><tr><th>Test case</th><th>FC</th><th>W</th><th>DC</th></tr></thead><tbody><tr><td>Duration</td><td>7 days</td><td>30 hr</td><td>10 days</td></tr><tr><td>Q<sub>NonSol</sub> (W/m<sup>2</sup>)</td><td>-75</td><td>0</td><td>-75</td></tr><tr><td>Q<sub>Sol</sub> (W/m<sup>2</sup>)</td><td>0</td><td>0</td><td>See text</td></tr><tr><td>τ<sup>w</sup> (Pa)</td><td>0</td><td>0.1</td><td>0</td></tr><tr><td>dθ/dz|<sub>0</sub> (K/m)</td><td>0.01</td><td>0.05</td><td>0.01</td></tr><tr><td>N<sup>2</sup><sub>0</sub> (s<sup>-2</sup>)</td><td>2 × 10<sup>-5</sup></td><td>1 × 10<sup>-4</sup></td><td>2 × 10<sup>-5</sup></td></tr></tbody></table>

Note. All three start with an initial homogeneous stable stratification $\partial\Theta/\partial z|_{0}$ (the associated $N_{0}^{2}$ is also given) and a surface temperature $\Theta_{surf} = 293.15$ K.

The diffusion is expressed implicitly because of the possibly large value of the diffusivity coefficient which would require small time steps. For the other terms, the Patankar trick is used to preserve the positivity of $\overline{\Theta}^{2}$ (Burchard, 2002; Lemarié et al., 2021; Patankar, 1980). This trick consists of expressing implicitly all the negative terms of the right-hand side (such as the dissipation term). The second right-hand side term is thus expressed explicitly if $\partial_{x}\Theta < 0$ and implicitly if $\partial_{x}\Theta > 0$. Boundary conditions for the temperature variance are zero at the bottom of the domain (Dirichlet condition), while at the surface a homogeneous Neumann condition is used for simplicity (no flux of temperature variance).

For every test case, we performed the simulations using $k-\epsilon$, $k\epsilon$, and KPP. The changes induced by $k\epsilon$, particularly the influence of the non-local term, will be analyzed by comparing with $k-\epsilon$. Concerning KPP, the simulations were done with and without its non-local term. The goal is to compare this term and its effect to the non-local term obtained in $k\epsilon$. The version of KPP used here is the original one described in Large et al. (1994).

## 2.5. LES

In order to validate the $k\epsilon$ model, we performed several LES. Practically, we use the Basilisk code (http://basilisk.fr, Popinet, 2020) to solve the three-dimensional Boussinesq equations in a small oceanic patch near the air-sea interface. We intend to explicitly compute the turbulent fluxes and the mean vertical profiles of temperature for buoyancy-driven convection and wind-driven turbulence. The domain is periodic in the horizontal direction. All variables are discretized at the cell center and are advected using the Bell-Collela-Glaz method. There is no explicit viscosity and no explicit diffusivity: both these terms are handled implicitly by the advection scheme. The surface forcing (wind and heat flux) is applied at the upper grid cell with a relaxation term. The bottom boundary condition is free slip for the velocity and inhomogeneous Neumann for the temperature (set to the initial stratification). The model is initialized with zero velocity and prescribed stratification for temperature (see next paragraph) to which we add a small random perturbation of magnitude $10^{-3}$ K. We use an adaptive time step adjusted with a CFL condition of 0.6. Averages are computed in a post-processing step: the overbar is interpreted here as a horizontal average and primes are deviations from this horizontal average. No temporal averaging is done.

## 2.6. The Three Test Cases

We evaluate the performance of the $k\epsilon$ parameterization with three different idealized surface forcing test cases: FC, wind-driven mixing (W), and diurnal cycle (DC). These cases are similar to the ones proposed by Van Roekel et al. (2018). The parameters for these three simulations are detailed in Table 2. For all three simulations, we take a linear equation of state $\rho = \rho_{0}[1 + \alpha(\Theta - \Theta_{0})]$ with $\rho_{0} = 1,025$ kg/m$^3$ the reference density, $\alpha = 2 \times 10^{-4}$ K$^{-1}$ the thermal expansion coefficient, and $\Theta_{0} = 293.15$ K the reference temperature. Rotational effects are not included. The LES domain size is $L_{x} = L_{y} = 100$ m, and $L_{z} = 50$ m (W, DC) or $L_{z} = 100$ m (FC). The grid resolution is isotropic (0.23 m) with $512 \times 512 \times 256$ (W, DC) or $512 \times 512 \times 512$ (FC) grid cells. The choice of the aspect ratio $L_{z}/L_{x}$ affects the convergence of the horizontal averaging, it is hence a trade-off between computation cost and quality of the resulting mean profiles. For the 1D simulations, the same vertical grids are used, and the time step is 360 s.

The temporal evolution for the DC test case of the net surface heat flux $Q_{0} = Q_{NonSol} + Q_{Sol}$ is made of $Q_{NonSol} = -75$ W m$^{-2}$ and a solar radiation $Q_{Sol}$ written as

$$Q_{Sol} = 235.62 \max[\cos(2\pi t/86400 - \pi); 0]. \qquad (24)$$

This solar radiation $Q_{Sol}$ penetrates through the water column according to a two-band exponential formulation with coefficients of a Jerlov type IB water mass (Paulson & Simpson, 1977). A graphical representation of $Q_{0}$ over a day is given in the results section in Figure 5b.
For FC and wind-driven mixing, the results will be compared with analytical solutions of the time evolution of $h_m$. For FC, the following evolution laws have been proposed (Marshall & Schott, 1999; Souza et al., 2020; Van Roekel et al., 2018)

$$h_m = \sqrt{2 Q_0 t/N_0^2} \quad (\text{Marshall & Schott, 1999}), \quad (25)$$

$$h_m = \sqrt{2.8 Q_0 t/N_0^2} \quad (\text{Van Roekel et al., 2018}), \quad (26)$$

$$h_m = \sqrt{3 Q_0 t/N_0^2} \quad (\text{Souza et al., 2020}), \quad (27)$$

with $Q_0$ the net surface heat flux, $N_0$ the initial buoyancy frequency, and $t$ the time. All three formulas have been confirmed by LES in their respective articles, so there is a priori no clue as to which one to choose. Van Roekel et al. (2018) shows that the prefactor 2.8 corresponds to the empirical rule of convection $\overline{w'\theta'}(h_m)/\overline{w'\theta'}_0 = -0.2$ (cf., Figure 1b and for example Large et al., 1994). Differences between the three formulas then likely indicate that the value -0.2 is not universal and depends on the intensity of the forcings, as suggested by Haghshenas and Mellado (2019). Here, for simplicity, we retain the prefactor 2.8 as the value $\overline{w'\theta'}(h_m)/\overline{w'\theta'}_0 = -0.2$ is the most commonly considered.

For wind-driven deepening, the following laws have been proposed (Kato & Phillips, 1969; Pollard et al., 1973; Price, 1979)

$$h_m = u_* (15 t/N_0^2)^{1/3} \quad (\text{Kato & Phillips, 1969}), \quad (28)$$

$$h_m = u_* (\sqrt{2} t/N_0)^{1/2} \quad (\text{Pollard et al., 1973}), \quad (29)$$

$$h_m = u_* (\sqrt{1.2} t/N_0)^{1/2} \quad (\text{Price, 1979}), \quad (30)$$

with $u_* = \sqrt{|\tau^w|/\rho_0}$ the surface wind friction velocity and $|\tau^w|$ the norm of the wind stress vector at the ocean surface. We decided to keep the expression of Price (1979) because it relies on constructive criticism of the formulas proposed by Kato and Phillips (1969) and Pollard et al. (1973) (see also Deleersnijder & Luyten, 1994).

Note that one limitation of these three test cases is that they do not take waves and Langmuir turbulence into account. Some suggestions for including these processes in $k-\epsilon$ have been proposed (e.g., for waves: Mellor and Blumberg (2004); e.g., for Langmuir: Axell (2002); Harcourt (2015)) and could be directly implemented in $k\epsilon$. Adding these features and analyzing the behavior of $k\epsilon$ in such test cases will be considered in future work.

## 2.7. Nondimensionalization

In order to compare the shape of the different profiles, variables are made dimensionless. The non-dimensional temperature is computed as $(\Theta - \Theta_{max})/\Theta^*$, with $\Theta_{max}$ the maximum of the temperature over the vertical, $\Theta^* = \overline{w'\theta'}_0/w^*$ a scaling of the temperature, $w^* = (-B_0h_m)^{1/3}$ a scaling of the velocity of the convective plumes (Marshall & Schott, 1999; Willis & Deardorff, 1974), and $B_0 = g\alpha Q_0/(\rho_0 c_p)$ the surface buoyancy flux. This normalization is valid only if $B_0$ is strictly negative and it is thus not used for the W and the DC test cases. The temperature flux is normalized by its surface value $\overline{w'\theta'}_0 = Q_0/(\rho_0 c_p)$. The depth is normalized by $h_m$. For the FC and the W test cases, $h_m$ is diagnosed as the depth of the maximum of $\partial_z\Theta$. For the LES, we rather take the mean depth of the 10 highest values because the profiles, which are instantaneous, can be noisy. For the DC test case, a consequence of the non-stationary surface forcing is that the profiles can exhibit non-canonical shapes (compared to the standard FC and W cases). To compute $h_m$, definitions based on integrated properties are then better suited than the ones based on thresholds or extrema. We choose to use the criterion proposed by Reichl et al. (2022) which is based on energetic principles. More precisely, the mixed layer depth is defined as the depth until which a given energy $\phi$ could homogenize the surface ocean. We take the value $\phi = 1 J/m^2$ and we perform
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Temporal evolution of the mixed layer depth *h*<sub>m</sub> for (a) free convection, (b) W, and (c) diurnal cycle test cases. The first 2 hr were removed in each plot to remove noisy *h*<sub>m</sub> which correspond to the initial mixed layer establishment.</figcaption></figure>

the calculation thanks to the python package "oceanmixedlayers" developed by the authors (Reichl et al., 2022, https://github.com/breichl/oceanmixedlayers).

**3. Results and Discussion**

**3.1. Mixed Layer Depth Evolutions**

The *k*–*ε* and the *k*et models exhibit comparable mixed layer depths *h*<sub>m</sub> in all three test cases (Figure 2). For W, both *k*–*ε* and *k*et give *h*<sub>m</sub> close to the analytical solution (Equation 30) whereas the LES exhibits *h*<sub>m</sub> which follows the analytical solution for half a day. After 13 hr, the LES predicts a deeper mixed layer (maybe due to the onset of internal waves breaking). For DC, both the *k*–*ε* and the *k*et models result in *h*<sub>m</sub> which is close to the LES and the analytical solution. To sum up, the *k*–*ε* and the *k*et models predict similar *h*<sub>m</sub>, and this *h*<sub>m</sub> reproduces well analytical solutions or the LES. The next sections will describe in more detail all three test cases to highlight the changes observed with the *k*et model in comparison to the *k*–*ε* scheme.

**3.2. FC Test Case**

In Figure 3a, we plot the dimensionless temperature flux profile computed with *k*–*ε*, *k*et, and the LES at the end of the 7 days of simulation for the FC test case. The *k*et flux is further decomposed into its downgradient (–*K*<sub>*t*</sub>*∂*<sub>*z*</sub>*Θ*) and its non-local (Γ<sub>*k*et</sub>) components (see Equation 18). The non-local term of KPP is also plotted. The *k*–*ε* profile
exhibits the classic pattern expected for an FC simulation: a linear decrease from the surface to the bottom of the mixed layer where it reaches a minimum which is approximately -0.2 times the surface flux. Note that even if the expression of the total flux changed drastically between k-ε and kεt, the kεt profile is still close to the k-ε profile. The non-local flux is positive (by definition), it is zero at the surface and at the bottom of the mixed layer. Hence, it does not add or remove any heat but rather redistributes heat within the mixed layer. This term is responsible for warming the upper part of the mixed layer and cooling the lower part of the mixed layer (the temperature equation is of the form $D_t\Theta = \dots - \partial_z w'\theta'$ and it is then the sign of $-\partial_z \Gamma_{ket}$ that is important to distinguish between cooling and warming). This is qualitatively the effect we expect from a coherent plume: plumes grow by entraining cold water near the surface, resulting in a warming of the upper part of the mixed layer, and then detrain in the environment, causing a cooling of the bottom part of the mixed layer.

In Figure 3b, we plot the dimensionless temperature profile of k-ε, ket, KPP, and the LES at the end of the 7 days of simulation. Dashed lines highlight the location $h_t$, the depth at which $\partial_z\Theta = 0$ in each case. The overall comparison with the LES is better with ket than with k-ε: while the k-ε model predicts $h_t = h_f = -0.86h_m$ (co-location $h_t = h_f$ results from the definition of a pure downgradient flux, see Equation 11 and Figure 3a), this co-location constraint is relaxed in the ket simulation, for which $h_t = -0.48h_m$, and is thus closer to the LES ($h_t = -0.49h_m$). The KPP scheme predicts $h_t = -0.19h_m$, whereas the KPP simulation without the non-local term $\Gamma_{KPP}$ gives $h_t = -0.83h_m$. Therefore, $\Gamma_{KPP}$ has the same expected effect of raising $h_t$ as the non-local term of kεt, even if none of the two KPP simulations (with or without $\Gamma_{KPP}$) give a satisfactory $h_t$ in comparison to the LES.

Figure 3c shows the temporal evolution of $h_t/h_m$ for the 7 days of simulation. In the LES, $h_t/h_m$ remains between -0.35 and -0.55 during the entire simulation. The k-ε value stabilizes around -0.85, the KPP simulation quickly stabilizes near -0.2 whereas KPP without the non-local term gives a continuous decrease of $h_t/h_m$ with values reaching -0.83 at the end of the 7 days. Finally, the ket curve converges to a value of -0.5 (comparable to the value of the LES). We conclude that the ket model performs better than k-ε and KPP for reproducing the shape of the temperature profile.

In Figure 3d, we plot the dimensionless temperature variance profile of k-ε, ket, and the LES at the end of the 7 days of simulation. Both schemes reproduce the global shape of the LES profile which consists of a local maximum at the surface and a local maximum at $h_m$. However, none of the two schemes give a good order of magnitude for these maxima (values are ≈ 100 times too small). We do not observe a drastic difference between the ket and k-ε profiles. However, if the ket temperature variance profile were equal to the k-ε one, that is, the equilibrium value, ket would reduce to k-ε by definition. In fact, the main differences between the two profiles are located in the inner part of the mixed layer (Figure 3e). In this zone, ket better follows the LES than k-ε. The difference between ket and the equilibrium value obtained by k-ε is due to either the variance memory $\partial_t\overline{\theta'^2}$ or the diffusion $D_{\overline{\theta'^2}}$ terms (cf., Equation 22). We performed two alternative runs of the ket model. One was realized by setting $\partial_t\overline{\theta'^2} = 0$ ("only diffusion") and another by taking $D_{\overline{\theta'^2}} = 0$ ("only memory"). The temperature variance profile of the "only diffusion" ket superimposes with the profile of the "normal" ket (Figures 3d and 3e). It would be tempting to say that the memory has no impact. This would imply that the profiles of the "only memory" ket would superimpose with the k-ε ones. However, this is not the case (Figures 3d and 3e). Therefore, the contributions of $\partial_t\overline{\theta'^2}$ and $D_{\overline{\theta'^2}}$ do not add up linearly: diffusion is the main contribution to the non-equilibrium part of $\overline{\theta'^2}$ in ket, with a moderate impact of the memory contribution.

Figure 3a shows a comparison between the non-local term of ket and the non-local term of KPP at the end of the 7 days of simulation. Qualitatively, these profiles share the property of vanishing at the surface and at the bottom of the mixed layer, and therefore act as a redistribution of heat in the mixed layer. Both are also positive, that is, corresponds to a temperature flux that follows the buoyancy effect (similar to coherent structures that are subjected to the buoyancy force). One difference is that $\Gamma_{KPP}$ has a single-mode shape whereas $\Gamma_{ket}$ presents a bi-modal shape. Analyses of the contribution of the different factors of $\Gamma_{ket} = \int_{h_m}^{h} \frac{\kappa}{\varepsilon} \beta_3 \overline{\theta'^2}$ (not shown) indicated that the mode close to the mixed layer bottom is mainly due to a maximum of $\overline{\theta'^2}$ whereas the mode closest to the surface is a result of a complex interaction of all the terms in the expression of the non-local term. Knowing that $\Gamma_{ket}$ presents a bi-modal shape could be of interest for adapting the KPP non-local term. For example, it would be
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results for the free convection test case. Dimensionless (a) temperature flux profile, (b) temperature profile, (d) temperature variance profile, and (e) zoomed-in temperature variance profile of $k-\epsilon$, $k\epsilon$, KPP, and the Large Eddy Simulations at the end of the 7 days of simulation. The KPP model was run with and without its non-local term $\Gamma_{KPP}$. In (b), dashed lines highlight the location $h_t$ of the zero of the gradient $\partial_z\Theta$, and (c) shows its temporal evolution $h_t/h_m$ for each simulation. In (d), the dashed box highlights the zoomed-in depths presented in (e). The two alternative versions of $k\epsilon$ presented in (d) and (e), “only diffusion” and “only memory,” are described in the text.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Results for the W test case. (a) Temperature profile and (b) horizontal velocity profile of $k-\varepsilon$, $k\varepsilon$ and the Large Eddy Simulations at the end of the 30 hr of simulation. For both cases, the $k-\varepsilon$ and $k\varepsilon$ curves are difficult to distinguish because they are almost superimposed.</figcaption></figure>

possible to write $\Gamma_{KPP}$ as a sum of two polynomials rather than one (Large et al., 1994). Nevertheless, one should be careful when doing a comparison of the two approaches because there are several possible definitions of $\Gamma_{k\varepsilon}$ as explained in Appendix A.

### 3.3. W Test Case

Figure 4 presents the temperature and horizontal velocity profiles of $k-\varepsilon$, $k\varepsilon$ and the LES at the end of the 30 hr of simulation for the W test case. In both cases, the $k-\varepsilon$ and $k\varepsilon$ are almost superimposed. Therefore, in a situation with no surface buoyancy flux, $k\varepsilon$ does not give any difference in comparison to $k-\varepsilon$. This is further demonstrated in Appendix B where it is shown that the part of $\Gamma_{k\varepsilon}$ associated with the non-equilibrium value of $(\partial_z \Theta)^2$ is zero. Since the scope of this paper is to describe the $k\varepsilon$ model and its differences with $k-\varepsilon$, no further analyses are thus made on the differences between $k-\varepsilon$ and the LES in the W test case.

### 3.4. DC Test Case

In order to get more insight into the memory effect mentioned earlier, we now focus on a DC case for which the forcing is no longer stationary. In DC, the $k\varepsilon$ model performs better than $k-\varepsilon$ for reproducing the LES temperature profiles (Figures 5c–5f). At the end of the night (6 hr, Figure 5c), so after 12 hr of FC, the $k\varepsilon$ scheme is better to predict the depth $h_t$ (depth of $\partial_z \Theta = 0$) separating the unstable and the stable part of the profile. At noon (Figure 5d), $k\varepsilon$ is better close to the surface where $k-\varepsilon$ is too unstable and too cold. The same behavior is observed just after sunset (18 hr, Figure 5e). At midnight (Figure 5f), the 10 first meters of the LES profile are nearly neutral. This is mainly reproduced by $k\varepsilon$ whereas the upper 10 m are unstable for $k-\varepsilon$.

The $k\varepsilon$ model better captures the stable/unstable fraction of the temperature profile. This is highlighted in Figure 6a which shows the percentage of the upper 30 m of the water column which is stably stratified ($\partial_z \Theta > 0$) for the 10 days of simulation. During daytime, the solar radiation increases this percentage up to 99% for the LES. This behavior is better reproduced by $k\varepsilon$ (up to 97%) than $k-\varepsilon$ (up to 95%). During nighttime, the free convective conditions decrease this percentage. The lowest value, obtained at the end of the night, decreases over the simulation as the mixed layer depth becomes closer to 30 m. After 10 days, this value is close to 50% for both $k\varepsilon$ and the LES whereas $k-\varepsilon$ gives a too unstable profile with only 20% of it which is stably stratified. Finally, $k\varepsilon$ better captures the instants of the abrupt decreases or increases between daytime and nighttime: $k-\varepsilon$ tends to obtain the abrupt decreases approximately 1 h too early and the sharp increase approximately 2 h too late (Figure 6b).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Results for the diurnal cycle test case over the ninth day of simulation. Subfigures (a) and (b) give an overview over the day by showing (a) eight temperature profiles in a single plot for *ket* and (b) the evolution of the net surface heat flux *Q*₀ (cf., Equation 24). The next subfigures present the comparison between *k*–*ε*, *ket* and the Large Eddy Simulations for (first column) the temperature profiles, and (second column) the temperature variance profiles. These comparisons are given at different times of the day (6, 12, 18, and 24 hr).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Temporal evolution of the percentage of the upper 30 m of the water column which is stably stratified ($\partial_z\Theta > 0$) in the diurnal cycle case. The evolution is shown (a) over the 10 days of simulation, and (b) zoomed in on the last 2 days.</figcaption></figure>

As for the FC case, *ket* does not reproduce better the maxima of the LES temperature variance profile than *k-ε* (Figures 5g–5j). However, the remaining part of the LES profile is better represented by *ket* than *k-ε*. As in the FC case, the “only diffusion” *ket* model (i.e., the *ket* model with variance memory set to zero, cf., Section 3.2) superimposes almost perfectly with *ket*. Therefore, the diffusion is again the main contributor to the non-equilibrium part of $\overline{\theta'^2}$.

**4. Conclusion**

The primary motivation behind this research was the need to improve the representation of oceanic convection processes in ocean models. A key achievement of this study is the analytical derivation of the non-local term within the *k-ε* model. With this strategy, we ensure that all components of the parameterization are consistent: this approach avoids potential redundancy that may arise when adding external components to existing parameterizations. We found that the non-local term is proportional to the temperature variance (high values of the variance correspond to a strongly inhomogeneous temperature field). This heterogeneity in the temperature field may in fact be related to the presence of intense coherent structures (convective plumes) that impact the turbulent fluxes. Compared to the standard *k-ε* for which the temperature variance equation is an equilibrium between the production and the dissipation terms, the temperature variance equation in *ket* contains a diffusion term and a tendency term ($\partial_t \overline{\theta'^2}$, that we also called “memory term” in this study).

To analyze the impact of both of these terms, and in order to assess the performance of *ket*, we have compared this new parameterization to *k-ε* and several LES in three test cases: free convection (FC), wind-driven mixing (W), and diurnal cycle (DC). For W, *ket* is similar to *k-ε*. This comparison suggests that shear turbulence does not trigger non-local processes. It strengthens the idea that non-local processes are made of convective plumes which are initiated by buoyancy fluxes. Indeed, *ket* differs from *k-ε* in the presence of a surface buoyancy flux (FC and DC). In these cases, the main effect of the non-local term is to modify the shape of the temperature profile and particularly the part of the mixed layer that is stably stratified. In FC, the LES temperature profile is stably stratified over the upper 40%–50% of the mixed layer, this is better reproduced by *ket* (50%) than *k-ε* (20%). For DC, *ket* also better reproduces the LES profiles, which present a cycle between very stratified profiles during daytime (almost 100% of the mixed layer is stably stratified) and profiles that resemble free convective ones during nighttime. In both FC and DC, we proved that diffusion is the main actor that drives the temperature variance out of equilibrium. More precisely, the effect of diffusion is to spread the temperature variance in the inner part of the mixed layer and results in a closer agreement to the LES profile.
The next step is to extend the derivation to include salinity. We are currently working on this approach: the main challenge is that the non-local term for salinity involves coupled equations with temperature, making the analytical derivation significantly more complex. Solving these coupled equations analytically is mathematically challenging and may require additional hypotheses. Another issue is also that the computational demands of implementing a coupled temperature-salinity non-local term within ocean models may increase. This can affect model efficiency and require adjustments in computational resources. Despite these difficulties, the extension of the non-local term derivation to salinity promises a more comprehensive and accurate representation of oceanic convection.

In the near future, our research plans also entail a systematic re-evaluation of all k-ε and kεt parameters. To achieve this, we will follow the framework proposed by Wagner et al. (2023) which employs an ensemble of LES in all oceanic convective regime (Legay et al., 2024; Moeng & Sullivan, 1994), in conjunction with Bayesian methods. Using this method, we will identify the key parameters in the closure and, more importantly, obtain an estimate and propose an error bar for all these parameters.

## Appendix A: Part of Γ<sub>ket</sub> Implicitly Present in the k-ε Model

For the sake of clarity, we recall here two expressions of the kεt model that will be discussed in this appendix

$$ \overline{w'\theta'} = -f_h \frac{k^2}{\epsilon} \partial_z \Theta + f_h^* \frac{k}{\epsilon} \beta_3 \overline{\theta'^2}, \quad (A1) $$

$$ \Gamma_{ket} = f_h^* \frac{k}{\epsilon} \beta_3 \overline{\theta'^2}. \quad (A2) $$

In Equation A2, part of Γ<sub>ket</sub> is associated with the equilibrium value of θ'². Indeed, in a k-ε model, we have (cf., Equation F1)

$$ \overline{\theta'^2} = -c_T \frac{k}{\epsilon} \overline{w'\theta'} \partial_z \Theta, \quad (A3) $$

which is the steady-state solution of the equilibrium between the production and the dissipation terms. By construction, if we use this equilibrium in the kεt model, then we obtain the k-ε model. If we insert Equation A3 in Equation A1, we indeed get a temperature flux w'θ' which is entirely downgradient and, by identification, we obtain a relationship between the stability function of k-ε and the stability functions of kεt which is the following

$$ \{f_h\}_{k-\epsilon} = \frac{\{f_h\}_{ket}}{1 + c_T \alpha_N \{f_h^*\}_{ket}}. \quad (A4) $$

This is interesting because it highlights the fact that, within the non-local term, the part of θ'² associated with its equilibrium value results in a downgradient contribution. In the kεt model, we obtain (cf., Equation 22)

$$ \overline{\theta'^2} = -c_T \frac{k}{\epsilon} \overline{w'\theta'} \partial_z \Theta + \frac{c_T}{2} \frac{k}{\epsilon} (D_{\overline{\theta'^2}} - \partial_t \overline{\theta'^2}). \quad (A5) $$

Inserting Equation A5 into Equation A1 gives

$$ \overline{w'\theta'} = -\{f_h\}_{k-\epsilon} \frac{k^2}{\epsilon} \partial_z \Theta + \frac{\{f_h^*\}_{ket}}{1 + c_T \alpha_N \{f_h^*\}_{ket}} \frac{c_T}{2} \frac{k^2}{\epsilon^2} \beta_3 (D_{\overline{\theta'^2}} - \partial_t \overline{\theta'^2}), \quad (A6) $$

and
$$ \hat{\Gamma}_{ket} = \frac{\{f_h^*\}_{ket}}{1 + c_T \alpha_N \{f_h^*\}_{ket}} \frac{c_T k^2}{2 \epsilon^2} \beta_3 (D_{\overline{\theta^2}} - \partial_t \overline{\theta^2}) \quad (A7) $$

is hence the effective non-local contribution of ket in comparison to the k-ε model. It is possible to implement the ket model with Equation A6 rather than Equation A1. However, such an implementation has proved to need very low time steps to be stable (in the order of just a few seconds). This is mainly because the term $(D_{\overline{\theta^2}} - \partial_t \overline{\theta^2})$ cannot be calculated explicitly and must be evaluated after solving the $\overline{\theta^2}$ equation. Quantities coming from different time steps are then mixed and identifying a clear path from time step n to time step n + 1 is not easy.

Keeping definition (Equation A2) is also an intentional choice because this expression exhibits some properties that are lost in Equation A7. Particularly, Equation A2 gives a $\Gamma_{ket}$ that is positive, that is, a turbulent temperature flux that follows the buoyancy effect (cold going down and hot going up) and hence close to what can be expected from a coherent plume that is subjected to the buoyancy force. Also, $\Gamma_{ket}$ is zero at the surface and at the base of the mixed layer, it therefore represents a redistribution of heat in the mixed layer and this again acts as a coherent plume. Thus, in practice, we keep an implementation based on Equation A1 which is described in Section 2.4 and the formula (Equation A7) will be used for comparing $k\epsilon$ to $k-\epsilon$. Particularly, when $\hat{\Gamma}_{ket}$ is close to zero, the $k\epsilon$ model has no sensible impact in comparison to the $k-\epsilon$ model. Profiles of $\hat{\Gamma}_{ket}$ for the FC and the W test cases are presented in Appendix B.

Equations A2 and A7 are two different splits downgradient/non-local of the turbulent temperature flux. These splits are arbitrary since both $\Gamma_{ket}$ and $\hat{\Gamma}_{ket}$ still contain dependencies on $\alpha_n$, which is directly related to $\partial_z \Theta$. We conclude that there is no unique way of defining the non-local term.

## Appendix B: Diagnostics of $\hat{\Gamma}_{ket}$ for FC and W

Based on Appendix A, we know that the part of $\Gamma_{ket}$ which is associated with the non-equilibrium value of $\overline{\theta^2}$ is given by

$$ \hat{\Gamma}_{ket} = \frac{\{f_h^*\}_{ket}}{1 + c_T \alpha_N \{f_h^*\}_{ket}} \frac{c_T k^2}{2 \epsilon^2} \beta_3 (D_{\overline{\theta^2}} - \partial_t \overline{\theta^2}), \quad (B1) $$

and the associated downgradient contribution to $w'\overline{\theta'}$ is

$$ \hat{DG}_{ket} = -\{f_h\}_{k-\epsilon} \frac{k^2}{\epsilon} \partial_z \Theta. \quad (B2) $$

In Figure B1, we evaluate $\hat{\Gamma}_{ket}$ and $\hat{DG}_{ket}$ for both the FC and the W test cases at the end of the simulations, and we compare them to $\Gamma_{ket}$ and $DG_{ket} = -\{f_h\}_{ket} \frac{k^2}{\epsilon} \partial_z \Theta$. For FC, an important proportion of $\Gamma_{ket}$ is explained by $\hat{\Gamma}_{ket}$, and so by the non-equilibrium values of $\overline{\theta^2}$. This is particularly the case in the inner part of the mixed layer ($-0.7 \ge z/h_m \ge -0.3$). This confirms what is observed in Figure 3e, that is, it is mostly in the inner part of the mixed layer that the temperature variance of $k\epsilon$ differs from its equilibrium value in $k-\epsilon$. For W, $\hat{\Gamma}_{ket}$ is close to zero and $k\epsilon$ thus brings no difference in comparison to a $k-\epsilon$ model (see again Figure 4). The temperature variance is close to its equilibrium value and the turbulent temperature flux $w'\overline{\theta'}$ is entirely explained by the downgradient term $\hat{DG}_{ket}$.
<figure><img src="image_9.png" /><figcaption>**Figure B1.** Decomposition of the turbulent temperature flux into $\overline{w'\theta'} = DG_{ket} + \Gamma_{ket}$ and $\overline{w'\theta'} = \overline{DG}_{kt} + \overline{\Gamma}_{kt}$ with $\overline{\Gamma}_{kt}$ the part of $\Gamma_{ket}$ which is associated with the non-equilibrium value of $\overline{\theta^2}$. (a) For the free convection test case at the end of the 7 days of simulation. (b) For the W test case at the end of the 30 hr of simulation.</figcaption></figure>

## Appendix C: Results of the $k-\epsilon-\overline{\theta^2}-\overline{w^2}$ Model

We detail here the results of the $k-\epsilon-\overline{\theta^2}-\overline{w^2}$ model which is a possible extension of the $k\epsilon$ model where the non-equilibrium is also considered for the $\overline{w^2}$ equation. By doing that, Equations 6 and 7 now form a system of eight equations with eight unknowns: $(\overline{u^2}, \overline{v^2}, \overline{u'v'}, \overline{u'w'}, \overline{v'w'}, \overline{u'\theta'}, \overline{v'\theta'}, \overline{w'\theta'})$. We solved this system with Mathematica and we obtained the following expressions:

$$ \overline{u'w'} = -f_m \frac{k^2}{\epsilon} \partial_z U, \qquad (C1) $$

$$ \overline{v'w'} = -f_m \frac{k^2}{\epsilon} \partial_z V, \qquad (C2) $$

$$ \overline{w'\theta'} = -f_h \frac{k^2}{\epsilon} \partial_z \Theta + f_h^* \frac{k}{\epsilon} \beta_3 \overline{\theta^2}, \qquad (C3) $$

which have the same shape as the ones found for the $k\epsilon$ model. Particularly, even if $\overline{w^2}$ is not in equilibrium anymore, the velocity fluxes $\overline{u'w'}$ and $\overline{v'w'}$ are still fully downgradient. The expressions of the stability functions $f_m$, $f_h$ and $f_h^*$ are:

$$ f_m = \frac{n_0 + n_2\alpha_M + n_3\alpha_T + n_4\alpha_W + n_5\alpha_W\alpha_N + n_6\alpha_W\alpha_M}{d_0 + d_1\alpha_N + d_2\alpha_M + d_3\alpha_N\alpha_M + d_5\alpha_M^2}, \qquad (C4) $$

$$ f_h = \frac{n_{2T}\alpha_M + n_{4T}\alpha_W + n_{5T}\alpha_W\alpha_N + n_{6T}\alpha_W\alpha_M}{d_0 + d_1\alpha_N + d_2\alpha_M + d_3\alpha_N\alpha_M + d_5\alpha_M^2}, \qquad (C5) $$

$$ f_h^* = \frac{n_{0T}^* + n_{1T}^*\alpha_N + n_{2T}^*\alpha_M}{d_0 + d_1\alpha_N + d_2\alpha_M + d_3\alpha_N\alpha_M + d_5\alpha_M^2}, \qquad (C6) $$

with $\alpha_N = \frac{k^2}{\epsilon}N^2$, $\alpha_M = \frac{k^2}{\epsilon}M^2$, $\alpha_T = \frac{k}{\epsilon}\beta_3^2\overline{\theta^2}$, and $\alpha_W = \frac{1}{k}\overline{w^2}$. Coefficients $n_i$, $n_{iT}$ and $d_i$ depends on the coefficients $c_i$ and $c_{iT}$; the expressions are given hereafter. Taking the values of the $c_i$ and $c_{iT}$ given in Table 1, the stability functions are approximately as follows
$$f_m = \frac{0.04693 - 0.00005303\alpha_M + 0.001996\alpha_T + 0.0896\alpha_W - 0.002994\alpha_W\alpha_N - 0.0001012\alpha_W\alpha_M}{1 + 0.03361\alpha_N + 0.01342\alpha_M + 0.00006267\alpha_N\alpha_M - 0.00001644\alpha_M^2}, \quad (C7)$$

$$f_h = \frac{0.0002651\alpha_M + 0.1681\alpha_W + 0.005649\alpha_W\alpha_N + 0.002952\alpha_W\alpha_M}{1 + 0.03361\alpha_N + 0.01342\alpha_M + 0.00006267\alpha_N\alpha_M - 0.00001644\alpha_M^2}, \quad (C8)$$

$$f_h^* = \frac{0.1120 + 0.003766\alpha_N + 0.001631\alpha_M}{1 + 0.03361\alpha_N + 0.01342\alpha_M + 0.00006267\alpha_N\alpha_M - 0.00001644\alpha_M^2}. \quad (C9)$$

Here are the expressions of the coefficients $n_i$, $n_{iT}$, and $d_i$:

$$n_0 = \frac{3c_4 - 2c_5}{6c_1}, \quad n_2 = \frac{-c_{4T}(3c_4 - 2c_5)(2 - 2c_{2T} - c_{4T})}{24c_1c_{1T}^2},$$

$$n_3 = \frac{(1 - c_3)(1 - c_{3T})(2c_{1T}c_5 + 3c_1(2 - 2c_{2T} - c_{4T}))}{6c_1^2c_{1T}^2}, \quad n_4 = \frac{2 - 2c_2 + c_5}{2c_1},$$

$$n_5 = \frac{-(1 - c_3)(2c_{1T}c_5 + 3c_1(2 - 2c_{2T} - c_{4T}))}{6c_1^2c_{1T}^2}, \quad n_6 = \frac{-c_{4T}(2 - 2c_2 + c_5)(2 - 2c_{2T} - c_{4T})}{8c_1c_{1T}^2},$$

$$n_{2T} = \frac{c_{4T}(3c_4 - 2c_5)}{12c_1c_{1T}^2}, \quad n_{4T} = \frac{1}{c_{1T}}, \quad n_{5T} = \frac{1 - c_3}{c_1c_{1T}^2},$$

$$n_{6T} = \frac{3c_1c_{4T}(2 - 2c_2 + c_5) + 2c_5c_{1T}(4 - 4c_2 + 3c_5)}{12c_1^2c_{1T}^2},$$

$$n_{0T}^* = \frac{1 - c_{3T}}{c_{1T}}, \quad n_{1T}^* = \frac{(1 - c_3)(1 - c_{3T})}{c_1c_{1T}^2}, \quad n_{2T}^* = \frac{c_5(1 - c_{3T})(4 - 4c_2 + 3c_5)}{6c_1^2c_{1T}},$$

$$d_0 = 1, \quad d_1 = \frac{1 - c_3}{c_1c_{1T}}, \quad d_2 = \frac{c_5(4 - 4c_2 + 3c_5)}{6c_1^2} - \frac{c_{4T}(2 - 2c_{2T} - c_{4T})}{4c_{1T}^2},$$

$$d_3 = \frac{c_5c_{4T}(1 - c_3)}{6c_1^2c_{1T}^2}, \quad d_5 = \frac{-c_5c_{4T}(4 - 4c_2 + 3c_5)(2 - 2c_{2T} - c_{4T})}{24c_1^2c_{1T}^2}. \quad (C10)$$

## Appendix D: Coefficients in the Second-Order Moment Equations

Coefficients $c_1, c_2, c_3, c_4, c_5, c_{1T}, c_{2T}, c_{3T}, c_{4T}, c_T$ used in Equations 3–5 are linked to the coefficients introduced by Canuto et al. (2001) through the following formulas:

$$c_1 = 1/\lambda, \quad c_2 = \alpha_1, \quad c_3 = 1 - \beta_5, \quad c_4 = 4/3 \alpha_1 - 4/5, \quad c_5 = \alpha_1 - \alpha_2, \quad (D1)$$
$$c_{1T} = \lambda_5/2, \quad c_{2T} = 3/4 \alpha_3, \quad c_{3T} = \gamma_1, \quad c_{4T} = \alpha_3/2, \quad c_T = 2 \lambda_8/(1 - \gamma_1).$$

## Appendix E: Expressions of the Main Tensors Under the Boundary Layer Approximation

After applying the boundary layer approximation, the tensors $P_{ij}, B_{ij}, S_{ij}, V_{ij}, Z_{ij}$ used in Equations 3–5 simplify to

$$P_{ij} = \begin{pmatrix} -2 \partial_z U \overline{u'w'} & -\partial_z U \overline{v'w'} - \partial_z V \overline{u'w'} - \partial_z U \overline{w'^2} \\ -\partial_z U \overline{v'w'} - \partial_z V \overline{u'w'} & -2 \partial_z V \overline{v'w'} & -\partial_z V \overline{w'^2} \\ -\partial_z U \overline{w'^2} & -\partial_z V \overline{w'^2} & 0 \end{pmatrix}, \quad (E1)$$
$$B_{ij} = \begin{pmatrix} 0 & 0 & \beta_3 \overline{u'\theta'} \\ 0 & 0 & \beta_3 \overline{v'\theta'} \\ \beta_3 \overline{u'\theta'} & \beta_3 \overline{v'\theta'} & 2 \beta_3 \overline{w'\theta'} \end{pmatrix}, \qquad (E2)$$

$$S_{ij} = \frac{1}{2} \begin{pmatrix} 0 & 0 & \partial_z U \\ 0 & 0 & \partial_z V \\ \partial_z U & \partial_z V & 0 \end{pmatrix}, \qquad (E3)$$

$$V_{ij} = \frac{1}{2} \begin{pmatrix} 0 & 0 & \partial_z U \\ 0 & 0 & \partial_z V \\ -\partial_z U - \partial_z V & 0 & 0 \end{pmatrix}, \qquad (E4)$$

$$Z_{ij} = \begin{pmatrix} \overline{u'w'} \partial_z U & \frac{1}{2} \overline{v'w'} \partial_z U + \frac{1}{2} \overline{u'w'} \partial_z V & \frac{1}{2} \partial_z U (\overline{w'^2} - \overline{u'^2}) - \frac{1}{2} \partial_z V \overline{u'v'} \\ \frac{1}{2} \overline{v'w'} \partial_z U + \frac{1}{2} \overline{u'w'} \partial_z V & \overline{v'w'} \partial_z V & \frac{1}{2} \partial_z V (\overline{w'^2} - \overline{v'^2}) - \frac{1}{2} \partial_z U \overline{u'v'} \\ \frac{1}{2} \partial_z U (\overline{w'^2} - \overline{u'^2}) - \frac{1}{2} \partial_z V \overline{u'v'} & \frac{1}{2} \partial_z V (\overline{w'^2} - \overline{v'^2}) - \frac{1}{2} \partial_z U \overline{u'v'} & -\overline{u'w'} \partial_z U - \overline{v'w'} \partial_z V \end{pmatrix}. \qquad (E5)$$

## Appendix F: The Algebraic System of 10 Equations of the GLS Formalism

For clarity, we give here the explicit writing of the 10 equations presented in Equations 6–8 and that are the basis of the GLS formalism:

$$\begin{aligned} 0 &= -c_1 \frac{\epsilon}{k} (\overline{u'^2} - \frac{2}{3}k) + (1-c_2) \left( -\frac{4}{3} \overline{u'w'} \partial_z U + \frac{2}{3} \overline{v'w'} \partial_z V \right) - \frac{2}{3} (1-c_3) \beta_3 \overline{w'\theta'} - c_5 \overline{u'w'} \partial_z U \\ 0 &= -c_1 \frac{\epsilon}{k} (\overline{v'^2} - \frac{2}{3}k) + (1-c_2) \left( -\frac{4}{3} \overline{v'w'} \partial_z V + \frac{2}{3} \overline{u'w'} \partial_z U \right) - \frac{2}{3} (1-c_3) \beta_3 \overline{w'\theta'} - c_5 \overline{v'w'} \partial_z V \\ 0 &= -c_1 \frac{\epsilon}{k} (\overline{w'^2} - \frac{2}{3}k) + \left( \frac{2}{3} - \frac{2}{3}c_2 + c_5 \right) \left( \overline{u'w'} \partial_z U + \overline{v'w'} \partial_z V \right) + \frac{4}{3} (1-c_3) \beta_3 \overline{w'\theta'} \\ 0 &= -c_1 \frac{\epsilon}{k} \overline{u'v'} - (1-c_2) \left( \overline{v'w'} \partial_z U + \overline{u'w'} \partial_z V \right) - \frac{1}{2} c_5 \left( \overline{v'w'} \partial_z U + \overline{u'w'} \partial_z V \right) \\ 0 &= -c_1 \frac{\epsilon}{k} \overline{u'w'} - (1-c_2) \overline{w'^2} \partial_z U + (1-c_3) \beta_3 \overline{u'\theta'} - \frac{1}{2} c_4 k \partial_z U - \frac{1}{2} c_5 \left( \overline{w'^2} \partial_z U - \overline{u'^2} \partial_z U - \overline{u'v'} \partial_z V \right) \\ 0 &= -c_1 \frac{\epsilon}{k} \overline{v'w'} - (1-c_2) \overline{w'^2} \partial_z V + (1-c_3) \beta_3 \overline{v'\theta'} - \frac{1}{2} c_4 k \partial_z V - \frac{1}{2} c_5 \left( \overline{w'^2} \partial_z V - \overline{v'^2} \partial_z V - \overline{u'v'} \partial_z U \right) \\ 0 &= -c_1 T \frac{\epsilon}{k} \overline{u'\theta'} - \left( 1 - c_{2T} - \frac{1}{2} c_{4T} \right) \overline{w'\theta'} \partial_z U - \overline{u'w'} \partial_z \Theta \\ 0 &= -c_1 T \frac{\epsilon}{k} \overline{v'\theta'} - \left( 1 - c_{2T} - \frac{1}{2} c_{4T} \right) \overline{w'\theta'} \partial_z V - \overline{v'w'} \partial_z \Theta \\ 0 &= -c_1 T \frac{\epsilon}{k} \overline{w'\theta'} - \overline{w'^2} \partial_z \Theta + (1-c_{3T}) \beta_3 \overline{\theta'^2} - \frac{1}{2} c_{4T} \left( \overline{u'\theta'} \partial_z U + \overline{v'\theta'} \partial_z V \right) \\ 0 &= -2 \overline{w'\theta'} \partial_z \Theta - \frac{2}{c_T} \frac{\epsilon}{k} \overline{\theta'^2} \end{aligned} \qquad (F1)$$
Appendix G: Coefficients of the Stability Functions for the GLS Formalism

Coefficients n₀, n₁, n₂, n₀T, n₁T, n₂T, d₀, d₁, d₂, d₃, d₄, d₅ of the GLS stability functions (Equations 12 and 13) have
the following definitions:

$$
n_0 = \frac{4 - 4c_2 + 3c_4}{6c_1},
$$

$$
n_1 = \frac{c_1 c_{1T} c_T (1 - c_{3T}) (4 - 4c_2 + 3c_4) - 2c_1 (1 - c_3) (2 - 2c_{2T} - c_{4T}) + 4c_{1T} (1 - c_3) (c_4 - c_5)}{6c_1^2 c_{1T}^2},
$$

$$
n_2 = \frac{-c_{4T}(4 - 4c_2 + 3c_4)(2 - 2c_{2T} - c_{4T})}{24c_1 c_{1T}^2},
$$

$$
n_{0T} = \frac{2}{3c_{1T}}, \quad n_{1T} = \frac{2(1 - c_3)}{3c_1 c_{1T}^2},
$$

$$
n_{2T} = \frac{c_1 c_{4T} (4 - 4c_2 + 3c_4) + 8c_5 c_{1T} (1 - c_2 + c_5) - 2c_4 c_{1T} (2 - 2c_2 + 3c_5)}{12c_1^2 c_{1T}^2},
$$

$$
d_0 = 1, \quad d_1 = \frac{7 - 7c_3 + 3c_1 c_T (1 - c_{3T})}{3c_1 c_{1T}},
$$

$$
d_2 = \frac{3c_5^2 + 6c_5(1-c_2) + 2(1-c_2)^2}{3c_1^2} - \frac{c_{4T}(2-2c_{2T}-c_{4T})}{4c_{1T}^2},
$$

$$
d_3 = \frac{c_5 c_{1T} (1 - c_3) (2 - 2c_2 + c_5)}{3c_1^3 c_{1T}^2} + \frac{c_1 c_{1T} c_T (1 - c_{3T}) (3c_5^2 + 6c_5(1 - c_2) + 2(1 - c_2)^2)}{3c_1^3 c_{1T}^2} + \frac{c_1 (1 - c_3) (3c_{4T} (1 - c_2 + c_5) - (1 - c_{2T}) (2 - 2c_2 + 3c_5))}{3c_1^3 c_{1T}^2},
$$

$$
d_4 = \frac{(1 - c_3)(4 - 4c_3 + 3c_1c_T(1 - c_{3T}))}{3c_1^2 c_{1T}^2},
$$

$$
d_5 = \frac{-c_{4T}(2 - 2c_{2T} - c_{4T})(3c_5^2 + 6c_5(1 - c_2) + 2(1 - c_2)^2)}{12c_1^2 c_{1T}^2}. \quad (G1)
$$

Taking the values of the cᵢ and cᵢT given in Table 1, a numerical estimate of these expressions gives the following
stability functions

$$
f_m = \frac{0.1067 + 0.01732\alpha_N - 0.0001205\alpha_M}{1 + 0.2398\alpha_N + 0.02872\alpha_M + 0.005154\alpha_N\alpha_M + 0.006930\alpha_N^2 - 0.00003372\alpha_M^2}, \quad (G2)
$$

$$
f_h = \frac{0.1120 + 0.003766\alpha_N + 0.0008871\alpha_M}{1 + 0.2398\alpha_N + 0.02872\alpha_M + 0.005154\alpha_N\alpha_M + 0.006930\alpha_N^2 - 0.00003372\alpha_M^2}. \quad (G3)
$$

Appendix H: The Algebraic System of Nine Equations of the kεt Parameterization

For clarity, we give here the explicit writing of the nine equations presented in Equations 6 and 7 and that are the
basis of the kεt parameterization:
$$
\begin{align*}
0 &= -c_1 \frac{\epsilon}{k} \left( \overline{u'^2} - \frac{2}{3}k \right) + (1-c_2) \left( -\frac{4}{3} \overline{u'w'} \partial_z U + \frac{2}{3} \overline{v'w'} \partial_z V \right) - \frac{2}{3} (1-c_3) \beta_3 \overline{w'\theta'} - c_5 \overline{u'w'} \partial_z U \\
0 &= -c_1 \frac{\epsilon}{k} \left( \overline{v'^2} - \frac{2}{3}k \right) + (1-c_2) \left( -\frac{4}{3} \overline{v'w'} \partial_z V + \frac{2}{3} \overline{u'w'} \partial_z U \right) - \frac{2}{3} (1-c_3) \beta_3 \overline{w'\theta'} - c_5 \overline{v'w'} \partial_z V \\
0 &= -c_1 \frac{\epsilon}{k} \left( \overline{w'^2} - \frac{2}{3}k \right) + \left( \frac{2}{3} - \frac{2}{3}c_2 + c_5 \right) \left( \overline{u'w'} \partial_z U + \overline{v'w'} \partial_z V \right) + \frac{4}{3} (1-c_3) \beta_3 \overline{w'\theta'} \\
0 &= -c_1 \frac{\epsilon}{k} \overline{u'v'} - (1-c_2) \left( \overline{v'w'} \partial_z U + \overline{u'w'} \partial_z V \right) - \frac{1}{2} c_5 \left( \overline{v'w'} \partial_z U + \overline{u'w'} \partial_z V \right) \\
0 &= -c_1 \frac{\epsilon}{k} \overline{u'w'} - (1-c_2) \overline{w'^2} \partial_z U + (1-c_3) \beta_3 \overline{u'\theta'} - \frac{1}{2} c_4 k \partial_z U - \frac{1}{2} c_5 \left( \overline{w'^2} \partial_z U - \overline{u'^2} \partial_z U - \overline{u'v'} \partial_z V \right) \\
0 &= -c_1 \frac{\epsilon}{k} \overline{v'w'} - (1-c_2) \overline{w'^2} \partial_z V + (1-c_3) \beta_3 \overline{v'\theta'} - \frac{1}{2} c_4 k \partial_z V - \frac{1}{2} c_5 \left( \overline{w'^2} \partial_z V - \overline{v'^2} \partial_z V - \overline{u'v'} \partial_z U \right) \\
0 &= -c_{1T} \frac{\epsilon}{k} \overline{u'\theta'} - \left( 1 - c_{2T} - \frac{1}{2} c_{4T} \right) \overline{w'\theta'} \partial_z U - \overline{u'w'} \partial_z \Theta \\
0 &= -c_{1T} \frac{\epsilon}{k} \overline{v'\theta'} - \left( 1 - c_{2T} - \frac{1}{2} c_{4T} \right) \overline{w'\theta'} \partial_z V - \overline{v'w'} \partial_z \Theta \\
0 &= -c_{1T} \frac{\epsilon}{k} \overline{w'\theta'} - \overline{w'^2} \partial_z \Theta + (1-c_{3T}) \beta_3 \overline{\theta'^2} - \frac{1}{2} c_{4T} \left( \overline{u'\theta'} \partial_z U + \overline{v'\theta'} \partial_z V \right) \tag{H1}
\end{align*}
$$

Appendix I: Coefficients of the Stability Functions of the kεt Parameterization

We give hereafter the expressions of the coefficients $n_0, n_1, n_2, n_{0T}, n_{1T}, n_{2T}, n_{0T}^*, n_{1T}^*, n_{2T}^*, d_0, d_1, d_2, d_3, d_4, d_5$ of the kεt stability functions (Equations 19–21). We point out that the expressions of the coefficient $n_{1T}$ and all the coefficients not multiplying $\alpha_n$ (i.e., $n_0, n_2, n_{0T}, n_{2T}, d_0, d_1$ and $d_5$) stay unchanged compared to the GLS ones (given in Appendix G).

$$
\begin{align*}
n_0 &= \frac{4 - 4c_2 + 3c_4}{6c_1}, & n_1 &= \frac{(1 - c_3)(2c_{1T}(c_4 - c_5) - c_1(2 - 2c_{2T} - c_{4T}))}{3c_1^2 c_{1T}^2}, \\
n_2 &= \frac{-c_{4T}(4 - 4c_2 + 3c_4)(2 - 2c_{2T} - c_{4T})}{24c_1 c_{1T}^2}, \\
n_3 &= \frac{(1 - c_3)(1 - c_{3T})(2c_{1T}(4 - 4c_2 + 3c_5) + 3c_1(2 - 2c_{2T} - c_{4T}))}{6c_1^2 c_{1T}^2}, \\
n_{0T} &= \frac{2}{3c_{1T}}, & n_{1T} &= \frac{2(1 - c_3)}{3c_1 c_{1T}^2}, \\
n_{2T} &= \frac{c_1 c_{4T}(4 - 4c_2 + 3c_4) + 8c_5 c_{1T}(1 - c_2 + c_5) - 2c_4 c_{1T}(2 - 2c_2 + 3c_5)}{12c_1^2 c_{1T}^2}, \\
n_{0T}^* &= \frac{1 - c_{3T}}{c_{1T}}, & n_{1T}^* &= \frac{(1 - c_3)(1 - c_{3T})}{c_1 c_{1T}^2}, \\
n_{2T}^* &= \frac{(1 - c_{3T})(3c_5^2 + 6c_5(1 - c_2) + 2(1 - c_2)^2)}{3c_1^2 c_{1T}}, \\
d_0 &= 1, & d_1 &= \frac{7(1 - c_3)}{3c_1 c_{1T}}, & d_2 &= \frac{3c_5^2 + 6c_5(1 - c_2) + 2(1 - c_2)^2}{3c_1^2} - \frac{c_{4T}(2 - 2c_{2T} - c_{4T})}{4c_{1T}^2}, \\
d_3 &= \frac{(1 - c_3)(3c_1 c_{4T}(1 - c_2 + c_5) + c_5 c_{1T}(2 - 2c_2 + c_5) - c_1(1 - c_{2T})(2 - 2c_2 + 3c_5))}{3c_1^3 c_{1T}^2},
\end{align*}
$$
$$d_4 = \frac{4(1-c_3)^2}{3c_1^2 c_{1T}^2}, \quad d_5 = \frac{-c_{4T}(2-2c_{2T}-c_{4T})(3c_5^2+6c_5(1-c_2)+2(1-c_2)^2)}{12c_1^2 c_{1T}^2}. \quad (11)$$

Taking the values of the $c_i$ and $c_{iT}$ given in Table 1, a numerical estimate of these expressions gives the following stability functions

$$f_m = \frac{0.1067 + 0.0001072\alpha_N - 0.0001205\alpha_M + 0.004673\alpha_T}{1 + 0.07843\alpha_N + 0.02872\alpha_M + 0.0003389\alpha_N\alpha_M + 0.001506\alpha_N^2 - 0.00003372\alpha_M^2}, \quad (12)$$

$$f_h = \frac{0.1120 + 0.003766\alpha_N + 0.0008871\alpha_M}{1 + 0.07843\alpha_N + 0.02872\alpha_M + 0.0003389\alpha_N\alpha_M + 0.001506\alpha_N^2 - 0.00003372\alpha_M^2}, \quad (13)$$

$$f_h^* = \frac{0.1120 + 0.003766\alpha_N + 0.003344\alpha_M}{1 + 0.07843\alpha_N + 0.02872\alpha_M + 0.0003389\alpha_N\alpha_M + 0.001506\alpha_N^2 - 0.00003372\alpha_M^2}. \quad (14)$$

## Data Availability Statement

All the codes used for the study are published on GitHub (Legay, 2024, https://github.com/legaya/James2024-ket/). The repository contains the three Jupyter Notebooks used for performing the 1D simulations and all the analyses, the version of the “oceanmixedlayers” package we used (Reichl et al., 2022, https://github.com/breichl/ocean-mixedlayers), the 1D model described in Section 2.4 as Fortran Modules, the Fortran codes needed for generating these modules, the files needed to perform the LES, and the LES results as netCDF files.

## References

Arakawa, A., & Schubert, W. H. (1974). Interaction of a cumulus cloud ensemble with the large-scale environment, Part I. *Journal of the Atmospheric Sciences*, 31(3), 674–701. https://doi.org/10.1175/1520-0469(1974)031(0674:IOACCE)2.0.CO;2

Axell, L. B. (2002). Wind-driven internal waves and Langmuir circulations in a numerical ocean model of the southern Baltic Sea. *Journal of Geophysical Research*, 107(C11), 25-1-25-20. https://doi.org/10.1029/2001JC000922

Burchard, H. (2002). In S. Bhattacharji, G. M. Friedman, H. J. Neugebauer, & A. Seilacher (Eds.), *Applied turbulence modelling in marine waters* (Vol. 100). Springer Berlin Heidelberg. https://doi.org/10.1007/3-540-45419-5

Burchard, H., & Baumert, H. (1995). On the performance of a mixed-layer model based on the $\kappa-\epsilon$ turbulence closure. *Journal of Geophysical Research*, 100(C5), 8523–8540. https://doi.org/10.1029/94JC03229

Burchard, H., & Bolding, K. (2001). Comparative analysis of four second-moment turbulence closure models for the oceanic mixed layer. *Journal of Physical Oceanography*, 31(8), 1943–1968. https://doi.org/10.1175/1520-0485(2001)031(1943:CAOFSM)2.0.CO;2

Burchard, H., & Petersen, O. (1999). Models of turbulence in the marine environment—A comparative study of two-equation turbulence models. *Journal of Marine Systems*, 21(1-4), 29–53. https://doi.org/10.1016/S0924-7963(99)00004-4

Canuto, V. M., Howard, A., Cheng, Y., & Dubovikov, M. S. (2001). Ocean turbulence. Part I: One-point closure model—Momentum and heat vertical diffusivities. *Journal of Physical Oceanography*, 31(6), 1413–1426. https://doi.org/10.1175/1520-0485(2001)031(1413:OTPIOP)2.0.CO;2

Canuto, V. M., Howard, A., Cheng, Y., & Dubovikov, M. S. (2002). Ocean turbulence. Part II: Vertical diffusivities of momentum, heat, salt, mass, and passive scalars. *Journal of Physical Oceanography*, 32(1), 240–264. https://doi.org/10.1175/1520-0485(2002)032(0240:OTPIVD)2.0.CO;2

Cheng, Y., Canuto, V. M., & Howard, A. M. (2002). An improved model for the turbulent pbl. *Journal of the Atmospheric Sciences*, 59(9), 1550–1565. https://doi.org/10.1175/1520-0469(2002)059(1550:AIMFTT)2.0.CO;2

Cheng, Y., Canuto, V. M., Howard, A. M., Ackerman, A. S., Kelley, M., Fridlind, A. M., et al. (2020). A second-order closure turbulence model: New heat flux equations and no critical Richardson number. *Journal of the Atmospheric Sciences*, 77(8), 2743–2759. https://doi.org/10.1175/JAS-D-19-0240.1

Deardorff, J. W. (1972). Theoretical expression for the countergradient vertical heat flux. *Journal of Geophysical Research*, 77(30), 5900–5904. https://doi.org/10.1029/JC077i030p05900

Deleersnijder, E., & Luyten, P. (1994). On the practical advantages of the quasi-equilibrium version of the Mellor and Yamada level 2.5 turbulence closure applied to marine modelling. *Applied Mathematical Modelling*, 18(5), 281–287. https://doi.org/10.1016/0307-904X(94)90336-0

Emanuel, K. A. (1991). A scheme for representing cumulus convection in large-scale models. *Journal of the Atmospheric Sciences*, 48(21), 2313–2329. https://doi.org/10.1175/1520-0469(1991)048(2313:ASFRCC)2.0.CO;2

Fearon, G., Herbette, S., Veitch, J., Cambon, G., Lucas, A. J., Lemarié, F., & Vichi, M. (2020). Enhanced vertical mixing in coastal upwelling systems driven by diurnal-inertial resonance: Numerical experiments. *Journal of Geophysical Research: Oceans*, 125(9), e2020JC016208. https://doi.org/10.1029/2020JC016208

Fox-Kemper, B., Ferrari, R., & Hallberg, R. (2008). Parameterization of mixed layer eddies. Part I: Theory and diagnosis. *Journal of Physical Oceanography*, 38(6), 1145–1165. https://doi.org/10.1175/2007JPO3792.1

Garanaik, A., Pereira, F. S., Smith, K., Robey, R., Li, Q., Pearson, B., & Van Roekel, L. (2024). A new hybrid mass-flux/high-order turbulence closure for ocean vertical mixing. *Journal of Advances in Modeling Earth Systems*, 16(1), e2023MS003846. https://doi.org/10.1029/2023MS003846

## Acknowledgments

We sincerely thank Gregory LeClaire Wagner, Ye Cheng, and one anonymous reviewer for their valuable time and helpful feedback. Their input has been instrumental in enhancing this manuscript. We thank Florian Lemarié for providing the 1D code. We thank Thierry Penduff and Pierre Brasseur for stimulating discussions on the mixed layer dynamics. The LES were performed on the GENCI HPC (allocation A0150112020). All other computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. This research was partially funded by the French Agence Nationale de la Recherche (ANR), project ANR-23-CE01-0009.
Garcia, J. R., & Mellado, J. P. (2014). The two-layer structure of the entrainment zone in the convective boundary layer. *Journal of the Atmospheric Sciences*, 71(6), 1935–1955. https://doi.org/10.1175/JAS-D-13-0148.1

Gaspar, P., Grégoris, Y., & Lefevre, J.-M. (1990). A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station papa and long-term upper ocean study site. *Journal of Geophysical Research*, 95(C9), 16–16193. https://doi.org/10.1029/JC095iC09p16179

Ghannam, K., Duman, T., Salesky, S. T., Chamecki, M., & Katul, G. (2017). The non-local character of turbulence asymmetry in the convective atmospheric boundary layer. *Quarterly Journal of the Royal Meteorological Society*, 143(702), 494–507. https://doi.org/10.1002/qj.2937

Gibbs, J. A., Fedorovich, E., & Van Eijk, A. M. J. (2011). Evaluating Weather Research and Forecasting (WRF) model predictions of turbulent flow parameters in a dry convective boundary layer. *Journal of Applied Meteorology and Climatology*, 50(12), 2429–2444. https://doi.org/10.1175/2011JAMC2661.1

Giordani, H., Bourdallé-Badie, R., & Madec, G. (2020). An eddy-diffusivity mass-flux parameterization for modeling oceanic convection. *Journal of Advances in Modeling Earth Systems*, 12(9), e02078. https://doi.org/10.1029/2020MS002078

Haghshenas, A., & Mellado, J. P. (2019). Characterization of wind-shear effects on entrainment in a convective boundary layer. *Journal of Fluid Mechanics*, 858, 145–183. https://doi.org/10.1017/jfm.2018.761

Hanjalić, K., & Launder, B. E. (1972). A Reynolds stress model of turbulence and its application to thin shear flows. *Journal of Fluid Mechanics*, 52(4), 609–638. https://doi.org/10.1017/S002211207200268X

Harcourt, R. R. (2015). An improved second-moment closure model of Langmuir turbulence. *Journal of Physical Oceanography*, 45(1), 84–103. https://doi.org/10.1175/JPO-D-14-0046.1

Kato, H., & Phillips, O. M. (1969). On the penetration of a turbulent layer into stratified fluid. *Journal of Fluid Mechanics*, 37(4), 643–655. https://doi.org/10.1017/S0022112069000784

Large, W. G., McWilliams, J. C., & Doney, S. C. (1994). Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. *Reviews of Geophysics*, 32(4), 363–404. https://doi.org/10.1029/94RG01872

Lazier, J. (2001). Deep convection. In *Encyclopedia of Ocean Sciences* (pp. 634–643). Elsevier. https://doi.org/10.1006/rwos.2001.0113

Legay, A. (2024). legaya/James2024-ket: 1.0 [Software]. Zenodo. https://doi.org/10.5281/zenodo.13121130

Legay, A., Deremble, B., Penduff, T., Brasseur, P., & Molines, J. (2024). A framework for assessing ocean mixed layer depth evolution. *Journal of Advances in Modeling Earth Systems*, 16(10), e2023MS004198. https://doi.org/10.1029/2023MS004198

Lemarié, F., Samson, G., Redelsperger, J.-L., Giordani, H., Brivoal, T., & Madec, G. (2021). A simplified atmospheric boundary layer model for an improved representation of air-sea interactions in eddying oceanic models: Implementation and first evaluation in NEMO (4.0). *Geophysical Model Development*, 14(1), 543–572. https://doi.org/10.5194/gmd-14-543-2021

Luyten, J., Pedlosky, J., & Stommel, H. (1983). The ventilated thermocline. *Journal of Physical Oceanography*, 13(2), 292–309. https://doi.org/10.1175/1520-0485(1983)013%3C0292:TVT%3E2.0.CO;2

Marshall, J., & Schott, F. (1999). Open-ocean convection: Observations, theory, and models. *Reviews of Geophysics*, 37(1), 1–64. https://doi.org/10.1029/98RG02739

Mellor, G. L., & Blumberg, A. (2004). Wave breaking and ocean surface layer thermal response. *Journal of Physical Oceanography*, 34(3), 693–698. https://doi.org/10.1175/2517.1

Mellor, G. L., & Yamada, T. (1982). Development of a turbulence closure model for geophysical fluid problems. *Reviews of Geophysics*, 20(4), 851–875. https://doi.org/10.1029/RG020i004p00851

Mironov, D. V., Gryanik, V. M., Moeng, C., Olbers, D. J., & Warncke, T. H. (2000). Vertical turbulence structure and second-moment budgets in convection with rotation: A large-eddy simulation study. *Quarterly Journal of the Royal Meteorological Society*, 126(563), 477–515. https://doi.org/10.1002/qj.49712656306

Moeng, C.-H., & Sullivan, P. P. (1994). A comparison of shear- and buoyancy-driven planetary boundary layer flows. *Journal of the Atmospheric Sciences*, 51(7), 999–1022. https://doi.org/10.1175/1520-0469(1994)051(0999:ACOSAB)2.0.CO;2

Moeng, C.-H., & Wyngaard, J. C. (1989). Evaluation of turbulent transport and dissipation closures in second-order modeling. *Journal of the Atmospheric Sciences*, 46(14), 2311–2330. https://doi.org/10.1175/1520-0469(1989)046(2311:EOTTAD)2.0.CO;2

Patankar, S. (1980). *Numerical heat transfer and fluid flow*. McGraw-Hill. https://doi.org/10.1201/9781482234213

Paulson, C. A., & Simpson, J. J. (1977). Irradiance measurements in the upper ocean. *Journal of Physical Oceanography*, 7(6), 952–956. https://doi.org/10.1175/1520-0485(1977)007(0952:IMITUO)2.0.CO;2

Pollard, R. T., Rhines, P. B., & Thompson, R. O. R. Y. (1973). The deepening of the wind-Mixed layer. *Geophysical Fluid Dynamics*, 4(4), 381–404. https://doi.org/10.1080/03091927208236105

Popinet, S. (2020). Basilisk flow solver and PDE library. Retrieved from http://basilisk.fr

Price, J. F. (1979). On the scaling of stress-driven entrainment experiments. *Journal of Fluid Mechanics*, 90(3), 509–529. https://doi.org/10.1017/S0022112079002366

Reffray, G., Bourdallé-Badie, R., & Calone, C. (2015). Modelling turbulent vertical mixing sensitivity using a 1-D version of NEMO. *Geoscientific Model Development*, 8(1), 69–86. https://doi.org/10.5194/gmd-8-69-2015

Reichl, B. G., Adcroft, A., Griffies, S. M., & Hallberg, R. (2022). A potential energy analysis of ocean surface mixed layers. *Journal of Geophysical Research: Oceans*, 127(7), e2021JC018140. https://doi.org/10.1029/2021JC018140

Reichl, B. G., & Hallberg, R. (2018). A simplified energetics based planetary boundary layer (ePBL) approach for ocean climate simulations. *Ocean Modelling*, 132, 112–129. https://doi.org/10.1016/j.ocemod.2018.10.004

Rodi, W. (1987). Examples of calculation methods for flow and mixing in stratified fluids. *Journal of Geophysical Research*, 92(C5), 5305. https://doi.org/10.1029/JC092iC05p05305

Siebesma, A. P., Soares, P. M. M., & Teixeira, J. (2007). A combined eddy-diffusivity mass-flux approach for the convective boundary layer. *Journal of the Atmospheric Sciences*, 64(4), 1230–1248. https://doi.org/10.1175/JAS3888.1

Souza, A. N., Wagner, G. L., Ramadhan, A., Allen, B., Churavy, V., Schloss, J., et al. (2020). Uncertainty quantification of ocean parameterizations: Application to the K-Profile-Parameterization for penetrative convection. *Journal of Advances in Modeling Earth Systems*, 12(12), e2020MS002108. https://doi.org/10.1029/2020MS002108

Stull, R. B. (1988). *An introduction to boundary layer meteorology* (Vol. 13). Springer Science & Business Media. https://doi.org/10.1007/978-94-009-3027-8

Thangam, S., Abid, R., & Speziale, C. G. (1992). Application of a new K-tau model to near wall turbulent flows. *AIAA Journal*, 30(2), 552–554. https://doi.org/10.2514/3.10952

Treguier, A. M., de Boyer Montégut, C., Bozec, A., Chassignet, E. P., Fox-Kemper, B., Hogg, A. M., et al. (2023). The mixed layer depth in the Ocean Model Intercomparison Project (OMIP): Impact of resolving mesoscale eddies. *EGUsphere*, 2023, 1–43. https://doi.org/10.5194/egusphere-2023-310
Troen, I. B., & Mahrt, L. (1986). A simple model of the atmospheric boundary layer; sensitivity to surface evaporation. *Boundary-Layer Meteorology*, 37(1-2), 129-148. https://doi.org/10.1007/BF00122760

Umlauf, L., & Burchard, H. (2003). A generic length-scale equation for geophysical turbulence models. *Journal of Marine Research*, 61(2), 235-265. https://doi.org/10.1357/002224003322005087

Umlauf, L., & Burchard, H. (2005). Second-order turbulence closure models for geophysical boundary layers. A review of recent work. *Continental Shelf Research*, 25(7-8), 795-827. https://doi.org/10.1016/j.csr.2004.08.004

Umlauf, L., Burchard, H., & Hutter, K. (2003). Extending the k-ω turbulence model towards oceanic applications. *Ocean Modelling*, 5(3), 195-218. https://doi.org/10.1016/S1463-5003(02)00039-2

Van Roekel, L., Adcroft, A. J., Danabasoglu, G., Griffies, S. M., Kauffman, B., Large, W., et al. (2018). The KPP boundary layer scheme for the ocean: Revisiting its formulation and benchmarking one-dimensional simulations relative to LES. *Journal of Advances in Modeling Earth Systems*, 10(11), 2647-2685. https://doi.org/10.1029/2018MS001336

Wagner, G., Hillier, A., Constantinou, N. C., Silvestri, S., Souza, A., Burns, K., et al. (2023). CATKE: A turbulent-kinetic-energy-based parameterization for ocean microturbulence with dynamic convective adjustment. *arXiv e-prints*. arXiv:2306.13204. https://doi.org/10.48550/arXiv.2306.13204

Warner, J. C., Sherwood, C. R., Arango, H. G., & Signell, R. P. (2005). Performance of four turbulence closure models implemented using a generic length scale method. *Ocean Modelling*, 8(1-2), 81-113. https://doi.org/10.1016/j.ocemod.2003.12.003

Wilcox, D. C. (1988). Reassessment of the scale-determining equation for advanced turbulence models. *AIAA Journal*, 26(11), 1299-1310. https://doi.org/10.2514/3.10041

Williams, R. G., Marshall, J. C., & Spall, M. A. (1995). Does Stommel's mixed layer "demon" work? *Journal of Physical Oceanography*, 25(12), 3089-3102. https://doi.org/10.1175/1520-0485(1995)025(3089:DSMLW)2.0.CO;2

Willis, G., & Deardorff, J. (1974). A laboratory model of the unstable planetary boundary layer. *Journal of the Atmospheric Sciences*, 31(5), 1297-1307. https://doi.org/10.1175/1520-0469(1974)031(1297:ALMOTU)2.0.CO;2

Zeierman, S., & Wolfshtein, M. (1986). Turbulent time scale for turbulent-flow calculations. *AIAA Journal*, 24(10), 1606-1610. https://doi.org/10.2514/3.9490

Zhou, B., Sun, S., Yao, K., & Zhu, K. (2018). Reexamining the gradient and countergradient representation of the local and nonlocal heat fluxes in the convective boundary layer. *Journal of the Atmospheric Sciences*, 75(7), 2317-2336. https://doi.org/10.1175/JAS-D-17-0198.1