Lukas Liechtenstein
email: lukas.liechtenstein@ec-lyon.fr

Fabien S Godeferd

Claude Cambon

Nonlinear formation of structures in rotating stratified turbulence

Keywords:

published or not. The documents may come

Introduction

Coherent structures are widely observed in instantaneous velocity fields in anisotropic statistically homogeneous turbulence (see e.g. [START_REF] Kimura | Particle dispersion in rotating stratified turbulence[END_REF]). The average shapes vary with respect to the statistically isotropic case depending on the nature of external distortions applied onto the flow. The structures can be visualized by plotting instantaneous Eulerian fields using an identification criterion chosen among many, for instance by plotting constant enstrophy surfaces as illustrated in figure 1. One can identify flattened objects in dominantly stratified turbulence ("pancakes", figure 1(a)) and vertically elongated structures in dominantly rotating turbulence ("cigars", figure 1(c)). Although the different scale elongations and aspect ratios for different cases can easily be identified at a glance, they can hardly be quantified merely as a function of the parameters used, such as the rotation rate or the stratification density gradient. We attempt to gain further insight into the evolution and scalings of these structures by establishing links between instantaneous structures and anisotropic statistics of the velocity field, especially considering Lagrangian and Eulerian second-order statistics. A key in this analysis is the role of nonlinear terms in Lagrangian versus Eulerian statistics, as previous studies have shown that structure formation is a nonlinear process (see e.g. [START_REF] Cambon | Energy transfer in rotating turbulence[END_REF] for the purely rotating case, or e.g. [START_REF] Riley | Direct numerical simulations of homogeneous turbulence in density-stratified fluids[END_REF][START_REF] Godeferd | Statistical modelling and direct numerical simulations of decaying stably-stratified turbulence: Part 2: Large and small scales anisotropy[END_REF] for the purely stratified one).

We consider here statistically homogeneous turbulence with vertical stable stratification characterized by the Brunt-Väisälä frequency N and vertical system rotation with frequency Ω in the Boussinesq approximation [START_REF] Boussinesq | Théorie analytique de la chaleur, mise en harmonie avec la thermodynamique et avec la théorie mécanique de la lumière[END_REF]. The system of equations in physical space is ∂u ∂t

+ u • ∇u -ν∇ 2 u = -∇p -f n × u + bn, (1)
∂b ∂t + u • ∇b -χ∇ 2 b = -N 2 n • u, (2)
∇ • u = 0 (3) g z x ρ N (a) Ω n u F c (b)
Figure 2. (a) Oscillating mechanism with stable stratification: a fluid particle in a negative density gradient ∂ρ/∂z under the influence of gravity g oscillates with frequency N . (b) Oscillating mechanism with system rotation: the Coriolis force F c acts on a fluid particle with a horizontal velocity u in a vertically rotating frame of reference with frequency Ω.

for the fluctuating velocity u and the buoyancy b, using a notation similar to that used by Cambon [START_REF] Cambon | Turbulence and vortex structures in rotating and stratified flows[END_REF], introducing the Coriolis parameter f = 2Ω. In addition to its relation with experimental studies for rotating stably stratified turbulence, our work is closely related to numerical and theoretical studies based upon equations ( 1)- [START_REF] Riley | Direct numerical simulations of homogeneous turbulence in density-stratified fluids[END_REF], with an emphasis on statistics and dynamics (see the recent reviews [START_REF] Staquet | Sommeria Internal gravity waves: from instabilities to turbulence[END_REF] and [START_REF] Cambon | Turbulence and vortex structures in rotating and stratified flows[END_REF]).

Apart from the compensated pressure gradient, the right-hand side terms containing the anisotropic parameters in (1) and (2) are linear. Dropping the nonlinear convective terms, they alone yield the linearized version of the Boussinesq approximation equations, which admit oscillating solutions due to physical mechanisms schematically shown in figure 2. Figure 2(a) describes this mechanism due to stratification, which generates buoyancy oscillations: fluctuations of density or vertical velocity on a fluid element in hydrostatic equilibrium will force a fluid element to regain its level by undergoing vertical oscillations. Turbulence with system rotation will be subjected to the Coriolis force F c , which forces horizontally moving fluid elements on circular trajectories, as shown in figure 2(b). Without pressure fluctuations, which ensure divergence-free properties, the oscillations in figures 2(a) and 2(b) stay strictly vertical and horizontal respectively. Complete linear solutions with the divergence-free properties of the flow allow more complex vertical/horizontal coupling, as studied in details by e.g. [START_REF] Hanazaki | Linear processes in unsteady stably stratified turbulence[END_REF] for the purely stratified case only. By combining the divergence-free velocity and buoyancy field into one vector v = (û (1) , û(2) , û(3) ) represented in a local reference frame in Fourier space, the so called Craya-Herring frame shown in figure 3(a), the initial fivecomponent problem (three velocity components, buoyancy, and pressure) reduces to a three-component solution (two velocity components plus buoyancy), as described in more detail in Appendix A. The unsteady parts of the solutions are internal propagating plane waves, also called inertio-gravity waves [START_REF] Cambon | Turbulence and vortex structures in rotating and stratified flows[END_REF]. A given initial velocity field u(x, t = 0) and buoyancy field b(x, t = 0) combined into the three-component vector v(k, t = 0)) therefore determines the general non viscous linear solution

v(x, t) = ǫ=0,±1 N ǫ e -iǫσt N -ǫ • v(k, t = 0) e ik•x d 3 k, (4)
within the Boussinesq approximation, by using the eigenmodes N ǫ , ǫ = 0, ±1 of the linear operator from equation (A.9). N 0 represents the non oscillating contribution to the velocity field in the linear limit, and N ±1 are the oscillating solutions, i.e. the inertio-gravity waves. Their dispersion relation

