Small-scale anisotropy in stably stratified turbulence

In this paper, small-scale statistics in stably stratified turbulence is studied theoretically and numerically. Expressions for the spectra of the velocity correlation, density fluctuation and buoyancy flux are derived using a perturbation method that generalizes the idea used by Ishihara et al for studying small-scale anisotropy in turbulent shear flows (2002 Phys. Rev. Lett. 88 154501). The expressions give not only the scaling in the inertial subrange but also estimates for anisotropy of the spectra at small scales. They were found to be in good agreement with a direct numerical simulation with 5123 grid points.

A turbulence system consists of a large number of degrees of freedom. The motion of each degree of freedom as well as its trajectory in phase space is sensitive to and affected by flow conditions, such as the initial and boundary conditions. Despite differences in flows, it is widely believed that, under certain conditions, the statistics of turbulence is universal, in the sense that it is insensitive to details of the flow conditions. It is also expected that the flow statistics can be described in terms of a few parameters. These expectations have been basically supported by experiments and direct numerical simulations (DNSs) of turbulence. Kolmogorov's hypotheses (K41) [1] are consistent with this idea of universality. According to K41, in fully developed turbulence at sufficiently high Reynolds numbers Re and at scales much smaller than the characteristic length scales L of the energy-containing eddies, there exists the so-called universal equilibrium range, where the statistics is insensitive to details of the flow conditions and can be described in terms of a few parameters, such as the energy dissipation rate per unit mass, the viscosity ν and the length scale .
This picture of turbulence is similar to that known in thermodynamics and statistical mechanics for a thermal equilibrium system, which also contains a large number of degrees of freedom. Such a system exhibits universal macroscopic relations that are insensitive to details of the differences between phase-space trajectories at the microscopic level, and can be described in terms of a few macroscopic state variables, such as pressure, mass density and temperature. In studies of thermal equilibrium or near-equilibrium systems, at least two classes of universal relations are known: (i) those characterizing the equilibrium state itself, such as Boyle-Charles' law and (ii) those characterizing the response to disturbances added to the equilibrium state (see e.g. [2]). As an example of the relation in (ii), consider a Newtonian incompressible fluid at rest in a thermal equilibrium state. If a disturbance or an external force is added to the fluid, the fluid is set into macroscopic motions that obey the following linear relation between the rate of strain tensor S ij and the viscous stress (momentum flux) tensor τ ij : where C ijmn is a fourth-order isotropic tensor that reflects the thermal equilibrium state of the fluid.
Even if the idea of universality, such as Kolmogorov's hypotheses, is asymptotically correct for Re → ∞, L/ → ∞, in real turbulence, Re and L/ can only be finite, so that the turbulence is not free from the effects (which we call here disturbances), which may cause deviations from the universal equilibrium state. For example, even if the statistics in the universal equilibrium range is locally isotropic for Re → ∞, L/ → ∞, in real turbulence with finite Re and L/ , the statistics cannot be isotropic in a strict sense owing to the non-universal disturbances. This fact, together with the analogy of turbulence with thermal equilibrium systems, encouraged us to examine the possible universality of turbulence in the sense of (ii) mentioned above.
One of the various possible disturbances that may cause deviations from the equilibrium state is the existence of a mean flow. Recently, Ishihara et al [3] (hereafter referred to as IYK) studied the effects of a homogeneous mean shear flow at small scales using a simple perturbation analysis. Their analysis suggests that the deviation Q ij (k) of the velocity correlation spectrum from the equilibrium state due to the mean shear can be given by at a large wavenumber k, where S ij is the rate of strain tensor of the mean shear flow, and the integration of Q ij (k) over the wave vector k gives the so-called turbulence Reynolds stress. Similar to equation (1), C ijmn is a fourth-order isotropic tensor that reflects the equilibrium state and is independent of S mn . However, unlike in equation (1), C ijmn in (2) depends on k and (2) holds only for sufficiently large k, i.e. small scales. The form of (2) agrees well with DNS results [3]. The scaling derived from (2) in the inertial subrange [3] is in agreement with the scalings obtained by Lumley [4] and Biferale et al [5] (see section 2). The equilibrium state is modified not only by the existence of the anisotropic and/or inhomogeneous mean flows, but also by external forces, such as buoyancy and magnetic fields. In this paper, we consider stably stratified turbulence as a typical example of turbulence subjected to anisotropic forces or disturbances. In section 2, we review the basic idea of IYK and consider how the equilibrium state responds to disturbances in a rather general framework. We then apply this idea to stably stratified turbulence and derive relations between the correlation spectra and the disturbance (or buoyancy force) in section 3. The results are compared with DNS data in section 4, where particular attention is paid to the angular dependence of the correlation spectra.

