Internal gravity waves in the energy and flux budget turbulence-closure theory for shear-free stably stratified flows
N. Kleeorin I. Rogachevskii I. A. Soustova Yu. I. Troitskaya O. S. Ermakova S. Zilitinkevich





gary@bgu.ac.il
http://www.bgu.ac.il/ gary










^1Department of Mechanical Engineering, Ben-Gurion University of the Negev, P. O. B. 653, Beer-Sheva
 8410530, Israel
 

^2Nordita, Stockholm University and KTH Royal Institute of Technology, 10691 Stockholm, Sweden
 

^3Institute for Atmospheric and Earth System Research, University of Helsinki, 00014 Helsinki, Finland
 

^4Institute of Applied Physics of the Russian Academy of Sciences, 603950 Nizhny Novgorod, Russia
 

^5Finnish Meteorological Institute, 00101 Helsinki, Finland
 

^6Faculty of Geography, Moscow State University, 119234 Moscow, Russia


^7Tyumen State University, 625003 Tyumen, Russia




We have advanced the energy and flux budget (EFB) turbulence closure
theory that takes into account a two-way coupling between internal
gravity waves (IGW) and the shear-free stably stratified turbulence.
This theory is based on the budget equation for the total (kinetic plus
potential) energy of IGW, the budget equations for the kinetic and
potential energies of fluid turbulence, and turbulent fluxes of potential temperature
for waves and fluid flow. The waves emitted at a certain level, propagate
upward, and the losses of wave energy cause the production of turbulence energy.
We demonstrate that due to the nonlinear effects more intensive waves produce
more strong turbulence, and this, in turns, results in strong damping of IGW.
As a result, the penetration length of more intensive waves is shorter than
that of less intensive IGW. The anisotropy of the turbulence produced by less
intensive IGW is stronger than that caused by more intensive waves.
The low amplitude IGW produce turbulence consisting up to 90 % of
turbulent potential energy. This resembles the properties of the observed high
altitude tropospheric strongly anisotropic (nearly two-dimensional) turbulence.




Internal gravity waves in the energy and flux budget turbulence-closure theory
for shear-free stably stratified flows
    S. Zilitinkevich^3,5,6,7
    December 30, 2023
=====================================================================================================================



§ INTRODUCTION


The classical theory of atmospheric flows is based on seminal papers by Rayleigh, Richardson, Prandtl, Kolmogorov, Obukhov and Monin (see, e.g., Refs. <cit.>). This theory implies that any turbulent flow can be considered as a superposition of the fully organised mean flow, and the fully chaotic turbulence characterised by the forward energy cascade from the larger eddies to smaller, resulting in the viscous energy-dissipation at the smallest eddies, with a pronounced inertial interval in the energy spectrum that is fully determined by the energy-dissipation rate <cit.>. The local characteristics of turbulence, in particular, turbulent fluxes that appear in the mean-flow equations, are generally controlled by the local features of the mean flow. It is also assumed that the turbulent flux of any quantity can be expressed as a product of the mean gradient of the quantity multiplied by a turbulent-exchange coefficient. This concept of down-gradient transport reduces the closure problem to determination of the turbulent-exchange coefficients, usually taken proportional to the turbulent kinetic energy and timescale (see, e.g., Refs. <cit.>). This has been formulated for neutrally stratified flows.

However, atmospheric flows often include, besides Kolmogorov turbulence, another type of motions, associated with the development of large-scale structures, e.g, large-scale coherent semi-organized structures (i.e., cloud cells and cloud streets) in turbulent convection <cit.> and internal gravity waves in stably stratified turbulence <cit.>.
The majority of efforts in development of the turbulence closure models for meteorological applications over half a century have been limited to ”mechanical closures" based on the sole use of the turbulent kinetic energy equation, disregarding turbulent potential energy, and daring only cautious corrections to the concept of down-gradient transport.

In stable stratification, such closures have resulted in the erroneous conclusion that shear-generated turbulence inevitably decays and that the flow becomes laminar in ”supercritical" stratifications (at Richardson numbers exceeding some critical value, see, e.g., <cit.>). Obvious contradictions of this conclusion via the well-documented universal presence of turbulence in strongly “supercritical” conditions typical of the free atmosphere and the deep ocean (see, e.g., Refs. <cit.>) were attributed to some unknown mechanisms and, in practical applications, mastered heuristically (see overviews in Refs. <cit.>). The decay of strongly stably stratified turbulence in direct numerical simulations (DNS) has been explained by the effect of diminishing the effective Reynolds number, which comes into play in the not-high-Reynolds-number flows in DNS, but remains insignificant in the very-high-Reynolds-number atmospheric flows.

These principal problems call for a revision of the traditional theory of atmospheric turbulence. The strongest motivation for the revision comes from the need to improve the modern numerical weather prediction, air quality, and climate models, in which turbulent planetary boundary layers couple the atmosphere with underlying land/water/ice surfaces. Stably stratified turbulence determines the vertical turbulent transport of energy and momentum and the turbulent diffusion of pollutants, aerosols, and other admixtures in the free atmosphere.

Numerous alternative turbulence closures in stratified turbulence have been formulated
using the budget equations for various turbulent parameters (in addition to the turbulent
kinetic energy) together with heuristic hypotheses and
empirical relationships (see reviews <cit.>).
Two-point turbulent closures have been developed as well (see reviews <cit.>),
which take into account a very detailed scale-by-scale and directional anisotropy,
that is almost lost in single-point closures.

Key ingredients of stably stratified turbulence are internal gravity waves.
In atmospheric and oceanic turbulence they have been a subject of intense research
(see, e.g., <cit.>).
In the atmosphere, internal gravity waves exist
at scales ranging from meters to kilometers, and are measured by direct probing
or remote sensing using radars and lidars.
The sources of internal gravity waves can be flows over complex terrain, strong wind
shears, convective and other local-scale motions underlying the stably stratified
layer, and wave-wave interactions <cit.>.

The different nature of fluctuations caused by turbulence and waves in stratified flows has been pointed out in
<cit.>.
The role of waves in turbulence closure models has been discussed in <cit.>
An additional negative term in the TKE budget equation (the rate of transfer of TKE into potential energy of wave-like motions) has been included in <cit.>. It was noted that with increasing stability, the wave-like motions dominate in comparison with velocity and buoyancy fluctuations of stratified turbulence, and fluctuations caused by waves suppress vertical mixing (see also <cit.>).

Analysis of the budgets of the wave kinetic energy and the wave temperature variance has been done in <cit.>.
They found significant buoyant production of the wave energy despite the strong static stability and energy transfer from waves to turbulence.
Different aspects related to the effects of
internal gravity waves (IGW) have been also discussed in <cit.>, where
it has been stressed that in the limit of small Froude number,
IGW only affect a poloidal part of the flow.
A strong toroidal cascade coexisting with weak IGW cascade
has been found as well <cit.>.

Numerous high resolution DNS of stably stratified turbulence
with IGW have been performed as well (see, e.g., <cit.>).
The role of IGW and their interaction with the
large-scale flow of vertically sheared horizontal winds has been studied in <cit.>.
It has been shown that most of the energy is concentrated along a
dispersion relation that is Doppler shifted by the horizontal winds. They pointed out that when uniform winds
are let to develop in each horizontal layer of the flow, waves whose phase velocity is equal to the horizontal
wind speed have negligible energy, which indicates a nonlocal transfer of their energy to the mean flow.
Scaling laws for mixing and dissipation in unforced stratified turbulence have been found in <cit.>.
Three regimes characterised by Froude number, namely (i) dominant waves, (ii) eddy-wave
interactions and (iii) strong turbulence have been observed in <cit.>.
An interaction between large-scale IGW and turbulent layer above
the pycnocline (where the density gradient is largest)
has been studied using DNS in <cit.>.
They have demonstrated that in the absence of IGW, the turbulence decays
and most of the turbulent energy is concentrated at the pycnocline center.
The turbulent eddies are collapsed
in the vertical direction and acquire the “pancake” shape.
The internal gravity waves significantly increase turbulent energy
<cit.>.