σ = (N 2 sin 2 θ + f 2 cos 2 θ) 1/2 (5)
as a function of the angle θ between the wave vector k and the vertical illustrates the importance of the ratio of rotation to stratification

α = f N . ( 6
)
The dynamics of the flow depends strongly on the non-dimensional parameter α, in addition to the non-dimensional parameters linked to the external forces: the Froude number Fr = q/N L and the Rossby number Ro = q/(f L), where q and L are velocity and length scales of the turbulence.

The case where σ is constant, i.e. α = 1, exhibits a very specific behaviour. It is not treated separately here, but nevertheless shown for comparison, as in figure 1(b), in which no specific structuring can be noted.

The linear solution (4) of the Boussinesq approximation gives us the possibility to analytically predict the velocity as a superposition of traveling plane waves and steady modes. As a result, long term statistical anisotropy for single-time second-order correlations can only appear by nonlinear mechanisms, unless it is already present in the initial conditions.

The definition of N ǫ , ǫ = 0, ±1, depends on f /N , and involves both velocity and buoyancy. On the one hand, they can be expressed in a fixed frame of reference. The contribution to the velocity field from N 0 is purely horizontal, while the wavy contribution from N ±1 can be horizontal, due to rotation only, as well as vertical, due to both rotation and stratification. On the other hand, the eigenmodes are more easily expressed in the Craya-Herring frame of reference (figure 3(a)). In this frame local to k (see Appendix), only two components appear for the velocity part, such that û(k) = û(1) (k)e 1 (k) + û(2) (k)e 2 (k), and a third component associated with buoyancy. This decomposition is generic for any divergence-free velocity field, and related to a classical toroidal-poloidal one in physical space.

From its definition, one computes the spectral vorticity as

ω(k) = ik û(1) (k)e 2 (k) -û(2) (k)e 1 (k) (7)
and derives û(1) (k) and û(2) (k) as a function of the vertical vorticity and vertical velocity (denoted with a subscript, and the orthogonal wave number is According to (4), turbulence with stable stratification and system rotation can be considered as a two-mode motion, one 'vortex' mode at ǫ = 0 and one 'internal wave' mode at ǫ = ±1. The vortex and internal wave mode are linearly independent, only linked through triadic exchange terms. Linear dynamics only affect the wave mode, and consist of periodic exchanges of energy between the wave part of the velocity (identical to the poloidal component only for pure stratification) and the buoyancy field, if stratification is present. In the case of pure rotation, only kinetic energy is concerned, and inertial waves exchange energy between poloidal and toroidal components of the velocity field.

k ⊥ = |k × n|) û(1) = i ω k ⊥ and û(2) = k k ⊥ û . (8)
We discuss two guidelines for a natural decomposition of the velocity field, which are not the same in a fundamental way. On the one hand, from the solution of the linear equations one can consider the velocity field as a composition of inertio-gravity waves and vortex modes. On the other hand, the velocity field can be divided into purely mathematical toroidal and poloidal parts. (See figure 3(b) for the differences between the two decompositions.) In the purely stratified case, the toroidal part is the non propagating vortex mode. For pure rotation, the toroidal mode is part of the propagating wave mode and the vortex mode does not exist. For cases in between, one part of the toroidal contribution is propagating, while the other part is non propagating.

In the following, all Eulerian fields are split using the toroidal/poloidal decomposition and discussed accordingly, as rotation and stratification have significantly different effects on these physically different velocity components. The article will move on to numerical methods in section 2, and then show results for Eulerian statistics as a function of time in section 3, for Lagrangian statistics in section 4, and for Eulerian structures and spectra in section 5. Finally, conclusions will be given in section 6.

Numerical Methods

We generate Lagrangian statistics using three different methods: a fully nonlinear direct numerical simulation (DNS); an analytical linear model using equation (4) for calculating two-time second-order Eulerian statistics (RDT/SCH); and a kinematic simulation (KS) based on a model using random Fourier modes for the generation of turbulent velocity fields, where the linear time evolution of the Fourier modes from equation ( 4) is exactly incorporated. In KS and DNS, trajectories are numerically computed and then used for calculating Lagrangian statistics. In the RDT/SCH, two-time second-order Eulerian statistics are put equal to their Lagrangian counterparts.

Eulerian statistics such as directional integral length scales or direction-dependent spectra are only evaluated from DNS, as KS is only designed to follow fluid elements and cannot modify by itself the input spectral distribution, which is chosen to be isotropic. The reversible anisotropy produced by the linear method and KS do not generate structures, as mentioned in section 1. We have not initiated KS with explicit initial anisotropy so as not to introduce an a priori additional anisotropy parameter, which would have been chosen arbitrarily.

Direct numerical simulation

Equations (1)-(3) are solved directly using a pseudo-spectral collocation method expressed in Fourier space following a classical scheme as e.g. in [START_REF] Rogallo | Numerical experiments in homogeneous turbulence[END_REF][START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF]. The computational domain uses periodic boundary conditions, and 256 3 points in physical space. Furthermore, the velocity field is completely de-aliased using a 2/3-truncation method in Fourier space. The initial energy spectrum is the narrow band distribution E(k) ∝ k 4 e -2(k/k i ) 2 , peaking at k i , around one tenth of the maximum wavenumber of the de-aliased field. The initial velocity vector directions are distributed randomly in space and therefore do not correlate with each other after initialization. From this initially random isotropic field we perform an isotropic precalculation and so provide a realistic velocity field as the starting point of the anisotropic runs. In fact, the velocity field becomes more "realistic" as we create higher order velocity correlations and built up isotropic energy dynamics during the precalculation. At the beginning of the anisotropic runs we allow for discontinuities in the statistical data derivatives due to the introduction of an anisotropic body force into the system.

