A V E N A I L L E

L G O S T I A U X 2 A N

D J S O M M E R

The role of fluctuations across a stable density interface

A statistical mechanics theory for a fluid stratified in density is presented. The predicted statistical equilibrium state is the most probable outcome of turbulent stirring. It results from a competition between turbulent transport and sedimentation of fluid particles by buoyancy effect. An approximate equipartition between kinetic and available potential energy is obtained. The slow temporal evolution of the vertical density profile is then related to the presence of irreversible mixing, which alters the global distribution of density levels. We propose a model in which the vertical density profile evolves through a sequence of statistical equilibrium states.

The theory is then tested with laboratory experiments in a two-layer stably stratified fluid forced from below by an oscillating grid. The turbulence produced by the grid spreads upwards in the lower region and is blocked at the interface between the dense (salty) and the light (fresh) water. The interface slowly moves upward by entraining fresh water in the turbulent region. Quantitative measurements of density fluctuations are made by planar laser induced fluorescence.

The density fluctuations across the interface are splitted in a "wave" part and a "turbulent" part. Temporal and spatial spectra of the wave part of the density field are well described by a previous theory due to Phillips.We argue that statistical mechanics predictions apply for the turbulent part of the density field. Assuming a two level global density distribution, the theory predicts a hyperbolic tangent shape for the mean vertical density profile, in agreement with experimental observations. The theory predicts that the interface thickness is inversely proportional to the Richardson number, which is also consistent with experimental observations. The density fluctuations obtained after removal of the wave part fit well with the statistical equilibrium theory in the interface region. However inside the mixed layer density fluctuations are instead controlled by a balance between eddy flux downward and dissipation by cascade to small scales. We report exponential tails for the density pdf in this region, similar to previously observed temperature pdf in high Rayleigh number convection experiments.

Introduction

A remarkable property of strongly stratified turbulent flows is their propensity to form density staircases with relatively thin interfaces separating regions of homogenised density [START_REF] Ruddick | The formation of layers in a uniformly stirred density gradient[END_REF][START_REF] Park | Turbulent mixing in stratified fluids: layer formation and energetics[END_REF]. Such staircases are routinely observed in the thermocline structure of the ocean. I that case, they often can be attributed to double diffusive convective instabilities involving temperature and salt [START_REF] Schmitt | Double diffusion in oceanography[END_REF]. However, staircases also occur in situations where the density variations are due solely to the temperature, as for instance in lakes [START_REF] Simpson | Temperature microstructure in a fresh water thermocline[END_REF], and we will restrict ourself in this paper to cases without double diffusion. The emergence of sharp density interfaces involves small scale mixing process, and yet their presence strongly affects the large scale transport of tracers such as pollutants or nutrients. Understanding what sets the interface shape and its evolution is therefore a challenging problem.

A satisfactory model for the interface shape remains elusive. [START_REF] Phillips | Turbulence in a strongly stratified fluid: is it unstable[END_REF]; [START_REF] Posmentier | The generation of salinity finestructure by vertical diffusion[END_REF] proposed a dynamical equation for the mean vertical density profile by modelling small scale turbulence with an eddy diffusivity depending on the local density gradient. The key assumption was to assume that the eddy diffusivity coefficient decreases sufficiently fast with the local density gradients above a given threshold, so that the eddy flux decreases with increasing density gradient. This approach explains the emergence of a sharp interface when stirring an initial stably (strongly) stratified linear density profile, but also predicts that the interface thickness becomes infinitely small (only limited by molecular diffusivity). By contrast, the density interfaces observed by [START_REF] Park | Turbulent mixing in stratified fluids: layer formation and energetics[END_REF] (with a permanent stirring) were found to be larger than the diffusive length scale. In order to overcome this difficulty, [START_REF] Balmforth | Dynamics of interfaces and layers in a stratified turbulent fluid[END_REF] proposed a model in which the eddy diffusion of density is coupled to a dynamical equation for the turbulent kinetic energy -also based on eddy viscosity. Under some assumptions on the forcing terms in the energy equation, [START_REF] Balmforth | Dynamics of interfaces and layers in a stratified turbulent fluid[END_REF] obtained solutions for the density profile showing the formation of homogeneous layers separated by sharp but finite interfaces.

Here we propose a radically different, but complementary approach: the emergence of a sharp but finite interface is interpreted as the most probable outcome of turbulent stirring between two regions of homogenised density, just as for the vorticity field in freely evolving two-dimensional turbulence [START_REF] Miller | Statistical mechanics of euler equations in two dimensions[END_REF][START_REF] Robert | Statistical equilibrium states for two-dimensional flows[END_REF]. We show that the most probable state results from a competition between turbulent transport that tends to increase the interface thickness, and buoyancy effects that tend to sharpen the interface by the vertical drift of fluid elements. This approach does not require any eddy diffusivity hypothesis: it relies on the assumption that the system sufficiently explores the phase space available under the energy conservation constraint. We will build upon the work of [START_REF] Tabak | Mixing in simple models for turbulent diffusion[END_REF]; [START_REF] Venaille | Modeling mixing in two-dimensional turbulence and stratified fluids[END_REF], which were to our knowledge the only attempts to generalise the equilibrium statistical mechanics approach to stratified turbulence.

We will use this statistical mechanics theory to interpret experimental observations of density fluctuations across a turbulent density interface. We will consider for that purpose an experimental idealisation of a mixed layer originally proposed by [START_REF] Rouse | Turbulent diffusion across a density discontinuity[END_REF]. An oscillating grid is localised at the bottom of a tank initially filled with a two layer stratified fluid. The grid-generated turbulence spreads in the homogeneous bottom layer and is blocked at the interface, where lighter fluid is slowly entrained in the turbulent layer. Since stirring is present only on the lower part, entrainment induces a slow depletion of the layer at rest and a progression of the mixed layer thickness. Here we distinguish two questions: i/ what sets the interface shape ? ii/ what sets the entrainment across the interface ?

Most previous studies have been focused on the second question, see e.g. [START_REF] Linden | Mixing in stratified fluids[END_REF]; [START_REF] Fernando | Turbulent mixing in stratified fluids[END_REF] for a review. [START_REF] Turner | The influence of molecular diffusivity on turbulent entrainment across a density interface[END_REF] addressed the role of molecular diffusivity on the entrainment velocity. Here we will assume that both Reynolds number and the Peclet numbers are very large. The key parameter of the problem is then the Richardson number Ri = g ′ L t /e c , where L t and e c are the integral length scale and kinetic energy of the turbulent velocity field, and g ′ the reduced gravity at the interface. Depending on this Richardson number, different flow regime have been identified close to the interface, see e.g. [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF] and references therein. Each of these regimes has led to a different entrainment model: i/ At moderate Richardson numbers, entrainment is mostly due to coherent vortices of the mixed layer impinging on the density interface [START_REF] Linden | The interaction of a vortex ring with a sharp density interface: a model for turbulent entrainment[END_REF][START_REF] Sullivan | The penetration of a density interface by heavy vortex rings[END_REF]. ii/ At larger Richardson number, entrainment is controlled by the breaking of Kelvin-Helmholtz instabilities [START_REF] Mory | A model of turbulent mixing across a density interface including the effect of rotation[END_REF]. iii/ At even larger Richardson number, entrainment is dominated by intermittent interfacial wave breaking (Fernando & Hunt 1997).

By contrast with these studies, the equilibrium statistical mechanics approach does not rely on a particular physical mechanism, but rather on the assumption that the flow sufficiently explores the phase space. We will show that this theory provides predictions for the interface shape, but that it must be supplemented by a model for irreversible mixing in order to describe entrainment across the interface. In real flows, there may exist dynamical regimes preventing the relaxation toward an equilibrium state, but the statistical mechanics approach provides at least an interesting limit case that has yet not been explored.

There are only few experimental and theoretical studies dealing with the vertical density profile across the interface, which is most often taken as given for theories on the entrainment velocity. In addition, most studies are limited to the estimate of the interface thickness. [START_REF] Crapper | The structure of turbulent density interfaces[END_REF]; [START_REF] Fernando | On the nature of the entrainment interface of a twolayer fluid subjected to zero-mean-shear turbulence[END_REF] found that the interface thickness did not vary significantly with the Richardson number, but [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF]; [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF] observed that the interface thickness decreases with increasing Richardson number. E & [START_REF] Hopfinger | On mixing across an interface in stably stratified fluid[END_REF] found a similar result, but with a finite asymptotic thickness depending on the Peclet number and the Reynolds number.

Here we present novel measurements of density fluctuations across the interface in these experiments, using planar laser induced fluorescence techniques, with index matching. This technique was already used by [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF]; [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF] in a similar experimental setting, or by [START_REF] Guyez | Turbulent mixing at a stable density interface: the variation of the buoyancy flux gradient relation[END_REF] in the case of two layers stably stratified Taylor-Couette experiments, but these studies did not focus on the role of density fluctuations.

As pointed out by [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF] , the density field of a snapshot and of a time (or spatial) averaged of the density field may be very different. We will address this important issue. In particular, we will show that the density interface can be decomposed into a large scale wave motion whose amplitude is well described by a theory due to [START_REF] Phillips | The dynamic of upper ocean[END_REF], and that removing this wave motion allows us to recover the turbulent density field predicted by statistical mechanics. [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF] already gave experimental evidence for Phillips theory concerning the interfacial wave amplitude. Here we confirm these observations and provide additional support for the theory. We note however that the Phillips approach is limited because it assumes that the velocity close to the interface is entirely due to the interface motion, while turbulent velocities in the mixed layer are in reality of the same order (or even larger) and may play an important role [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF]). We will also show that the mean density profile which is a priori prescribed in the Phillips theory can be accounted for by the statistical mechanics theory.

Importantly, the experiments presented in this paper are performed with a permanent forcing mechanism, which allows us to sustain a quasi-stationary turbulent field. We will interpret this turbulent field as an effective heat bath that sets the level of fluctuation for the density field. We will neglect the direct effect of buoyancy on the turbulent eddies. Although this may be true sufficiently far form the density interface, the interplay between turbulence and stratification may be important in real flows [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF][START_REF] Mcdougall | Measurements of turbulence in a zero-mean-shear mixed layer[END_REF]). In addition, we consider an experimental setting with only one turbulent layer. However, previous work do not show qualitative differences between the single stirred and the double stirred case [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF]. Note that the evolution of the vertical density profiles in the presence of a permanent forcing differs from the studies by [START_REF] Linden | Mixing across a density interface produced by grid turbulence[END_REF]; [START_REF] Whitehead | Turbulent mixing of two-layer stratified fluid[END_REF], who considered the evolution of the mean vertical density profile between isolated mixing events, obtained either by dropping an horizontal grid in the flow, or by moving horizontally a vertical rod, and waiting for the turbulence to decay between each mixing event.

