F F I O L I

P A D A V I D S O N

S B D A L Z I E L 2 A N

D N S W A M I N A T H A N

The evolution of a stratified turbulent cloud

Localized regions of turbulence, or turbulent clouds, in a stratified fluid are the subject of this study, which focuses on the edge dynamics occurring between the turbulence and the surrounding quiescent region. Through laboratory experiments and numerical simulations of stratified turbulent clouds, we confirm that the edge dynamics can be subdivided into materially-driven intrusions and horizontally-travelling internal wave-packets. Threedimensional visualizations show that the internal gravity wave-packets are in fact largescale pancake structures that grow out of the turbulent cloud into the adjacent quiescent region. The wave-packets were tracked in time, and it is found that their speed obeys the group speed relation for linear internal gravity waves. The energetics of the propagating waves, which include waveforms that are inclined with respect to the horizontal, are also considered and it is found that, after a period of two eddy turn-over times, the internal gravity waves carry up to 16% of the cloud kinetic energy into the initially quiescent region. Turbulent events in nature are often in the form of decaying turbulent clouds, and it is therefore suggested that internal gravity waves radiated from an initial cloud could play a significant role in the re-organization of energy and momentum in the atmosphere and oceans.

Introduction

Turbulent flows in the presence of strong stable stratification occur in the atmosphere, oceans and lakes. Turbulence episodes in practice are often localized and without a continuous source of energy as a result of, for example, local shear instability or internal gravity wave breaking (see the review of Riley & Lelong 2000, and references therein). Localized regions of turbulence have been observed in the stratified oceanic thermocline by [START_REF] Grant | Some observations of the occurrence of turbulence in and above the thermocline[END_REF], according to whom "the most striking conclusions from this record are that there is turbulence in the thermocline and that it is in patches". In the review of [START_REF] Gregg | Diapycnal mixing in the thermocline: a review[END_REF] oceanic data is presented showing that mixing in the stratified ocean presents large scale intermittency, which supports the idea of turbulence occurring in localized patches. Such patchiness of the turbulent activity on scales greater than the largest eddies has also been found by [START_REF] Mahrt | Intermittency of atmospheric turbulence[END_REF] when analyzing atmospheric aircraft data. Motivated by these observations, we consider the problem of a localized and decaying region of turbulence sitting in an otherwise unperturbed and stably stratified environment. We call such a localized turbulent region a turbulent cloud ; we consider stratified turbulent clouds, where the cloud interior and the surroundings are initially linearly stratified. Stratified turbulent clouds are governed by inhomogeneous dynamics, in contrast to the dynamics of homogeneous stratified turbulence. The latter class of problems has been covered extensively, with theoretical (see [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]), numerical (see [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF][START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Waite | Stratified turbulence dominated by vortical motion[END_REF] and experimental approaches (see [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]. The inhomogeneous problem, on the other hand, has received relatively little attention, but its manifestation in nature makes it an interesting problem to study.

Turbulent clouds present a turbulent/non-turbulent interface between the cloud and its quiescent surroundings. A number of experimental studies on this interface and its dynamics were published some years ago (see [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF][START_REF] Thorpe | On the layers produced by rapidly oscillating a vertical grid in a uniformly stratified fluid[END_REF][START_REF] Ivey | Boundary mixing in a stratified fluid[END_REF][START_REF] Silva | Experiments on collapsing turbulent regions in stratified fluids[END_REF]. In all these experiments a grid is oscillated, either vertically or horizontally, in a stratified fluid creating a localized turbulent region, which after some time collapses under the effect of gravity with the formation of laterally spreading intrusions. Turbulent mixing in the intrusion interiors results in smaller absolute density gradients compared to the ambient linear stratification; these horizontal density differences lead to horizontal pressure gradients that push the intrusions forwards into the quiescent region (see Figure 1). The displacement of dye lines in the initially quiescent region is a signature of internal wave emission from the turbulent cloud, which was observed by some authors (see [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF][START_REF] Silva | Experiments on collapsing turbulent regions in stratified fluids[END_REF]. [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF] interpret the waves as being horizontally-travelling columnar modes excited by the motion of the intrusions. Columnar modes are known to propagate in front of single intrusive gravity currents generated by the release of mixed fluid of intermediate density in a stratified container (see [START_REF] Manins | Intrusion into a stratified fluid[END_REF]. Density-driven intrusions with their associated staircase-like density profiles have been observed in lakes (see e.g. [START_REF] Caldwell | Thermal microstructure on a lake slope[END_REF][START_REF] Hogg | Mixing of an interflow into the ambient water of Lake Iseo[END_REF] and in the ocean (see [START_REF] Armi | Some evidence for boundary mixing in the deep ocean[END_REF].

In a recent numerical study, [START_REF] Basak | Dynamics of a stratified shear layer with horizontal shear[END_REF] considered the evolution of a horizontal shear layer in a fluid with a vertical stable density gradient. The horizontal shear layer initialized with broadband fluctuations, rapidly undergoes transition to turbulence with the formation of a turbulent cloud surrounded by two opposite-moving laminar streams. The edge dynamics are considered and intrusions spreading into the surrounding streams are observed. The internal wave field is highlighted using the horizontal divergence -∂u z /∂z (u z is the vertical component of velocity), although this tool cannot identify any horizontal waves which are present in the problem (horizontal wave modes, as we will see, do not possess any vertical velocity). When using horizontal divergence as a signature of internal wave emission, inclined wave beams are found to propagate in the quiescent region. Turbulent wakes are another example of localized turbulent region and several laboratory experiments and numerical simulations have treated the dynamics of stratified turbulent wakes including the radiated internal wave field [START_REF] Gilreath | Experiments on the generation of internal waves in a stratified fluid[END_REF][START_REF] Lin | Internal waves generated by the turbulent wake of a sphere[END_REF][START_REF] Bonneton | Internal waves produced by the turbulent wake of a sphere moving horizontally in a stratified fluid[END_REF][START_REF] Spedding | Turbulence, similarity scaling and vortex geometry in the wake of a towed sphere in a stably stratified fluid[END_REF][START_REF] Diamessis | A spectral multidomain penalty method model for the simulation of high reynolds number localized incompressible stratified turbulence[END_REF][START_REF] Abdilghanie | The internal gravity wave field emitted by a stably stratified turbulent wake[END_REF]. Stationary lee waves generated by the sphere or body moving with respect to the stratified fluid are observed, in addition to "random internal waves" [START_REF] Gilreath | Experiments on the generation of internal waves in a stratified fluid[END_REF] that are observed to propagate with an inclination angle with respect to the horizontal plane [START_REF] Bonneton | Internal waves produced by the turbulent wake of a sphere moving horizontally in a stratified fluid[END_REF][START_REF] Diamessis | A spectral multidomain penalty method model for the simulation of high reynolds number localized incompressible stratified turbulence[END_REF][START_REF] Abdilghanie | The internal gravity wave field emitted by a stably stratified turbulent wake[END_REF]. Gravitational collapse of the wake could give rise to intrusions growing in the spanwise direction and forcing horizontal internal gravity waves. However, due to the little mixing occurring within the wake, clear intrusions have not been clearly observed in stratified wake experiments and simulations, which lead some authors to discount this mechanism as a significant source of IGW generation (see [START_REF] Lin | Internal waves generated by the turbulent wake of a sphere[END_REF] and even to put into question the fact that high-Fr stratified turbulent wakes eventually collapse under gravity (see [START_REF] Abdilghanie | The internal gravity wave field emitted by a stably stratified turbulent wake[END_REF]. The absence of layered intrusions in stratified wakes makes the problem quite distinct from the turbulent cloud problem under consideration, where intrusions are known to exist and are a prominent feature of the edge dynamics.

In the context of homogeneous stratified turbulence, it has been know for a long time (see [START_REF] Lin | Wakes in stratified fluids[END_REF][START_REF] Kimura | Diffusion in stably stratified turbulence[END_REF] that initially isotropic turbulence in a stratified fluid rapidly becomes anisotropic, with the formation of pancake-like structures in its interior. The formation of the pancakes is sometimes known as layering and its causes are yet unclear, with a number of mechanism having been put forward over the years, for example the decorrelation instability of [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] and the zig-zag instability (see [START_REF] Billant | Experimental evidence for a new instability of a vertical columnar vortex pair in a stronly stratified fluid[END_REF]. It should be noted that both the zig-zag and decorrelation instability mainly affect the potential vorticity component of the flow (or vortical mode, see [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF]; no layering mechanisms involving the internal wave component have been proposed thus far. In the case of the zig-zag instability acting on a columnar vortex pair, it is not clear how this layering mechanism would manifest itself in turbulence that is initially isotropic.

There are a number of analogies between stratified and rotating turbulent flows: both classes of turbulence display anisotropic structures in their interior (the structures are elongated along the rotation axis in rotating turbulence), and both stratified and rotating fluids sustain wave propagation in their interior. In a rotating fluid the restoring force is provided by the Coriolis force, which is responsible for oscillations in the interior of the fluid known as inertial waves [START_REF] Lighthill | Waves in fluids[END_REF]. Recently, [START_REF] Davidson | On the evolution of eddies in a rapidly rotating system[END_REF] proposed a theory according to which an isolated eddy in a rotating environment deforms into a columnar structure as a result of linear inertial wave propagation. This result could be relevant also for turbulent clouds composed of a multitude of eddies. The authors also present the results from experiments on a turbulent cloud in a rotating system showing that columnar structures emerge from the cloud and grow into the quiescent region (see Figure 2). The structures are indeed found to be created by the linear superposition of inertial waves, as demonstrated by comparing their growth rate to the linear group speed of inertial waves. These results provide motivation for studying similar internal gravity wave mechanisms in the context of stratified turbulent clouds.

Having considered the previous work on turbulent clouds in stratified fluids, we have indentified a number of open questions regarding this flow. Firstly, is there a manifestation of pancake formation or layering in this inhomogeneous problem? If so, what causes the formation of pancakes in the flow? Additionally, what is the complete internal wave dispersion pattern from the turbulent cloud? Finally, are the internal waves linear and how much energy do they carry away from the cloud? We attempt to answer these questions using laboratory experiments and numerical simulations in conjunction.

The remainder of this paper is structured as follows: in §2 results from laboratory experiments on stratified turbulent clouds are presented, and in §3 direct numerical simulations of the problem are described and discussed. A discussion follows in §4 comparing and contrasting the experiments and the DNS, and exploring in more detail the analogy between stratified and rotating turbulent clouds. Finally a brief set of concluding remarks is given in §5. Throughout the paper we sometimes use the following abbreviations: DNS for direct numerical simulation, IGW for internal gravity wave and PV for potential vorticity.

An experiment on turbulent clouds in a stratified fluid

We have performed a variant of the experiments of [START_REF] Thorpe | On the layers produced by rapidly oscillating a vertical grid in a uniformly stratified fluid[END_REF] and [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF]. [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF] created a cloud of turbulence in a stratified fluid by continuously oscillating a grid. The turbulence initially occupies only a small part of the tank, but then spreads in the horizontal direction, with the formation of horizontal intrusions as evidenced by dye initially inserted in the turbulent cloud. We have repeated their experiment except that, instead of continuously oscillating the grid, we oscillate just once in the horizontal direction and restrict this motion to the right-half of the tank, creating an initial cloud of turbulence which is then free to evolve (this is similar to the experiments of [START_REF] Liu | Collapse of a turbulent front in a stratified fluid. 1 Nominally two-dimensional evolution in a narrow tank[END_REF]). The tank we used is 105 cm long, 26 cm wide and 52 cm deep, filled with a linearly stratified salt solution to a depth of 45 cm. The linear stratification was established using the standard double-bucket method [START_REF] Fortuin | Theory and application of two supplementary methods of constructing density gradient columns[END_REF]. To create turbulence in the tank we used a rake of vertical bars so as to avoid any externally-imposed vertical lengthscale (the vertical scale is instead selected by the dynamics); the rake is composed of vertical bars with a bar size b = 6.4 mm and bar spacing M = 31.8 mm. The amplitude of the initial oscillation was ∼ 40 cm, which defines the initial horizontal extent of the turbulent cloud, and three rake speeds were utilized: U ≈ 8, 14, 28 cm s -1 . The rake was set into motion from the right-end of the tank, it was then traversed at a uniform speed up to the centre of the tank, and then traversed uniformly back to its original position.

Decaying turbulence experiments in a non-stratified fluid were performed with grids of the same bar and mesh size driven at similar speeds U by [START_REF] Dalziel | Decay of rotating turbulence: some particle tracking experiments[END_REF] and [START_REF] Staplehurst | Structure formation in rotating turbulence[END_REF] and they found that, at the beginning of the power law region u 2 ∼ t -n , typical velocities and lengthscales of the turbulence are u ≈ 0.05U and l ≈ b. We take these as estimates of u and l at the beginning of the decay also in our stratified experiments, since the initial Froude number is Fr 0 1 in all runs. Experiments were carried out with different initial linear stratifications, with Brunt-Väisälä frequencies N = g/ρ 0 |dρ/dz| (ρ 0 is a reference density taken as the density of the water close to the free surface) in the range N = 0.6-1.4 rad s -1 . The initial Froude number based on turbulence quantities, Fr 0 = u/N l, lies in the range Fr 0 = 0.8-3.6. The turbulence was visualized using two techniques: in some experiments pearlessence flakes, which highlight regions of strong shear (Savaş 1985), were distributed throughout the salt solution (this was achieved by adding the pearlessence to the double-bucket system before filling the tank), while in others fluorescein (a green fluorescent dye) was added to the initial turbulent cloud. The flow was illuminated with a vertical light sheet and images of reflected light intensity captured by a digital camera, viewing in the plane shown in the schematic of Figure 3. Dye marks fluid motion and hence our fluorescein experiments highlight the deformation of the cloud edge due to the turbulent advection of fluid. In contrast, pearlessence is distributed throughout the tank, and therefore shows both the turbulent motion of the fluid initially in the cloud and any internal waves propagating in the quiescent region adjacent to the cloud. Pairs of nominally identical experiments were carried out, using in one case pearl essence flakes and in the other case fluorescein. This results in two sets of images for the same experimental conditions, which can be compared and contrasted. In particular, by "subtracting off" the dyed region shown by the fluorescein images from the pearlessence images, a signature of the wave dispersion from the turbulent cloud is obtained.

Figure 4 shows fluorescein and pearlessence images obtained under the same experimental conditions. In both sets of images horizontal layered structures are seen to emerge from the edge of the turbulent cloud and grow into the quiescent region. From the time of the second frame shown onwards, the horizontal extent of the structures is different across the two visualizations, with the pearlessence structures extending far beyond the fluorescein structures. It would therefore appear that the two types of structures visualized are not the same physical structure. By considering in the pearlessence images only the structures that have extended beyond the dyed region of fluid, we can say that these elongated horizontal structures must be associated with the internal wave field. We are therefore in the presence of layered structures created by the superposition of internal gravity waves. The vertical lengthscale associated with the pearlessence and the fluorescein structures are roughly the same, and hence it is likely that, as suggested by [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF], the intrusions marked by the dye excite horizontally-travelling waves. The presence of horizontal fluid intrusions as well as horizontally-travelling internal waves in the stratified turbulent cloud problem is different to what happens in the rotating cloud problem, where only vertically-travelling waves are seen to emerge from the turbulent cloud [START_REF] Davidson | On the evolution of eddies in a rapidly rotating system[END_REF]). We will return to this point in the discussion section.

In summary, we have performed experiments with stratified turbulent clouds that show the appearance of horizontally-elongated wave-packets propagating ahead of horizontal fluid intrusions. The tools we have used are essentially qualitative tools consisting of different flow visualization techniques. The physical picture obtained from these simple experiments motivated us to study the turbulent cloud problem further using quantitative tools. Our approach has been to use direct numerical simulations to study the wave dispersion from the cloud more in detail, particularly to check the possible linearity of the waves and consider their energetics. The results we obtained are presented below.

Direct numerical simulations

Numerical methods

Motivated by our experimental inspection of stratified turbulent clouds, we performed numerical simulations of a localized region of turbulence in a stratified fluid. The numeri-cal methodology we employed 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 u is the velocity vector, u z the vertical component of velocity, p the perturbation pressure from the hydrostatic pressure distribution, ν the kinematic viscosity, D the density diffusivity and b the buoyancy defined as b = -gρ ′ /ρ 0 (where ρ ′ is the perturbation density from the linear density profile and gravity is acting in the vertical z direction).

The Brunt-Väisälä frequency N is taken as a constant input parameter in the simulations. Equation (3.3) for the buoyancy is obtained either from the energy equation, if the stratification is due to temperature differences, or from the equation for the concentration of a constituent (e.g. salt in water), if the stratification is due to variations of constituent concentration in the fluid [START_REF] Kundu | Fluid mechanics[END_REF]. Equations (3.1)-( 3.3) are solved using a pseudo-spectral method based on Rogallo's algorithm (see [START_REF] Yeung | Lagrangian characteristics of turbulence and scalar transport in direct numerical simulations[END_REF].

The numerical domain is a cube of size l box , and the boundary conditions are periodic in all three directions. 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 decaying stratified simulations were initialized with a fully-developed field of turbulence obtained from a homogeneous isotropic pre-computation in a constant-density fluid. The pre-computation was initialized with random Fourier modes satisfying a prescribed energy spectrum of the form E(k) = Ak 4 exp(-k 2 /a 2 ). During the pre-computation this artificial velocity field is allowed to evolve until the velocity derivative skewness S reaches an approximately steady value at S ≈ -0.5 [START_REF] Davidson | Turbulence. An introduction for scientists and engineers[END_REF]). If k max is the maximum resolved wavenumber, and η and l the Kolmogorov and integral lengthscales at the end of the pre-computation, we only retained the pre-computations with k max η 1.5 and l box /l > 10. This ensures that the small scales are resolved and that the large scales are not too 'boxed-in' (see [START_REF] Pope | Turbulent flows[END_REF]. The fully-developed turbulence thus obtained is then confined to a central slice of the cube by applying a high-order exponential filter, which suppresses velocity fluctuations outside the central region. The width of the filter is 0.2l box . We use this localized turbulent field, an example of which is given in Figure 5, to initialize the stratified simulations. The buoyancy field is initialized with b = 0 everywhere in the domain so that at the beginning of the stratified simulations the density profile is linear throughout the periodic cube, i.e. no mixing has initially taken place within the cloud.

In a single run, Lagrangian particles were added to the simulation. Lagrangian particles consist in idealized fluid elements whose position can be found by integrating dr/dt = u(t, r), where u is the fluid velocity at r. A second-order predictor-corrector scheme is used to find the position of the particles. Since the Eulerian velocity field is only available at the discrete grid nodes, interpolation has to be used to find the fluid velocity at the particle positions. This interpolation is the dominant source of error in the extraction of Lagrangian data [START_REF] Yeung | Lagrangian statistics from direct numerical simulations of isotropic turbulence[END_REF]; in the numerical code fourth-order cubic splines are used to accurately interpolate Eulerian velocity field data. Table 1. Simulation parameters for different runs considered. In all cases the subscript 0 on a certain quantity corresponds to the initial value of the quantity at t = 0. The quantity T0 corresponds to the initial integral timescale T0 = l0/u0. The Reynolds number is Re0 = u0l0/ν. The buoyancy Reynolds number is given at the start of each run as R0 = Re0Fr 2 0 since the turbulence is initially isotropic (we do not give this quantity for the case F1.5 since the Froude number is initially high, far away from being a strongly stratified case). The Kolmogorov scale at the start of the stratified runs is η0 = (ν 3 /ǫ0) 1/4 .

Eulerian visualizations

We now present numerical simulation results, starting from Eulerian visualizations of the flow field. The simulations we performed are listed in table 1, together with the relevant initial parameters and durations of the decaying turbulence simulations. The simulations are named based on their initial Froude number, for example case F0.1 corresponds to a simulation with Fr 0 = 0.1. In listing the initial quantities of the different simulations we do not distinguish between horizontal and vertical directions, as one normally would do in stratified turbulence simulations. The reason behind this is that the initial conditions for all stratified runs are isotropic within the cloud, coming from a pre-computation in a constant-density fluid. For all cases the Schmidt number is taken as Sc = 1.

Figure 6 shows iso-surfaces of u x , the velocity component in the x-direction, for case F0.1 at different times throughout the simulation. The structure of the u x field rapidly becomes pancake-like, with the growth of horizontal flattened structures out of the turbulent cloud and into the quiescent region. In some cases the pancakes detach from the cloud forming isolated structures. The pancakes maintain their coherence for a time in- dye used in our experiments. The other iso-surfaces given in Figure 7 are for u x , u y and b; these quantities are not materially conserved and their fields will be re-arranged both by turbulent advection and by any internal waves present in the simulation. The Π field remains compact at the late time and its extent does not differ greatly from the initial extent of the turbulent cloud. On the other hand, the other three fields present large-scale pancakes extending up to the edge of the domain. Using similar reasoning as in §2, it can be concluded that the large-scale pancakes highlighted by the plots in Figures 6 and7 are internal waves propagating in the quiescent region. These pancake structures are nothing but the three-dimensional equivalent of the horizontal structures evidenced in the experiments with pearlessence (see Figure 4). As pointed out by [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF], there could be motions in a stratified flow that do not possess potential vorticity but which are not internal waves. It is possible that such advective motions could be causing the large-scale structures observed in our runs, which would invalidate our interpretation of them being the result of the superposition of internal waves. To dispel any doubts regarding the nature of the pancakes, we seeded the initial turbulent cloud with Lagrangian particles. This was done for a single case, run F0.1, and the results are presented in the next section.

Lagrangian visualizations

Inertia-less particles were randomly distributed throughout the initial turbulent cloud for case F0.1. The displacement of the particles at a certain time can be compared with the extent of the pancake structures from Figures 6 and7. If the particles are displaced up to the edge of the pancakes then it is likely that these structures are caused by advection of fluid from the cloud. By contrast, if the pancakes lie ahead of the particles, internal waves are probably at the origin of these large scale motions (as the motions are beyond the region of action of turbulent advection).

The particle positions at two times are shown in Figure 8. The two lowest iso-surfaces of u x showing the pancake structures have been superimposed on the particle positions for visual aid. It is clear from the figure that at the early time of N t = 9.6 the particles are confined to a region close to the initial turbulent cloud, whereas the u x field propagates ahead in the quiescent region. This is even more evident at the later time N t = 19.6. These results confirm that the pancake motions in this inhomogeneous problem are created by the superposition of internal waves.

To justify our use of potential vorticity as a tracer in §3.2 we plot the particle positions superimposed on two iso-surfaces of Π at a late time in Figure 9. The Eulerian and Lagrangian fields match closely, but there are some discrepancies at the edge of the cloud which could be due to the choice of the iso-surface levels. This result suggests that fluid advection that does not contain potential vorticity is negligible in our simulations.

Are the pancakes composed of linear waves?

Having confirmed that large-scale pancake structures are created in the turbulent cloud problem by the superposition of internal waves, a natural question that arises is: are these waves linear or non-linear? Because the structures created by the waves can be thought of as wave-packets travelling at the group speed, another way of posing this question would be: do the pancake-shaped wave-packets travel at the group speed as predicted by linear theory?

To answer this question our approach is to track individual pancakes from u x isosurface plots and compare their speed with the expected linear group speed. As can be appreciated by inspection of Figure 6, the pancake structures elongate primarily in the x-direction because the forcing of the intrusions is predominantly in the direction of the inhomogeneity, i.e. the x-direction. Therefore the tracking works by focusing in on an individual pancake-shaped wave-packet -corresponding to an iso-surface with a specific isovalue -and finding the maximum extension along the x-axis of its isosurface. The maximum x-extension of the wave-packet corresponds to the wave-packet front. This procedure is repeated for successive timesteps, and hence plots of wave-packet front displacement against time can be computed. Such plots are presented in Figure 10 for a collection of pancake-shaped wave-packets taken from the low-Froude number simulations, runs F0.03, F0.1 and F0.2. The wave-packet displacement grows linearly with time, meaning that the group speed (given by the slope of the curves) is constant, a result which is consistent with linear wave propagation.

The expression for the group velocity of internal gravity waves can be written as (see [START_REF] Kundu | Fluid mechanics[END_REF])

c g = ± N k z k 3 (k z e h -k h e z ), (3.4)
where k z and k h are the vertical and horizontal wavenumbers so that the wavevector k is given by k = k h e h + k z e z , and k = |k|. We are in the presence of horizontallytravelling wave-packets meaning that the wavevector must be purely vertical with k h = 0. An expression for the group velocity according to linear theory is therefore c g = ±(N/k z )e h = ±(N λ z /2π)e h , where λ z is the vertical wavelength. We therefore would expect the horizontal internal waves created by the initially isotropic turbulent cloud to have a group speed c g ∼ N l 0 , if linear theory applies. To check if the wave-packets are indeed composed of linear internal gravity waves, Figure 10 time non-dimensionalized by multiplication by N . If the horizontally-travelling waves are linear the curves should collapse on a single straight line and indeed they do. So far we have verified the proportionality of the group speed c g to N and l 0 , but we have not confirmed the value of the constant that ensures the exact expression c g = N/k z is enforced. To check that c g = N/k z in the simulations, we considered vertical slices at an x location in the initially quiescent region, with snapshots taken after the passage of the wave-packets. There is a clear vertical wavelength associated with the u x field in such vertical slices (see Figure 11) and we extracted the dominant vertical wavenumber k z by taking the two-dimensional FFT of these fields. It was hence possible to calculate c g = N/k z and compare it with the group speed obtained from the wave-packet tracking. The results are presented in table 2; more than one wave-packet was tracked for every run so an average wave-packet group speed was calculated for every simulation. For all three low-Fr runs the wave-packet c g compares favourably to the expected group speed from linear theory.

In summary, then, we have tested the linearity of the internal waves composing the pancake structures, and the results shown that the waves do indeed act as linear internal gravity waves. However, as can be appreciated from Figure 6, the horizontal internal waves are of finite-amplitude and we therefore refer to them as quasi-linear waves.

Wave dispersion pattern

It is of interest to study the internal wave dispersion pattern from the turbulent cloud. We have already confirmed the presence of horizontally-travelling waves (with frequency ω ≈ 0) in the simulation, but we expect a spectrum of frequencies in the radiated wave field since the initial turbulent cloud contains wavevectors k in all directions. Internal gravity waves with frequencies in the range 0 ω N are expected according to the linear dispersion relation ω = N k h /k. To highlight off-axis waves we use ∂u z /∂z, the negative horizontal divergence, which has been used previously by a number of authors Run N kz at |û 2D x |max cg = N/kz Wave-packet cg F0.03 44.5 3.5 12.7 12.9 F0.1 44.5 5 8.9 9.1 F0.2 22.2 4.5 4.9 4.8

Table 2. Comparison of group speed as expected from linear theory with measured wave-packet speed. The wave-packet cg was averaged over all wave-packets tracked from an individual run.

Only low-Fr runs are considered.

as a signature of the wave field emitted by a turbulent cloud [START_REF] Diamessis | A spectral multidomain penalty method model for the simulation of high reynolds number localized incompressible stratified turbulence[END_REF][START_REF] Basak | Dynamics of a stratified shear layer with horizontal shear[END_REF][START_REF] Spedding | Turbulence, similarity scaling and vortex geometry in the wake of a towed sphere in a stably stratified fluid[END_REF]. It should be emphasized, however, that the ∂u z /∂z field cannot show horizontal internal waves with k = k z e z that contain purely horizontal velocity (by continuity k • u = 0 for linear internal gravity waves), and hence the u x field is used to show these zero-frequency waves.

Figure 12 contains x-z slices of the u x and ∂u z /∂z fields at two times for run F0.1. As previously seen in the three-dimensional visualizations of §3.2, the u x field shows the presence of horizontal internal waves in the initially quiescent region. The horizontal waves are clearly visible from the time N t = 10.9 and create a close-to-uniform layering, similarly to the experimental slices of Figure 4. At this early time, diagonal waves highlighted by ∂u z /∂z have not propagated far from the initial turbulent region, which could be because, from equation 3.4, their horizontal group speed c gh = N k 2 z /k 3 = (N/k z ) sin 3 Θ is lower than the (horizontal) group speed of purely horizontally-propagating waves c g = N/k z . By the later time N t = 32.2 both waveforms (horizontal and diagonal) have filled the computational domain. The typical angle of inclination of the diagonal waves is Θ = 50-55 • from the vertical, corresponding to a wave frequency ω = N cos Θ = 0.57-0.64N . The u x field at the sides of the turbulent cloud is now a complicated superposition of horizontal and diagonal wave beams, since also the waves with ω > 0 have an associated u x field. It is evident from figure 12(c) that also the cloud interior has become layered in u x , though through very different (non-linear) dynamics compared to the largely linear dynamics on the sides of the turbulent cloud.

We have found evidence for the presence of two frequency bands for the internal gravity wave radiation from the turbulent cloud: horizontal internal gravity waves with ω ≈ 0 and inclined waves of finite frequency ω ≈ 0.6N . The second type of IGW is consistent with waves emitted from a well-mixed turbulent patch [START_REF] Dohan | Internal waves generated from a turbulent mixed region[END_REF][START_REF] Taylor | Internal gravity waves generated by a turbulent bottom ekman layer[END_REF] and with IGWs emitted from a stratified turbulent wake (see [START_REF] Bonneton | Internal waves produced by the turbulent wake of a sphere moving horizontally in a stratified fluid[END_REF], but towards the lower end of the spectrum observed in a recent numerical study of wave radiation from a stratified turbulent wake at high Re (see [START_REF] Abdilghanie | The internal gravity wave field emitted by a stably stratified turbulent wake[END_REF]. The horizontal waves are consistent with columnar modes emitted in front of intrusions (see [START_REF] Manins | Intrusion into a stratified fluid[END_REF][START_REF] Maurer | Gravity-driven intrusions in stratified fluids[END_REF]. Apart from the continuously displaced dye lines, which are a signature of columnar wave activity, observed by [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF] and [START_REF] Silva | Experiments on collapsing turbulent regions in stratified fluids[END_REF], these horizontal waves have not previously been observed or characterized in stratified turbulent cloud experiments or simulations.

In part, this is because many previous studies [START_REF] Linden | The deepening of a mixed layer in a stratified fluid[END_REF][START_REF] Dohan | Internal waves generated from a turbulent mixed region[END_REF][START_REF] Taylor | Internal gravity waves generated by a turbulent bottom ekman layer[END_REF] considered turbulent patches that extended to cover the entire horizontal extent of the experimental/numerical domain, which directly rules out horizontal wave propagation. Horizontally localized turbulent clouds have been considered before by [START_REF] Basak | Dynamics of a stratified shear layer with horizontal shear[END_REF], but their use of ∂u z /∂z as a signature of internal gravity waves means that they only observed inclined waves in their simulations. In figure 18 of their paper, vertical slices of ρ ′ , ω y and ω x show the growth of layered structures out of the turbulent cloud; the authors interpret these structures as being solely due to intrusions, but they could represent both intrusions and horizontal waves travelling in front of them. It is interesting to note that in the same figure the authors show ω z (a first order approximation to the potential vorticity) and this quantity remains confined to the turbulent cloud, hinting to the fact that the horizontal structures to the sides of the cloud could in fact be internal gravity waves. The absence of a clear collapse event and of ensuing intrusions in stratified turbulent wakes [START_REF] Lin | Internal waves generated by the turbulent wake of a sphere[END_REF][START_REF] Abdilghanie | The internal gravity wave field emitted by a stably stratified turbulent wake[END_REF]) could be the reason why horizontal IGWs have not been observed before in this context.

Wave energetics

The simulations were performed in the absence of forcing and the energy in the domain should therefore decay in time. A look at the kinetic energy in the box as a function of time shows that E K does not decrease monotonically with time but rather presents largeamplitude oscillations (see Figure 13(a)). These oscillations are in anti-phase with the oscillations on the potential energy plot (the potential energy per unit mass at a point is

E P = b 2 /2N 2 ),
which suggests a process of exchange of kinetic and potential energy and hence the presence of gravity waves in the system. This is confirmed by the evolution of the total energy in the box, E tot = E K + E P , which presents a classical viscous (and diffusive) decay. The period of the large-amplitude oscillations in Figure 13(a) is T ≈ 0.07 giving a frequency ω ≈ 89.8 ≈ 2N , which is higher than the upper bound on the frequency of internal gravity waves in the system ω = N . This apparent contradiction is due to the fact that if we are in the presence of IGWs, u x , u y , u z ∝ exp(iωt) and the kinetic energy will be E K = 1/2(u 2

x + u 2 y + u 2 z ) ∝ exp(2iωt), and so we can interpret the oscillations as being IGWs with ω = N . It therefore appears that high-frequency waves are persistent in our runs and influence the total kinetic energy evolution. They are so visible on the plot of E K because these waves are all created at t = 0 (and hence are in phase and will interfere constructively) as the density field ρ ′ (t = 0) = 0 adjusts to the turbulent velocity field imposed on it. After reaching their maximum amplitude these high-frequency waves gradually dissipate through the action of viscosity. Internal gravity waves with ω = N have zero group speed and correspond to trapped modes that remain confined to the initial turbulent cloud; they should be distinguished from propagating internal waves with non-zero group speed that were identified in §3.5.

To find the amount of kinetic energy in the propagating waves we integrated the kinetic energy in the regions ahead of the turbulent cloud as identified by the potential vorticity field. A threshold was set on |Π| (typically taken as 0.1% of |Π| max ) below which only propagating waves were assumed to exist. The result of this integration is given in Figure 13(b) for run F0.1. It is clear that the kinetic energy in the propagating waves is a small but appreciable fraction of the kinetic energy left behind in the cloud E cloud = E K -E wave , the maximum value of E wave /E cloud being ≈ 0.12 = 12% for run F0.1. The evolution of E cloud , similarly to the total kinetic energy in the box, presents oscillations created by internal waves with ω = N , confirming that there are trapped modes that remain localized within the cloud. The radiated kinetic energy as a function of time does not present any visible oscillations and, after an initial transient, reaches an approximately steady value at t = t peak . This suggests that the horizontal internal waves with ω ≈ 0 contain most of the radiated energy, while the off-axis waves do not contribute appreciably to the energetics. Such a view is consistent with the visualizations in Figure 6 and with the experimental observations of Figure 4, as in both cases only horizontal structures are observed.

In figure 13, the cumulative kinetic energy dissipation t 0 ǫ dt and cumulative potential energy dissipation t 0 χ dt are given, where ǫ = 2νS ij S ij is the kinetic energy dissipation rate and χ = (D/N 2 )∇b • ∇b is the rate of diapycnal mixing, and both quantities are volume averaged. Both ǫ and χ are confined to the turbulent cloud, as is confirmed by plots of these quantities averaged over vertical y-z planes as a function of x (not shown). By equipartition the kinetic energy of the propagating waves should equal their potential energy, which allows a direct comparison of E wave to be made not only with t 0 ǫ dt but also with t 0 χ dt. Hence, if we compare the kinetic energy in the propagating waves with the cumulative kinetic and potential energy dissipation we see that at N t = N t peak = 5.3, E wave is about 50% of the cumulative dissipation and 118% of the cumulative mixing occurring within the cloud, and at the end of the run these values have dropped to 13% and 27% of the integrated dissipation and mixing respectively. These figures show that internal gravity wave radiation is an efficient source of energy dispersion from a cloud of turbulence, especially at the beginning of its evolution.

The percentage of kinetic energy held in the propagating waves relative to the kinetic energy in the cloud is given in table 3 at t = t peak . Also the maximum value of E wave /E cloud is given for each decaying run. The fraction of radiated kinetic energy appears to be constant throughout the runs at t = t peak but the maximum value of this quantity is an inverse function of the Froude number, with the lowest-Fr case presenting the greatest wave energy percentage, equal to ∼ 16%. This is a significant portion of the total E K and means that propagating internal waves generated by a localized turbulence event could be important in the context of meteorology and oceanography as they could alter momentum and energy budgets in different regions of the atmosphere or oceans.

Discussion

Stratified vs rotating turbulence

It is of interest to revisit the analogy between stratified and rotating turbulent flows in light of the results presented in the previous sections. [START_REF] Davidson | On the evolution of eddies in a rapidly rotating system[END_REF] considered the linear initial-value problem of an eddy sitting in a rotating system and found that the component of angular momentum in the rotating reference frame along the direction of the rotation axis remains confined to the cylinder parallel to the rotation axis and circumscribing the eddy. The authors use this result to show that isolated eddies in a rotating fluid will deform as a result of linear internal wave propagation that preferentially channels energy along the rotation axis giving rise to columnar structures. Experimental results are also presented showing that columnar wave-packets obeying linear theory emerge from a rotating turbulent cloud, thus corroborating the theory. In essence, in the rotating turbulent cloud problem, free and spontaneous emission of inertial waves from the cloud results in the formation of columnar structures. At first sight the analogy between rotating and stratified flows appears to carry through to the stratified turbulent cloud problem, where we have shown that quasi-linear internal waves form horizontal pancake-like structures. However, as evidenced by the experimental fluorescein visualizations, non-linear fluid advection also results in the emergence of horizontal layered structures or intrusions from the turbulent cloud. The intrusion height roughly scales with the vertical wavelength obtained from the pearlessence visualizations. It is therefore likely that, as put forward by [START_REF] Browand | The behavior of a turbulent front in a stratified fluid: experiments with an oscillating grid[END_REF], the intrusions push on the fluid immediately ahead of them and thereby create a steadily moving front that propagates into the quiescent region. This means that in the stratified turbulent cloud problem, forced emission of internal gravity waves (with the forcing provided by the intrusions' motion) results in the creation of pancake structures. The analogy, which in any case was not a formal one, therefore seems to break down: the rotating case is a free emission problem, while the stratified case is a forced problem.

Intrusions in experiments and simulations

Moving now to a comparison between the experimental and numerical results, when inspecting the potential vorticity iso-surfaces in Figure 7 one can notice that the structures are not very pronounced and the Π field remains localized even at late times. Recall that we use potential vorticity as a tracer in the same way as we used fluorescein as a tracer in the experiments; therefore the numerical intrusions shown by Π are far less pronounced than the experimental intrusions, to the point that one has doubts whether materially- driven intrusions are at all present in the numerical results. To clarify this it is important to note that at high Reynolds numbers the relevant timescale for the potential vorticity is the advection timescale, or eddy turnover time, T = l/u (see [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF]. Now because the initial Froude number is low in all the runs considered thus far, the typical initial advection velocity u 0 will be much lower than the group speed of horizontal waves c g ∼ N l 0 , meaning that internal gravity waves will evolve much faster than the potential vorticity field. This explains why in Figure 7 simulation with Fr 0 = 0.2 again up to 18 eddy turnover times to check for the presence of intrusions in the simulation. In Figure 14 vertical slices of the potential vorticity are presented at 2 y-locations, for run F0.2ext and for the high-Froude number run F1.5, after the same number of eddy-turnover times. Isolated intrusions have emerged from the turbulent cloud in both cases; the intrusions are less ordered than in the experiments, possibly due to the finite width of the turbulent slice, which is significantly deformed at this late time. In both runs with different Fr 0 the intrusions are about to reach the domain boundaries at the time t ≈ 16.5T 0 , which confirms that the eddy turnover time is indeed the correct timescale to describe the motion of the intrusions. The height of the intrusions seems to differ slightly across the two runs, with the F1.5 intrusions, for which the initial Froude number is higher, having a greater height than the F0.2 intrusions. A similar variation in the intrusion height was observed in the experiments, with growing intrusion heights with increasing Froude number. To quantify this we measured a vertically-averaged intrusion height from the fluorescein images (it was not attempted to do so in the simulations because of the lack of order of the numerical intrusions). A plot of the intrusion height h normalized by the initial turbulence lengthscale l 0 (estimated as the bar size b) as a function of the initial turbulent Froude number Fr 0 = u 0 /N l 0 is presented in Figure 15. The line in figure 15 is a power law with the following functional form, h/l 0 = 6.1(Fr 0 ) 0.313 , which provides a plausible fit to the data even though the experimental error bars are large. We now try to rationalize this behaviour of the intrusion height. The main physical ingredient of the analysis is that intrusions form as a result of the collapse of mixed pockets of fluid within the turbulent cloud. In this discussion we refer to collapse as the re-stratification of the density field that occurs in mixed fluid regions with ρ tot ≈ const, with the associated horizontal spreading of fluid in the form of intrusions. We identify the onset of collapse as the moment when inertial forces that mix the density field approximately equal the buoyancy forces, i.e. when Fr ∼ 1. Now, because Fr 0 1 in all our experiments, the turbulence has time to mix the density field and consequently pockets of mixed fluid will grow. This process lasts until the Froude number has decayed to a value of order unity, when collapse occurs. Hence the height of the intrusions will depend on Fr 0 , with greater Froude numbers allowing more time for the turbulence to mix the density field until it feels the effect of the stratification. To make any progress we need to estimate the initial evolution of the turbulent velocity and lengthscale to be able to estimate the time at which the turbulent Froude number becomes of order of unity. Because Fr 0 1 in all the runs there will be a period of quasithree-dimensional evolution and we might therefore use the estimates for energy decay in homogeneous isotropic turbulence, neglecting as a first approximation the inhomogeneity of the turbulent cloud. In a recent review [START_REF] Davidson | The minimum energy decay in quasi-isotropic grid turbulence[END_REF] starts from the well-established expression du 2 /dt = -Au 3 /l and finds the evolution of the kinetic energy and the integral lengthscale in decaying homogeneous isotropic turbulence. The evolution of u 2 and l can be written as: where u 0 and l 0 are the velocity scale and integral lengthscale of the turbulence at an initial time, taken as t = 0. Combining these two expressions a well-known result is obtained, notably that the integral timescale T of the turbulence grows linearly in time:

T T 0 = l/l 0 u/u 0 = 1 + A n u 0 t l 0 (4.4)
where T 0 is the initial integral timescale. Crucially this expression is also equal to the ratio between the initial Froude number Fr 0 and the Froude number at time t, that is T /T 0 = Fr 0 /Fr = [1 + (A/n)(u 0 t/l 0 )]. When F r = 1, equation 4.4 gives Fr 0 = [1 + (A/n)(u 0 t/l 0 )], which is the quantity of which the power is taken in equation 4.2. Therefore it is also true that at collapse l l 0 = (Fr 0 ) 1-n/2 (4.5)

Up to this moment three-dimensional mixing takes place within the cloud, thus the vertical scale of the mixed fluid regions will be set by l. We have said that intrusions originate from such mixed regions, and hence it follows that the ensuing intrusion height is h ∼ l. A corollary of this argument is that h ∼ l b = u/N (where l b is the buoyancy scale) or, equivalently, h ∼ l oz (where l oz = (ǫ/N 3 ) 1/2 is the Ozmidov lengthscale). This follows from the fact that l oz = (ǫ/N 3 ) 1/2 ∼ (u 3 /lN 3 ) 1/2 = Fr 1/2 l b ∼ l b , when Fr ∼ 1, where we have used the estimate ǫ ∼ u 3 /l and the fact that the turbulence is initially isotropic. Therefore equation 4.5 becomes:

h l 0 ∼ (Fr 0 ) m (4.6)
which is an expression for the intrusion height (the substitution m = 1 -n/2 has been performed). [START_REF] Davidson | The minimum energy decay in quasi-isotropic grid turbulence[END_REF] reports a range of decay exponents, n = 1.10-1.37. The corresponding range for the exponent m in equation 4.6 is m = 0.315-0.45. The exponent obtained from the experimental data is m = 0.313, which is at the boundary of the range predicted by this semi-empirical model. The positive match with the experimental trend suggests that our physical picture is plausible, and that turbulent mixing inside the cloud causes the formation of the laterally-spreading intrusions. In a different context, [START_REF] Spedding | Vertical structure in stratified wakes with high initial froude number[END_REF] predicts the vertical scale of a stratified turbulent wake to be L v ∼ Fr 1/3 0 , i.e. a power law with a very similar exponent to the one we obtained. The main assumption is that the vertical scale is set by the Ozmidov lengthscale at collapse, which as seen above is equivalent to assuming Fr ∼ 1 as we have done.

Conclusions

The evolution of an inhomogeneous cloud of stratified turbulence, and the surrounding flow, is dominated by the growth of pancake structures, and these are composed of intrusions and horizontally-propagating internal gravity waves that are forced by the intrusions' motion. series of events leading to layering in this problem is thought to be as follows: (a) mixing within the turbulent cloud creates regions of ρ ≈ const, (b) when these mixed regions reach a critical vertical scale l ∼ u/N , they collapse under the action of gravity forming intrusions spreading to the sides of the cloud, (c) horizontal waves are forced by the intrusions and propagate at their level, enlarging the layer (as in Figure 4). Such a layering mechanism is quite distinct from the previously proposed layering mechanisms of stratified turbulence, notably the decorrelation mechanism (see [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF]) and the zig-zag instability (see [START_REF] Billant | Experimental evidence for a new instability of a vertical columnar vortex pair in a stronly stratified fluid[END_REF]. These mechanisms mainly affect the PV component of the flow, whereas the mechanism here proposed involves both potential vorticity motions and internal gravity waves. Because mixing within homogeneous stratified turbulence is not uniform but rather occurs in intermittent patches, this layering mechanism need not be restricted to the inhomogeneous problem under study and could be at work in a larger class of stratified turbulent flows.

Figure 1 .

 1 Figure 1. Sketch of motion of intrusions emerging from a stratified turbulent cloud. Return currents are established in between intrusions, with the resulting density profile across the intrusions given on the right of the image (reproduced from Browand et al. (1987) with permission).

Figure 2 .

 2 Figure 2. Image of the flow in rotating turbulent cloud experiment by Davidson et al. (2006). The turbulent cloud was initially confined to the top of the image; image taken at a time Ωt = 10 after initiation of the turbulence (Ω is the system rotation rate about the vertical axis).

Figure 3 .

 3 Figure 3. Experimental set-up

Figure 4 .

 4 Figure 4. Flow visualization of edge of stratified turbulent cloud at three different times. Top: fluorescein visualization, bottom: pearlessence visualization. Experimental parameters: U = 14.5 cm s -1 , N = 1.4 rad s -1 . The initial integral timescale T0 is estimated as T0 = l/u = 0.88 s with u = 0.73 cm s -1 and l = 0.64 cm.

Figure 5 .

 5 Figure 5. Initial condition for stratified DNS: iso-surfaces of ux from fully developed turbulence confined to a vertical slice.

Figure 6 .Figure 7 .

 67 Figure 6. Iso-surfaces of constant ux with every colour corresponding to a different isovalue. Case F0.1 at four different times: (a) N t = 4.7 (b) N t = 8.3 (c) N t = 13.5 (d) N t = 19.8.

Figure 8 .

 8 Figure 8. Particle positions superimposed on two lowest iso-surfaces of ux for case F0.1 at two times: (a) N t = 9.6 (b) N t = 19.8. Every black dot corresponds to a Lagrangian particle.

Figure 9 .

 9 Figure 9. Particle positions superimposed on two lowest iso-surfaces of Π from Figure 7. Case F0.1 at time N t = 19.8.

Figure 10 .

 10 Figure 10. Displacement of the front of 13 different wave-packets tracked in simulations F0.03, F0.1, F0.2 as a function of time. Different symbols correspond to different individual wave-packets. In (b) ∆x and t are normalized by the initial integral lengthscale l0 and N -1 respectively.

Figure 11 .

 11 Figure 11. Vertical slice in the quiescent region after wave passage. In (a) the black line at 0.8L shows the location of the slice, while the grey area corresponds to the initial extent of the turbulent cloud. The snapshot in (b) was taken at a time N t = 13.5 and the shading shows the ux field on the slice.

Figure 12 .

 12 Figure 12. Vertical slices of ux field in (a)-(c) and duz/dz field in (b)-(d) at two different times (the slice is at y = 0.1l box ). Case F0.1 at times: (a)-(b) N t = 10.9 (c)-(d) N t = 32.2.

Figure 13 .

 13 Figure 13. (a) Plots of the evolution with time of EK , EP and Etot for simulation F0.1 (all energies are found by averaging their value throughout the entire box) (b) Evolution of energetics related quantities for case F0.1: kinetic energy in turbulent cloud E cloud = EK -Ewave, kinetic energy in propagating waves Ewave (found by integrating ahead of PV structures), and cumulative kinetic and potential energy dissipation.

Figure 14 .

 14 Figure 14. Vertical slices of Π field for (a)-(b) run F0.2ext at t/T0 = 16.5 (c)-(d) run F1.5 at t/T0 = 16.6. Slices are x-z planes taken at (a) y = 0.1l box (b) y = 0.9l box for case F0.2ext and (c) y = 0.1l box , (d) y = 0.4l box for case F1.5.

Figure 15 .

 15 Figure 15. Intrusion height h normalized by initial integral lengthscale l0 as a function of initial Froude number F r0 = u0/N l0.

Table 3 .

 3 Simulation Fr 0 N t peak Ewave/E cloud at t = t peak (Ewave/E cloud )max Energetics-related quantities for the propagating waves. The fractions Ewave/E cloud are given in percentages.

	F0.2	0.2	4.8	10.6%	10.9%
	F0.1	0.1	5.3	10.2%	11.6%
	F0.03	0.032	6.5	10.1%	16.2%

† Email address for correspondence: aplm2@cam.ac.uk

The authors would like to thank Prof P.K. Yeung of Georgia Tech who provided the pseudo-spectral code used for the computations. A.M. gratefully acknowledges financial support received from EPSRC.


A Maffioli

G Brethouwer

E Lindborg

Mixing efficiency in stratified turbulence

 can be re-interpreted in this light because the observed variation of Γ in their study as well as in datasets from recent oceanic and atmospheric measurements occurs at a Froude number of order unity, close to the transition value F r = 0.3 found in our simulations.

Introduction

The manifestation of mixing in the atmosphere and oceans takes different forms from phyto-plancton blooms in the oceans to the creation of well-mixed fronts in the atmosphere. Mixing in the atmosphere and oceans is a key factor to consider when estimating global energetics. Central to the attempt of quantifying mixing in the oceans is the concept of eddy diffusivity, defined as K ρ = B/N 2 [START_REF] Osborn | Oceanic fine structure[END_REF] where B = -bu z is the buoyancy flux, b = -ρ ′ g/ρ 0 being the buoyancy acceleration related to the density perturbation ρ ′ from the background stratification and u z being the vertical velocity perturbation. The idea is that the buoyancy flux can be modelled as B = -bu z = K ρ d b/dz, where the mean buoyancy gradient is d b/dz = N 2 . [START_REF] Osborn | Estimates of the local rates of vertical diffusion from dissipation measurements[END_REF] inspected the turbulent kinetic energy equation in the presence of buoyancy and assuming steady-state conditions introduced the flux Richardson number Ri f = B/(B + ǫ k ), which is the ratio of buoyancy flux to turbulence production by the mean flow and can be thought of as a mixing efficiency. A related quantity is the mixing coefficient Γ = B/ǫ k , so that K ρ = Γ ǫ k /N 2 . A constant mixing efficiency η = Ri f = 0.17 was assumed by [START_REF] Osborn | Estimates of the local rates of vertical diffusion from dissipation measurements[END_REF] leading to a mixing coefficient Γ = Ri f /(1-Ri f ) = 0.2, a value which has been widely used in oceanographic applications ever since. [START_REF] Salehipour | Diapycnal diffusivity, turbulent Prandtl number and mixing efficiency in Boussinesq stratified turbulence[END_REF] suggested the use of ǫ p instead of the buoyancy flux B when calculating the mixing coefficient, because the buoyancy flux contains both reversible and irreversible exchanges of kinetic and available potential energy, while we are interested in the irreversible conversion of available potential energy into background potential energy, quantified by ǫ p . Throughout the paper we therefore use the definitions Γ = ǫ p /ǫ k for the mixing coefficient and η = ǫ p /(ǫ k + ǫ p ) for the mixing efficiency. Direct support for an eddy diffusivity of the form K ρ = ǫ p /N 2 = Γ ǫ k /N 2 was provided by the work of [START_REF] Lindborg | Vertical dispersion by stratified turbulence[END_REF] who derive an analytical expression for the mean square particle displacement 1/2 δz 2 , which increases linearly in time, the constant of proportionality being K ρ .

The parameters that could conceivably affect mixing in stratified turbulence are the buoyancy Reynolds number Re b , the turbulent Froude number F r and the Reynolds number Re. Taking advantage of the estimate ǫ k ∼ u 3 /ℓ, we define these as

Re b = ǫ k νN 2 , F r = ǫ k N u 2 , Re = u 4 νǫ k , (1.1)
These definitions differentiate themselves from the more classical definitions using a turbulent lengthscale ℓ. Throughout the paper u and ℓ are considered to be horizontal scales so that the results generalize to the case of strongly stratified turbulence and its anisotropic conditions. In this case then F r represents a horizontal Froude number, as opposed to the vertical Froude number F r v = u/(N ℓ v ) based on a vertical lengthscale ℓ v . As shown by [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF], F r v = O(1) in the limit of low F r and therefore F r v does not further influence the dynamics. As pointed out for example by [START_REF] Ivey | Density stratification, turbulence, but how much mixing?[END_REF] the three parameters listed in equation (1.1) are not independent since Re b = Re F r 2 . In problems of mixing through a density discontinuity a relevant parameter is the bulk Richardson number Ri b = g∆ρℓ/(ρ 0 u 2 ) where ∆ρ is the density jump across the interface and ρ 0 is a reference density. If we substitute ∆ρ/ℓ with the continuous stratification density gradient |dρ/dz| it is clear that Ri b ∼ N 2 ℓ 2 /u 2 ∼ F r -2 using the definition of N = -(g/ρ 0 )dρ/dz. Hence problems with a density interface can also effectively be characterized by the parameters given in (1.1).

Since the seminal work of [START_REF] Osborn | Oceanic fine structure[END_REF]; [START_REF] Osborn | Estimates of the local rates of vertical diffusion from dissipation measurements[END_REF] it has been found that mixing efficiency and mixing coefficient are not constants but vary in a certain parameter range. [START_REF] Barry | Measurements of diapycnal diffusivities in stratified fluids[END_REF] found a variation of the mixing coefficient in their grid stirring experiments for buoyancy Reynolds numbers Re b = ǫ k /(νN 2 ) > 300 that was well described by Γ ∝ Re [START_REF] Lozovatsky | Mixing efficiency in natural flows[END_REF]. Recently, [START_REF] Mater | The quest for an unambiguous parameterization of mixing efficiency in stably stratified geophysical flows[END_REF] have suggested that a multi-parameter approach is more suitable to describe mixing processes in stratified turbulence.

Classical parametrizations of mixing have focused on the bulk Richardson number Ri b , which is closely related to the turbulent Froude number since Ri b ∼ F r -2 . In experiments of mixing across a density interface by [START_REF] Turner | The influence of molecular diffusivity on turbulent entrainment across a density interface[END_REF]; [START_REF] Kato | On the penetration of a turbulent layer into stratified fluid[END_REF] the focus was on the entrainment velocity u e across the interface. This is closely related to the mixing efficiency and it can be shown that η = Ri b (u e /u) ∼ F r -2 (u e /u) (see [START_REF] Turner | Buoyancy effects in fluids[END_REF]. In the limit of weak stratification, corresponding to high values of F r, [START_REF] Turner | The influence of molecular diffusivity on turbulent entrainment across a density interface[END_REF] found that the normalized entrainment velocity u e /u reaches a constant, implying the result η ∝ F r -2 for weakly stratified turbulence. This result has also been suggested by [START_REF] Holford | Turbulent mixing in a stratified fluid[END_REF] who state that the buoyancy field behaves as a passive scalar in the limit of high F r. At the other end of the spectrum, strong stratification and low F r lead to an entrainment velocity u e /u ∝ Ri -1 b ∼ F r 2 in the experiments by [START_REF] Kato | On the penetration of a turbulent layer into stratified fluid[END_REF]. This means a constant mixing efficiency η in the limit of strongly stratified turbulence, which has been confirmed by several more recent experiments [START_REF] Park | Turbulent mixing in stratified fluids: layer formation and energetics[END_REF][START_REF] Oglethorpe | Spontaneous layering in stratified turbulent Taylor-Couette flow[END_REF][START_REF] Olsthoorn | Vortex-ring-induced stratified mixing[END_REF].

Hence two parametrizations of mixing exist: a more classical one based on Ri b and therefore essentially on the turbulent Froude number F r and a more recent one based on Re b . We now turn to the evidence from scaling of the equations of motion to try and shed light on this debate.

The Boussinesq set of equations for a linearly stratified fluid is given by:

∂u ∂t + u • ∇u = - 1 ρ 0 ∇p + b e z + ν∇ 2 u, (1.2) ∂b ∂t + u • ∇b = -N 2 u z + D∇ 2 b. (1.3)
We consider first the limit of weak stratification and high F r where horizontal and vertical lengthscales can be assumed to be approximately equal as can be done for horizontal and vertical velocity scales. From the buoyancy equation ( 1.3) a simple balance between advection u • ∇b and the background stratification term N 2 u z leads to a scaling for the buoyancy as b ∼ N 2 ℓ. Comparing this to the advection term in the Boussinesq momentum equation

(1.2) u • ∇u ∼ u 2 /ℓ it is clear that their ratio is (N 2 ℓ)/(u 2 /ℓ) = N 2 ℓ 2 /u 2 ∼ F r -2 .
Hence in the limit of F r ≫ 1, the buoyancy term can be neglected to leading order in the momentum equation and buoyancy effectively behaves as a passive scalar in the presence of a mean scalar gradient N 2 as pointed out by [START_REF] Holford | Turbulent mixing in a stratified fluid[END_REF]. It is possible from equations (1.2)-( 1.3) to form equations for the turbulent kinetic energy and potential energy. At high Reynolds number and high Peclet number, the kinetic and potential energy dissipation rates reach a finite and positive limit that is of the same order as the advection term in these energy equations:

ǫ k ∼ u • ∇ |u| 2 2 ∼ u 3 ℓ , (1.4) ǫ p ∼ u • ∇ b 2 2N 2 ∼ b 2 u N 2 ℓ ∼ N 2 uℓ. (1.5)
From these scalings it follows that the mixing coefficient

Γ = ǫ p /ǫ k ∼ (N 2 uℓ)/(u 3 /ℓ) = N 2 ℓ 2 /u 2 ∼ F r -2 in the limit of high F r.
This result is analogous to that obtained when considering Turner's experiment. The Prandtl number P r = ν/D does not affect the analysis, because we assume that the Reynolds number Re and Peclet number P e = Re P r are both sufficiently high that the kinetic and potential energy dissipation rates become independent of viscosity ν and diffusivity D, respectively. The strongly stratified turbulence theory was developed by several researchers in this field; its two main conditions are F r ≪ 1 and Re b ≫ 1 and these conditions both have to be met within the theoretical framework (developed by [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]). The scaling analyses behind this work have proven to be able to describe the layered large-scale appearance of strongly stratified turbulence, observed previously by many authors and reviewed extensively by [START_REF] Riley | Fluid motions in the presence of strong stable stratification[END_REF]. A consensus has emerged recently that for values Re b > 10 stratified turbulence has the expected form of the horizontal energy spectrum

E h (k h ) ∼ ǫ 2/3 k k -5/3 h
and that the vertical lengthscale is approximately equal to ℓ v = u/N . This has been shown in many numerical works in the last ten years [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF] 2007) tested many of the predicted results of the strongly stratified turbulence theory and plotted, for example, the anisotropy of the dissipation S diss against Re b in their forced stratified DNS simulations. This quantity is defined as S diss = ν (∂u x /∂z) 2 + (∂u y /∂z) 2 /ǫ k and represents the ratio of dissipation due to vertical gradients to the overall dissipation. The authors of this work found that S diss rapidly goes towards the isotropic value of 4/15 ≈ 0.267 [START_REF] Taylor | Statistical theory of turbulence[END_REF] as the buoyancy Reynolds number is increased. Similar results were obtained in decaying stratified turbulence by [START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF]; [START_REF] Maffioli | Dynamics of stratified turbulence decaying from a high buoyancy Reynolds number[END_REF] who both found values close to S diss ≈ 0.41 at times when the stratified turbulence in their DNS was most vigorous and Re b = 5 -9. Furthermore, as discussed extensively in [START_REF] Smyth | Length scales of turbulence in stably stratified mixing layers[END_REF]; [START_REF] Smyth | Differential diffusion in breaking Kelvin-Helmholtz Billows[END_REF], in the case of stratified mixing layers a good condition for determining if there is active three-dimensional turbulence at the small scales is that Re b > 20. The issue of whether the isotropic limit of the dissipation terms will be reached in strongly stratified turbulence is still being studied and constitutes an active area of research. Note however that the isotropic limit is not necessary for our main scaling analysis to be applicable because the scaling is based on the classical assumption that the dissipation rates reach a finite and positive limit at high Reynolds number, irrespective of the values of viscosity and diffusivity.

We now proceed to the main part of this paper that is concerned with estimating the dependency of Γ on the important physical parameters of the problem. Considering our novel scaling analysis at high F r and the above arguments, it is our contention that in stratified turbulence the mixing efficiency is described by the Froude number only. The buoyancy Reynolds number Re b should not play any role. To test this theoretical reasoning, we have performed DNS of stratified turbulence across a large range of Froude numbers and Reynolds numbers.

Numerical methodology

We perform direct numerical simulations (DNS) of turbulence in a linearly stratified fluid, with constant Brunt-Väisälä frequency N . The equations which are solved directly are the Boussinesq set of equations (1.2)-(1.3) together with the continuity equation ∇ • u = 0. We have included a body force f in the Boussinesq momentum equation (1.2) to ensure that the turbulence reaches statistical stationarity. These equations 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 (for more details see [START_REF] Rogallo | Numerical experiments in homogeneous turbulence[END_REF].

The turbulence is kept in a steady state through the use of forcing. We use random forcing that is uncorrelated in time, in a way that the physical location of the body force keeps changing in time. Most of the simulations are at high F r, for which we use isotropic forcing in spherical shells of radius k max f = 5. The forcing is mostly concentrated at smaller wavenumbers with a forcing power spectrum

P (k) = α 2 exp(-(k -k peak f ) 2 )
with k peak f = 3. The value of α is determined at every timestep to ensure that the forcing

power P f = all k (û * • f + 1/2| f | 2 ∆t
) is kept to a constant value. Since f ∝ α this results in a quadratic equation for α, which we solve at every timestep. Of the two roots for α we choose the one with minimum absolute value, which therefore minimizes the magnitude of the force; we find this technique reduces the oscillations in time for most quantities of interest, such as ǫ k , ǫ p and RMS values of the velocity components. At low Froude numbers F r < 0.2, we found that isotropic forcing leads to non-stationary solutions, with growth of energy in the shear modes (modes with k x = k y = 0) that quickly dominate the overall kinetic energy. We therefore utilize vortical forcing concentrated in the vertically rotational modes with k z = 0. This type of forcing in modes with k z = 0 and finite and small k h takes a long time to leak energy to the shear modes, with finite and small k z and k h = 0, and hence a quasi-stationary state can be reached at intermediate times.

Vortical forcing has been used by numerous authors when studying stratified turbulence at low F r (see [START_REF] Waite | Stratified turbulence dominated by vortical motion[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF][START_REF] Augier | Stratified turbulence forced with columnar dipoles: numerical study[END_REF].

Most simulations are run on cubic domains except for the three simulations with strongest stratification for which we take advantage of the anisotropy of the lengthscales and use rectangular domains with L z < L h . In all cases, since we expect isotropic dissipation at high Re, the grid spacing is constant in all directions and it is chosen to ensure k max η = 1.5. The Prandtl number in all cases is P r = 1. The list of all simulations that have been performed for this study is given in table 1 together with the important physical parameters for each run. The strongly stratified runs are performed on increasingly large grids in order to meet the condition Re b > 10 as F r is decreased. As a result our most strongly stratified run with F r = 0.02 and Re b = 17 requires a grid of 4096 2 × 1024 collocation points. The values quoted are time averages over the stationary period of each simulation as are all the values plotted in the next section, in which we discuss the results of the various stratified DNS performed.

Results of DNS runs

Variation of mixing coefficient at constant Re b

We have performed simulations at constant buoyancy Reynolds number and at different Froude numbers. This was achieved by using different number of grid points for each simulation and hence varying also the Reynolds number while keeping Re b ≈ 1000. This value of Re b is well within the energetic-regime of [START_REF] Barry | Measurements of diapycnal diffusivities in stratified fluids[END_REF]; [START_REF] Shih | Parameterization of turbulent fluxes and scales using homogeneous sheared stably stratified turbulence simulations[END_REF] so that the mixing coefficient is expected to not be constant. As shown in figure 1 indeed find a variation of Γ across the five runs but this can now not be due to variations in Re b , which is kept constant, but rather is due to a changing F r. The values of Γ span an order of magnitude as do the values of F r in the simulations. At F r = 0.29 we have Γ = 0.51, a high value compared to the often quoted Γ = 0.2 value [START_REF] Osborn | Estimates of the local rates of vertical diffusion from dissipation measurements[END_REF]. At the other end of the F r-range, F r = 2.9 (but at the same Re b ), we have Γ = 0.05. Such a low value of the mixing coefficient is a result of the weak stratification meaning that there is not much of a density difference for the turbulent flow to mix.

Mixing coefficient at varying Reynolds numbers

The kinetic and potential energy dissipation rates reach the finite and positive limits given in equations (1.4)-(1.5) at high-Re λ , the Reynolds number based on the Taylor microscale, defined as Re λ = uλ/ν with λ = (u 2 / < (∂u x /∂x) 2 >) 1/2 . This was shown by [START_REF] Donzis | Scalar dissipation rate and dissipative anomaly in isotropic turbulence[END_REF] for turbulence with a passive scalar and is expected to be true also for stratified turbulence. If Re λ is not high enough both the non-dimensional dissipations become a function of the Reynolds number, that is ǫ k ℓ/u 3 = f (Re λ ) and ǫ p ℓ/(ub 2 /N 2 ) = g(Re λ ), and these functions are given in [START_REF] Donzis | Scalar dissipation rate and dissipative anomaly in isotropic turbulence[END_REF] in the case of passive scalar advection. Because the result Γ ∝ F r -2 follows from the dissipation scalings in equations (1.4)-(1.5), at finite Reynolds number we expect instead Γ = h(F r, Re λ ) and the mixing efficiency now depends not only on the Froude number but also on the Taylor scale Reynolds number.

To confirm these ideas we consider first the non-dimensionalized dissipation rates given in figure 2 for four runs at constant F r ≈ 3 but performed on successively larger grids so as to increase Re λ from low values to Re λ ≈ 240. This value is considered to be a good high value above which both dissipations reach a constant limit when nondimensionalized (see [START_REF] Donzis | Scalar dissipation rate and dissipative anomaly in isotropic turbulence[END_REF][START_REF] De Bruyn Kops | Classical scaling and intermittency in strongly stratified Boussinesq turbulence[END_REF]. In making the dissipations non-dimensional we use the standard isotropic definition of the turbulent lengthscale, ℓ = (π/u 2 ) (E(k)/k)dk based on the 3-D energy spectrum E(k), expected to be valid at this high F r. We find a non-dimensionalized kinetic energy dissipation that is remarkably close to the form found by [START_REF] Donzis | Scalar dissipation rate and dissipative anomaly in isotropic turbulence[END_REF] in their collection of DNS (the exact same constants as in their paper are used in the curve shown in figure 2(a)). Also the nondimensionalized potential energy dissipation is similar to the prediction of [START_REF] Donzis | Scalar dissipation rate and dissipative anomaly in isotropic turbulence[END_REF] for the non-dimensionalized dissipation of passive scalar variance, confirming that buoyancy behaves as a passive scalar at high Froude numbers. The next step is to check that, for these four simulations at small to moderate Re λ , also the mixing coefficient Γ changes, approaching a constant at the highest Re λ . The steady-state value of Γ obtained for the four runs is given in figure 3, and there is clearly a significant variation, with Γ decreasing for increasing values of Re λ and, possibly, the approach of a constant value Γ ≈ 0.033 at the highest Re λ = 242.

Mixing coefficient at high Re λ

We now consider only the runs with Re λ > 200, for which the dissipation rates ǫ k and ǫ p have approximately reached their respective finite and positive limits.

In figure 4 we show Γ as a function of F r for these runs. If we focus on the high-F r behaviour we see that indeed Γ ∝ F r -2 for F r > 1, which confirms the scaling arguments delineated in §1. These scaling arguments have as a first result E P /E K ∼ E P /u 2 ∝ F r -2 and this is also confirmed by the simulations for F r > 1 as shown in figure 5. Of course our simulations are in a stationary state, which makes them easier to analyse, while time-evolving simulations such as those of [START_REF] Shih | Parameterization of turbulent fluxes and scales using homogeneous sheared stably stratified turbulence simulations[END_REF]; [START_REF] Salehipour | Diapycnal diffusivity, turbulent Prandtl number and mixing efficiency in Boussinesq stratified turbulence[END_REF] have an extra degree of freedom. It may be that some definition of the Reynolds number may be important to describe the decay of the turbulence, or the "age of a mixing event" as described by [START_REF] Smyth | Length scales of turbulence in stably stratified mixing layers[END_REF], in this case. However, at F r > 1 the relevant Reynolds number should not be the buoyancy Reynolds number. This follows from a simple argument, the Ozmidov scale ℓ oz ∼ F r 3/2 ℓ > ℓ ≈ L box for F r > 1. The buoyancy Reynolds number is a measure of the dynamic range from the Kolmogorov scale through to ℓ oz since it can be written as Re b = (ℓ oz /η) 4/3 . However, the Ozmidov scale is now larger than ℓ, the largest physical scale in the problem, and for this reason Re b does not have a physical influence on the stratified turbulence being simulated.

In the limit of low F r and strong stratification, we expect the mixing coefficient to reach a constant value as in the experiments by [START_REF] Park | Turbulent mixing in stratified fluids: layer formation and energetics[END_REF]; [START_REF] Oglethorpe | Spontaneous layering in stratified turbulent Taylor-Couette flow[END_REF]; [START_REF] Olsthoorn | Vortex-ring-induced stratified mixing[END_REF] as long as Re b > 10. We have chosen this criterion because there is evidence that for Re b = O(10) the turbulence starts showing the characteristics of strongly-stratified turbulence as discussed in §1. In figure 4 we see that Γ has a maximum around F r = 0.3 and then drops as the stratification is increased. Our results confirm that Γ indeed approaches a constant value Γ ≈ 0.33 as F r → 0. The ratio E P /u 2 is also predicted to be a constant of order unity by the strongly stratified turbulence theory (see [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF]. From figure 5 this ratio appears to slowly go to a constant value E P /u 2 ≈ 0.15 for F r < 0.3.

Discussion and conclusions

There is a revealing relationship linking F r and Re b in the limit of high Reynolds number when the dissipative scales are isotropic and ǫ k = 15ν < (∂u x /∂x) 2 >: equation (4.1) to find that the corresponding turbulent Froude number is F r = 0.45, which is close to the value of F r = 0.3 at which we observed the drop in Γ (see figure 4). We can take this a step further and try to estimate the value of F r at which Γ starts varying in ocean and atmosphere field data. [START_REF] Davis | The modification of bottom boundary layer turbulence and mixing by internal waves shoaling on a barrier reef[END_REF] report turbulence Reynolds numbers greater Re ℓ = 1500 in their ocean measurements. Using Re ℓ = uℓ/ν ∼ u 4 /(νǫ k ) = 15Re 2 λ we can make the estimate Re λ ≈ √ 15Re = 150, not much larger than the value from the DNS of [START_REF] Shih | Parameterization of turbulent fluxes and scales using homogeneous sheared stably stratified turbulence simulations[END_REF]. The variation in Γ was observed for Re b > 100 or for F r > 15Re b /Re λ = 0.26, again close to our transition Froude number.

We finally consider the atmospheric data of [START_REF] Lozovatsky | Mixing efficiency in natural flows[END_REF]. From the data of the VTMX experiment in [START_REF] Monti | Observations of flow and turbulence in the nocturnal boundary layer over a slope[END_REF], night-time values for the horizontal turbulent velocity were u ≈ 0.3 m/s, while a value for the dissipation can be inferred from the compensated spectra in Figure 4(b) of [START_REF] Lozovatsky | Mixing efficiency in natural flows[END_REF] as ǫ k ≈ 1.8 × 10 -3 m 2 /s 3 . Using ν = 1.4 × 10 -5 m 2 /s for air at 10 • C one obtains a value for the Taylor microscale as λ = 0.1 m and Re λ ≈ 2000. This finally gives a transition Froude number as F r = 0.39 corresponding to Re b = 40,000. So for completely different transition values based on buoyancy Reynolds number (to which we can add Re b = 1000 from the present DNS) we have similar value of turbulent Froude number close to F r ≈ 0.3 at which the mixing coefficient starts dropping considerably. The variation in Γ does therefore appear to be a Froude number effect. It remains to clearly show Γ ∝ F r -2 in field data, which is essential in our eyes to finish resolving this open issue in stratified turbulent mixing. The first evidence in this direction comes from the work of [START_REF] Wells | The relationship between flux coefficient and entrainment ratio in density currents[END_REF] who find Γ ∼ F r -2 at high F r in the case of oceanic gravity currents.

In conclusion, we have presented results from direct numerical simulations of constant-N forced stratified turbulence covering almost 3 orders of magnitude in F r and a vast range of Re b . The simulations at high Re λ show a clear behaviour of the mixing coefficient

.

  In their DNS of stratified shear flow,[START_REF] Shih | Parameterization of turbulent fluxes and scales using homogeneous sheared stably stratified turbulence simulations[END_REF] found a constant mixing coefficient Γ ≈ 0.2 for 7 < Re b < 100 but a varying mixing coefficient going as Γ ∝ Re -1/2 b for Re b > 100. Ocean field measurements by Davis & Moninsmith (2011); Walter et al. (2014) have found similar variations of Γ ∝ Re -1/2 b at high Re b > 100. Atmospheric boundary layer measurements within the vertical transport and mixing experiment (VTMX) also have a similar variation of Γ with Re b but now at Re b > 10 4 , suggesting a completely different bound on the buoyancy Reynolds number (see

Figure 1 .

 1 Figure 1. Mixing coefficient as a function of Froude number at a constant Re b ≈ 1000. Each point corresponds to one of the five simulations labelled R1kF2.9 through to R1kF0.3.

Figure 2 .Figure 3 .

 23 Figure 2. (a) Kinetic energy dissipation and (b) potential energy dissipation, both non-dimensionalized, for 4 runs with F r ≈ 3. The relevant DNS runs are R1kF2.9, 256F2.9, 512F3 and 1024F3.1. The curves represent the fit proposed by Donzis et al. (2005) with the same constants in the case of ǫ k ℓ/u 3 , and the same value for D but a slightly different value for C in the case of ǫpℓ/(ub 2 /N 2 ) (we use C = 0.31).

Figure 4 .

 4 Figure 4. Mixing coefficient as a function of Froude number (log-log plot shown in inset).

Figure 5 .

 5 Figure 5. Ratio of potential energy to horizontal kinetic energy as a function of Froude number (log-log plot shown in inset).

Table 1 .

 1 List of DNS runs performed: relevant non-dimensional parameters and type of forcing.

	Run	Nx = Ny Nz	F r	Re b Re λ Re forcing
	R1kF2.9	96	96 2.90 1010	42	120	iso
	R1kF1.6	192	192 1.64	990	74	370	iso
	R1kF0.9	384	384 0.94	980	129 1120	iso
	R1kF0.5	768	768 0.52	960	229 3480	iso
	R1kF0.3	1536	1536 0.29	990	423 11930	iso
	1024F0.7	1024	1024 0.70 2340 266 4730	iso
	1024F1.6	1024	1024 1.58 10430 250 4180	iso
	1024F3.1	1024	1024 3.10 37370 242 3900	iso
	1024F5.9	1024	1024 5.86 133430 241 3880	iso
	1024F12	1024	1024 11.97 537250 237 3750	iso
	256F2.9	256	256 2.85 4190	88	520	iso
	512F3	512	512 3.03 13370 148 1460	iso
	R200F0.14	1024	1024 0.141 200	390 10130 vort
	R57F0.09	1024	1024 0.091	57	319 6760	vort
	R14F0.04	1024	512 0.044	14	324 7000	vort
	R15F0.03	2048	512 0.035	15	432 12420 vort
	R17F0.02	4096	1024 0.020	17	805 43180 vort

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

We would like to thank P.K. Yeung for providing us with the DNS code used for this study. We thank Pierre Augier for helpful discussions on mixing efficiency. We are indebted to the Leverhulme Trust network on Waves & Turbulence headed by Peter Davidson for the many discussions it generated between its participants. We acknowledge PRACE for awarding us via the StratForce project access to the Fermi cluster at CINECA in Italy. Computational resources at PDC in Sweden were made available by SNIC. G.B. acknowledges financial support by the Swedish Research Council (grant number 621-2013-5784).

as Γ ∝ F r -2 for F r > 1, confirming the scaling analysis for weakly stratified turbulence. In the opposite limit, F r < 1 a peak in Γ is found at F r ≈ 0.3 with a high value of mixing coefficient Γ = 0.51, significantly larger than the Γ = 0.2 value that is commonly used in oceanographic applications. The mixing coefficient then drops to values around Γ = 0.33; these are still high values and are due to the presence of a strong stratification, which high-Re b turbulence is able to mix efficiently. To what extent our results can be generalized to all types of flows is an open question. We find it very likely that the general dependence of Γ on F r in the two limits F r ≪ 1 and F r ≫ 1 will show some degree of universality.