Eulerian temporal statistics are calculated during the calculation, as saving of the velocity field at all time steps is too memory consuming. Only the last velocity field is saved, to calculate spectra and isosurfaces in a post-processing step.

Trajectories in DNS are obtained by solving the fluid particle motion equation; the velocity at the location of the particle is computed from the known Eulerian velocity field using a standard spatial interpolation scheme of 6-th order with Lagrange polynomials. Time marching is done with an intermediate time step to ensure the stability condition for the particle equation. The raw data of the trajectories are stored and processed at a later stage to calculate Lagrangian statistics.

Run Resolution Duration

ν f N Re λ Initial PE A1 256 3 30T N/2π 1/500 0 8π 39 E pot = E pol B1 512 3 6.5T N/2π 1/1000 π 8π 51 E pot = E pol C1 256 3 30T N/2π 1/500 8π 8π 38 E pot = E pol D1 256 3 32.5T f /2π 1/500 10π π 41 E pot = E pol E1 256 3 24T f /2π 1/500 8π 0 42 E pot = 0 B2 256 3 8T N/2π 1/1500 π 10π 76 E pot = 0 C2 256 3 8T N/2π 1/1500 10π 10π 80 E pot = 0 D2 256 3 8T f /2π 1/1500 10π π 83 E pot = 0 Table 1.
Sets of parameters of the DNS runs presented.

Kinematic simulation

The kinematic simulation (KS) models a turbulent velocity field as a superposition of random Fourier modes so that the field is automatically incompressible. No dynamical equation is involved in the basic version of KS, although oscillations physically related to internal waves can be explicitly introduced for each mode using the linear solution from equation ( 4) [START_REF] Nicolleau | Turbulent diffusion in stably stratified non-decaying turbulence[END_REF]. Due to the statistical axisymmetry of the here treated velocity field the choice of random wave vectors is restricted to one random Fourier mode to represent all azimuthal angles for one realization. A carefully chosen discretization is used for representing modes at different polar angles (directional asymmetry) and for different wave vector lengths k. Due to the lack of complete representation, different realizations (≈ 50) ensure good average statistics. By utilizing an initial velocity field û(i) (k, 0) in the Craya-Herring frame of reference, the velocity components at time t = 0 for a point x in real space can be written down as a discrete Fourier sum

u i (x, 0) = ℜ N k n=1 M θ m=1 û(1) (k mn , 0)e (1) i (k mn ) + û(2) (k mn , 0)e (2) i (k mn ) e ikmn•x (9)
with the definitions of e (i) i and û(i) i given in Appendix A. This formulation implies the existence of exactly one k for each resolved wave vector length k n and each polar angle θ m , but with a random azimuthal angle φ representing all azimuthal modes. The total number of resolved modes is therefore N k M θ .

By assuming that the linear evolution is comparatively faster than de-correlation due to the nonlinear evolution, one can deduce a quasi-linear time evolution of the velocity field where the time evolution of û(i) (k mn , t) is given by the first two components of a discrete version of equation ( 4) as

u i (x, t) = ℜ N k n=1 M θ m=1 e ikmn•x û(1) (k mn , t)e (1) i (k mn ) + û(2) (k mn , t)e (2) i (k mn ) (10) Run M θ N k Duration of run f N k max /k i Initial PE F 50 400 200T N/2π 0 16π 1000 E pot = 0 G 50 400 200T N/2π 2π 16π 1000 E pot = 0 H 50 400 200T N/2π 16π 16π 1000 E pot = 0 K 50 400 200T f /2π 16π 2π 1000 E pot = 0 L 50 400 200T f /2π 16π 0 1000 E pot = 0
v(k mn , t) = ǫ=0,±1 N ǫ m e -iǫσmt N -ǫ m • v(k mn , 0) , (11)
noting σ m the dispersion frequency associated with wavevector k mn . KS can then be used as a Lagrangian model of turbulent diffusion by solving the fluid particle motion equation

ẋ(t) = u(x(t), t), (12)
where the velocity u is found at any time T and any point x in real space by using Fourier sums of the random wave vectors as in equation [START_REF] Vincent | The spatial structure and statistical properties of homogeneous turbulence[END_REF].

The time integration of equation ( 12) is done by discretization of the trajectories of the fluid elements by constant time steps and extrapolating the new position of a Lagrangian fluid element by an Adams-Bashforth-Moulton predictor-corrector algorithm of the fourth order. The spectral energy distribution, which is fixed throughout time by the algorithm, is k 3 for k min < k < k i and k -2 for k i < k < k max . The values for the gradients of the spectra are chosen in approximation to values from the DNS. The values for k min , k i and k max are 0.1, 1 and 1000 respectively. Table 2 gathers the parameters of the KS data.

Linear model: RDT/SCH

