## Abstract

The buoyancy subrange of stably stratified turbulence is defined as an intermediate range of scales larger than those in the inertial subrange. This subrange encompasses the crossover from internal gravity waves (IGWs) to small-scale turbulence. The energy exchange between the waves and small-scale turbulence is communicated across this subrange. At the same time, it features progressive anisotropization of flow characteristics on increasing spatial scales. Despite many observational and computational studies of the buoyancy subrange, its theoretical understanding has been lagging. This article presents an investigation of the buoyancy subrange using the quasi-normal scale elimination (QNSE) theory of turbulence. This spectral theory uses a recursive procedure of small-scale modes elimination based upon a quasi-normal mapping of the velocity and temperature fields using the Langevin equations. In the limit of weak stable stratification, the theory becomes completely analytical and yields simple expressions for horizontal and vertical eddy viscosities and eddy diffusivities. In addition, the theory provides expressions for various one-dimensional spectra that quantify turbulence anisotropization. The theory reveals how the dispersion relation for IGWs is modified by turbulence, thus alleviating many unique waves' features. Predictions of the QNSE theory for the buoyancy subrange are shown to agree well with various data.

## 1. Introduction

Most geophysical and planetary flows are affected by stable stratification that gives rise to internal gravity waves (IGWs). These flows are also turbulent, and the interaction between turbulence and waves has significant impact upon the transport of momentum and scalars on different scales [1,2]. Understanding and quantification of this interaction is one of the important unresolved problems of geophysical fluid dynamics and fluid dynamics in general [3,4]. Based upon the observations in the atmospheric boundary layer, Mahrt *et al.* [5] identified three flow regimes that differed by relative importance of the effects of turbulence and stable stratification. The differences could be quantified in terms of the gradient Richardson number, *Ri*, defined as *Ri*=*N*^{2}/*S*^{2}, where *N* is the Brunt–Väisälä frequency defined below and *S* is the mean shear. The regimes were classified as weakly stable, transitional and very stable for , and , respectively. Using the results for the eddy viscosities and eddy diffusivities obtained from the quasi-normal scale elimination (QNSE) theory [6,7], one can establish correspondence between the ranges of *Ri* and the strength of stratification on different scales as measured by the ratio *k*/*k*_{O}, where *k* is a wavenumber, *k*_{O}=(*N*^{3}/*ϵ*)^{1/2} is the Ozmidov wavenumber and *ϵ* is the rate of the viscous dissipation. More specifically, the ranges [0,0.05], [0.05,0.3] and [0.3,) for the Richardson number *Ri* would roughly correspond to (,3], [3,0.3] and [0.3,0) for *k*/*k*_{O}.

Mahrt *et al.* [5] based their regime classification upon the magnitudes of the vertical eddy viscosity and eddy diffusivity. An alternative classification can use the one-dimensional spectra of the kinetic and potential energies. For instance, by analysing velocity and temperature distributions in free atmosphere, Alisse & Sidi [8] diagnosed layers where the vertical spectra of the horizontal velocity and temperature fluctuations scaled as either or , *k*_{z} being the vertical wavenumber. The former was identified with the dynamical regime of Kolmogorov turbulence, while the latter was likened to the IGW-dominated regime in free troposphere and stratosphere [9], where this spectrum was referred to as *canonical* because of its propensity to develop as a universal distribution nearly invariant for all seasons, meteorological conditions and geographical locations throughout the atmosphere. The transition between the two regimes takes place close to the Ozmidov scale, .

Bouruet-Aubertot *et al.* [10] explored classification of flow regimes using vertical and horizontal one-dimensional spectra obtained from oceanic observations. Similar to the atmospheric observations, the Ozmidov scale, *k*^{−1}_{O}, was identified with the transition from the Kolmogorov (referred to as inertial-convective) to the stratified inertial turbulent, or buoyancy, or saturated IGW subrange. For brevity, the latter will be referred to as the buoyancy subrange. As in the case of atmospheric turbulence, the one-dimensional vertical spectra in the Kolmogorov and buoyancy subranges follow the and distributions, respectively [11,12]. The behaviour of the horizontal spectra in these subranges differs from their vertical counterparts as they were suggested to simply preserve their Kolmogorov distribution throughout the buoyancy subrange [13–15].

Note that the and transition in the lower atmosphere was observed much earlier by Shur [16]. These observations stimulated the development of early theories of stratified turbulence [17–23]. A review of some later work was given by Galperin & Sukoriansky [2] who also summarized the application of the QNSE theory to stably stratified turbulence in the limits of weak and strong stratification.

For the case of arbitrary stratification, QNSE is a complicated spectral theory whose main derivations are quite difficult to follow. In the limit of weak stratification, however, its mathematics simplify significantly and the theory becomes completely analytical. This limit is important for understanding the transition from the Kolmogorov to the buoyancy-dominated flow regime and modification of various one-dimensional spectra. Furthermore, a fundamental theory of this transition sheds light upon the universality of the canonical spectrum whose physical nature has not been fully understood. Other important outcomes of this theory are improved understanding and quantification of the energy partition between internal waves, turbulence and small-scale dissipation in atmospheric and oceanic flows. The purposes of this paper are precisely the development of the analytical theory of the transitional Kolmogorov–buoyancy subrange and studies of various properties of turbulence and waves within this subrange.

The analysis will be based upon the QNSE theory developed for the case of arbitrary stratification [6]. In this study, the derivations will be performed in the limit of weak stratification. All the basic assumptions of the theory remain unchanged in this limit as well as all the basic steps of calculations. The theory uses the energy balance in the assumption of stochastically steady state but does not consider the prognostic energy equation and ensuing energy fluxes.

The time-dependent energy equation is considered in a conceptually different spectral theory by Cambon and co-workers [24,25]. These authors develop a fully axisymmetric multi-modal extension of the eddy-damped quasi-normal Markovian (EDQNM) approximation [26] for flows with stable stratification. For rotating and stratified flows, this approach was further developed in Cambon & Scott [27], Cambon [28,29], Godeferd & Staquet [30], Cambon *et al.* [31] and other publications using both EDQNM and direct numerical simulation (DNS). The extensive body of work of the Lyon group has been summarized in the recent book by Sagaut & Cambon [32]. One of the central points of their approach is the use of the Craya–Herring decomposition in which the incompressibility condition, expressed as the orthogonality of the velocity vector and its wavevector **k**, is used via choosing a local coordinate system normal to **k** which encloses the velocity vector. This vector is further decomposed into toroidal and poloidal components, with the former being orthogonal to the gravity vector **g**, i.e. belonging in the horizontal plane and representing the vortical mode. The temperature mode at wavevector **k** (recast as the buoyancy) is formally joined with the velocity vector as the third component of a formal vector . The solenoidal part of thus includes the toroidal and poloidal velocity components, while the pseudo-dilatational part, in the language of compressible flows, is the buoyancy scalar. Linear and nonlinear operators are then expressed in terms of .

Over the years, the Craya–Herring and related wave-vortex decompositions have proved to be a powerful tool for studies of anisotropic turbulent flows [32–36]. They are closely related to the eigenmode expansions in terms of the gravity wave and potential vorticity modes [37]. Among important achievements of the EDQNM formalism combined with the Craya–Herring decomposition are demonstration of spectral anisotropization of turbulence and collapse of the energy of vertical fluctuations under the action of stable stratification. The former manifests itself in energy concentration in the vertical cone in Fourier space and concurrent formation of horizontal layering in physical space, while the latter results in the establishing of a two-component flow [25,32]. As shown in §8, the QNSE theory yields similar results.

Unlike EDQNM, the main focus of QNSE is on the ‘dressed’ viscosities and diffusivities, or, equivalently, the renormalized velocity and temperature Green functions. They are computed by applying the successive coarsening procedure to the coupled system of momentum and temperature equations transformed using the Feynman diagram technique [6]. This approach yields different renormalizations for the viscosity and the diffusivity and generates different ‘dressed’ viscosities and ‘dressed’ diffusivities in the horizontal and vertical directions. In the limit of weak stratification, all ‘dressed’ variables and spectra can be calculated analytically. As will be shown in §8, tendencies to spectral anisotropization emerging from these calculations agree with the EDQNM predictions outlined in Godeferd & Cambon [25] and Sagaut & Cambon [32].

The paper is organized as follows. Section 2 formulates basic equations, and outlines the assumptions and methods of QNSE. Section 3 elaborates the limit of weak stable stratification. Sections 4 and 5 provide detailed calculation of the effective viscosities and diffusivities in the limit of weak stable stratification. These sections are of a technical nature and a reader who is not interested in formal mathematical derivations can proceed directly to §§6 and 7. Section 8 discusses a broad range of spectral characteristics of stably stratified turbulence. Finally, §10 provides discussion and conclusions.

## 2. Model outline

The QNSE, as well as other spectral theories [25], is formulated for a three-dimensional turbulent flow with imposed vertical stabilizing temperature gradient. The flow is governed by the momentum, temperature and continuity equations in the Boussinesq approximation,
2.1
2.2
2.3
where is the unit vertical vector directed upward, *P* is the pressure, *ρ* is the constant reference density, *ν*_{0} and *κ*_{0} are the molecular viscosity and diffusivity, respectively, *β* is the thermal expansion coefficient, *g* is the gravity acceleration directed downward, *Θ* is the mean potential temperature and *T* describes the fluctuations of *Θ*.

In the spectral domain bounded by the Kolmogorov wavenumber, , space–time Fourier transforms of the velocity and temperature fields are 2.4 and 2.5 where and are the Fourier modes of velocity and temperature, respectively, and is a four-dimensional vector in Fourier space. A detailed derivation of the model has been given elsewhere [38,6,7], and only a brief review of its main points will be given here.

QNSE operates with Fourier-transformed governing equations. These equations are strongly nonlinear as the Reynolds number, *Re*, measuring the ratio of nonlinear to viscous terms in the momentum equation is large on large scales. Straightforward application of perturbation methods to obtain large-scale solutions to these equations would yield divergent results. The situation is different when the smallest scales are considered, however, where viscous processes dominate and *Re*=*O*(1) [39]. The smallness of *Re* prompts the exploration of the renormalized perturbation methodology operating with ‘dressed’, or effective, or eddy viscosity and eddy diffusivity rather than their ‘bare’ (i.e. molecular) values [40,41]. This methodology allows one to gradually coarse-grain the flow domain by recursive elimination of small shells of small-scale modes and calculate compensating effective eddy viscosity and eddy diffusivity accounting for the transport processes on the eliminated scales [42,43]. This approach is a cornerstone of the QNSE theory [6,38] whose quintessence is the calculation of the scale-dependent ‘dressed’ viscosities and diffusivities using a coupled system of the momentum and temperature equations. They are written in the form of the Langevin equations that can be used to derive expressions for the second-order moments and various spectra.

To proceed with the development of the QNSE theory, one needs to derive a self-contained equation for the velocity . This is achieved by excluding the pressure term in Fourier-transformed equation (2.1) using the incompressibility equation. The resulting equation still contains the temperature, , given by the Fourier-transformed equation (2.2),
2.6
where
2.7
is the ‘bare’ temperature Green function. The temperature can be excluded from the momentum equation by the use of a formal solution to the temperature equation (2.6) in the form of the Langevin equation with the renormalized temperature Green function ,
2.8
and
2.9
where *κ* is the effective diffusivity, is the stochastic forcing evoked by vertical velocity fluctuations,
2.10
and accounts for nonlinear stochastic stirring of a mode by all other modes. This leads to a self-contained canonical equation for ,
2.11
where is a large-scale stochastic force caused by large-scale instabilities responsible for maintaining statistically steady state,
2.12
and
2.13
*δ*_{βμ} is the Kronecker delta-symbol, is the ‘bare’ tensorial velocity Green function,
2.14
and
2.15
*N*=(*βg*(d*Θ*/d*z*))^{1/2} is the Brunt–Väisälä frequency, *ϕ* is the angle between the vector **k** and the vertical and is the ‘bare’ auxiliary Green function,
2.16
The tensorial form of the Green function is the result of the action of stable stratification. The presence of the complex poles in the denominator of the factor given by equation (2.15) points to the presence of IGWs that coexist with turbulence. Their properties will be elaborated later. The domain of definition of the canonical equations extends from 0 to *Λ*, *Λ*=*k*_{d}.

The recursive scale elimination starts with subdivision of the domain of definition on two sub-domains, *D*^{<}=(0,*Λ*−*ΔΛ*] and *D*^{>}=(*Λ*−*ΔΛ*,*Λ*], where *ΔΛ* is infinitesimal. Correspondingly, ‘slow’ () and ‘fast’ () velocity and temperature modes are defined such that *k*∈*D*^{<} for **u**^{<},*T*^{<} and *k*∈*D*^{>} for **u**^{>},*T*^{>}. Any vector mode *u*_{α} and scalar mode *T* can be represented as *u*_{α}=*u*^{<}_{α}+*u*^{>}_{α} and *T*=*T*^{<}+*T*^{>}. For the fast modes, *Re*=*O*(1). Ensemble averaging of the fast modes over the shell *D*^{>} shrinks the domain of definition from *D*=*D*^{<}∪*D*^{>} to *D*^{<}. Renaming the new dissipation cutoff, *Λ*−*ΔΛ*, into *Λ* renders *Λ* variable and associates it with *dynamic* dissipation cutoff. The averaging generates *O*(*ΔΛ*) corrections to the viscosity and diffusivity which account for the transport processes on now unresolved scales. Emerging from this procedure, ‘dressed’ viscosity and diffusivity are flow-dependent, in contrast to their ‘bare’ constant molecular values appearing in the governing equations (2.1) and (2.2) with all scales resolved.

At the next step, the shrunk domain of definition is again subdivided onto two sub-domains, (0,*Λ*−*ΔΛ*] and (*Λ*−*ΔΛ*,*Λ*], and new groups of slow and fast variables are identified. The procedure of elimination of the fast variables from the shell *ΔΛ* is repeated, yielding new corrections to the viscosity and diffusivity, and so on. Successive repetition of this procedure results in a cyclic process, each iteration of which modifies ‘dressed’ viscosity and diffusivity and shrinks the spectral domain by additional *ΔΛ*. The decrease of *Λ* and increase of eddy viscosity are such that the Reynolds number based upon the parameters specified in the vicinity of *Λ* remains *O*(1), thus signalling that the successive scale elimination procedure is mathematically sound.

The ensemble averaging in the shell *D*^{>} relies upon representation of the velocity modes in that shell via Langevin equation,
2.17
where is the *renormalized* Green function as it contains renormalized viscosity and diffusivity resulting from elimination of small scales up to *Λ*. These viscosity and diffusivity will be calculated later. Note that in all renormalized Green functions, zero indexes are suppressed. The random force is generated by nonlinear interactions as elaborated in the earlier studies [40,44–47]. Numerous attempts to derive this force from first principles have been unsuccessful so far [40,48], and so its form needs to be postulated. Physically, accounts for stirring of a mode *k* by all other modes; it is solenoidal, zero-mean, white noise in time, isotropic and homogeneous in time and space. Its correlation function is thus given by
2.18
where the forcing amplitude *D* is proportional to the energy dissipation rate *ϵ*. As elaborated in Sukoriansky *et al.* [6], the proportionality factor is not adjustable; it is calculated from the energy balance considerations.

As shown in Sukoriansky *et al.* [38,6], an important demand to **f** is that 〈*f*_{α}(**p**)*f*_{β}(**q**)*f*_{σ}(**k**−*p*−**q**)〉=0 for such vector triads that *p*,**q** and **k**−**p**−**q**∈*D*^{>}. This property alone suffices to develop a rigorous self-contained mathematical procedure of successive averaging. This force does not need to be Gaussian although a Gaussian field would meet the above requirement. Generally, could be characterized as quasi-normal. The combination of the quasi-normal forcing and eddy damping represented by the eddy viscosity and eddy diffusivity places QNSE in the class of quasi-normal eddy-damped theories of turbulence [26,32,40,49].

Derivations reveal that not only the Langevin equation (2.17) acquires tensorial form under the action of stable stratification but in addition, corrections to eddy viscosity and eddy diffusivity differ in vertical and horizontal directions. This behaviour is recognized by introducing a set of new variables instead of traditional *νk*^{2} and *κk*^{2},
2.19
and
2.20
where ; *k*_{x}=*k*_{1}; *k*_{y}=*k*_{2}; ; and *k*_{z}=*k*_{3} are the horizontal and vertical wavenumbers, respectively, and
2.21
and
2.22
The initial values on the right-hand side of these equations are *ν*_{0}=*ν*=*ν*_{h}=*ν*_{z}, *κ*_{0}=*κ*=*κ*_{h}=*κ*_{z}, *δν*_{z}=0 and *δκ*_{z}=0. For convenience, the pairs *ν*_{h}, *ν*_{z} and *κ*_{h}, *κ*_{z} will be referred to as horizontal and vertical eddy viscosities and eddy diffusivities, respectively. The dependence on will be understood in further derivations as dependence on three variables, *ω*,*k* and *k*_{3}, i.e. the vertical coordinate, *k*_{3}, is singled out and the expressions for the Green functions become
2.23
2.24
2.25

It is well known that viscosity and diffusivity are well defined when there is a spectral gap between the eliminated ‘subgrid’ scales *k*>*Λ* that generate effective viscosity and diffusivity [39,50] and ‘resolvable’ scales, *k*<*Λ*, upon which viscosity and diffusivity act. In real flows, however, such a gap does not exist, and so one has to introduce two-parametric viscosity and diffusivity that reveal a cusp-like behaviour on the smallest resolvable scales *k*≃*Λ* [51]. Note in this respect that all renormalized Green functions within QNSE are, in fact, two-parametric as they depend on the dynamic dissipation cutoff, *Λ*, and vector . The QNSE theory uses an important simplification known as the distant interaction approximation in which the limit *k*/*Λ*→0 is taken, and only the terms up to *O*(*k*^{2}) are retained [43,47]. De facto, this approximation introduces a spectral gap between the resolvable and subgrid scales and accordingly, renormalized viscosities and diffusivities are taken as functions of *Λ* only. In doing so, the dependence on *k* and, thus, the cusp-like behaviour of the viscosities and diffusivities are lost. The degree of inaccuracy introduced by this approximation cannot be estimated rigorously at the present time. However, its utility has been demonstrated by Smith & Woodruff [47], who showed that it provides an accurate description of the effective viscosity in isotropic three-dimensional turbulence. In addition, Smith & Woodruff [47] provided some heuristic arguments in its justification. The distant interaction approximation will be used throughout this study. For brevity, the parameter *Λ* will be omitted in the arguments of all Green functions.

The corrections to viscosity and diffusivity generated by elimination of the shell *D*^{>} are calculated from the corrections to the inverse Green functions [6],
2.26
and
2.27
where
2.28
and *U*_{αβ}(*Ω*,*q*,*q*_{3}) is the spectral tensor. Evaluation of this tensor using equations (2.17), (2.18), (2.23) and (2.15) yields
2.29
where
2.30

The integrals in (2.26) and (2.27) are calculated in the distant interaction approximation and the limit *ω*→0 as only the long-time behaviour is considered. In §§4 and 5, it will be shown that
2.31
and
2.32
where *Δν*, *Δδν*_{z}, *Δκ* and *Δδκ*_{z} are all *O*(*ΔΛ*). In the limit *ΔΛ*→0, one obtains a system of four coupled ordinary differential equations for *ν*, *δν*_{z}, *κ* and *δκ*_{z} as functions of *Λ*. This system of equations will be derived and solved in the limit of weak stable stratification.

## 3. The limit of weak stable stratification

From the governing equations, one infers that there are two competing processes in the system at a scale *Λ*^{−1}: nonlinear interactions and buoyancy. Their relative importance is measured by the ratio of nonlinear (*ν*(*Λ*)*Λ*^{2}) to buoyancy (*N*) frequencies, giving rise to the spectral Froude number,
3.1

Later, it will be shown that for , *ν*(*Λ*) approximately obeys the Richardson law, *ν*(*Λ*)∼*ϵ*^{1/3}*Λ*^{−4/3}. Substitution of this expression in (3.1) yields
3.2
where
3.3
is the Ozmidov wavenumber. This wavenumber delineates a scale, , at which the turbulence eddy turnover time is equal to the time scale of internal wave, *N*^{−1}. Turbulence dominates on smaller scales, while larger scales are strongly influenced by the waves. Both and *Λ*/*k*_{O} can be used as non-dimensional parameters characterizing the strength of stratification. Equation (2.15) indicates that in the case of weak stratification, asymptotic expansions only contain even powers of *N* such that the first terms in such expansions are small. The derivations of the QNSE theory simplify considerably when only terms are retained. Some of these simplifications can be traced to the term *B* in equation (2.30), where the last term on the right-hand side is small and can be neglected in the limit of weak stratification. This reduces the complexity of the integrands in equations (2.26) and (2.27). Other simplifications follow from the fact that stratification-induced anisotropic corrections *δν*_{z} and *δκ*_{z} to *ν* and *κ*, respectively, are small. This observation will be used later to simplify the derivations.

Alternatively to , expansions can also be carried out in powers of (*Λ*/*k*_{O})^{−4/3}. Note that the early theories of stably stratified turbulence mentioned in §1 dealt with *O*[(*k*/*k*_{O})^{−4/3}] corrections to the Kolmogorov theory, i.e. de facto they considered the limit of weak stratification.

