Evolution of small-scale turbulence at large Richardson numbers

The theory of stratified turbulent flow developed earlier by the authors is applied to data from different areas of the ocean. It is shown that turbulence can be amplified and supported even at large gradient Richardson numbers. The cause of that is the exchange between kinetic and potential energies of turbulence. Using the profiles of Brunt–Väisälä frequency and vertical current shear given in Forryan et al. (2013), the profiles of the kinetic energy dissipation rate are calculated. The results are in reasonable agreement with the experimental data.

1 Introduction

At present, it is well established that the processes in the upper mixed layer of the ocean and inland waters play a significant role in both the development of global climate models and the creation of regional weather forecast models (e.g. Hostetler et al., 1993; Ljungemyr et al., 1996; Tsuang et al., 2001; Mackay, 2006). Small-scale and mesoscale processes effectively interact with each other and provide an energy sink for currents and waves of larger scales. Such processes as wind wave breaking and surface and subsurface shear flows can cause turbulent mixing and the resulting fine-structure formation with areas of sharp gradients of temperature and salinity. Whereas the mechanisms of generation of small-scale turbulence are understood reasonably well, the problem of its interaction with other types of motions and long-time support is less clear. Earlier works were based on Miles' instability condition $R_{\mathrm{i}}<\mathrm{1}/\mathrm{4}$, where R_i is the gradient Richardson number. However, in many observations, turbulence exists in a quasi-stationary regime at much larger R_i, up to 10 and more. In some works, this was explained by the presence of a fine structure of currents, with thin layers of strong shear (Smyth et al., 2013). A more general description is based on the semi-empirical K−ε equations (Burchard, 2002; Burchard and Bolding, 2002; Mellor and Yamada, 1982) showing that the developed turbulence can be amplified and supported under a softer condition R_i<1 (Monin and Yaglom, 1965). Similar equations were used in more specific models of formation of the upper turbulent layer (Ostrovsky and Soustova, 1969) and of the action of internal waves on turbulence (Ivanov et al., 1983; Strang and Fernando, 2001; Stretch et al., 2001). However, even that is insufficient to explain many observations in the ocean and atmosphere, where the turbulence is observed at significantly larger R_i (Forryan et al., 2013; Avicola et al., 2007; Galperin et al., 2021). New theoretical models have also appeared for describing nonstationary turbulent processes in the atmosphere and ocean. They are based on a spectral approach, confirming, in particular, the absence of a critical Richardson number for describing turbulent-wave processes in a stratified fluid (Sukoriansky et al., 2003, 2005b, a; Galperin and Sukoriansky, 2020; Galperin et al., 2021, 2007).

In Ostrovsky and Troitskaya (1987), within the framework of a kinetic approach, a closed nonstationary model of turbulence interacting with a variable current was suggested. It includes mutual transformation between kinetic and potential energies of turbulence. The latter is associated with the density fluctuations occurring in stratified turbulence (Monin and Ozmidov, 1981). The theory suggested in Ostrovsky and Troitskaya (1987) is based on the equation for the variable probability distribution function of fluid velocity and density. This approach reduces the uncertainty of standard semi-empirical K−ε schemes and naturally includes the potential energy of turbulence in the model. As a result, small-scale turbulence can be supported at a non-zero level by the average shear at any finite values of the gradient Richardson number, without a threshold. Later this approach was further developed in Zilitinkevich et al. (2007), Zilitinkevich et al. (2013), and Soustova et al. (2020) in application to atmospheric turbulence, where the energy and flux budget (EFB) model was added to the theory. The correspondence between this theory and the K−ε model as well as the proper parametrization using the turbulent Prandtl number are discussed in recent works (Gladskikh et al., 2023; Soustova et al., 2020).