The paper is organised as follows. The equilibrium statistical mechanics theory is introduced and discussed in the second section. It is shown in a third section that one must add to this equilibrium theory a model for irreversible mixing in order to describe entrainment across the interface. The experimental setting is presented in a fourth section. The fifth section contains experimental observations on interfacial waves, and statistical mechanics predictions for the mean vertical density profile. It is also shown in this section that the density fluctuations far from the interface are intrinsically out of equilibrium, and their pdf display exponential tails. We conclude and summarise the main results in the sixth section.

2. Turbulent density interfaces as a statistical equilibrium state 2.1. Boussinesq equations and their conservation laws Consider a flow in the Boussinesq approximation, taking place in a domain V. At each time t and each point x = (x, y, z) ∈ V the system is described by the reduced density b = g (ρ -ρ 0 ) /ρ 0 and the velocity field u. Here ρ is the fluid density, g the gravity and ρ 0 a reference density. Note that with our convention the density b is the opposite of the buoyancy, and this field will be simply referred to as the density field in the remaining of this paper. The velocity field is non-divergent: ∇u = 0.

(2.1)

The flow dynamics is expressed by the advection and molecular diffusion of density b

∂ t b + u∇b = κ∆b, (2.2)
and by the momentum equation

∂ t u + u∇u = -∇P -be z + ν∆u + F, (2.3)
where F is a mechanical forcing and k, l is the vertical unit vector.

We consider here a case without forcing and dissipation (F = 0, κ = ν = 0). If the velocity field remains differentiable within the domain V where the flow takes place, the total energy of the flow

E = V 1 2 u 2 + bz dx (2.4)
is a dynamical invariant. In addition, as a consequence of the density advection equation, the global distribution (i.e. histogram) of density levels

g(σ) = 1 |V| V δ(b -σ)dx (2.5)
is also conserved. Note that g(σ) can be related to the sorted density profile b s (z), obtained by allowing each fluid particle to settle down to its position at rest (lower potential energy), using -g(b s )db s = Adz, where A is the horizontal domain area. Therefore the function g is proportional to the inverse of the derivative of this sorted profile, g(b s ) = -A(db s /dz) -1 . Boussinesq equations admit other dynamical invariants called Casimir functionals, which are related to the conservation of Ertel potential vorticity, see e.g. [START_REF] Salmon | Lectures on geophysical fluid dynamics[END_REF]. However, in this paper, we will not take into account these invariants to compute the equilibrium state. The fact that in the absence of rotation, statistical equilibria are not affected by dynamical invariants related to potential vorticity conservation has been recently discussed in detailed in the framework of the equilibrium statistical mechanics of the shallow water system [START_REF] Renaud | Equilibrium states and energy partition of the shallow water model[END_REF].

Equilibrium states

Here we proceed by analogy with statistical mechanics of two dimensional Euler or quasigeostrophic flows originally proposed by [START_REF] Miller | Statistical mechanics of euler equations in two dimensions[END_REF]; [START_REF] Robert | Statistical equilibrium states for two-dimensional flows[END_REF], see Sommeria ( 2001); [START_REF] Majda | Nonlinear dynamics and statistical theories for basic geophysical flows[END_REF]; [START_REF] Bouchet | Statistical mechanics of two-dimensional and geophysical flows[END_REF]; [START_REF] Lucarini | Mathematical and physical ideas for climate science[END_REF] for recent reviews.

We define a microscopic state of the system as a given fine grained density field b(x, y, z) and a non-divergent velocity field u(x, y, z) = (u, v, w). These microscopic fields are relevant phase space variables because they satisfy a Liouville theorem, stating that the flow in phase space is non-divergent. Then Liouville theorem can be written formally as

V dxdydz δ∂ t b δb + δ∂ t u δu + δ∂ t v δv + δ∂ t w δw = 0 . (2.6)
The proof of such a Liouville theorem for the velocity field (u, v, w) is a classical result for 3D Euler dynamics. It can for instance be obtained by decomposing the velocity field into Fourier modes, see e.g. [START_REF] Bouchet | Statistical mechanics of two-dimensional and geophysical flows[END_REF] and references therein. The generalisation to the Boussinesq system is straightforward, since the only difference with 3d Euler is the presence of an additional term linear in b in the momentum equation, as well as an additional equation describing the pure advection of density b by the nondivergent velocity field.

A macroscopic state of the system is defined by the probability density field ρ(x, σ, v) that describes the probability to measure a given density level b = σ and a given velocity value u = v in the vicinity of the point x.

The conservation laws can be expressed as constraints on these probability field. Indeed, the energy defined in Eq. (2.4) and the global distribution of density levels defined in Eq.

(2.5) can be expressed as a functionals of ρ,

E[ρ] = V dx +∞ -∞ dv +∞ -∞ dσ ρ v 2 2 + σz , (2.7) H σ [ρ] = V dx +∞ -∞ dv +∞ -∞ dσ ρ . (2.8)
Finally, the probability density field ρ is normalised at each point x:

+∞ -∞ dv +∞ -∞ dσ ρ(x, σ, v) = 1 . (2.9)
Just as in the case of the vorticity field in 2D turbulence, we anticipate that a typical microscopic state b, u picked at random for a given set of constraint is characterised by large scale variations superimposed with wild fluctuations at small scales. This is what motivates the introduction of the probability field ρ to describe the field at a macroscopic level: the large scale flow will be obtained by computing the mean quantities

b(x) ≡ +∞ -∞ dv +∞ -∞ dσσρ(x, σ, v) , u(x) ≡ +∞ -∞ dv +∞ -∞ dσvρ(x, σ, v) , (2.10)
while the local fluctuations around this mean flow will be statistically described by the distribution ρ(x, σ, v). Classical counting arguments allow us to estimate the number of microscopic states associated with a given macroscopic field ρ, and to show that the most probable state is the one that maximises the mixing entropy (2.11) among all the states that satisfy the constraints of the problem

S = - V dx +∞ -∞ dv +∞ -∞ dσ ρ ln ρ ,
E[ρ] = E and H σ [ρ] = g(σ)
, see e.g. [START_REF] Miller | Statistical mechanics of euler equations in two dimensions[END_REF]; [START_REF] Robert | Statistical equilibrium states for two-dimensional flows[END_REF]. This variational problem can be written in the compact form:

S (E, g(σ)) = max ρ {S[ρ] | E[ρ] = E, H σ [ρ] = g(σ)} , (2.12)
where S (E, g(σ)) is the equilibrium entropy, and where the maximum is searched over all the density probability fields ρ which are normalised at each point x.

Determination of the equilibrium

The first step to find the equilibrium state is to compute critical points of the variational problem, i.e. to find the field ρ such that first variations of the mixing entropy around this state do vanish, given the constraints of the problem. One needs for that purpose to introduce the Lagrange multipliers β and γ(σ) associated respectively with the energy constraint (2.7) and with the constraints of the global density distribution (2.8), and then compute first variations with respect to ρ: δS -βδE + γ(σ)δH σ dσ = 0 .

(2.13) Equation (2.13) with the normalisation constraint Eq. (2.9) yield the following necessary and sufficient condition for ρ to be a critical point of the variational problem: σ) .

