Mathilde ; Cavelier

B.-J Gréa

A Briard

L Gostiaux

The subcritical transition to turbulence of Faraday waves in miscible fluids

Keywords: Parametric instability < Instability, Transition to turbulence < Instability, Faraday waves < Waves/Free-surface Flows Faraday instability, mode selection, wavebreaking

We study the development and the breaking process of standing waves at the interface between two miscible fluids of small density contrast. In our experiment, a subharmonic wave is generated by a time-periodic vertical acceleration via the Faraday instability. It is shown that its wavelength may be selected not only by the linear process predicted by the Floquet theory and favoring the most unstable modes allowed by the tank geometry, but also by a nonlinear mode competition mechanism giving the preference to subcritical modes. Subsequently, as the standing wave amplitude grows, a secondary destabilization process occurs at smaller scales and produces turbulent mixing at the nodes. We explain this phenomenon as a subcritical parametric resonance instability. Different approaches derived from local and global stability analysis are proposed to predict the critical wave steepness. These theories are then assessed against various numerical and experimental data varying the frequencies and the amplitudes of the forcing acceleration.

Introduction

The triggering of waves at the interfaces between fluids of different densities by vertical vibrations is a well known phenomenon first observed by [START_REF] Faraday | On the forms and states of fluids on vibrating elastic surfaces[END_REF] and extensively reviewed in [START_REF] Miles | Parametrically forced surface waves[END_REF]. It constitutes a classical example of parametric instability in fluids and it has greatly helped the understanding of pattern formation in nonlinear systems [START_REF] Edwards | Patterns and quasi-patterns in the faraday experiment[END_REF][START_REF] Kudrolli | Patterns and spatiotemporal chaos in parametrically forced surface waves: a systematic survey at large aspect ratio[END_REF][START_REF] Godrèche | Hydrodynamics and nonlinear instabilities[END_REF][START_REF] Kahouadji | Numerical simulation of supersquare patterns in faraday waves[END_REF]. In this context, a considerable number of studies successfully characterized the instability onset using linear Floquet theory. We mention here only the major contributions of [START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF] and [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF]. In addition, the development of new weakly-nonlinear approaches was of paramount importance to predict the saturation amplitudes of the Faraday waves [START_REF] Douady | Experimental study of the faraday instability[END_REF][START_REF] Zhang | Pattern formation in weakly damped parametric surface waves driven by two frequency components[END_REF][START_REF] Chen | Amplitude equation and pattern selection in faraday waves[END_REF][START_REF] Skeldon | Can weakly nonlinear theory explain faraday wave patterns near onset?[END_REF], to disentangle the multi-modal interactions leading to spatiotemporal chaos [START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF][START_REF] Meron | Low-dimensional chaos in surface waves: Theoretical analysis of an experiment[END_REF][START_REF] Gollub | Spatiotemporal Chaos in Interfacial Waves[END_REF] or to evidence the bifurcations and hysteresis phenomena [START_REF] Rajchenbach | Faraday waves: their dispersion relation, nature of bifurcation and wavenumber selection revisited[END_REF][START_REF] Périnet | Hysteretic faraday waves[END_REF].

By contrast, the Faraday instability in the turbulent regime has been less studied. On the one hand, this is possibly due to a lack of theoretical tools as strong nonlinearities still † Email address for correspondence: benoit-joseph.grea@cea.fr (b): Time series images from the camera zooming on one wavelength and presenting two oscillation periods of the primary wave. This illustrates the different stages of the wavebreaking with first a 'blurring' of the interface at the node followed by a 'roll-up'. This case corresponds to the b5 experiment which parameters are detailed in table 1. (c): Visualisation of the interface at wavebreaking in the direct numerical simulation DNSd3 (the parameters are given in table 2). The reference frame as well as the acceleration direction are also indicated.

lie 'almost entirely outside the realm of available analytical techniques', as commented by [START_REF] Miles | Parametrically forced surface waves[END_REF]. On the other hand, there are few Faraday experiments dedicated to the subject as most of them are conducted in small apparatus with high viscosity fluids to better control the dissipation process [START_REF] Bechhoefer | An experimental study of the onset of parametrically pumped surface waves in viscous fluids[END_REF]. Besides, the turbulent regime has been investigated, in particular for miscible fluids with small density contrast (see [START_REF] Zoueshtiagh | Experimental and numerical study of miscible faraday instability[END_REF][START_REF] Amiroudine | Mixing generated by faraday instability between miscible liquids[END_REF] where it is observed that the turbulent mixing layer driven by vertical vibrations grows and eventually saturates. This indeed occurs as the natural frequencies of the system decrease with the enlargement of the layer and are no longer parametrically excited by the periodic forcing. By retaining only the nonlinear interactions of turbulence with the mean flow [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF], the final size of the turbulent mixing layers can be predicted analytically [START_REF] Gréa | What is the final size of turbulent mixing zones driven by the faraday instability[END_REF]. This prediction has been recently confirmed experimentally in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF].

Concerning more specifically the transition to turbulence in the Faraday problem, it is known since [START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF] and [START_REF] Meron | Parametric excitation of multimode dissipative systems[END_REF]

that chaotic behaviours often

The subcritical transition to turbulence 3 appear for parameters in the vicinity of the neutral branch intersections of the stability diagram. The experiments presented in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] postulate several scenarios of transition to turbulence. For instance, due to the large dimensions of the tank allowing viscous effects to be negligible, harmonic and subharmonic modes can interact to generate mixing at small scales as also reported in numerical simulations [START_REF] Briard | Harmonic to subharmonic transition of the faraday instability in miscible fluids[END_REF]. However this experimental campaign also evidences that turbulence can result from the breaking process of a single Faraday mode. This phenomenon is illustrated in figure 1 (see also online movies) showing a growing subharmonic primary Faraday wave subjected to a destabilization process occurring at the nodes and rapidly producing turbulent mixing. The objective of this work is to investigate and explain this mechanism.

The breaking of Faraday waves at free surfaces is known to appear at the wave crest and leads to the formation of jets (as shown by [START_REF] Jiang | Period tripling and energy dissipation of breaking standing waves[END_REF][START_REF] Wright | Numerical studies of two-dimensional faraday oscillations of inviscid fluids[END_REF][START_REF] Longuet-Higgins | Vertical jets from standing waves; the bazooka effect[END_REF][START_REF] Kalinichenko | Breaking of faraday waves and jet launch formation[END_REF]. This comes from the modulation of the primary wave interacting with its first temporal harmonics. More generally, the crests of the waves in the ocean are also subject to destabilization. This topic due to its importance is well documented, shedding light on the many instability mechanisms which can develop [START_REF] Banner | Wave breaking in deep water[END_REF][START_REF] Kiger | Air-entrainement mechanisms in plunging jets and breaking waves[END_REF]). Yet, the breaking process of free surface waves differs sensitively from the observations reported in figure 1. By contrast, our problem presents very close similarities with the destabilization of standing waves described by [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF], also in the context of miscible fluids with small density variations. While the primary wave in Thorpe's experiment is not generated by vertical vibrations but by lateral plungers, a vortex is still produced at the wave node as in figure 1. More precisely, several phases can be identified for the instability with a 'blurring' of the interface preceding its 'roll-up'. [START_REF] Kalinichenko | Developement of a shear instability in nodal zones of a standing internal wave[END_REF] has also observed and investigated experimentally the breaking process of a Faraday wave between miscible or immiscible fluids. In particular, he reported that the secondary instability starts for wave steepness ka ∼ 0.4, with k the wavenumber and a the amplitude of the Faraday wave. The Rayleigh-Taylor type instability does not seem to play a role in the process as the acceleration induced by the primary wave displacement is not sufficient to invert the gravity. In place, both [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] and [START_REF] Kalinichenko | Developement of a shear instability in nodal zones of a standing internal wave[END_REF] convenes that a sort of Kelvin-Helmholtz instability, 'although not in a simple form', is at work. Due to the strong time dependence of this configuration, evaluating locally the Richardson number at the node cannot be sufficient to assess the importance of the shear instability. Additionally, in the context of internal gravity waves, the role of subharmonic secondary parametric instabilities in the breaking process has been explored [START_REF] Mcewan | Parametric instability of internal gravity waves[END_REF][START_REF] Bouruet-Aubertot | Breaking of standing internal gravity waves through two-dimensional instabilities[END_REF][START_REF] Benielli | Excitation and breaking of internal gravity waves by parametric instability[END_REF][START_REF] Staquet | Internal gravity waves: From instabilities to turbulence[END_REF][START_REF] Sutherland | Internal gravity waves[END_REF][START_REF] Yalim | Parametrically forced stably stratified flow in a three-dimensional rectangular container[END_REF]. Can this mechanism also apply to Faraday waves?

This paper is organized as follows: we give in section 2 a brief description of the experiments and numerical simulations used for this study. In section 3, we analyze the characteristics of the primary Faraday wave, emphasizing in particular the mode selection mechanism. The section 4 is dedicated to the wavebreaking process with two theoretical approaches proposed and shedding light on the importance of a subharmonic secondary instability. We then detail our methodology in order to measure the wavebreaking amplitudes in section 5. Finally, the analysis and discussion of the results in view of the theoretical predictions are provided in section 5.3.

Generalities

This work, dedicated to the wavebreaking of Faraday waves, relies on several experiments already presented in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] and initially designed to study the turbu-4 Cavelier et al.

lent mixing driven by vertical vibrations. We first detail the configuration used and the parameters considered. Next, we present the direct numerical simulations which allow us to explore an even broader range of parameters and to identify how the transition to turbulence takes place.

Experimental set-up and parameters