The solution (4) yields the exact linear solution, provided the initial velocity field is known. It is here called rapid distortion theory (RDT) as its derivation is closely related to the short-time approximation of the Boussinesq equations [START_REF] Kimura | Particle dispersion in rotating stratified turbulence[END_REF]. As in classical RDT, not only the velocity field may be computed at each time directly from (4), but one can obtain as well analytical expressions for the single-point statistics of the turbulent field, upon integration of two-point correlations of û over the complete spectral space. This yields the kinetic energy spectra, but also the potential energy, and eventually, if needed, two-point two-time spectra of the horizontal or vertical velocity components, as well as those of the poloidal and toroidal components. Unfortunately, this method only applies to Eulerian correlation spectra, but can hardly be used as such for obtaining Lagrangian ones, without further hypotheses. Recent works have demonstrated the validity of one additional assumption to Taylor's temporal method of integration for single particle dispersion (as in equation ( 19)). It is a simplified Corrsin hypothesis (SCH), which allows to replace the Lagrangian velocity correlations by second order correlations statistics obtained from RDT, which are evaluated at two times. This gives an analytical expression for calculating one-particle Lagrangian displacement correlations using two-time velocity correlations ( [START_REF] Kaneda | Suppression of vertical diffusion in strongly stratified turbulence[END_REF][START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF]). Details of the method are explained in [START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF], which also describes the general results for the stratified/rotating case with arbitrary initial partition of potential and kinetic energies. As an example, let us write down only the vertical displacement correlation ∆ 33 =< (∆x 3 ) 2 > at time t for a case with equipartitioned energy (E pot (0) = E pol (0) = E tor (0)) with dispersion relation σ:

∆ 33 (t) = 2E kin (0) 1 0 1 -cos 2 θ 1 -cos σt σ 2 d(cos θ), (13)
where only the initial kinetic energy E kin (0) appears, in addition to the dispersion relation σ. In the integral [START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF], the directivity of the flow is felt through the presence of the trigonometric functions in the dispersion relation, which alone determines different evolutions of single particle dispersion with this model.

Temporal evolution of Eulerian statistics in DNS

Three DNS runs (B2, C2, D2) are presented: one dominantly stratified case (α = 0.1, red); one dominantly rotating case (α = 10, blue); and a case with constant σ α = 1 (green), as in this case the dispersion relation reduces to σ = 1. The parameters for the different runs are shown in Table 1.

All DNS are freely decaying as anisotropy in freely decaying DNS can evolve more naturally than in forced DNS. Some test runs (not presented) where done for forced DNS with various methods of forcing, in the large or small scales domain, or only in the horizontal or vertical direction. The type of forcing influences so much the anisotropy, which is a priori not known, that a directionally "natural" forcing can hardly be done. Furthermore, to show the exchange of energy between the poloidal mode and the potential energy, the runs B2, C2, and D2 are initialized with zero potential energy. Figure 4 shows the energy decay of the runs. The decomposition of the energy in E tor , E pol and E pot in figures 4(b), (c) and (d) add up to the total energy in figure 4(a). The decay rate at the end of the run is around t -0.5 for the very rapidly rotating case D2 and the strongly stratified case B2, so significantly slower than in isotropic turbulence, which roughly decays with t -1.4 . The decay rate for the case with constant σ C2 approaches unity. In the cases B2 and D2 the dispersion law depending on θ alters the energy exchanges in an anisotropic way and thereby slows down the decay of energy, an important feature in anisotropic turbulence. Note that for the total energy in figure 4(a) no oscillations are found in all three cases, as the exchange of energy between the potential mode and the poloidal mode is added up. The purely kinetic energy composed of only toroidal and poloidal part oscillates, as part of the total energy is stored as potential energy. The potential energy in figure 4(d) for C2 is significantly lower than for B2 and C2, as stratification is comparatively weak. The poloidal and the potential parts exchange energy with linear mechanisms, which can clearly be seen by the complementary nature of the oscillations with a frequency N for B2 and C2 in figures 4(c) and (d). The oscillations are damped by linear effects (angular phase mixing) and also by nonlinear transfer mechanisms, which explain the oscillations for C2 in the toroidal energy in figure 4(b).

To sum up, the mechanism of the energy equilibrium with initially zero potential energy for the three cases can be explained in the following way: for significant stratification, there exists an immediate linear rapid exchange of energy between the poloidal and potential parts, which tends to distribute the available energy evenly on both parts by exhibiting oscillations of frequency N . This exchange is damped by linear (phase mixing) and nonlinear (phase mixing, triadic exchange term) mechanisms, which also transport energy to the toroidal mode. For weak stratification and significant rotation, the mechanism is the same, but the amount of available potential energy is smaller and the period of the oscillations is longer due to comparatively smaller values of N . Nevertheless, for this case there is a linear exchange of energy between the toroidal and poloidal mode, keeping them at the same level at all times.

For the case α = 1, the energy does a double linear exchange, first from poloidal to potential energy, then from toroidal to poloidal energy. Therefore, oscillations can be seen in all three modes, though they are smallest for the toroidal part and largest for the potential part.

Anisotropy in the energy due to stratification is obvious, if one compares the poloidal with the toroidal part. However, this anisotropy is generated rather by linear mechanisms due to the initialization of zero potential energy and so contributes only marginally to anisotropic structure formation. When looking at energy evolution of cases with initial equidistribution of energy between the three modes (not shown), anisotropy due to nonlinear effects is relatively weak, with a slightly higher kinetic energy in the horizontal direction.

To quantify the scales of structures, it is better to look at anisotropic statistical measures, such as the integral length-scales

L n ii (t) = 1 < u i u i > ∞ 0 < u i (x)u i (x + r n ) > dr n (14)
for different velocity components u i and different separation directions n (vertical and horizontal). In figure 5(a) we show the correlation length scales for the horizontal velocity u x in the horizontal and the vertical direction. With dominant stratification, the elongation in the horizontal (increase in L x xx , dashed red curve) as well as the reduction in the vertical (decrease in L z xx , solid red curve) length scales clearly illustrates the pancake like structures observed for B2 in figure 1(a). Even more dramatically the scale elongation of L z xx (solid blue curve) for D2 illustrates the formation of vertically elongated structures in figure 1(c) for turbulence with dominant rotation, while the horizontal structures (dashed blue curve) do not evolve in size.