ρ (σ, x, v) = β 2π 3/2 e -β v 2 2 ρ b (z, σ), ρ b (z, σ) ≡ e -βσz+γ(σ) +∞ -∞ dσ e -βσz+γ(
(2.14)

The values of the Lagrange multipliers β and γ(σ) are (implicitly) determined by the expression of the constraints E[ρ] = E and H σ [ρ] = g(σ), given by Eq. (2.7) and Eq.

(2.8) , respectively. Relaxation equations towards these equilibria (for a given global distribution of density level, and for either a prescribed energy and a prescribed inverse temperature) are proposed in [START_REF] Venaille | Modeling mixing in two-dimensional turbulence and stratified fluids[END_REF], with an application to restratification problems. Three remarkable properties are satisfied by the equilibrium state.

First, the probability density field (2.14) is expressed as a product of the probabilities for density σ and velocity v, which means that b and v are two independent quantities at equilibrium.

Second, the predicted velocity distribution is Gaussian, isotropic and homogeneous in space: The local kinetic energy e c of the equilibrium state is therefore homogeneous in space, with

e c = 1 2 dσ dv v 2 ρ = 3 2 β -1 .
(2.15)

The inverse of β defines an effective "temperature" of the turbulent field. According to Eq. (2.15), this temperature is proportional to the variance of the velocity fluctuations. Third, the distribution of density levels ρ b (z, σ) depends only on the vertical coordinate z. Using Eq. (2.14) and the definition of the averaged density b(z) given in Eq. (2.10) , the expression of the distribution of density levels can be written as

ρ b (σ, z) = ρ b (σ, 0)e -β σz- z 0 dz ′ b(z ′ ) , (2.16)
where ρ b (σ, 0) must be determined using the constraints H σ [ρ] = g(σ) given Eq. (2.8). Any moments of the density distribution can be computed using b n = σ n ρdσ.

(2.17)

Applying this expression to the second moment, we can easily show from Eq. (2.14) that the density variance is proportional to the vertical gradient of the mean density, with a coefficient of proportionality given by the inverse temperature:

db dz = β b 2 -b 2 .
(2.18)

To conclude, the kinetic energy of the equilibrium state is homogeneous in space, and the variance of the density fluctuations is proportional to this kinetic energy, with a coefficient proportional to the vertical mean density gradient. It is shown in Appendix A that Eq.

(2.18) implies equipartition between kinetic energy and available potential energy in a low energy limit, which allows to interpret the widely reported mixing efficiency coefficient of 0.25 as a consequence of the rapid relaxation of the system towards statistical equilibrium.

Sharp interfaces as a Fermi-Dirac distribution

Let us consider the particular case of an initial state composed of two density levels in equal proportion. This would for instance be the global distribution of a tank initially filled with two uniform density layers of equal thickness. Note that one can always choose the reference density such that the two levels are symmetric around 0, writing:

σ ∈ {-Σ, Σ} . (2.19)
We introduce the probability p to measure the level Σ at height z:

ρ b = p(z)δ (σ -Σ) + (1 -p(z)) δ (σ + Σ) (2.20)
According to equation (2.14) and using the fact that the two density level are in equal proportions to eliminate γ(σ), we obtain p(z) = e -βΣz e -βΣz + e βΣz , (2.21) (to avoid a constant of integration, we have chosen the origin of the z-axis at the interface between the two layers at rest). This is reminiscent of a Fermi-Dirac distribution. Indeed, for the two level system, the incompressibility constraint plays the same role as the exclusion principle for the statistics of Fermions. Following this analogy, the density field is a collection of fluid particles with energy e p = σz, with a Fermi level ε f = 0, in thermal contact with a heat bath characterised by the inverse temperature β. Using Eq. (2.21) and Eq. (2.15), the mean density profile b = Σp -Σ(1 -p) can be expressed in term of the density jump ∆b = 2Σ and the local kinetic energy e c :

b(z) = - ∆b 2 tanh z ∆h , ∆h ≡ 4e c 3∆b , (2.22)
where we have introduced the interface thickness ∆h. This tanh profile was previously obtained by [START_REF] Tabak | Mixing in simple models for turbulent diffusion[END_REF] using similar arguments. However, [START_REF] Tabak | Mixing in simple models for turbulent diffusion[END_REF] did not relate the inverse temperature β to an effective turbulent heat bath.

According to Eq. (2.22), the Richardson number based on the interface thickness Ri ∆h ≡ 3∆h∆b/(4e c ) is always equal to one. This means that the interface thickness is the typical vertical length of an overturning of the stable interface provided that an order one fraction of the kinetic energy inside the overturning region is transferred into potential energy.

The interface thickness can also be expressed in term of a global Richardson number based on the layer thickness H.

∆h = 2H Ri H , Ri H ≡ H∆b 2e c /3 . (2.23)
Note that the term 3e c /2 is the square of a typical turbulent velocity. The interface is sharp when ∆h ≪ H, i.e. whenever Ri H ≫ 1. In that case the potential energy associated with an overturning of the density field at the domain scale H is much larger than the total kinetic energy. By contrast, the density field becomes homogeneous in the limit for which the total kinetic energy is much larger than the potential energy associated with an overturning of the density field at the domain scale H, in which case Ri H ≪ 1.

More generally, it is clear from Eq. (2.16) that the mean density profile and the statistics of density levels result from a competition between turbulent transport and buoyancy repelling. Turbulent transport is related to the turbulent kinetic energy, which increases when the inverse temperature β decreases. Buoyancy repelling is related to the difference between values of the density levels σ initially present in the flow. In the high energy limit (βσ → 0), the density statistics are homogeneous in space, and there is no stratification.

2.4. Turbulent velocity fluctuations as a "heat bath"

The previous results have been obtained in a microcanonical framework, assuming that the total energy is conserved. However the total energy is not actually conserved. Given a given cut-off length scale, kinetic energy is indeed transferred to scales smaller than this cut-off, no matter how small this cut-off length scale. This has two consequences: i/ in the presence of viscosity, the fluctuations will be dissipated, no matter how small is the viscosity ii/ even in the absence of viscosity, the velocity field may become non-differentiable, hence breaking energy conservation. This contrasts sharply with equilibrium states of the 2D Euler equations, in which case small scale vorticity fluctuations do not contribute to the total energy, which belongs entirely to a large scale flow structure. Therefore a forcing term is needed in order to maintain a statistically steady state. But in that case the system is out-of equilibrium. However, we propose here a phenomenological interpretation of the results obtained in previous subsections by assuming that the observed density field remains close to an equilibrium state, even when the velocity field is strongly out of equilibrium. We assume that just as in the equilibrium case above, the kinetic energy is homogeneous in space and is related to the (dynamical) 'temperature' through Eq. (2.15). This temperature is set by a balance between forcing and dissipation in the momentum equation. Since β is given, the relevant statistical ensemble is the canonical one, and the equilibrium state is the minimiser of the free energy

F [ρ] ≡ -S[ρ] + βE[ρ] .
(2.24)

According to Eq. (2.14), the total mixing entropy can be expressed as the sum of the mixing entropy associated with the density field and the mixing entropy associated with the velocity field: (2.25) and according to Eq. (2.15) the total energy is

S = S b - 3 2 ln β, S b ≡ - V dx +∞ -∞ dσ ρ b ln ρ b ,
E[ρ b ] = E p [ρ b ] + 3 2β , with E p [ρ b ] = V dx +∞ -∞ dσ σzρ b . (2.26)
Since β is given, minimising the functional F defined Eq. ( 2.24) is equivalent to minimising

F b [ρ b ] ≡ -S b [ρ b ] + βE p [ρ b ] . (2.27)
Finally, we obtain a variational problem in the canonical ensemble:

F (β, g(σ)) = min ρ b {-S b [ρ b ] + βE p [ρ b ] | H σ [ρ b ] = g(σ)} , (2.28)
which means that we look for the probability density field ρ b that minimises a free energy while conserving the global distribution of density levels.

To conclude, the density field characterised by its potential energy and its global distribution of density levels can be considered as a subsystem in thermal contact with an effective heat bath provided by the turbulent velocity field. Note that forcing is required to maintain this turbulent velocity field, but we still assume that no forcing and no dissipation is present in the dynamics of the density field. The effect of including dissipation in the density dynamics is discussed in the next section.

Entrainment and irreversible mixing

In the previous section, it was assumed that the global distribution of density levels is conserved. There is in that case no temporal evolution of the vertical mean density profile once the equilibrium state is reached. This is because the statistical mechanics approach does not take into account irreversible mixing through turbulent cascade. This irreversible mixing process changes the global distribution of density over time. If this global distribution of density levels evolves on a sufficiently slow time scale, one may assume that this evolution occurs through a sequence of equilibrium states. We show in Fig. 1 the sequence of vertical density profiles in the case of a two level configuration b ∈]-Σ, Σ[ with decreasing values of Σ for a fixed inverse temperature β = 3/(2e c ) (prescribed by a turbulent heat bath). This shows a trend towards complete homogenisation of the density field with a persistence of the density interface.

Let us assume that the density field is anti-symmetric (in a statistical sense) with respect to an interface located at z = 0, which is the case for instance if the initial condition is made of two layers of homogeneous fluid with equal depth H. We define the entrainment velocity as the relative temporal variation of the averaged density in the lower layer:

U e = - H 2 d t b b , with b ≡ 1 H 0 -H dz b . (3.1)
This definition is consistent with Eq. (1) in [START_REF] Turner | The influence of molecular diffusivity on turbulent entrainment across a density interface[END_REF].

A simplified model in the two-level case

Let us now come back to the two level configuration b ∈] -Σ, Σ[ in order to devise a simple model for irreversible mixing, assuming that the dynamics goes through a sequence of equilibria. It amounts to find a dynamical equation for the level Σ(t).

We propose to model the temporal evolution of the density level Σ as a simple linear relaxation process towards the averaged density in the lower layer introduced Eq. (3.1):

∂ t Σ = -s (Σ -b ) , (3.2)
where s is a mixing rate, i.e. the inverse of a relaxation time, which can be interpreted as a typical stretching time or a cascade rate. The physical motivation for this model is that the interface acts as a barrier for irreversible mixing, in such a way that the turbulence tends to homogenise the fluid independently in each layer. Considering that b is given by the equilibrium profile Eq. (2.22), using Eq. (2.23), and taking the limit of a sharp interface Ri H ≫ 1, Eq. (3.2) becomes

∂ t Σ Σ = s log 2 Ri H + o Ri -1 H . (3.3) Still in the high Richardson limit (Ri H ≫ 1), we get ∂ t b / b = ∂ t Σ/Σ + o Ri -1
H , which, using Eq. ( 3.3) and Eq. (3.1), yields

U e = s H Ri H log 2 2 + o Ri -1 H . (3.4)
Assuming that the velocity field is not affected by stratification, the straining rate s can be obtained on dimensional ground as ∝ e 1/2 c /L t where e c and L t are the turbulent kinetic energy and length scale in the absence of stratification. We recover in that case the classical result U e ∼ Ri -1 H initially proposed by [START_REF] Rouse | Turbulent diffusion across a density discontinuity[END_REF] who assumed first the presence of a sharp interface between two homogeneous layers, and second that the increase of the potential energy is proportional to the energy production by mechanical stirring. However, many experimental observations suggest that for very large values of Ri H the power relation between entrainment velocity and Richardson number is steeper that -1. Different arguments have been proposed to account for these observed power law, see e.g. [START_REF] Fernando | Turbulent mixing in stratified fluids[END_REF] and references therein. The simplified model presented above translates the scaling for the entrainment velocity U e ∼ Ri -n H into a scaling for the cascade rate s ∼ Ri -n+1 H .

To conclude, we have devised a toy model in which the mean vertical profile evolves through a sequence of equilibrium states described by a tanh profile until the flow is fully homogenised. The main caveat of the model is that it assumes a two level distribution for the global density distribution, while mixing through turbulent cascade leads to the creation of a continuum of density levels between its extremal values. This aspect will be discussed in more details in the experimental part of this paper.

Comparison with a model based on turbulent diffusion

In order to appreciate the difference between the model proposed in the previous subsection and other approaches based on an effective turbulent diffusivity, it is instructive to consider the simple case of the temporal evolution of an initial step function through the heat equation ∂ t b = K∂ zz b with an homogeneous diffusion coefficient K and no density flux at the upper and lower boundary ∂ z b| z=0,H = 0. The sequence of vertical density profiles from the initial condition to the final homogeneous state is shown on the right panel of Fig. 1. We clearly see that the route toward complete homogenisation is different The value of the density level Σ(t) is decreasing from 1 (plain green line) to 0 (dashed red line). b) purely diffusive relaxation of an initial step function towards an homogeneous density profile. The diffusion coefficient is homogeneous in space.

in the diffusive case and in the quasi-equilibrium case, for which the interface thickness remains quasi-constant through the homogenisation process.

