Parking selection in vehicle routing

Théo Le Colleter, Fabien Lehuédé, Olivier Péton IMT Atlantique, LS2N, UMR CNRS 6004, Nantes, France Dorian Dumez, HEC Montréal, CIRRELT, Montréal, Canada Email: {fabien.lehuede;olivier.peton}@imt-atlantique.fr, dorian.dumez@hec.ca This paper presents a variant of the vehicle routing problem considering the delivery of customers in cities with a combination of walking and driving. The objective is first to offer a better model for actual delivery problems where drivers are usually assumed to have a single parking location to deliver each customer. Second, we show that serving some customers through walking trips is more efficient in congested areas. We introduce the Park-and-Loop Routing Problem with Parking Selection (PLRP-PS), which extends the work of [START_REF] Coindreau | Vehicle routing with transportable resources: Using carpooling and walking for on-site services[END_REF][START_REF] Cabrera | The Doubly Open Park-and-loop Routing Problem[END_REF]. We investigate the case where parking locations should be selected among a large set of parking areas in a variant of Large Neighborhood Search (LNS).

The park-and-loop routing problem with parking selection

We consider a set of customers N with non-unitary demand q i , ∀i ∈ N who have to be served using a homogeneous vehicle fleet K. We consider that each vehicle starts from and return to a depot 0. It is assigned to a driver who will perform deliveries by walking from a parking location to the customers. We consider a large set P of parking locations. Each parking location i ∈ P is associated with a parking time pt i . For each pair of nodes i, j in {0} ∪ P , we assume a driving time t d ij and distance d d ij are determined, while a walking time t w ij and distance d w ij are defined for each i, j ∈ P ∪ N . A vehicle route starts from the depot, visit a sequence of parking locations in P , and returns to the depot. The sum of customers demands served by a vehicle route must not exceed the vehicle capacity Q d and a vehicle route duration must not exceed a maximum value h max . A walking-trip starts from a parking location in P , visits a sequence of customers, and returns to the parking location. It must not serve more than the walking capacity denoted Q w . The overall walking distance of each driver should not exceed a given limit w max . At a parking location, a walking-trip can start after the vehicle is parked.

Then the driver may perform several walking-trips from his vehicle.

The PLRP-PS consists of designing vehicle routes and walking-trips, respecting vehicle capacities, maximum route duration and the maximum walking distance of each driver, such that each customer is served with its demand by exactly one driver, lexicographically minimizing the number of vehicles and the sum of driving times, walking times and parking times. This problem is illustrated on Figure 1

Solution approach

To solve the PLRP-PS, the proposed algorithm is based on a variant of LNS called Small and Large Neighborhood Search (SLNS). We investigate different strategies that can be used in the SLNS repair operators to select parking locations.

Small and Large Neighborhood Search

The principle of SLNS where first introduced in [START_REF] Dumez | Hybridizing Large Neighborhood Search and Exact Methods for Generalized Vehicle Routing Problems with Time Windows[END_REF]. It combines ingredients from the seminal paper of [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF] and the recent paper of [START_REF] Christiaens | Slack induction by string removals for vehicle routing problems[END_REF]. In this algorithm, most iterations perform small destroy and repair operations on the current solution, using in particular, fast list heuristics to re-insert customers. On a regular basis, large destructions are performed, involving additional operators to bring diversification in the search.

Classical destroy and repair operators where tailored to the problem and extend to integrate parking selection. To create the vehicle routes and walking trips, we introduce three types of insertions, namely insertion in a new route, insertion through a new parking location and insertion through a served parking location, which allows customers to be inserted in existing walking-trips.

Insertion strategies for parking selection

Considering a huge number of parking locations may be available in a city, independently of the number of customers, we are interested in reducing the possibilities that are evaluated when selecting a parking location. Let us presume that a customer has been selected and we aim at evaluating its insertion in a new route or through a new parking location. We propose different parking selection strategies. A parking selection strategy is composed of a parking selection criterion and parking selection method.

We consider two criteria to select parking locations from each customer: The first criterion is the walking time between the parking and the customer (radius-r criterion).

The second criterion is the ranking of the parking in the list of parking ordered by nondecreasing distance to the considered customer (k-nearest criterion). In addition, two approaches involve these criteria to reduce the set of evaluated parkings: The first one filters the parking set from which each customer can be served a priori, keeping the locations below a given threshold on the chosen criterion. In the second approach, for each customer a sorted list of parking locations is designed based on the selected criterion.

The locations are then evaluated according to this order, with a given probability to stop after each location evaluation.

Experiments and case study

The efficiency of the approach is demonstrated on the Park-and-Loop Routing Problem (PLRP), with 13 new best solutions found on an existing benchmark solved by [START_REF] Coindreau | Vehicle routing with transportable resources: Using carpooling and walking for on-site services[END_REF] and also solved by [START_REF] Cabrera | The Doubly Open Park-and-loop Routing Problem[END_REF]. Realistic instances were generated based on open data from the city of Nantes. The proposed parking selection strategies are compared and tuned on these instances. Figure 2 present experiments on problem variants, where the average driving, walking and parking times per customer is represented for different problem sizes. In this figure, the CVRP supposes parking at each customer, the PLRP allows drivers to perform walking trips while parked at customers, the PLRP-PS(C+DA) adds parking at delivery areas, while the PLRP-PS(DA) includes only delivery areas. Comparing the CVRP and the PLRP-PS(DA), we find that combining walking and driving in urban delivery can save 19% of working time on average compared to the classical vehicle routing approach.

In addition, parking at customers saves walking time when this is possible, and serving several customers from one parking is even more profitable in this case.

  , with three vehicles routes (green, grey and blue) starting from the black depot on the left. Walking paths are in red. Red points indicate customers and the dark blue points indicate the selected parking locations. Light blue points indicate unused delivery areas.

Figure 1 :

 1 Figure 1: Best found solution for instance "nantes-50-1" with 50 customers. Map created using Open Street Map ©: www.openstreetmap.org/copyright, travel time and distances with Open-RouteService ©: www.openrouteservice.org

Figure 2 :

 2 Figure 2: Objective function value decomposition for each type of time and per customer


Pascale Bouruet-Aubertot

C Koudella

C Staquet

K B Winters

Particle dispersion and mixing induced by breaking internal gravity waves

Keywords: Particle dispersion, Mixing, Internal gravity waves, Wave breaking