The energy- and flux-budget (EFB) theory of turbulence closure for stably stratified dry atmospheric flows has been recently developed in Refs. <cit.>. In accordance with wide experimental evidence, the EFB theory shows that high-Reynolds-number turbulence is maintained by shear in any stratification, and the ”critical Richardson number", treated over decades as a threshold between the turbulent and laminar regimes, factually separates two turbulent regimes: the strong turbulence, typical of atmospheric boundary layers, and the weak three-dimensional turbulence, typical of the free atmosphere or deep ocean and characterized by sharp decrease in heat transfer in comparison to momentum transfer. The principal aspects of the EFB theory have been verified against scarce data from the atmospheric experiments, direct numerical simulations, large-eddy simulations (LES) and laboratory experiments relevant to the steady state turbulence regime.

In stably stratified turbulence, large-scale internal gravity waves cause additional vertical turbulent flux of momentum and additional productions of turbulent kinetic energy (TKE), turbulent potential energy (TPE) and turbulent flux of potential temperature <cit.>. For the stationary, homogeneous regime, the EFB theory in the absence of the large-scale IGW yields universal dependencies of the flux Richardson number, the turbulent Prandtl number, the ratio of TKE to TPE, and the normalised vertical turbulent fluxes of momentum and heat on the gradient Richardson number, Ri (see Refs. <cit.>). Due to the large-scale IGW, these dependencies lose their universality.  The maximal value of the flux Richardson number (universal constant  0.2-0.25 in the no-IGW regime) becomes strongly variable in the turbulence with large-scale IGW. In the vertically homogeneous stratification, the flux Richardson number increases with increasing wave energy and can even exceed 1.
The large-scale internal gravity waves also reduce anisotropy of turbulence. Indeed, in contrast to the mean wind shear, which generates only horizontal component of the turbulent kinetic energy, IGW generate both horizontal and vertical components of TKE. These waves increase the share of TPE in the turbulent total energy (TTE = TKE + TPE). A well-known effect of IGW is their direct contribution to the vertical transport of momentum. Depending on the direction (downward or upward), IGW either strengthen or weaken the total vertical flux of momentum. Predictions from this theory <cit.> are consistent with available data from atmospheric and laboratory experiments, DNS and LES.

In the theory discussed in <cit.>,
the stably stratified turbulence is produced
by a large-scale wind shear. This theory takes only
into account a one-way coupling corresponding to
the effect of large-scale IGW with random phases
on stably stratified turbulence,
while the feedback effect of the
turbulence on IGW has been ignored in <cit.>.
In view of applications, the theory discussed in <cit.> describes only
turbulence in the lower troposphere (up to the altitudes about 1 – 1.5 km).

The goal of the present study is to investigate the two-way coupling between large-scale IGW and shear-free stably stratified turbulence. This implies that the turbulence is produced solely by dissipation of IGW propagating in stratified turbulent flows.
In the analysis, we use the budget equations for the kinetic and
potential energies of both, fluid turbulence and large-scale IGW with random phases.
We also apply the budget equations for turbulent heat flux and momentum.
We demonstrate that due to the nonlinear effects the penetration length of the more intense IGW is less than that for the less intensive IGW (with lower wave energy). The low amplitude IGW produce turbulence consisting up to 90 % of potential energy.
The results of the present study describe only
the upper troposphere (located at the altitudes about 10 – 15 km),
see, e.g., <cit.>, and references therein.

This paper is organized as follows.
In Section II we outline properties of the internal gravity waves
propagating in a fluid in the absence of turbulence. We also discuss here
the energy budget equations for IGW.
In Section III we formulate governing equations for the energy and flux budget
turbulence-closure theory for stably stratified turbulence with large-scale IGW.
In Section IV we study the effects of large-scale IGW on turbulence for the steady-state
regime.
The two-way coupling between turbulence and large-scale IGW is analysed in Section V.
Finally, conclusions are drawn in Section VI.




§ LARGE-SCALE IGW IN THE STABLY STRATIFIED FLOWS


In this study we focus on the effect of large-scale internal
gravity waves on the stably stratified turbulence.



 §.§ Linear IGW in the absence of turbulence


First we outline properties of the internal
gravity waves propagating in a fluid in the absence of turbulence and neglecting dissipations.
These waves are described by the following equations:

    ∂Ṽ^ W∂ t + (Ṽ^ W·∇) Ṽ^ W   =    - ∇P̃^ Wρ_0 + β Θ̃^ W e ,
     
    ∂Θ̃^ W∂ t + (Ṽ^ W·∇) Θ̃^ W   =     - β^-1 N^2  Ṽ^ W_j e_j,

where Ṽ^ W, Θ̃^ W  and P̃^ W are the velocity, potential temperature and pressure characterizing IGW, e is the vertical unit vector, β=g/T_0 is the buoyancy parameter,  g is the acceleration due to gravity,  N= (β |∇_z Θ|)^1/2 is the Brunt-Väisälä frequency, Θ is the potential temperature of fluid defined as Θ= T  (P_0/ P)^1-1/γ, where T is the absolute temperature, T_0  is its reference value, P is the fluid pressure, P_0  is its reference value,  γ=c_p/c_v =1.41 is the specific heats ratio, and ρ_0 is the fluid density. The potential temperature Θ̃^ W for waves is defined in a similar way.
Equations (<ref>) and (<ref>) are written in the Boussinesq approximation with the continuity equation,
div Ṽ^ W =0.

Solution of linearized equations (<ref>) and (<ref>) is given by

    Ṽ^ W   =   ( e -  k_h k_z  k_h^2)   V_∗^ W(z)  cos[φ( r) - ω t],
     
    Θ̃^ W   =   N^2(z) ωβ  V_∗^ W(z)  sin[φ( r) - ω t],
     
    P̃^ W   =    - ρ_0  (k_z  ω k_h^2)   V_∗^ W(z)  cos[φ( r) - ω t],

(see, e.g., <cit.>), where  φ( r) is the wave phase, k= k_h +  e k_z  is the wave vector; k_h = (k_x, k_y)  is the horizontal wave vector,
V_∗^ W(z) is the wave velocity amplitude,
and the frequency of IGW is given by

    ω = N(z) k_h  k .

Propagation of IGW in an inhomogeneous medium is determined in the approximation of geometrical optics by the following Hamiltonian equations:

    ∂ r∂ t = ∂ω∂ k,
    
    ∂ k∂ t = - ∂ω∂ r,

(see, e.g., <cit.>), where r is the radius-vector of the center of the wave packet. Since the Brunt-Väisälä frequency N=N(z) depends on the vertical coordinate, Eq. (<ref>) yields k_h= const and the vertical component of the wave vector depends on z, i.e., k_z=k_z(z).  Equations (<ref>) and (<ref>) for the IGW with a fixed frequency allow us to determine k(z) as

    k(z)=k_0  N(z)  N(z_0),