## 4. Computation of eddy diffusivities

As elucidated in the study of Sukoriansky *et al*. [6], evaluation of the integrals in (2.26) and (2.27) consists of the following three steps:

— contour integration over the frequency

*Ω*;— invocation of the distant interaction approximation; and

— integration over the spherical shell

*D*^{>}, which is composed of multiplication by*ΔΛ*and angular integration over the surface of the sphere with the radius*Λ*.

In all calculations, simplifications related to the symmetry and solenoidality properties of the projection operator, *P*_{αβ}(**k**)=*P*_{βα}(*k*) and *k*_{α}*P*_{αβ}(**k**)=0, are routinely enforced. Similarly enforced is the continuity equation, *u*_{β}*k*_{β}=0. Because the integral (2.26) is factored with , any term in that integral proportional to *k*_{β} will be set to zero automatically.

The temperature integral is simpler of the two, and it is convenient to calculate it first. The integrand in equation (2.27) in the limit *ω*→0 is
4.1
The distant interaction approximation implies that the expression on the right-hand side of (2.27) is computed up to the order of *O*(*k*^{2}). The factor *k*_{α}*k*_{β} in front of the integral is already *O*(*k*^{2}), and so only the *O*(1) terms should be retained in the integrand (4.1). Setting *k*=0 and *k*_{3}=0 in *G*_{T} gives
4.2
The use of equation (2.29) renders factorization of ,
4.3
Because equation (4.3) includes scalar factors dependent on *Ω* and tensorial factors independent of *Ω*, it is convenient to represent it as a sum,
4.4
where
4.5
4.6
4.7
4.8