Density interfaces are sometimes fitted with error functions, see e.g. [START_REF] Crapper | The structure of turbulent density interfaces[END_REF]; [START_REF] Linden | Mixing across a density interface produced by grid turbulence[END_REF]; [START_REF] Whitehead | Turbulent mixing of two-layer stratified fluid[END_REF]. The error function is the solution of the heat equation for a constant diffusion coefficient, in the case of an initial step function in an unbounded domain. In the case a bounded domain, this error function is a good fit for the density profile as long as the interface thickness remains much smaller than the domain size. With a proper rescaling of the z axis and of the density axis, the error function and the hyperbolic tangent functions are hardly discernible. We note however that the physical mechanisms underlying the choice of one function rather than the other to fit experimental data are drastically different. Indeed, the choice of an error function result from a model based on a local turbulent diffusivity hypothesis. By contrast, there is no such assumption required for the choice of a tanh profile: the density profile is interpreted in that case as the equilibrium state of a two level system, which results from the competition between turbulent transport and buoyancy.

Two-layer stratified fluid forced by an oscillating grid

Experimental setting

A tank with horizontal cross section 40 × 40 cm 2 is filled with a layer of dense fluid below a layer of light fluid, see Fig. 2-a. The density is homogeneous in each layer, and each layer depth is initially around 40 cm . This experimental setting is similar to the one described in [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF].The novelty comes from measurements techniques. We use Planar LASER Induced Fluorescence (PLIF) with index matching between both layer in order to observe quantitatively density fluctuations. The lower layer contains water, salt and rhodamine. The upper layer contains water and ethanol, such that the optical index is the same in each layer. The density difference is imposed by the concentration in salt and ethanol in both layers.

Turbulence is generated by an horizontal grid oscillating vertically at 5 Hz, with a grid mesh of 10 cm (including the 2 cm thickness of the grid bars), and an amplitude of 8 cm (crest to crest). The same forcing is used for all the experiment. Turbulence properties have been characterised using PIV measurements in a case without stratification. We observed in that case an exponential decay of the kinetic energy, with an e-folding depth L t = 10 cm interpreted as the integral length scale of turbulence, which is of the order of the grid mesh, see Fig. 2-b and Appendix B. We see that typical turbulent velocities close to the density interface are of the order of U ∼ 1 cm.s -1 . This corresponds to a Reynolds number Re = L t U/ν ≈ 10 3 associated with moderate turbulence close to the interface at the beginning of the experiment. This means that viscous effect may be important once the interface has moved up by around 10 cm. Note however that deeper in the mixed layer, the Reynolds number increases by two order of magnitude. The Peclet number is P e = L t U/κ ≈ 10 6 , with κ the salt diffusivity, also of the same order for alcohol and Rhodamine.

The main control parameter is the density jump ∆ρ/ρ varying from 0.01% to 0.8%, see Tab. 1. The initial interface height denoted H is slightly different from one experiment to another. Therefore, each experiment is characterised by two non-dimensional parameters, namely, the bulk Richardson number Ri and the Richardson number based on the interface height Ri H , respectively defined by

Ri = L t ∆b 2e t /3 , and Ri H = H∆b 2e c /3 , (4.1)
where H is the initial interface height, e c the turbulence kinetic energy measured at z = H in the homogeneous case and L t the integral length scale of turbulence in the homogeneous case. We see Tab. 1 that the bulk Richardson number varies from 1 to 150. In practice, a well defined, sharp interface was only observed for Ri ≥ 10. For lower bulk Richardson number, the density field did not reach a quasi-stationary state presenting a turbulent density interface. For this reason, we will mostly focus on experiments characterised by Ri > 10 in order to test statistical mechanics predictions.

A snapshot of the density field is shown Fig. 3-a. The density field is observed in a 25 × 25 cm 2 frame centred 10 cm above the initial density interface at the beginning of the experiment, in the central part of the tank, see Fig. 2. For each density snapshot, light adsorption is compensated, as well as the presence of possible imperfections in the LASER sheet, as dark bands due to bubbles or dust in the optical path.

Each experiment is performed during 800 second, during which a snapshot of the density field is recorded every second. The grid oscillation starts after t = 5 seconds, and stopped at t = 700 seconds. Then the relaxation to rest is observed. The temporal evolution of the x-averaged density is shown in Fig. 4, we see qualitatively that the interface associated with the x-averaged density remains sharp, and its thickness increases with decreasing Richardson numbers.

Relation with the statistical mechanics model

Several assumptions are required to interpret this experiment in the framework of the equilibrium statistical mechanics theory introduced in section 2. The source of kinetic energy is localised at the grid position in the experiment, implying a vertical decay of the kinetic energy. This contrasts with the statistical equilibrium stating that the kinetic energy (or the effective temperature) is homogeneous in space. Our working hypothesis is that prediction from equilibrium statistical mechanics for the density field may be applied near the density interface, by considering that the effective temperature of the equilibrium state is provided by velocity fluctuations that would be observed at the interface in the absence of stratification.

A second difficulty is that there must be sufficient mixing in phase space for the system to reach the equilibrium state. It is clear that the interface motion at large scale is dominated by waves, for which nonlinear effects driving this mixing are inhibited. We shall therefore assume that the turbulent fluctuations involved in the statistical equilibrium are limited to small scales, after exclusion of a "wave motion " of the density interface.

A third limitation for the applicability of the theory is the irreversible dissipation of density fluctuations through turbulent cascade, changing the global distribution of density levels with time. As discussed in section 3, this effect is actually related to entrainment across the interface. We will assume that the global distribution of density levels evolves on a time scale longer than the one required for the system to reach the equilibrium state.

A fourth limitation of the present experiment comes from the asymmetric forcing: only the lower layer is turbulent. A direct consequence of this asymmetric forcing is that the turbulent layer is actually penetrating into the quiescent layer: contrary to the homogeneous case discussed in previous section, the interface is shifting vertically in the experiment. Testing equilibrium statistical mechanics predictions requires in that case to be in the reference frame moving with the interface, assuming that the fluid evolves through a sequence of stationary states. It is thus necessary to define precisely what is the interface height in the experiment.

Definition of the interface height

One major difficulty associated with the definition of the interface height and thickness stems from the fact that instantaneous density profiles may be very different from temporal or spatial averages, see e.g. [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF]. In order to define the interface height at each time t and location x , [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF][START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF] considered iso-density contours parametrised by

α ≡ ρ -ρ min ρ max -ρ min , (4.2)
and defined the interface height h(x, t) as the height of the contour α = 1/2 As noticed by [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF]; [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF], this method cannot be applied when interface overturns. This occurs for instance with wave breaking. As an example, one clearly sees on the snapshot of Fig. 3-a that iso-density contours do not define singledvalued functions for the interface height h(x, t).

We propose and discuss in this paper a method to define the location of a corrugated interface, which is well defined even in the presence of overturning events in he density field. This method allows us to distinguish interfacial waves form turbulent fluctuations around the interface. In order to find the interface height h(x, t) on the density field snapshot Fig. 3-a, a "sorted density field" is computed Fig. 3-b: at each horizontal point x, the vertical sequence of n z pixels is sorted so that the density is decreasing with increasing height z. As a consequence of this sorting procedure, iso-density contours are always single-valued functions of the horizontal x coordinate. The interface h(x, t) is then defined as the height of the intermediate density contour α = 1/2 of the sorted density field, which is represented as a thick black line in Fig. 3-a,b.

Estimate of the entrainment velocity and test of PLIF calibration

The spatial and temporal variability of the interface height h(x, t) will be explored in next section. Here we consider the slow temporal evolution of the x-averaged interface height h(t) shown Fig. 5-a. We assume that there is a time scale separation between a fast temporal variability of the interface height and a slow penetration of the turbulent layer into the layer at rest. The corresponding entrainment velocity U e is estimated by considering a linear fit of the interface height h(t) between t = 100 and t = 300 seconds for each experiment.

The PLIF calibration is checked on Fig. 5-a by comparing the temporal evolution of the x-averaged interface h(t) with the estimate h 0 < b 0 > / < b >, where < b > (t) is the averaged density in the turbulent layer, below the interface, and h 0 , < b 0 > are the initial interface height and density in the lowest layer. Despite the temporal variability of h(t) associated with interfacial waves, there is a good agreements between the estimate of the interface elevation h(t) using the sorting algorithm, and the estimate of the interface elevation h 0 < b 0 > / < b > using mass conservation. A systematic drift is observed only in the strongly stratified case, probably due to light absorption that was not completely cancelled with our data analysis procedure for this particular experiment.

The entrainment coefficient E = U e /(2e c /3) 1/2 is defined as the ratio between the entrainment velocity and the rms turbulent velocity. This rms turbulent velocity were estimated at the interface height z = h, using PIV measurements in a case without stratification. The variations of the entrainment coefficient E with the bulk Richardson number Ri defined Eq. (4.1) is plotted Fig. 5-b. We find a power law E ∼ Ri -n with the exponent n between 1 and 3/2, consistently with previous observations in the same range of Richardson numbers, either in similar experimental setting [START_REF] Fernando | Turbulent mixing in stratified fluids[END_REF], or in Taylor-Couette experiments [START_REF] Guyez | Turbulent mixing at a stable density interface: the variation of the buoyancy flux gradient relation[END_REF]). In the remaining of this paper, we do not further investigate entrainment mechanism, which would require more detailed measurements of the velocity field close to the interface in order to better characterise vertical density fluxes. We rather assume that the system evolves slowly through a sequence of statistically steady states, and we focus on the properties of the density fields associated with these states.

Characterisation of interfacial waves, interface shape and density fluctuations across the interface

On the one hand, it is clear from Fig. 3 that the interface is corrugated and present small scale structures, consistent with the statistical mechanics approach. On the other hand, low frequency oscillations of the interface height are clearly visible in Fig. 4. These low frequency oscillations cannot be explained with the statistical mechanics approach. We show in the first subsection that some properties of the observed interfacial waves are well described by a heuristic theory due to [START_REF] Phillips | The dynamic of upper ocean[END_REF]. The characterisation of interfacial waves will allow us to propose a criterion to distinguish a "wave part" and a "turbulent part" for the fluctuations of the density field close to the interface. The "turbulent part" of the density field will be considered in the second subsection in order to test statistical mechanics predictions for the interface shape. Finally, we show in the last subsection that density fluctuations within the mixed layer and sufficiently far from the interface are much larger than expected from equilibrium statistical mechanics arguments, and present exponential tails.

Interfacial waves

The spectrum of the interface elevation in a two layer fluid subject to an external forcing was predicted by [START_REF] Phillips | The dynamic of upper ocean[END_REF] with heuristic arguments. In this framework, the predicted spectrum does not depend on the forcing mechanism. The only input of the theory is the interface thickness ∆h. We give in the following some experimental evidence for Phillips' theory, assuming that the interface thickness is the one predicted by statistical mechanics, and also discuss limitation of this approach.