The experimental set-up is now briefly introduced since the details can be found in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF]. We fill a cuboidal tank of inner length W = 94.6 cm, width D = 11 cm, and height H = 67 cm, with salt and fresh water (see figure 1). The salt water density takes the values ρ 1 = 1030, 1060 or 1090 kg m -3 while, for the fresh water, we get ρ 2 = 998 kg m -3 . This corresponds to various Atwood numbers, expressing the density contrast, A = (ρ 1ρ 2 )/(ρ 1 + ρ 2 ) ∈ {0.015; 0.03; 0.045}. The heavier salty water layer is initially placed at the bottom. It is separated from the lighter fresh water by a thin diffuse interface of thickness δ = 0.5 -1.5 cm located at half the height of the tank. This thickness may vary due to the filling procedure of the tank. The values of δ can be measured either by the initial image from the camera or by the vertical density profiles obtained from a probe before the experiment starts. An hexapod oscillates the tank along the z-vertical direction (for the horizontal directions, x corresponds to the length W and y is along the width D of the tank). This generates a well-controlled time dependent vertical acceleration of intensity G(t) = G 0 (1 + F cos ωt). Here G 0 = 9.81m s -2 is the usual gravitational acceleration, ω the frequency, and F the forcing parameter. This forcing parameter is related to the vertical displacement amplitude of the hexapod, a h , as F = a h ω 2 /G 0 . In the experiments, the acceleration does not change sign since F < 1, although the displacement amplitude of the vessel can be as large as a h = 45 cm. Table 2: Label (series and number) and parameters in physical units (Atwood number, forcing parameter and frequency) taken for the direct numerical simulations presented in this work. The cases DNSa, d and e correspond to the wavebreaking detection. The series DNSb and c are dedicated to the competition between mode (2, 0) and (3, 0) where the selected mode appears underlined. The parameter r expresses the initial amplitude ratio (r = 0 corresponding to a pure (2, 0) mode). The initial amplitude of the interface perturbation, the y-spanwise perturbation amplitude 1 for DNSf together with the interface thicknesses δ are also detailed. The computation domain is of cubic size with length W = 94.6 cm or 2W for the series labelled with * . All the DNSs have a 1024 3 grid resolution.

Series Number

A F ω [rad s -1 ] km,n [m -1 ] Mode (m, n) δ [cm]

Series Number

A F ω [rad s -1 ] Mode (m, n) r / 1 [cm] δ = 3σ [cm]
We select in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] the experiments with sharp initial interfaces and developing a single Faraday wave. The cases exhibiting different modes appearing simultaneously are not considered. Therefore, our study is based on 18 experiments shown in table 1, and grouped by values of A.

The primary standing Faraday waves observed in these experiments are characterized by a horizontal wavenumber, k m,n = k 2

x + k 2 y , associated with the mode index in the x and y directions, respectively m = k x W/π, and n = k y D/π. So for instance m = 2 corresponds to a wavelength equalling the width W of the tank. As the primary wave is "two-dimensional", this implies a zero mode index n = 0. It can be seen in table 1 that configurations with odd and even modes have been investigated.

Direct numerical simulations

In this work, we also provide direct numerical simulations (DNS) in order to explore a broader panel of parameters and to investigate the inner mechanisms of wavebreaking.

These simulations solve numerically the Navier-Stokes equations under the Boussinesq approximation. They express the dynamics of the incompressible fluid velocity U(x, t) and the concentration of heavy fluid C(x, t). Here for miscible fluids with small Atwood number, the dimensionless concentration C(x, t) ∈ [0 1] is related to the density as ρ(x, t) = ρ 2 + (ρ 1ρ 2 )C(x, t). In the reference frame attached to the container, this leads to the classical system of equations

∂ t U + U • ∇U = -∇Π -2AG(t)Ce z + ν∇ 2 U, (2.1a)
∂ t C + U • ∇C = D∇ 2 C, (2.1b) ∇ • U = 0.
(2.1c)

In (2.1a-2.1c), Π refers to a reduced pressure and ν, D are the kinematic viscosity and molecular diffusion coefficients respectively. This set of equations constitutes our theoretical framework in order to predict the wavebreaking. It also describes reasonably well the flow dynamics at large scale in the experiments despite the variations of the viscosity and diffusion coefficients, < 20%, between fresh and salt water.

The simulations are performed in a triply periodic cubic box of size W (or 2W ) using the code already described in [START_REF] Briard | Harmonic to subharmonic transition of the faraday instability in miscible fluids[END_REF][START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF]. Therefore, we do not seek at reproducing the tank's walls which do not play a direct role in the wavebreaking phenomenology (although the walls can play a decisive role in the final transition to turbulence). The code is based on a pseudo-spectral collocation method with two-third rule dealiasing. The time advancement is realized through a third-order low-storage strongstability-preserving Runge-Kutta scheme, with implicit viscous terms. All the simulations use a 1024 3 grid box with a slab-decomposition on 1024 cores. Due to the vertical periodicity, a thin penalization layer is applied to freeze the velocity and concentration fields at the top and bottom of the computational domain. This method is extensively described in appendix B of [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF]. Several tests varying the width of the penalization band have been conducted in order to ensure that this vertical treatment has no impact on the dynamics of the interface. In all the simulations presented, the amplitude of the Faraday wave is less than half of the vertical height of the non-penalized domain.

The simulations parameters are presented in table 2. In order to have a well resolved flow field, the viscosity and diffusion are fixed at ν = D = 2.26 × 10 -6 m 2 s -1 . This corresponds roughly to twice the real viscosity of water but largely overestimates the molecular diffusion (the Schmidt number is 1 instead of 700 for salt water). Due to the dimensions of the tank, this limitation does not prevent to properly capture at least the first stages of secondary instabilities developing on the Faraday wave.

We now detail the initial conditions taken in the simulations. While the initial velocity is U = 0 in the simulations, the initial concentration profile is taken two-dimensional (2D) of the form (for the DNSa, d and e series and outside the penalization band)

C(x, z) = 1 2 1 + tanh z -ξ(x) σ , with ξ(x) = sin(k m,0 x).
(2.2)

The parameter σ in (2.2) sets the initial width δ of the interface (δ ≈ 3σ). The function ξ(x) indicates the initial perturbed interface position of sinusoidal shape with wavelength k and of small amplitude = 1.5 cm. Therefore, the initial interface is slightly more diffused and has a larger amplitude in the simulations compared to the experiments. This is to have at least 10 grid points across the interface layer and to ensure grid convergence of the simulations. Without ambiguity, the initial wavenumber k m,0 for DNSa, d, and e also corresponds to the observed wavenumber at later times indicated in table 2 and characterizing the subharmonic Faraday wave. This is due to our choice for the forcing frequency taken nearly as twice the value of the dispersion relationship of an inviscid interface ω = 2 AG 0 k m,0 .

In order to explore more broadly the effect of mode selection, we also propose simulations DNSb and c with an initial interface position defined as

ξ(x) = [r cos(k 3,0 x) + (1 -r) cos(k 2,0 x)].
(2.3)

Here, the parameter r thus expresses the initial ratio amplitude between the modes (3, 0) and (2, 0). These simulations are conducted in a computational domain twice the size of the tank, 2W , in order to allow the development of odd modes otherwise forbidden due to the periodic boundary conditions. Also, the interface thicknesses δ are doubled to keep at least 10 grid points across the interface while the viscosity and diffusion coefficients are multiplied by 4 still to ensure grid convergence of these simulations.

In the simulation series DNSa-e, the flow remains two dimensional even after the secondary instability starts. In order to study the full transition to turbulence, we consider simulations DNSf where the interface position is slightly perturbed in the spanwise direction y. Introducing the normalized white noise function f , the interface position is given by ξ(x, y) = sin(k m,0 x) + 1 f (y).

(2.4) In practice, the y disturbance amplitude is set such that 1 = 10 -2 . However, we have also tested various simulations varying the 1 parameter, not presented as exhibiting the same phenomenology as DNSf. The breaking of the spanwise symmetry invariance in DNSf can also be produced by the lateral boundary layers in the experiments. Therefore, various simulations mimicking the lateral boundary layers were also conducted using the penalization method introduced in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF]. These simulations (not presented) give also similar results as DNSf which will be discussed in section 5.3.3.

Mode selection mechanism of the Faraday wave

In this part, we discuss the primary wave characteristics and try to figure out which linear or nonlinear mechanism eventually selects the dominant wavelength of the instability in the experiments.

Linear theory

It is well-known since [START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF] that, when modelling the Faraday instability, the amplitudes, η k , for the interface modes of wavenumber k are ruled by a Cavelier et al.

Mathieu equation ηk + 2γ(k) ηk + Ω 2 (k)(1 + F cos ωt)η k = 0.
(3.1)

In (3.1), we define the inviscid frequency of the diffuse interface Ω and the viscous damping term γ, both of which depend on the horizontal wavenumber k. Note that the decoupling of each inviscid mode is only true in the limit of small damping. However, the full analysis of this problem can be performed using the method proposed by [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF].

The inviscid frequency Ω(k) within the deep water approximation is thus a growing function of k and can be evaluated for a given vertical density profile, see for instance [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] for piecewise linear profile

Ω(k) = AG 0 k 1 + kδ/2 1/2 . (3.2)
For small wavenumbers, i.e. kδ 1 with δ the thickness of the interface, the classical dispersion relationship for an interface within the deep water approximation Ω = √ AG 0 k is recovered. In the large wavenumber limit, kδ 1, the interface mode reduces to Ω(k) = 2AG 0 /δ corresponding to the local buoyancy or Brunt-Väisälä frequency at the interface.

Mode (2, 0) (3, 0) (4, 0) (5, 0) (6, 0) (8, 0) γ w = 9.75 × 10 -3 1.08 × 10 -2 1.15 × 10 -2 1.22 × 10 -2 1.28 × 10 -2 1.36 × 10 -2 γ b = 1.16 × 10 -4 2.60 × 10 -4 4.62 × 10 -4 7.22 × 10 -4 1.04 × 10 -3 1.85 × 10 -3 γ δ = 1.74 × 10 -3 2.61 × 10 -3 3.48 × 10 -3 4.35 × 10 -3 5.22 × 10 -3 7.27 × 10 -3

Table 3: Values for the damping coefficients, γ w , γ b , γ δ , in s -1 and evaluated for the largest wavelengths developing in the experiment. We assume here that the Atwood number is A = 0.03 and the thickness of the interfacial layer is δ = 0.5 cm. Here, the top boundary is taken as a wall to evaluate γ h (the values would be nearly the same for a free surface).

The viscous dissipation term, γ(k), expressing the small interfacial mode damping, can have different origins. The damping coming from the bulk flow for a sharp interface takes the form γ b (k) = 2νk 2 (Lamb 1945;Landau & Lifshitz 2013). However, due to the velocity gradients, a significant damping can also occur within the thin layer separating the two fluids. Assuming a piecewise linear vertical density profile, [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] obtained the expression γ δ (k) = AG 0 νk 2 /Ω 2 δ ≈ νk/δ for kδ 1. In this linear theory, we wish also to account for the damping generated by the boundary layers at the various walls (top, bottom and laterals) existing in the experiments. The boundary layer widths in the experiment can be evaluated using δ w = (2ν/Ω) 1/2 . This gives values around δ w ∼ 1 -2 mm using the parameters of the experiments, showing that the boundary layer widths are much smaller than the characteristics wavelengths of the instability and the size of the tank. In this condition, [START_REF] Keulegan | Energy dissipation in standing waves in rectangular basins[END_REF]; [START_REF] Miles | Surface-wave damping in closed basins[END_REF] have derived an expression for the damping of free surface waves in a rectangular basin due to the laminar boundary layers. This result has also been generalized to our problem by [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] as detailed in appendix B and leads to the following expression for the damping coefficient The neutral curves (thick plain lines) have a slight dependence on the Atwood number explaining that they are not completely superimposed. The symbols correspond to the parameters taken in the experiment in table 1. The shapes indicate the Atwood number and the colors reveal which mode is eventually selected.