The purpose of this paper is to analyze diapycnal mixing induced by the breaking of an internal gravity wave -the primary wave -either standing or propagating. To achieve this aim we apply two different methods. The first method consists of a direct estimate of vertical eddy diffusion from particle dispersion while the second method relies upon potential energy budgets [Winters, K.

Introduction

In stably stratified fluids such as the ocean or the atmosphere, turbulent motions are responsible for irreversible mixing. As stratification effects limit vertical motions, vertical exchanges are strongly inhibited. As a consequence vertical mixing is weak. This is attested by values of vertical eddy diffusivity which are of a few orders of magnitude lower than horizontal diffusion coefficients, especially in the ocean interior. Nevertheless, in strongly stratified regions, the evaluation of vertical mixing is crucial for two main reasons. The first reason is to allow the prediction of exchanges within chemical or biological processes, namely in the oceanic thermocline the exchanges of nutrients, oxygen, CO 2 and other gases. The second reason is to provide parameterizations of eddy diffusivity in large scale circulation models, which do not resolve these motions. This is a central issue since in the stratified ocean, vertical mixing is one of the basic processes which sets the intensity of the general circulation [START_REF] Toole | Turbulent mixing in the ocean: intensity, causes, and consequences[END_REF]. Consequently, predictions of circulation models appear to be strongly sensitive to the parameterization of vertical eddy diffusivity [START_REF] Bryan | Parameter sensitivity of primitive equation ocean general circulation models[END_REF][START_REF] Cummins | Sensitivity of the GFDL ocean general circulation model to a parameterization of vertical diffusion[END_REF].

It is well-known that the intensity of vertical mixing strongly depends on the mechanism responsible for the onset of turbulence [START_REF] Gregg | Diapycnal mixing in the thermocline: a review[END_REF]. Therefore, process-oriented studies are useful in order to help quantify this issue. In this paper, we adopt this approach and focus on mixing induced by an internal gravity wave of high frequency (frequency of the same order as the buoyancy frequency). In the ocean, such high frequency waves are generated, for example, by means of flow interaction with small scale bathymetric variability. Our study describes the evolution of high frequency waves in a linear stratification. This study could be of relevance in the oceanic thermocline: there, shear instability is an important process which promotes mixing and it is more than likely that internal gravity waves initiate this process [START_REF] Toole | Turbulent mixing in the ocean: intensity, causes, and consequences[END_REF]. We are interested here in the evolution of a primary wave toward breaking. The dynamics of this wave, either standing or propagating, have been previously studied and are reported in [START_REF] Bouruet-Aubertot | Breaking of standing internal gravity waves through two-dimensional instabilities[END_REF][START_REF] Polzin | Statistics of the Richardson number: mixing models and finestructure[END_REF] for standing waves and in [START_REF] Koudella | Breaking progressive internal gravity waves: two-and three-dimensional numerical experiments[END_REF] and [START_REF] Koudella | Ondes internes de gravité en fluide stratifié: instabilités, turbulence et vorticité potentielle[END_REF] for propagating waves.

The major difficulty in mixing studies is to accurately estimate the mixing rate, that is the rate at which potential energy irreversibly increases through mixing. This rate is the turbulent diapycnal flux. This flux is most often compared either to the dissipation rate of energy in the flow or, assuming the mean (horizontally averaged) density gradient varies smoothly over the integral scale of turbulent events, to the mean density gradient. The former ratio defines the mixing efficiency while the latter defines the turbulent diapycnal diffusivity (or simply the eddy diffusivity). In the present paper, we shall use an exact method to compute the diapycnal flux, following [START_REF] Winters | Available potential energy and mixing in stratified fluids[END_REF] valid whatever the flow conditions. Predictions by the model of [START_REF] Osborn | Oceanic finestructure[END_REF] are recovered using this method when the flow is statistically stationary and the density gradient is approximately constant over the integral scale of turbulence. We shall address two main questions, which we present within the context of the oceanic thermocline.

The eddy diffusivities inferred from particle dispersion and from the diapycnal diffusive flux: how do they compare?

The eddy diffusivity can also be inferred directly from particle dispersion, in a way similar to the tracer experiments made in situ. We aim here at providing a careful comparison between this estimate of the eddy diffusivity and that inferred from the diapycnal flux. This is of particular interest regarding the interpretation of in situ measurements. Indeed, regardless of the turbulent processes that yield mixing, the eddy diffusivity should not depend upon the method that is used to compute it, either from microstructure measurements or by passive scalar diffusion. Recently this uncertainty has been lifted with the dispersion experiments of [START_REF] Ledwell | Evidence for slow mixing across the pycnocline from an open-ocean tracer-release experiment[END_REF] who estimated the effective diffusion of dye after a few months. The latter estimate appeared to be in good agreement with estimates from microstructure measurements. Later, careful comparison between the two methods in a smaller water reservoir has been conducted by [START_REF] Wuest | Comparison of diapycnal diffusivity measured by tracer and microstructure techniques[END_REF] in lakes. They reported an agreement between the two methods within a factor of 2.

Can we relate vertical mixing with the process of breaking of a large scale primary wave?

This question involves both the rate of mixing and the spatial distribution of mixing. In the ocean it is usually assumed that mixing efficiency is constant. This is one of the main assumption of the Osborn model which relates the eddy diffusion coefficient to the rate of dissipation of kinetic energy [START_REF] Osborn | Estimates of the local rate of vertical diffusion from dissipation measurements[END_REF]. Our study aims at determining whether the temporal intermittency of breaking has any effect on the value of the mixing efficiency. Secondly, because of the irreversible processes associated with breaking, it is expected that the background stratification will be modified in a specific way. This question deals with the relationship between the observed density structure and the flow dynamics. In the ocean a consequence of mixing is that density profiles often show a step-like structure. Our purpose is thus to study the role of internal wavebreaking with respect to the step-like structure mentioned by [START_REF] Eriksen | Measurements and models of fine structure, internal gravity waves, and wave breaking in the deep ocean[END_REF]. For this aim, we shall mainly investigate an argument proposed by [START_REF] Garrett | Oceanic mixing by breaking internal gravity waves[END_REF], which relates this step-like structure to the instability of thin layers of high vertical shear, this vorticity concentration resulting from internal waves of high vertical wavenumber.

The paper is organized as follows. In Section 2 we describe the numerical models for Eulerian and Lagrangian simulations. An overview of the wave dynamics is provided in Section 3. Section 4 concerns numerical experiments of particle dispersion. Dispersion coefficients of sets of particles released in the flow while the wave is breaking are estimated and we discuss these results in the light of previous works. Section 5 is devoted to mixing. We first introduce the analysis proposed by [START_REF] Winters | Available potential energy and mixing in stratified fluids[END_REF] which relies on the discrimination between background and available potential energies. This analysis allows a calculation of the diffusion coefficient from the diffusive flux across the isopycnals, i.e. the constant density surfaces [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF]. We then compare this evaluation with that from particle dispersion. The energetics of mixing is discussed next and we provide a parameterization of the vertical eddy diffusivity versus the turbulent Froude number of the flow. Mixing results in a weak distortion of the density profile, which we discuss in relation with oceanic measurements at small scales. Comparisons between the analysis of the twoand three-dimensional simulations of the same primary wave are eventually conducted and conclusions are drawn in Section 6.

Numerical models for Eulerian and Lagrangian simulations

Resolution of the Boussinesq equations

We solve the two-and three-dimensional Navier-Stokes equations in the Boussinesq approximation. In this approximation, the mean density, as well as the total density, is close to a constant density ρ 0 (associated with the state of hydrostatic equilibrium) so that buoyancy solely comes into play through the buoyancy force and incompressibility can be assumed. The initial density profile, ρ lin , is linearly decreasing with the vertical coordinate y so that the buoyancy frequency, N 2 lin = (-g/ρ 0 ) dρ lin /dy, is constant (g is the acceleration due to gravity). We use Cartesian coordinates (x, y) in two dimensions and (x, y, z)inthree dimensions, where y is the vertical coordinate, directed upwards. The components of the velocity field, u, are denoted as (u, v, w).

In two dimensions, the incompressibility condition is taken into account by introducing a stream function ψ, which is related to the two velocity components u and v by u=∂ y ψ, v=-∂ x ψ. We solve the dimensional Boussinesq equations in terms of the vorticityψ and the reduced density fluctuations

ρ * =(g/ρ 0 )ρ ′ ∂ t ψ + J ( ψ, ψ) = ∂ x ρ * + ν 2 ψ
(1)

∂ t ρ * + J(ρ * ,ψ) =-N 2 lin ∂ x ψ + ν Pr ρ * (2)
where the nonlinear advective terms are written by means of the Jacobian J. Note that we decompose the total buoyancy into the sum of a fluctuating part, ρ * , and a mean part, N 2 lin , and solve the equation for the fluctuating part. This is done for numerical reasons since the reduced density, ρ * ,v erifies the required boundary conditions as opposed to the total buoyancy (see below). The numerical computations are performed within a square domain in a vertical plane.

In three dimensions, we solve the full equations of motions, within the Boussinesq approximation

∂ u ∂t + u •∇ u =-∇p ′ -ρ * i y + ν ∇ 2 u (3) ∂ρ * ∂t + u •∇ρ * = N 2 lin v + ν Pr ∇ 2 ρ * (4) ∇• u = 0 (5)
where p ′ is the pressure fluctuation divided by ρ 0 and i y is the unit vector along the vertical direction orientated upward. We shall only consider propagating waves in this case, within a cubic domain. In both sets of equations, the friction is modeled as an ordinary diffusion, using a Laplacian term. ν is the kinematic viscosity and Pr=ν/κ is the Prandtl number, which varies within [0.72, 10] but most often taken equal to unity. Time will be displayed in units of the buoyancy period T=2π/N lin . Boundary conditions are periodic in all (two or three) directions for velocity and density fluctuations in the case of propagating waves. In contrast, for two-dimensional standing waves, normal velocities at the walls are set to zero as well as density fluctuations at the horizontal boundaries. Assuming that these functions can be expressed as Fourier series, these two cases are compatible with the use of a pseudo-spectral method (e.g. [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]. In the case of standing waves, symmetry properties arise as a consequence of the boundary conditions: each field is either an odd or an even function of each spatial coordinate which permits one to develop the field into either sine or cosine series. Integration in time is performed using an explicit third order Adams-Bashforth scheme. Resolution is equal to 513 2 for two-dimensional standing waves; for propagating waves, the resolution is 512 2 in two dimensions and 256 3 in three dimensions. Before each run, the viscosity ν is adjusted so that the equation for the total energy is satisfied within 1%. The value of ν depends upon the amplitude of the primary wave, for a given spatial resolution. Note that through dimensional analysis three nondimensional numbers arise which are a Froude number, which characterizes the importance of nonlinear terms with respect to buoyancy effects, a Reynolds number and a Prandtl number. Parameters of the runs as well as these three nondimensional numbers are summarized in Table 1.

The results of our calculations can be easily transposed to oceanic and atmospheric flows having the same Froude number and Reynolds number as those of our calculations. We look here at the evolution of a large scale primary wave, of wavelength equal to λ=2π in the x and y directions, and which is evolving in a background linear stratification, of buoyancy frequency N lin . The numerical domain has an area equal to λ 2 in the case of propagating waves and to λ 2 /4 in the case of standing waves. Comparisons with oceanic or atmospheric flows can be made through a similarity in Froude number, Fr, and Reynolds number, Re. Fr and Re are defined from the primary wave parameters and background stratification: Fr=U/(N lin λ)a n dRe=Uλ/ν, where U is the amplitude of the horizontal a Wavevector, k, is taken equal to (1, 1). A is the amplitude of the stream function of the primary wave, λ the wavelength of the wave, dt the time step (in units of N -1 ), ν the viscosity and Pr is the Prandtl number. Two dimensionless numbers, the Froude number (Fr) and the Reynolds number (Re), are given in the last two columns. Applications of our calculations to oceanic or atmospheric flows can be made provided that the Froude number is the same. We take the velocity of the wave and the wavelength of the wave as characteristic velocity and length for the calculation of Fr and Re. velocity of the primary wave. As an example of application of our simulations to an oceanic situation, let us consider an internal wave propagating in the thermocline, with Fr=0.04 and Re=10 700 (as for runs P1 and S1 in Table 1). These parameters may characterize a wave with wavelength λ=20 m and frequency N/ √ 2 evolving in a background stratification N=2, 5×10 -3 rad s -1 . Note that the Reynolds number of this oceanic wave is slightly higher than that of the numerical simulation (Re=3×10 4 compared to Re=10 700). However, this difference in Reynolds numbers has no effect on the dynamics which is controlled by nonlinear wave-wave interactions and is therefore prescribed by the value of the Froude number. We checked indeed that the rate of decrease of the energy of the wave does not vary with the Reynolds number [START_REF] Bouruet-Aubertot | Breaking of standing internal gravity waves through two-dimensional instabilities[END_REF]. This means in particular that the product of the enstrophy and the molecular viscosity is constant. We discuss further the implications for mixing in Section 5.2.4.

Numerical method for particle trajectories

Trajectories for weightless particles have been calculated for the case of a two-dimensional propagating wave only, by solving the equation of advection

d dt X( a, t) = V( X( a, t), t) (6)
where X( a, t) and V( X( a, t), t) denote the position and velocity of the particle being at position a at t=0.

Subroutines for particle trajectories have been provided by Babiano and are described in [START_REF] Babiano | Single particle dispersion, Lagrangian structure function and Lagrangian energy spectrum in two-dimensional incompressible turbulence[END_REF]. A Runge-Kutta method of order 2 is chosen for time implementation and spatial interpolation is performed by using cubic spline following the method introduced by De Boor (1978). Note that owing to periodic boundary conditions, the matrix of the basis of spline functions is calculated in spectral space, once and for all, at the initial time. This contributes to the efficiency of the method (further details are provided in [START_REF] Zouari | Expériences numériques lagrangiennes à partir de modèles euériens[END_REF].

Wavebreaking processes

Initial conditions

The initial condition is a simple monochromatic internal gravity wave of large scale, which is either propagating or standing in a vertical plane. This wave is referred to as the primary wave.

In two dimensions, the propagating primary wave at t=0 is prescribed as

ψ(x, y, t = 0) = A sin(k x x + k y y), ρ * (x,y,t = 0) = A k x N 2 lin ω sin(k x x + k y y) (7)
while in three dimensions, it is prescribed as u(x,y,z,t = 0) = Ak y cos(k x x + k y y), v(x, y, z, t = 0) =-Ak x cos(k x x + k y y), w(x,y,z,t = 0) = 0 (8)

ρ * (x,y,z,t = 0) = A k x N 2 lin ω sin(k x x + k y y) (9)
Thewavevector(k x , k y ) is equal to (1, 1) in two dimensions and to (k x , k y , k z )=(1, 1, 0) in three dimensions. A is the amplitude of the stream function (as well as, in the present case,

where k x =k y =1, the amplitude of the velocity field, whether two-or three-dimensional).

The frequency ω is inferred from the dispersion relationship

ω = N lin cos θ ( 10
)
where θ is the angle of the wavevector with the horizontal direction. Since the initial condition ( 7) is a solution of the nonlinear Boussinesq equations, an additional perturbation must be introduced at the initial state to destabilize the primary wave. In two dimensions, a white noise of uniform distribution is thus superimposed upon the velocity field at t 0 ,of amplitude 10 -5 A. In three dimensions, the perturbation is made of a superposition of linear internal waves with random phases, of the same amplitude as in two dimensions. In the standing case (i.e. for an internal wave in a closed domain) the exact solution of the nonlinear Boussinesq equations can only be computed as an expansion as a function of the Froude number (e.g. [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF][START_REF] Bouruet-Aubertot | Breaking of standing internal gravity waves through two-dimensional instabilities[END_REF]. Therefore, for simplicity, the initial condition we have chosen is the solution of the linearized Boussinesq equations

ψ(x, y) = A sin(k x x)sin(k y y), ρ * (x, y) =-A k x N 2 lin ω cos(k x x)sin(k y y) (11)
This standing wave is a superposition of four waves making the same angle with the horizontal. It can also be viewed as a Green-Taylor vortex in a vertical plane, having the same vorticity field. In this standing case, the actual solution will progressively depart from the linear approximation because of nonlinear interactions, so that no small amplitude white noise needs to be added.

As in all such stability problems, the time of breaking will depend on the amplitude of the noise perturbation for propagating waves and on the degree of imbalance from the nonlinear solution for standing waves.

Overview of the propagating and standing wave evolutions

It has been well-known since the two-dimensional linear stability analysis of [START_REF] Drazin | On the instability of an internal gravity wave[END_REF] and [START_REF] Mied | The occurrence of parametric instabilities in finite-amplitude internal gravity waves[END_REF] that a statically stable gravity wave is always unstable in a plane defined by its wave vector and the gravity vector, whatever its initial amplitude. This instability is of the parametric type and is the analogue, but for a primary wave, of the parametric subharmonic instability (PSI) isolated by [START_REF] Mccomas | Resonant interaction of oceanic internal waves: the dynamic balance of internal waves[END_REF] in the Garrett-Munk spectrum. Energy is thus transferred toward motions of smaller vertical scale and lower frequencies. The linear stability analysis of [START_REF] Klostermeyer | Two-and three-dimensional parametric instabilities in finite amplitude internal gravity waves[END_REF] and [START_REF] Lombard | Instability and breakdown of internal gravity waves. 1. Linear stability analysis[END_REF] extends Drazin and Mied's study to primary waves of arbitrary amplitude: a statically stable (small amplitude wave) is most unstable to two-dimensional perturbations (in the plane defined by the wave vector and the gravity vector) but a statically unstable wave is most unstable to three-dimensional perturbations. It is noteworthy that the physical mechanism underlying the wave instability is the same whatever be its amplitude, i.e. of the parametric type. In the limit of a primary wave of vanishingly small amplitude, the resonant interaction theory is recovered. The numerical studies of [START_REF] Bouruet-Aubertot | Breaking of standing internal gravity waves through two-dimensional instabilities[END_REF][START_REF] Polzin | Statistics of the Richardson number: mixing models and finestructure[END_REF], [START_REF] Koudella | Breaking progressive internal gravity waves: two-and three-dimensional numerical experiments[END_REF] and [START_REF] Koudella | Ondes internes de gravité en fluide stratifié: instabilités, turbulence et vorticité potentielle[END_REF], while confirming these predictions, focus upon the nonlinear development of the instability, the breaking process and subsequent turbulent regime. They show that the instability of a high frequency internal gravity wave leads to transfers of energy to secondary waves of lower frequencies and that these secondary waves eventually break through shear instability. Our purpose here is to briefly describe the mechanisms leading a small amplitude (statically stable) primary wave to wavebreaking, owing to their possible role on mixing processes. Both two-and three-dimensional calculations will be described. In two dimensions, the evolution of the wave is illustrated by density and vorticity fields at different times, in Fig. 1 for propagating waves (run P1) and in Fig. 2 for standing waves (run S2). The initial amplitudes are the same.

The dynamics is first dominated by the oscillation of the primary wave (Figs. 1a and2a). For propagating waves, constant vorticity lines are straight lines which make an angle π/4 to the vertical, while for standing waves the vorticity field is that of a Green-Taylor vortex, as already noted. A weakly nonlinear regime first sets in during which a primary instability slowly grows through a mechanism of resonant triad interactions. As a consequence, an organized perturbation progressively develops, as attested by the deformation of the initially smooth isopycnals (Figs. 1c and2c). This perturbation actually results from a local parametric instability, in agreement with theoretical predictions as discussed above. Indeed, the perturbation is a wave oscillating with a frequency equal to half the primary wave frequency. The frequency of the perturbation is deduced from the mean angle of either the isopycnal or the vorticity bands with the horizontal. The vorticity bands make an approximate angle of 21 • to the horizontal in both cases very close to the theoretical value of 20.7 • .However, the mechanism producing the parametric instability is quite different in these two cases and depends on the nature of the wave, either propagating or standing. Basically the mechanism of parametric instability can be described by a parametric oscillator when the scale of the primary wave is large compared to that of the perturbation: the primary wave acts as a parametric forcing upon the perturbation, since its oscillation leads to a periodic variation of the frequency of the small scale perturbations which are advected by the large scale motion. From the dispersion relation ω = N cos(θ) = -(g/ρ 0 )[d( ρ(y) +ρ(x, y, t))/dy]cos(θ), two possible mechanisms of modulation of the frequency of the perturbation can be inferred: the first one is a variation of the buoyancy frequency induced by the primary wave (through d ρ/dy while the second one is a variation of the orientation, θ , of the wave vector of the perturbation advected by the large scale motion. The latter effect dominates for standing waves in the center of the domain where the perturbation is being advected into a rocking motion by the primary wave. Therefore, the perturbation is localized along a straight diagonal band centered at the extremum in stream function of the primary wave. In contrast the former effect is very likely responsible for the growth of a wave packet for propagating waves. In this case the growth of the perturbation arises in locations of positive vorticity of the primary wave [START_REF] Koudella | Ondes internes de gravité en fluide stratifié: instabilités, turbulence et vorticité potentielle[END_REF].

The amplitude of this perturbation grows exponentially until the density field locally overturns. This configuration appears not to be immediately dynamically unstable. In fact breaking seems to be initiated by the instability of the shear associated with the perturbation rather than by a convective instability. Indeed, the latter instability is very likely inhibited in the vertical plane due to the large scale shear induced by the primary wave [START_REF] Deardorff | Gravitational instability between horizontal plates with shear[END_REF]. The occurrence of wavebreaking through shear instability is characterized by the onset of fine scale structures of high intensity in the vorticity field. The latter instability first occurs in the middle of the domain for standing waves (Fig. 2), but, for propagating waves, manifests itself in ubiquitous parallel sets of unstable shear layers (Fig. 1e andf). For standing waves, local overturnings have extended over the whole domain during the latest stage and are still occurring despite the high decrease of the primary wave energy (Fig. 2g). This coincides with the extension of shear instabilities to the whole domain (see Fig. 2h). For propagating waves by contrast, no more overturning occurs within the flow (Fig. 1g).

Time evolution of the maximum of the vertical density gradient (Fig. 3c andd) and the kurtosis of the vorticity, Ω 4 (x, y) dx dy/ Ω 2 (x, y) dx dy 2 , where Ω(x, y) refers to the vorticity (Fig. 3a andb) are also displayed. These quantities provide evidence of changes in the dynamics, from quasi-periodic to turbulent, and of the occurrence of instabilities. Indeed positive values of the maximum of the vertical density gradient indicate local overturning while high values of the kurtosis of the vorticity denote shear instability, through which the primary wave breaks. During the quasi-periodic stage the kurtosis has a constant value of 1.5 for standing waves and 2 for propagating waves (a Gaussian field would have a kurtosis of 3 for comparison). A sudden increase occurs at the onset of wavebreaking, detected by positive values of dρ/dy (see Fig. 3c at t=32T and Fig. 3d at t=65T). For increasing times peaks of high values of the kurtosis of the vorticity are noted as long as intermittent shear instabilities occur (Fig. 3a andb). Eventually a decrease is observed during the latest stage of relaminarization of the flow. These results reveal a significant difference between standing and propagating waves: wavebreaking events occur over a much longer time interval in the case of standing waves than for propagating waves. This results from the spatial organization of the primary instability initiating breaking which is confined in the middle of the domain for a standing wave while more widely distributed in shear layers for propagating waves (see Figs. 1 and2). This difference in spatial organization of overturnings explains why energy is more rapidly dissipated for propagating waves than for standing waves. However, though the patterns of the instability differ, a unique relationship is obtained between the spatially averaged diffusivity and the turbulent Froude number of the flowwhateverthenatureofthewave(seeSection 5.2.3).

Prior to the analysis of mixing it is important to address the question of the relevance of two-dimensional numerical simulations to describe accurately the wave breaking process and subsequent wave turbulence. To this aim we compare the behavior of the propagating primary wave to that of the same wave subjected to a three-dimensional perturbation (see Section 3.1). This run is referred to as P10 in Table 1.

In three dimensions, the most unstable perturbations are two-dimensional and of parametric type, in agreement with the predictions by [START_REF] Klostermeyer | Two-and three-dimensional parametric instabilities in finite amplitude internal gravity waves[END_REF] and [START_REF] Lombard | Instability and breakdown of internal gravity waves. 1. Linear stability analysis[END_REF]. The analysis of run P10 reveals that shear instability occurs next, in the vertical plane of the primary wave, similar to the two-dimensional case. Convective instability eventually arises, after the onset of shear instability, leading to fully three-dimensional dynamics. However, when the wave is of a fairly small amplitude, the growth rate of the convective instability is so much smaller with respect to that of the shear instability that the mixing properties are appropriately described by a two-dimensional computation in this case.

It is worthwhile to compare the breaking of this primary wave of small amplitude with that of a statically unstable wave (e.g. [START_REF] Lombard | Instability and breakdown of internal gravity waves. 1. Linear stability analysis[END_REF][START_REF] Fritts | Vorticity dynamics in a breaking internal gravity wave. Part 2. Vortex interaction and transition to turbulence[END_REF]. In our case, the amplitude of the wave does not increase and wavebreaking occurs because the wave loses its coherence owing to the instability process. The mechanisms leading to breaking are mostly two-dimensional, as just discussed, and the overturning of isopycnals resulting from these instabilities may eventually trigger a three-dimensional buoyancy induced instability. By contrast, a statically unstable wave is most unstable to three-dimensional instability (the latter wave may result from upward propagation in the atmosphere, for instance). The resulting breaking is fully three-dimensional in this case (see [START_REF] Staquet | Mixing by internal gravity waves[END_REF] for a discussion). We recall, however, that the instability mechanism deriving its energy from the primary wave is of the parametric type, regardless of the amplitude of the wave [START_REF] Klostermeyer | Two-and three-dimensional parametric instabilities in finite amplitude internal gravity waves[END_REF].

Particle dispersion

This section deals with the breaking of a two-dimensional propagating wave. We determine the eddy diffusivity from the dispersion of fluid particles resulting from breaking.

Accuracy of the numerical method for tracking particle trajectories

The subroutines that compute the trajectories of particles have been previously tested for accuracy [START_REF] Babiano | Single particle dispersion, Lagrangian structure function and Lagrangian energy spectrum in two-dimensional incompressible turbulence[END_REF][START_REF] Zouari | Expériences numériques lagrangiennes à partir de modèles euériens[END_REF]. We therefore briefly describe only two tests that we performed.

We first considered the behavior of particles in a stationary unstratified flow. The initial condition is a Green-Taylor vortex (ψ(x, y) ∼ sin(k x x)sin(k y y)), which is a stationary unstable solution of the Navier-Stokes equations in a homogeneous fluid. Two straight lines of particles were introduced at the initial state and we checked that particle trajectories remained on streamlines and described closed circles. The precision of the method is equal to 6×10 -13 after 120 eddy turnover times.

We then considered a single propagating wave, which is an exact solution of the nonlinear equations. For a periodic wave of wavevector (k x , k y ), particle trajectories follow straight segments, verifying the equation k x x+k y y=constant (with k x , k y =1). The deviation of particle trajectories from the analytic solution is of the order of the accuracy of the computer after 15 000 time steps (i.e. 6.75 wave periods).

Results

Particles were released in the flow during the turbulent stage in order to provide an estimate of vertical dispersion induced by wavebreaking. The time of release of particles is the time of the first local maximum of the kurtosis of the vorticity. This time corresponds to the onset of shear instability and wavebreaking.

We focus on the case of horizontal lines of particles because this distribution of particles appeared to be the most relevant for diapycnal mixing. A distribution into clouds of particles will be mentioned for comparison. Three lines of particles have been released in the flow at time t=36.62, in units of buoyancy periods (run P1), that is just after the beginning of wavebreaking. Each line is made of 3072 particles (i.e. one particle every dx/6). The time evolution of these sets of particles is displayed in Fig. 4 at successive times. The corresponding vorticity fields are shown as well (Fig. 2, second column). The vertical spreading of the particles is first highly variable along the horizontal direction as shown by the very different shapes of the lines, which are either horizontally elongated or vertically dispersed (see Fig. 4c at t=36.62+1.6). These different behaviors result from the spatial inhomogenities of the velocity field. This is clearly evident on the second line of particles at t=36.62+1.6 (Fig. 4c): the left half of the line is dispersed in a zone of small intense eddies corresponding to an unstable region (see Fig. 4d) while the right half is strained under the action of long vortices of small vertical scales, corresponding to a zone of stability. For increasing time, t>36.62+4.8, the vertical dispersion of the lines becomes homogeneous over the horizontal direction. This horizontally homogeneous behavior starts when the width of the lines becomes larger than the typical vertical extent of these small scale intense eddies so that the effects of these small scale structures of the velocity field are averaged. The vertical spreading of each line is summarized in Fig. 5 where vertical profiles of particle concentration are plotted. These vertical profiles are obtained by horizontally averaging the number of particles per grid cell (Fig. 5). The vertical profile evolves progressively from three delta functions to three Gaussian like shapes. As expected the maximum of particle concentration decreases as the lines vertically disperse and the width of the peak increases characterizing the dispersion.

Our first attempt to evaluate a dispersion coefficient, K, is to infer it from the equation of evolution of the vertical profiles of particle concentration. This method appears to be unreliable due to the spiky shape of the profiles. A preliminary fit of these curves to a Gaussian function is no more helpful since some of the curves are highly non symmetric with respect to the maximum concentration. Preferably the dispersion coefficient is inferred from two quantities characterizing the variance of the vertical distribution of particles [START_REF] Hunt | Diffusion in the stably stratified atmospheric boundary layer[END_REF]. The first quantity, denoted Z2 0 , is the mean square of the vertical distance from the center of gravity for all the particles of a given line: Z2

0 (t) = 1/N N i=1 [y i (t) -y G (t)] 2
, where y i (t) is the vertical position of the particle i and y G (t) the vertical position of the center of gravity of the set of particles at time t; the second quantity, denoted Z2 sq , is the mean square vertical deviation from the initial distribution located at y 0 for a given line: Z2 sq (t) = 1/N N i=1 [y i (t)y 0 ] 2 . Time evolutions of these two quantities are displayed in Fig. 6. A logarithmic plot (not shown here) enhances the distinction between the two regimes, the first regime dominated by shear effects and the second one exhibiting a linear law that is Z2 0 (t) and Z2 sq (t) ∼ t. The first regime is characterized by a rapid and fluctuating increase of the variance. This behavior is observed for each line of particles though noticeable differences are seen. These differences result from the local influence of small scale eddies as already mentioned. Later as the width increases the influence of small scale eddies is averaged. As a result an almost linear regime is reached for Z2 0 (t) and Z2 sq (t). We note as well, residual oscillations resulting from wave motions superimposed on the mean evolution. We would like to stress that the spatial distribution of the particles is crucial to the concentration becoming statistically homogeneous horizontally before the end of the turbulent regime. Indeed, for an initial distribution of fluid particle into discrete clouds, no such behavior is obtained. The main explanation is that, because of the small horizontal extend of the cloud (∼10 dy), the time evolution of the set of particles is dictated by motions of small scale structures and does not become independent of these structures before the end of the wave turbulence regime.

Theoretical predictions for Z(t) 2 have been derived by [START_REF] Pearson | A statistical model of fluid-elements motions and vertical diffusion in a homogeneous stratified turbulent flow[END_REF], for a stably stratified statistically stationary turbulence. The behavior of Z2 (t) depends on whether fluid particles exchange their density: when there is no diffusion, the model predicts that Z(t) 2 ceases to grow after a time of the order of the buoyancy period and reaches an asymptotic value, which is a decreasing function of the Froude number of the turbulence (defined as the ratio between buoyancy time scale and velocity time scale). When there is a small diffusive change of density, the mean-square vertical displacement eventually displays a slow linear growth instead of reaching an asymptotic value. The growth rate is proportional to the eddy diffusivity. In the present case, though breaking is not a stationary process, the statistical properties of the flow remain about the same during a few buoyancy periods and this is attested for instance by saturated kinetic and potential energy spectra remaining at the same level during breaking [START_REF] Bouruet-Aubertot | Stratified turbulence produced by internal wave breaking[END_REF]. The predictions of the model may therefore be tested. The trajectories of the particles in our computations are driven by the velocity field only but this velocity field is that of fluid parcels that are subjected to density diffusion. Therefore, the particles can be interpreted as fluid particles. The predictions of the model should thus be tested for the case where particles exchange density. Fig. 6 reveals a good agreement between these predictions and our numerical finding. Indeed, Z(t) 2 displays a slow growth, after four buoyancy periods.

We compute the dispersion coefficient K from the slope of Z2 (t) [START_REF] Hunt | Diffusion in the stably stratified atmospheric boundary layer[END_REF], where Z=Z 0 or Z sq . In our simulations, because of residual oscillations, K is calculated from the slope of the least square line (see Fig. 6). Values for K between 0.0025 and 0.0034 are obtained. To allow for comparisons with other works, we divide K by the molecular diffusivity of our numerical model. This ratio defines the Cox number

(K = (1/2)(d/dt) ( Z2 (t) 2 ); see
Cox = K κ ( 12
)
where κ=ν/Pr. The Cox number is indicative of the increase in diffusion due to turbulence since for a periodic wave, i.e. for a laminar regime, Cox≈1. The values obtained for K yield the range 12-16.3 for the Cox number, corresponding to a variation of the order of 30%. This variation could very likely be reduced by considering a forced primary wave, a case for which particle experiments could be conducted over a longer time interval.

Comparisons with other works

We first compare our results with other two-dimensional numerical simulations of internal gravity waves [START_REF] Schilling | Vertical mixing of passive scalars owing to breaking gravity waves[END_REF] and Kelvin-Helmholtz waves in atmospheric jets [START_REF] Schilling | Particle dispersion due to dynamical instabilities in the lower atmosphere[END_REF]. Their values for K are much higher than ours, the effective dispersion coefficient lying within the range 0.3-0.7 m 2 s -1 . Also, the dispersion coefficient is highly variable, both spatially and temporally (see Fig. 6 of [START_REF] Schilling | Particle dispersion due to dynamical instabilities in the lower atmosphere[END_REF] for instance). These results illustrate the difficulty of evaluating K in flows which are not in statistically steady states such as the turbulent regime resulting from a wave breaking.

Dispersion by freely decaying stratified turbulence was recently studied by [START_REF] Kimura | Diffusion in stably stratified turbulence[END_REF] using three-dimensional direct simulations of 128 3 resolution. The time evolution of Z2 sq showed a smoother evolution since only slight oscillations could be noted for the highest values of N (Fig. 1 of [START_REF] Kimura | Diffusion in stably stratified turbulence[END_REF]). This results from the buoyancy flux which is more significant in wave-turbulence flows as opposed to stratified turbulence.

We eventually compare our results for the vertical eddy diffusivity with those obtained through tracer release experiments in lakes [START_REF] Wuest | Comparison of diapycnal diffusivity measured by tracer and microstructure techniques[END_REF]. These authors found high discrepancies between boundary and internal diapycnal mixing: values of eddy diffusivity lie within 7×10 -6 and 8×10 -5 m 2 s -1 near the boundaries while these in the interior are up to 4×10 -6 m 2 s -1 . Corresponding Cox number (i.e. K/κ) are [50, 571] and 28, respectively (with κ=1×10 -7 and 4×10 -7 m 2 s -1 ). Our value of 16 obtained for the wave of the smallest amplitude is in reasonable agreement with the value obtained in the interior.

Mixing

Background

We have applied the analysis of mixing in stably stratified fluids proposed by [START_REF] Winters | Available potential energy and mixing in stratified fluids[END_REF] and [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF]. The method relies on the distinction between two kinds of potential energies, available and background. The available potential energy is the reversible part of potential energy associated with waves and periodically exchanged with kinetic energy; the background potential energy is the irreversible part of potential energy associated with fluid mixing. When all fluid motions have ceased, the potential energy solely consists in the background potential energy. It is therefore the minimum potential energy of the fluid. These ideas were proposed by [START_REF] Lorenz | Available potential energy and the maintenance of the general circulation[END_REF] and applied by [START_REF] Thorpe | Turbulence and mixing in a Scottish loch[END_REF] to experimental data in one dimension. [START_REF] Winters | Available potential energy and mixing in stratified fluids[END_REF] and [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF] have derived a mathematical framework that formalizes these concepts in three spatial dimensions.

We first provide the expression of the diapycnal diffusive flux from which the eddy diffusivity is inferred. We next derive potential energy budgets, following [START_REF] Winters | Available potential energy and mixing in stratified fluids[END_REF].

An exact expression for the diapycnal diffusive flux

As recalled in Introduction, when the vertical mean density gradient varies smoothly over the scale of turbulent motion, the diapycnal flux is assumed to depend linearly upon that gradient. The proportionality constant is the diapycnal eddy diffusivity. The determination of the diapycnal flux is difficult from a practical point of view because it amounts to extracting a diffusive contribution from a turbulent field. Models are therefore usually employed to compute it (see e.g. [START_REF] Gregg | Diapycnal mixing in the thermocline: a review[END_REF] and [START_REF] Toole | Turbulent mixing in the ocean: intensity, causes, and consequences[END_REF] for a review). A new method has been recently proposed by [START_REF] Winters | Diascalar flux and the rate of fluid mixing[END_REF] to compute that flux exactly. The key ingredient is a stable density profile obtained by reordering the fluid parcels adiabatically so that the heaviest fluid particle has the lowest altitude, following the method described by [START_REF] Thorpe | Turbulence and mixing in a Scottish loch[END_REF]. The expression for the diapycnal flux is

φ d (y, t) = κ (∇ρ) 2 I -dρ * /dy ( 13
)
where κ is the molecular diffusivity, the index I refers to an average along an isopycnal and ρ * is the sorted profile. Note that the diapycnal flux thus defined is always positive. Also, its expression is formally similar to that proposed by [START_REF] Osborn | Oceanic finestructure[END_REF], but using an average along the isopycnals and a sorted density profile in place of a horizontal average and a mean density gradient. Following these lines, the eddy diffusivity should be defined as

K d (y, t) = φ d (y, t) -dρ * /dy (14)
Note that K d is always positive, whatever be the flow dynamics.

Budgets of background potential energy and available potential energy

The background potential energy per unit volume E b is simply the potential energy of the background (or sorted) density profile

E b (t) = g V V yρ * (y, t) dV (15)
It can also be expressed in terms of the total density field (since the sorting process is an adiabatic redistribution of mass)

E b (t) = g V V y * ( x, t)ρ( x, t) dV ( 16
)
where y * ( x, t) is the vertical position, in the state of minimum potential energy, of the fluid element located before sorting at position x at time t.

The energy released during the sorting process is the available potential energy. It is equal to the difference between the total potential energy and the background potential energy

E a = E p -E b (17)
that is

E a (t) = g V V (y -y * ( x, t))ρ( x, t) dV (18)
These equations provide a general definition for the available potential energy, valid for whatever be the flow conditions. For linearized fluid motion, in which case the sorted density profile remains close to some initial or reference profile because density fluctuations about this profile are small, the available potential energy can be approximated by the variance of the density fluctuations. Note that closed analytical expressions of the local available potential energy have been derived by [START_REF] Holliday | On potential energy density in an incompressible, stratified fluid[END_REF].

The volume averaged potential energy E p satisfies the equation

d dt E p = Φ a -Φ i + Ėp surf adv + Ėp surf dif (19)
while the equation for the background potential energy is

d dt E b = Φ d + Ėb surf adv + Ėb surf dif (20)
We use the following notation: For both propagating and standing waves, any change of E b can result either from diapycnal diffusion within the fluid interior or from advective or diffusive processes at the boundaries. Note that for a closed system with adiabatic boundaries, E b can only increase under the action of diapycnal diffusion within the fluid interior since surface fluxes vanish.

• Φ d = (g/H )
The evolution equation for the available potential energy per unit volume is derived from Eqs. ( 19) and (20), using Eq. ( 17)

d dt E a = Φ a -(Φ d -Φ i ) + ( Ėp surf adv -Ėb surf adv ) + ( Ėp surf dif -Ėb surf dif ) (21)
Thus, surface fluxes cancel if there is no density overturning at the boundaries that is if y * =y.I nt h efluid interior, the time rate of change of available potential energy depends fundamentally on the relationship between the buoyancy flux Φ a and the turbulent diapycnal flux Φ d -Φ i .

Results

Temporal evolution of the Cox number: comparison with the values from particle dispersion

The purpose of this section is to provide evidence of the correlation between mixing processes and the flow dynamics and to compare the value of the maximum Cox number obtained from mixing analysis and from particle dispersion.

Time evolution of the Cox number, K d /κ, is displayed in Fig. 7a for standing waves (run S2) and Fig. 7b for propagating waves (run P1). Both the spatially averaged Cox number and its maximum value are plotted in Fig. 7. During the first stage, which is dominated by the periodic oscillation of the primary wave, the Cox number is equal to 1+o(Fr 2 ), where Fr is the Froude number of the primary wave (defined in Section 2.1); it has therefore a value close to 1. The onset of wavebreaking results in an increase of the Cox number. As long as overturnings occur, its evolution is characterized by the occurrence of peaks of high values. A decrease is eventually observed for propagating waves during the latest stage of the flow when most of the energy of the primary wave has been dissipated through breaking (Fig. 7). This temporal evolution displays the same qualitative behavior as those of dρ/dy and of the kurtosis of the vorticity presented in Fig. 3. More precisely maximum values are associated with strong wavebreaking events which are denoted by high positive values of dρ/dy and high value of the kurtosis of the vorticity (Fig. 3a andc at t=65andFig.3b anddatt=32). Comparisons between propagating and standing waves of same amplitude show that slightly higher values of the Cox number are obtained for propagating waves. The reason is that the Cox number depends upon the energy of the primary wave (and not upon its amplitude), during the initial quasi-periodic regime, and the energy of the propagating wave per unit area is twice the energy of the standing wave per unit area.

Comparison with values of the Cox number obtained from particle dispersion reveals a good agreement with the present values of this parameter. Values obtained from particle dispersion are within 12 and 16.3 for K/κ, for times within 40 and 46T. These values well compare with the value of the maximum Cox number at time t=40T which is of 11. This is one of the main results of the paper. It states that the vertical diffusivity can equivalently be inferred from particle dispersion experiments or from the direct computation of the diapycnal flux. However, it would be valuable to perform simulations of a forced wave for further comparison between the two methods.

The energetics of mixing

We describe in this section the budgets of background and available potential energy and a time-averaged mixing efficiency is defined.

The budget of background potential energy.

The temporal evolutions of the change in the background potential energy, E b (t)-E b (0), normalized by the change in the total wave energy over the duration of the calculation, E wave is represented in Fig. 8a andb for standing waves and propagating waves, respectively. E wave is defined by E wave (t=t end ) -E wave (t=0), where E wave (t) is the sum of the kinetic and available potential energy at time t. This normalized change in background potential energy will be referred to as E * b in the following. The time integral of the sum of Φ d and Ėb surf dif is also plotted in Fig. 8b. For standing waves the time integral of the fluctuating diffusive flux, Φ d -Φ i , is displayed in Fig. 8c. The time integrated fluxes in these figures are also divided by E wave .

The overall evolution of E * b displays two main stages, whatever the propagating or standing nature of the primary wave. E * b is first nearly constant during the quasi-periodic regime. By the end of this regime (t∼60T for standing waves and t∼28T for propagating waves), t h ei n c r e a s ei nE * b is less than 5% of the total variation in the wave energy (Fig. 8a andb). E * b starts growing significantly from the onset of wavebreaking (at t∼65T for run S2), seeFig.8aandatt∼35T for run P1, see Fig. 8b). As shown by Eq. ( 20 potential energy. For propagating waves, this sum describes only the mean evolution of this energy since the advective surface flux Ėb surf adv adds an oscillatory component. This is clearly illustrated in Fig. 8b. These oscillations, however, become significant only when wavebreaking begins. We now describe in more detail the behavior of E * b for the propagating and standing wave cases separately.

The density difference between the upper and lower boundaries is maintained constant in our study; thus, if the fluid were at rest (in which case N=N lin ), the vertical flux of matter associated with diffusive processes should be made up for near the boundaries in order not to mix the fluid there. It follows that, because of the constant value of N lin , the surface diffusive flux should exactly oppose the diapycnal flux within the fluid volume: Ėb surf dif =-Φ i .We find that this relation actually holds all along the fluid evolution for the propagating case. It follows that the mean evolution of E b for propagating waves can be interpreted as solely due to the fluctuating diapycnal flux Φ d dt -Φ i dt. During the quasi-periodic regime, theincreaseofE b is fairly small because the fluctuating diapycnal flux varies like A 2 , where A is the amplitude of the periodic wave. As wavebreaking occurs, the diapycnal flux Φ d increases (while Φ i retains a constant value), which leads to a significant growth in E b . During the regime subsequent to breaking, the fluctuating diapycnal flux nearly vanishes because no more breaking occurs and E b ceases to grow (Fig. 8b).

For standing waves, the surface advective flux vanishes and, in addition to the constancy of N, ρ ′ issettozeroatthey-boundaries. This implies that, at the y-boundary, the role of the surface diffusive flux is to compensate any change in density at this boundary. It follows that, during the quasi-periodic regime, this flux should oppose the diapycnal flux. Fig. 8a shows that this is indeed the case: the background potential energy does not grow during the quasi-periodic regime. When breaking occurs, E * b displays a rapid growth but soon reaches a plateau. The reason is very likely that breaking mostly occurs in the center part of the domain and increases the diapycnal flux there, while the surface diffusive flux remains nearly unchanged. Time is needed for this diapycnal flux to also occur at the boundaries, and this accounts for the delay in the growth of the surface diffusive flux, as manifested by the plateau eventually reached by E * b . The surface diffusive flux thus acts as a local heat flux for the upper boundary and cooling for the lower boundary. E b resumes its growth from ∼105T due to a new burst of overturning activity; this is attested in Fig. 3 by the occurrence of a local maximum of the kurtosis of the vorticity. Eventually E b decreases since diffusion at the boundaries dominates because of less turbulent activity (see Fig. 3). It is expected that when all the energy put into the system, that is the wave energy, has been dissipated, E b will recover its initial value since the density difference across the domain height ρ top -ρ bottom remains constant. This is only the case for standing waves for which surface diffusion at the upper and lower boundaries causes a decrease in background potential energy which will eventually compensates for the increase in background potential energy resulting from mixing.

The budget of available potential energy.

The time evolution of E a is displayed in Fig. 9. The result of breaking is a steep decrease of E a (Fig. 9a andb). Before a detailed description we would like to underline a difference between standing and propagating waves which is crucial when mixing is estimated from the buoyancy flux. The volume averaged E a is oscillating in time in the case of standing waves (Fig. 9a) while this is not the case for propagating waves (Fig. 9b). This result can best be seen in terms of the buoyancy flux: for a propagating wave this flux vanishes when integrated over the whole domain, whereas for a standing wave the buoyancy flux oscillates when volume averaged. We find, in contrast to E b , that the evolution of E a is only slightly modified during breaking by surface fluxes. For propagating waves for which the time integrated buoyancy flux is small (see Fig. 9b), the evolution of E a follows that of -(Φ d -Φ i ) dt. In this case the rate of decrease of E a is close in absolute value to the fluctuating (turbulent) diffusive flux. Note that the rate of decrease of E a is not exactly zero during the periodic stage since Φ d dt> Φ i dt (see Eq. ( 27)). By contrast, for standing waves, the budget of the different terms contributing to the E a evolution is very different since the time integrated buoyancy flux is an important component (Fig. 9a). In this case the turbulent diffusive flux can not be inferred directly from the rate of decrease of E a and requires careful time averaging when fluxes Φ d and Φ a are not available. 22)) during the turbulent regime for runs P1 (Pr=0.72), P6 (Pr=1), P7 (Pr=5)andP8(Pr=10). The dependency for a linear wave, 1/(1+Pr) is also plotted (dashed line).

Mixing efficiency.

A relevant quantity regarding the energetics of mixing is the mixing efficiency which allows one to calculate the energy used to mix the fluid given the input of energy in the system. In what follows we focus on mixing within the interior of the fluid. This point of view is justified by the fact that breaking occurs within the fluid interior in localized overturning events of scales much smaller than the domain height. We define a time-integrated mixing efficiency during each regime, quasi-periodic and turbulent. This quantity consists in the ratio of time integrated energy fluxes, namely the fluctuating diapycnal flux divided by the instantaneous change in the wave energy:

η(t) = t t 0 (Φ d -Φ i ) dt [ E wave ] t t 0 (22)
where t 0 refers either to the initial time when η is calculated during the quasi-periodic regime or to the end of the quasi-periodic regime when η is calculated during the turbulent regime.

We have computed η for different values of the Prandtl number (0.72, 1, 5 and 10), in the case of propagating waves. Similar values of mixing efficiency are obtained for standing waves. η has a constant value during each regime. During the quasi-periodic regime, we have checked that the dependency of the mixing efficiency upon the Prandtl number follows that of a periodic wave, namely η=1/(1+Pr). The dependency of η upon the Prandtl number during wavebreaking is displayed in Fig. 10. The Prandtl number dependence for a periodic wave (1/(1+Pr)) is also displayed for comparison. The evolution of η(Pr) for Pr<1.5 follows that of a periodic wave showing that wave diffusion is dominant in these computations. Then, for increasing Prandtl number, η(Pr) differs from 1/(1+Pr) and the decrease of η is much weaker. This evolution is consistent with the fact that the dynamics is more turbulent for increasing Prandtl number. Moreover, it appears that η(Pr) seems to reach an asymptotic value in a similar way as turbulent flows, in which mixing efficiency does not depend on the molecular Prandtl number. This asymptotic value of η would compare quite well with the usual maximum value of about 0.2-0.3 (if we refer to the value of 0.25 which is obtained for Pr=10) obtained in the ocean as in laboratory experiments [START_REF] Linden | Mixing in stratified fluids[END_REF]. However, the effect of the restriction of our simulations to a vertical plane is to be addressed in order to assess this result. This point is further discussed in Section 5.2.5 in which a comparison with preliminary analysis of three-dimensional calculations is conducted.

Parameterization of mixing versus the wave parameters and background stratification

We focus on the evolution of the spatially averaged Cox number ( Cox( x, t) x ), referred to as Cox(t) , versus the turbulent Froude number. The Cox number is defined by Eq. ( 12). 23)) for all runs of Table 1 (except runs P6-P10, at all discrete times). The dotted line is the straight line: Cox∼1.12+2.78(Fr t -Fr wave ) 2 , where Fr wave =Fr t (0).

We define the turbulent Froude number as

Fr t = √ 2E Z N (23)
where E Z is the instantaneous enstrophy. This definition of Fr t is equivalent to the definition proposed by [START_REF] Ivey | On the nature of turbulence in a stratified fluid. Part I. The energetics of mixing[END_REF] for the turbulent Froude number, Fr t =(u/L c )/N, where u denotes the RMS of the velocity of turbulent motions and L c is the scale of the most energetic events. Indeed, with the kinetic energy spectra of both propagating and standing waves behaving as k -3 (Bouruet- [START_REF] Bouruet-Aubertot | Stratified turbulence produced by internal wave breaking[END_REF][START_REF] Koudella | Breaking progressive internal gravity waves: two-and three-dimensional numerical experiments[END_REF], it follows that the enstrophy E Z behaves as k -1 so that the value of Fr t is controlled by the largest, and thus most energetic, scales. Note that the associated gradient Richardson number 1/Fr 2 t is close to the minimum Richardson number. In Fig. 11 is displayed Cox(t) at all discrete times of the runs of Table 2, versus (Fr t -Fr wave ) 2 ; Fr wave is the initial Froude number of the wave and is defined by √ 2E Z (t = 0)/N. Note that the Froude number of the wave comes into play since the value of the Cox number during the initial quasi-periodic stage depends upon the wave amplitude. Fig. 11 shows that the behavior of Cox(t) can be described using the law Cox(t) ∼1.12 + 2.79 × (Fr t -Fr wave ) 2 ( 24)

with the correlation coefficient of 0.978. This law can be accounted for by using the fact that, for run P1 for instance, ǫ and Φ d behave very closely (not shown). Assuming a linear dependency between these two quantities, of the form Φ d =λǫ+µ, one gets Cox = λ ′ Fr 2 t + µ ′ with λ ′ =λN 2 ν/(2Φ i )a n dµ ′ =µ/Φ i (Φ i having a constant value, as already noted).

Eq. ( 24) is valid whatever be the nature of the wave, standing or propagating. This shows that the averaged eddy diffusivity does not depend upon the details of the spatial organization of the instability (see Section 3.2). It is useful to note that this relationship is similar to that obtained for a two-dimensional stratified shear layer [START_REF] Staquet | Mixing in a stably-stratified shear layer: two-and three-dimensional numerical experiments[END_REF]. This similarity is consistent with the process of wavebreaking to eventually occur through shear instability.

Interestingly the relationship between vertical eddy diffusivity and the squared Froude number is analogous to mixing models based upon a linear dependency between the dissipation rate of turbulent kinetic energy and the finescale shear level of the internal wave field in the ocean [START_REF] Polzin | Statistics of the Richardson number: mixing models and finestructure[END_REF]. This finescale shear level is dictated by the value of the shear at the scale which limits the Garrett and Munk spectral domain (constant spectrum of the vertical shear) and the buoyancy subrange (k -1 z vertical shear spectrum). More generally this parameterization is similar to empirical parameterization used in large scale models for instance (e.g. [START_REF] Pacanowski | Parameterization of vertical mixing in numerical models of tropical oceans[END_REF].

The structure of the background density profile resulting from mixing

To get a direct view of mixing, we display in Fig. 12 (for run S2) and in Fig. 13 The relative change δρ b (y, t) reaches at most 5% for the case of highest amplitude. As shown by Figs. 12 and 13, the main feature of δρ b is the succession of positive and negative values. When the slope of δρ b (y, t) is positive, the buoyancy frequency is smaller than the initial value and conversely for positive values as seen from the buoyancy frequency profile. This effect, though very small, is the signature of the tendency of the flow to organize into layers. This layering, as we shall refer to, is a direct consequence of the breaking process since this process results from the formation of horizontal shear layers that ultimately become unstable. This flow organization is observed both for standing and propagating waves, although differences can be noted. Indeed, while the amplitude of the distortion of the background density profile is maximum at half depth for standing waves, this amplitude is fairly constant with depth for propagating waves. This feature which results from the difference in the spatial structure of the primary instability of the wave (see vorticity fields Figs. 1d and2d) has no effect on the spatially averaged diffusivity (see previous section and Fig. 11).

Corresponding profiles of the Cox number (Figs. 12 and 13) show a similar organization into layers, of alternatively higher and smaller eddy diffusivities. Note that the same behavior is observed in mixing layers where maximum eddy diffusivities are found in the zone of maximum shear [START_REF] Staquet | Mixing in a stably-stratified shear layer: two-and three-dimensional numerical experiments[END_REF]. This suggests that low amplitude (statically stable) breaking internal gravity waves could be modeled as a succession, along the vertical direction, of shear layers separated by quieter regions.

It might appear surprising that the changes of the background density profiles are fairly weak, whereas a high mixing efficiency is obtained. The reason is that mixing efficiency only relates to the wave energy, i.e. to the energy input to the system, without reference to the total amount of potential energy. The wave energy is much smaller than the energy that would be required to fully mix the fluid in which it propagates. In other words, instabilities will occur on a vertical scale that is much smaller than the scale of the density gradient. Efficient mixing may thus occur.

Therefore, the reason for this weak distortion of the background density profile is that we study the effect of mixing over a fairly small time interval, of a few tens of buoyancy periods, corresponding to the time of decrease of the energy of the unforced wave. It would be interesting to study the effects of mixing over longer times for a forced wave. In order to do this we estimate the time required to get a step like structure of length scale λ v .W e assume that the flow is forced and that a steady state has set up during wavebreaking. The rate of input of energy is thus opposite to the rate of decrease of the wave energy during wavebreaking, that is

ǫ i =- dE wave dt (25)
T h ev a r i a t i o no fe n e r g yi nt h es y s t e m , E wave , required for an increase in background potential energy, E b , depends upon the mixing efficiency, η, as follows:

E b = η E wave ( 26
)
where

E wave = ǫ i T (27)
for a constant rate of input of energy, ǫ i . Note that since molecular diffusion almost cancels surface diffusion, η is close to the mixing efficiency defined in Section 5.2.4, that is

η = (Φ d -Φ i ) dt E wave (28)
The variation in E b associated with the evolution from a linear background density profile to a step-like profile of scale λ v is

[ E b ] T 0 = H λ v g ρ 0 d ρ dz λ 2 v 12 ( 29
)
where E b is expressed per unit of mass and with H being the height of the domain. Note that the variation in E b is smaller within a factor of (λ v /H) 2 than that to get a full mixing of the fluid. This leads to the following expression of T versus the scale of the layers:

T ∼ H λ v g ρ 0 d ρ dz λ 2 v 12ηǫ i (30)
T is estimated for run P1, that is with the following numerical values:

ǫ i =6×10 -4 [L 2 ][T 3 ],
λ v /H=16 and η=0.25. We get T=3160 in units of N -1 . This result can be easily transposed to an oceanic internal gravity wave of frequency N/ √ 2 such that the Froude number, Fr=U/(Nλ), is equal to 0.04, that of run P1. Note that the time T depends upon the rate of decrease of the energy, ǫ i , and the mixing efficiency. Therefore, since the rate of decrease of the energy does not depend upon the Reynolds number, provided that the Reynolds number is sufficiently high (see Section 2.1), our result still holds for higher Reynolds number and the similarity in Reynolds number is not a compulsory condition. The application of our results to an oceanic wave of wavelength equal to 20 m evolving in a linear stratification of N=2.5×10 -3 rad s -1 gives a time T of about 14 days. Thus, once the secondary waves of lower frequencies have become unstable, a time T∼14 days is required to fully mix the layers. Note that this time would only be of 2 days if the Froude number is increased by a factor of 2. These estimates show that even for the lowest wave amplitude the time required for reaching a step-like structure is consistent with oceanic timescales. Moreover, they can be related with the scenario proposed by [START_REF] Eriksen | Measurements and models of fine structure, internal gravity waves, and wave breaking in the deep ocean[END_REF]: wavebreaking occurs in regions of high shear associated with low frequency waves. Here, in our calculations, these regions of high shear are waves of lower frequency (compared to the primary wave frequency) whose energy arise from the parametric instability of the primary wave.

Comparison with a preliminary three-dimensional calculation

The purpose of this section is to determine how mixing is affected by the restriction of the dynamics to a vertical plane. To this aim we analyze mixing induced by the same propagating wave as in run P1 but whose dynamics is resolved in three dimensions, using a high resolution (see [START_REF] Koudella | Breaking progressive internal gravity waves: two-and three-dimensional numerical experiments[END_REF]. The resolution of the three-dimensional simulation is 256 3 (as opposed to 512 2 for run P1) and the Reynolds number is twice smaller compared to run P1. Therefore, for consistency, we compare the three-dimensional case to its two-dimensional analogue with same resolution and Reynolds number.

The time evolution of the total energy is displayed in both cases in Fig. 14a. Very similar evolution is obtained except that wavebreaking occurs slightly earlier in two dimensions and this feature is attributed to the lowest degree of freedom in this case. We now focus upon the mixing properties of both flows.

The temporal evolution of E b , relative to the total change in the wave energy during the flow evolution, is plotted in Fig. 14b for both two-and three-dimensional simulations. Very similar curves are displayed. The effect of the smaller Reynolds number (compared to run P1) should be noted: the increase of E b resulting from the pure molecular diffusion is of 15% at the onset of breaking while it is only of 0.5% for run P1 (see Fig. 14). At the end of the turbulent regime the relative change in the background potential energy is of about 55% in both cases. Therefore, the mixing properties of both flows, in terms of the global increase of potential energy associated with mixing, are similar.

Regarding mixing efficiency, very close values are also obtained in the two-and threedimensional cases. During wavebreaking the mixing efficiency reaches a value of 0.55 in both cases. We explain this similarity by the fact that the dynamics of the primary wave, which is of low amplitude, remains quasi-two-dimensional in the three-dimensional calculations (except for a weak convective instability after the onset of shear instability), in agreement with linear stability analysis. The high value of the mixing efficiency during the breaking regime also indicates that the increase in the background potential energy is strongly contributed by purely molecular diffusive effects, which act similarly in two and three dimensions. It is very likely that the mixing efficiency during wavebreaking will decrease as the amplitude of the wave will increase because the dynamics become fully three-dimensional. This point is under present investigation.

Conclusion

In this paper we provide estimates of diapycnal mixing induced by breaking internal gravity waves that are either standing or propagating. We focus on the case of a primary wave of fairly high frequency (N/ √ 2) and of small amplitude, that is, which is statically stable. This primary wave of high frequency is unstable toward perturbations of lower frequencies. These perturbations which are secondary waves of half the primary wave frequency take their energy from the primary wave through a mechanism of resonant interactions. These secondary waves of lower frequency then become unstable and break. This process of destabilization of a high frequency wave provides insights on how a high frequency internal gravity wave can have an important role in transferring energy towards lower frequency motions. Note that these processes of destabilization of a primary wave are not altered by the presence of a shear flow of much lower frequency provided that the primary wave travels in a plane normal to the lower frequency shear flow [START_REF] Bouruet-Aubertot | Numerical experiments on internal gravity waves in an accelerating shear flow[END_REF]. This study might be of interest in the ocean where waves of fairly low frequency are energetic and whose breakdown is responsible for significant mixing. We analyze the results of two-dimensional numerical simulations and we compare these results with those obtained from a three-dimensional simulation.

We show how our results can be transposed to an oceanic internal gravity wave provided that the Froude number of the oceanic wave is within the range of the Froude numbers here considered. Indeed the dynamics itself depends only upon the Froude as this has been shown previously [START_REF] Bouruet-Aubertot | Breaking of standing internal gravity waves through two-dimensional instabilities[END_REF]. It results that even if the similarity in Reynolds number is not be achieved (the Reynolds number being slightly higher in the ocean), the rate of decrease of the energy (kinetic+potential) does not depend on the Reynolds number.

We first provide a direct estimate of vertical eddy diffusion from particle dispersion. A local dispersion coefficient is inferred from the time evolution of sets of particles which are released in the flow during the turbulent regime. We show that after an early evolution dominated by shear effects, a diffusion law is eventually reached. The vertical dispersion coefficient is computed and found to be fairly independent of the location of the particles in the flow.

The second method for computing the vertical (actually diapycnal) eddy diffusivity relies on potential energy budgets [START_REF] Winters | Available potential energy and mixing in stratified fluids[END_REF], following ideas previously introduced by [START_REF] Lorenz | Available potential energy and the maintenance of the general circulation[END_REF]. In this method the diapycnal eddy diffusivity is inferred from the diffusive flux across the isopycnals of the background density profile. A good agreement with the value inferred from particle dispersion is obtained. However, further calculations but of a forced wave would be valuable for a more detailed comparison. This agreement is of particular importance regarding the interpretation of in situ values either inferred from microstructure measurements or from tracer dispersion. The range of variation of the maximum over time of the Cox number (i.e. the eddy diffusivity relative to the molecular diffusivity) is [9,262] and, for the spatially averaged Cox number [5,60]. The former range is consistent with recent measurements of vertical diffusion in the interior of the ocean, Cox∼100 [START_REF] Ledwell | Evidence for slow mixing across the pycnocline from an open-ocean tracer-release experiment[END_REF].

We next define a mixing efficiency in which the increase of potential energy resulting from mixing is compared to the wave energy loss. The mixing efficiency during wavebreaking is Prandtl number dependent which attests of the contribution of pure molecular diffusion to mixing. This contribution of molecular diffusion decreases as the Prandtl number increases and we find significant variation in mixing efficiency depending on the dynamics. Therefore, mixing efficiency can not be assumed to be constant, especially for high Prandtl number (Pr=10). In this case (Pr=10) the mixing efficiency for weakly nonlinear waves is 0.09 while the mixing efficiency for wavebreaking is 0.25. This value of 0.25 is in agreement with the 'canonical value' of 0.2-0.3 usually taken for mixing efficiency in the oceanic thermocline. However, further comparisons would require to study the behavior of mixing efficiency versus the Prandtl number in the case of higher Reynolds number simulations.

Mixing results in a weak distortion of the background density profile,madeofasuperposition of less stratified layers and more stratified layers. This structure results from the process of wavebreaking which takes place through unstable shear layers. This structure is very similar to that observed at small scale in the ocean, though the distortion obtained in the numerical simulations is very weak. One reason for this weakness lies in the fact that the wave is unforced. We support this argument by estimating the time that would be required for a full mixing of a layer of a given vertical thickness, for a forced wave; this time is found to be consistent with oceanic timescales. Comparisons with shear layer behavior [START_REF] Staquet | Mixing in a stably-stratified shear layer: two-and three-dimensional numerical experiments[END_REF] suggest that mixing produced by wavebreaking occurs in a way similar to that of a succession of shear layers.

We then show that the Cox number varies like the square of the turbulent Froude number. One interest in this finding is that this law is analogous to mixing models in the Ocean interior in which the eddy diffusion coefficient is related to the finescale shear of the internal wave field [START_REF] Polzin | Statistics of the Richardson number: mixing models and finestructure[END_REF]. Also, this law is identical to an empirical parameterization used in large scale models, in which the Richardson number of the flow comes into play.

Eventually we provide a comparison between two-and three-dimensional calculations of the evolution of the same wave. We show that mixing processes behave very closely. This results from the fact that the breaking of a wave of small amplitude is mostly two-dimensional. The effect of the wave amplitude upon mixing processes in the three-dimensional case is under present investigation.

Fig. 1 .

 1 Fig. 1. Propagating waves (run P1): density fields (first column) and vorticity fields (second column) in the vertical (x, y) plane at different times: (a and b) t=32.6; (c and d) t=36.6; (e and f) t=37.9 and (g and h) t=46.1. Time is expressed in units of the buoyancy period, T. The values of the extrema of the density and the vorticity are indicated on each color scheme.

Fig

  Fig. 1. (Continued).

Fig. 2 .

 2 Fig. 2. Standing waves (run S2): density fields (first column) and vorticity fields (second column) in the vertical (x, y) plane at different times: (a and b) t=49.6; (c and d) t=65; (e and f) t=72 and (g and h) t=95. Time is expressed in units of the buoyancy period, T. The values of the extrema of the density and vorticity are indicated on each color scheme.

Fig

  Fig. 2. (Continued).

Fig. 3 .

 3 Fig. 3. Time evolution of the maximum over space of the kurtosis of the vorticity and of the vertical density gradient: (a and b) for propagating waves (run P1) and (b and d) for standing waves (run S2). The dotted line represents the limit between stable and unstable stratification.

Fig. 4 .

 4 Fig. 4. Time evolution of three lines of particles released in the flow after breaking (run P1). An extended domain of size [-2π,4π] 2 , instead of [0, 2π ] 2 is visualized to give a realistic view of the dispersion. (a) Particles at the time of release, t=36.62; (b) vorticity field at same time (values of the vorticity within -11.73 and 7.99); (c) particles at t=36.62+1.6; (d) vorticity field at same time (values of the vorticity are within -9.60 and 9.77); (e) particles at t=36.62+3.2; (f) vorticity field at same time (values of the vorticity are within -9.10 and 9.12); (g) particles at t=36.62+8 and (h) vorticity field at same time (values of the vorticity are within -7.50 and 8.50). Time is expressed in units of the buoyancy period, T.

Fig. 4 .

 4 Fig. 4. (Continued).

Fig. 5 .

 5 Fig. 5. Vertical profiles of the concentration in particles (i.e. horizontal average of the number of particles per grid cell: c(y, t)), at different times. t 0 denotes the time of release of the particles. The Gaussian fit is represented with a dotted line. (a) c(y, t=t 0 ); (b) c(y, t=t 0 +1.6); (c) c(y, t=t 0 +3.2) and (d) c(y, t=t 0 +8). Time is expressed in units of the buoyancy period, T.

Fig. 6 .

 6 Fig. 6. Vertical dispersion vs. time: (a) time evolution of Z 0 (t) 2 = N i=1 (y i (t)y G (t)) 2 for line 1 (upper line); (c) same for line 2 (middle line); (e) same for line 3 (lower line); (b) time evolution of Z sq (t) 2 = N i=1 (y i (t)-y G (t)) 2 for line 1; (d) same for line 2 and (f) same for line 3. The least square line, calculated during the late evolution, is indicated by a dotted line. The dispersion coefficient, K, is equal to half the slope of this line, its value is given in terms of the molecular diffusivity.

  y top y bottom φ d (y) dy is the volume averaged diapycnal flux. H is the vertical extend of the fluid container. • Φ i =-(gκ/H )[ ρ y=y topρ y=y bottom ] is the rate of increase of E p within the interior of the fluid volume due to pure molecular diffusion only. Φ i will be referred to as the laminar diapycnal flux and the difference Φ d -Φ i as the fluctuating diapycnal flux. Note that, in the present case where a constant density gradient is imposed to the fluid, Φ i retains the same constant value throughout the evolution of the flow. • Φ a = (g/V ) V vρ dV is the volume averaged buoyancy flux and represents the reversible rate of exchange between kinetic and potential energy. • Ėp surf adv and Ėb surf adv are the rate of change of E p and E b due to density advection across the boundaries. These terms vanish for standing waves. • Ėp surf dif and Ėb surf dif are the rate of change of E p and E b due to diffusion at the boundaries.

Fig. 7 .

 7 Fig. 7. Time evolution of the maximum over the numerical domain of the maximum Cox number (black line) and of the spatially averaged Cox number (gray line): (a) standing waves (run S2) and (b) propagating waves (run P1).

  ), the change in the background potential energy is controlled by surface fluxes, Ėb surf dif and Ėb surf adv ,a n d by diffusive processes within the fluid volume associated with the diapycnal flux, Φ d .For standing waves where the surface advective flux vanishes, the sum of the time-integrated diffusive surface flux and diapycnal flux describe exactly the evolution of the background

Fig. 8 .

 8 Fig. 8. The budget of the background potential energy. (a) Time evolution of the change in the background potential energy E b (t)-E b (0), relative to the total change in the primary wave energy during the flow evolution E wave : (E b (t)-E b (t=0))/ E wave for standing waves (run S2); (b) same as part (a) (gray line) but for propagating waves (run P1). Also plotted is (Φ d + E b surf dif ) dt/ E wave (black line) and (c) energy arising from mixing within the interior: E m (t) = t t=0 (Φ d -Φ i ) dt/ E wave (run S2).

Fig. 9 .

 9 Fig. 9. The budget of the available potential energy. (a) Run S2. Available potential energy, E a (black line); the gray line is the opposite of the time integrated fluctuating diapycnal flux Φ d -Φ i added to E a (0) and (b) same but for propagating waves (run P1); the time integrated buoyancy flux t 0 (Φ a ) dτ is also represented (thin black line).

Fig. 10 .

 10 Fig.10. Mixing efficiency (defined by Eq. (22)) during the turbulent regime for runs P1 (Pr=0.72), P6 (Pr=1), P7 (Pr=5)andP8(Pr=10). The dependency for a linear wave, 1/(1+Pr) is also plotted (dashed line).

Fig. 11 .

 11 Fig. 11. Volume averaged Cox number vs. the turbulent Froude number (defined by Eq. (23)) for all runs of Table 1 (except runs P6-P10, at all discrete times). The dotted line is the straight line: Cox∼1.12+2.78(Fr t -Fr wave ) 2 , where Fr wave =Fr t (0).

  (for run P1) the relative change in the background density profile δρ b (y, t) = [ρ b (y, t)ρ b (y, t 0 )/ρ b (y, t)]; the change in the background density profile is weak, as we shall see, and this is why this relative change is displayed. The vertical profile of the buoyancy frequency, N 2 b =-(g/ρ 0 )(dρ b /dy) and of the Cox number are also presented on the figures. Two times are shown: the time at which the Cox number reaches an absolute maximum value and the end of the calculation.

Fig. 12 .

 12 Fig. 12. Standing waves (run S2): (a) vertical profile of the relative change in the background density profile, [ρ b (y, t)-ρ b (y, t=0)]/ρ b (y, t=0) (in %) at the time the Cox number reaches its absolute maximum (t=110); (b) buoyancy frequency profiles, N 2 b (y, t) at t=110; (c) Cox number profile at t=110; (d) same as part (a) but at the end of the calculation (t=134); (e) same as part (b) but at t=134 and (f) same as part (c) but at t=134. Time is expressed in terms of the buoyancy period, T.

Fig. 13 .

 13 Fig. 13. Propagating waves (run P1): (a) vertical profile of the relative change in the background density profile, [ρ b (y, t)-ρ b (y, t=0)]/ρ b (y, t=0) (in %) at the time the Cox number reaches its absolute maximum (t=42); (b) buoyancy frequency profiles, N 2 b (y, t) at t=42; (c) Cox number profile at t=42; (d) same as part (a) but at the end of the calculation (t=64); (e) same as part (b) but at t=64 and (f) same as part (c) but at t=64. Time is expressed in terms of the buoyancy period, T.

Fig. 14 .

 14 Fig. 14. Comparison between a two-and a three-dimensional simulation of a propagating breaking gravity wave of low amplitude (runs P1 and P10 in Table 1). (a) Total energy vs. time and (b) time evolution of the change in the background potential energy E b (t)-E b (0), relative to the total change in the primary wave energy during the flow evolution E wave :(E b (t)-E b (t=0))/ E wave .

Table 1

 1

Table of

 of

		runs a
	Run	A	dt	ν	Pr	Resolution	Fr=U/(Nλ)	Re=U/(Nλ)
	Propagating waves
	P1 P2 P3 P4 P5 P6 P7 P8 P9 P10	0.256 0.384 0.5 0.62 0.75 0.256 0.256 0.256 0.256 0.256	0.004 0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.008 0.012	1.5×10 -4 1.5×10 -4 3×10 -4 3×10 -4 3×10 -4 1.5×10 -4 1.5×10 -4 1.5×10 -4 3×10 -4 3×10 -4	0.72 0.72 0.72 0.72 0.72 1 5 10 0.72 0.72	512 2 512 2 512 2 512 2 512 2 512 2 512 2 512 2 256 2 256 3	4.07×10 -2 6.11×10 -2 7.96×10 -2 9.87×10 -2 1.19×10 -1 4.07×10 -2 4.07×10 -2 4.07×10 -2 4.07×10 -2 4.07×10 -2
	Standing waves
	S1 S2 S3	0.192 0.256 0.4	0.001 0.001 0.001	1.25×10 -4 1.5×10 -4 1×10 -3	1 1 1	513 2 513 2 513 2	3.05×10 -2 4.08×10 -2 6.37×10 -2

Acknowledgements

This work has been initiated during a visit of KBW at ENS Lyon, thanks to a CNRS Associate Research position. We thank A. Babiano for having provided subroutines for Lagrangian particles trajectories. We acknowledge as well the two anonymous reviewers for useful comments which helped to improve the manuscript. Calculations have been performed on the Cray C98 of IDRIS (CNRS computer centre) under contracts no. 970580 and 980580.