Although anisotropy and scale elongation can qualitatively be compared to the visualized structures in figure 1, the ratio of scales does not up to now reach asymptotic values.

To find quantitative aspect ratios of vorticity isosurfaces, vorticity length scales are investigated. They are defined similarly to equation ( 14), but by correlating vorticity instead of velocity. Velocity correlation length scales indicate a large scale, in particular, when looking at an energy spectrum, the scale which is cascading with k -1 . Similarly, dimensionally, the vorticity correlation length scales are characteristic of scales where the energy spectrum is decaying with k -3 . In isotropic turbulence, this scale will be found in the viscous range. In anisotropic turbulence, as shown in section 5, this needs not be the viscous scales, but is nevertheless difficult to locate.

However, the direct link between these length scales and vorticity isosurfaces is useful to gain insight into anisotropic coherent structures. In figure 5(b) we show L z ωzωz (solid lines) as well as L y ωxωx (dashed lines). The most surprising feature of these curves is that they are complementary, if one swaps the values of the parameters N and f . This is not at all true for velocity correlation lengths, as seen in figure 5(a). As the nonlinear mechanisms creating these anisotropies are very different, this might only be a coincidence. Moreover an asymptotic value seems to be reached for L z ωzωz for dominant stratification as well as L y ωxωx for dominant rotation. This is in agreement with an asymptotic vertical Froude number, or an asymptotically fixed width of layers for stratification or fixed diameter of vortices for rotation. The horizontal diameter of the pancakes as well as the vertical length of the cigars seems to grow linearly.

Rapid Distortion (linear) Theory (RDT) applied to initial isotropic data with equipartition of poloidal and potential energy (the latter condition just to simplify without essential lack of generality) strictly conserves isotropy of any single-time double correlation. For example, the isotropic relationships are not altered for the integral length scales mentioned above. They are equal in horizontal and vertical directions, in contrast to results from full nonlinear statistical spectral closures which are in excellent agreement with high resolution DNS [START_REF] Godeferd | Detailed investigation of energy transfers in homogeneous stratified turbulence[END_REF][START_REF] Godeferd | Statistical modelling and direct numerical simulations of decaying stably-stratified turbulence: Part 2: Large and small scales anisotropy[END_REF]. The differences between single-time linear and nonlinear behaviour indicate a strong influence of nonlinearity in the formation of structures. Accordingly, it is shown in DNS that the structures take considerably more time to form, when increasing the parameters of stratification N , the Brunt-Väisälä frequency and Ω, the system rotation frequency, and therefore increasing the dominance of linear terms.

∆ x (t, t ′ ) ∆ y (t, t ′ ) ∆ z (t, t ′ )
x(X, t)

x(X, t ′ ) s Figure 6. Schematic plot of single particle dispersion in three-dimensional space for one fluid element.

Lagrangian structures

Already linear models break isotropy, if one looks at Lagrangian or two-time Eulerian statistics. In this section we use results from RDT evaluated at two times (linear), Lagrangian statistics from KS (random Fourier modes) and DNS (nonlinear), details of which are given in section 2. We analyze the most basic Lagrangian quantity, the single particle dispersion ∆x, schematically shown in figure 6 for one fluid element. The Lagrangian position of a fluid element labeled by the initial position X

x(t) = x(X, t) (15)
has a Lagrangian velocity related to the Eulerian velocity field u(x, t) by

V (t) = u(x(X, t), t). ( 16
)
Therefore the position of the particle advected by the Lagrangian velocity field can be written down as

ẋ(t) = u(x(X, t), t) (17)
which subjects the Lagrangian position of the particle x(t) to feedback by itself and consequently to a nonlinear evolution. By integration along the trajectory of a fluid element, one gets mean displacements along each i-th direction

∆x i (t, t ′ ) = x i (t) -x i (t ′ ) = t t ′ ẋi (s)ds. ( 18
)
which, as covariances, give single particle dispersions [START_REF] Taylor | Diffusion by continuous movements[END_REF] ∆

ii (t, t ′ ) =< ∆x i (t, t ′ ) > 2 = t t ′ ds ′ t t ′ < ẋi (s) ẋi (s ′ ) > ds. ( 19
)
in all three space directions, i = 1, 2, 3. For the analytical linear method, the Lagrangian velocity correlations in equation ( 19) are replaced by their Eulerian counterparts derived from Rapid Distortion Theory (RDT) [START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF], following the simplified Corrsin hypothesis (SCH) (see section 2).

For studying anisotropy, it is interesting to compare the single particle dispersion in the horizontal and the vertical directions separately. The vertical single particle dispersion ∆ 33 (0, t) =< (x 3 (t 0 ) -x 3 (t)) 2 > is shown in figure 7 tN/2π RDT/SCH, the KS for runs F to L and the DNS for runs A1 to E1. The similarity of the three figures is surprising, when one considers the differences in the methods with which the results have been obtained. Three kinds of regimes may be identified: a ballistic regime (∆ 33 ∝ t 2 ) at small times, valid for all the values of α and all models, then at longer times a weak Brownian regime for rotation dominant cases, such that ∆ 33 ∝ t, and finally an oscillating plateau for cases with non-zero stratification, even if rotation is dominant, at a time of around t = N/2π. The plateau illustrates a confinement of vertical displacement for fluid particles, which scales with any non-zero value of N .

t 2 t c) DNS Nonlinear α = 0 α = 1/8 α = 1 α = 8 α = ∞
Nevertheless, there are differences between the linear, the random and the nonlinear evolutions. The single-time Reynolds stresses are isotropic in RDT as well as in KS, and therefore the maximum vertical dispersion scales with N and q 2 . In DNS, a mean value leads to a similar plateau, though scaling with N and 3u 2 3 instead of q 2 , which is different due to the anisotropy of the Reynolds stress tensor in the Eulerian velocity field generated by DNS.