γ w ≈ ν Dδ w = √ νΩ √ 2D . (3.3) 1.5 2.0 2.5 3.0 3.5 4.0 ω/ AG 0 k 2,0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 F (2,0) (3,0) (4,0) (5,0) (6,0) (7,0) (8,0) (9,0) A = 0.015 A = 0.030 A = 0.045 A = 0.015 A = 0.030 A = 0.045
Here, (3.3) thus expresses the dominant contribution of the lateral walls (in the plane z-x) to the damping. We gather in table 3 the numerical values of the damping coefficients originating from the bulk, the interfacial layer separating the fluids and from the boundary layers at the wall. It can be shown that these values do not vary more than several percents if we account for the viscosity contrast between fresh and salt water. Therefore, it clearly indicates that the dissipation occurs essentially in the viscous layers at the walls as γ w is larger than the other contributions. The damping γ w scales indeed like ν 1/2 in (3.3) while it is linear in ν for the dissipation γ b or from the interfacial layer γ δ . [START_REF] Bechhoefer | An experimental study of the onset of parametrically pumped surface waves in viscous fluids[END_REF] have extensively discussed this aspect and they suggest to use fluids with high viscosity in order to better control the dissipation in experiments dedicated to the study of the instability threshold. By contrast, our study is focused on the wavebreaking mechanism explaining why we favor the use of low viscosity fluids. Note that for larger wavenumber k, the contributions from the interface layer regain in importance and cannot be neglected.

The stability diagram corresponding to the first subharmonic tongue is represented in figure 2. It is plotted for the different large scale modes of the tank and derived using the damping γ w from (3.3) and γ δ (the latter contribution being smaller). The neutral curves of (3.1) are computed using the method proposed by [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] and also used in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] assuming different Atwood number values A and an initial interface thickness δ = 1 cm. For a given mode k, the minimum forcing F th able to destabilize the interface occurs at the frequency corresponding to the first subharmonic resonance, Ω(k) = ω/2. The classical asymptotic theory of the Mathieu equation, in the limit of small damping, allows the derivation of the threshold as F th = 8γ/ω (see [START_REF] Rajchenbach | Faraday waves: their dispersion relation, nature of bifurcation and wavenumber selection revisited[END_REF] for instance). The threshold F th varies very weakly for the different modes presented in figure 2. Indeed the contribution due to the damping from the viscous layer at the walls scales like γ w ∼ ω 1/2 leading to a decrease of F th at larger ω. However this effect is compensated at larger k by the contribution from the damping at the interface scaling like γ δ ∼ ω 2 .

The parameters taken in the experiments with F 0.3, also indicated in figure 2, are situated in unstable regions well above the viscous thresholds determined by the linear Floquet theory. As a consequence, at least two or more modes can be simultaneously subharmonically unstable in these experiments.

Linear or nonlinear mode selection?

In this section, we investigate the mechanisms leading to the mode selection of the primary wave. As shown in figure 2 and due to the large acceleration forcing F , several modes can be linearly unstable and play a role in the interface dynamics. Surprisingly, a single mode, corresponding nearly always to the smallest unstable wavelength, emerges from this process; there is a clear tendency to favor the modes pertaining to the right unstable tongues in figure 2 (the mode reported in table 1 is also indicated by the color of the symbol in figure 2). In addition, the selection mechanism does not apparently discriminate the even or the odd modes of the tank as both can be observed in the experiments.

One would expect the modes with the largest linear growth rates to be selected first: this is why only the modes in the first subharmonic band are considered here, as the higher resonance regions exhibit much smaller amplification rates. For a given mode, the Floquet theory shows that the maximum amplification occurs for parameters close to the subharmonic resonance frequency located at the center of the instability tongue. However, the results in figure 2 reveal that in many cases, the selected mode does not have the largest linear growth rate. Even more, some of the observed modes are hardly unstable and should have very small growth rate from linear theory (such as EXPa1 for instance). This statement stands even if we account for some experimental uncertainties in term of Atwood number (±0.001) or initial interface width (±0.5 cm). It can be shown that these effects only slightly modify the instability tongues of figure 2. In particular, a larger interface thickness would left-shift a bit the instability tongues of figure 2 as the natural frequencies Ω are decreased (the damping dominated by the viscous layer at the wall remains unchanged). Besides, we have checked that the hexapod movement is well controlled and remains sinusoidal. Therefore, it is unlikely to have spurious forcing frequencies in the system which may change the linear stability of the problem.

The initial perturbation of the interface may also play a role in the mode selection mechanism. A large initial amplitude on a given mode can explain why it appears even if it has not the largest growth rate during the linear phase. This would suggest that an initial condition at small scales is at work in the experiments although we have not observed such disturbance or could not identify a source able to generate it. Besides, this cannot shed light on the appearance of linearly stable modes.

By contrast, Faraday experiments with immiscible fluids have revealed the ability of nonlinearities to select modes and generate transient chaotic regimes (see for instance [START_REF] Ciliberto | Pattern competition leads to chaos[END_REF][START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF]. Using weakly nonlinear approaches, Meron & Procaccia The subcritical transition to turbulence 11 (1986) and [START_REF] Meron | Parametric excitation of multimode dissipative systems[END_REF] have already detailed how the mode suppression phenomenon can occur. Considering two modes close to the first subharmonic resonance, the nonlinear cubic coupling terms in each amplitude equation have a sign determined by the detuning parameter ∆ = Ω 2 ω 2 -1 4 . Here, the detuning parameter for a given mode expresses the departure of the forcing frequency from the first subharmonic resonance. Therefore, when two modes are in competition, their respective detuning parameters ∆ generally have an opposite sign because the forcing frequency ω lies between the two subharmonic resonance frequencies (see figure 2). This explains the mode suppression since one mode can develop, even being linearly stable, by pumping energy from the other one. The theory indeed shows that the vanishing mode, ∆ < 0, is supercritical as the nonlinear coupling damps the instability while the dominant one, ∆ 0, is subcritical as it is being reinforced by the nonlinearities. As the wave amplitude grows, the frequency of the wave tends to diminish [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF]) as well as the detuning parameter of each modes [START_REF] Godrèche | Hydrodynamics and nonlinear instabilities[END_REF]. Hence, this left-shifts the instability bands of figure 2 and favors the subcritical modes at smaller wavelength.

Numerical analysis of the mode competition

At this stage, the mode suppression due to a nonlinear coupling between competing modes can explain the mode selection evidenced in figure 2. Besides, the mode amplitude is no longer negligible compared to its wavelength when the selection process is at work, suggesting that the nonlinear effects are an important ingredient to account for. We wish to assess further this hypothesis by the mean of direct numerical simulations (DNS) with well-characterized initial conditions. Two series of DNS have been performed using 1024 3 grid points (Series DNSb and c in Tab. 2), with A = 0.03. The frequency ω and the forcing F taken in the simulations are also represented in the phase diagram of figure 3. It is important to stress here that the phase diagram does not account for the wall damping as simulations are performed in a triply periodic box. The two series of DNS start from the same location in the phase diagram (point A) with F = 0.8 and ω = 3.07 (or equivalently ω/ AG 0 k 2,0 = 2.2). This corresponds to parameters with the two unstable modes having nearly the same exponential growth rate as ∼ e λωt . Indeed, the Floquet exponent λ takes the value 0.09 for mode (3, 0) and 0.07 for mode (2, 0). We fix the forcing parameters in one series and vary the amplitude ratio r. In the other series, the relative amplitude is set at r = 0.5 and we decrease the forcing frequency ω in order to explore more deeply the (2, 0) subharmonic instability tongue. The simulations are stopped when the wavebreaking occurs and we report which of the (2, 0) or (3, 0) mode prevails at this moment in figure 3. This procedure is performed both visually and by computing the Fourier modes of the interface.

The simulations clearly evidence the mode suppression phenomenon to the benefit of the modes with the smallest wavelength. The results reported in figure 3 show the dominance of mode (3, 0) even starting from a small initial amplitude (transition occurs at r = 0.1) or developing in a region where it is linearly stable (for small ω). The phenomenon can be scrutinized in more detail on the snapshots extracted from the two series in the figure 4 where the mode (3, 0) emerges from cases with initial r = 0.25, 0.5 or with the frequencies ω = 2.6, 2.8. Indeed, in the last row of figure 4, one can observe that mode (2, 0) prevails only for ω = 2.4 (ω/ AG 0 k 2,0 = 1.72) and for larger ω, say ω = 2.6 (or ω/ AG 0 k 2,0 = 1.86) mode (3, 0) is visible despite being linearly stable. As a consequence, the mode competition greatly enhances the sensitivity to initial conditions in the experiment. As importantly, this process breaks the symmetry of the primary wave as can be seen in both the experiments (figure 1) and the simulations (figure 4). Some specificities of the mode suppression in our Faraday experiment with miscible 12 Cavelier et al. accounting for viscosity and diffused interface δ = 1.8 cm (as no walls are present in the DNS, the damping coefficient γ ≈ γ δ and critical thresholds F th are very small). The symbols' color indicates which mode emerges from the simulation, the mixed colors express that both modes can be observed. Two series of DNSb, c (see table 2) are presented here starting from point A. In the DNSb group, the initial amplitude ratio r between modes (3, 0) and (2, 0) is set at r = 0.5 and we decrease the forcing frequency ω. In the DNSc group, the frequency and forcing are fixed and we vary r. The point corresponding to DNSe is also placed.

1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 ω/ AG 0 k 2,
fluids are now addressed. We have not observed oscillations between two specific modes as in [START_REF] Ciliberto | Pattern competition leads to chaos[END_REF][START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF] or similarly [START_REF] Yalim | Parametrically forced stably stratified cavity flow: complicated nonlinear dynamics near the onset of instability[END_REF] in the context of a stable stratification. This is notably because, at large forcing parameter, the interface irreversibly grows allowing continuously new modes to be destabilized. The modes can change in our experiment as already reported in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF]. But it always corresponds to a one way transition from large to small wavelength for interface modes.