Dispersion relation for interfacial waves

Phillips (1977) considered a sharp density interface separating two homogeneous layers of height H/2 ≫ ∆h, with a density jump ∆b = g∆ρ/ρ. The buoyancy frequency inside the interface can be estimated as N ∼ ∆b ∆h .

(5.1)

Let us first assume that the interface elevation is a monochromatic wave characterised by the frequency ω and the wavenumber modulus K. Let us consider in addition that the wavenumber modulus is such that K ≪ 2π/(∆h), so that the interface can be considered infinitely thin at lowest order, and K ≥ 2π/H, so that the limit of deep water can be considered. Under these assumptions, the dispersion relation is

ω = ∆bK 2 for N ≪ ω ≤ ∆bπ H . (5.2)
The presence of interfacial waves is revealed Fig. 6 by the spatial-temporal spectrum of the interface elevation, denoted ĥ(f, k), where h(x, t) is the observed interface elevation defined with the sorting algorithm introduced in subsection 4.3.

For each frequency f = ω/2π, the spectra shown in Fig. 6 have been normalised by their maximum value over the horizontal wavenumbers k. We see that significant contributions to the spectrum are always located inside the region delimited by the dispersion relation Eq. (5.2). Note that the interface variations are measured on a line in the x direction and not in the horizontal (x, y) plane, so any wave numbers |k| < 2ω 2 /∆b may correspond to an interfacial wave, according to Eq. (5.2). This is why the spatio-temporal spectra of Fig. 6 are not merely peaked around the dispersion relation.

Prediction of the interface elevation amplitude for a monochromatic wave

For a given interfacial wave with wavenumber modulus K, the interface elevation amplitude is denoted a K . We assume that the flow around the interface is entirely due to the potential flow associated with the interface deformation. This allows us to estimate the velocity field close to the interface as

U K ∼ a K ω .
(5.3)

Given that the flow in the mixed layer is strongly turbulent, this hypothesis may be questioned, and we will provide further discussion on this point at the end of this subsection. When the interface is infinitely sharp (∆h = 0), the horizontal velocity field due to the variation of the interface elevation is discontinuous across the interface, with a velocity jump given by ∆U K ∼ a K ω. Let us now consider that the interface is characterised by a small but non-zero thickness (∆h = 0) Using the estimate of the velocity jump obtained in the limit of an infinitely sharp interface, the vertical gradient of the horizontal velocity field is estimated as

∂ z U ∼ ∆U K ∆h ∼ a K ω ∆h .
(5.4)

Defining the local (or gradient) Richardson number inside the interface as

Ri loc = N 2 (∂ z U ) 2 ∼ ∆b∆h a 2 K ω 2 ∼ ∆h a 2 K K , (5.5)
a sufficient condition for stability of the flow inside the thin interface is Ri loc > 1/4 [START_REF] Miles | On the stability of heterogeneous shear flows[END_REF]. The key idea of [START_REF] Phillips | The dynamic of upper ocean[END_REF] is then to assume i/ that the flow is actually unstable whenever Ri loc < 1/4; ii/ that this instability eventually leads to wave breaking, which limits the growth of the wave amplitude a K ; iii/ that this is the dominant mechanism to extract energy from the interfacial wave; iv/ that an external mechanism constantly supplies energy to the wave. Then a steady state can be reached, and the saturated wave amplitude a K is such that the condition of criticality Ri loc = 1/4 is satisfied. Injecting Eq. (5.2-5.4) in Eq. (5.5), this condition for criticality yields

a 2 K ∼ ∆h K .
(5.6)

Spatial power spectrum of the interface elevation

Let us now assume that the interface is an (isotropic) collection of waves with wavenumber (k, l) (and wavenumber modulus K = √ k 2 + l 2 ), and that these waves do not interact. Let us write Ψ(k, l) the spatial power spectrum of the interface elevation. The variance of interface elevation at wave number K is related to the power spectrum through

a 2 K = √ k ′2 +l ′2 >K dk ′ dl ′ Ψ(k ′ , l ′ ) ∼ K 2 Ψ(k, l).
(5.7)

Injecting then Eq. (5.6) in Eq. (5.7) yields

Ψ(k, l) ∼ ∆h K 3 .
(5.8)

The experimental spatial power spectra are obtained by measuring the interface elevation along a line in the x direction. Let us call Ψ x (k) the power spectrum of the interface elevation along this direction. It is given by

Ψ x (k) ≡ +∞ -∞ dl Ψ(k, l) ∼ ∆h k 2 .
(5.9)

The experimental observation of the interface elevation power spectrum Ψ x (k)/∆h is shown Fig. 7-a. The spectral amplitudes have been normalised by the interface thickness predicted Eq. (2.22) with statistical mechanics arguments, i.e. by ∆h ∼ e c /∆b. We see that the Phillips prediction of a -2 slope for this spectrum is consistent with the behaviour of the experimental spectrum at low wave numbers in Fig. 7-a. However, some care must be taken to interpret these spatial spectra for large wavenumbers. Indeed, the meaning of the small scale spatial fluctuations of the interface h(x, t) defined with the sorting algorithm is not clear. These small scales may be dominated by the presence of turbulent fluctuations in the density field rather than by interfacial waves. For these reasons, we will consider in the remaining the interface h interp (x, t), obtained for each time t by fitting the interface elevation h(x, t) with a third order polynomial. Finally, rescaling the spatial power spectrum shown Fig. 7-a by the height ∆h allows us to obtain a reasonable collapse of the three different experiments, consistently with Eq. (5.9).

Temporal power spectrum of the interface elevation

Let us now introduce Φ(ω) the temporal power spectrum of the interface elevation. The variance of interface elevation at frequency ω is

a 2 K = +∞ ω dω ′ Φ(ω ′ ) ∼ ωΦ(ω).
(5.10) Injecting Eq. (5.10) in Eq. (5.6), and using the dispersion relation Eq. ( 5.2) yields

Φ(ω) ∼ ∆b∆h ω 3 .
(5.11)

Considering the statistical mechanics prediction in Eq. (2.22) for the interface thickness ∆h, one gets (5.12) This means that the amplitude of the interface displacement frequency spectrum is independent from the density jump ∆b. It only depends on the local kinetic energy. To our knowledge, this simple but important consequence of an interface thickness scaling as the inverse of the Richardson number (Ri ∼ L t ∆b/e c or Ri H ∼ H∆b/e c ) has not been discussed previously.

Φ(ω) ∼ e c ω 3 .
Experimental observations of the interface elevation temporal power spectrum are shown Fig. 7-b. The agreement between the theoretical predictions and experimental results is good. For frequency higher than the gravest linear mode of interfacial waves (represented as vertical dashed lines), the spectrum slope is consistent with the -3 slope predicted by Phillips theory. Note that the maximum observed frequencies are always larger than the buoyancy frequency N = ∆b/∆h. Perhaps more strikingly, no rescaling have been used to plot the spectra, and yet they all collapse on the same curve in the regime where the -3 slope is observed. Since the interface elevation is roughly similar for all the experiment , the kinetic energy e c around the interface within the mixed layer is not expected to vary significantly from one experiment to another. The collapse of all the experiments on the same curve confirms therefore the prediction of Eq. (5.12). This equation was obtained under the assumption that the interface thickness varies as the inverse of the Richardson number, consistently with statistical mechanics predictions. The collapse of the spectra on a single curve at high frequency is therefore an indirect test of these statistical mechanics prediction for the interface thickness of the mean vertical density profile.

Previous experimental observations of a -3 slope of the temporal power spectrum of the interface elevation were provided by [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF]. This slope is also consistent with the experiments by [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF]. However, [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF] found that the amplitude of the frequency power spectra of the interface elevation scaled as Ri -1 . Their scaling amounted to an interface thickness decreasing as Ri -2 . By contrast, [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF] found an interface thickness that was not varying significantly with the Richardson number. Our result is intermediate between both cases.

Finally, we observed in the experiments at high Richardson numbers the presence of a well identified peak in the frequency spectrum. The frequency of this peak were slightly smaller than the gravest linear interfacial mode represented as vertical dashed lines on Fig. 7-b. These sloshing frequencies may probably be attributed to nonlinear interactions between interfacial waves, and may also be associated with the presence of solitons. To our knowledge, there were no previous observation of such sloshing dynamics in similar experiments. Understanding this phenomenon will require further work.

Consistency of the approach

A key assumption of Phillips theory is that close to the interface, the velocity field U K at scale K is entirely due to the potential flow created by the variations of the interface elevation. In the experiment, the energy source for the waves is the turbulent velocity field. Let us call U turb,K the rms turbulent velocity at scale K, close to the interface, within the turbulent mixed layer. This velocity contains the contribution of all wave numbers lager than K. Phillips approach is consistent at scales such that U K ≫ U turb,K .

In order to obtain a simple estimate for U turb,K , we assume that the turbulent velocity field is unaffected by stratification, and well described by the phenomenology of three dimensional homogeneous isotropic turbulence. Of course, this hypothesis is too simplistic, but it can be used as a lowest order estimate. A much more detailed analysis and discussion on the coupling between turbulent and stratification is provided in Fernando & Hunt (1997). For the sake of simplicity, we assume in the following that there is no turbulent motion at scales larger the integral length scale L t †. According to previous notations, the velocity field at the integral length scale is U turb,K ∼ e 1/2 c . At scales K ≫ 2π/L T , assuming an inertial range, the turbulent velocity is obtained by dimensional analysis: U turb,K ∼ ǫ 1/3 K -1/3 , with ǫ ∼ e 3/2 c /L t the energy dissipation rate. Following [START_REF] Carruthers | Velocity fluctuations near an interface between a turbulent region and a stably stratified layer[END_REF], we assume that temporal fluctuations of the velocity field at a given point are given by the random advection of turbulent eddies by the integral scale eddies. This yields ω ∼ e 1/2 c K, and U turb,K ∼ ǫ 1/3 e 1/6 c ω -1/3 . Injecting Eq. (5.12-5.10-5.2) in Eq. ( 5.3), we obtain the estimate U K ∼ e 1/2 c for the velocity close to interface associated with the variation of the interface elevation. In other words, whatever the scale K such that 2π/∆h ≪ K ≤ 2π/H, the potential flow created by the interface elevation variations is of the order of the rms turbulent velocity e 1/2 c . We stress that this result relies on the assumption h ∼ e c /∆b, which was done based on the statistical mechanics result Eq. (2.22). We see that the condition for consistency U turb,K ≪ U K is valid for sufficiently small scales (K ≫ 2π/L t ), or for sufficiently high frequency (ω ≫ 2πe

1/2 c /L t ). The cut-off frequency f c = e 1/2
c /L t is of the order of 0.1 s -1 in the experiment. The maximum observed frequency is f max = 1 s -1 on Fig. 7-b. The criterion f ≫ f c is therefore only marginally satisfied in the range of the observed -3 slope.