The oscillations in the linear method are very regular as well as being damped at a lower rate. The damping with this method is strictly linear due to phase mixing. Any nonlinear triadic interaction is excluded, which explains the regularity and slow decay of the oscillations. Furthermore, there is no phase mixing for the case α = 1, which means that the oscillations for this case for the RDT/SCH are undamped. Accordingly, the ascent of the one-particle dispersion seen for the case α = 1 in DNS is a nonlinear effect. As the dynamics is strictly linear in KS, we attribute the similar evolution of the one particle dispersion to the nonlinearity in the fluid particle motion equation. This "nonlinear" tendency of a reduced confinement is confirmed by simulations of both KS and DNS with higher Froude numbers. The transition from the ballistic to the Brownian regime is observed at a time of around t = f /2π, as the cases α = 8 and α = ∞ show. At later times, in DNS (figure 7(c)), the case α = ∞ increases its rate of vertical dispersion, probably a numerical artefact.

The horizontal dispersion is shown in figure 8. Again, the similarity is outstanding, though this time figures 8(a) and (b), from the linear method and the random model, seem to yield more similar results. As for the vertical dispersion, in all three cases the horizontal one-particle dispersion ∆ 11 (0, t) =< (x 1 (t 0 ) -x 1 (t)) 2 > shows a t 2 law at small times. With rotation one finds a stage with Brownian (∆ 11 ∝ t) behaviour at about 0.1 integral time scales, which, depending on the amount of stratification, returns to a t 2 time evolution. Contrarily to single particle dispersion in isotropic turbulence, a Brownian behaviour expected at very large times is not observed in any of the cases. The differences between the three figures is mainly the scaling, done here with the initial integral length scale and initial mean velocity. Furthermore, the Brownian transient dispersion of the DNS data is not as long as in the other two cases. This might be a low Reynolds number effect.

For pure rotation, the Brownian regime seems not to go back to a ballistic regime for the linear data. Furthermore, the vertical diffusivity is exactly twice the horizontal one [START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF], a ratio which is approximately recovered for the nonlinear data.

Eulerian structures and statistics

The most detailed second order statistics are given by two-point velocity correlations, or the related spectral tensor. For any anisotropic homogeneous turbulent flow, this spectral tensor can be reduced to an optimal set of three components, one of which is the kinetic energy density tensor e(k), the other two providing a very fine grained description of the anisotropic structure of the flow (details in [START_REF] Cambon | Energy transfer in rotating turbulence[END_REF]). In our case, axisymmetry means that e depends on the polar angle θ only, and its departure from a spherical equidistribution characterizes directional anisotropy. Hence, the angular distribution of e is a key information which characterizes cigar-or pancake-shaped structures in physical space [START_REF] Godeferd | Detailed investigation of energy transfers in homogeneous stratified turbulence[END_REF][START_REF] Cambon | Energy transfer in rotating turbulence[END_REF][START_REF] Godeferd | Statistical modelling and direct numerical simulations of decaying stably-stratified turbulence: Part 2: Large and small scales anisotropy[END_REF]. Equivalent information is contained in the distribution of the poloidal and toroidal kinetic energy density spectra, presented hereafter, since e(k, θ) = e tor (k, θ) + e pol (k, θ).

k x k y k z û(k(n, m)) k n θ m A (a) θ k z (b)
We analyze here the anisotropy of the velocity field by using energy spectra as a function of wave number k and the azimuthal angle θ measured with the vertical (polar) Fourier direction [START_REF] Godeferd | Statistical modelling and direct numerical simulations of decaying stably-stratified turbulence: Part 2: Large and small scales anisotropy[END_REF]. The energy is integrated over cuts of spherical shells along the polar angle direction. In a discrete flow field, for n shells of width ∆k and radius k n divided into m equal polar sectors of angle ∆m and angle θ m we can write down an expression for the directional kinetic energy spectrum

E(k n , θ m ) = m 2 θm+∆θ/2 θm-∆θ/2 cos θdθ -1 |k|∈In,θ k ∈Jm û * (k) • û(k). ( 20
)
It corresponds to an integration of the blue surface in figure 9(a) over a shell width ∆k, which is the piecewise equivalent to the 4πk 2 weighting used to obtain spectra from spectral densities. The spectra represent classically integrated spherical spectra, but instead of interpreting them along a ray from the center of the sphere, they represent spectra along a side of cones with different opening angles, schematically sketched in figure 9 By using a Craya-Herring type decomposition of energy in a toroidal and poloidal mode (details in Appendix A or e.g. in [START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF]) it is observed that most of the anisotropy in dominantly stratified cases is found in the toroidal mode, since θ increases from top to bottom for the directional energy spectrum curves. The poloidal mode distributes the turbulent kinetic energy in a much more isotropic way. Furthermore, the largest anisotropy can be found at intermediate wave numbers, where one would expect an inertial range for higher Reynolds numbers. The energy in the polar direction clearly dominates, illustrating a quasi two-component flow in the horizontal direction. Linear dynamics, e.g. from Rapid Distortion Theory (RDT) do not manage to create such anisotropic effects. The toroidal and the poloidal mode are linearly independent, linked through a nonlinear triadic exchange term. Linear dynamics only affect the poloidal mode periodically exchanging energy between the poloidal part of the velocity and the buoyancy field. RDT applied to initial isotropic data with equipartition of poloidal and potential energy strictly conserves isotropy of any single-time second-order correlation. Dominantly rotating cases show similar anisotropy in toroidal and poloidal parts of the energy with a maximum at the smallest scales, θ decreasing from top to bottom. The isotropy at large scales is much more conserved than in dominantly stratified cases. However, at small scales the energy density is much higher in the polar direction, indicating a high degree of correlation for the velocity in the real-space vertical direction. This indicates a trend toward a two-dimensional three-component flow field correlated in the vertical direction. Linear inertial waves exchange kinetic energy between the poloidal and the toroidal part of the velocity field and therefore cannot explain the anisotropy at large k. The role of nonlinear interactions with little involvement of the poloidal mode and the related angular energy drain (rather than a classical cascade) is confirmed, suggesting a simplified scheme supported by both DNS results and statistical models [START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF].