The more complex transitions evidenced in figure 14 of [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF], for instance, refer to modes pertaining not to the same instability band or being of different nature.

The irreversible mixing produced by the rapid breaking of the primary waves also explains this aspect.

Another difference with past Faraday immiscible experiments conducted in a shallow basin is that in our case the dominant waves correspond to those with the smallest wavelength as already discussed. Noticeably this point clearly agrees with the nonlinear theory of mode suppression. It can be shown that within the deep water approxima- 2) with ω varying and A = 0.030, F = 0.8, r = 0.5, = 1.5 cm. We put the slices of width W of the concentration field at four instants starting from the initial interface at t = 0 and ending when the wavebreaking occurs; pure fluids remain in white while mixed fluid with C = 0.5 ± 0.15 is in black. tion, the subcritical modes are indeed those with small wavelength, see [START_REF] Rajchenbach | Faraday waves: their dispersion relation, nature of bifurcation and wavenumber selection revisited[END_REF] for details.

In this part, it has been evidenced that the mode selection of the primary wave may result from a complex nonlinear mode competition process. When this is the case, the subcritical mode is eventually selected. In the following, we try to explain the breakdown of the Faraday waves.

Modelling the breakdown of Faraday waves

We now present two heuristic models dedicated to the breakdown of Faraday waves initiating the transition to turbulent mixing. The objective is to evaluate the critical wave steepness at which the breakdown may occur. By emphasizing two simple approaches, we aim at exploring various frameworks for the breakdown and at disentangling the inner mechanisms responsible for the instability.

Both models, although relying on different assumptions, suggest that the breakdown results from a subharmonic secondary instability at small scales. Therefore, one key ingredient in these approaches is to account for the unsteadiness of the primary wave. This aspect differs from secondary instability analysis relying on a frozen base flow used for instance in the context of Kelvin-Helmholtz instability [START_REF] Salehipour | Turbulent diapycnal mixing in stratified shear flows: the influence of prandtl number on mixing efficiency and transition at high reynolds number[END_REF].

The first approach, hereafter referred as global, is based on the fact that the breakdown of the Faraday waves changes the monotony of horizontally averaged density profiles (this point is more specifically detailed later in section 5.2). We therefore seek at identifying Figure 5: Frameworks applied to model the wavebreaking of the primary wave (a) and detailed in section 4. For the global approach (b), we consider the stability of the horizontally averaged concentration profiles. For the local approach (c), we study the development of small perturbations of wavenumber k wb at the node of the primary wave.

the conditions for a small disturbance to develop around the mean profile characterizing a Faraday wave of given amplitude (see figure 5). By contrast, the second model proposed relies on the local analysis of small perturbations at the interface node of the primary wave (see also figure 5).

The global approach

A simple model equation

In order to derive a simple model for the breakdown of Faraday waves, the concentration and velocity fields are decomposed into a mean and a fluctuating part as C = C + c and U = U + u. A mean quantity * is obtained by averaging along the horizontal x, y directions. The system (2.1) is classically averaged also in order to find the equations for the mean flow and its fluctuations. It can be directly shown that the mean velocity is zero, U = 0, due to the symmetries and the incompressibility condition. The mean vertical concentration profile C(z, t) evolves principally due to the vertical buoyancy flux wc as

∂ t C + ∂ z wc = D∂ 2 zz C. (4.1)
In this global approach, the primary Faraday wave is thus embedded in the mean vertical concentration profile C(z, t) but also has fluctuation components satisfying (4.1).

Further simplifications are now made in order to obtain a more tractable model. We seek a 'rapid' secondary instability occurring at small scales and located at z = 0. In this context, the rapid acceleration theory initially developped by [START_REF] Hunt | Rapid distortion theory and the problems of turbulence[END_REF] and applied to turbulent mixing layer in [START_REF] Gréa | The rapid acceleration model and the growth rate of a turbulent mixing zone induced by Rayleigh-Taylor instability[END_REF] provides a convenient framework for expressing the dynamics of small scale disturbances. We thus discard all the dissipative and the non linear terms, except for the coupling between the fluctuations and the mean concentration field. In addition, the small scale disturbance only sees a uniform mean concentration gradient at z = 0 determined by the length L = -1/∂ z C(0, t). Within this quasi-homogeneous approximation, the small scale fluctuating quantities and pressure p are determined by

∂ t u = -∇p -2AG(t)ce z , (4.2a) ∂ t c = w L(t) , (4.2b) ∇ • u = 0. (4.2c)
One recognizes in (4.2) the classical equations for an internal gravity wave, except it is driven by a time varying acceleration and mean density gradient. When it is uniform across the layer, this mean gradient can be evaluated from the mixing layer length [START_REF] Andrews | A simple experiment to investigate twodimensional mixing by rayleightaylor instability[END_REF] and previously used in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] for the fully turbulent regime. However, this property is lost when the flow takes the form of a single laminar wave. In this case, the inverse concentration gradient is maximum at z = 0, and we will discuss in a next section different ways to evaluate it.

L int = 6 +∞ -∞ C(1 -C)dz introduced by
These waves depends on their orientations but for this heuristic model we focus only on waves with a wavevector in the horizontal plane. Differently oriented modes are thought to be less relevant in the secondary instability partly because they are less likely to modify the mean concentration profile; the feedback of the fluctuations to the mean concentration profile is indeed controlled by the vertical buoyancy flux term wc which is weaker for vertically oriented modes. Eliminating w in (4.2), we obtain the following Mathieu-like equation (see a more detailed derivation in [START_REF] Gréa | What is the final size of turbulent mixing zones driven by the faraday instability[END_REF])

c + L L ċ + 2AG(t) L c = 0. (4.3)
The concentration fluctuations c are therefore fully determined by L expressing the amplitude of the primary Faraday wave. We need to determine the condition on L for which the perturbations may develop. Indeed, the rise of the perturbation foreshadows the breaking process of the primary wave and the onset of turbulence. More precisely, the equation (4.3) exhibits the buoyancy frequency defined as Ω B = (2AG 0 /L) 1/2 and a damping term L/L expressing the variations of Ω B as the mixing zone width L evolves. As will be seen in section 4.2, the Ω B frequency is relevant for the secondary instability because the shear at the nodes of the primary wave is directly driven by the wave amplitude.

We now detail the implications of this model equation regarding the wavebreaking phenomenology as observed in figure 1.

The subcritical nature and the criterion for the wavebreaking

The stability of the model equation (4.3) has been extensively discussed in [START_REF] Gréa | What is the final size of turbulent mixing zones driven by the faraday instability[END_REF] and [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF] allowing the prediction of the final widths of the turbulent mixing zones. This saturation criterion evaluates when the subharmonic instability stops or equivalently when the inner frequencies of the layer are no longer in

Ω 2 0B /ω 2 , F ) (with Ω 0B = (2AG 0 /L 0 ) 1/2
). The two colored areas indicate the first subharmonic tongues with β = 0 and β = 0.7 respectively. The dashed blue lines correspond to the approximation (4.4) while the red dashed dotted ones correspond to (B 9). resonance with the forcing. Thus, it did not seek to account for the unsteadiness of L. Conversely, we wish to interpret the wavebreaking as the development of a secondary subharmonic instability at small scales. For a small disturbance characterized by the buoyancy frequency Ω B , we thus aim at finding when it becomes parametrically unstable as the result of both the forcing and the primary wave oscillations. The enlargement of the primary wave amplitude determines not only the instability threshold, but also the later amplification of the secondary instability growth rate. This explains why the secondary instability rapidly develops at the interface and sheds light on the subcritical nature of this secondary instability. This peculiarity is indeed well-known for non linear Mathieu systems such as (4.3) as detailed in [START_REF] Soliman | Indeterminate sub-critical bifurcations in parametric resonance[END_REF] or [START_REF] Godrèche | Hydrodynamics and nonlinear instabilities[END_REF].

In order to derive an analytic criterion for the wavebreaking, we consider the inverse mean concentration gradient L having the following simple form L(t) = L 0 (1+β cos(ωt)). The length L is thus expected to be proportional to the amplitude η p of the primary Faraday wave, L(t) ∼ |η p | in the laminar phase. Here, the proportionality coefficient depends on the shape of the nonlinear primary wave. Also, the parameter β expresses the relative amplitude of the Faraday wave oscillations while L 0 is the mean over one oscillation period. This expression does not account for the primary mode growth, which is assumed small over an oscillation period. It also expresses that for a subharmonic instability, L oscillates at the frequency ω while η p oscillates at the frequency ω/2. However, the higher temporal harmonics of L(t) or η p for the primary wave, are discarded. In addition, it is important to note that the primary Faraday mode amplitude is in phase with the acceleration G(t).

In this context, we use the Floquet analysis detailed in appendix B to find the secondary subharmonic instability onset. The results are represented in the stability diagram of figure 6 exhibiting the instability tongues for β = 0 and β = 0.7. Note that the case β = 0 corresponds indeed to the classical stability diagram of a Mathieu equation. In this representation, the right-hand sides of the neutral branches (solid black lines of figure 6) determine the critical amplitude of the primary wave and the beginning of the secondary instability as L 0 grows. This gives a critical threshold L crit which should be close to the one characterizing the wavebreaking L wb if the instability develops quickly (we still have L crit L wb ).

An analytical approximation for the critical threshold L crit is also derived in appendix B leading to

L crit ≈ (4 -2F ) 1 + β/2 2AG 0 ω 2 , for F 1 and β 1. (4.4)
As shown in figure 6 (blue dashed line), the criterion (4.4) slightly underestimates L crit at small F and large β while being reasonably correct for the parameters taken in the experiment. However, it becomes very bad at large F , even leading to negative values for F 2. Despite being more complicated, a better approximation can be derived corresponding to (B 9) and shown by the red dashed dotted lines of figure 6.

In (4.4), we see that the forcing F together with the movement of the primary wave characterized by β contribute to the destabilization of the primary wave. In particular, even without acceleration forcing (F = 0), the secondary instability can be triggered by the primary wave oscillations. This also leads to striking differences in terms of growth rates. For instance, for the subharmonic resonance at F = 0.7 we find a Floquet exponent λ = 0.09 for β = 0 (also corresponding to the growth of the primary wave), while λ = 0.2 for the case with β = 0.7. Therefore, the acceleration induced by the primary wave increases sensibly (but not drastically) the growth rate of the secondary modes. This effect explains why the breakdown occurs in a time scale much shorter than the growth of the primary wave.

The local approach