Small-scale statistics in turbulent shear flows
In this section, we consider an incompressible turbulent velocity field u(x, t) that obeys the Navier-Stokes equations where ν is the kinematic viscosity, p the pressure, F the external force and the density is assumed to be unity. In fully developed turbulence at high Reynolds numbers, the time dependence of u in the Eulerian frame of reference is dominated by the so-called sweeping effect of large eddies. This effect is not to be confused with the non-linear interactions between small eddies of similar size. To avoid confusion, it is convenient to introduce a co-ordinate system moving with the velocity of a fluid particle at a certain time, as in K41. Let v be the velocity in the moving frame of reference, defined as v(r, t) ≡ u(r + x 0 + u 0 t, t + t 0 ) − u 0 with u 0 ≡ u(x 0 , t 0 ), and let us consider the statistics in a local domain D = {(r, t) | r L and |t| T }, where L and T are the characteristic length and time scales of the energy-containing eddies, respectively, or those of the mean flow. In terms of v, (3) can be written as where G = G(r, t) = F (r + x 0 + u 0 t, t + t 0 ). If we decompose the field into the mean and fluctuating parts, such as v = v +ṽ with ṽ = 0, then (5) gives where E ≡ −M +G, and the brackets denote a statistical mean or average taken over an appropriate ensemble. Equation (6) shows that the dynamics ofṽ are governed by three terms: (a) the convective term (ṽ · ∇)ṽ, (b) the viscous term ν∇ 2ṽ and (c) the term E representing effects of the mean flow and external forces. Here, we ignore the pressure term, since it can be eliminated by the use of the incompressibility condition (4).
To obtain some idea of the relative magnitude of these terms, we assume here for the sake of simplicity that (i) the external force is negligible (G = 0), i.e. E = −M, (ii) the mean shear rate S ij ≡ ∂ v i /∂r j is almost constant in the local domain D and (iii) the statistics are almost homogeneous in the domain D. The mean field v may then be expanded for small r as v i = S ij r j + · · · and where the higher-order terms in r have been omitted. Hereinafter, we use the summation convention for repeated indices.
Let v be the characteristic velocity of eddies of size in D or, equivalently, the characteristic velocity |v| at r ∼ , where the symbol ∼ denotes equality to an order of magnitude. Then, dimensional considerations suggest that the order of magnitude of the three terms in (6) associated with the motions of these eddies are given by where S is an appropriate norm of the tensor S ij . Equation (8) gives M/[(ṽ · ∇)ṽ] ∼ S /v , so that |M| will be much smaller than |(ṽ · ∇)ṽ| provided that δ ≡ S /v 1. According to K41, v obeys the scaling v ∼ 1/3 1/3 in the inertial subrange L η, where η is the dissipation length scale at which v 2 η /η ∼ νv η /η 2 , i.e. Re ≡ v /ν ∼ 1. Then (8) yields M/[(ṽ · ∇)ṽ] ∼ δ ∼ S 2/3 / 1/3 . This implies that the effect of M decreases with the size faster when compared with the effect of the convective term (ṽ · ∇)ṽ and there may exist a range of ( E ≡ ( 1/3 /S) 3/2 , in the present case) for which the magnitude of E is much smaller when compared with that of the nonlinear convective term. A similar result may hold true for turbulence under various kinds of external forces, such as stably stratified turbulence (see section 3), rotating turbulence, magnetohydrodynamic turbulence subjected to a uniform magnetic field and β-plane turbulence in two dimensions. The ideas discussed in this section are also applicable to these types of turbulence.
In terms of the characteristic time scales τ N , τ ν and τ E associated with the terms (a), (b) and (c), respectively, (8) can be written as To take into account the incompressibility condition, it is convenient to introduce wave vector space. Let Q ij be the velocity correlation spectrum tensor, defined as Recalling that the Navier-Stokes dynamics with S ij = 0 is compatible with the local homogeneities and isotropies at small scales L, i.e. in the wavenumber range k ≡ |k| k L ≡ 1/L, we write Q as where Q (0) ij represents the spectrum in the absence of mean shear and external forces, and is isotropic in the high wavenumber range k k L , The above discussion leading to (8) and (9) gives This implies that in the high wavenumber range for which δ(k) 1, the effects of the term M are weak when compared with those of the convective term, so that it may be treated as a perturbation when estimating Q ij (k) on the basis of (3) or, equivalently, (6). Since E in (6) is linear in S ij (see (7)), it is tempting to assume that Q is of the form Here, from symmetry considerations, it is shown that C ijmn can be written without loss of generality as for almost stationary turbulence, where a and b are isotropic functions of k, i.e. they depend on k only through k [3]. This form has been confirmed to be in good agreement with DNS results [3]. Quantitative estimates of a(k) and b(k) have been obtained both in DNSs [3,6] and experiments [7]; these are also in good agreement. If we assume that K41 holds, then [4], who discussed the influence of viscosity, buoyancy, magnetic fields and elasticity on the scaling, and also with Biferale et al [5], who presented a prediction that generalized the scaling for all orders of structure functions and all SO(3) sectors. Note that (i) equation (11) with (12) applies not only for the similarity range, but is valid as long as δ(k) 1, and (ii) equation (11) gives not only the scaling in the inertial subrange, but also gives the anisotropy (i.e. the dependence on the components i, j and on the direction of k) of Q ij (k) in terms of only two scalar functions a(k) and b(k), unlike the expressions in Lumley [4] and Biferale et al [5]. The form (11) with (12) is also consistent with the Lagrangian renormalized approximation (LRA) [8]; a theoretical estimate of a(k) and b(k) has been obtained on the basis of the LRA [6].