Furthermore, weak rotation with dominant stratification as well as weak stratification with dominant rotation is not detectable. As the terms containing the anisotropy parameters are linear in the Boussinesq system of equations, the mixed cases with domination of stratification or rotation are very similar to the pure cases. Effectively, the influence of rotation in figures 10(a) and (d), as well as the influence of stratification in figures 10(c) and (f) can not be detected. The case α = 1 shows results similar to the isotropic case (black curves in all figures), meaning that in this case the inertio-gravity waves lead to a compensation of the two mechanisms of anisotropy.

Conclusions and open issues

There is a paradox: in a stratified fluid, the layering or pancake dynamics, observed from velocity snapshots or single-time Eulerian statistics, is a nonlinear phenomenon, as are the columnar structures in rotating turbulence. How can this be partly predicted by trajectories created with linear dynamics with KS, or two-time statistics calculated in the RDT approximation?

In linearly generated flow fields, this kind of anisotropy is missing in the Eulerian velocity fields, but their Lagrangian statistics do exhibit anisotropic features. We put forward two means by which the anisotropy of the velocity field, if absent, can be mimicked in the Lagrangian statistics.

First, two-time second-order velocity statistics indeed show isotropy breaking, which is amplified by a time integration, as performed in previous works about dispersion [START_REF] Taylor | Diffusion by continuous movements[END_REF][START_REF] Kaneda | Suppression of vertical diffusion in strongly stratified turbulence[END_REF]. The latter work also incorporates a procedure for turning Eulerian two-point two-time correlation spectra into Lagrangian ones [START_REF] Cambon | Turbulent diffusion in rapidly rotating flows with and without stable stratification[END_REF][START_REF] Kaneda | Suppression of vertical diffusion in strongly stratified turbulence[END_REF], hence exemplifying a possible link between the two. Second, the fluid particle trajectory itself may be a source of nonlinearity which can introduce irreversible anisotropic features into Lagrangian statistics. The intrinsic nonlinearity in the Lagrangian equation for the particle position, ẋ(X, t) = u(x(X, t), t), might be additionnally fed by the nonlinearity of the advection term in the Navier-Stokes equation. Using the results of Taylor [START_REF] Taylor | Diffusion by continuous movements[END_REF], this argument can be extended to the two-time statistics of RDT/SCH. Strictly speaking, the trajectory argument does not hold per se for the purely analytical linear model, although RDT appears to be informative enough for Lagrangian one-particle anisotropy, and this comes in support to the simplified Corrsin hypothesis.

Let us discuss the first point in more detail. Two-time Eulerian and Lagrangian double velocity correlations incorporate more information than single-time (only Eulerian) double correlations. Isotropy breaking by the presence of unsteady wave and steady vortex modes, arising from a purely linear operator (e.g. our equation ( 4)), is reflected significantly by two-time velocity correlations. This anisotropy is then magnified by the Taylor time-integration when calculating Lagrangian diffusivities in the "RDT+SCH" model.

If single-time double correlations are considered, no significant anisotropy can result from the RDT linear solution; isotropy is broken in single-time third-order correlations, so that anisotropy for Eulerian statistics is triggered by nonlinear transfer terms mediated by triple correlations. Accordingly, nonlinear dynamics are essential for a correct anisotropic description of single-time Eulerian statistics, and for the related formation of organized structures as exhibited by DNS snapshots.

In our observations of the linear regime, dispersive waves appear more efficient than "soft" organised structures to render the Lagrangian diffusivity anisotropic. From this viewpoint, the trapping and confinement of trajectories by dispersive waves can be similar to the effect of a robust vortex. If we focus on vertical rms displacement, this efficiency depends only on the dispersion law: confinement in the vertical direction with plateau (figure 7) is found in the stratified case, but not in the purely rotating case. Hence Lagrangian anisotropy is found to be maximum in the stratified case, since vertical motion is completely dominated by waves inducing a plateau, in contrast with the steady "vortex" part of the horizontal motion which is responsible for free dispersion in horizontal directions. In the pure rotating case, the entire motion is dominated by inertial waves, which reduce the growth of Lagrangian displacement in agreement with a Brownian-type law, but do not stop it since no plateau appears. In the latter case, anisotropy only consists of a ratio 2 of vertical to horizontal diffusivities.

If one now considers organised structures, for which nonlinear dynamics are needed, they are not robust enough to significantly affect Lagrangian diffusivity (for instance one would expect a plateau for horizontal rms displacement, if the "cigars" were really strong vertical vortices in rotating flows). This competition between linear effects of anisotropic dispersive waves and nonlinear effects of structure formation could be different when studying other Lagrangian statistics such as two-particle dispersion. This question shall be adressed in a future work by computations of two-particle dispersion by KS and DNS, as well as by the development of a two-particle two-time RDT velocity correlation model.