The distribution of density levels

Now that we have characterised the properties of the interface elevation, we focus in this subsections on the statistical properties of density fluctuations. The temporal evolution of the distribution of density levels ρ b (z, σ, t) is shown Fig. 8. Each plotted distribution is obtained by building a normalised histogram of density levels for each depth z, using a sequence of 100 images separated by one second, and each successive plots of Fig. 8 are separated by 100 seconds. The maximal density levels decrease slowly with time in the lower layer, with a concomitant increase of the interface height. An experiment with weak stratification is shown on Fig. 8-d. In that case the flow is rapidly fully homogenised, the interface is not well defined, and the density field does not evolves through a sequence of stationary states.

Around the density interface, the density distributions of Fig. 8 is closed to a double peaked function. However, the fluctuations of the interface elevation have not be removed to obtain these statistics. The doubled peaked function for the density distribution around the interface is therefore mostly due the rapid motion of interfacial waves around a slowly evolving mean interface height. The sloshing dynamics of the interface only affects the the density field close to the interface. Far from the interface, the motion of the density field is not affected by the variations of the interface elevation and the method used to build the histograms of Fig. 8 is relevant to describe density fluctuations.

We see that separating the part of the density statistics due to the wave motion of the interface from the actual turbulent fluctuations is difficult in practice, and we will propose in the following a rudimentary decomposition of the density fields into "waves" and "turbulence" close to the interface.

The mean vertical profile after removing the effect of large scale interfacial waves

We present here an experimental test of the statistical mechanics predictions for the mean density profile around the interface. As a starting point, we note that the global distribution in the experimental of density levels is close to a double delta function, for Richardson numbers from 10 to 150, since the ratio of the interface width with the layer depth is much smaller than one. Equilibrium statistical mechanics theory predicts in that case a tanh shape for the x-averaged density profile, see Eq. (2.22) in subsection 2.3.

In order to test the statistical mechanics prediction for the density interface, it is necessary to cancel the spurious effect of sloshing dynamics on density statistics close to the interface. We rebuild for that purpose the histograms of density levels obtained initially Fig. 8 by considering a frame of reference following the interface height for each image. The underlying assumption is that the wave motion of the interface is dominated at lowest order by spatial mode with typical length scale larger that the horizontal size of the image.

Assuming that the system is in a quasi-stationary state on time intervals of 200 seconds after t = 100 s for each experiment, the x-averaged vertical profile of density levels b exp (z) is fitted with the function

b f it (z) = ∆b f it 2 1 -tanh z -h ∆h f it .
(5.13)

According to the statistical mechanics prediction in Eq. (2.22), all experimental vertical mean density profile should collapse on the same curve if the vertical axis is redefined by z * = (zh)/∆h f it and if the density is rescaled as b * = b/∆b f it Remarkably, all the vertical profile collapse on a curve that is very close to the predicted tanh relation (black curve) on Fig. 9-a. We note that the fit is better above the interface than below the interface †. A possible reason is that the tanh-profile corresponds to an equilibrium state for a two level system, while irreversible mixing through turbulent cascade leads to the creation of new intermediate density levels close to the interface. The effect of irreversible mixing through turbulent cascade is discussed in more detailed in the next subsection.

Variation of the interface thickness with the Richardson number

According to the statistical mechanics prediction Eq. (2.22), the interface thickness and the Richardson number Ri H introduced Eq. (4.1) are related through ∆h = 2H/Ri H . We check Fig. 9-b that the observed interface width ∆h f it is inversely proportional to the Richardson number Ri H . However, the coefficient of proportionality is larger than the one predicted by the theory: we observe ∆h f it ≈ 3H/Ri H . A possible reason for this discrepancy may be attributed to our estimate of the rms kinetic energy in the Richardson number Ri H . Indeed, e c were estimated by considering the rms turbulent kinetic energy measured in a case without stratification at the height of the interface. Assuming that the factor ∆h f it /∆h ≈ 3/2 may be attributed to the modification of turbulent properties due to the stratification, and that the statistical mechanics predictions for the interface height are correct, we define an effective Richardson number Ri f it H such that ∆h f it = 2H/Ri f it H , and we use this Richardson number to estimate the effective energy e f it c = (3/2)H∆b/Ri f it H . This effective energy can then be used to define the effective † As explained in subsection 3.2, a fit with the error function would be as good as the fit with the tanh function: the two function would be indiscernible on Fig. 2.22-b.

temperature close to the interface

β f it = 3 2e f it c = 1 ∆h f it ∆b f it .
(5.14)

Since the pioneering work of [START_REF] Crapper | The structure of turbulent density interfaces[END_REF]; [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF], the variations of the interface thickness with the Richardson number has remained highly debated, partly because different experimental settings and measurement techniques have been used, partly because different definitions for the interface thickness have been considered. [START_REF] Crapper | The structure of turbulent density interfaces[END_REF]; [START_REF] Wolanski | Turbulent entrainment across stable density step structures[END_REF]; [START_REF] Fernando | On the nature of the entrainment interface of a twolayer fluid subjected to zero-mean-shear turbulence[END_REF] found that the interface thickness was independent from the Richardson number based not he turbulence length scale (Ri = ∆bL t /(2e c /3)). By contrast, [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF] distinguished a static thickness h s measured after stopping mechanical stirring from a dynamical thickness h d measured in the presence of turbulence. They observed that the static thickness h s was independent from the Richardson number, while (h d -h s )/h s ∼ Ri -1 . Using PLIF measurements, [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF] observed that the mean interface thickness was decreasing with the Richardson number, but with a different scaling (h s ∼ L t Ri -2 ) and that the interfacial wave amplitude was decreasing as h d ∼ L t Ri -1 for sufficiently large Ri. However, [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF] found h d ∼ L t using a similar method (but better spatial resolution) to determine the interface thickness. [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF][START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF] defined the interface thickness ∆h(x, t) as the height difference between two prescribed iso-density contours (the height difference between the contours α = 0.2 and α = 0.8, where α is defined Eq. (4.2)). Just as in the case of the interface height h(x, t), this method can not be applied to a density field with a strongly corrugated interface, presenting overturning events everywhere.

Here we have presented two different experimental results supporting a scaling of the interface thickness with the inverse of the Richardson number. First, the amplitude of the frequency spectra do collapse on the same curve at high frequency, which is predicted with Phillips theory and the additional assumption that the interface thickness scales as Ri -1 . Second, the fit of the interface shape obtained in a frame of reference following the x-averaged interface elevation also yields a similar scaling. This second test is not fully satisfactory: indeed, in the presence of a perfectly thin interface ∆h ≈ 0, with interfacial waves of wavelength smaller than the windows of observation and wave amplitude scaling as Ri -1 , our procedure to obtain the interface thickness would lead to a scaling ∆h ∼ Ri -1 . We expect that the scaling ∆h ∼ Ri -1 predicted by statistical mechanics is valid for large but moderate Richardson numbers (Ri ∼ 10), when the interface is permanently breaking, while it is not valid for very high Richardson numbers, when the interface is only breaking intermittently. In this case, the interface acts as a mixing barrier that prevent mixing in physical space, and in phase space.

Vertical profile of the variance of density levels

Statistical mechanics predicts not only the mean vertical density profile but also the presence of density fluctuations across the interface. These fluctuations can be related to the mean density profile and to the effective temperature of the turbulent flow through Eq. (2.18). We assume that the inverse temperature is given by β f it = ∆h f it ∆b f it -1 defined Eq. (5.14). We also assume that the mean density profile b(z) is well described by the tanh profile defined Eq. (5.13). Then Eq. (2.18) yields

b 2 -b 2 (∆b f it ) 2 = 1 2 cosh z -h ∆h f it -2
.

(5.15)

The observed vertical variation of the variance of density fluctuations is plotted in Fig. 9-c for different experiments (corresponding to different Richardson numbers). The vertical coordinate z * = (zh)/∆h f it has been rescaled by the interface thickness for each experiment, and we consider the frame of reference following the interface elevation, just as in Fig. 9-a. The variance of density fluctuations is rescaled by ∆b f it 2 . The thin black line is the statistical mechanics prediction given by Eq. (5.15), This theoretical prediction is qualitatively correct sufficiently close to the interface. Far above the interface, the statistical mechanics theory overestimates the fluctuations : these fluctuations are absent in the upper layer since the turbulence is located in the lower layer and around the interface. Far below the interface, within the mixed layer, equilibrium theory underestimate the density fluctuations: filaments of light fluids entrained in the mixed layer are not stirred as much as would be predicted by the equilibrium theory: the density field is strongly out of equilibrium in this region.

Although the observed vertical profiles of density fluctuations close to the interface are close to the predicted vertical profile, their amplitude remains smaller than the amplitude predicted by the equilibrium theory. The main reason for this discrepancy is that the theory does not take into account irreversible mixing of density, which tends to decrease the density fluctuations. In addition, we clearly see on Fig. ( 8) that the density distributions close to the interface contains a continuum of density levels between the extremal values. The presence of these intermediate density levels which were not initially present in the two layer density field are also evidence for irreversible mixing by turbulent cascade. This effect can not be captured in the framework of the simple two level system. Modelling the combined effect of irreversible mixing (through turbulence cascade) and of the relaxation toward equilibrium has partly been addressed by [START_REF] Venaille | Modeling mixing in two-dimensional turbulence and stratified fluids[END_REF] and will be the object of future work To conclude, one can distinguish two regions for the density field: i/ The region close to the density interface, which may be interpreted as an equilibrium state once the effect of interfacial waves and interface increase due to entrainment are removed. The theory predict correctly the mean vertical density profile, but overestimate the fluctuations, and the observed density distribution is different from the initially postulated two level distribution. ii/ The region far from the interface (z ≫ ∆h), which is strongly out of equilibrium since the observed variance of density fluctuations is much larger than the one predicted by the equilibrium theory. The aim of the next section is to describe in more details the properties of density fluctuations in this out-of equilibrium region.