where z=z_0  is the hight where the IGW source is located and k_0=k(z=z_0).
Equation (<ref>) determining the z-dependence of k(z)=(k_z^2 + k_h^2)^1/2  implies that the vertical wave numbers, k_z(z), change when the IGW propagates through the stably stratified flow.
When the IGW source is located at the surface (for z_0=0), the vertical wave number is
k_z(z)=k_h [N^2(z)/ω^2 - 1]^1/2. When the IGW source is located at the upper boundary of the layer under consideration (for z_0=H), the vertical wave number is given by k_z(z)=-k_h [N^2(z)/ω^2 - 1]^1/2.
It follows from these expressions that k(z) = k_h N(z)/ω.



 §.§ Budget equation for total energy of IGW


Let us study the effect of large-scale IGW on the stably stratified turbulence.
We assume that the wavelengthes and periods of the large-scale IGW are
much larger than the turbulence spatial scales and timescales, respectively.
This allows us to treat the large-scale IGW with respect to turbulence
as a kind of mean flows with random phases.
We neglect small molecular dissipation of IGW for large Reynolds and Péclet numbers.
At the low frequency part of the IGW spectra, we limit our analysis to frequencies
essentially exceeding the Coriolis frequency.

We consider the flow fields as a superposition of three components: mean fluid velocity,
U(t,z), fluctuations of fluid velocity, u(t, r), and random large-scale IGW velocity field,  V^ W(t, r), i.e., the total velocity is v= U +  u +  V^ W.
To determine the random large-scale wave velocity field,  V^ W(t, r) ≡⟨ v⟩ - ⟨ v⟩,
besides the ensemble averaging over turbulent fluctuations (denoted by the angular brackets), we also perform an averaging in time over the IGW period (denoted by overbar).
We also use similar decompositions for the total potential temperature,  Θ_ tot =Θ + θ + Θ^W, and for the total pressure,  P_ tot=P + p + P^ W, where
P and Θ are the mean potential temperature and pressure, respectively,
p and θ are fluctuations of the potential temperature and pressure, respectively, and
the wave fields are determined as Θ^ W(t, r) = ⟨Θ_ tot⟩ - ⟨Θ_ tot⟩ and P^ W(t, r) = ⟨ P_ tot⟩ - ⟨ P_ tot⟩.
We assume that the mean fields depend on z-coordinate and time, while the large-scale wave fields depend on all three coordinates and can be represented as an ensemble of wave packets with narrow frequency range and random phases.

In this study we consider a shear-free turbulence, and for simplicity we assume that the mean fluid velocity is zero. We consider low-amplitude approximation for the large-scale IGW and neglect the wave-wave interactions, but take into account the interactions between turbulence and the large-scale IGW.
Equations for the large-scale IGW in stably stratified turbulence are given by

    ∂ V^ W_i ∂ t   =    - ∇_i P^ Wρ_0 + β Θ^ W e_i - ∇_j τ_ij^ W,
     
    ∂Θ^ W∂ t   =    - β^-1 N^2   V^ W_j e_j - ∇_j F_j^ W ,

where the interaction between the turbulence and the large-scale IGW is described by the effective Reynolds stress tensor, τ_ij^ W = ⟨ v_i v_j ⟩ - ⟨ v_i v_j ⟩,
and the effective flux of potential temperature, F_i^ W = ⟨ v_i Θ_ tot⟩ - ⟨ v_i Θ_ tot⟩.
Equations (<ref>)–(<ref>) are mean-field equations,
where the divergence of the Reynolds stress
in the mean momentum equation, and the divergence of the turbulent
heat flux in the mean potential temperature equation
determine effects of turbulence on the mean fields.
To derive Eqs. (<ref>)–(<ref>) for IGW, we obtain two system of equations:
the first system is obtained by the ensemble averaging of the exact momentum
and potential temperature equations over turbulent fluctuations (denoted by the angular brackets);
the second system is obtained by the additional time averaging of the first system of
mean-field equations over the IGW periods (denoted by overbar).
Equations (<ref>)–(<ref>) are obtained by the subtraction of the second system
of equations from the first system.

Using Eqs. (<ref>)–(<ref>), we derive budget equation for the total wave energy,
E^ W=E_K^ W+E_P^ W,

    ∂ E^ W∂ t + ∇·Φ^ W = -D^ W ,

where E_K^ W=( V^ W)^2/2 is the turbulent kinetic wave energy, E_P^ W=(β/N)^2  (Θ^ W)^2/2 is the turbulent potential wave energy, Φ^ W = ρ_0^-1 F^ W  P^ W is the flux of the total wave energy, E^ W, and
D^ W is the dissipation rate of large-scale IGW given by

    D^ W   =     - β^2  N^2(F_j^ W ∇_j Θ^ W - F_j^ W Θ^ W ∇_j ln N^2 )
    
           - τ_ij^ W ∇_j V^ W_i,

is the dissipation rate of the large-scale IGW that is determined by the work of turbulence caused by the interaction with the large-scale IGW. In the gradient approximation, the dissipation rate of the large-scale IGW is positive, D^ W>0. Indeed, as has been shown in <cit.> (see Eq. (3.72) in <cit.>), the flux Φ^ W of the total wave energy, E^ W, is given by Φ^ W = V^ g E^ W, where V^ g=±ω ( k/k^2 -  k_h/k_h^2) is the group velocity of the large-scale IGW.

Let us obtain a steady-state solution of budget equation (<ref>)
for the total wave energy considering two cases:

(i) Non-dissipative large-scale IGW, i.e., the dissipation rate of the total wave energy vanishes, D^ W=0. In this case steady-state solution of Eq. (<ref>)
reads: ∇ · ( V^ g E^ W)=0.
Since k_h= const and ω= const, we find that the vertical profile of the wave amplitude in the solution for Ṽ^ W(t, r) [see Eq. (<ref>)] is given by:

    |V_∗^ W(z)|=|V_0|   [N^2(z)/ω^2 - 1]^-1/4 ,

where V_0 is the constant of integration.
This profile is in agreement with that obtained in <cit.>.

(ii) Dissipative large-scale IGW.
In this case the dissipation rate of the total wave energy is

    D^ W= K_H(1 +  Pr__T) k^2 E^ W ,

which is derived for a narrow wave packet, and
a homogeneous and isotropic turbulence (see Appendix <ref>), where
Pr__T=K_M/K_H is the turbulent Prandtl number, K_M is the eddy viscosity
and K_H is the eddy diffusivity.
Thus, the steady-state solution of Eq. (<ref>)
reads: ∇ · ( V^ g E^ W)=-D^ W, which can be rewritten as

    d Φ^ W_z  dz   =    - σ_g Φ^ W_z ,

where Φ^ W_z= V^ g_z E^ W and

    σ_g    =    K_H(1 +  Pr__T) k_h^3 (N^2 - ω^2)^-1/2(N ω)^4 .

Equation (<ref>) yields

    |V_∗^ W(z)|=|V_0|   |N^2(z)/ω^2 - 1|^-1/4exp(- τ_g  2) ,