In addition to the above mentioned "paradox", another important result of our study is the strong angular dependence of anisotropic spectra obtained from high resolution DNS, down to the smallest scales. The specific anisotropization mechanism corresponds to a nonlinear energy drain, confirmed by DNS to be angle-to-angle, rather than shell-to-shell as for isotropic cascade. This dynamics agrees with results from existing anisotropic statistical models [START_REF] Godeferd | Detailed investigation of energy transfers in homogeneous stratified turbulence[END_REF][START_REF] Cambon | Turbulence and vortex structures in rotating and stratified flows[END_REF][START_REF] Godeferd | Statistical modelling and direct numerical simulations of decaying stably-stratified turbulence: Part 2: Large and small scales anisotropy[END_REF]. In addition to this confirmation, our new results are that angular dependence is essentially restricted to the toroidal mode in stably stratified turbulence, and that the quasi-isotropic shape of poloidal energy scales like k -2 . The physical interpretation of both results is an open question, to be adressed with the help of statistical models. Finally, our work has raised the issue of how to produce deterministic initial data for the computations, that can be physically relevant as sources of motion, especially for wave propagation [16].

This system of differential equations can be rewritten as a matrix equation where, by standard linear algebra, ξ ǫ = v • N ǫ * , which can easily be obtained by orthogonality of N ǫ . From the diagonalized system of equations, one can deduce the time dependent term, so getting a general solution for v which has the following form:

v(k, t) = ǫ=0,±1
N ǫ e -iǫσt N -ǫ • v(k, 0) (A.12)

Figure 1 .

 1 Figure 1. Enstrophy isosurfaces from freely decaying DNS with a resolution of 256 3 , at 20% of the maximum enstrophy in the box, colored by the vertical velocity. The Taylor microscale based Reynolds number is around 50. Adimensional parameters: (a) α = 0.1, Fr λ = 0.1, (b) α = 1, Fr λ = 0.1 and (c) α = 10, Ro λ = 0.05, where α = f /N , Fr = u/(N λ) and Ro = u/(f L).

Figure 3 .

 3 Figure 3. (a) The unit vectors in the local Craya-Herring frame of reference. (b) The decomposition of a Fourier velocity vector in its poloidal and toroidal part compared to the decomposition into its wave and vortex part.

Figure 4 .

 4 Figure 4. Energy evolution of DNS runs B2, C2, D2: (a) total energy; (b) toroidal part; (c) poloidal part; (d) potential energy.

Figure 5 .

 5 Figure 5. DNS time evolution of the directional correlation length scales: (a) velocity: dashed lines: L x xx ; solid lines: L z xx (b) vorticity: dashed lines: L y ωxωx ; solid lines: L z ωzωz

Figure 7 . 1 0 ( 1 -

 711 Figure 7. (a) Linear prediction (computed from ∝ 1 0 (1 -cos 2 θ) 1-cos σt σ 2 d(cos θ)) as defined in section 2, (b) random Fourier modes model (runs F to L), and (c) nonlinear simulation (runs A1 to E1) of one-particle vertical dispersion for different α = f /N . The two top curves (blue and gray) are scaled with f instead of N as it is the dominant parameter. Vertical range for all three plots is the same.

Figure 8 . 1 0

 81 Figure 8. (a) Linear prediction (computed as ∝ 1 0 (1 + cos 2 θ) 1-cos σt σ 2 d(cos θ), (b) random Fourier modes model, and (c) nonlinear simulation of one-particle horizontal dispersion for different α = f /N . Vertical range for all three plots is the same.

Figure 9 .

 9 Figure 9. (a) Integrating surface for angular dependent spectra. (b) Direction of integration of an angular spectrum along k.

  (b). This representation of the flow field allows us to analyze anisotropy in the flow at all length scales and is shown for the toroidal mode in figures 10(a), (b) and (c) and for the poloidal mode in figures 10 (d), (e) and (f) for DNS runs B2, C2 and D2.

Figure 10 .

 10 Figure 10. Poloidal (top) and toroidal (bottom) angular spectra for DNS runs B2 (left), C2, (middle) and D2 (right). Black curves show an isotropic spherical spectrum for comparison. Vertical range for all six plots is the same.

  σ s = N sin θ and σ r = f cos θ. The eigenvalues of the matrix are 0 and ±iσ, where σ = σ 2 s + σ 2 r and the normalized eigenvectors are This basis of eigenmodes can now be used to express v as v = ǫ=0,±1 ξ ǫ N ǫ , (A.11)

Table 2 .

 2 Sets of parameters of the KS runs presented.

Acknowledgments

High resolution DNS computations have been performed thanks to computational time provided in equal part by IDRIS (computing center of the Centre National de la Recherche Scientifique) and CCRT (computing center of the Commissariat à l' Énergie Atomique).

Appendix A. Craya-Herring frame

Divergence-free flows can be characterized as a two-component flow by defining a local (as a function of k) orthonormal reference frame (Craya-Herring frame, see figure 3(a)) with the first direction horizontal, the third direction along k and the last direction orthogonal to the first two

This decomposition is essentially the division of the toroidal (mode 1) and poloidal (mode 2) part, so the velocity vector with two components can then be presented as

with its toroidal û(1) (k) and poloidal û(2) (k) modes. Following e.g. [START_REF] Cambon | Turbulence and vortex structures in rotating and stratified flows[END_REF] we introduce a third component û(3) as an imaginary component representing the square root of the potential energy

and so construct a hermitian vector field

represented in the Craya-Herring frame of reference. This leads to an expression for the total energy density of the fluid

In components of the new reference frame the linearized equations ( 1) and ( 2) take the following form

∂ t û(2) + f cos θû (1) + iN sin θû (3) = 0 (A.7)

-iN ∂ t û(3) + N 2 sin θû (2) = 0 (A.8)