In this paper, the kinetic model of turbulence is used to describe the evolution and structure of the upper turbulent layer with the parameters taken from in situ observations. Particular attention is paid to the cases of the large Richardson number and the role of turbulent potential energy in explaining the observation data. As an example, we use some data from the paper (Forryan et al., 2013) that provided a relatively detailed set of measurements for three cruises taken in 2006–2009 in different areas of the world ocean: North Atlantic (cruise D3406, June–July 2006, and cruise D321, July–August 2007) and Southern Ocean (cruise JC29, November–December 2008). These experiments aimed to study turbulent mixing in the presence of a stratified shear flow associated with mesoscale motions such as eddies and fronts. With the given profiles of current shear and buoyancy frequency taken from Forryan et al. (2013), the theory developed in Ostrovsky and Troitskaya (1987) yields the results that satisfactorily agree with the measurements of the turbulent dissipation rate given in Forryan et al. (2013). The details of the measurements can be found in Forryan et al. (2013) and references therein.

Figure 1 Profiles of the Richardson number for three cruises calculated from Fig. 2 of Forryan et al. (2013).

Figure 2 Temporal variation of kinetic (a) and potential (b) energy for the conditions of cruise JC29. From top to bottom: $z=\mathrm{20},\mathrm{30},\mathrm{50},\mathrm{100}$, and 180 m. Here, $K_{\mathrm{0}}=P_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{6}}$ m² s⁻².

2 Basic equations

The general equations obtained in Ostrovsky and Troitskaya (1987) (see also Gladskikh et al., 2023), are shown in the Appendix. Here they will be used for the particular case of known profiles of horizontal current shear $\partial \langle u_x\rangle /\partial z=V_z$, where 〈u_x〉=V(z) is the ensemble-average horizontal velocity, and the average density is 〈ρ(z)〉=R(z). Here, z is the vertical coordinate. As a result, we have a system of two equations for the kinetic energy of turbulence K(t,z) and potential energy P(t,z) per unit volume. The latter is related to the density fluctuations $\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle$:

$$\begin{array}{}\text{(1)}& P={\displaystyle \frac{\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle {\mathit{g}}^{\mathrm{2}}}{\mathrm{2}{N}^{\mathrm{2}}{R}^{\mathrm{2}}}}.\end{array}$$

Here, g is the gravity acceleration, R(z) is defined above, and N²(z) is the squared Brunt–Väisala frequency. Note that here the z direction is chosen downwards. Under the above conditions, the general Eq. (18) of Ostrovsky and Troitskaya (1987) and Eq. (6) of Gladskikh et al. (2023) are reduced to two equations for K:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial K}{\partial t}}& =V_z^{\mathrm{2}}L\sqrt{K}-N^{\mathrm{2}}L\sqrt{K}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{3}P}{K}}(\mathrm{1}-G)\right)\\ & -{\displaystyle \frac{G{K}^{\mathrm{3}/\mathrm{2}}}{L}}+{\displaystyle \frac{\mathrm{5}}{\mathrm{3}}}{\displaystyle \frac{\partial }{\partial z}}\left(L\sqrt{b}{\displaystyle \frac{\partial K}{\partial z}}\right),\\ {\displaystyle \frac{\partial P}{\partial t}}& =N^{\mathrm{2}}L\sqrt{K}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{3}P}{K}}(\mathrm{1}-G)\right)\\ & -{\displaystyle \frac{D{K}^{\mathrm{1}/\mathrm{2}}}{L}}P+{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{L\sqrt{K}}{N^{\mathrm{2}}{R}^{\mathrm{2}}}}{\displaystyle \frac{\partial \left(N^{\mathrm{2}}{R}^{\mathrm{2}}P\right)}{\partial z}}\right).\end{array}\end{array}$$

