A M A F F I O L I

P A D A V I D S O N

Dynamics of stratified turbulence decaying from a high buoyancy Reynolds number

We present direct numerical simulations of unforced stratified turbulence with the objective of testing the strongly stratified turbulence theory. According to this theory the characteristic vertical scale of the turbulence is given by ℓ v ∼ u h /N , where u h is the horizontal velocity scale and N the Brunt-Väisälä frequency. Combined with the hypothesis of the energy dissipation rate scaling as ǫ ∼ u 3 h /ℓ h , this theory predicts inertial range scalings for the horizontal spectrum of horizontal kinetic energy and of potential energy, according to E(k h ) ∝ k -5/3 h

. We begin by presenting a scaling analysis of the horizontal vorticity equation from which we recover the result regarding the vertical scale, ℓ v ∼ u h /N , highlighting in the process the important dynamical role of large-scale vertical shear of horizontal velocity. We then present the results from decaying DNS which show a good agreement with aspects of the theory. In particular, the vertical Froude number is found to reach a constant plateau in time, of the form F r v = u h /(N ℓ v ) = C with C = O(1) in all the runs. The derivation of the dissipation scaling ǫ ∼ u 3 h /ℓ h at low Reynolds number in the context of decaying stratified turbulence, highlights that the same scaling holds at high R = ReFr 2 h ≫ 1 as well as at low R ≪ 1, which is known (see [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF] but not sufficiently emphasized in recent literature. We find evidence in our DNS of the dissipation scaling holding at R = O(1), which we interpret as being in the viscous regime. We also find ǫ k ∼ u 3 h /ℓ h , ǫ p ∼ u 3 h /ℓ h (with ǫ = ǫ k + ǫ p ), in our high-resolution run at earlier times corresponding to R = O(10), which is in the transition between the strongly stratified and the viscous regimes. The horizontal spectra of horizontal kinetic energy collapses in time using the scaling E h (k h ) = C 1 ǫ 2/3 k k -5/3 h and the horizontal potential energy spectra is well described by E p (k h ) = C 2 ǫ p ǫ -1/3 k k -5/3 h . The presence of an inertial range in the horizontal direction is confirmed by the constancy of the energy flux spectra over narrow k h -ranges. However, the vertical energy spectra are found to differ significantly from the expected E h (k v ) ∼ N 2 k -3 v scaling, showing that F r v is not of order unity on a scale-by-scale basis, thus providing motivation for further investigation of the vertical structure of stratified turbulence.

Key words: Authors should not enter keywords on the manuscript, as these must be chosen by the author during the online submission process and will then be added during the typesetting process (see http://journals.cambridge.org/data/relatedlink/jfmkeywords.pdf for the full list)

Introduction

The atmosphere and oceans present strong stable density stratification over significant portions of their height and depth. These systems are also affected by the Earth's rotation, which however has a much smaller influence than the stratification. In terms of the flow dynamics, as highlighted by [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF], "when turbulence is generated in the presence of stable density stratification, it often occurs locally and without a continuous source of energy" as a result of local Kelvin-Helmholtz (KH) instability or breaking of internal gravity waves (IGWs). The study of decaying stratified turbulence is therefore of most interest. Comparatively few previous studies of decaying stratified turbulence exist, both in terms of experiments [START_REF] Fincham | Energy dissipation and vortex structure in freely decaying stratified grid turbulence[END_REF][START_REF] Praud | Decaying grid turbulence in a strongly stratified fluid[END_REF] and numerical simulations [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Bartello | Sensitivity of stratified turbulence to the buoyancy Reynolds number[END_REF]. On the other hand, forced stratified turbulence has received considerable attention particularly in the case of numerical studies where the achievement of statistical stationarity is desirable (see e.g. [START_REF] Smith | Generation of slow large scales in forced rotating stratified turbulence[END_REF][START_REF] Laval | Forced stratified turbulence: successive transitions with Reynolds number[END_REF][START_REF] Waite | Stratified turbulence dominated by vortical motion[END_REF][START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF][START_REF] Waite | Stratified turbulence at the buoyancy scale[END_REF][START_REF] Almalkie | Kinetic energy dynamics in forced, homogeneous, and axisymmetric stably stratified turbulence[END_REF][START_REF] Augier | Stratified turbulence forced with columnar dipoles: numerical study[END_REF]. Most of the recent numerical studies of forced stratified turbulence employ purely horizontal forcing to avoid direct excitation of IGWs, and force modes with k z = 0 to avoid selecting a vertical scale, two conditions which are hardly met in nature (neither KH instability or IGW breaking lead to this kind of forcing). Moreover, the forcing directly excites a certain horizontal scale, which will therefore not be dynamically selected.

The strongly stratified turbulence theory has been developed by a number of researchers in the first decade of this century. The first step has been the work of [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF] who postulate that the vertical lengthscale of the turbulence should be ℓ v = u h /N and find a new set of reduced order equations. In this paper, the vertical spectrum of horizontal kinetic energy is predicted to be E h (k v ) ∼ N 2 k -3 v ; this form of the vertical spectrum had classically been found by [START_REF] Dewan | Saturation and the "universal" spectrum for vertical profiles of horizontal scalar winds in the atmosphere[END_REF] who interpreted it as being due to saturation of internal gravity waves. In the following years, [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF] makes the hypothesis that horizontal kinetic energy and potential energy are approximately equipartitioned and uses ǫ k ∼ ǫ p ∼ u 3 h /ℓ h (ǫ k is the viscous dissipation while ǫ p is the dissipation due to diffusion of buoyancy) in his analysis to find that the horizontal spectrum of horizontal kinetic energy should scale as E h (k h ) ∼ ǫ 2/3 k k -5/3 h and the horizontal spectrum of potential energy as E p (k h ) ∼ (ǫ p /ǫ 1/3 k )k -5/3 h . Atmospheric spectra have been explained using the strongly stratified turbulence theory by [START_REF] Riley | Stratified turbulence: a possible interpretation of some geophysical turbulence measurements[END_REF]. The strongly stratified turbulence regime is subject to the conditions Fr h = u h /(N ℓ h ) ≪ 1 and R = Re Fr 2 h ≫ 1 (Re = u h ℓ h /ν is the Reynolds number based on horizontal scales). The latter condition is equivalent to the buoyancy Reynolds number condition, Re b = ǫ k /(νN 2 ) ≫ 1, as discussed by [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF].

At this point it is necessary to consider the recent literature behind the condition R = Re Fr 2 h ≫ 1 and the difference between R and the buoyancy Reynolds number Re b = ǫ k /(νN 2 ). The importance of the buoyancy Reynolds number in geophysical turbulence and in stratified mixing layers has been emphasized by [START_REF] Smyth | Length scales of turbulence in stably stratified mixing layers[END_REF]. Taking a more theoretical standpoint, [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF] consider the importance of diffusive effects on strongly stratified turbulence and identify the vertical diffusion component as a major contributor to viscous effects, but only if the condition R ≫ 1 is not met. Following this paper, an important contribution was brought by [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF] who identify the local gradient Richardson number Ri g as an important non-dimensional parameter in stratified turbulence. In particular they find that the local Richardson number has to be small in some areas of the flow in order for turbulence to develop, with the formation of KH instabilities. In their article [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF] also suggest that Ri g ∼ 1/R describes well their simulations and a wider range of stratified turbulence problems. The route to the development of strongly stratified turbulence starting from Taylor-Green vortices has been further pursued by [START_REF] Hebert | Predicting turbulence in flows with strong stable stratification[END_REF] who confirm that, in the horizontal plane of maximum vertical shear, Ri g ≈ 1/R. They also find that, in this horizontal plane, R ∼ Re b , leading to Ri g ∼ νN 2 /ǫ k . Moreover, the buoyancy Reynolds number Re b has more recently emerged (see [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF][START_REF] Almalkie | Kinetic energy dynamics in forced, homogeneous, and axisymmetric stably stratified turbulence[END_REF][START_REF] Bartello | Sensitivity of stratified turbulence to the buoyancy Reynolds number[END_REF] as another measure of the relative importance of viscous effects in strongly stratified turbulence. Again the condition Re b ≫ 1 is necessary for strongly stratified turbulence without viscous effects. However, the two non-dimensional parameters R and Re b are only equivalent when the dissipation scaling ǫ k ∼ u 3 h /ℓ h holds. This brings us back to the main focus of the paper, the conditions ǫ k ∼ u 3 h /ℓ h and ǫ p ∼ u 3 h /ℓ h on which much of the strongly stratified turbulence theory resides (see [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]. These conditions stem from the evolution equation of total energy in purely decaying stratified turbulence, dE tot /dt = -ǫ where ǫ = ǫ k + ǫ p and E tot = E k + E p , E k is the kinetic energy and E p the potential energy, and all quantities are volume averages. If, by virtue of the strong anisotropy in the velocity components and lengthscales of stratified turbulence, we set E k ∼ u 2 h and choose a typical energy cascade timescale to be τ ∼ ℓ h /u h , the dissipation scaling is recovered after invoking equipartition, E k ∼ E p , and using the convective timescale τ also for the dissipation of buoyancy and potential energy. This dissipation scaling has been used by many authors (see [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF][START_REF] Riley | Stratified turbulence: a possible interpretation of some geophysical turbulence measurements[END_REF] to estimate the horizontal lengthscale as ℓ h ∼ u 3 h /ǫ. The only evidence of this scaling holding in decaying stratified turbulence comes from [START_REF] Hebert | Predicting turbulence in flows with strong stable stratification[END_REF], but this study is limited to a horizontal plane where very large vertical shear is present and does not consider the entire simulation domain. The equivalent form ǫ ∼ u 3 /ℓ is well established in homogeneous, isotropic turbulence, with extensive experimental and numerical evidence [START_REF] Yeung | Universality of the Kolmogorov constant in numerical simulations of turbulence[END_REF][START_REF] Davidson | Turbulence. An introduction for scientists and engineers[END_REF]). The dissipation scaling ǫ ∼ u 3 h /ℓ h was predicted to be valid also at low R by [START_REF] Godoy-Diana | Vertical length scale selection for pancake vortices in strongly stratified viscous fluids[END_REF]; [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF], a fact that is sometimes overshadowed in the strongly stratified turbulence literature.

Concerning spectra in the horizontal direction, [START_REF] Bartello | Sensitivity of stratified turbulence to the buoyancy Reynolds number[END_REF] find that the spectra are independent of stratification and of buoyancy Reynolds number only if Re b > 10. They then find a k -5/3 h -scaling of the horizontal spectra for their most stratified case satisfying this condition on Re b . It is important to note that their definition of buoyancy Reynolds number is Re (ℓ v /ℓ h ) 2 , which is different to the above definition based on ǫ k .

The aim of this paper is to understand the applicability of the four predictions of the strongly stratified turbulence theory:

F r v ∼ 1, ǫ ∼ u 3 h /ℓ h , horizontal spectra scaling as E h (k h ) ∼ ǫ 2/3 k k -5/3 h and E p (k h ) ∼ (ǫ p /ǫ 1/3 k )k -5/3 h and vertical spectra scaling as E h (k v ) ∼ N 2 k -3 v
. Our first step is to consider the full vorticity equation and provide a scaling analysis in the limit of high R, as presented in §2. We then continue with presenting direct numerical simulations (DNS) of stratified turbulence, starting with the numerical methodology and description of runs in §3. The main results are presented in §4, where we consider the validity of Fr v ∼ 1 in our runs, the dissipation scaling for ǫ k and ǫ p and the horizontal spectra of horizontal kinetic energy and potential energy. This section ends with a look at the vertical spectra of the stratified turbulent flows that have been simulated. A discussion section follows in §5 in which the anisotropy of the RMS velocity components is considered at high R and finally concluding remarks are given in §6.

2. Scaling analysis of the vorticity equation in the limit of high R

We start by presenting a scaling analysis of the vorticity equation. This has the advantage of not requiring any assumption regarding the scaling of the pressure term, as this term is eliminated when we take the curl of the momentum equation. In this analysis we assume from the outset that ℓ v ≪ ℓ h so that layers are present in the flow, whereas [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF] take the aspect ratio α = ℓ v /ℓ h as a free parameter and take Fr h ≪ 1 (a condition we do not require a priori ).

As mentioned above, our starting point is the vorticity equation that is obtained by taking the curl of the Boussinesq momentum equation. We divide the vorticity equation into two equations, one valid for the horizontal component of vorticity ω h = ω x e x +ω y e y , and the other valid for the vertical component of vorticity ω z :

∂ω h ∂t + u • ∇ω h = ω • ∇u h -e z × ∇b + ν∇ 2 ω h (2.1) ∂ω z ∂t + u • ∇ω z = ω • ∇u z + ν∇ 2 ω z (2.2) ∇ h • u h + ∂u z ∂z = 0 (2.3) ∂b ∂t + u • ∇b = -N 2 u z (2.4)
where u is the velocity, ω = ∇ × u the vorticity, b = -ρ ′ g/ρ 0 the buoyancy (with g = -ge z being the gravitational acceleration), ν the kinematic viscosity and N = -(g/ρ 0 )dρ/dz the Brunt-Väisälä frequency that quantifies the density stratification. The total density is thus ρ = ρ 0 + ρ+ρ ′ , the first term being a reference density, the second being the linear background density profile and the third being the perturbation density from the linear profile. We have added mass conservation and the buoyancy equation without diffusive effects to this set as they are necessary for the analysis.

We note in passing that the vertical vorticity equation simplifies to the advection equation Dω z /Dt = 0 if we neglect viscous effects and take the vortex stretching to be small compared to horizontal advection of ω z . This in turn is the form of the potential vorticity (PV) conservation equation for the linear PV component if we neglect viscous and diffusive effects [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF].

We now focus on the integral scale contributions to (2.1)-(2.4), and make the main assumption of the scaling and set ℓ v ≪ ℓ h ; this assumption is based on the observation in laboratory experiments, numerical simulations and indeed in the atmosphere, that turbulence in a stratified fluid brings about the formation of pancake eddies or layers with vertical scales much smaller than their horizontal scales. Now, from the continuity equation we see that u v ∼ u h ℓ v /ℓ h . Therefore the scaling for the horizontal vorticity is

ω h ∼ u h /ℓ v + u v /ℓ h ∼ u h /ℓ v (1 + (ℓ v /ℓ h ) 2 ) ≈ u h /ℓ v .
The scaling for the vertical vorticity is ω z ∼ u h /ℓ h . We take time to scale as a horizontal advective timescale, t ∼ ℓ h /u h , which amounts to filtering out any internal gravity waves (unless their timescale is very low and of the order of ℓ h /u h ). Finally, the scaling for the buoyancy b is taken as b ∼ N 2 ℓ v ; this follows directly from (2.4) using continuity, u v /ℓ v ∼ u h /ℓ h .

We now focus on the horizontal vorticity equation and provide scalings for every term in the equation:

Advection : u • ∇ω h ∼ u h u h ℓ h ℓ v + u v u h ℓ 2 v = u 2 h ℓ h ℓ v (2.5) Vortex stretching : ω • ∇u h ∼ u h ℓ v u h ℓ h + u h ℓ h u h ℓ v = u 2 h ℓ h ℓ v (2.6) Buoyancy term : ∇b × e z ∼ N 2 ℓ v ℓ h (2.7) Viscous diffusion : ν∇ 2 ω h ∼ ν u h ℓ 2 h ℓ v + u h ℓ 3 v = νu h ℓ 3 v 1 + ℓ v ℓ h 2 ≈ νu h ℓ 3 v (2.8)
It is clear that the first two terms -advection and vortex stretching -have the same order of magnitude under the assumptions made. We are interested in the form the equations assume at high Reynolds number, in the inviscid limit. Therefore we form a stratified Reynolds number by calculating the ratio of the advective term and the viscous term and set this term to be much greater than unity, so that the viscous term can be neglected:

advection viscous diffusion ∼ u 2 h /ℓ h ℓ v νu h /ℓ 3 v = u h ℓ 2 v νℓ h ≫ 1 (2.9)
Finally, because we are interested in the dynamics strongly influenced by stratification, we equate the buoyancy term with the advection term. This leads to:

N 2 ℓ v ℓ h ∼ u 2 h ℓ h ℓ v (2.10) from which Fr v = u h N ℓ v ∼ 1 (2.11)
We have thus recovered the F r v = O(1) scaling proposed by many authors [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF][START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF]. It follows from this result that Fr h = u h /N ℓ h ∼ ℓ v /ℓ h ≪ 1 using our assumption of layered flow. So indeed the horizontal Froude number has to be low for the scalings to hold. If we further substitute the vertical scale obtained in this way, ℓ v ∼ u h /N , into the expression of (2.9) we have:

advection viscous diffusion ∼ u h (u h /N ) 2 νℓ h = u h ℓ h ν u h N ℓ h 2 = Re Fr 2 h = R (2.12)
Evidently we require R = Re Fr 2 h ≫ 1 to ensure that the viscous terms may be neglected.

Numerical methodology

We now discuss direct numerical simulations of decaying stratified turbulence. The numerical methodology we employ consists in solving the Boussinesq set of equations directly for a linearly stratified fluid. The Boussinesq set of equations can be written as:

∇ • u = 0, (3.1) ∂u ∂t + u • ∇u = - 1 ρ 0 ∇p + b e z + ν∇ 2 u, (3.2) ∂b ∂t + u • ∇b = -N 2 u z + D∇ 2 b, (3.3)
where p is the perturbation pressure from the hydrostatic pressure distribution, D is the buoyancy diffusivity and the remaining symbols have already been defined. Equations (3.1)-( 3.3) are solved using a pseudo-spectral method based on Rogallo's algorithm [START_REF] Rogallo | Numerical experiments in homogeneous turbulence[END_REF]. Time advancement is carried out using a second-order Runge-Kutta predictor-corrector integration scheme, while the viscous and diffusive terms are integrated exactly by using suitable integrating factors. De-aliasing of the non-linear terms is performed using a combination of truncation and phase-shifting. The simulations are performed in rectangular domains with L x = L y > L z , with vertical to horizontal length ratios of 1/6, with the exception of one run at higher resolution where the box aspect ratio is 1/4. The resolution of the simulations is chosen so that the grid spacing is isotropic, i.e. ∆x = ∆y = ∆z. We would like to capture dynamics associated with both potential vorticity and with internal gravity waves. Previous decaying simulations with Taylor-Green vortices (see [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF] were skewed towards vortical motions. However, internal gravity waves have been found to be important in creating organized structures as shown recently by [START_REF] Maffioli | The evolution of a stratified turbulent cloud[END_REF] for localized patches of turbulence. The runs are therefore initialized with isotropic random-phase velocity fields satisfying:

E(k) = U 2 k 4 exp(-k 2 /k 2 p ) (3.4)
where E(k) = 1/2 < û * • û > is the 3D energy spectrum and the < . . . > averaging is over spherical shells at radius k. For all runs we take U 2 = 0.22 and k p = 5 giving a narrow spectrum peaking at √ 2k p ≈ 7. We choose the viscosity, and hence the initial Reynolds number, in order to ensure that our fields are well resolved. We use the condition k max η > 1 (k max is the maximum resolved wavenumber after truncation of the highest modes and η = (ν 3 /ǫ k ) 1/4 is the Kolmogorov lengthscale) as the criterion for good resolution of the dissipation scales (as delineated by [START_REF] Eswaran | An examination of forcing in direct numerical simulations of turbulence[END_REF][START_REF] De Bruyn Kops | Direct numerical simulation of laboratory experiments in isotropic turbulence[END_REF], and this condition is met at all times throughout the evolution of the turbulent field. This is a severe requirement considering that η, after an initial drop as more small-scale turbulence is generated, steadily grows throughout a decaying run. Some trial-and-error was required to obtain the desired value for ν as the minimum in η was not known a priori. A typical plot of k max η throughout run R80F0.2 (see table 1) is given in figure 1; by the end of the run k max η = 2.9.

The various runs are listed in table 1 together with non-dimensional parameters at the peak in kinetic energy dissipation when small-scales have formed. From this moment onwards the simulation approaches "real" decaying stratified turbulence and hence values for the important parameters are given here and not at the start of the simulation. The horizontal Froude number for all cases is Fr h,ss = O(1) so that initially stratification is not dominant but it will become so as the turbulence decays further. This is the situation that is thought to occur in the atmosphere and oceans with vigorous turbulent outbursts that then decay until Fr h ≪ 1 and stratification dominates as discussed by [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF]. In order to obtain Fr h,ss = O(1) the initial Froude number was also of order one, which justifies the use of isotropic initial conditions. We use a single grid resolution for the initial four runs -denoted R1.8kF1.8, R330F0.6, R80F0.2, R26F0.1 following their "stratified Reynolds number" and horizontal Froude number when small scales have formed. A higher resolution run is also performed, R49F0.1. All runs have R ss ≫ 1 and, importantly, R > 10; this is the condition to have strongly stratified turbulence, together of course with Fr h ≪ 1, which will be attained as turbulence decays. Also the buoyancy Reynolds number is high when small-scales have formed, and all runs have Re b,ss ≫ 1 except run R26F0.1, where Re b,ss = 4.2. All cases are run up to when R is of order one. Hence we cover the entire evolution of strongly stratified turbulence. The only exception is the high-resolution run R49F0.1, which was run for a smaller amount of physical time reaching a final R end = 6.3.

The buoyancy field is initialized with b = 0 everywhere. This results in high-frequency internal gravity waves with ̟ = N being generated at t = 0 that affect the u z and b fields as the buoyancy field reacts to having large vertical velocities initially. The horizontal lengthscale ℓ h presented in what follows is an integral lengthscale obtained from the transverse velocity correlations:

ℓ x = 1 < u 2 y > r0 0 < u y (x)u y (x + re x ) > dr, (3.5) ℓ y = 1 < u 2 x > r0 0 < u x (x)u x (x + re y ) > dr, (3.6) ℓ h = 1 2 (ℓ x + ℓ y ), (3.7)
where r 0 is the horizontal separation at which the velocity correlation first becomes equal to zero and < . . . > now represents a volume average throughout the computational domain. The vertical lengthscale ℓ v is obtained in a similar fashion using the transverse correlations of horizontal velocity:

ℓ (x) z = 1 < u 2 x > r0 0 < u x (x)u x (x + re z ) > dr, (3.8) ℓ (y) z = 1 < u 2 y > r0 0 < u y (x)u y (x + re z ) > dr, (3.9) ℓ v = 1 2 (ℓ (x) z + ℓ (y) z ).
(3.10) The spectra presented throughout the paper are one-dimensional spectra found by summation over "planar shells":

E (x) h (k h ) = |kx|=k h ky kz 1 2 (û x û * x + ûy û * y ), (3.11) E (y) h (k h ) = kx |ky|=k h kz 1 2 (û x û * x + ûy û * y ), (3.12) E h (k h ) = 1 2 (E (x) h (k h ) + E (y) h (k h )), (3.13) E h (k v ) = kx ky |kz|=kv 1 2 (û x û * x + ûy û * y ), (3.14)
where * denotes complex conjugation. In these definitions we use only horizontal velocity components as most of the strongly stratified theory has been developed considering only horizontal velocities and taking E k ∼ u 2 h as done by [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]. The potential energy spectra are computed in an analogous fashion by summation of b * b /(2N 2 ).

Results

Visualization of flow field and emergence of strongly stratified turbulence

The turbulence simulations were initialized with isotropic velocity fields and zero density perturbations. When the stratification is abruptly switched on at t = 0 the turbulence takes some time to adjust to the vertical restoring force caused by the linear background density profile. The global effect is that < u 2 z > slowly drops while buoyancy < b 2 > is generated. As is evident from figure 2, this process takes the form both of high-frequency exchanges of kinetic and potential energy (the oscillations on E p and 1/2 < u 2 z > correspond to ̟ = N internal gravity waves as observed by [START_REF] Maffioli | The evolution of a stratified turbulent cloud[END_REF], in their simulations of a turbulent cloud with stratification) and of a slower non-linear process on a timescale similar to the initial turnover time τ 0 .

The effect of the stratification on the turbulence structure is apparent if we turn to visualisations of vertical and horizontal slices through the parallelepiped-shaped domains (see figures 3 and 4). At the peak in dissipation the buoyancy field is highly-turbulent and approximately isotropic (see figure 3(a)), but as time progresses and the turbulence dies Froude numbers at the same non-dimensional time, we see from figure 5 that the layering is more apparent as Fr h,ss is lowered (the instantaneous Fr h in the four images is decreasing from figure 5(a) to (d)). In addition, we have the interesting feature that the high-Fr h,ss runs (R1.8kF1.8 and R330F0.6) still appear vigorously turbulent at this stage, while the low-Fr h,ss cases (R80F0.2 and R26F0.1) appear to be in a viscous regime. This observation can be explained if we consider the value of R at this time for the different runs: R = 8.9 for case R1.8kF1.8, R = 6.3 for case R330F0.6, R = 3.1 for case R80F0.2 and R = 1.6 for case R26F0.1. Hence runs R1.8kF1.8 and R330F0.6 still marginally satisfy R ≫ 1 at this time and still constitute strongly stratified turbulence, while runs R80F0.2 and R26F0.1 have R ∼ 1 and they are in a transition phase between the highly-turbulent and the viscous regimes of stratified turbulence.

Having obtained a visual appreciation of the stratified turbulence under investigation, we now characterize it quantitatively, looking in particular at the energy dissipation rate, at correlation lengthscales and energy spectra.

Vertical Froude number

One of the main results of our scaling analysis in §2 is the recovery of the vertical scale relation ℓ v ∼ u h /N as found in the original scaling analysis of the momentum equation by [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF]. Of course, this result holds only for high buoyancy Reynolds numbers, or for high values of the parameter R = ReFr 2 h . However, our runs were intentionally stopped when R falls to order unity allowing us to check if ℓ v scales as the buoyancy scale. When R = O(1) we expect there to be a transition from the scaling ℓ v ∼ u h /N to a viscous scaling. This means that, if we plot Fr v as a function of time, we expect it to reach a constant value while R ≫ 1. Considering our visualisations in Figure 5, the stratified turbulence in our DNS runs appears to change character and be less vigorous at values of R < 6.3. We therefore substitute the condition R > 7 for the more stringent (but not quantifiable) condition R ≫ 1, and plot the vertical Froude number as a function of time only for times at which this condition is met. This is done in Figure 6 and it is clear that Fr v = const for all runs and the aysmptotic value is similar in all cases, at Fr v = 0.34-0.37. We also show the full evolution of Fr v (inset panel in figure 6); this shows a slow variation of Fr v at late times (especially true for the low-Fr 0 case F0.16) as the stratified turbulence is more affected by vertical viscous diffusion. This is accordance with a transition at R = O(1) between strongly stratified turbulence and a viscously-dominated regime.

All in all the simulations provide convincing evidence that, after an initial transition from random-phase velocity fields to fully-developed stratified turbulence, the vertical Froude number is of order one. This is one of the central claims of the strongly-stratified turbulence theory and, to our knowledge, it had not been verified in the context of freely decaying simulations.

Dissipation scaling

It has been suggested by numerous authors that in stratified turbulence both the viscous dissipation ǫ k and the diffusive dissipation ǫ p = (D/N 2 ) ∂b/∂x i ∂b/∂x i scale with horizontal length and velocity scales:

ǫ k = A k u 3 h ℓ h , A k = const; ǫ p = A p u 3 h ℓ h , A p = const. (4.1)
Both constants A k and A p are expected to be finite and of order one.

In figure 7(a) we show the total energy dissipation rate ǫ = ǫ k +ǫ p averaged throughout the simulation domain and normalized by u 3 h /ℓ h as a function of time for the four runs at the lower resolution. According to the relations (4.1) the quantity ǫℓ h /u 3 h should reach a constant value of order one. We see that this is indeed the case at late times, however both R and the buoyancy Reynolds number are close to 1 at these times and therefore we expect viscosity to be playing a role here. To clarify further that our four lower resolution simulations are not in the strongly stratified turbulence regime when the curves flatten out, we apply the same bound used in the previous section to test F r v ∼ 1. The portion of the four curves presenting a close to constant non-dimensional dissipation are all clearly beyond the bound R = 7 (which corresponds to a buoyancy Reynolds number Re b ≈ 2).

To understand this apparently contradictory result, we consider the opposite limit to the one treated in the above scaling analysis in §2, namely the limit of strong viscous dissipation R ≪ 1. In this limit, experiments on stratified turbulence (Fincham et al. 1996) show that vertical diffusion dominates the viscous dissipation and that the flow is still highly anisotropic and composed of pancake structures with ℓ v ≪ ℓ h . Given that this flow regime is dominated by the action of viscosity in the vertical direction, we take the vertical lengthscale to obey a viscous diffusion process ℓ v ∼ √ νt, similarly to Godoy-Diana et al. (2004). We also make the assumption that the total dissipation is well described by ǫ ∼ u 3 h /ℓ h , which is what has to be later verified in the scalings. Hence the sum of kinetic and potential energy obeys the evolution equation

dE tot /dt = -Au 3 h /ℓ h . Now, invoking equipartition, E k ∼ E p , with E k ∼ u 2
h and E p ∼ b 2 /N 2 , we arrive at the evolution equation for the horizontal kinetic energy

du 2 h dt = -A k u 3 h ℓ h , u 2 h ∼ t -n , (4.2)
where the power-law behaviour is expected as the flow is still turbulent, despite being in the viscous regime. Hence we have

-n u 2 h t = -A k u 3 h ℓ h (4.3) from which u h /ℓ h = n/(A k t)
ǫ k ℓ h u 3 h ǫpℓ h u 3 h R > 7 R < 7 (b)
let ǫ k ≈ ν(∂u h /∂z) 2 ∼ νu 2 h /ℓ 2 v ∼ u 3
h /ℓ h and we have recovered the dissipation scaling for the viscous dissipation. Considering that we are interested in the case Pr = ν/D = 1, which is the case in the DNS runs, the same analysis clearly also leads to ǫ p ∼ u 3 h /ℓ h and therefore the total dissipation ǫ presents the same scaling. In conclusion we can say with some confidence that this re-derivation and figure 7(a) further strengthen the scaling proposed for the viscous regime (see [START_REF] Godoy-Diana | Vertical length scale selection for pancake vortices in strongly stratified viscous fluids[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF].

We now consider the high-resolution run R49F0.1. In figure 7(b) the full evolution of this run is presented. It is relatively clear that after an initial transient associated with the stratified turbulence becoming fully-developed, the curve flattens out at around t/τ 0 = 11.5. Again, the boundary demarcating the region in which R > 7 from the later times at which R < 7 is given, and it occurs close to the end of the simulation. There is a period of about 3τ 0 ahead of this boundary during which both ǫ k and ǫ p scale with horizontal quantities alone and, when normalized, reach a constant plateau. Despite the fact that the dissipation scaling is presented earlier than the low-resolution runs, at this time we are already in the transition between the strongly stratified regime and the viscous regime. At t/τ 0 = 11.5, R = 10 but the buoyancy Reynolds number Re b = 2.1 as can be seen in figure 8 showing the concurrent evolution of R and Re b for runs R80F0.2 and R49F0.1. This value of Re b is too low as can be seen more clearly when calculating the ratio of dissipation due to vertical gradients of horizontal velocity to the total dissipation, ǫ k,hv /ǫ k = 0.69 at this time. The turbulent motion is starting to be dominated by vertical shear and we are approaching the viscous regime. So we have been unable to verify the dissipation scaling in the strongly stratified regime. Our partial result do however support the claim that the dissipation scaling ǫ ∼ u 3 h /ℓ h is a robust result, holding across the various regimes of stratified turbulence.

Horizontal spectra of horizontal kinetic energy

According to [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF] the horizontal spectrum of horizontal kinetic energy (which is very similar to the spectrum of full kinetic energy, given that the vertical component of velocity is small) should scale as

E h (k h ) ∼ ǫ 2/3 k k -5/3 h
for strongly stratified turbulence. Let us now see if this is the case for our decaying runs. Figure 9 shows one-dimensional horizontal spectra of horizontal kinetic energy compensated by multiplication by ǫ -2/3 k k 5/3 h for the four DNS runs. These compensated spectra should present a plateau at a constant value of Ẽh (k h ) = ǫ -2/3 k k 5/3 h E h (k h ) over the stratified inertial range if the Lindborg scaling applies. We have chosen to plot the spectra using a linear axis for the y-axis as often log-log plots show apparent plateaus where there still is significant variation in the compensated spectra. We can see that the compensated spectra present about half a decade of constancy for runs F1.3, F0.6 and F0.3, starting from a wavenumber k h ≈ 10. We have chosen the time instants in figure 9 as they present the clearest stratified inertial range and these times correspond to R ≈ 15 for the three runs with a clear plateau. This is in line with previous knowledge [START_REF] Bartello | Sensitivity of stratified turbulence to the buoyancy Reynolds number[END_REF] that Re (ℓ v /ℓ h ) 2 = 10 appears to be the lower bound for the strongly stratified horizontal spectra to be observed. In our DNS simulations this condition actually would translate to Re F r 2 h /F r 2 v ≈ 8Re F r 2 h = 8R > 10 or R > 1.3 using the value Fr v = 0.35 reached by all our runs. We find in our runs that R ≈ 15 is needed for the stratified inertial range to be observed. The extent of the stratified inertial range should be from the horizontal integral scale to the Ozmidov scale, and this range is determined solely by the Froude number since k oz /k ℓ h ∼ Fr -3/2 h [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]. Hence the lower the Froude number the more wavenumbers the inertial range covers and this is consistent with our results with run F0.3 presenting the longest inertial range.

There is a limit, however, to how low we can push the Froude numbers in a decaying context as exemplified by run R26F0.1 where no stratified inertial range is visible. We interpret this as the effect of the vertical dissipation of large-scale horizontal motions which becomes important at R = O(10). Clearly, since R = Re Fr 2 h , this occurs at earlier times for runs that start from low Fr ss . Now, the stratified cascade with its inertial range inevitably sets in after the dissipation peak and formation of smallest scales. As a result it appears that, for the case R26F0.1, there is not enough time for the spectrum to become fully developed before vertical diffusion dominates. To add some numerical evidence to this argument, all runs present a peak in dissipation at t = 4.4τ 0 and we have R = 15 at a later time highly dependent on Fr ss . For the case R26F0.1 this time is t = 6.6τ 0 (for the other runs the times are given approximately in figure 9), and hence 2 initial eddy turnover times are not sufficient for the stratified cascade to fully develop. On the other hand run R80F0.2 has about 4.5 turnover times before diffusion starts becoming important so that an inertial range can be established.

The bottom plot in figure 9 corresponds to run R49F0.1 at an early time. At this instant, there is no plateau observable on the horizontal spectrum of horizontal kinetic energy, but rather a dip at k h ≈ 20 and a bump at k h ≈ 80, similarly to what observed by [START_REF] Augier | Stratified turbulence forced with columnar dipoles: numerical study[END_REF].

The compensated spectra in figure 9 are for a given time instant. It is interesting to observe the evolution of the spectra in time to check if the scaling with the kinetic energy dissipation is correct. Results for the one-dimensional horizontal spectra of horizontal KE over a range of times are presented in figure 10 for cases R330F0.6 and R80F0.2. It is clear that the curves collapse well over two to three integral timescales, showing that the dependence on the kinetic energy dissipation is met in our simulations. Moreover, the constant describing the compensated horizontal spectra is C 1 = 0.44 ± 0.02 for the four runs excluding run R26F0.1. For times beyond the ones considered in figure 10 the spectra are viscously dominated without much of an inertial range, highlighting that the buoyancy Reynolds number is low at late times.

It is important to note that k -5/3 h horizontal spectra are presented well before the dissipation scaling ǫ ∼ u 3 h /ℓ h is shown, meaning that the spectra result is applicable to the strongly stratified turbulence regime. For instance, the high-resolution run has ǫ k,hv /ǫ k = 0.419 at the time at which the spectra are shown, which is not too far from the isotropic value of ǫ k,hv /ǫ k = 4/15 = 0.267. The -5/3 form of the horizontal kinetic

k h k h k h k h Ẽh (k h ) Ẽh (k h ) Ẽh (k h ) Ẽh (k h ) 10 0 10 1 10 2 10 3 0 0.2 0.4 0.6 0.8 1 k h Ẽh (k h ) Figure 9. Compensated horizontal spectra of horizontal kinetic energy, Ẽh (k h ) = ǫ -2/3 k k 5/3 h E h (k h ).
Top left: run R18.kF1.8 at t = 20.4τ0 (R = 14.6, Fr h = 0.16), top right: run R330F0.6 at t = 12.3τ0 (R = 15.9, Fr h = 0.12), middle left: run R80F0.2 at t = 8.8τ0 (R = 15.6, Fr h = 0.096), middle right: run R26F0.1. at t = 5.4τ0 (R = 20, Fr h = 0.11). The bottom panel is the run R49F0.1 with its horizontal spectrum of horizontal kinetic energy at t = 5.6τ0 (R = 29.7, Fr h = 0.11).

energy spectra is a statement that the dissipation scaling holds on a scale-by-scale basis, or v 3 h /s h ∼ ǫ k where v h and s h are local scales in the cascade. It is therefore surprising that this scaling holds in the inertial range in our DNS runs but, at the same times, not at the integral scales, which are feeding the downscale cascade of energy in the horizontal. At high enough R and buoyancy Reynolds number it is expected that both results will be

(k h ) = ǫ -2/3 k k 5/3 h E h (k h
) for (a) R330F0.6 from t = 11τ0 up to t = 13.6τ0, (b) R80F0.2 from t = 7.4τ0 up to t = 9.6τ0.

recovered, but in our simulations the horizontal-transfer-to-dissipation balance is found scale-by-scale but not in a volume averaged sense.

Horizontal spectra of potential energy

The horizontal spectra of potential energy have the prediction is

E p (k h ) ∼ ǫ p ǫ -1/3 k k -5/3 h
, which is found from cascade arguments (analogous to those that lead to the form for the horizontal spectrum of horizontal kinetic energy) and on insisting that the timescale of the buoyancy fluctuations is the same as the convective timescale. It follows from this that the form of the horizontal spectrum is a Corrsin-Obukhov spectrum, which governs passive scalar advection. Figure 11 shows five compensated horizontal spectra of potential energy, one for every run at the same times as for the horizontal kinetic energy spectra in figure 9. Except for the high-resolution run R49F0.1 presenting a clear inertial range behaviour, the results are less convincing than the KE spectra. Nonetheless the spectra present a limited inertial range as shown more clearly by the time evolution of the potential energy spectra for runs F0.6 and F0.3 in figure 12. The constant appears to be C 2 = 0.61 ± 0.03 across the runs, except for case R26F0.1 where, as for the kinetic energy spectra, there is no clear plateau.

Horizontal energy flux spectra

We have calculated the horizontal energy flux spectra Π k (k h ) and Π p (k h ) by summing the transfer spectra over cylinders of radius k h in Fourier space; here k h is a two-dimensional wavenumber, k h = k 2

x + k 2 y . For more details on the definition of kinetic energy and potential energy transfer spectra the reader is referred to [START_REF] Augier | Spectral analysis of the transition to turbulence from a dipole in stratified fluid[END_REF]. The horizontal energy flux spectra highlight the energy transfers in the horizontal direction at a certain wavenumber k h . If Π k (k h ) > 0 then there is downscale transfer of kinetic energy to higher wavenumbers than k h and if Π k (k h ) < 0 there is an upscale transfer at this wavenumber. Moreover, a region of constant Π k or Π p defines an inertial range (see [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]. In figure 13 we show Π k and Π p together with the corresponding 2-D horizontal spectra for runs R80F0.2 and R49F0.1. The times at which the energy flux spectra are shown are towards the end of the period over which the 1-D spectra present a collapse over time (see Figure 10). The energy spectra are compensated to show a k -5/3 h -range, which is very limited in run R80F0.2 and slightly more pronounced h Ep(k h ). The runs and corresponding times are the same as those given in figure 9.

in the high-resolution run R49F0.1. In all cases the energy flux spectra are constant over a narrow range, which correlates well with the k -5/3 h -region in the 2-D spectra (except for Π k in run R80F0.2 which is flat starting from wavenumbers k h < 10, where there is clearly no k -5/3 h -region). In both runs, the region Π p = const extends to higher wavenumbers than the region Π k = const. The ranges are narrow due to the limited dynamic range in the runs, however the fact that Π k and Π p even present a region of constancy provides confirmation that we indeed have an inertial range in the horizontal direction. The kinetic and potential energy fluxes are negative at small wavenumbers showing that there is a small upscale transfer at large scales. However, in both the runs under consideration, the transfer of kinetic energy to the shear modes k h = 0 is negligible at the times shown, and it is actually downscale for run R80F0.2 from this mode. At other times, there is a transfer of kinetic energy towards the k h = 0 mode but this mode never dominated the dynamics in our simulations.

As is apparent from figure 13 we have Π k < ǫ k and Π p < ǫ p in the inertial range. If we form total transfers Π = Π k + Π p then the ratio of the total transfer to the total dissipation is Π/ǫ ≈ 0.53 in the inertial range for both runs under consideration. Hence only about half of the dissipation rate can be attributed to downscale transfer of energy in the horizontal in our decaying stratified turbulence. This means that the remaining half of the dissipation rate must be accounted for by energy transfers in the vertical direction or by dissipation acting at scales of the order k h ≈ 30 or lower, presumably due to vertical shear at these horizontal scales. The cascade of energy in the horizontal has therefore been confirmed in our DNS of stratified turbulence but perhaps it should be put into the perspective that it accounts for only roughly half the total energy dissipation in the flow. Further simulations with a greater dynamic range are needed to confirm this statement.

Vertical energy spectra

Vertical spectra of horizontal kinetic and potential energy are predicted to scale as

E h (k v ) ∼ E p (k v ) ∼ N 2 k -3 v
according to the strongly stratified turbulence theory. This particular form for E h (k v ) means that the vertical Froude number is of order unity on a scale-by-scale basis. It is important to test whether the vertical spectra in our decaying runs obey this scaling or not. In figure 14(a), vertical spectra of horizontal kinetic energy for all four cases at specific times are presented, the spectra having been compensated so that a horizontal plateau would show accordance with the scaling above. There is no evidence of such a plateau in any of the DNS runs, the spectra being considerably shallower than k -3 v before the rapid viscous decay at high wavenumbers. Similar vertical spectra were found by Bartello & Tobias in their decaying simulations (Tobias 2014).

It should be pointed out that the range of vertical wavenumbers at which the k -3 vspectrum is predicted is narrow, ranging the scales from the buoyancy scale u h /N to the Ozmidov scale ℓ oz = ǫ k /N 3 , which scales like Fr -1/2 h (see [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]. Given the moderately low Froude number in our runs, which is F r h = O(0.1) during the time when R = O(10), the vertical spectra results are not so surprising. Limited support has The vertical spectra in our DNS runs at intermediate to late times are actually rather similar to their horizontal equivalents, as shown in figure 14(b) for both kinetic and potential energy. The vertical and horizontal spectra present closely matching trends at varying wavenumber meaning that the vertical spectra are closer to k -5/3 v than k -3 v over the stratified inertial range. Beyond a transition wavenumber that corresponds to the Ozmidov scale ℓ oz it is expected that there is a return to an isotropic Kolmogorov spectrum of the form E [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF][START_REF] Augier | Spectral analysis of the transition to turbulence from a dipole in stratified fluid[END_REF]. Our shallow vertical spectra may be the result of an early transition to isotropic conditions due to the moderately low Fr h . We will discuss this point further in the discussion that follows in §5.

h (k v ) ∼ ǫ 2/3 k k -5/3 v (see

Discussion

The fact that the vertical spectra do not present a N 2 k -3 v -dependency has an important implication: the vertical Froude number is not of order unity on a scale-by-scale basis in our simulations. We will now see that this is actually not so surprising if we consider that our scaling argument from which we found Fr v ∼ 1 implies that Fr h ≪ 1. Let us consider scale-by-scale Froude numbers Fr s h = v h /N s h and Fr sv = v h /N s v , where v h , s h = π/k h and s v = π/k v are local scales in the cascade. The reason for using a factor of π and not 2π when relating a scale to its corresponding wavenumber is that s h and s v are scales and not wavelengths, i.e. a region of correlated positive buoyancy such as a yellow structure in figure 3 presents scales s h and s v while if we consider also an adjacent blue structure we will have a wavelength. We found some evidence from our horizontal spectra that at intermediate scales in the inertial range we should have ǫ ∼ v 3 h /s h . This leads to

Fr s h = v h N s h ∼ (ǫs h ) 1/3 N s h = ǫ 1/3 k 2/3 h π 2/3 N (5.1)
so it is clear that Fr s h steadily grows from its integral scale value Fr h as we move to smaller scales. Run F0.3 presents the clearest inertial range of E h (k h ), see figure 9, but on substituting the values of dissipation at this time in (5.1), we see that Fr s h > Fr v for k h > 80. This value compares well with the wavenumber at the end of the stratified inertial range, estimated to be k h = 50, as can be seen in figures 9 and 10. However, this value of k h is only about 10 vertical wavenumbers after the peak in E h (k v ) so it appears that the condition Fr s h ≪ 1 is not well satisfied for most of the vertical scales and Fr sv can therefore not be of order unity for these scales. The typical argument reads that the horizontal Froude number becomes unity at the Ozmidov scale with k oz = π N 3 /ǫ k ≈ 470 for run F0.3 at t = 8.8τ 0 (the time of figure 9). But, by using the value of Fr v from our simulations, the estimate for the transition wavenumber is greatly reduced (the use of the total dissipation also reduces this value, but by a much smaller amount), making it practically impossible to observe a vertical inertial range as predicted by the strongly stratified turbulence theory.

In the future, greater computer capabilities will allow Fr h to be pushed much lower, while keeping R > O(10), thus allowing Fr s h to be small over a large range of scales.

The scaling analysis of §2 should therefore be valid over these scales, for which Fr sv ∼ 1 as a result. It follows that the level of anisotropy decreases as we go down the cascade:

v z v h ∼ s v s h = Fr s h Fr sv ∼ Fr s h ∼ ǫ 1/3 k 2/3 h N (5.2)
where v z is a local vertical velocity scale and we have used continuity in the first step. Therefore the distribution of vertical kinetic energy across horizontal scales is given by

v 2 z ∼ ǫ 1/3 k 2/3 h N 2 ǫ k h 2/3 ∼ ǫ 4/3 k 2/3 h N 2 , (5.3)
which grows as we move to smaller and smaller scales up to a transition wavenumber that is of the order of k oz . This means that, as discussed previously by [START_REF] Riley | Recent progress in stratified turbulence[END_REF]; [START_REF] Davidson | Turbulence in rotating, stratified and electrically conducting fluids[END_REF], the vertical velocity fluctuations peak at the Ozmidov scale, ℓ oz , rather than at the integral scale ℓ h . From (5.3) we therefore have

< u 2 z >∼ v 2 z (k h = k oz ) ∼ ǫ/N from which the following result is obtained, < u 2 z > u 2 h ∼ ǫ/N u 2 h ∼ u h N ℓ h = Fr h (5.4)
rather than < u 2 z >= Fr 2 h u 2 h , the naive deduction obtained by applying continuity without realizing that vertical velocity is a small scale quantity in stratified turbulence.

We have attempted to test equation (5.4) in our simulations. However, large oscillations with frequency equal to twice the Brunt-Väisälä frequency on the < u 2 z >-curve, appear to cause a mismatch between < u 2 z > /u 2 h and F r h . This can be observed in figure 15. By removing the mode k z = 0 from the vertical 1-D spectra of u 2 z , one isolates the dynamics that is independent of this k z = 0 wavenumber, at which the oscillations are expected according to the IGW dispersion relation. The summation of the 1-D spectra of u 2 z from the first non-zero vertical wavenumber to the largest vertical wavenumber is presented together with the evolution of F r h in time in figure 15. The match of functional forms is good for both low-F r h,ss runs considered: run R80F0.2 and the high-resolution run R49F0.1.

Another observation that we may add to support this relation is that visualisations of u z appear to be isotropic as seen in figure 16, at variance with the other fields, u x , u y and b, which invariably present the characteristic layering of stratified turbulence (see figures 3-5). This is suggestive of u z being a smaller scale quantity, as at small scales we have a return to isotropy. It is hoped that further simulations of stratified turbulence at lower Fr h will confirm these preliminary findings.

Conclusions

We have presented results from a scaling analysis and direct numerical simulations of stratified turbulence starting from a high value of R and a Froude number of order one. The alternative scaling analysis based upon the vorticity equations, leads simply and directly to the well-known scaling result for the vertical lengthscale of the turbulence as ℓ v = u h /N . In doing so it highlights the importance of horizontal vorticity in stratified turbulence, and in particular of vertical shearing of horizontal velocity at large scales (which can be interpreted physically as pancake vortices sliding one on top of the other). A frequent assumption in the stratified turbulence literature is ǫ k ∼ u 3 h /ℓ h and ǫ p ∼ u 3 h /ℓ h , which can be useful when interpreting the results of forced simulations. The results of our runs have shown its validity, yet only for the transition regime between the strongly stratified turbulence and the viscous regimes. This result does suggest that the dissipation scaling is robust, holding across the regimes of stratified turbulence. The horizontal spectra are largely in line with the Lindborg forms over a limited inertial range. The presence of an inertial range is confirmed by the energy flux spectra Π k and Π p being approximately constant over a narrow k h -range. The k -5/3 h -spectra are observed from early times so that we are still in a state close to strongly stratified turbulence. The reason why the spectra present k -5/3 h -ranges at times when the dissipation does not show ǫ ∼ u 3 h /ℓ h is not understood, and is a point that deserves further attention. The constants C 1 and C 2 for the horizontal kinetic and potential energy spectra are different in our runs, while previous results had them as being very close with C 1 = C 2 ≈ 0.5 as found by [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]; Lindborg & Brethouwer (2007). More importantly, the vertical spectra E h (k v ) and E p (k v ) are very far from k -3 v and are closer to k -5/3 v over a range of wavenumbers. This result is thought to be due to the not sufficiently low horizontal Froude numbers in the cases under consideration. It is possible that by pushing horizontal Froude numbers down to very low values Fr h = O(10 -3 ) the form N 2 k -3 v for the vertical spectra will be recovered and the scaling Fr v ∼ 1 will be found to hold over a range of

Figure 1 .

 1 Figure 1. Evolution of kmaxη in time for run R80F0.2

Figure 2 .Figure 3 .Figure 4 .Figure 5 .

 2345 Figure 2. Evolution of box-integrated kinetic energy components and potential energy for case R26F0.1. The horizontal velocity uh = < 1/2(u 2 x + u 2 y ) > so that u 2h is the total horizontal kinetic energy and 1/2 < u 2 z > is the vertical kinetic energy. The potential energy in our runs is given by Ep =< b 2 > /2N 2 .

Figure 6 .

 6 Figure 6. Vertical Froude number evolution at times for which R > 7 for DNS runs. Full Fr v evolution given in inset.

Figure 7 .Figure 8 .

 78 Figure 7. (a) Total energy dissipation ǫ = ǫ k + ǫp normalized by u 3 h /ℓ h as a function of time for four runs at different initial Fr . The vertical black segment indicates, for each run, the bound between times at which R > 7 and later times at which R < 7. (b) The viscous dissipation and diffusive dissipation normalized analogously for run R49F0.1 with the same bound indicated on the plot.

Figure 10 .

 10 Figure 10. Horizontal spectra of horizontal kinetic energy compensated such that Ẽh (k h ) = ǫ

Figure 11 .

 11 Figure 11. Compensated horizontal spectra of potential energy, Ẽp(kh) = ǫ -1 p ǫ 1/3 k k 5/3

Figure 12 .h

 12 Figure 12. Compensated horizontal 1-D spectra of potential energy with Ẽp(kh) = ǫ -1 p ǫ 1/3 k k 5/3 h Ep(k h ) for (a) F0.6 from t = 11τ0 up to t = 13.1τ0, (b) F0.3 from t = 7.4τ0 up to t = 8.8τ0.

Figure 13 .h

 13 Figure 13. Energy flux spectrum in the horizontal direction together with corresponding horizontal 2-D spectrum (compensated): (a) horizontal kinetic energy spectrum Ẽh (k h ) = ǫ -2/3 k 5/3 h E h (k h ) and kinetic energy flux spectrum for run R80F0.2, (b) potential energy spectrum Ẽp = ǫ -1 p ǫ 1/3 k k 5/3 h Ep(k h ) and potential energy flux spectrum for run R80F0.2 (both (a) and (b) are at t/τ0 = 9.3), (c) horizontal kinetic energy spectrum and horizontal kinetic energy flux spectrum for run R49F0.1, (d) potential energy spectrum and potential energy flux spectrum for run R49F0.1 (both (c) and (d) are at t/τ0 = 7.5). The horizontal wavenumber is shifted by 0.5 to the right to show the transfers to the shear modes with k h = 0.

Figure 14 .

 14 Figure 14. (a) Compensated vertical spectra of horizontal kinetic energy, Ẽh (kv) = N -2 k 3 v E h (kv). (b) Horizontal and vertical 1-D spectra of horizontal kinetic energy (main plot) and potential energy (inset) for case R1.8kF1.8 at t = 27.2τ0.

Figure 15 .

 15 Figure 15. Evolution of the non-dimensionalized vertical velocity variance and the horizontal Froude number in (a) run R80F0.2, (b) run R49F0.1. The solid grey curve (cyan online) corresponding to (< u 2 z > | kz =0 /u 2 h ) represents the variance of uz (non-dimensionalized by u 2 h ), with the removal of the contributions to the variance of the wavenumber kz = 0. The solid dark curve (red online) corresponds to the horizontal Froude number. The dashed curve (green online) with the large and periodic oscillations is the variance of u 2 z non-dimensionalized by u 2 h considering all vertical wavenumbers.

Figure 16 .

 16 Figure 16. Visualisation of the uz field in run R80F0.2 at t/τ0 = 10, corresponding to Froude number Fr h = 0.09 and R = 13.3.

Table 1 .

 1 RunNx = Ny Nz Ress Fr h,ss Rss Re b,ss R end t end /τ0 DNS runs including horizontal and vertical resolutions and parameters at the peak in kinetic energy dissipation (denoted by the subscript ss, standing for "small-scales"). The runs are denoted by their value of Rss and Fr h,ss . The initial eddy turnover time is τ0 = ℓ0/u0.

	R1.8kF1.8	1536	256 570 1.76 1770 790	1.0 117.3
	R330F0.6	1536	256 860 0.62 330 170	1.5	65.8
	R80F0.2	1536	256 1410 0.24	80	30	1.7	40.3
	R26F0.1	1536	256 1820 0.12	26	4.2	0.7	36.5
	R49F0.1	2048	512 2220 0.15	49	13	6.3	15.7

  giving a timescale for this turbulent flow as t = nℓ h /(A k u h ) ∼ ℓ h /u h (the usual eddy turn-over time). Substituting this into our diffusion length ℓ v we have ℓ v ∼ νℓ h /uh = ℓ h / √ Re, a scaling for ℓ v first proposed by[START_REF] Godoy-Diana | Vertical length scale selection for pancake vortices in strongly stratified viscous fluids[END_REF] to describe the viscous regime of stratified turbulence. The final step is then to

		1.2
		1					R1.8kF1.8
							R330F0.6
		0.8					R80F0.2
							R26F0.1
	ǫℓ h u 3 h	0.6
		0.4
		0.2
		0
		0	20	40	60	80	100	120
					t/τ0
				(a)
		0.35
		0.3
		0.25
		0.2
		0.15
		0.1
		0.05
		0	0	5	10		15
					t/τ0

† Email address for correspondence: maffioli@mech.kth.se

scales. The computational requirements for such a run are very high, but they could lead to the resolution of an essential open issue in stratified turbulence.

The authors would like to thank Prof P.K. Yeung of Georgia Tech for valuable discussions and for providing the pseudo-spectral code used for the simulations.