The previous global approach has the main advantage of being relatively simple through relying on the horizontal averaging process. The drawback, however, is the loss of track of the physical mechanism responsible for the secondary instability. Also, it assumes that the secondary instability mode only results from the interaction with the mean component of the primary wave. This assumption may appear excessive and it motivates us to propose a complementary method. We therefore perform a stability analysis of the flow generated at the node of the primary wave.

The equations driving the evolution of an interfacial perturbation amplitude, η, at the node of the primary wave, are given in the inviscid limit by (see the derivation in appendix C)

η -2iAk wb U η + AG(t)k wb -k 2 wb U 2 -iAk wb U η = 0. (4.5)
In (4.5), k wb represents the wavevector modulus of the perturbation, and U is the oscillating tangential velocity induced by the primary wave at the node. Not surprisingly, we recover the equation in the Boussinesq limit first derived by [START_REF] Kelly | The stability of an unsteady Kelvin-Helmholtz flow[END_REF] and expressing the dynamics of an oscillating sheared interface. In fact, supposing the instability am-Cavelier et al.

plitude and wavelength are small, respectively kη 1 and κ = k wb /k 1, this allows us to neglect the interface tilting at the node and to assume the quasi-homogeneity of the perturbation. Note that often at least two vortices may appear in the experiments at the node of the primary wave during the breaking process suggesting the validity of the homogeneity assumption.

Many works have been dedicated to the stability analysis of (4.5) in the case of G constant and U oscillating at a single frequency ω. In addition to the parametric resonant modes [START_REF] Kelly | The stability of an unsteady Kelvin-Helmholtz flow[END_REF] identified, [START_REF] Lyubimov | Developement of a steady relief at the interface of fluids in a vibrational field[END_REF] and [START_REF] Khenner | Stability of planeparallel vibrational flow in two-layer system[END_REF] have shown and derived a criterion for the existence of Kelvin-Helmholtz type modes at the interface. These latter modes generate a longwave instability observed in most experiments and which may exhibit frozen wave patterns at high forcing frequency [START_REF] Wolf | Dynamic stabilization of the interchange instability of a liquid-gas interface[END_REF][START_REF] Wunenburger | Frozen wave induced by high frequency horizontal vibrations on a co2 liquid-gas interface near the critical point[END_REF][START_REF] Yoshikawa | Oscillatory kelvinhelmholtz instability. part 2. an experiment in fluids with a large viscosity contrast[END_REF][START_REF] Gaponenko | Dynamics of the interface between miscible liquids subjected to horizontal vibration[END_REF][START_REF] Lyubimov | Viscosity effect on the longwave instability of a fluid interface subjected to horizontal vibrations[END_REF]Gréa & Briard 2019). Noticeably, frozen waves are completely steady structures analogous to the inclined equilibrium positions of a strongly and horizontally oscillated pendulum.

The local analysis emphasizes the importance of the shear in the breakdown process of the primary wave. Although the works of [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] and [START_REF] Kalinichenko | Developement of a shear instability in nodal zones of a standing internal wave[END_REF] have early recognized this aspect, the nature of the instability, either Kelvin-Helmholtz (KH) or parametric resonance (PR) type, has not been fully identified in this Faraday wave context. The results of [START_REF] Kelly | The stability of an unsteady Kelvin-Helmholtz flow[END_REF] and [START_REF] Khenner | Stability of planeparallel vibrational flow in two-layer system[END_REF], do not apply to our specific configuration where the acceleration oscillates at frequency ω while the shear velocity is subharmonic with frequency ω/2. We thus reconsider the problem of (4.5) by taking a primary wave of the form η p (t) = a cos(ωt/2) leading to

U (t) = - ωa 2 sin(ωt/2), and U 2 = ω 2 a 2 8 (1 -cos ωt) = 1 2 AG 0 k 1 + 4∆ a 2 (1 -cos ωt), (4.6)
where in the last expression, the subharmonic resonance condition ω 2 = 4AG 0 k/(1+4∆) for an inviscid interface is used. In order to study the stability of (4.5), it is further convenient to introduce the new variable Y defined as η = Y e t 0 iAk wb U (t )dt giving (at leading order in A)

Ÿ + AG(t)k wbk 2 wb U 2 Y = 0. (4.7) Due to the change of variable expression and that U oscillates at ω/2, the response in η will be also subharmonic. However, for small Atwood number, η ≈ Y and the response can be also nearly synchronous. Indeed, by substituting the expression for U (t) into (4.7), we obtain a simple Mathieu equation on the form

Ÿ + (P + Q cos(ωt)) Y = 0, (4.8) with P = AG 0 k κ - 1 2(1 + 4∆) κ 2 (ka) 2 and Q = AG 0 k κF + 1 2(1 + 4∆) κ 2 (ka) 2 .
In figure 7, we show the stability diagram of (4.8) in a κ-ka plane (the subharmonic resonance condition for the primary wave is again used). The instability tongue corresponding to the KH type modes appears for P 0 which stands as the classical criterion for the inviscid KH instability, k wb U 2 AG 0 [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF]. The parametric resonance zones start for P 0 but have also a continuation in the opposite half-plane. Remarkably, the instability zones exhibit a very weak dependence on κ for κ 1. Therefore, at given κ and as the primary wave amplitude grows, the perturbation passes through the successive instability zones, first the PR types then lately the KH one. The growth rates can be computed with the Floquet exponent and show a maximum approximately in the middle of each zone. The KH and PR1 growth rates are larger compared to the other instability zones. Therefore, the breakdown of the Faraday wave is expected to occur when the wave steepness, ka, lies in the instability KH or PR1 zones. This local theory is inviscid which explains that the growth rates are higher at large κ. Of course, the viscosity and the thickness δ of the interface should play a role and moderate this aspect.

Similarly to the global approach, we can propose an approximation for the critical wave steepness corresponding to the onset of the PR1 instability. Using the asymptotic expression for the neutral curves of the Mathieu equation in the limit P → -∞, Q → +∞ detailed in [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], we obtain

ka crit ≈ 1 3 (1 + 4∆)(1 + F ). (4.9)
This simple expression (4.9) corresponds indeed to the plateau (it does not depend on κ) separating the PR1 and PR2 bands in the small perturbation wavelength limit κ 1 as shown in figure 7 (black dotted line).

At this stage, two theoretical approaches have been proposed to study the breaking process of Faraday waves. Before assessing the validity of these approaches, it is first necessary to detail how we detect the wavebreaking in both the experiments and simulations.

Cavelier et al.

Data analysis of the experiments and simulations

In this section, we detail the methodology used in order to measure the primary wave amplitudes and the inverse mean density gradients. We then investigate various strategies to detect the instants at which the wavebreaking occurs. In the last part we discuss the validity of the global and local wavebreaking theories against the data from the experiments and DNS.

Primary wave amplitude and inverse mean density gradient

The global and local approaches detailed in section 4 rely on different scales characterizing the primary wave. These correspond respectively to the inverse mean concentration gradient L(t) at the rest interface position z = 0 and the wave amplitude η p (t) (see also figure 5). It is therefore important to link these two quantities in order to compare the global and the local theories in the experiments. As already remarked, we can expect the proportionality between the inverse mean gradient and the wave amplitude, L(t) = κ|η p (t)|. This assumes that the shape of the primary wave is frozen which also determines the coefficient κ. By considering the position ξ p (x, t) = η p (t) sin kx of a sharp sinusoidal interface, as for the standing wave in appendix C, we can easily show that κ = π (the wave amplitude is supposed larger than the interface thickness δ). This choice is also motivated by the expression for the finite amplitude standing wave profile given by [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] in deep water approximation and showing that the higher order corrections are negligible even at moderate wave steepness.

The global theory assumes an inverse density gradient oscillating as L(t) = L 0 (1 + β cos ωt) while in the local theory the primary wave amplitude evolves as η p = a cos ωt/2. By expanding |η p (t)| in Fourier series and truncating at leading order, we simply get L 0 = 2κa/π and β = 2/3. This allows to express the results from the global theory in terms of critical steepness as for the local theory. In particular, the critical wave steepness for the global theory expressed by (4.4) in the limit of small F , using the resonance condition and taking κ = π (sinusoidal interface) gives

ka crit = 3 8 (2 -F )(1 + 4∆).
(5.1)

We now detail how to measure the inverse local gradient L and the wave amplitude η p in practice in the experiments. The inverse local gradient is a difficult quantity to obtain directly from the mean concentration profile. The latter, resulting from the postprocessing of the images of the camera as detailed in [START_REF] Briard | The turbulent faraday instability in miscible fluids[END_REF], can be a bit noisy, particularly when the secondary instabilities start at the node. The quantity |η p | can be obtained either from the crest-to-crest amplitude of the wave on the raw camera images, or from an arbitrary threshold on the mean concentration (here taking the height where 0.1 C 0.9). Again this method suffers from being sensitive to small variations in the mean concentration profile. Therefore, we find it more convenient to measure these quantities indirectly using the integral length L int previously introduced in section 4.1.1. By assuming the shape of the concentration profiles frozen, we can deduce all the characteristic lengths from L int which for a sinusoidal interface gives L int = 2.4|η p | = 0.76L.

To ensure that these relations apply in our problem, we have measured the amplitude |η p | at a local maximum just before the wavebreaking using the threshold method on C, the inverse mean concentration gradient L and the integral length L int . In figure 8a-b, we see that the sinusoidal profile is a good fit of the interface position and provides a good evaluation of the mean concentration profile. If the amplitude is too small though, the fit is less satisfactory probably because the interface thickness should be accounted for in the evaluation of the mean concentration profile. The method shows that the inverse mean concentration gradient is maximum at the center of the layer z = 0, in accordance with the fact that the wavebreaking starts at the nodes of the primary wave. Besides, we recover the expected relation between L and the wave amplitude |η p |. Therefore, the different lengths L and |η p | can be correctly evaluated from L int . In figure 8c, we further check that the relation between L int and |η p | holds on multiple experimental data before but close to the wavebreaking. The correlation is thus satisfactory and gives us confidence in our method to extract the amplitudes and the inverse concentration gradients necessary to explore the wavebreaking phenomenon.

Wavebreaking detection

In this part, we detail the strategy in order to identify the instants corresponding to the wavebreaking in the experiments and the simulations.