Stably stratified turbulence
We consider here a stably stratified turbulence velocity field that obeys the Boussinesq approximation. It can be written as after appropriate normalizations, where u, ρ and p are the fluctuating velocity, density and pressure with zero mean, respectively; ν and κ are the kinematic viscosity and molecular diffusivity, respectively; e 3 is the unit vector in the x 3 -direction, which is anti-parallel to gravity and N is the Brunt-Väisälä frequency given by Here, g is the acceleration due to gravity, ρ 0 the reference density andρ(x 3 ) a stably stratified background-density profile with dρ/dx 3 (<0) being constant. This system contains not only the velocity and pressure fields but also the fluctuating density field. It is not difficult to generalize the ideas of the previous section to include the effect of the density field. Let us define the density spectrum P(k, t) and the buoyancy-flux spectrum tensor B i (k, t) as follows: Hereafter, we consider statistically quasi-stationary states and omit writing time t in the spectrum tensors.
Corresponding to (10), let us write where X (0) denotes the equilibrium spectrum in the absence of buoyancy, i.e. N = 0, and X denotes the deviation from X (0) . We assume here that the equilibrium spectra for k k L are isotropic, which implies that Let v and ρ be the characteristic velocity and amplitude of the density fluctuations, respectively, associated with eddies of size ∼ 1/k. Then, similar to (9), we have where τ N is the characteristic time scale associated with the nonlinear convective terms (but in the moving frame of reference as in (5)) and τ E is the linear interaction between u and ρ (the N terms) in (13) and (15). In wave vector space, we have τ N /τ E ∼ δ(k) ≡ N/[kv(k)] for eddies of size 1/k, so that in the high wavenumber range, for which δ(k) 1, the effect of the buoyancy force is weak when compared with that of the convective term and can be treated as a perturbation to the dynamics (13)- (15) with N = 0. Then, by following the idea leading to (11), we try to write expansions for the deviation spectra in powers of δ(k) or N = Ne 3 , e.g.

