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.
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.  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=N2/S2, 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/kO, where k is a wavenumber, kO=(N3/ϵ)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/kO.
Mahrt et al.  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  diagnosed layers where the vertical spectra of the horizontal velocity and temperature fluctuations scaled as either or , kz 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 , 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.  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−1O, 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 . These observations stimulated the development of early theories of stratified turbulence [17–23]. A review of some later work was given by Galperin & Sukoriansky  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 . 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  for flows with stable stratification. For rotating and stratified flows, this approach was further developed in Cambon & Scott , Cambon [28,29], Godeferd & Staquet , Cambon et al.  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 . 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 . 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 . 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  and Sagaut & Cambon .
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 , 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) . 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Θ/dz))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 Λ, Λ=kd.
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. , 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 νk2 and κk2, 2.19 and 2.20 where ; kx=k1; ky=k2; ; and kz=k3 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 k3, i.e. the vertical coordinate, k3, 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≃Λ . 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(k2) 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 , who showed that it provides an accurate description of the effective viscosity in isotropic three-dimensional turbulence. In addition, Smith & Woodruff  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 , 2.26 and 2.27 where 2.28 and Uαβ(Ω,q,q3) 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 Λ/kO 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 (Λ/kO)−4/3. Note that the early theories of stably stratified turbulence mentioned in §1 dealt with O[(k/kO)−4/3] corrections to the Kolmogorov theory, i.e. de facto they considered the limit of weak stratification.
4. Computation of eddy diffusivities
— 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(k2). The factor kαkβ in front of the integral is already O(k2), and so only the O(1) terms should be retained in the integrand (4.1). Setting k=0 and k3=0 in GT 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 S1(Ω,q,q3)
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 R1 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 S2(Ω,q,q3)
This term possesses seven poles, three of which are given by equation (4.9), and the rest, owing to the factor B(Ω,q,q3), 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 R2 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(k2) 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 k3 were set to zero owing to the presence of the O(k2) factor, kαkβ, in front of the integral in (2.27). The contractions of this factor with and are, thus, O(k2) 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 q3 only (these terms are due to the residues R1 and R2) 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 qh 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, qh=(q1,q2,0), and vertical, q3, 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 Sd=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 k2 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|,k3−q3) and B=B(Ω,q,q3).
Because A and B are both small, the term S4 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 .
Similar to the case of the temperature equation (4.4), only the scalar functions Si 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 S1
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 R1(k,k3,q,q3) (for brevity, the arguments of Ri 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 S2
The product of and A produces, in addition to the three poles of S1, 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 R2 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 S3
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, S3 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(k2) 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 Ri 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 Ri are calculated up to O(k2) and O(k), respectively. The O(k2) 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 k2 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. , 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.  showed that νn≃0.2D1/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  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. , 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 kO being the Ozmidov wavenumber. Similarly to , the ratio k/kO 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. , 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  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 . 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  or lateral stirring by mesoscale fluctuations . Some studies suggest that Γ is not a constant, and its value may vary by an order of magnitude .
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/kO, 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/kO→0, they revert to the classical dispersion relation for linear waves (something akin to the transition from turbulence to ‘wavulence’ ). Unlike the latter case, however, in which ω depends only on the angle θ when N is constant, for moderate or k/kO, ω 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 , p. 108). In line with these findings, DNSs by Pham et al.  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/kO 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 kt, 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 32kO. 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→kt. 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.  attempted to establish an objective criterion differentiating between turbulence and waves in stably stratified flows. Their conclusion was that, ‘in agreement with Stewart , “there is probably no really clear-cut distinction between turbulence and waves”, at least when both are coexisting at the same scales’ , 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 . 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, Reb is viewed as a measure of turbulence activity [53,75,81–83]. It follows from (7.6) that for kd≃kt, Reb≃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: Reb<7 was associated with a regime dominated by the molecular processes; a regime with 7<Reb<100 was called transitional; and a regime with Reb>100 was referred to as energetic (see also ). 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/kO≪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, Prt=νz/κz, is a growing function of the ratio k/kO, which measures the strength of stratification. As shown in Sukoriansky et al. [6,7] and Galperin et al. , this result provides evidence of the absence of the critical gradient Richardson number, Ricr, in strong stratification (the gradient Richardson number is defined as Ri=N2/S2, S being the mean shear, and Ricr is associated with complete suppression of turbulence, i.e. flow laminarization). The conclusion about the absence of Ricr agrees with abundant data. To name only a few examples, Izakov , 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 , p. 16 concludes that ‘the turbulence never vanishes at large Richardson number partly because of shear generation by the submeso motions’. Mack & Schoeberlein  were unable to identify Ricr in the oceanic data. In simulations by Fritts et al. , 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 Ricr has been found restrictive in many Reynolds stress models (RSM), leading to insufficient mixing in predicted flows . To alleviate this shortcoming, many recent RSM were designed Ricr-free [92–96]. This new development in modelling of stably stratified turbulent flows was highlighted in the recent overviews by Fernando & Weil  and McFarlane .
The absence of Ricr 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 . 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, Reb. 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  and Brethouwer et al.  indicate that turbulent mixing is never completely quashed, even under very strong stratification, for as long as Reb≫1. Along with this, one needs to keep in mind that Reb 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 Reb is very small  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 , we define by E1(k1) a one-dimensional longitudinal spectrum of the velocity component u1, while E1(k2) and E1(k3) are the transverse spectra of u1 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 U11(ω,k) is the horizontal component of the spectral tensor given by equation (2.29), and CB=0.214. Equation (8.1) quantifies the pioneering measurements by Shur  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’ . 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, CB=0.214, is in agreement with DNS  (CB≃0.1), large eddy simulations  (CB≃0.2) and laboratory experiment  (CB≃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≡k3/kO, and F(x) is a universal function. The universal scaling (8.2) was first proposed by Gargett et al.  based upon observations in the upper ocean. Later, this scaling was confirmed by Gregg et al.  using PATCHEX and PATCHEX north datasets. A similar scaling arises from the recent observations by Ledwell et al.  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.  for k3/kO<1. To make a meaningful comparison for k3/kO>1, one needs to take account of the regime classification by Shih et al.  (see also ). The theory agrees very well with the data from PATCHEX north for which Reb≃125 and kd/kO≃37, i.e. the flow is in the energetic regime. For the PATCHEX data, on the other hand, Reb≃14 and kd/kO≃7, i.e. the Kolmogorov inertial subrange is practically absent, and no agreement with the theory for k>kO can be expected.
The spectrum (8.1) has been observed not only in the ocean but also in the atmosphere. Vanzandt  and Smith et al.  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 . 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 CB 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. , Weinstock , Hines , Gardner  and Dewan  emphasized the importance of nonlinear interactions for obtaining a universal power law spectrum. Weinstock  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 E1(k1), E1(k2), E2(k1) and E2(k2) 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 k1 and k3, 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  and Sagaut & Cambon  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 k1 and k3, i.e. Etot(k1)=E1(k1)+E2(k1)+E3(k1) and Etot(k3)=E1(k3)+E2(k3)+E3(k3), 8.5 and 8.6
Clearly, Etot(k1) decreases and Etot(k3) 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 k3 for every energy component and with the results obtained from EDQNM by Godeferd & Cambon  (fig. 1 in their paper) and by Sagaut & Cambon  (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 .
Note that the direction k3 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  or vortical  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 , Holloway & Hendershott , Huang et al.  and Lee & Smith  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 , 8.7 where ϵθ is the dissipation rate of , and CCO=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 CCO≃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 CT=0.62, CN1=0.406 and CN2=0.134, and the first term on the right-hand side is the one-dimensional Corrsin–Obukhov spectrum . The spectrum ET(k) is related to the spectrum of the turbulence potential energy, Ep(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=CTCN1Γ+CN2≃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.  also gave γ3=0.2. The energy equipartition between the wave parts of the kinetic and potential energies has been discussed by Staquet . 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  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 E1(k3) 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, CM1=0.097 and CM2=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ΘCM1(ϵp/ϵ)+CM2. The first addend in this expression is small compared with the second one, and can be neglected, leaving γ1≃CM2. 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 kh only and dominates on relatively large scales.
The existence of the spectrum (8.15) is not indisputable. A study by Klymak & Moum  suggests that the inertial subrange with the Corrsin–Obukhov spectrum extends to horizontal wavenumbers kh smaller than kO by some one or two orders of magnitude. Further discussion along these lines can be found in Riley & Lindborg .
Bouruet-Aubertot et al.  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, kh. For large values of kh, these spectra exhibited robust Corrsin–Obukhov regime with the characteristic k−5/3h spectral slope. For relatively small kh, 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Θ/dz 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 , 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  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  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 , 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 §6a, 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−1O, 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 . 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  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.  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−1O 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. , fig. 2. This result supports the validity of the EBO model and the physical reasoning behind it.
In §7c, 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  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 §7a. 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 .
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 N2 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.  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, Prt, 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−1O). In agreement with the conclusions reached in Ledwell et al. , 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.  for isotropic three-dimensional case and by Sukoriansky et al. [151,152] for two-dimensional turbulence with inverse cascade and sign-changing renormalized viscosity.
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
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.