Here, L is the outer scale of turbulence and G is the anisotropy parameter tending to 1 for strongly anisotropic turbulence with a small vertical scale compared to the horizontal scale. For details, see Ostrovsky and Troitskaya (1987). Here we consider the model of locally isotropic turbulence, for which G∼0.5. The parameter L can be taken from empirical data (Rodi, 1980) or found from the turbulence spectrum (Forryan et al., 2013; Lozovatsky et al., 2006). C and D are empirical constants. The terms $C{K}^{\mathrm{3}/\mathrm{2}}/L$ and $D{K}^{\mathrm{1}/\mathrm{2}}P/L$ in Eq. (2) define the dissipation rates of kinetic and potential energy, respectively.

Before analyzing the full system (2), we note that some significant conclusions can be made from a reduced, local ODE system following from Eq. (2) after neglecting the last, diffusive terms in these equations:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial K}{\partial t}}& =V_z^{\mathrm{2}}L\sqrt{K}-N^{\mathrm{2}}L\sqrt{K}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{3}P}{K}}(\mathrm{1}-G)\right)\\ & -{\displaystyle \frac{C{K}^{\mathrm{3}/\mathrm{2}}}{L}},\\ {\displaystyle \frac{\partial P}{\partial t}}& =N^{\mathrm{2}}L\sqrt{K}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{3}P}{K}}(\mathrm{1}-G)\right)-{\displaystyle \frac{D{K}^{\mathrm{1}/\mathrm{2}}}{L}}P.\end{array}\end{array}$$

Figure 3 Profiles of kinetic (blue) and potential (green) energies for JC29 at t=40 000 s.

Figure 4 Profiles of the turbulent kinetic energy dissipation rate for JC29. Green – interpolated data of Forryan et al. (2013). Blue – theory.

The coordinate z is now a parameter in these equations. In particular, they define a stationary distribution (a stable equilibrium point on the phase plane of variables K and P):

$$\begin{array}{}\text{(4)}& K_{\mathrm{st}}\left(z\right)={\displaystyle \frac{V_z^{\mathrm{2}}{L}^{\mathrm{2}}}{\mathrm{2}C}}f\left(R_{\mathrm{i}}\right),P_{\mathrm{st}}\left(z\right)={\displaystyle \frac{V_z^{\mathrm{2}}{L}^{\mathrm{2}}}{D}}-K_{\mathrm{st}},\end{array}$$

where $R_{\mathrm{i}}=V_z^{\mathrm{2}}/N^{\mathrm{2}}$ is the Richardson number, and