Comparison with DNS
To examine the theoretical conjectures obtained in the previous section, we performed a DNS of decaying stratified turbulence that obeys (13)- (15). The fluctuating velocity u, density ρ and pressure p were assumed to be periodic with a period 2π in each of the three Cartesian directions. The background densityρ appears in (13)-(15) only through the constant N as given in (16). The DNS was based on an alias-free spectral method with 512 3 grid points and a maximum wavenumber of k max = 241. A fourth-order Runge-Kutta method is used for the time evolution. We set N = 1 and ν = κ = 2.89 × 10 −4 . There have been extensive studies of stably stratified turbulence using numerical simulations (see e.g. [9]- [17]). In the present DNS, (i) the number of grid points was as large as 512 3 , (ii) N was as small as 1 (this was because we are interested in the high wavenumber range for which δ(k) 1) and (iii) the natural viscosity and diffusivity were used (i.e. no turbulence models such as hyper-viscosity/diffusivity were used). To the authors' knowledge, the present DNS is unique in having these three characteristics.
In the DNS, Q ij (k, t), P(k, t) and B(k, t) were given by whereû(k, t) andρ(k, t) (k = k 1 e 1 + k 2 e 2 + k 3 e 3 and k 1 , k 2 , k 3 = 0, ±1, ±2, . . .) are the Fourier transforms of u(x, t) and ρ(x, t), respectively, and Re denotes the real part. We assumed that the ensemble average f K of f summed over a wave vector domain, K, can be approximated by the summation over domain K in one realization, provided that the number of wave vector modes in K is not too small. Therefore, in practice, we used only one realization in the following analysis. We used a statistically quasi-stationary state obtained from a preliminary DNS of forced turbulence without a density field for the initial velocity field of the present DNS. A log-log plot of the initial kinematic energy spectrum E k (k) versus the wavenumber k showed a slope close to that of the k −5/3 spectrum near k = 10 (figure omitted). The initial density field was generated randomly under the constraint (kinetic energy) : (potential energy) = 2 : 1, within each shell k − 1/2 k < k + 1/2 in the wave vector space, i.e. E P (k) ≈ (1/2)E Q (k). At the initial instant, u and ρ were almost uncorrelated, i.e. E B (k) ≈ 0.
The DNS was performed up to Nt = 5.0. After a certain initial transient period, E B (k) showed a persistent feature of E B (k) < 0 for small k and E B (k) > 0 for large k. The feature is consistent with previous numerical studies, although the critical wavenumber k c at which E B (k) changes the sign is different in simulations under different conditions (initial conditions; external forces; turbulence models/natural viscosity and diffusivity; and resolutions); in the present DNS, k c k E , where k E = N 3/2 −1/2 , while in the simulation in [15], k c k E (see [12,14,15,17,18] and references therein for more details on the sign of E B (k)). Hereafter, we will focus on the scaling and angular dependence of the spectra, for which the theoretical predictions were given in section 4.
The magnitude of E B (k), as well as those of E Q (k) and E P (k), decayed with time in the present DNS. The spectra E Q (k), E P (k) and E B (k) at different Nt values were quite similar except during the initial transient period. In the following analysis of the DNS data, we used the representative spectra at Nt = 3.0 after the transient period. Some parameter values that characterize the field at Nt = 3.0, as well as the initial field at Nt = 0.0, are listed in table 1, where the total kinematic energy E Q , the total potential energy E P , the energy dissipation rate and the integral length scale L 0 are defined by where u 2 = 2E Q /3. The Taylor microscale λ, the Kolmogorov microscale η and the Taylor microscale Reynolds number R λ are given by According to K41, the nonlinear interaction between eddies of similar size is dominant in the inertial subrange k E k k ν , k κ , where k E = N 3/2 −1/2 . In our run, k E , k ν and k κ were estimated to be k E = 4.2 and k ν = k κ = 219 at Nt = 3.0.  (d) Strictly speaking, there may be modifications to Q, P and B from the equilibrium spectra in the range owing not only to buoyancy, but also to non-stationarity, in particular those of and χ, as discussed for a turbulent shear flow in [6]. However, the modifications are expected to be small provided that τ N τ t , where τ t is the characteristic time scale of the energy or scalar dissipation rate and is given by τ t = /(d /dt) or τ t = χ/(dχ/dt); τ N is the time scale associated with interactions between eddies of similar size. This may be the case in the wavenumber range k k t , where k t is the wavenumber at which τ N ∼ τ t . The use of the K41 scaling gives τ N ∼ −1/3 k −2/3 . In the present DNS, k t was estimated to be about 1.5 at Nt = 3.0.
Furthermore, the deviation Q, P and B due to the non-stationarity can be shown to be isotropic and of the form to leading order for small τ N /τ t , where q t and p t are appropriate functions of k. In the following, where our main interest is in the anisotropy of the small scales, we ignore the possible effects of non-stationarity. Figures 1(a), (b), (c) and (d) show E B (k), E q1 (k), E q2 (k) and E p1 (k), respectively, together with the isotropic spectra E Q (k) and E P (k) obtained from the present DNS. Equations (30)-(32) and the Kolmogorov scaling (39)-(44) suggest that, in the inertial subrange, to leading order for small δ(k). The spectra (at k ≈ 10) in figure 1 are consistent with, or at least quite close to, these estimates, although the ranges over which the DNS values follow the scaling laws implied by (i)-(iii) are limited, especially for E q1 , presumably because the Reynolds number of the present DNS, where the Taylor micro-scale Reynolds number was R λ ≈ 170, was only moderate.
Next we consider the angular dependence of the spectra. Before comparing the theoretical estimates and DNS values, recall that the θ dependences (27)-(29) and (36)-(38) of the spectra Q ii (k), q 1 (k), q 2 (k), P(k), p 1 (k) and B 3 (k) are obtained without assuming Kolmogorov scaling, i.e. they are free from the Kolmogorov scaling assumption. To examine a theoretical conjecture for the θ dependence, let us define the normalized function X( , θ) as where is the domain [k 0 , k 1 ], dk denotes an integration with respect to k over and X(k, θ) is the average of X(q) over q satisfying k − k/2 q < k + k/2 and θ − θ/2 θ < θ + θ/2 with q 3 = q cos θ and k = 1, θ = 3 • . Figure 2 shows the DNS values of Q ii ( , θ) and P( , θ) for different values of (i.e. k 0 and k 1 ). If X(k) is an isotropic function of k, then X( , θ) must be θ-independent and equal to 1. The θ dependence of the DNS values Q ii ( , θ) and P( , θ) in figure 2 is consistent with (27) and (28) since the values are almost constant and close to 1; however, they have some small θ-dependent deviations from the constant. The large fluctuations observed in Q ii ( , θ) and P( , θ), especially for with small k 0 and k 1 , are presumably due to the lack of the number of wave vector modes in .
To compare the theory and DNS results regarding θ dependences in more detail, let us first consider (29) and (36), according to which When comparing with the DNS values, it is convenient to use X( , θ) rather than X(k, θ)/ X(k) k to avoid large fluctuations in the denominator. Equation (45) yields  Next, consider (37) and (38), according to which where Figure 4 shows DNS values of q A 2 ( , θ) and p A 1 ( , θ) for various values of , together with the theoretical predictions of 5 4 (3 cos 2 θ − 1) obtained from (37) and (38). The DNS values and the theoretical conjecture are again in good agreement with each other, especially for large wavenumbers, similar to figure 3.