Exponential shape of the tracer distribution far from the interface

Sufficiently far from the interface, the fluid motion is not influenced by the motion induced by the interfacial waves, and the distribution of density levels observed in Fig. 8 (obtained without change of reference frame) can directly be interpreted as turbulent fluctuations.

We see on Fig. 10-a the distribution of density levels plotted for each experiment at five different depth corresponding to five prescribed values of the density variance (relative to the density variance at the interface for each experiment). The density levels on the horizontal axis are normalised by ∆b f it . According to Eq. (5.13), ∆b f it is the maximal value of the x-averaged vertical profile of density for each experiment. Since the initial density jump is ∆b 0 > ∆b f it , density levels larger than ∆b f it can be observed.

The depth and normalised density variance corresponding to each plotted density distribution are shown Fig. 10-b. We see on Fig. ( 8) that a given value of the density variance is associated with a well defined depth when sufficiently close to the interface, but that the depth associated with a given value of the variance are more scattered far from the interface.

Strikingly, the density distributions of Fig. 10-a do collapse qualitatively well, and are characterised by exponential tails, with an e-folding depth that decreases at increasing distances from the interface.

Exponential tails in the distribution of a tracer in turbulent flow have been previously reported either in the case of an isolated source discharging the tracer into an infinite (unconfined) medium [START_REF] Duplat | A nonsequential turbulent mixing process[END_REF], or in the case a confined medium with large scale inhomogeneities due the injection of the tracer at the boundaries. This is for instance the case in convection experiments, where exponential tails in the temperature distribution have been reported for very high Reynolds numbers [START_REF] Castaing | Scaling of hard thermal turbulence in rayleigh-bénard convection[END_REF]. In these convection experiments, one can consider that the temperature in the bulk is statistically homogeneous. The statistically steady distributions result from a competition between turbulent cascade that tends to dissipate temperature fluctuations and the tracer fluxes at the boundaries that inject fluctuations in the bulk [START_REF] Pumir | Exponential tails and random advection[END_REF]. In this approach, the density is considered as a passive tracer As far as the distribution of density levels is concerned, the mixed layer in the present stably stratified experiment is analogous to the mixed layer in the convection experiment at very high Reynolds number. The only difference is that the source of density fluctuations in the stably stratified experiment comes only from the interface at the top of the mixed layer, while the injection of density fluctuations comes from both the upper and the lower layer in convection experiments. This asymmetry in the injection of density fluctuations explains why only one tail of the density distribution present an exponential shape in the case of the stably stratified experiment.

Conclusion

We have proposed a statistical mechanics interpretation of the formation of sharp but highly corrugated density interface between region of homogeneous density in the presence of turbulence, building upon previous work by [START_REF] Tabak | Mixing in simple models for turbulent diffusion[END_REF], which generalises the Miller-Robert-Sommeria approach for the vorticity in two-dimensional turbulence to the density in three dimensional stratified turbulence. The statistics of the density field is predicted as the most probable outcome of turbulent stirring. An effective "heat" bath is provided by the turbulent velocity field. The temperature of this "heat bath" is proportional to the turbulent eddy kinetic energy.

In the case of a system initially composed of two homogeneous layers with a stable density interface, the theory predicts a tanh-shape for the mean vertical density profile. This equilibrium density profile is interpreted as the result of a competition between turbulent transport that tends to smooth out the interface and buoyancy forces that tend to sharpen the interface by the sorting of the fluid elements by density. For large Richardson numbers, buoyancy takes over turbulent transport and the interface is thin. More precisely, the interface thickness is inversely proportional to the Richardson number.

The equilibrium theory alone cannot describe entrainment across the interface, which would eventually lead to complete homogenisation of the density field. In the case of the experiments presented in this paper, turbulence generation is limited to the lower layer. The interface is progressively drifting upward by the entrainment of fresh fluid, and the lower layer density decreases due to this mixing process. Entrainment is related to irreversible mixing of density levels through turbulent cascade, a process which changes the global distribution of density levels. We propose a simplified model which keeps a two level system at statistical equilibrium, accounting for the effect of entrainment by a progressive decrease of the lower layer density.

The advantage of the statistical mechanics approach over previous models is that it provides a prediction for finite interface thickness in the presence of a source of turbu-lence without assumption about turbulent diffusivity. In addition, it does not rely on a particular mechanism at stake close to the interface (wave breaking, Kelvin-Helmholtz instability,...) that may depend on the Richardson number [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF]. The only assumption is that the system sufficiently explores the phase space.

In order to test the equilibrium theory, we assumed the existence of a statistically stationary state, and obtained statistics of the density field close to the interface using a reference frame following the horizontally averaged interface height. This method allowed us to get rid of the out-of equilibrium effects of the entrainment and of the gravest spatial modes of the interfacial gravity waves (which were found to dominate the temporal fluctuations of the interface elevation). We found that the shape of the observed mean vertical density profile is well fitted by a hyperbolic tangent function, as predicted by the equilibrium theory in the case of a two level system. Furthermore the interface thickness was found to be inversely proportional to the Richardson number between Ri ∼ 10 and Ri ∼ 100, as expected form the theory. We do not know whether these predictions would remain valid for higher Richardson numbers (still for a turbulent flow). Indeed, the interface then becomes locally sharper than the resolution of our measurements, with spatially and temporally intermittent mixing events across this interface, see also [START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF]. The statistical theory also provides a good prediction for the vertical profile of the rms density fluctuation variance at the interface. The value of this variance is however smaller than predicted by theory. The discrepancy is probably due to the dissipation of density fluctuations by the turbulent cascade.

Below the interface, the measured density fluctuations are by contrast much stronger than the statistical equilibrium prediction, which is very low in relation with the quasiuniform mean density, according to Eq. (2.17). These density fluctuations are strongly out-of equilibrium: they are transported from the interface by turbulent transport processes. In this region, we observed exponential tails in the probability distribution function of density levels. Such exponential tails can be attributed to isolated filaments that are entrained from the interface and then stirred through turbulent cascade in the lower layer. The density then behaves as a passive scalar, and the shape of its pdf results from a competition between turbulent cascade and the injection of light density filaments entrained from the interface. This situation is analogous to convection experiment at high Reynolds numbers, where exponential tails in the density pdf have been also previously reported.

Density fluctuations discussed above correspond to rather small scales, for which turbulent motion prevails. At larger scale the interface fluctuates as internal waves, for which we were able to check the dispersion relation. We observe that the spatial and temporal spectra of these waves can be interpreted by a theory due to [START_REF] Phillips | The dynamic of upper ocean[END_REF] which states that the wave amplitude at each frequency is such that the induced flow is close to criticality for shear instability at the interface. In the range of frequency much larger than the gravest mode, and much smaller than the buoyancy frequency characterising the density gradient inside the interface, this theory predicts a -3 slope for the interface elevation frequency power spectrum. We confirm this result, which also agrees with previous observations by [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF]. The value of the predicted spectral energy depends on the interface thickness as an input. We observed that the amplitude of the waves in this range of frequency is independent from the Richardson number, which corresponds to the statistical equilibrium prediction of an interface thickness varying as the inverse of the Richardson number. This scaling is different than the one obtained by [START_REF] Hannoun | Turbulent mixing at a shear-free density interface[END_REF], and we do not know whether it would still be valid for higher Richardson number (beyond 100), as explained above. Note finally that the Phillipp's approach is limited to frequencies for which the turbulent velocities are negligible with respect to the potential flow associated with the variations of the interface elevation. A more detailed study of the coupling between interfacial waves and turbulence has been addressed by Fernando & Hunt (1997);[START_REF] Mcgrath | Turbulence, waves and mixing at shear-free density interfaces. Part 2. Laboratory experiments[END_REF]. We focused in this paper to the visualisation of the density field using PLIF technique with index matching. The only PIV measurement presented in this paper were performed in a case without stratification. The underlying assumption was that turbulence properties were not much affected by the presence of stratification (below the interface). We chose to use PLIF only in order to obtain quantitative measurements of turbulent fluctuations in the density field (the presence of particle used the PIV measurements alters the quality of the visualisation). However, we hope to address in future for a more detailed study of the buoyancy fluxes close to the interface using Simultaneous PIV and PLIF measurement, which has recently been performed in the context of entrainment in gravity currents [START_REF] Odier | Entrainment and mixing in a laboratory model of oceanic overflow[END_REF].

To conclude, equilibrium statistical mechanics allows us to interpret qualitatively the formation of thin but corrugated and turbulent density interface between regions of homogeneous density. The approach is limited by out-of equilibrium effect such as irreversible mixing and interfacial wave excitation. We believe that combining these different phenomenological approaches will lead to fruitful models in the context of turbulent mixing across a density interface, or more generally in any stratified turbulence problem where density fluctuations play an important role, as for intense in gravity currents [START_REF] Odier | Fluid mixing in stratified gravity currents: The prandtl mixing length[END_REF][START_REF] Odier | Understanding and modeling turbulent fluxes and entrainment in a gravity current[END_REF]. It would be also interesting to check the relevance of equilibrium statistical mechanics at the interface of two immiscible fluids. On the one hand, this would avoid the issue with irreversible mixing through turbulence cascade, but on the other hand, surface tension effects may be influential.

7. Appendix A: Equipartition and mixing efficiency of 0.25 in the low energy limit

Eq. (2.18) allows to derive an interesting side result concerning energy equipartition in the low energy limit. This result may in turn be used to predict a mixing efficiency coefficient, which is a measure of the fraction of the energy injected in a stratified fluid that is actually used to irreversibly increase the potential energy of the flow.

Let us consider a stratified fluid initially at rest, characterised by its density profile b bg (z) (the index "bg" stands for "background"). The potential energy of this state if denoted E bg,p ≡ H 0 dz bz. Let us then assume that the fluid is isolated, that a given amount of energy ∆E is injected in the system, and that a statistical equilibrium is reached on a time scale much shorter than the time scale for viscous dissipation and for irreversible mixing of density. Let us call E c and E p the kinetic and the potential energy of the equilibrium state. The available potential energy of the equilibrium state is The second equality is obtained by noting that the equilibrium state is statistically invariant on the horizontal, and that the horizontal integral amounts to an ensemble average. For the equilibrium state, the potential energy ∆E p may be qualified as available since one recover b = b bg just by setting the effective temperature to 0 (β = +∞). Since each fluid particle conserves its density in the absence of dissipation, a fluid particle of density b at height z and time t can be seen as a fluid particle initially at height z bg (b), displaced by ξ(x, y, z) = z -z bg (b(x, y, z)) in the vertical direction. Let us assume that this displacement ξ is sufficiently small, which is ensured by considering a low energy limit. At lowest order, the density fluctuation defined by b ′ (x, y, z) = b(x, y, z) -b bg (z) is proportional to the vertical fluid particle displacement: b ′ = ξ∂ z b bg , and the available potential energy can be expressed (still at lowest order) as