The wavebreaking phenomenon can be defined as a local overturning of the interface at the node. This can be directly observed in the images from the camera by noticing the appearance of vortices giving birth to a mushroom-like structure (see figure 9). Since this visual criterion is somewhat subjective, we also propose an automated detection procedure using a simple algorithm based on the Thorpe displacement (introduced by [START_REF] Thorpe | Turbulence and mixing in a scottish loch[END_REF]). Indeed, this displacement characterizes the distance a parcel of fluid has to move vertically in order to be in stable equilibrium with the surrounding water. Hence, a completely stable profile would have zero displacements whereas the existence of some local overturns generates non zero Thorpe displacements. Accordingly, we compute each 1D vertical transect the displacements δ T = z *z, where z is the position of a fluid parcel on the instantaneous concentration field, and z * is the position of the same parcel on the vertically sorted concentration field (see figure 9). We detect when this displacement |δ T | exceeds a certain threshold (here |δ T | > 5 -12 px = 5.5 -13.2 mm) chosen a bit smaller than the initial interface widths. For a given experiment, this algorithm provides almost the same image (and thus same time and value for L) as the one chosen by eyes only (see figure 9). However, it is not possible to detect from this method the beginning of the 'blurred' region as defined in [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] and indicating the secondary instability onset leading to the wavebreaking at the nodes.

The Thorpe displacement can also be used to measure the vortex size at wavebreaking and thus gives the ratio κ = k wb /k necessary for the local theory. To determine k wb , we take the maximum displacement δ T evaluated at the image given by the Thorpe displacement method. By construction, this measurement cannot be smaller than the arbitrary threshold chosen for the wavebreaking detection. In practice, δ T exceeds by more than three times this value. As a summary, we illustrate in figure 10 the whole procedure allowing for the wavebreaking detection and the measurements of the wave amplitudes and inverse concentration gradients.

The time evolution of the integral length L int in a simulation revealing the subharmonic oscillations (thus L int oscillates at frequency ω) and the growth of the primary wave are presented in figure 10a. The instant corresponding to wavebreaking is determined by evaluating the Thorpe displacement. This corresponds to the apparition of the vortices at the nodes. In addition, we show the evolution of the integral length L int,s computed Cavelier et al.

from the sorted concentration profiles as in [START_REF] Briard | Harmonic to subharmonic transition of the faraday instability in miscible fluids[END_REF]. This quantity expresses the evolution of irreversible mixing by distinguishing the available and the background potential energies (see also for instance [START_REF] Winters | Available potential energy and mixing in density-stratified fluids[END_REF][START_REF] Peltier | Mixing efficiency in stratified shear flows[END_REF][START_REF] Davies Wykes | Efficient mixing in stratified flows: Exerimental study of a rayleigh-taylor unstable interface within an otherwise stable stratification[END_REF]. Noticeably, we observe that the growth of the irreversible mixing starts just after the wavebreaking time. In figure 10, we verify further that this process is due to wavebreaking as the curve of L int,s gets detached from the interface width L int,d expressing the thickening of the interface by pure diffusion. However, the slow evolution of L int,s would make it difficult to build a wavebreaking detection criterion from it. In figure 10b, we also present the different mean concentration profiles before and after the wavebreaking time and renormalized by L int . This confirms that the mean concentration profiles can be considered as frozen and well represented by a sinusoidal interface before the wavebreaking. This aspect is important in order to link the inverse concentration gradient to the wave amplitude and to compare the global and the local approaches. At and after the wavebreaking, the mean concentration profiles are drastically distorted, and exhibit inversions of the stratification due to the roll-up of the interface. The predictions from the local and global theories are also plotted in figure 10a showing a good agreement with the measured amplitude at the wavebreaking. Again this is expected to have L wb > L crit for both theories. We discuss more thoroughly this point in the next section.

Results and discussion

Critical steepness values

We now analyze the experimental and numerical measurements of wavebreaking in order to assess the global and local theories presented in section 4. The instants corresponding to the wavebreaking are detected using the procedure detailed in section 5.2. They are usually close but not necessarily exactly located at a maximum of the primary wave. Therefore, we perform a linear interpolation between two successive maxima in order to evaluate the amplitude a at the wavebreaking (thus deduced from the integral scale L int as explained in section 5.1). In figures 11 and 12, the wavenumber ratio κ and the wave steepness ka corresponding to the wavebreaking are reported and superimposed to the stability curves obtained from the local and the global approaches.

The wave steepness measured at the wavebreaking detection is roughly ka ∼ 0.75 from the experiments and the DNS. These values are thus located inside the parametric resonance instability band above the thresholds provided by the local and global approaches. Some points are slightly below the critical steepness provided by the global theory at ∆ = 0.04 in figure 11d and 12d but it corresponds to negative detuning cases which have a lower critical steepness value. For the moderate F parameters investigated in this work, both theories predict roughly the same critical steepness value, around ka crit ∼ 0.5. This result clearly gives a strong credit to a wavebreaking process due to a secondary subharmonic instability appearing when the primary wave reaches a critical amplitude: it is therefore subcritical. Conversely, we have conducted a simulation (DNSe case, see table 2 and figure 3) with parameters close to the neutral curve and where the primary wave reaches a saturation amplitude below the critical steepness value. This kind of wave does not experience wavebreaking as expected from the theory.

A tendency emerges from figure 11, that the critical wave steepness slightly decreases at large wavenumber ratio κ. This observation stands despite the difficulty to evaluate κ in the experiments but is less obvious in the DNS results of figure 12. Although the neutral curves depend very weakly on κ (they do not depend on it with the global theory), the instability growth rates evaluated from Floquet theory increase when the small perturbation wavelength κ becomes large in the local theory. This result is characteristic of classical inviscid theories such as for the KH instability. It is therefore not surprising to see that the wavebreaking occurs earlier, i.e. at smaller ka, when the secondary instability is initiated at smaller wavelengths. Also, one cannot exclude the possibility that the unstable modes from the higher parametric bands (PR2 for instance) play a role in the destabilization process of the primary wave. More precisely it may explain the blurring of the interface already observed by [START_REF] Kalinichenko | Developement of a shear instability in nodal zones of a standing internal wave[END_REF] well before the wave breaking at ka ∼ 0.4. These PR bands are however expected to be seriously damped by the viscosity in addition of having a lower growth rate. Conversely, the growth rates of the PR1 modes as computed in the local theory are very large compared to the primary Faraday wave ones. We find indeed Floquet exponents λ, varying between 2 and 8, corresponding to an amplification 500 over one oscillation period (by contrast the primary wave has Floquet exponents around 0.08). This is due to the additional forcing generated by the primary wave and explains why the breaking process develops over one or two oscillation periods. It seems also justified to neglect in the theory the growth of the primary wave amplitude for these parameters. By contrast, the growth rates evaluated by the global approach are much more modest (∼ 0.2). This point is probably due to the assumption 2).

linking the local mean concentration gradient to the amplitude of the primary wave detailed in section 5.1. Indeed, as soon as the secondary instability develops, the shape of the mean concentration profile is distorted due to the perturbation feedback. This aspect is well demonstrated in the simulation result of figure 10. It eventually leads to a local change of sign in the stratification which would in return considerably increase the instability growth rate. However this effect is not accounted for in the present global approach which still remains effective at predicting the secondary instability onset.