$$\begin{array}{}\text{(5)}& \begin{array}{rl}f\left(R_{\mathrm{i}}\right)& =\mathrm{1}-(\mathrm{4}-\mathrm{3}G)R_{\mathrm{i}}+\\ & {\left[\mathrm{1}+R_{\mathrm{i}}^{\mathrm{2}}(\mathrm{4}-\mathrm{3}G{)}^{\mathrm{2}}+R_{\mathrm{i}}(\mathrm{4}-\mathrm{6}G)\right]}^{\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

It is noteworthy that, at R_i→∞, f has a non-zero limit $f_{\mathrm{\infty }}=\mathrm{6}(\mathrm{1}-G)/(\mathrm{4}-\mathrm{3}G)>\mathrm{0}$ (it is 1.2 for G=0.5). Hence, the turbulent energy remains finite at large Richardson numbers. General features of this solution and the turbulent Prandtl number following from it are discussed in Ostrovsky and Troitskaya (1987) and Gladskikh et al. (2023). In what follows, the applicability of these simple solutions will be verified by comparison with the solutions of the full Eq. (2). In the experiments, the dissipation rate of kinetic turbulent energy ε is commonly measured as a characteristic of turbulence. Within the semi-empirical approach, it is defined as (Kolmogorov, 1941)

$$\begin{array}{}\text{(6)}& \mathit{\epsilon }={\displaystyle \frac{C{K}^{\mathrm{3}/\mathrm{2}}}{L}}.\end{array}$$

It is necessary to choose the empirical constants in the above equations and in Eq. (6). There exists a broad literature discussing these values for different laboratory and on-site conditions. In the semi-empirical models they are commonly defined by scaling. We choose the outer turbulence scale based on the results of the spectral approach (Forryan et al., 2013; Galperin et al., 2007), in which the minimal wave number for the energy-carrying spectrum approximated by empirical functions is about 2 cpm. Here we take L=0.58 m. The range of the constant C is also wide in the literature. Since we are mainly interested in the quasi-equilibrium regime when the shear source is balanced with dissipation of turbulent energy, we use the data of Rodi (1980) to take $C=D=\mathrm{0.09}$. Anyway, we are mainly concerned with the order of the obtained values; indeed, in the data of Forryan et al. (2013) considered below, the spread of the data is up to an order of magnitude.

In what follows, we solve systems (2) and (3) using the Wolfram Mathematica 13 program and compare them with each other and the data of in situ measurements.

Figure 5 Profiles of kinetic (blue) and potential (green) energies for D306 at t=40 000 s.

Figure 6 Profiles of the turbulent kinetic energy dissipation rate for D306. Green – interpolated data of Forryan et al. (2013). Blue – theory.

3 Application of the local model

As mentioned, here we apply the theory to the data of the three cruises described in Forryan et al. (2013). First, we digitized the red curves in Fig. 2 of that paper (as mentioned, there is a large dispersion of real data, but we naturally use the mean profiles). Then, using the given profiles of V_z and N², we calculated the Richardson number as shown in Fig. 1.

Note that, for a water density of 1000 kg m⁻³, the static pressure (dbar) coincides with depth in meters, so that in our calculations we use the depth, neglecting small differences in density. As seen from Fig. 1, the Richardson number exceeds 1 within all ranges of the available data, and its maximum lies in the range of 10–23.

Figure 7 Profiles of kinetic (blue) and potential (green) energies at t=40 000 for D321.

Figure 8 Profiles of the turbulent kinetic energy dissipation rate for D321. Green – interpolated data of Forryan et al. (2013). Blue – theory.

3.1 Cruise JC29

Now, using the interpolation of digitized data for N² and V_z given in Fig. 2 of Forryan et al. (2013), we solved Eq. (3) with the initial conditions $K(\mathrm{0},z)=K_{\mathrm{0}}\exp (-\mathrm{0.1}z)$ and $P(\mathrm{0},z)=P_{\mathrm{0}}\exp (-\mathrm{0.01}z)$ (we remind the reader that, in the local model, z is a given parameter). The values K₀ and P₀ varied from 10⁻⁶ to 10⁻⁵ m² s⁻² with slight changes in transient processes but with the same asymptotic values of K and P at large times. Figure 2 shows the solutions for several depths covering the range shown in Forryan et al. (2013).

Figure 9 Cruise JC29: comparison of profiles of kinetic (a) and potential (b) energies obtained from (3) (local) and (2) (full).

Figure 10 Cruise D306: comparison of profiles of kinetic (a) and potential (b) energies obtained from (3) (local) and (2) (full).

Figure 11 Cruise D321: comparison of profiles of kinetic (a) and potential (b) energies obtained from (3) (local) and (2) (full).

Here, the constant levels of energy are established in several hours. It is natural to use the asymptotic values for comparison with the measurement data. In the subsequent plots we use the log–lin presentation, following Forryan et al. (2013). Figure 3 shows the depth profile of the energies at large times. Note that, according to the second Eq. (4), the asymptotic values of P and K follow each other.

Using this solution, we calculate the turbulent dissipation rate (Eq. 6) and compare it with the data of Forryan et al. (2013) after digitizing both and interpolating them by smooth functions. The result is shown in Fig. 4. Here the difference between theory and measurements is mainly within a half-order. Considering the large spread of experimental data, this is a rather good agreement.

3.2 Cruise D306

To save space, for another two cruises, we show only the asymptotic profiles of the corresponding values at large times. For cruise D306, the values of kinetic and potential energies are of the same order as for cruise JC29 (Fig. 5).

Figure 6 shows the turbulence dissipation rate.

Here again, the difference between the theory and data of Forryan et al. (2013) is within a half-order.

3.3 Cruise D321

The corresponding dependencies for cruise D321 are given in Figs. 7 and 8.

Here again, one can see good agreement between the theory and the mean measured profile.

4 Comparison with the full system

The above results were obtained neglecting vertical turbulence diffusion. To verify this approximation, we solved the full system (2) with the same parameters and initial conditions, adding boundary conditions for fluxes of kinetic and potential energy:

$$\begin{array}{}\text{(7)}& FK={\displaystyle \frac{\sqrt{K}\partial K}{\partial z}},FP={\displaystyle \frac{\sqrt{K}\partial P}{\partial z}}.\end{array}$$

They are given to be compatible with the initial conditions at initial points z₀ from which the plots start in Fig. 2 of Forryan et al. (2013) and tend to zero values at the deepest points. Then, the solutions for K and P were compared with those of the local system (3). The figures below show such a comparison for the asymptotic values. In what follows, the above local solutions are shown in blue and the solutions of Eq. 2 with diffusion in orange.

In all three cases, the local and full models are practically identical. Evidently, this means the closeness of data for the dissipation rate that is a function of K. Hence, for the vertical scales of average values, vertical diffusion can be neglected, and one can use the simplified local Eq. (3).

5 Discussion and conclusions

In this paper we demonstrated that including the potential energy of turbulence (associated with density fluctuations in the presence of stratification) in the semi-empirical Reynolds-type equations of a turbulent flow allows us to explain the existence and evaluate the parameters of small-scale turbulence at large Richardson numbers. Application of these equations to the results of Forryan et al. (2013), where the measurements of profiles of the buoyancy frequency, current shear and dissipation rate of turbulent energy are shown together for three different areas of the ocean, provide not only qualitative but also reasonable quantitative agreement between the theoretical and experimental data. We have also shown that the contribution of turbulent diffusion to the level of turbulent pulsations is insignificant in the above case.

For further progress, more specific experimental data sets are desirable. Indeed, in the above calculations, we used the average data for buoyancy frequency, velocity shear, and rate of kinetic energy dissipation plotted in red in Fig. 2 of Forryan et al. (2013). However, the same figure shows significant, up to 1 order of magnitude, data scatter from each cruise, obtained in different locations and on different days. The authors do not specify the confidence intervals of the data, and no correlation between the different curves for the shear and buoyancy frequency curves is known. Still, it is possible to evaluate the maximal span of the results based on the extreme curves in each plot. For that, we digitized the leftmost and rightmost curves for the shear and buoyancy, interpolated them, and used the resulting functions to calculate the limits of the theoretical dissipation rate using Eqs. (4) and (5). For cruise 306, this is shown in Fig. 12.

Figure 12 Maximal (black) and minimal (gray) profiles of the turbulent kinetic energy dissipation rate for D306 calculated from the maximal possible scatter of data for N2 and S given in Forryan et al. (2013). Red – the average value shown in blue in 6.

Note that the data scattering for this value shown in Forryan et al. (2013) lies within the maximal theoretical limits, which implies that, knowing the real data for a specific location and time of measurement, we would reasonably predict the corresponding depth dependence for the turbulent kinetic energy dissipation rate. Similar results take place for other cruises.

Note in conclusion that the kinetic approach used in the equations used above allows us to naturally include the potential energy under consideration. Considering a large variation of empirical parameters given in different sources (Smyth et al., 2013; Liu et al., 2017; You et al., 2003), for further study it seems important to find more experimental data allowing us to apply the theory. We also plan to extend the present approach to description of time-varying turbulence in the field of internal waves (e.g., Moum et al., 2022).

Appendix A: Dynamical equations for a turbulent stratified flow

Here we briefly outline the general system of equations for a turbulent stratified flow obtained in Ostrovsky and Troitskaya (1987) and developed in Soustova et al. (2020) and Gladskikh et al. (2023). Without dwelling on the details which are described in these works, here we briefly outline the main points of the model. It begins by introducing the variable probability distribution function:

$$\begin{array}{}\text{(A1)}& f(\mathit{v},\mathit{\lambda },\mathit{r},t)=\langle \mathit{\delta }(\mathit{u}-\mathit{v})\mathit{\delta }(\mathit{\rho }-\mathit{\lambda })\rangle ,\end{array}$$

where δ is the Dirac delta function and the angular parentheses denote the ensemble averaging. Using this together with the Navier–Stokes equations for u and ρ and supposing a Gaussian distribution function, the authors of Ostrovsky and Troitskaya (1987) obtained expressions for the average fluxes of turbulent energy, momentum, and mass. They are the same as in the common K−ε theory, except for the mass flux, with the form

$$\begin{array}{}\text{(A2)}& \langle {\mathit{\rho }}^{\prime }{u}_i^{\prime }\rangle =-LV\left({\displaystyle \frac{\partial \langle \mathit{\rho }\rangle }{\partial x_i}}+g_i{\displaystyle \frac{\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle }{V^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}}}-{\displaystyle \frac{\mathit{g}{\mathit{\beta }}_i}{V^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}}}\right),\end{array}$$