### (a) Frequency integration

#### (i) Contribution from *S*^{1}(*Ω*,*q*,*q*_{3})

This term possesses three poles, 4.9 where 4.10 and 4.11

Two of these poles, *Ω*_{1} and *Ω*_{2}, lie in the upper half-plane, while *Ω*_{3} belongs in the lower half-plane. The frequency integral *R*^{1} is computed using the residue at the pole *Ω*_{3}, yielding
4.12
Using the smallness of *δν*_{z} and *δκ*_{z} compared with *ν* and *κ*, this expression can be simplified to read
4.13

#### (ii) Contribution from *S*^{2}(*Ω*,*q*,*q*_{3})

This term possesses seven poles, three of which are given by equation (4.9), and the rest, owing to the factor *B*(*Ω*,*q*,*q*_{3}), are
4.14
Four of these seven poles, *Ω*_{1}, *Ω*_{2}, *Ω*_{4} and *Ω*_{6}, are located in the upper half-plane, while the remaining three belong in the lower half-plane. The frequency integral *R*^{2} is calculated using the residues at the latter poles. In the limit of weak stratification, the result is
4.15

The result of the frequency integration in (2.27) can be written as 4.16

### (b) Distant interaction approximation

Recall that in this approximation, only the terms up to *O*[(*k*/*Λ*)^{2}] are retained in thus facilitating isolation of *O*(*k*^{2}) terms that generate corrections to effective viscosities and diffusivities and as additional benefit, leading to simplification of the angular integration. This approximation has already been applied to the integrand (4.1) in which *k* and *k*_{3} were set to zero owing to the presence of the *O*(*k*^{2}) factor, *k*_{α}*k*_{β}, in front of the integral in (2.27). The contractions of this factor with and are, thus, *O*(*k*^{2}) and cannot be further simplified. The respective expressions are
4.17
and
4.18

### (c) Integration over the spherical shell *D*^{>}

The integration of (4.16) over a spherical shell *D*^{>} consists of multiplication of the integrand by −*ΔΛ* (the negative sign is stipulated by contraction of the domain when scales are eliminated) and integration over the surface of a *d*-dimensional sphere, *Σ*_{d}, with a radius *Λ* (for convenience, some of the calculations are performed in *d*-dimensional space; *d*=3). In neutrally stratified isotropic turbulence, the angular integration is straightforward and yields sums of products of the Kronecker *δ*-symbols [38,43]. In flows with stable stratification, however, the integrand (4.16) is anisotropic and the procedure of angular integration requires further development. Inspection of this integrand reveals that it is composed of the products of terms, some of which depend on *q*_{3} only (these terms are due to the residues *R*^{1} and *R*^{2}) while others, produced by the tensorial functions and , involve a polynomial dependence on the scalar products of the type *k*_{α}*q*_{α}. The latter terms preserve the isotropy in the horizontal plane. Taking advantage of this, one can perform the angular integration over the surface of a *d*-dimensional sphere, *Σ*_{d}, by, firstly, integrating over the surface of a *d*−1-dimensional sphere *Σ*_{d−1} (for *d*−1=2, such a sphere is a circle in the horizontal plane) and, secondly, integrating over the vertical semicircle:
4.19
Note that in the integration over *Σ*_{d−1}, the horizontal radius *q*_{h} is a function of the vertical coordinate given by
4.20
where *θ* is the angle between the vector **q** and the vertical axis; *q* is replaced by *Λ* reflecting the fact that *D*^{>} is a spherical shell of the radius *Λ*.

Owing to the horizontal isotropy of the integrand (4.16), the integration over *Σ*_{d−1} can be performed analytically after representing vector **q** as a sum of its horizontal, **q**_{h}=(*q*_{1},*q*_{2},0), and vertical, *q*_{3}, projections,
4.21
where
4.22
The integration over a *d*−1-dimensional spherical surface of radius *Λ* can be accomplished by using the following auxiliary relationships valid for the horizontal components of the vector **q**:
4.23
4.24
4.25
where *S*_{d}=2*π*^{d/2}/*Γ*(*d*/2) is the surface area of a *d*-dimensional unit sphere, and *α* through *δ* are indexes in the horizontal plane.

Integrations over the circle *Σ*_{d−1} and the angle *θ* can be performed analytically, and the result is structurally identical to equation (2.32). The correction to is composed of two parts, one of them is proportional to *k*^{2} and the other to . They provide *O*(*ΔΛ*) corrections to *κ* and *δκ*_{z} denoted as *Δκ* and *Δδκ*_{z}, respectively. Upon dividing these corrections by *ΔΛ*, taking the limit *ΔΛ*→0, and substituting *d*=3, one obtains differential equations for *Λ*-dependent *κ* and *δκ*_{z},
4.26
and
4.27
These equations are coupled and their right-hand sides involve *ν* and *δν*_{z} such that they require two additional equations. These equations will be derived in §5.

## 5. Computation of the effective viscosities

Equations for *ν* and *δν*_{z} will be derived in this section, using exactly the same steps and algorithms as those used above to obtain the effective diffusivities.

### (a) Frequency integration

It is convenient to bring the factor *P*_{αμθ}(**k**) inside the integral in (2.26) and represent the integrand as
5.1

Substitute equations (2.14) for the Green function and (2.29) for the velocity correlator in (5.1) and represent the integrand as a factorized sum,
5.2
where the tensorial terms are
5.3
5.4
5.5
5.6
and the scalar terms, in the limit *ω*→0, are
5.7
5.8
5.9
5.10
where , *A*=*A*(−*Ω*,|**k**−**q**|,*k*_{3}−*q*_{3}) and *B*=*B*(*Ω*,*q*,*q*_{3}).

Because *A* and *B* are both small, the term *S*^{4} can be neglected, leading to a simplified equation (5.2),
5.11
This is a significant simplification as the neglected terms produced the most complicated expressions in the full solution [6].

Similar to the case of the temperature equation (4.4), only the scalar functions *S*^{i} are *Ω*-dependent; they generate a set of poles in the frequency integration in (2.26) whose contributions can be calculated using standard contour methods. The tensorial terms, on the other hand, only affect the angular integration. The contour integration is performed next.

#### (i) Contribution from *S*^{1}

This term produces three poles,
5.12
where *ω*_{q} is given by equation (4.10) and
5.13

Two of these poles, *Ω*_{1} and *Ω*_{3}, are located in the upper half-plane, while *Ω*_{2} belongs in the lower half-plane. The frequency integral *R*^{1}(*k*,*k*_{3},*q*,*q*_{3}) (for brevity, the arguments of *R*^{i} will be suppressed thereon) is computed using the residue at pole *Ω*_{2},
5.14
For weak stratification and using distant interaction approximation, this equation yields
5.15

#### (ii) Contribution from *S*^{2}

The product of and *A* produces, in addition to the three poles of *S*^{1}, two more complex poles in the upper half-plane,
5.16
where
5.17

Out of the five poles given by (5.12) and (5.16), four are located in the upper and one in the lower half-plane. The frequency integral *R*^{2} is computed using the residue at pole *Ω*_{2} and using the weak stratification assumption,
5.18
The distant interaction approximation yields further simplification to (5.18),
5.19

#### (iii) Contribution from *S*^{3}

The function *B* given by (2.30) has the same poles as *A* and its conjugate *A** except that now **k**=0 in the argument of *A*,
5.20
Thus, *S*^{3} has seven poles as given by equations (5.12) and (5.20); four of these poles are located in the upper and three in the lower half-plane. The frequency integral is computed using the residues at the lower half-plane poles. In the limit of weak stratification, the result is
5.21
Further simplification of this equation using the distant interaction approximation yields
5.22

Upon completing frequency integration, the integrand *I*_{αβ} takes the form
5.23

### (b) Distant interaction approximation

Similarly to §4, the integrand (5.23) is simplified by retaining only terms up to *O*[(*k*/*Λ*)^{2}]. Unlike the case of the temperature equation, however, where the integral was factored with an *O*(*k*^{2}) term and the integrand only needed to be evaluated at the zero order, in the present case, the integral in (2.26) is only factored with *P*_{αμθ}(**k**), which is *O*(*k*). Thus, each of the factors and *R*^{i} in (5.23) should be computed up to *O*(*k*). Computational techniques used in this section are similar to those used in §4 for calculation of the effective diffusivities.

The integrand (5.23) can be represented as
5.24
where
5.25
and and *R*^{i} are calculated up to *O*(*k*^{2}) and *O*(*k*), respectively. The *O*(*k*^{2}) expansions of do not depend on and, thus, cannot be simplified further. Full expressions for these terms are given in appendix A.

Integration over the spherical shell *D*^{>} consists of the same two steps as described in §4 for the case of eddy diffusivities: firstly, one performs integration over the horizontal plane *Σ*_{d−1}, and secondly, the result is integrated over the angle *θ* from 0 to *π*. The integration can be carried out analytically, and the result is structurally identical to equation (2.31). The correction to consists of two parts, one of which is proportional to *k*^{2} and the other to . They provide *O*(*ΔΛ*) corrections to *ν* and *δν*_{z} denoted as *Δν* and *Δδν*_{z}, respectively. Upon dividing these corrections by *ΔΛ*, taking the limit *ΔΛ*→0, and substituting *d*=3, one obtains differential equations for *Λ*-dependent *ν* and *δν*_{z},
5.26
and
5.27
Along with (4.26) and (4.27), equations (5.26) and (5.27) form a closed coupled system of ODEs, which can be solved analytically.

For practical use, one needs to relate the forcing amplitude *D* to observable parameters. This can be done via the use of the energy balance equation. In the neutral case, the Green function *G*_{αβ} reverts to the scalar (auxiliary) Green function and, as shown by Sukoriansky *et al*. [6], can be used to calculate the three-dimensional energy spectrum,
5.28
This spectrum can be used in the energy dissipation equation,
5.29
where *ν*_{n} is the effective viscosity in neutral flows. This equation quantifies the energy loss by the slow modes (0<*k*<*Λ*) as acted upon by the fast modes (*k*>*Λ*) represented by the effective viscosity *ν*_{n}. Sukoriansky *et al.* [6] showed that *ν*_{n}≃0.2*D*^{1/3}*Λ*^{−4/3}. Substituting this expression and (5.28) in (5.29), one finds
5.30
which yields *ν*_{n}(*Λ*) in the form of the Richardson law [46,52],
5.31
Clearly, in the neutral case, equations (5.26), (5.27), (4.26) and (4.27) yield *ν*→*ν*_{n}, *δν*_{z}→0, *κ*→*αν*_{n} and *δκ*_{z}→0, where *α* is the inverse turbulent Prandtl number in neutral stratification. Yakhot & Orszag [43] estimated that *α*≃1.4.

It is convenient to use *ν*_{n} for normalization of the solution to the system (5.26), (5.27), (4.26) and (4.27). Because the variables *ν*_{h}, *ν*_{z}, *κ*_{h} and *κ*_{z} are better suited for practical applications, one can use equations (2.19) and (2.20) to make necessary transformation and write the solution in the form
5.32
5.33
5.34
5.35

A general remark needs to be made regarding the forcing in the Langevin equation (2.17). This forcing is assumed isotropic in neutrally stratified flows. Although this assumption may be violated in flows with stable stratification, it is asymptotically valid in the limit of weak stratification considered in the present study. This approximation does not preclude the system from acquiring anisotropic features via anisotropization of the velocity and temperature Green functions. The next order correction can be computed but, even in this approximation, as shown in the following sections, good agreement with data has been achieved.

## 6. Anisotropization of transport properties

Equations (5.32)–(5.35) provide a fully analytical description of the onset of flow anisotropization owing to the action of stable stratification. Derived from nearly first principles, this description provides detailed and important insights into the physics of the small-scale mixing in flows combining turbulence and waves. As emphasized by Ivey *et al.* [53], the requirement of high accuracy is crucial for today's complex numerical models of oceanic and atmospheric circulations. One of the most important theoretical results is that with decreasing , the horizontal effective viscosity and diffusivity increase over their neutral values, while their vertical counterparts exhibit an opposite tendency. In quantitative terms, taking , one finds that all effective viscosities and diffusivities differ from their neutral isotropic values by less than 2 per cent while for , the deviations do not exceed 10 per cent.

Thus far, the wavenumbers *k* and *Λ* denoted different entities. The following sections, however, feature *Λ* only. Reverting to a conventional nomenclature for a running wavenumber, *Λ* will be renamed as *k* from now on.

### (a) Scaling with the Ozmidov wavenumber

The spectral Froude number in equations (5.26), (5.27), (4.26) and (4.27) is defined according to equation (3.1). In flows with weak stratification, this definition could also be based upon *ν*_{n} because the difference between *ν* and *ν*_{n} is small. The latter definition, combined with (5.31), gives
6.1
*k*_{O} being the Ozmidov wavenumber. Similarly to , the ratio *k*/*k*_{O} can be used as a non-dimensional parameter characterizing the strength of stratification. The case of weak stratification corresponds to the limit . Using (5.31), equations (5.32)–(5.35) can be recast in the form
6.2
6.3
6.4
6.5
Note that for identified in §1 as a range of weak stable stratification, the anisotropic corrections to the isotropic values of the effective viscosity and diffusivity do not exceed about 30 per cent.

The normalized effective viscosities and diffusivities are shown in figure 1. This figure clearly demonstrates gradual anisotropization of the main flow characteristics, which leads to the increase of the horizontal and decrease of the vertical effective viscosity and diffusivity compared with their neutral values. The decrease in the vertical direction is much stronger than the increase in the horizontal. The anisotropization of stably stratified turbulence has been well documented in experiments and numerical simulations [35,54–56]. It is further exacerbated in flows that feature horizontal vis-à-vis vertical shear [57–60].

Figure 1 shows that for , the vertical effective viscosity and diffusivity become negative. Recall that equations (6.2)–(6.5) are valid only in weak stratification. As was shown in Sukoriansky *et al.* [6], in the full solution, *ν*_{z} and *κ*_{z} are still decreasing but always remain positive. In fact, in strong stratification, *ν*_{z}/*ν*_{n}, *ν*_{h}/*ν*_{n} and *κ*_{h}/*ν*_{n} attain constant values, while *κ*_{z} becomes scale-independent [6] and can be expressed by the Ellison–Britter–Osborn (EBO) equation [61–64],
6.6
where the mixing efficiency *Γ* is a constant estimated in the QNSE theory at about 0.4 [2]. This value is somewhat larger than the widely accepted in the oceanographic literature estimate of 0.2 [61,65–67]. However, it is well in line with the atmospheric data where *Γ* has been found in the range between 0.33 and 1 [68–72] and DNS [73–75]. Note that *Γ* derived from observations may be affected by factors such as double diffusion processes [76] or lateral stirring by mesoscale fluctuations [77]. Some studies suggest that *Γ* is not a constant, and its value may vary by an order of magnitude [62].

## 7. Turbulence and internal waves

### (a) Dispersion relation for internal waves in presence of turbulence

The Langevin equation for the velocity (2.17) describes a linear, forced, stochastic oscillator whose eigenfrequencies are given by the secular equation
7.1
Straightforward calculation yields a solution to (7.1),
7.2
which describes decaying in time waves. The frequencies of these waves are given by the real part of (7.2),
7.3
or, in terms of the ratio *k*/*k*_{O},
7.4
where is the frequency of the linear internal waves. Equations (7.3) and (7.4) are valid for arbitrary stratification and in fact are dispersion relations for internal waves when they coexist with turbulence. These equations quantify turbulence-caused internal wave frequency shift. In the limit of strong stratification, i.e. when or *k*/*k*_{O}→0, they revert to the classical dispersion relation for linear waves (something akin to the transition from turbulence to ‘wavulence’ [78]). Unlike the latter case, however, in which *ω* depends only on the angle *θ* when *N* is constant, for moderate or *k*/*k*_{O}, *ω* depends on the magnitude of the wavevector as well. This dependence modifies some of the properties of linear waves. For instance, because *ω* is no longer a homogeneous function of the components of **k**, the orthogonality between **k** and the group velocity, ∂*ω*/∂**k**, is lost (see [3], p. 108). In line with these findings, DNSs by Pham *et al.* [74] demonstrate that the linear internal wave theory does not work for waves generated by turbulence.

In contrast to the situation with viscosities and diffusivities, the effect of turbulence on waves is much stronger in weak rather than in strong stratification. Being intuitively clear, this is also evident from the presence of a term proportional to in equation (7.3) rather than its inverse in equations (5.32)–(5.35). Dispersion relation (7.4) shows that in some range of scales, waves coexist with turbulence although waves' frequencies are attenuated by turbulent scrambling. On small scales, the scrambling completely overwhelms orderly fluctuations owing to waves, as the increasing or *k*/*k*_{O} cause the expressions in the curly brackets of the respective equations (7.3) and (7.4) to eventually become negative. The criterion *ω*=0 provides a threshold at which waves vanish. The corresponding critical Froude number, , can be associated with the threshold of internal wave disappearance under the action of turbulence. A similar criterion can be derived in terms of the critical wavenumber *k*_{t},
7.5
Because the waves are attenuated by turbulence at relatively large *k* at which turbulence is close to isotropic, one can use equations (6.2)–(6.5) to find out that *κ*_{z}/*ν*_{n}−*ν*_{z}/*ν*_{n}≃0.4 and *κ*_{h}/*ν*_{n}−*ν*_{h}/*ν*_{n}≃0.4, giving
7.6
The maximum threshold wavenumber is attained along the directions *θ*=±*π*/2, where it is approximately equal to 32*k*_{O}. This threshold is delineated by a vertical line in figure 1. Clearly, the waves are completely suppressed by turbulence on scales where turbulence dominates and where effective viscosities and diffusivities are close to their isotropic neutral values.

The dispersion relation (7.4) normalized with *ω*_{0} is shown in figure 2. This figure clearly demonstrates the dependence *ω*(*k*) and gradual suppression of waves as *k*→*k*_{t}. Internal waves exist only for the wavevectors sweeping the surface of the cylindrical body shown in figure 2, whose boundary at *ω*=0 is approximated by equation (7.6) quite accurately.

In terms of the spectral Froude number, criterion (7.6) reads 7.7 This equation provides further evidence that waves can survive down to small scales whereby turbulence becomes almost isotropic. As discussed earlier, for , all effective viscosities and diffusivities deviate from their neutral isotropic values by no more than 2 per cent.

These results have important theoretical implications, as they establish a quantitative hedge between turbulence and waves. The difference between these two fluid dynamics phenomena is usually smeared in flows where they coexist. Jacobitz *et al.* [79] attempted to establish an objective criterion differentiating between turbulence and waves in stably stratified flows. Their conclusion was that, ‘in agreement with Stewart [80], “there is probably no really clear-cut distinction between turbulence and waves”, at least when both are coexisting at the same scales’ [79], p. 11. A mathematically sharp delineation between these two entities can be stipulated by the fact that only waves have a dispersion relation. In flows on a *β*-plane, for example, turbulence and waves coexist on virtually all scales, which is evidenced by broadening of Rossby waves' spectral peaks with increasing wavenumbers [2]. One has yet to investigate whether or not turbulence causes degeneration of the Rossby wave dispersion relation. In stably stratified flows, such degeneration does take place, thereby carving a spectral domain populated by near-Kolmogorov turbulence with no waves. As will be shown in §9, the distinction between turbulence and waves can also be established based upon their effects on transport processes.

### (b) The buoyancy Reynolds number

An important characteristic of stably stratified turbulence is the buoyancy Reynolds number measuring the distance between the Kolmogorov and Ozmidov wavenumbers,
7.8
Sometimes, *Re*_{b} is viewed as a measure of turbulence activity [53,75,81–83]. It follows from (7.6) that for *k*_{d}≃*k*_{t}, *Re*_{b}≃100. This is an interesting result showing that if the waves impinge upon the viscous dissipation processes, i.e. the inertial subrange of isotropic turbulence is not broad enough, the anisotropization may impact scales proximate to the dissipation range. This result is consistent with some observational studies [84,85] where anisotropization of velocity fluctuations on the scales of dissipation was detected in flows with . DNSs by Shih *et al.* [86,83] are also in line with these observations. These DNSs were used to identify three flow regimes differing by their turbulence intensity: *Re*_{b}<7 was associated with a regime dominated by the molecular processes; a regime with 7<*Re*_{b}<100 was called transitional; and a regime with *Re*_{b}>100 was referred to as energetic (see also [53]). The QNSE theory adds another dimension to this classification by showing that the physics governing the hierarchy of these regimes is intimately related to turbulence–internal wave interaction.

### (c) The absence of the critical Richardson number

The dispersion relation (7.4) indicates that the internal waves populating large scales (*k*/*k*_{O}≪1) become practically indistinguishable from the linear waves. The flow itself, however, may remain far from laminar if judged by the horizontal (isopycnal) and vertical (diapycnal) transfers of momentum and heat, which, as discussed in §6, are strongly anisotropic. In addition, different asymptotic behaviours of *ν*_{z} and *κ*_{z} reflect, on the one hand, the property of the waves to mix momentum but not a scalar, and reveal, on the other, that the turbulent Prandtl number, *Pr*_{t}=*ν*_{z}/*κ*_{z}, is a growing function of the ratio *k*/*k*_{O}, which measures the strength of stratification. As shown in Sukoriansky *et al.* [6,7] and Galperin *et al.* [87], this result provides evidence of the absence of the critical gradient Richardson number, *Ri*_{cr}, in strong stratification (the gradient Richardson number is defined as *Ri*=*N*^{2}/*S*^{2}, *S* being the mean shear, and *Ri*_{cr} is associated with complete suppression of turbulence, i.e. flow laminarization). The conclusion about the absence of *Ri*_{cr} agrees with abundant data. To name only a few examples, Izakov [88], p. 356 mentions that ‘direct measurements of turbulent pulsations of temperature and wind speed …show that turbulent pulsations are always present in the troposphere and lower stratosphere…’. By analysing turbulence variability in very stable atmospheric boundary layers, Mahrt [89], p. 16 concludes that ‘the turbulence never vanishes at large Richardson number partly because of shear generation by the submeso motions’. Mack & Schoeberlein [90] were unable to identify *Ri*_{cr} in the oceanic data. In simulations by Fritts *et al.* [91], pp. 4–5, it was found that ‘…an initial is not a reliable predictor of the absence of turbulence in flows exhibiting significant scale interactions’. The presence of *Ri*_{cr} has been found restrictive in many Reynolds stress models (RSM), leading to insufficient mixing in predicted flows [87]. To alleviate this shortcoming, many recent RSM were designed *Ri*_{cr}-free [92–96]. This new development in modelling of stably stratified turbulent flows was highlighted in the recent overviews by Fernando & Weil [97] and McFarlane [98].

The absence of *Ri*_{cr} is intrinsic to the QNSE theory while in RSM this feature needs to be introduced ad hoc in order to comply with data.

Evidence of preservation of the turbulent nature of a flow in a strongly stratified environment can be inferred from studies of the potential vorticity-connected non-propagating mode that renders stratified turbulence radically different from wave turbulence [32]. This mode is related to the toroidal mode in the Craya–Herring decomposition, and retains its three-dimensionality under the action of any stratification. As will be shown in the next section, flow three-dimensionality is also evident in the QNSE-derived one-dimensional spectra.

It can be argued that since *Ri* does not include molecular viscosity, it cannot characterize the importance of the molecular processes. An appropriate parameter for this purpose would be the buoyancy Reynolds number, *Re*_{b}. These two parameters are independent of each other and, if used in tandem, could provide more complete characterization of different aspects of stably stratified turbulence. For instance, recent investigations by Billant & Chomaz [99] and Brethouwer *et al.* [100] indicate that turbulent mixing is never completely quashed, even under very strong stratification, for as long as *Re*_{b}≫1. Along with this, one needs to keep in mind that *Re*_{b} is only a measure of the width of the Kolmogorov inertial range and not a characteristic of flow laminarization. In fact, energy transfer to horizontal motions and flow layering continue even when *Re*_{b} is very small [101] such that forced stably stratified turbulence may never look like a laminar flow.

## 8. Energy and temperature spectra

### (a) Vertical spectrum of the kinetic energy

In anisotropic problems, the traditional three-dimensional kinetic energy spectrum provides limited information on flow energetics. More detailed information can be gained from various one-dimensional spectra. Among those, the vertical spectrum of the horizontal kinetic energy is one of the most important for oceanic and atmospheric applications. Following the nomenclature in Monin & Yaglom [102], we define by *E*_{1}(*k*_{1}) a one-dimensional longitudinal spectrum of the velocity component *u*_{1}, while *E*_{1}(*k*_{2}) and *E*_{1}(*k*_{3}) are the transverse spectra of *u*_{1} in the horizontal and vertical directions. In the QNSE theory, all spectra can be calculated analytically in the limit of weak stratification [2,6], e.g.
8.1
where *U*_{11}(*ω*,**k**) is the horizontal component of the spectral tensor given by equation (2.29), and *C*_{B}=0.214. Equation (8.1) quantifies the pioneering measurements by Shur [16] in the atmospheric boundary layer that pointed to the transition from the to slope in the kinetic energy spectrum. Shur's work stimulated the rise of early theories of stably stratified turbulence [19–23,103]. Some authors attribute the range with the spectrum to the Kelvin–Helmholtz instability and refer to it as ‘IGW turbulence’ [104]. Other studies associate the spectral transition with the horizontal layering of the flow and formation of the slow manifold [105–107] that take place on scales of the order of the Ozmidov scale. The numerical value of the coefficient, *C*_{B}=0.214, is in agreement with DNS [108] (*C*_{B}≃0.1), large eddy simulations [109] (*C*_{B}≃0.2) and laboratory experiment [110] (*C*_{B}≃0.2).

Equation (8.1) can be used to derive an expression for the spectrum of the vertical shear, , in the universal or canonical form used in physical oceanography [111,112],
8.2
where *x*≡*k*_{3}/*k*_{O}, and *F*(*x*) is a universal function. The universal scaling (8.2) was first proposed by Gargett *et al.* [111] based upon observations in the upper ocean. Later, this scaling was confirmed by Gregg *et al.* [12] using PATCHEX and PATCHEX north datasets. A similar scaling arises from the recent observations by Ledwell *et al.* [113] in the Antarctic Circumpolar Current west of Drake Passage. The QNSE theory gives the canonical spectrum a solid theoretical footing by demonstrating that it can be derived analytically from nearly first principles.

Figure 3 shows remarkable quantitative agreement between the theoretical equation (8.2) and the observational data by Gregg *et al.* [12] for *k*_{3}/*k*_{O}<1. To make a meaningful comparison for *k*_{3}/*k*_{O}>1, one needs to take account of the regime classification by Shih *et al.* [83] (see also [53]). The theory agrees very well with the data from PATCHEX north for which *Re*_{b}≃125 and *k*_{d}/*k*_{O}≃37, i.e. the flow is in the energetic regime. For the PATCHEX data, on the other hand, *Re*_{b}≃14 and *k*_{d}/*k*_{O}≃7, i.e. the Kolmogorov inertial subrange is practically absent, and no agreement with the theory for *k*>*k*_{O} can be expected.

The spectrum (8.1) has been observed not only in the ocean but also in the atmosphere. Vanzandt [114] and Smith *et al.* [9] showed that the observed vertical spectra of the horizontal velocity preserve universal form for all seasons, meteorological conditions and geographical locations throughout the troposphere, stratosphere, thermosphere and mesosphere. This form has been referred to as a canonical gravity wave spectrum [115]. In many studies, this spectrum has been identified with the saturated spectrum of IGWs and its exponent obtained based upon parametrization of breaking waves and dimensional considerations (see [115–125], and references therein). An often used estimate of *C*_{B} throughout the atmosphere is [70,115,126,127], which is close to the QNSE-derived value of 0.2.

Other studies related the canonical spectrum to the buoyancy range turbulence [128,129]. Studies by Smith *et al.* [9], Weinstock [130], Hines [119], Gardner [131] and Dewan [132] emphasized the importance of nonlinear interactions for obtaining a universal power law spectrum. Weinstock [130] discussed the connection between the saturation amplitude of each wave and the collective dissipation caused by all other waves. In the QNSE theory, turbulence and internal waves are considered as one entity in strongly nonlinear, anisotropic regime that gives rise to universal scalings.

### (b) Energy anisotropization

Even though the energy transfers are not calculated in the QNSE theory explicitly, modifications of various one-dimensional spectra provide characterization of energy redistribution between different velocity components and different directions in Fourier space. Using a definition similar to (8.1) and the approximation of weak stratification, one finds that *E*_{1}(*k*_{1}), *E*_{1}(*k*_{2}), *E*_{2}(*k*_{1}) and *E*_{2}(*k*_{2}) are weakly affected by stratification and remain near-Kolmogorovian. This result is in agreement with some recent observational and numerical evidence [13–15].

Consider now the energy of the vertical fluctuations as a function of *k*_{1} and *k*_{3},
8.3
and
8.4

These equations demonstrate that the energy of the vertical fluctuations has a tendency to degenerate with increasing *N* and/or increasing scales. This tendency corresponds to the formation of a two-component flow as the third, vertical, component collapses. Godeferd & Cambon [25] and Sagaut & Cambon [32] provided a thorough investigation of this behaviour.

It is instructive now to compare the QNSE-derived expressions for the one-dimensional spectra of total energy in the directions *k*_{1} and *k*_{3}, i.e. *E*_{tot}(*k*_{1})=*E*_{1}(*k*_{1})+*E*_{2}(*k*_{1})+*E*_{3}(*k*_{1}) and *E*_{tot}(*k*_{3})=*E*_{1}(*k*_{3})+*E*_{2}(*k*_{3})+*E*_{3}(*k*_{3}),
8.5
and
8.6

Clearly, *E*_{tot}(*k*_{1}) decreases and *E*_{tot}(*k*_{3}) increases with increasing strength of stratification pointing to the tendency to energy accumulation in the vertical conical shape in Fourier space. It is consistent with energy accumulation in the vicinity of *k*_{3} for every energy component and with the results obtained from EDQNM by Godeferd & Cambon [25] (fig. 1 in their paper) and by Sagaut & Cambon [32] (fig. 7.5 in their book) who identify it with the toroidal cascade. Combined with the collapse of the energy of the vertical fluctuations, this cascade leads to the formation of the vertically sheared horizontal flows [32,106] and, eventually, vertically layered and decorrelated horizontal flows.

One needs to keep in mind that equations (8.3)–(8.5) are all derived in the approximation of weak stratification and can only be used to illustrate various tendencies in energy redistribution under the action of stable stratification. They should not be extended to scales where their right-hand sides become negative. In the case of strong stratification, full analysis shows that the energy remains positive at all scales and arbitrary stratification [6].

Note that the direction *k*_{3} plays a special role in the dynamics of stably stratified turbulence not only as a result of the energy accumulation in its vicinity but also because of the absence of waves. The horizontal velocity mode corresponding to this direction is known as non-propagating [32] or vortical [33] and can be identified with a slow manifold.

The accumulation of energy in the slow manifold is a typical feature of anisotropic turbulence with dispersive waves. The process of turbulence anisotropization under the action of extra strains that support specific linear waves has been explained within the analysis of triad interactions. Holloway [133,134] applied this analysis to internal waves, while, for example, Rhines [135], Holloway & Hendershott [136], Huang *et al.* [137] and Lee & Smith [138] analysed Rossby waves. Essentially, the necessity to satisfy the frequency resonance condition hampers the efficiency of the triad interactions and reduces their ability to transfer energy across the spectrum. The ensuing anisotropic transfer results in energy accumulation in zero eigenfrequency modes and gives rise to the steepening of the corresponding one-dimensional spectra. Even though the slow manifold is wave-free, this energy redistribution is facilitated by turbulence–wave interaction in other modes and carried out by nonlinear interactions. It is thus suggestive that the characteristic steep, one-dimensional, canonical spectrum of stably stratified turbulence (8.1) is not a spectrum of gravity waves, but rather is a manifestation of the slow manifold.

### (c) Vertical spectrum of the potential energy

Using the temperature Langevin equation (2.8), one can show that when the temperature is a passive scalar, the three-dimensional spectrum of one half of the temperature variance, , obeys the Corrsin–Obukhov law [43],
8.7
where *ϵ*_{θ} is the dissipation rate of , and *C*_{CO}=1.03 is the Corrsin–Obukhov constant. The value of this constant is in fair agreement with laboratory measurements and DNSs that show some spread but point to *C*_{CO}≃1 [139,140].

In the case of weak stratification, equations (6.2)–(6.5) can be used to compute the vertical spectrum of the temperature variance,
8.8
where *C*_{T}=0.62, *C*_{N1}=0.406 and *C*_{N2}=0.134, and the first term on the right-hand side is the one-dimensional Corrsin–Obukhov spectrum [2]. The spectrum *E*_{T}(*k*) is related to the spectrum of the turbulence potential energy, *E*_{p}(*k*), by
8.9
Furthermore, the potential energy dissipation rate, *ϵ*_{p}, is related to *ϵ*_{θ} by *ϵ*_{p}=(*αg*/*N*)^{2}*ϵ*_{θ}. In the case of strong stratification, *ϵ*_{p} and *ϵ* are related by
8.10
where *Γ* is the mixing efficiency (see equation (6.6)) [2,74,75]. Using these relationships, equation (8.8) can be re-written in terms of the potential energy spectrum,
8.11
where *γ*_{3}=*C*_{T}*C*_{N1}*Γ*+*C*_{N2}≃0.23. Similar to the vertical spectrum of the horizontal kinetic energy (8.1), the spectrum (8.11) combines the turbulent, Corrsin–Obukhov part and the contribution due to stable stratification. Remarkably, the latter contributions in the horizontal components of the kinetic energy and in the potential energy turn out to be approximately the same, . Numerical simulations by Bouruet-Aubertot *et al.* [108] also gave *γ*_{3}=0.2. The energy equipartition between the wave parts of the kinetic and potential energies has been discussed by Staquet [104]. Various expressions for the temperature spectra were developed previously [22,141] although they did not account for flow anisotropy.

Stratospheric observations generally confirm the scaling (8.11) extrapolated to the case of strong stratification,
8.12
with the proportionality coefficient between 0.1 and 0.3 [142–146]. Cot [147] confirmed the proportionality of the vertical spectra of the horizontal kinetic and potential energies in the troposphere and the stratosphere. Noting observational difficulties, he nevertheless estimates the ratio of these spectra at about 2.5–3. The QNSE equations (8.1) and (8.12) value this ratio at 2 (note that *E*_{1}(*k*_{3}) needs to be doubled to account for the contributions from both horizontal kinetic energy components), in fair agreement with other estimates.

### (d) Horizontal spectrum of the potential energy

In the case of weak stratification, equations (6.2)–(6.5) can be used to derive an analytical expression for the horizontal spectrum of the temperature variance,
8.13
where *C*_{Θ}=0.87, *C*_{M1}=0.097 and *C*_{M2}=0.535. This expression can be re-written in a form similar to equation (8.8) but for the horizontal spectrum of the temperature variance,
8.14
or the horizontal spectrum of the potential energy,
8.15
where *γ*_{1}=*C*_{Θ}*C*_{M1}(*ϵ*_{p}/*ϵ*)+*C*_{M2}. The first addend in this expression is small compared with the second one, and can be neglected, leaving *γ*_{1}≃*C*_{M2}. Similar to the vertical spectrum of the potential energy, its horizontal counterpart is also composed of the Corrsin–Obukhov part and the stratification contribution that depends on *N* and *k*_{h} only and dominates on relatively large scales.

The existence of the spectrum (8.15) is not indisputable. A study by Klymak & Moum [148] suggests that the inertial subrange with the Corrsin–Obukhov spectrum extends to horizontal wavenumbers *k*_{h} smaller than *k*_{O} by some one or two orders of magnitude. Further discussion along these lines can be found in Riley & Lindborg [15].

Bouruet-Aubertot *et al.* [10] provided characterization of temperature variance spectra using high-resolution temperature measurements in deep water made during the Rockall Channel Studies. Based upon these observations, the temperature variance spectra were computed as functions of the horizontal wavenumber, *k*_{h}. For large values of *k*_{h}, these spectra exhibited robust Corrsin–Obukhov regime with the characteristic *k*^{−5/3}_{h} spectral slope. For relatively small *k*_{h}, however, the spectra underwent steepening and seemed to approach the slope. Figure 4 compares the low wavenumber ranges of the measured spectra with the wave-dominated subrange of the theoretical expression (8.14). These comparisons are of particular interest because the theoretical expression depends on d*Θ*/d*z* only, a parameter that can be estimated with much better accuracy than the dissipation rates of the kinetic and potential energies, *ϵ* and *ϵ*_{θ}, which were not measured independently. Owing to this uncertainty, comparisons in the turbulence-dominated subrange were not sought. In the buoyancy subrange, however, the agreement between the data and the theoretical expression is quite good not only with regard to the spectral slope but also for the amplitude. This agreement is significant and can be viewed as confirming the functional dependence in (8.14), including the value of the numerical coefficient *γ*_{1}.

## 9. Transport of momentum and scalar

What are the implications of the changes introduced by stable stratification in spectral characteristics of turbulence for the transport of momentum and scalar? One could rephrase this question to align it with the issues raised by Weinstock [130], i.e. is there a connection between the effective diffusivity and dissipation and the gravity wave spectrum? And its reverse, what can be inferred from the observed spectra for the diffusivity? Assuming that ‘the diffusivity depends on the scale being diffused’, Weinstock [130] introduced a formal cutoff at the diffusive scale in the integral that relates the diffusivity to the spectrum in his model. Ultimately, this approach yields a scale-dependent diffusivity. For the experimental determination of the diffusivity coefficient, however, Weinstock [130] reverted to the scale-independent EBO expression similar to (6.6) in which the dissipation rate, *ϵ*, is measured.

It is well known that linear internal waves make no contribution to diapycnal scalar transport [4], p. 166. It is important to verify that wave-dominated large-scale modes comply with the requirement of the diminishing diapycnal transport in the QNSE theory. In the limit of weak stratification, equations (6.2)–(6.5) predict that the effective viscosities and diffusivities preserve their scale-dependent, Richardson-law-like structure. The limit of strong stratification requires the invocation of the full QNSE solution. In §6*a*, it was shown that in this limit, the QNSE asymptotic expression for the diapycnal diffusivity coincides with the EBO expression (6.6). This expression can be rearranged in a form congruent with the Richardson law,
9.1
Mirroring the classical Richardson diffusion law, the QNSE theory implies that the effective viscosity and diffusivity at a scale *k*^{−1} result from contributions by all scales smaller than *k*^{−1}. Within this framework, equation (9.1) has the following interpretation. When the diffusion scale *k*^{−1} is in turbulence-dominated range extending from up to a fraction of *k*^{−1}_{O}, all scales up to *k*^{−1} contribute to the diapycnal diffusivity. Larger, wave-dominated scales make no contribution to *κ*_{z} and so it becomes scale-independent being sustained by the turbulent scales only [1]. This result demonstrates that *κ*_{z} has correct asymptotic behaviour in the limit of strong stable stratification and provides a new interpretation of the EBO model. It also clarifies the model by Weinstock [130] by showing that the scale-dependency of *κ*_{z} should only be confined to the turbulence subrange.

The scale-independence of *κ*_{z} in strongly stratified flows is an important result that can be verified against data. Studies of the mixing in the oceanic pycnocline by Ledwell *et al.* [149] indicate that the diapycnal diffusivity remains approximately constant and largely scale-independent on the time scales from six to 24 months even though the vertical spread of the tracer during that time was considerable. From the physical viewpoint, because the vertical diffusion is carried out within turbulent patches dominated by three-dimensional overturning turbulence, the vertical scale of these patches measured by *k*^{−1}_{O} is a pertinent scale for the eddy diffusivity (rather than the vertical spread of a tracer). The calculation of the evolution of an initial tracer concentration profile to a vertical profile after 1-year integration using a one-dimensional diffusion equation with a constant, scale-independent eddy diffusivity estimated according to the EBO model was in good agreement with the measurements reported by Ledwell *et al*. [113], fig. 2. This result supports the validity of the EBO model and the physical reasoning behind it.

In §7*c*, it was argued that the absence of the critical Richardson number reflects the fact that in flows with strong stratification, the diapycnal momentum flux, being carried out by waves, remains finite while its scalar counterpart diminishes. Fernando & Weil [97] hypothesize that such behaviours of diapycnal momentum and scalar fluxes can be used to develop a demarcation line between the scales dominated by either turbulence or waves. The EBO model (9.1) provides even sharper differentiation between turbulence- and wave-dominated scales. Recall that another criterion to separate turbulence and waves subranges was considered in §7*a*. That criterion referred to the scale at which the internal wave dispersion relation degenerates and the waves disappear completely.

## 10. Discussion and conclusions

The main result of this paper is a derivation of an analytical theory of transition from the classical Kolmogorov to buoyancy-dominated, stratified turbulence. This result is achieved in the framework of the QNSE theory, which is maximally proximate to first principles. Considering turbulence and waves as one entity, the theory unifies spectral characteristics, effective diapycnal and isopycnal viscosities and diffusivities, rates of dissipation of the kinetic and potential energies and flow anisotropization. A theoretical framework capable of comprehensive analysis of turbulence and internal waves has been sought since the inception of early theories of stably stratified turbulence in late 1950s–early 1960s; the need for such a theory was reiterated by Weinstock [130].

The theory quantifies both the decrease of vertical viscosities and diffusivities and the increase of their horizontal counterparts compared with the neutral case. In addition, it details how turbulence affects IGWs, quantifies waves' frequency attenuation and demonstrates that turbulence can cause degeneration of the internal wave dispersion relation. It is of importance that this degeneration coincides with the scales at which turbulence begins to undergo anisotropization.

Studies of one-dimensional vertical spectra of kinetic and potential energies reveal that they are composed, respectively, of the sum of the Kolmogorov or Corrsin–Obukhov parts and parts generated by contributions from stable stratification, the latter distinguished by the characteristic dependence on *N*^{2} and spectral exponents with the −3 slope. For the kinetic energy, these steep slopes are attributed not to the direct effect by IGWs but to the formation of a slow manifold facilitated by the frequency resonances of wavevector triads. The horizontal spectrum of temperature fluctuations appears to have a similar structure. Such two-term functional representation enables one to extend the universal scaling first proposed by Gargett *et al.* [111] for the vertical shear spectrum to the vertical and horizontal spectra of the potential energy. It is of importance that the levels of both kinetic and potential energies in the buoyancy subranges are independent of turbulence characteristics and determined solely by the Brunt–Väisälä frequency. This result may explain the enigmatic energetic universality of the celebrated vertical spectra of the kinetic energy in the atmosphere and of the Garrett and Munk spectrum of the IGWs in the ocean.

It was found that the spectral estimates based upon the exact results obtained in the approximation of weak stratification held on scales where stratification is strong. It is possible, however, that on these scales, some of the spectra may return to their Kolmogorov shape albeit with different coefficients. Also, it is likely that on such scales, the effects of the Earth rotation cannot be neglected. Because the current version of the QNSE theory does not include rotational effects, studies of spectral modifications on large scales are beyond the scope of this paper. They will be addressed in future research.

Along with the spectral estimates, another important outcome of the QNSE theory is the re-evaluation of the transport of momentum and scalars in flows that combine turbulence and waves. Reflecting the fact that internal waves can transfer momentum but not scalars, the effective vertical viscosity remains finite while the vertical diffusivity diminishes on large scales where stratification is strong. On the other hand, both horizontal viscosity and diffusivity increase on large scales compared with the case of neutral stratification. Such behaviour of transport coefficients is quite different from the neutral case and points to: (i) growing vertical Prandtl number, *Pr*_{t}, with stratification strength; (ii) the failure of the Reynolds analogy; (iii) anisotropization rather than degeneration of turbulent transport and ensuing unlikeliness of flow laminarization in arbitrarily strong stratification; and (iv) the absence of the critical Richardson number.

Another outcome of the QNSE theory is a new interpretation of the EBO diapycnal mixing model according to which only turbulent scales contribute to diapycnal diffusivity in flows that combine turbulence and waves. Applied in the oceanographic context, this result stipulates that the breaking internal waves contribute to the diapycnal scalar mixing by generating turbulence dominated scales of the order of *O*(*k*^{−1}_{O}). In agreement with the conclusions reached in Ledwell *et al.* [113], the prevailing in the mid-depth oceanic interior values of *N* and *ϵ* yield the diapycnal diffusivity that appears to be an order of magnitude smaller than what is necessary to maintain the overturning circulation.

Concluding, we re-emphasize that the QNSE theory allows one to calculate anisotropic, scale-dependent viscosities and diffusivities and, thus, the response functions. In the EDQNM framework, one can calculate turbulent fluxes while making certain approximations regarding the response functions. Both theories are mutually complementary and, as was shown earlier, produce complementary results (e.g. collapse of the vertical fluctuations and energy accumulation in the slow, non-propagating mode). One can expect that future progress can be achieved by combining these approaches. This direction has already been explored by Dannevik *et al.* [150] for isotropic three-dimensional case and by Sukoriansky *et al.* [151,152] for two-dimensional turbulence with inverse cascade and sign-changing renormalized viscosity.

## Acknowledgements

Some of the materials of this paper were reported at the Second International Conference on Turbulence Mixing and Beyond and we are grateful to the conference organizers and in particular Dr Snezhana Abarzhi for giving us the opportunity to present these results. We also thank Dr Pascale Bouruet-Aubertot for providing us with the data from studies of the bottom boundary layer in the Rockall Channel. Many thanks also go to an anonymous reviewer who provided a thorough report that helped us to clarify and improve the presentation. Partial support of this research by the ARO grant nos W911NF-05-1-0055 and W911NF-09-1-0018, ONR grant nos N00014-07-1-1065 and N00014-10-1-0519, and the Israel Science Foundation grant no. 134/03 is gratefully acknowledged.

## Appendix A. Expressions for in equation (5.25)

A 1 A 2 A 3

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales I’.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.