Results and discussions
In this paper, we have assumed that, in fully developed turbulence at sufficiently small scales and sufficiently high Reynolds numbers, there exists an equilibrium state, whose statistics is governed primarily by the inherent dynamics, consisting of the convective, pressure and viscous effects, and is insensitive to the flow conditions at large scales. The equilibrium state is generally  modified in response to possible external disturbances, such as the existence of a mean field, magnetic fields, buoyancy and so on. In section 3, by extending the ideas of IYK on the analysis of small-scale anisotropy in turbulent shear flows, we derived expressions to modify the velocity correlation spectrum Q ij (k), the density spectrum P(k) and the buoyancy-flux spectrum B i (k) due to the buoyancy effect in stably stratified turbulence. The expressions give not only the scaling in the inertial subrange, but also estimates of the anisotropy of the spectra in terms of a few scalar functions. The theoretical conjectures were compared with a DNS and found to be in good agreement. Expression (11) with (12) for turbulent shear flow as well as the expressions (24)-(26) for the spectra Q, P and B in stably stratified turbulence are consistent with the LRA. Regarding (11), a theoretical estimate of the two universal functions a(k) and b(k) in the inertial subrange was obtained based on the LRA [6]. An estimate of the universal functions α 1 , α 2 , α 3 , β 1 , β 2 and γ in (39)-(44) for the inertial subrange can also be obtained from the LRA, at least in principle. Such a theoretical estimate remains to be pursued. A preliminary analysis suggests that the k −3 spectra in (42)-(44) can be modified into one that is proportional to k −3 log k, as is the case for the enstrophy cascade range in two-dimensional turbulence [19,20].
The scaling (39)-(44) can also be obtained from simple dimensional arguments under certain weak assumptions. However, the possibility of anomalous scaling, which is not obtained by simple dimensional arguments, cannot be excluded a priori. In fact, (i) a closure approximation under an appropriate linearization with respect to small perturbations in the velocity correlation spectrum, as well as (ii) a simple passive vector model with pressure-like effects, which are  generalizations of the well-known Kraichnan's passive scalar model [21,22], suggest the possibility of anomalous anisotropic scaling of the velocity correlation spectrum in turbulence that is almost isotropic [23]. Clarification of the possible anomalous scaling in turbulence with and without mean flow, buoyancy, etc is an interesting challenge that requires further theoretical studies as well as experiments and DNSs at high Reynolds numbers.