A large range of κ values has been found in the experiments and simulations, which comes equally from the primary and secondary wavelength measurements. Noticeably, the simulations provide larger (up to 110 in figure 12) partly because the primary wavelength is large. Therefore, an important question raised is how the wavenumber k wb is selected. Neither the global theory (which does not depend on k wb ) nor the local theory (which relies on the analysis of an infinitely thin interface) can explain this aspect. When the interface thickness reaches comparable size with the secondary instability wavelength, the natural frequency becomes bounded and the instability growth rate is expected to be limited. By analogy, the pure KH instability (represented by the Taylor-Goldstein equation) has a maximum growth rate around k wb δ ∼ 1 at low Richardson number. This also depends on the shear and density profiles considered for the analysis (see for instance [START_REF] Taylor | Effect of Variation in Density on the Stability of Superposed Streams of Fluid[END_REF][START_REF] Hazel | Numerical studies of the stability of inviscid stratified shear flows[END_REF][START_REF] Caulfield | Multiple linear instability of layered stratified shear flow[END_REF]. The same process for the PR mode seems at work in our experiments as suggested by the measurements of k wb . This has been also evidenced theoretically by [START_REF] Poulin | Parametric instability in oscillatory shear flows[END_REF] on oscillatory uniform shear layer configurations. However, the local analysis accounting for the interface thickness is rendered complex by the coupling between the internal layer modes induced by the time varying shear and acceleration.

Exploring further the forcing parameter effect

Despite performing well for the parameters corresponding to the experiments, the local and global approaches differ sensitively in their dependence on the forcing parameter. While the critical wave steepness grows linearly in F in the local theory as shown by (4.9), the F dependence in the global criterion is more complex, even exhibiting a diminution of the wave steepness at small F as stated by (5.1). In order to investigate this dependence, we show in figure 13 the results from the simulations corresponding to the DNSd cases (see table 2). These simulations are conducted using the same initial conditions and frequency ω, varying only the forcing parameter F . The results show a critical wave steepness globally evolving linearly with . This trend seems reasonably well predicted by the local approach although underestimating the wave breaking at large F . By contrast, the global approach clearly underestimates the critical wave steepness. As an explanation, it should be stated that both the global and local approaches become limited when F is large because the growth of the primary wave amplitude a cannot be neglected over one oscillation period. In the DNSd cases of figure 13, the wavebreaking occurs during the first oscillation for F 3, when the sign of gravity is inverted. Therefore at large F the secondary instability is expected to change nature as becoming triggered more by the growth of the primary wave than its oscillations. In that respect, the wavebreaking process is no longer parametric but becomes of KH type similar to the one appearing in Rayleigh-Taylor instability [START_REF] Birkhoff | Helmholtz and Taylor instabilities[END_REF][START_REF] Daly | Numerical study of two fluid rayleightaylor instability[END_REF][START_REF] Baker | Singularity formation during rayleightaylor instability[END_REF].

The final transition to turbulence

So far, we have not demonstrated that the secondary instability developing at the node of the primary wave triggers the full transition to a turbulent regime. Indeed, the simulations with a 2D interface initialization are unable to develop the cascading process 2. The wavebreaking is detected at ωt = 14.5 for both simulations. leading to turbulence. We thus provide as an example the simulation DNSf where the 2D initialization is perturbed by a small random white noise in the spanwise y direction.

In figure 14, we compare at different times DSNf with the simulation DNSd3 using the same parameters (see table . 2) but having a simple 2D initial condition.

Both simulations evolve similarly until wavebreaking (here at ωt = 14.5). This means also that the critical amplitude is not modified by the small spanwise y perturbation, hardly visible since its amplitude is two orders of magnitude smaller than that of the x 2D perturbation. This feature has been reproduced on other simulations at different ε 1 parameters or in simulations where the spanwise invariance is broken by the presence of lateral walls (not presented for the sake of conciseness). After wavebreaking, the simulations differ sensitively as DNSf (or simulations with the spanwise invariance broken) exhibits a rapid transition to turbulent mixing. The process is so violent that turbulence does not stay bounded to the node of the Faraday wave and spreads around the whole layer. The images of figure 14 at ωt = 25.7 (see also online movie of DNSf) reveal that the final mechanism leading to turbulent mixing is related to the merging and the stretching of the secondary vortices at the node forced by the oscillations of the Faraday wave.

In order to identify the specificity of this transition, we evaluate the Reynolds and bulk Richardson numbers as classically introduced for the Kelvin-Helmholtz instability [START_REF] Caulfield | Layering, instabilities, and mixing in turbulent stratified flows[END_REF]. Taking for the interface half width σ and for the half shear velocity U = aω/2, we find Re = U σ/ν = 313 at the wavebreaking ωt = 14.5. For the bulk Richardson number we obtain Ri = AG 0 σ/U 2 = 0.023. While the low value of the Richardson number is consistent with the existence of a strong shear instability, the Reynolds number is unexpectedly small for a mixing transition. This is despite being probably underestimated due to the blurring of the interface before the wavebreaking. Perhaps the fact that the transition does not result from a convective secondary instability inside the vortices as often observed for the Kelvin-Helmholtz instability [START_REF] Salehipour | Turbulent diapycnal mixing in stratified shear flows: the influence of prandtl number on mixing efficiency and transition at high reynolds number[END_REF] may explain this aspect. By contrast, we can consider a global Reynolds number based on the width and velocity of the layer. This gives Re = a 2 ω/ν ∼ 10 4 and agrees with a mixing transition criteria of Re > 10 4 proposed by [START_REF] Dimotakis | Turbulent mixing[END_REF].

Conclusion

In this work, we have studied experimentally, numerically and theoretically the wavebreaking mechanism leading to turbulence of growing Faraday waves at the interface between miscible fluids of small density contrast. The Floquet linear theory reveals that several subharmonic modes can be simultaneously unstable due to the quantification induced by the tank geometry and the large forcing accelerations considered. It is shown that, by accounting for the viscous damping at the walls and from the bulk flow, the parameters in our experiments are located well above the instability threshold, suggesting that the inviscid theory provides a good first approximation for this problem. We have evidenced that the mode selection of the Faraday wave results not only from a linear process but also a nonlinear competition favoring the modes with smaller wavelength. More precisely, the natural frequencies of the system decrease as the mode amplitude grows. Therefore, the subcritical modes are those with positive detuning which are favored when the primary wave amplitude has reached a critical value. By contrast, those with negative detuning are supercritical and become damped as being less in parametric resonance with the forcing frequency. This mode competition phenomenology also explains the sensitivity to initial conditions and the symmetry breaking of the dominant mode which has been observed in our experiments.

The idea of considering the nonlinear mode interaction as a linear process but with a time evolving reference state was first successfully introduced to describe the saturation of turbulent Faraday waves. In addition to use it for describing the mode competition phenomena, we also employ this approach to explain the breakdown of the primary wave as the appearance of a secondary subharmonic instability at small wavelengths. This theory, referred to here as 'global' as resulting from a horizontal averaging process, reveals that the secondary instability is principally due to the oscillations of the primary wave, thus explaining why it develops very rapidly compared to the growth rate of the primary Faraday wave. We then also propose a criteria giving the critical steepness of the primary wave at which the wavebreaking is expected to appear. Due to its simplicity, this theory cannot explain why the secondary instability occurs at the node of the primary wave nor how it depends on its wavelength. This is why we study the flow in a local frame attached to the node of the primary wave leading to a local approach. This reveals Cavelier et al.

the importance of the shear in the development of the secondary instability leading to the wavebreaking. It is shown from a stability analysis that the unstable modes can be of either Kelvin-Helmholtz or of parametric resonant type, the latter occurring earlier during the growth of the primary wave.

In order to assess the global and local theories, we evaluate the primary wave amplitudes in the experiments and in the numerical simulations, together with the wavenumbers associated with the wavebreaking phenomenon. This is done using the Thorpe displacement indicating a local overturning of the interface, as well as extracting the background density profiles evolving due to an irreversible mixing. The results reveal the subcritical nature of the wavebreaking detected roughly for wave steepnesses ka ∼ 0.75. These values are consistently larger than ka ∼ 0.4 reported by [START_REF] Kalinichenko | Developement of a shear instability in nodal zones of a standing internal wave[END_REF] but associated with the earlier 'blurring' stage of the instability. Remarkably, both theories indicate that this should correspond to a parametric resonant subharmonic instability developing when the primary wave amplitude reaches a critical value. Therefore the present mechanism for Faraday waves presents similarities with the breaking process of internal gravity waves. However, when the forcing parameter is increased, our approaches come to a limitation as the wavebreaking process changes in nature resulting more from the amplitude growth of the primary wave than its oscillations. In this case, the phenomenon becomes similar to the secondary vortices appearing in the classical Rayleigh-Taylor instability.

Simulations with 2D initial conditions perturbed along the spanwise direction evidence that the final transition to turbulence originates from the secondary instability at the node of the Faraday wave. The stretching and merging of the secondary vortices driven by the oscillations of the primary wave is an efficient mechanism to produce mixing, differing sensitively with the transition scenarios observed in the context of the Kelvin-Helmholtz instability.

A =                       
                       , ( B 5)
and

B =                        1 + F 2 0 F 2 0 0 1 -F 2 0 F 2 . . .
F 2 0 1                        . ( B 6)
with an = (n -1 + α)(n + α), bn = (n + 1 + α)(n + α). For the diagonal terms in matrix A, we have dn = (n + α) 2 for all n > 0 and for n = 0, we have d ± 0 = α 2 ± β 2 a0. In order to get a simple analytic approximation for the neutral branch for However, this expression is not completely satisfactory to approximate the neutral branch corresponding to the onset of the subharmonic instability for larger β and F as shown in figure 6. A better expression can be found with N = 1 which leads to The equations describing the field generated by an insviscid interfacial wave of small amplitude can be found in many classical textbooks (see for instance Sutherland 2010). We recall briefly the procedure, expressing a 2D incompressible velocity disturbance by its stream function as (up, wp) = (-∂zψp, ∂xψp). Seeking modal solutions on the form ψp(x, z, t) = Ap(t) ψp(z)e ikx , and for the interface deformation ξp(x, t) = ηp(t)e ikx , we thus obtain the degenerate Rayleigh equation : We turn our attention to continuity of pressure on each side:

- At this stage, we wish to perform the stability analysis of this flow in the vicinity of the node, x = 0, z = 0. We thus rescale the dimensions with the typical wavenumber k wb of the secondary instability. Therefore, in the small perturbation wavelength limit of κ = k wb /k 1, the flow induced by the primary wave simply reduces to an horizontal interface subjected to a discontinuous tangential velocity ξp = 0, up = ± ηp = ±U and wp = 0.

(C 7)

C.2. Secondary instability

We perform the linear stability analysis of the base flow coming from the primary wave and defined by (C 7). We consider the small velocity perturbation (u, w) = (-∂zψ, ∂xψ). Again, seeking modal solutions on the form ψ ± (x, z, t) = A ± (t) ψ± (z)e ik wb x , and for the interface deformation ξ(x, t) = η(t)e ik wb x . We obtain the degenerate Rayleigh equation

?Figure 1 :

 1 Figure 1: The breaking of a Faraday wave in the FARAMIX experiment. (a): Visualization showing the tank geometry and the configuration (also presented in Briard et al. (2020)).(b): Time series images from the camera zooming on one wavelength and presenting two oscillation periods of the primary wave. This illustrates the different stages of the wavebreaking with first a 'blurring' of the interface at the node followed by a 'roll-up'. This case corresponds to the b5 experiment which parameters are detailed in table 1. (c): Visualisation of the interface at wavebreaking in the direct numerical simulation DNSd3 (the parameters are given in table2). The reference frame as well as the acceleration direction are also indicated.

Figure 2 :

 2 Figure 2: Stability diagram for (3.1) in a non dimensional frequency ω/ AG 0 k 2,0 and forcing F plane. The colored regions correspond to the first subharmonic instability band associated with the different modes of the tank (the mode number is indicated in the figure). The diagram is obtained using the damping coefficient γ = γ δ + γ w at three different Atwood numbers A and considering an interface thickness δ = 1 cm.The neutral curves (thick plain lines) have a slight dependence on the Atwood number explaining that they are not completely superimposed. The symbols correspond to the parameters taken in the experiment in table 1. The shapes indicate the Atwood number and the colors reveal which mode is eventually selected.

Figure 3 :

 3 Figure3: Parameters of the DNS (symbols) in the stability diagram ω-F . The colored areas correspond to the subharmonic instability tongues for the modes (2, 0) and (3, 0) accounting for viscosity and diffused interface δ = 1.8 cm (as no walls are present in the DNS, the damping coefficient γ ≈ γ δ and critical thresholds F th are very small). The symbols' color indicates which mode emerges from the simulation, the mixed colors express that both modes can be observed. Two series of DNSb, c (see table 2) are presented here starting from point A. In the DNSb group, the initial amplitude ratio r between modes (3, 0) and (2, 0) is set at r = 0.5 and we decrease the forcing frequency ω. In the DNSc group, the frequency and forcing are fixed and we vary r. The point corresponding to DNSe is also placed.

Figure 4 :

 4 Figure 4: Mode selection in 6 DNS of figure 3. (a) Cases corresponding to DNSc (see table 2) with r varying and A = 0.030, ω = 3.07, F = 0.8. The amplitude of the interface deformation is = 3 cm. (b) Cases corresponding to DNSb (see table2) with ω varying and A = 0.030, F = 0.8, r = 0.5, = 1.5 cm. We put the slices of width W of the concentration field at four instants starting from the initial interface at t = 0 and ending when the wavebreaking occurs; pure fluids remain in white while mixed fluid with C = 0.5 ± 0.15 is in black.

Figure 6 :

 6 Figure 6: Stability analysis of (4.3) with L(t) = L 0 (1 + β cos(ωt)) and represented in the plane (Ω 2 0B /ω 2 , F ) (with Ω 0B = (2AG 0 /L 0 ) 1/2). The two colored areas indicate the first subharmonic tongues with β = 0 and β = 0.7 respectively. The dashed blue lines correspond to the approximation (4.4) while the red dashed dotted ones correspond to (B 9).

Figure 7 :

 7 Figure 7: Stability curve of (4.8) for F = 0.7 and ∆ = 0 in the κ-ka plane or the P -Q plane (insert) corresponding to the classical representation of the Mathieu equation. The blue colored areas show the Kelvin-Helmholtz (KH) and parametric resonance (PR) instability regions. The dashed curve corresponds to P = 0. The area corresponding to P < 0 is located above the dashed curve in the κ-ka representation. The critical wave steepness value indicated by the black dotted curve corresponds to criterion (4.9).

Figure 8 :

 8 Figure 8: (a) Visualisation of the interface in experiment EXPa7 (see table 1) at the amplitude maximum just before the wavebreaking and compared to a sinusoidal profile (red line). (b) Mean concentration profile from experiment and sinusoidal interface. (c) int plots as a function of |η p | at a maximum amplitude just before the wavebreaking for all the experiments of table 1. The values for |η p | are evaluated from the crest-to-crest distance of the wave measured directly on the images or using the mean concentration profiles. The two arrows correspond to the EXPa7 case shown in (a) and (b).

Figure 9 :

 9 Figure 9: Procedure for the wavebreaking detection. Left: Visual criterion from the calibrated camera image. Right: The Thorpe displacements δ T evaluated by sorting the concentration field in each vertical transect of the calibrated image. The colorbar indicates the displacement in [mm].

Figure 10 :

 10 Figure 10: Time evolution of the mixing zone width and the mean concentration profiles at different instants extracted from the DNSd2 (see parameters in table. 2). (a) Evolution of the lengths L int and L int,s computed from the mean and sorted concentration profiles, respectively, as a function of ωt. The dashed line corresponds to the integral scale L int,d expressing the thickening of the interface by diffusion only. The star symbol at ωt = 17.75 indicates the wavebreaking detected by the Thorpe displacement. The horizontal lines correspond to the theoretical wavebreaking predictions, here converted in terms of integral length. (b) Mean concentration profiles at different times corresponding to the local maxima of L int in (a) and renormalized by the integral mixing zone width L int . The inserted images illutrate the state of the interface at the same instants.

Figure 11 :

 11 Figure 11: Parametric instability bands (colored regions) and experimental data (symbols) corresponding to the wavebreaking in a κ-ka representation. The parameters for the experiments are detailed in table 1. The left (a, c) diagrams correspond to the local theory (PR1 instability band). The right (b, d) ones are associated with the global theory. We place the theoretical instability zones with ∆ = 0 in the (a, b) figures and ∆ = 0.04 for the (c, d) figures both with F = 0.3 and 0.7. The symbol colors indicate successively (a) the forcing parameter F , (b) the Atwood number A, (c) the detuning ∆, and (d) the primary wave mode in the experiments.

Figure 12 :

 12 Figure 12: Same as figure 11 but for the DNS data corresponding to the series DNSa (see table2).

Figure 13 :

 13 Figure 13: Wave steepness at wavebreaking as a function of the forcing parameter F for the DNSd cases (symbols) with ω = 4.29 rad s -1 and A = 0.045 (see table 2). The dashed and dotted curves correspond respectively to the local and global criteria. Inserts: Visualisation of the concentration field in DNSd3 and 11 just after the wavebreaking.

Figure 14 :

 14 Figure 14: Visualization of the interface at different times of DNSd3 (left) and DNSf (right) corresponding to the parameters in table 2. The wavebreaking is detected at ωt = 14.5 for both simulations.

  The flow induced by the primary wave at the node

  = e ±kz on each side of the interface. Discarding the second order terms, the condition expressing the interface dynamics is wp = ikψp = of ψp across the interface and ikAp = ηp The equation for the horizontal momentum writes (again discarding the quadratic terms) ∂tup = -∂ 2 tz ψp = -

  this gives ψ+ = e -k wb z and ψ-= e k wb z . The condition at z = 0 is thereforew = ik wb ψ = Dξ Dt = ξ ± U ik wb ξ, , the continuity of pressure gives p + -ρ1G(t)η = p --ρ2G(t)η. (C11)

  For a standing wave, we thus have ξp(x, t) = ηp(t) sin kx, wp(x, t) = ηp(t)e

	ρ1G(t)ξp + ρ1	∂ 2 tz ψp ik	= -ρ2G(t)ξp + ρ2	∂ 2 tz ψp ik	,	(C 4)
	leading to
	ηp + AG(t)kηp = 0.			(C 5)

∓kz sin kx, ψp(x, t) = -ηp(t) k e ∓kz cos kx and up(x, t) = ± ηp(t)e ∓kz cos kx. (C 6)

Acknowledgment

The authors wish to acknowledge Y. Atik and L. Brosset from the Gaztransport Technigaz (GTT) Motion Analysis and Testing Laboratory for accessing the facility and technical support during the experiments in the context of the FARAMIX project.

Declaration of Interests

The authors report no conflict of interest.

Appendix A. The damping of a gravity waves in a rectangular tank

In this appendix, we evaluate the contributions from the walls to the damping of an interfacial wave in a rectangular tank. The geometry corresponds to the experiment detailed in section 2.2. We consider a wave with horizontal wavenumber k and a small amplitude a, oscillating at the frequency Ω (the damping is thus supposed to be very small). The rest interface position is located in the middle of tank and we study the two cases depending if the top boundary is a wall or a free surface. This result has already been established by [START_REF] Keulegan | Energy dissipation in standing waves in rectangular basins[END_REF] in the context of a free surface wave and generalized by [START_REF] Thorpe | On standing internal gravity waves of finite amplitude[END_REF] for two fluids configurations. The objective here is to disentangle the various contributions from the walls, in particular to show that the damping comes principally from the lateral boundaries in our experiment.

The subcritical transition to turbulence

31

We assume potential flow away from the walls. If the viscous boundary is very thin, the potential φ for the wave can be expressed by

= -aG 0 Ω cosh kH/2 cosh k(z -H/2) cos kx cos Ωt, for z > 0.

Hence, the 2D velocity field (for z < 0) is given by u 0 = ∇φ

sin kx cos Ωt,

For z > 0 one obtains the same expression but using z -H/2 in the cosh.

Classically (Landau & Lifshitz 2013), the damping coefficient can be evaluated from the total mechanical energy of the system as

where E mech is the sum of kinetic and potential energy, and the symbol * corresponds the mean over one oscillation period of the wave. In an oscillating system, the mean mechanical energy is provided by twice the mean kinetic energy. Therefore, using that kinetic energy comes principally from the potential flow, we have

In the limit of small Atwood number, the mechanical energies in the upper and lower parts of the tank are the same so we can restrict the integration over the lower part of the tank. Also, taking the mean over a period of time of cos 2 (Ωt) brings a factor half, we get

In order to compute the damping coefficient, we determine the mean value of the energy dissipation due to the walls. To this aim, it is necessary to evaluate the velocity field in the Stokes layers which should match the potential solution u 0 , w 0 oscillating at frequency Ω away from the walls. Following Landau & Lifshitz (2013), we use the fact that the oscillations of a viscous liquid around a solid body is equivalent to the oscillations of a solid body immersed in a viscous liquid. Hence, each wall can be assimilated to a plate oscillating in its own plane.

The mean mechanical energy dissipated in the layers adjacent to the wall is equal to the mean kinetic energy dissipated in those layers. For this, we use that the wavelength and dimensions of the tank are large compared to the Stokes layer thickness δ w = 2ν/Ω. Contribution from the vertical walls at y = 0 and y = D First, we consider the vertical wall at y = 0 with z 0. The fluid velocities in the viscous layer are given by

sin kx e -y/δw cos(Ωty/δ w )cos Ωt ,

cos kx e -y/δw cos(Ωty/δ w )cos Ωt .

(A 6) So the velocity is zero at the wall and we recover the potential solution u 0 away from it.

Cavelier et al.

The wall contribution to the mean dissipation of mechanical energy can be expressed by

Taking the square of the y-derivative, we get at leading order in 1/δ w ,

Therefore, the full contribution to the mean mechanical energy dissipation of the two vertical walls along x, summing the bottom and the upper part, leads to Ėmech = 4 Ė1 . So the damping coefficient γ w1 corresponding to the two lateral walls at y = 0 and y = D is simply given by

Vertical walls x = 0 and x = W Similarly, near the vertical wall at x = 0 and z 0, the velocity in the viscous layer is given by

The amount of energy dissipated on this wall is

Again considering the upper part of the tank and the opposite wall at x = W we obtain a dissipation of 4 Ė2 . This gives the following damping coefficient

Bottom and upper walls at z = -H/2, H/2 Considering the bottom wall at z = -H/2, the fluid velocities in the viscous layer are

At leading order in 1/δ w we obtain the energy loss as

The subcritical transition to turbulence

33

The same amount of energy is dissipated on the upper walls z = H/2. So the damping coefficient accounting for the bottom and the upper wall is

Summing all the contributions from the wall we finally obtain the damping coefficient γ w = γ w1 + γ w2 + γ w3 .

Appendix B. Stability analysis of the global approach equation

In this part, we detail the Floquet analysis of (4.3). After choosing

with τ = ωt and Ω 2 0B = 2AG 0 /L 0 . The Floquet theorem states that the solution of (B 1) is on the form

with λ the real Floquet exponent characterizing the growth rate of the instability and α = 0 or 1/2 depending if the instability is harmonic or subharmonic. By setting λ = 0, we can determine the neutral branches of the instability. Inserting the solution in (B 1), we get

Y n e i(n+α)τ = 0, (B 3)

Reorganizing the sums, we obtain for each n

This constitutes a generalized eigenvalues problem of the form AX = sBX where X is constructed from the real and imaginary parts of the vector Y . Following [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] to express the condition for c being real and truncating the solution, we restrict the problem to 0 n N leading to (2N + 2) × (2N + 2) matrix sizes for A and B. Furthermore, we focus on the subharmonic instability corresponding to α = 1/2. In this case, the pentadiagonal matrices A and B are written as

Appendix C. Derivation of the local approach equations

In this part, we detail the derivation of (4.5), first recalling the flow field generated by the primary wave and then the equations for a perturbation at the node. By combining the previous conditions, we obtain the equation already derived by [START_REF] Kelly | The stability of an unsteady Kelvin-Helmholtz flow[END_REF] and here expressed in the Boussinesq limit η -2iAk wb U η + AGk wb -k 2 wb U 2 -iAk wb U η = 0. (C 13)

It should be stressed that this equation is derived in the context of a small primary wave amplitude, kη 1, and small perturbation wavelength, κ 1, in order to consider the primary wave solution corresponding to (C 7).