where τ_g=∫_0^z σ_g(z')  dz'. In the absence of turbulence (K_M, K_H → 0), the parameter τ_g tends to 0, so that Eq. (<ref>) coincides with
Eq. (<ref>). It follows from Eq. (<ref>) that in the absence of turbulence in the vicinity of a resonance, ω =N(z=z_r), the amplitude of the large-scale IGW tends to infinity. This implies that the low amplitude approximation does not valid, and the nonlinear effects (e.g., the wave braking) should be taken into account. For instance, the wave braking can cause an additional production of turbulence. On the other hand, in the presence of turbulence the infinite growth of the wave amplitude does not occur if the first and the second spatial derivatives of the Brunt-Väisälä frequency vanish at the surface z=z_r <cit.>.



§ ENERGY AND FLUX BUDGET (EFB) MODEL FOR TURBULENCE WITH LARGE-SCALE IGW


In this section we formulate the budget equations for stably stratified turbulence
with large-scale IGW.



 §.§ Budget equations for turbulence with large-scale IGW


In the framework of energy and flux budget turbulence  model, we use
the budget equations for the turbulent kinetic energy (TKE), E_K=⟨ u^2 ⟩/2, the intensity of the potential temperature fluctuations, E_θ=⟨θ^2 ⟩/2, and the vertical turbulent flux, F_z = ⟨ u_z θ⟩, of potential temperature
accounting for large-scale IGW:

    DE_K  Dt + ∇_z  Φ_K   =   β  F_z - ε_K
    - τ_ij^ W ∇_j V^ W_i + β V^ W_z Θ^ W,
    
    
    D E_θ Dt + ∇_z  Φ_θ   =    - F_z  ∇_z Θ - ε_θ - F_j^ W ∇_j Θ^ W,
     
    D F_z  Dt + ∇_z  Φ_F    =   β ⟨θ^2
    ⟩ - 1ρ_0 ⟨θ ∇_z p ⟩ - ⟨ u_z^2 ⟩ ∇_z  Θ
           - ε_z^(F) - τ_j3^ W ∇_j Θ^ W
    - F_j^ W ∇_j V^ W_z,

(see, e.g., Ref. <cit.>),
where D / Dt = ∂ /∂ t +  U·∇.
The terms Φ_K, Φ_θ and Φ_F include
the third-order moments. In particular, Φ_K = ρ_0^-1⟨ u_z   p⟩ + (1/2) ⟨ u_z   u^2 ⟩ determines the flux of E_K;
Φ_θ = ⟨ u_z  θ^2 ⟩/2 determines the flux of E_θ and Φ_F = ⟨ u_z^2 θ⟩ + ρ_0^-1 ⟨θ  p ⟩ / 2 determines the flux of F_z.
The terms ε_K=E_K/t_T, ε_θ=E_θ/(C_p   t_T) and
ε^(F)=F_z/(C_F   t_T) are the dissipation rates of the turbulent kinetic energy E_K, the intensity of the potential temperature fluctuations E_θ and the vertical turbulent heat flux F_z. These  dissipation rates are expressed using the Kolmogorov  hypothesis <cit.>, where t_T=ℓ_0 /E_K^1/2 is the turbulent dissipation timescale, ℓ_0 is the integral scale of turbulence, and C_p  and C_F are dimensionless constants. The turbulent potential energy (TPE), E_P, is defined as E_P = β^2 E_θ / N^2.

The third term in the right hand side of Eq. (<ref>) contributes to the traditional
gradient turbulent flux (proportional to - ∇_z  Θ) of potential temperature,
while the first and the second terms in the right hand side
of Eq. (<ref>) describe a non-gradient contribution to the vertical turbulent flux of
potential temperature. In stably stratified flows
the gradient and non-gradient contributions to the vertical turbulent flux of
potential temperature have opposite signs. This implies that
the non-gradient contribution decreases the traditional gradient turbulent flux.

The budget equations for the turbulent kinetic energies,
E_α = ⟨ u_α^2⟩/2, along the x, y
and z directions can be written as follows:

    DE_α Dt + ∇_z  Φ_α   =   δ_α 3 β  F_z - ε_αα^(τ) + 1  2 Q_αα
    - τ_α j^ W ∇_j V^ W_α ,

where α=x,y,z, the term ε_αα^(τ) = E_α/3t_T
is the dissipation rate of the turbulent kinetic energy components,
Φ_α= ρ_0^-1⟨ u_α  p⟩ + (1/2) ⟨ u_z   u_α^2 ⟩ determines the flux of E_α, the term Q_αα
is correlations between the fluctuations of the pressure, p, and the small-scale velocity shears: Q_ij = ρ_0^-1 (⟨ p ∇_i u_j⟩ + ⟨ p ∇_j u_i⟩).
In Eq. (<ref>) we do not apply the summation convention for the double Greek indices.
Equations (<ref>)–(<ref>) are obtained by averaging over the ensemble of turbulent fluctuations and over the period of large-scale IGW. These equations without the large-scale IGW terms can be found in <cit.> (see also <cit.>).

Equations (<ref>)–(<ref>) contain production terms caused
by the large-scale IGW (see <cit.>).
In particular, the terms - τ_ij^ W ∇_j V^ W_i
+ β V^ W_z Θ^ W in the right hand side of Eq. (<ref>)
determine the production rate of TKE, and the term F_j^ W ∇_j Θ^ W
in Eq. (<ref>) contributes to the production rate of TPE, while the terms
- τ_j3^ W ∇_j Θ^ W
- F_j^ W ∇_j V^ W_z in Eq. (<ref>) describe
the production rate of the vertical turbulent flux, F_z, of potential temperature.
To close the system of the budget equations (<ref>)–(<ref>), one needs to determine
the Reynolds stress for the wave fields, τ_ij^ W, and the turbulent flux of potential temperature
for the wave fields, F^ W, which have been derived in <cit.>:

    τ_ij^ W   =    - C_τ  t_T (τ_im∇_m V^ W_j + τ_jm∇_m V^ W_i ) ,
     
    
    F_i^ W   =    - C_F   t_T (τ_im∇_m Θ^ W + τ_i3^ W∇_z Θ + F_m∇_m V^ W_i ) .

These quantities are caused by interactions between the large-scale IGW and turbulence. They are determined by subtracting of the ensemble-averaged equations (but not averaged over the IGW period) from exact equations for these quantities, assuming that ω t_T ≪ 1.
To derive Eqs. (<ref>) and (<ref>), we
assume that the effective dissipation rate, ε_ij^(τ,W), of
the Reynolds stress for the wave fields, τ_ij^ W, and
the effective dissipation rate, ε_i^(F,W),
of the turbulent flux of potential temperature
for the wave fields, F^ W, are expressed as
ε_ij^(τ,W) =
τ_ij^ W/(C_τ  t_T) and ε_i^(F,W) =
F_i^ W/(C_F   t_T), respectively, where C_τ is a dimensionless constant.
In Eqs. (<ref>) and (<ref>) we also
omit the terms which are quadratic in wave amplitude,
because they do not contribute to the correlations
τ_ij^ W ∇_j V^ W_i
and  τ_ij^ W ∇_j Θ^ W
entering in Eqs. (<ref>)–(<ref>) <cit.> (see also <cit.>).

If there is no separation of timescales between the
turbulent scales and the IGW scales, there are additional correction factors
in the right hand side of Eqs. (<ref>) and (<ref>). In particular,
the right hand side of Eq. (<ref>) is multiplied by the factor
[1 + (C_τω t__T)^2]^-1, while
the right hand side of Eq. (<ref>) is multiplied by the factor
[1 + (C_F ω t__T)^2]^-1.
However, since the free constants C_τ and C_F are small, i.e.,
C_τ^2 ≈ C_F^2 ≈ 10^-2 (see Sect. IV),
these correction factors are of the order of 1 even when there is
no separation of the timescales.
On the other hand, any mean-field theory formally requires
a separation of scales for its validity.



 §.§ Summary of assumptions and steps of derivations


In our analysis, we use budget equations for the one-point second moments
for the following reasons.
We develop a mean-field theory and do not study small-scale structure
of turbulence. In particular, we study  large-scale long-term dynamics, i.e.,
we consider effects in the
spatial scales which are much larger than the integral scale of turbulence
and in timescales which are much longer than the turbulent timescales.
We limit our analysis to the geophysical flows, in which the vertical variations
of the mean fields are much larger than the horizontal variations, so that the
terms associated with the horizontal gradients in the budget equations for turbulent
statistics can be neglected. For instance, in typical atmospheric flows, the vertical
scales are much smaller than the horizontal scales, so that the mean vertical velocity
is much smaller than the horizontal velocity.

We have made the following assumptions related to internal gravity waves:


  * The periods and wavelengths of IGW are larger than the turbulent correlation time
and the integral scale of turbulence.


  * We assume that the large-scale wave fields
can be represented as an ensemble of wave packets with narrow
frequency range and random phases.


  * We restrict the analysis to a low-amplitude approximation for the large-scale IGW
and neglect the wave-wave interactions, but take into account the interactions
between turbulence and the large-scale IGW. We leave an
account for the wave-wave interactions for further study.


  * The ensemble of the large-scale IGW has a power-law spectrum
and is isotropic in the horizontal plane.


  * We neglect small molecular dissipation of IGW considering turbulence with
large Reynolds and Péclet numbers.
At the low frequency part of the IGW spectra, we limit our analysis to frequencies
essentially exceeding the Coriolis frequency.



We apply the energy and flux budget turbulence closure model, which
assumes the following:


  * The characteristic times of variations of
the turbulent kinetic energy, the turbulent potential energy,
and the vertical turbulent flux F_z of potential temperature
are much larger than the turbulent timescale.
This allows us to obtain a steady-state solution
for the budget equations of TKE, TPE and F_z
for a homogeneous stratified hydrodynamic turbulence.
These budget equations include production terms caused by the large-scale IGW.


  * We neglect the divergence of the fluxes of TKE, TPE and F_z
(i.e., we neglect the divergence of third-order moments).
In the present study, we restrict the analysis to the effects of the
large-scale IGW on the second-order statistics and leave
the third-order moments (the fluxes of energies and
the fluxes of momentum and heat fluxes) for further study.


  * Dissipation rates of TKE, TPE and F_z are expressed using
the Kolmogorov  hypothesis <cit.> (see also <cit.>), i.e.,
ε_K=E_K/t_T, ε_θ=E_θ/(C_p   t_T) and
ε^(F)=F_z/(C_F   t_T).


  * We assume that the term ρ_0^-1 ⟨θ ∇_z p ⟩ in Eq. (20) for the vertical turbulent flux of potential
temperature is parameterised by C̃_θ β ⟨θ^2
⟩ with C̃_θ < 1. This implies that β ⟨θ^2
⟩ - ρ_0^-1 ⟨θ ∇_z p ⟩ = C_θ β ⟨θ^2
⟩ with the positive coefficient C_θ = 1 -C̃_θ that is less than 1.
The justification of this assumption is discussed
in <cit.>.


  * The considered stratified hydrodynamic turbulence with large-scale IGW
is shear-free and isotropic in the horizontal plane.


In the next section we will use Eqs. (<ref>)–(<ref>) to study effects of the large-scale IGW on turbulence.



§ EFFECTS OF LARGE-SCALE IGW ON TURBULENCE: STEADY-STATE REGIME


In this section we study the effects of large-scale internal gravity waves on turbulence,
using the steady state version of Eqs. (<ref>)–(<ref>). Solving the system of these equations,
we obtain the turbulent kinetic energy E_K and the vertical turbulent flux F_z
of the potential temperature:

    E_K    =   2 C_F  3 (ℓ_0 N_0)^2  Ri__ W Ŵ ,
     
    
    F_z    =    -K_H  N_0^2 β ,

where N_0=N(z=z_0) and the coefficient of turbulent (eddy) diffusivity is

    K_H    =   (2 C_F  3)^3/2 ℓ_0^2   N_0 (Ŵ - 1) ( Ri__ W Ŵ)^3/2 .

Here we have introduced the following two key parameters characterizing
effects of the large-scale IGW on turbulence:

(i) the wave Richardson number:

    Ri__ W = N_0^2   H^2 (Γ(μ)   Pr_T^(0))   E^ W ;


(ii) the parameter Ŵ:

    Ŵ≡- τ_ij^ W ∇_j V^ W_i E_K/t_T = C_τ Ψ(Q,μ)  E^ W ℓ_0^2  E_K   H^2 ,

where E^ W is the total wave energy, H is the hight of the layer where the large-scale IGW propagate,
the dimensionless functions Γ(μ) and Ψ(Q,μ) are given below, μ
is the exponent of the energy spectrum of the large-scale IGW and Q=[N(z)/N(z_0)]^2 is the dimensionless lapse rate.

The reason for using these parameters for characterization of the effects of large-scale IGW
on turbulence is as follows.
We consider a shear-free turbulence. The turbulence is produced only by large-scale IGW, so that
the classical gradient Richardson number, Ri_g=N^2/S_ wind^2 tends to infinity,
because the mean wind shear, S_ wind, vanishes.
The most appropriate parameter in this case is the effective Richardson number
Ri__ W associated with the amplitude of the wave.
For instance, the large wave Richardson number implies a low-amplitude wave.

Note that the Froude number Fr = U/(L N_0) used in fluid dynamics is
related to the wave Richardson number Ri__ W as Fr =  Ri__ W^-1/2,
where the velocity U = √(E^W) is related to the wave energy E^ W and the scale
L = H  (Γ(μ)   Pr_T^(0))^-1/2 is proportional
to the hight H of the layer in which the large-scale IGW propagate.
We use the notion of the wave Richardson number following tradition of the atmospheric physics
and meteorology, where different kinds of Richardson numbers (e.g., the gradient Richardson number
and the flux Richardson number) are used.

Another important parameter, Ŵ, is defined as the ratio
of the TKE production rate caused by the internal gravity waves to the dissipation rate of TKE.
This parameter describes the efficiency of the turbulence production by waves.
It depends on the wave Richardson number, Ŵ=Ŵ( Ri__ W),
and it can be interpreted as the squared ratio of the wave shear, √(E^ W) / H,
and the turbulent shear, √(E_K) / ℓ_0.
The function Ŵ( Ri__ W) is determined by the cubic algebraic equation:
Ŵ^3 + B_1 Ŵ^2 + B_2 Ŵ + B_3 =0, with coefficients B_k given in Appendix <ref>.

The function Γ(μ) in Eq. (<ref>) depends on the exponent μ of the IGW energy spectrum.
In particular, for a power law, k^-μ, energy spectrum of IGW existing in the range of the wavenumbers,
H^-1≤ k ≤ L_W^-1, the function Γ(μ) is given by: Γ(μ)=|C_μ|
(H/L_W)^3-μ for 1<μ<3; Γ(μ)=4 ln(H/L_W) for μ=3 and
Γ(μ)=C_μ for μ >3, where C_μ=2(μ-1)/(μ-3) <cit.>.
The function Ψ(Q,μ) in Eq. (<ref>) that is related to the parameter Ŵ,
is given by

    Ψ(Q,μ) = 2  3[1 + 3 (Q-1)A_z] Γ(μ) .


Now using the steady-state version of Eqs. (<ref>)–(<ref>),
we determine various dimensionless parameters versus the wave Richardson number.
Assuming for simplicity constant Brunt-Väisälä frequency which implies that Q=1,
we obtain the vertical anisotropy parameter, A_z = E_z/E_K, and the ratio of turbulent potential and kinetic energies, E_P/ E_K:

    A_z    =   5 - Ŵ  (1- A_z^∗)   (1- Ri_ f^∗)  5 + 2 Ŵ (1- A_z^∗)   (1- Ri_ f^∗) ,
    
    E_P E_K   =    C_P  [Ŵ(1 + 1/ Pr__T^(0))-1],

where E_P = (β/N)^2 E_θ is the turbulent potential energy, A_z^∗ and Ri_ f^∗ are the vertical anisotropy parameter and the flux Richardson number in the limit of very large gradient Richarson number in the absence of the IGW, respectively, and Pr__T^(0) is the turbulent Prandtl number
for a nonstratified turbulence (at zero gradient Richardson number).
Equation (<ref>) determines the vertical anisotropy parameter A_z.
Since A_x + A_y + A_z=1 and A_x = A_y, we obtain that A_x = A_y = (1- A_z)/2,
where the horizontal anisotropy parameters are defined as A_x = E_x/E_K and A_y = E_y/E_K.
The turbulent viscosity, K_M=2 C_τ A_z E_K^1/2ℓ_0, is given by

    K_M    =   (2 C_F  3)^1/2 2 C_τ ℓ_0^2 N_0 [ Ri__ W Ŵ( Ri__ W)]^1/2 ,

so that the turbulent Prandtl number is

    Pr__T = 3 C_τ C_F [ Ri__ W Ŵ( Ri__ W) Ŵ -1] ( C_A - Ŵ C_A + 2 Ŵ) ,

where C_A=5   (1- A_z^∗)^-1  (1- Ri_ f^∗)^-1.
Equations (<ref>)–(<ref>) and (<ref>) allow us to determine the nondimensional ratio F_z^2/ E_K E_θ:

    F_z^2  E_K E_θ   =   2 C_F  3 C_P [(Ŵ -1)^2  Ri__ W Ŵ( Ri__ W)]
    
          ×[Ŵ(1 + 1  Pr__T^(0))-1]^-1 .


Let us discuss the choice of the dimensionless empirical constants in the developed theory. There are two well-known universal constants: the limiting value of the flux Richardson number Ri_ f^∗=0.25 for an extremely strongly stratified turbulence (i.e., at infinite gradient Richardson number)
and the turbulent Prandtl number Pr__T^(0)=0.8 for a nonstratified turbulence (at zero gradient Richardson number).
The vertical anisotropy parameter for an extremely strongly stratified turbulence is A_z^∗ =0.03
and the constants C_F=C_τ/ Pr__T^(0) = 0.125, where C_τ is the coefficient determining the turbulent viscosity  (K_M=2 C_τ A_z E_K^1/2ℓ_0) for a non-stratified turbulence. The constant C_P=0.417 describes the deviation of the dissipation timescale of E_θ=⟨θ^2 ⟩/2 from the dissipation timescale of TKE. The constants C_F, C_P and A_z^∗ are determined from numerous meteorological observations, laboratory experiments and LES (see details in <cit.>).
The results essentially depends on the empirical constant C_θ, e.g.,
the constant C_θ=1/15 is chosen to get a small parameter A_z^∗ to reproduce a quasi-two-dimensional turbulence for an extremely strongly stratified turbulence.

In Fig. <ref> we show the dependence of the vertical anisotropy parameter A_z= E_z/E_K on the wave Richardson number Ri__ W for different values of the parameter C_θ.
The anisotropy parameter depends strongly on the empirical constant C_θ, e.g., it becomes negative, A_z<0, when C_θ<1/15.
This indicates that the system does not reach a steady-state.






The ratio of the turbulent potential energy, E_P, to the total turbulent energy E=E_K+E_P versus the wave Richardson number Ri__ W is shown in Fig. <ref>. It is seen here that for large wave Richardson number (i.e., for the low amplitude IGW), the turbulent potential energy is about 90 % of total turbulent energy. Without waves, the ratio E_P/E is less than 0.2, while in the presence of the large-scale IGW this ratio can reach 0.9.

The non-dimensional turbulent kinetic energy E_K/(ℓ_0 N_0)^2 versus Ri__ W [see Eq. (<ref>)] is shown in Fig. <ref>. It is clearly seen in this figure, that the turbulent kinetic energy decreases rapidly with the increase of Ri__ W, and E_K is less than the Ozmidov energy scale, (ℓ_0 N_0)^2.
Note that the Ozmidov length scale, √(ε_K/N_0^3),
is well known as a rough threshold between
anisotropic scales and isotropic ones, where
ε_K is the dissipation rate of TKE.



In Figs. <ref> and <ref> we show the turbulent Prandtl number Pr__T and the non-dimensional squared potential temperature flux F_z^2/ E_K E_θ versus the wave Richardson number Ri__ W, respectively. It can be seen in these figures that the turbulent Prandtl number increases with decrease of the wave energy E^ W, becomes quite large, and the heat transfer becomes weaker, i.e., the ratio  F_z^2/ E_K E_θ decreases with decrease of the wave energy.







§ TWO-WAY COUPLING BETWEEN TURBULENCE AND LARGE-SCALE IGW


In this section we consider the two-way coupling of stably-stratified turbulence and large-scale IGW.
The large-scale IGW emitted at a certain level, propagate
upward, and the losses of wave energy cause the production of turbulence energy.
Equation (<ref>) for the total wave energy,
E^ W, reads:

    ∂ E^ W∂ t + ∇_z (V^ g  E^ W) = -γ_d E^ W,

where the damping rate, γ_d, of the large-scale IGW is given by

    γ_d = C_F  (1 +  Pr__T^(0))  Ψ(Q,μ)  ℓ_0  H^2  E_K^1/2 ,

V^ g = (1-μ^-1)   N_0   H   f(Q) is the group velocity of the large-scale IGW and the function f(Q) is

    f(Q) = (Q-1)^1/2 Q ∫_0^π/2(1 + cos^2 ϑ Q-1)^1/2 sin^2 ϑ d ϑ .


The steady-state version of Eq. (<ref>) can be written in the form of the nonlinear Byger equation, ∇_z I = - κ  I, where I=V^ g  E^ W and κ=γ_d/V^ g.
This equation can be also rewritten as

    ∇_z  Ri__ W = √( Ri__ W) H_ eff Ŵ( Ri__ W)  ,

where the effective damping length scale H_ eff is defined as

    H_ eff = C_μ H (H ℓ_0)^2 ,

and C_μ=2(μ-1)/(μ-3) depends on the exponent μ of the energy spectrum of the large-scale IGW.
Equation (<ref>) allows us to obtain the spatial profile of the wave Richardson number.
In particular, the function Ri__ W(z), as the result of
solution of Eq. (<ref>) is shown in Fig. <ref>.
As follows from Eq. (<ref>) and Fig. <ref>, the effective damping length scale H_ eff is much larger than the equilibrium height H. This implies that the large-scale IGW penetrate almost the entire atmosphere height and generate weak turbulence. These waves are significantly attenuated at z = 10 H_ eff.



Equation for the total turbulent energy E=E_K + E_P is

    ∂ E ∂ t + ∇_z Φ = Π - ε,

where Π = γ_d E^ W is the production rate of the total turbulent energy E caused by the damping of the large-scale IGW, and ε=E/(C_P t_T) is the dissipation rate of the total turbulent energy. Equations (<ref>) and (<ref>) yield the budget equation for the sum, E+E^ W of the turbulent total energy and the wave total energy:

    ∂(E + E^ W) ∂ t + ∇_z (Φ + V^ g  E^ W)= - E  C_P t_T .

This equation describes energy exchange between turbulence and large-scale internal gravity waves.

Let us analyze the two-way interaction between large-scale IGW and stably-stratified turbulence.
In Fig. <ref> we show the vertical profile of the normalised turbulent kinetic energy, E_K(z)/(ℓ_0 N_0)^2. It decreases rapidly with height reaching the value 2 × 10^-3 at z = 10 H_ eff.
The vertical profile of the normalized squared turbulent flux F_z^2/ E_K E_θ of potential temperature (see Fig. <ref>) is similar to the vertical profile of the turbulent kinetic energy, shown in Fig. <ref>.







The vertical profile of A_z is shown in Fig. <ref>, which demonstrates that A_z decreases with height.
This implies that anisotropy of turbulence increases with height, and this effect is more significant for C_θ=  1/15. For small wave Richardson numbers (i.e., for intensive waves), turbulence is almost isotropic, but for large heights the turbulence anisotropy becomes more significant compared to the case of a sheared stably stratified turbulence.
In particular, when hight z varies from z= 0.1   H_ eff to z = 50 H_ eff, the parameter A_z decreases in 10 times (see Fig. <ref>).
The latter implies formation of a ”pancake" structure in turbulent velocity field for large z (see, e.g., <cit.>).





The turbulent Prandtl number and the ratio of potential to the total energy E_P/E shown in Figs. <ref> and <ref>,
increase with height as the wave becomes less intensive.
It is seen in Fig. <ref> that for the upper layers, the turbulent potential energy is larger than the  turbulent kinetic energy. This implies that the temperature fluctuations dominate over the velocity fluctuations.
In particular, the ratio E_p/E increases up to 0.85 at z=50 H_ eff.
For comparison, in a shear-produced stably stratified turbulence without IGW, the ratio of turbulent potential energy to total turbulent energy, E_p/E, reaches 0.15 at very large gradient Richardson numbers (see Fig. 7 in <cit.>).
On the other hand, in a shear-produced stably stratified turbulence with IGW (where only the one-way
coupling is taking into account), the maximum ratio E_p/E
increases up to 0.45 at very large gradient Richardson numbers
(see Fig. 5 in <cit.>).



§ CONCLUSIONS




  * Within the new EFB turbulence closure model a system of equations describing two-way interactions between internal gravity waves (IGW) and turbulence has been derived. This system includes the budget equation for the total (kinetic plus potential) energy of IGW, and the equations for the kinetic and potential energies of turbulence, turbulent heat fluxes for waves and flow for arbitrary stratification. The general physical picture in two-way coupling between IGW and turbulence is following: waves emitted at a certain level (in case if there is no refraction), propagate upward. The losses of wave energy cause the production of turbulence. Therefore, waves transfer energy, while turbulence causes its losses. We have shown that more intensive waves penetrate in a shorter distance, while less intense IGW penetrate in larger distances. This is caused by the nonlinear effects, where the more intense IGW produce more strong turbulence that results in more intensive damping of IGW.

  * The analysis of the effects of IGW on the anisotropy of turbulence has shown that for less intensive waves the turbulence anisotropy is stronger. Low amplitude waves produce anisotropic turbulence with a low energy, and the total turbulent energy consists up to 90 % of potential energy.  This property resembles one observed in high altitude tropospheric nearly two-dimensional turbulence.

  * We also have demonstrated that the kinetic energy of turbulent fluctuations has the Ozmidov energy scale
    (which is the product of the Brunt-Väisälä frequency and turbulence integral scale). When turbulence is produced only by waves, the gradient Richardson number tends to infinity, because it is a shear-free turbulence (for which the wave shear is much larger than the wind shear). The most appropriate parameter in this case is the effective Richardson number associated with the amplitude of the wave.
    
 



The detailed comments on our manuscript by the two anonymous referees
which essentially improved the presentation of our results,
are very much appreciated.
NK and IR acknowledge support from
the Israel Science Foundation governed by the
Israel Academy of Sciences (grant No. 1210/15).
IR also acknowledges the Research Council of Norway under the FRINATEK (grant  No. 231444)
and the National Science Foundation (grant No. NSF PHY-1748958).
IS and YT acknowledge support from Russian Foundation for Basic Research (RFBR grants
No. 18-05-00292, No. 18-05- 00265, No. 17-05-41117, No. 19-05-00366).
SZ acknowledges support from the Academy of Finland (grant ClimEco No.
314 798/799, FMI-FI), Russian Science Foundation (RSF grant No. 15-17-20009,
IPFAN-RU), Russian Foundation for Basic Research (RFBR grants No. 18-05-60299,
IPFAN-RU, and No. 18-05-60126, INM/MSU-RU), Russian Science Foundation
(RSF grant No. 14-17-00131, Tyumen State University-RU).
NK and IR acknowledge the hospitality of Nordita and
Institute for Atmospheric and Earth System Research (INAR) of
University of Helsinki. IR also acknowledges the hospitality of
the Kavli Institute for Theoretical Physics in Santa Barbara.








§ DERIVATION OF EQ. (<REF>)



Equation (<ref>) follows from Eq. (<ref>). In particular,
substituting Eqs. (<ref>)–(<ref>) and (<ref>)–(<ref>) into Eq. (<ref>),
performing averaging, and taking into account that for the linear IGW
E_K^ W = E_P^ W = (1/2) E^ W,
we arrive at expression (<ref>) for the dissipation rate D^ W
of the total wave energy E^ W.

A simple and approximate derivation of Eq. (<ref>) is as follows.
The dissipation rate of the wave kinetic energy E_K^ W is
D^ W_K = - K_M  V^ W·Δ V^ W,
which in k space reads 2 K_M k^2 E_K^ W.
The dissipation rate of the wave potential energy E_P^ W is
D^ W_P= - (β^2/N^2) K_H Θ^ W ΔΘ^ W,
which in k space reads 2 K_H k^2 E_P^ W.
Taking into account that for the linear IGW,
E_K^ W = E_P^ W = (1/2) E^ W,
we arrive at the expression (<ref>) for the dissipation rate D^ W=D^ W_K+D^ W_P
of the total wave energy E^ W.



§ THE FUNCTION Ŵ( RI__ W)



The function Ŵ( Ri__ W) is determined by the following cubic algebraic equation:
Ŵ^3 + B_1 Ŵ^2 + B_2 Ŵ + B_3 =0, where


    B_1    =    - [1 - 2 C_θ C_P  + C_A (C_θ C_P  (1 + 1/ Pr__T^(0)) - C_F  5)
    + 2  3 Ri__ W]  [ C_F - 2 C_θ C_P  (1 + 1/ Pr__T^(0))]^-1,
     
    
    B_2    =   [C_A (1 + C_θ C_P) + 2 - C_A  3  Ri__ W]
    [ C_F - 2 C_θ C_P  (1 + 1/ Pr__T^(0))]^-1,
     
    
    B_3    =   C_A  3  Ri__ W [C_F - 2 C_θ C_P  (1 + 1/ Pr__T^(0))]^-1 ,


and C_A= 5 (1- A_z^∗)^-1(1- Ri_ f^∗)^-1.




MY71
A. S. Monin and  A. M. Yaglom, Statistical
Fluid Mechanics (MIT Press, Cambridge,
Massachusetts, 1971), Vol. 1.

MY75
A. S. Monin and A. M. Yaglom, Statistical Fluid
Mechanics  (MIT Press, Cambridge, Massachusetts,
1975), Vol. 2.

K41
A. N. Kolmogorov, Doklady AN SSSR 32, No. 1, 19 (1941).

Z91
S. S. Zilitinkevich, Turbulent Penetrative Convection  (Avebury Technical, Aldershot, 1991).

EB93
D. Etling and R. A. Brown, Boundary-Layer Meteorol. 65, 215 (1993).

AZ96
B. W. Atkinson and J. Wu Zhang, Rev. Geophys. 34, 403 (1996).

Z02
S. S. Zilitinkevich, Quart. J. Roy. Meteorol. Soc. 128, 913 (2002).

EKRZ02
T. Elperin, N. Kleeorin, I. Rogachevskii and S. Zilitinkevich,
Phys. Rev. E 66, 066305 (2002);
Boundary-layer Meteorol. 119, 449 (2006).

B74
T. Beer, Atmospheric Waves (Wiley, New York, 1974).

GH75
E. E. Gossard and W. H. Hooke, Waves in the Atmosphere (Elsevier, New York, 1975).

M01
Yu. Z. Miropol'sky, Dynamics of Internal Gravity Waves in the Ocean
(Springer, Dordrecht, 2001).

N02
C. J. Nappo, An Introduction to Atmospheric Gravity Waves
(Academic Press, London, 2002).

BU09
O. Bühler, Waves and Mean Flows (Cambridge Univ. Press, New York, 2009).

S10
B. R. Sutherland, Internal Gravity Waves (Cambridge Univ. Press, Cambridge, 2010).

C61
S. Chandrasekhar, Hydrodynamic and
Hydromagnetic Stability, (Dover Publications
Inc., New York, 1961), Sect. 2.

M86
J. Miles, Phys. Fluids 29, 3470 (1986).

SF01
E. J. Strang and H. J. S. Fernando, J. Phys. Oceanogr. 31, 2026
(2001).

O01
Y. Ohya, Boundary-Layer Meteorol. 98, 57 (2001).

BAN02
R. M. Banta,  R. K. Newsom,  J. K. Lundquist, Y.
L. Pichugina,  R. L. Coulter and  L. Mahrt,
Boundary-Layer Meteorol. 105, 221 (2002).

PAR02
E. R. Pardyjak,  P. Monti and  H. J. S. Fernando,
J. Fluid Mech. 459, 307 (2002).

MO02
P. Monti,  H. J. S. Fernando,  M. Princevac, W.
C. Chan, T. A. Kowalewski and  E. R. Pardyjak,
J. Atmos. Sci. 59, 2513 (2002).

LA04
J. S. Lawrence, M. C. B. Ashley,  A. Tokovinin
and T. Travouillon, Nature 431, 278 (2004).

M14
L. Mahrt, Annu. Rev. Fluid Mech. 46, 23 (2014); Boundary-Layer Meteorol. 135, 1 (2010).

CAN09
V. M. Canuto, In.: Lect. Notes Phys., vol  756, ”Interdisciplinary aspects of turbulence", W. Hillebrandt and F. Kupka (Eds.) pp. 107-160 (Springer, Berlin, 2009).

WT03
W. Weng and P. Taylor, Boundary-Layer Meteorol. 107, 371  (2003).

UB05
L. Umlauf and H. Burchard H., Continental Shelf Research 25, 725 (2005).

C01
C. Cambon, Eur. J. Mech. B (fluids) 20, 489510 (2001).

SC18
P. Sagaut and C. Cambon, Homogeneous Turbulence Dynamics (Springer, New York,
2018).

IOS83
A. V. Ivanov,  L. A. Ostrovsky,  I. A. Soustova,  L. Sh. Tsimring, Dynamics of Atmosheres and Ocean
7, 221 (1983).

OT87
L. A. Ostrovsky and  Yu.I. Troitskaya, Izv. Akad. Nauk. SSSR, Fiz. Atmos. Okeana 23, 1031 (1987).

GM79
C. Garrett and W. Munk, Annu. Rev. Fluid Mech. 11, 339 (1979).

HM06
K. R. Helfrich and W. K. Melville, Annu. Rev. Fluid Mech. 38, 395 (2006).

JDS06
Q. Jiang, J. D. Doyle and R. B. Smith, J. Atmos. Sci. 63, 617 (2006).

S07
R. B. Smith, J. Atmos. Sci. 64, 594 (2007).

PZ14
R. Plougonven and  F. Zhang, Rev. Geophys. 52, 33 (2014).

SN15
J. Sun, C. J. Nappo, L. Mahrt, D. Belusíc, B. Grisogono, et al.,
Rev. Geophys. 53, 956 (2015).

YM15
E. Yigit and A. S. Medvedev, Advances in Space Res. 55,  983 (2015).

C99
G. Chimonas, Boundary-Layer Meteorol. 90, 397 (1999); ibid 102, 139 (2002).

SS02
C. Staquet and  J. Sommeria, Annu. Rev. Fluid Mech. 34, 559 (2002).

FA03
D. C. Fritts and  M. J. Alexander, Rev. Geophys. 41, 1003 (2003).

JRF05
F. G. Jacobitz, M. M. Rogers and  J. H. Ferziger, J. Turbulence 6, 1  (2005).

JSG03
L. H. Jin, R. M. C. So and T. B. Gatski,
J. Fluid Mech. 482, 207 (2003).

BP04
H. Baumert and H. Peters, J. Phys. Oceanogr. 34, 505 (2004);
Ocean Sci. 5, 47 (2009).

ST04
C. Staquet, Surv. Geophys. 25, 281 (2004).

FE81
J. J. Finnigan and F. Einaudi, Quart. J. Roy. Met. Soc. 107, 807 (1981).

FEF84
J. J. Finnigan, F. Einaudi and D. Fua, J. Atmosph. Sci. 41, 2409 (1984).

EF84
F. Einaudi, J. J. Finnigan and D. Fua, J. Atmosph. Sci. 41, 661 (1984).

F88
J. J. Finnigan, J. Atmosph. Sci. 45, 486 (1988);
Boundary-Layer Meteorol. 90, 529 (1999).

EF93
F. Einaudi and J. J. Finnigan, J. Atmosph. Sci. 50, 1841 (1993).

LI06
E. Lindborg, J. Fluid Mech. 550, 207 (2006).

CM15
P. Clark di Leoni and P. D. Mininni, Phys. Rev. E 91, 033015 (2015).

PR18
A. Pouquet, D. Rosenberg, R. Marino and C. Herbert, J. Fluid Mech. 844, 519 (2018).

DOZ13
O. A. Druzhinin, L. A. Ostrovsky and S. S. Zilitinkevich,
Nonlin. Proc. Geophys. 20, 977 (2013).

DO15
O. A. Druzhinin and L. A. Ostrovsky,
Nonlin. Proc. Geophys. 22, 337 (2015).

ZKR07
S. S. Zilitinkevich,  T. Elperin,  N. Kleeorin
and I. Rogachevskii, Boundary-Layer Meteorol.
125, 167 (2007).

ZKR08
S. S. Zilitinkevich,  T. Elperin,  N. Kleeorin,
I. Rogachevskii, I. Esau, T. Mauritsen and M.
Miles, Quart. J. Roy. Met. Soc. 134, 793 (2008).

ZKR09
S. S. Zilitinkevich,  T. Elperin,  N. Kleeorin,
V. L'vov and I. Rogachevskii, Boundary-Layer Meteorol. 133, 139 (2009).

ZKR10
S. S. Zilitinkevich, I. Esau, N. Kleeorin, I. Rogachevskii
and R. D. Kouznetsov, Boundary-Layer Meteorol.
135, 505 (2010).

ZKR13
S. S. Zilitinkevich, T. Elperin, N. Kleeorin, I.
Rogachevskii, and I.  Esau,
Boundary-Layer Meteorol. 146, 341 (2013).

T73
J. S. Turner, Buoyancy Effects in Fluids
(Cambridge University Press, Cambridge, 1973).

W62
S. Weinberg, Phys. Rev. 126, 1899 (1962).

KF94
J. C. Kaimal and  J. J. Fennigan, Atmospheric Boundary Layer Flows (Oxford
University Press, New York, 1994).

CM93
V. M. Canuto and  F. Minotti, J. Atmos. Sci. 50, 1925
(1993).

CCH02
Y. Cheng, V. M. Canuto and A. M. Howard, J., Atmosph.
Sci. 59, 1550 (2002).