where $V=\sqrt{\langle {u^{\prime }}^{\mathrm{2}}\rangle }$, and L, as above, is the outer scale of turbulence, g is the gravity acceleration, and β_i are the components of the vector

$$\begin{array}{}\text{(A3)}& \mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\int \mathrm{d}{\mathit{r}}_{\mathrm{1}}{\displaystyle \frac{\partial }{\partial \mathit{r}}}{\displaystyle \frac{\mathrm{1}}{\left|\mathit{r}-{\mathit{r}}_{\mathrm{1}}\frac{\partial }{\partial z_{\mathrm{1}}}\langle {\mathit{\rho }}^{\prime }(\mathit{r},t){\mathit{\rho }}^{\prime }({\mathit{r}}_{\mathrm{1}},t)\rangle \right|}},\end{array}$$

which characterizes the effect of pressure fluctuations arising from random displacements of a particle in a stratified fluid.

Equation (A2) for the mass flux includes the summand $g_i\frac{\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle }{V^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}}-\frac{\mathit{g}{\mathit{\beta }}_i}{V^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}}$, which, as shown below, leads to some significant differences from results obtained within the framework of known gradient models. The physical meaning of the additional terms in the expression mentioned above is related to the dependence of the force acting upon a random displacement of a liquid particle in a stratified medium on the shape of the liquid volume, i.e., on the ratio of the characteristic scales L_z and L_r. As shown in Ostrovsky and Troitskaya (1987), the components of the vector β have the form ${\mathit{\beta }}_x={\mathit{\beta }}_y=\mathrm{0},{\mathit{\beta }}_z=\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle R$ for a statistically homogeneous field of density fluctuations. Here, R is the anisotropy parameter:

$$\begin{array}{}\text{(A4)}& R\simeq \left\{\begin{array}{ll}\mathrm{1},& L_z\ll L_r,\\ \simeq {\left(\frac{L_r}{L_z}\right)}^{\mathrm{2}}& L_z\gg L_r,\end{array}\right.\end{array}$$

where L_z and L_r are the vertical and horizontal scales, respectively, of the density field correlation.

As a result, one obtains the equations for the mean values of velocity, density, turbulent kinetic energy $K=\mathrm{3}{V}^{\mathrm{2}}/\mathrm{2}$, and variance of density pulsations $\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle$:

$$\begin{array}{}\text{(A5)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \langle u_i\rangle }{\partial t}}+\langle u_j\rangle {\displaystyle \frac{\partial \langle u_i\rangle }{\partial x_j}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial \langle p\rangle }{\partial x_i}}+g_i{\displaystyle \frac{\langle \mathit{\rho }\rangle -{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}}\\ & ={\displaystyle \frac{\partial }{\partial x_j}}\left(L\sqrt{K}\left({\displaystyle \frac{\partial \langle u_i\rangle }{\partial x_j}}+{\displaystyle \frac{\partial \langle u_j\rangle }{\partial x_i}}\right)\right),\end{array}\text{(A6)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \langle \mathit{\rho }\rangle }{\partial t}}+\langle u_i\rangle {\displaystyle \frac{\partial \langle \mathit{\rho }\rangle }{\partial x_i}}\\ & =\mathrm{2}{\displaystyle \frac{\partial }{\partial x_i}}L\sqrt{K}\left({\displaystyle \frac{\partial \langle \mathit{\rho }\rangle }{\partial x_i}}+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}K{\mathit{\rho }}_{\mathrm{0}}}}\left(g_i\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle +\mathit{g}{\mathit{\beta }}_i\right)\right),\end{array}\text{(A7)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial K}{\partial t}}+\langle u_i\rangle {\displaystyle \frac{\partial K}{\partial x_i}}-L\sqrt{K}{\left({\displaystyle \frac{\partial \langle u_i\rangle }{\partial x_j}}+{\displaystyle \frac{\partial \langle u_j\rangle }{\partial x_i}}\right)}^{\mathrm{2}}\\ & -{\displaystyle \frac{\mathit{g}}{{\mathit{\rho }}_{\mathrm{0}}}}L\sqrt{K}\left({\displaystyle \frac{\partial \langle \mathit{\rho }\rangle }{\partial z}}+{\displaystyle \frac{\mathrm{3}\mathit{g}}{\mathrm{2}K{\mathit{\rho }}_{\mathrm{0}}}}\left(\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle +{\mathit{\beta }}_z\right)\right)+{\displaystyle \frac{C{K}^{\mathrm{3}/\mathrm{2}}}{L}}\\ & ={\displaystyle \frac{\mathrm{5}}{\mathrm{3}}}{\displaystyle \frac{\partial }{\partial x_i}}\left(L\sqrt{K}{\displaystyle \frac{\partial K}{\partial x_i}}\right),\end{array}\text{(A8)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle }{\partial t}}+\langle u_i\rangle {\displaystyle \frac{\partial \langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle }{\partial x_i}}\\ & -\mathrm{2}{\displaystyle \frac{\partial \langle \mathit{\rho }\rangle }{\partial x_i}}L\sqrt{K}\left({\displaystyle \frac{\partial \langle \mathit{\rho }\rangle }{\partial x_i}}+\left(g_i\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle -\mathit{g}{\mathit{\beta }}_{\mathit{i}}\right){\displaystyle \frac{\mathrm{3}}{\mathrm{2}K{\mathit{\rho }}_{\mathrm{0}}}}\right)\\ & +{\displaystyle \frac{D{K}^{\mathrm{1}/\mathrm{3}}}{L}}\langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle ={\displaystyle \frac{\partial }{\partial x_i}}L\sqrt{K}{\displaystyle \frac{\partial \langle {{\mathit{\rho }}^{\prime }}^{\mathrm{2}}\rangle }{\partial x_i}}.\end{array}\end{array}$$

In an incompressible fluid considered here, ∇u=0. The potential energy of fluctuations is determined from the last Eq. (A8). Equation (A6) is a particular case of this system for the given average current and density stratification. In general, such effects as internal wave damping by turbulence can be included in the solution as well. On the other hand, the turbulence “breakdown” phenomenon, in which, in certain phases of the wave, the velocity shear cannot maintain a nonzero level of turbulent energy obtained using the common semi-empirical equations (Ivanov et al., 1983), does not exist here. This is also confirmed by numerical calculations using the parametrization based on the model above, given in the work (Gladskikh et al., 2023).