∆E p = 1 2 dxdydz b ′2 ∂ z b bg = 1 2 dz b 2 -b 2 ∂ z b bg . (7.2)
Again, the second equality is obtained by noting the horizontal integration amounts to an ensemble average. Injecting then Eq. (2.18) in Eq. ( 7.2) gives

∆E p = 1 2β dz ∂ z b ∂ z b bg . (7.
3)

In the low energy limit, we get at lowest order ∂ z b ≈ ∂ z b bg , which yields ∆E p = V /(2β), where V is the volume where the flow takes place. The inverse temperature β is related to the total kinetic energy ∆E c = V e c through Eq. (2.15), which yields ∆E p = ∆E c /3. This expresses equipartition of the energy between the available potential energy and the three degrees of freedom of the kinetic energy.

Finally, the ratio of the available potential energy with the total energy injected in the system is

η ≡ ∆E p ∆E p + ∆E c = 1 4 . (7.4)
Let us now assume that once the equilibrium state is reached, the density fluctuations are smoothed out on each horizontal plane due to the combined effect of direct turbulent cascade and molecular diffusivity. Let us also assume that the rate of kinetic energy dissipation is equal to the rate of dissipation for the density variance at each height. These hypothesis ensure that the rhs of Eq. (2.18) remains constant through the flow evolution, and that the profile b remains the equilibrium state throughout the flow evolution. At sufficiently large time, once the fluctuations around b are irreversibly mixed, the flow is at rest and this b becomes the new background density profile. The increase of potential energy ∆E p defined by Eq. (7.1) accounts therefore for the irreversible increase of potential energy †. The mixing efficiency defined as the irreversible increase of potential energy normalised by the total energy injected in the system is then simply given by Eq. (7.4).

To conclude, a mixing efficiency coefficient of 1/4 can be interpreted as a result of energy equipartition at equilibrium under the following assumption: i/ a low energy limit ii/ a time scale for the dissipation of horizontal density fluctuations much larger than the time scale to reach the equilibrium state iii/ a rate of dissipation of the variance of density fluctuation at each height equal to the rate of dissipation of the kinetic energy. If these conditions are not fulfilled, then the mixing efficiency should be smaller than 0.25.

Mixing efficiency coefficients between 0.2 and 0.3 are widely used in modelling context. However, experiments and simulations and observations seem to show that there is no universal mixing efficiency in stably stratified turbulence [START_REF] Peltier | Mixing efficiency in stratified shear flows[END_REF]; [START_REF] Ivey | Density stratification, turbulence, but how much mixing?[END_REF]. Our result suggest that the value 0.25 may be interpreted as a limit case in the framework of equilibrium statistical mechanics.

Appendix B : Turbulence properties in the homogeneous case

We performed one experiment without stratification in order to have a reference flow that we analysed using PIV measurements. This allowed us to obtain some characteristics of the flow in the turbulent layer only, assuming that these flow properties would not be much different in the presence of stratification provided that one consider the dynamics sufficiently below the interface. In particular, this methods allows to estimate the energy flux due to the grid forcing, and the decay of the turbulence strength with altitude.

Let us call L t (z) the integral length scale of turbulence. Sufficiently close to the source, the integral length scale of turbulence if given by the grid mesh of the oscillating grid L t = L f . Sufficiently far from the source, the only length scale of the problem is the distance from the source and L t ∼ z. The temporal evolution of the kinetic energy may be modelled as (8.1) where the effect of turbulence is modelled as effective viscosity, with a turbulent energy flux

∂ t e c = ∂ z (ν t ∂ z e c ) - c d L t e 3/2 c ,
F ec = -ν t ∂ z e c , ν t = ae 1/2 c L t . (8.2)
The sink of energy is given by a Kolmogorov dissipation term. Sufficiently close to the grid, L t = L f is a constant and stationary kinetic energy profile is exponential:

e c = e 0 c exp - z L t , L t = a c d 1/2 L (8.3)
with e 0 c the kinetic energy at the source located in z = 0. Sufficiently far from the source L t ∼ z and stationary kinetic energy profile is given by a power law e c ∼ z -2 , see e.g. [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF].

  Figure 1. a) Sequence of statistical equilibrium states in the case of a two level global density distribution (b ∈ {-Σ(t), Σ(t)}), at fixed temperature β (prescribed by a turbulent heat bath).The value of the density level Σ(t) is decreasing from 1 (plain green line) to 0 (dashed red line). b) purely diffusive relaxation of an initial step function towards an homogeneous density profile. The diffusion coefficient is homogeneous in space.

Figure 2 .

 2 Figure 2. a) Experimental setting. The tank is filled below with a layer of water with salt and rhodamine, and above with water, ethanol and rhodamine. The LASER sheet illuminates the centre of the tank. The density field is observed in the 25 × 25 cm 2 window represented with a black line. b) Vertical variation of the turbulent kinetic energy in the 25 × 25 cm 2 window when there is no stratification (PIV measurements).

Figure 3 .

 3 Figure 3. a) Snapshot of the density field at the centre of the tank, normalised between 0 (light fluid, blue color) and 1 (dense fluid, red color) . b) Same density field, but for each value of x the vertical density profile is sorted with denser fluid below, see subsection 4.3. The black line is the interface defined as the height of the intermediate density level (the contour α = 1/2) in the sorted field. The white line is a fit of this interface height with an order 2 polynomial.

  Fig 4 for three different experiments associated with decreasing Richardson numbers from panel a to c. It always takes around 20 seconds before the turbulence reaches the interface. Then the averaged density of the lower layer decreases while the interface height increases slowly. In the three experiments presented

Figure 4 .

 4 Figure 4. Temporal evolution of the density fields averaged in the horizontal direction for 3 different experiments, see Tab. 1: a) ∆b = 0.13 m.s -2 (EXP4) b) ∆b = 0.033 m.s -2 (EXP1) c) ∆b = 0.011 m.s -2 (EXP7). The black line represents the interface (see text). The density is normalised such that it varies between 0 and 1 at the initial time t = 0 for each experiment.

Figure 5 .

 5 Figure 5. a) mass conservation and validation of PLIF measurements. The thin plain line represents the temporal evolution of the x-averaged interface height h(t) normalised by the initial height h0, and obtained with the sorting algorithm. The "+" symbols represent the temporal evolution of the averaged density < b > in the mixed layer, normalised by the initial density < b0 > in this layer. b) Entrainment coefficient as a function of the average Richardson number for each experiment. The entrainment coefficient is defined as the ratio between the interface velocity and the turbulent velocity. The interface velocity is obtained by a linear fit of the curve of the left panel over the interval t = 100 -300 s. The turbulent velocity around the interface height is estimated by using the PIV measurements performed in the case without stratification.

Figure 6 .

 6 Figure 6. a) Interface height spectrum for EXP1 (∆b = 0.03 m.s -2 ). For each frequency f , the spectrum is normalised by its maximal value. The black line is the dispersion relation for interfacial waves given Eq. (5.2). b) Same plot for EXP4 (∆b = 0.13 m.s -2 ).

Figure 7 .

 7 Figure 7. a) Spatial power spectrum. Blue diamonds: ∆b = 0.011 m.s -2 (EXP7); Green circles: ∆b = 0.03 m.s -2 (EXP1); Red squares: ∆b = 0.13 m.s -2 (EXP4); b) Temporal power spectra of the x-averaged interface height time series. The vertical dashed lines represent the location of the gravest linear interfacial mode for each experiment, given by ω0 = π∆b/L, where L is the lateral extension of the tank.)

Figure 8 .

 8 Figure 8. Temporal evolution of the distribution of density levels: a) EXP4; b) EXP1; c) EXP7; d) EXP6. For each experiment the density levels are normalised from 0 (light fluid, blue color) to 1 (dense fluid, red color), and a logarithmic scale has been chosen to visualise fluctuations far from the interface.

Figure 9 .

 9 Figure 9. a) Vertical profile of the x-averaged density. The thin black line is the tanh profile predicted by equilibrium statistical mechanics in the two level case. Red circle: EXP7; black triangle: EXP1; magenta triangle: EXP2; green star: EXP3; blue diamond: EXP4. b) Variation of the interface thickness ∆h f it with the Richardson number. c) Vertical variations of the variance of the density distribution.The thin black line in the statistical mechanics prediction Eq. (5.15).

Figure 10

 10 Figure 10. a) Probability distribution function of the density at increasing distances from the interface. Each symbol corresponds to a different experiment (o: EXP7, ⋄: EXP2, * : EXP4). For each experiment, 5 different depth (associated with 5 different colors on the plot) are considered. Those 5 depth are chosen such that the density variances of the pdf normalised by the density variance at the interface are the same for a given color (the relation between depth and the prescribed values of density variance is shown on panel b). For each experiment, the density levels on the x axis are normalised by the maximum value of the x-averaged density, denoted ∆b f it . b) Variation of the height (vertical axis) associated with prescribed values of relative density variance ( b 2 -b 2 / b 2 -b 2

  ∆E p ≡ dxdydz (b -b bg ) z = dz b -b bg . (7.1)

Table 1 .

 1 Parameters

	Experiment EXP1 EXP2 EXP3 EXP4 EXP5 EXP6 EXP7
	∆ρ/ρ (%)	0.33	0.3	0.8	1.3	0.1	0.4	0.11
	Ri	28	28	84	144	0.9	3.5	12
	RiH	127	125	370	640	3.5	14	55

for the different experiments. The density jump ∆ρ/ρ is estimated as the beginning of each experiment. See Eq. (4.1) for the definition of Ri and RiH.

† This assumption can not be fully valid in our experiment, since the presence of large scale flow structure filling the whole domain are often reported in confined turbulent flows).

† Note that Eq. (7.2) is derived from Eq. (7.1) by assuming that each fluid particle conserves its density; these equations are no more equivalent once the fluctuations of density have been smoothed out at each height z (in which case the available potential energy vanish). However, the result ∆Ep = ∆Ec/3 remains valid since it is obtained at equilibrium.