2.1. Experimental set-up and parameters

The experimental set-up is now introduced briefly since the details can be found in Briard et al. (Reference Briard, Gostiaux and Gréa2020). We fill a cuboidal tank of inner length $W=94.6 \ \text {cm}$, width $D=11 \ \text {cm}$, and height $H=67 \ \text {cm}$, with salt and fresh water (see figure 1). The salt water density takes the values $\rho _1=1030$, $1060$ or $1090 \ \text {kg} \ \text {m}^{-3}$, while for the fresh water we get $\rho _2=998 \ \text {kg} \ \text {m}^{-3}$. This corresponds to various Atwood numbers expressing the density contrast: $\mathcal {A}=(\rho _1-\rho _2)/(\rho _1+\rho _2) \in \{0.015, 0.03, 0.045 \}$. The heavier salt water layer is initially placed at the bottom. It is separated from the lighter fresh water by a thin diffuse interface of thickness $\delta =0.5 - 1.5 \ \text {cm}$ located at half the height of the tank. This thickness may vary due to the filling procedure of the tank. The values of $\delta$ 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.

A 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 \omega t)$. Here $G_0=9.81\ \text {m}\ \text {s}^{-2}$ is the usual gravitational acceleration, $\omega$ is the frequency, and $F$ is the forcing parameter. This forcing parameter is related to the vertical displacement amplitude of the hexapod, $a_h$, as $F= a_h \omega ^{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.

We select in Briard et al. (Reference Briard, Gostiaux and Gréa2020) 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 $\mathcal {A}$.

Table 1. Label (series and number), Atwood number, forcing parameter and frequency considered for the experiments in this work. The wavenumbers and mode types corresponding to the primary Faraday wave are also indicated. The initial interface thickness $\delta$ is measured either by a probe when available or directly from the camera (labelled with $*$).

The primary standing Faraday waves observed in these experiments are characterized by a horizontal wavenumber, $k_{m,n}=\sqrt {k_x^{2}+k_y^{2}}$, associated with the mode index in the $x$ and $y$ directions, respectively, $m=k_x W /{\rm \pi}$ and $n=k_y D /{\rm \pi}$. 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.

2.2. 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 ${{\boldsymbol {U}}} ({{\boldsymbol {x}}},t)$ and the concentration of heavy fluid $C({{\boldsymbol {x}}},t)$. Here for miscible fluids with small Atwood number, the dimensionless concentration $C({{\boldsymbol {x}}},t) \in [0 \ 1]$ is related to the density as $\rho ({{\boldsymbol {x}}},t)= \rho _2+(\rho _1-\rho _2)\,C({{\boldsymbol {x}}},t)$. In the reference frame attached to the container, this leads to the classical system of equations

(2.1a)$$\begin{gather} \partial _t {{\boldsymbol{U}}}+{{\boldsymbol{U}}} \boldsymbol{\cdot} \boldsymbol{\nabla} {{\boldsymbol{U}}}={-}\boldsymbol{\nabla} \varPi -2 \mathcal{A}\,G (t)\,C {\boldsymbol{e}}_z + \nu\,\nabla ^{2} {{\boldsymbol{U}}}, \end{gather}$$

(2.1b)$$\begin{gather}\partial _t C +{{\boldsymbol{U}}} \boldsymbol{\cdot} \boldsymbol{\nabla} C=\mathcal{D}\,\nabla ^{2} C, \end{gather}$$

(2.1c)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} {{\boldsymbol{U}}}=0. \end{gather}$$

In (2.1a)–(2.1c), $\varPi$ refers to a reduced pressure, and $\nu, \mathcal {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 $2 W$) using the code already described in Briard et al. (Reference Briard, Gréa and Gostiaux2019,  Reference Briard, Gostiaux and Gréa2020). Therefore we do not seek to reproduce 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-thirds rule dealiasing. The time advancement is realized through a third-order low-storage strong-stability-preserving Runge–Kutta scheme, with implicit viscous terms. All the simulations use a $1024^{3}$ grid box with a pencil decomposition on either 1024 or 2048 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 described extensively in appendix B of Briard et al. (Reference Briard, Gostiaux and Gréa2020). 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 $\nu =\mathcal {D}=2.26 \times 10^{-6}\ \text {m}^{2}\ \text {s}^{-1}$. This corresponds to roughly 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 us from properly capturing at least the first stages of secondary instabilities developing on the Faraday wave.

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, DNSd and DNSe correspond to the wavebreaking detection. The series DNSb and DNSc are dedicated to the competition between modes $(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 $\epsilon$ of the interface perturbation and the $y$-spanwise perturbation amplitude $\epsilon _1$ for DNSf, together with the interface thicknesses $\delta$, 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 DNS have a $1024^{3}$ grid resolution.

We now detail the initial conditions taken in the simulations. While the initial velocity is ${{\boldsymbol {U}}}=0$ in the simulations, the initial concentration profile is taken as two-dimensional (2-D) of the form (for the DNSa, DNSd and DNSe series and outside the penalization band)

(2.2)\begin{equation} C(x,z)=\frac{1}{2} \left ( 1+\text{tanh} \left[ \frac{z-\xi(x)}{\sigma}\right]\right), \quad \text{with} \ \xi(x)=\epsilon \sin (k_{m,0} x). \end{equation}

The parameter $\sigma$ in (2.2) sets the initial width $\delta$ of the interface ($\delta \approx 3 \sigma$). The function $\xi (x)$ indicates the initial perturbed interface position of sinusoidal shape, with wavelength $k$ and of small amplitude $\epsilon =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 $20$ grid points across the interface layer and to ensure grid convergence of the simulations.

Without ambiguity, the initial wavenumber $k_{m,0}$ for DNSa, DNSd and DNSe 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 as nearly twice the value of the dispersion relationship of an inviscid interface $\omega =2 \sqrt {\mathcal {A} G_0 k_{m,0}}$.

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

(2.3)\begin{equation} \xi(x) = \epsilon [r \cos(k_{3,0} x) + (1-r) \cos(k_{2,0} x)]. \end{equation}

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 $\delta$ are doubled to keep at least 10 grid points across the interface, while the viscosity and diffusion coefficients are still multiplied by $4$ to ensure grid convergence of these simulations.

In the simulation series DNSa–DNSe, 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

(2.4)\begin{equation} \xi(x,y)=\epsilon \sin (k_{m,0} x)+\epsilon_1\,f (y). \end{equation}

In practice, the $y$ disturbance amplitude is set such that $\epsilon _1=10^{-2} \epsilon$. However, we have also tested various simulations varying the $\epsilon _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 Briard et al. (Reference Briard, Gostiaux and Gréa2020). These simulations (not presented) give results similar to those for DNSf, which will be discussed in § 5.3.3.

3.1. Linear theory

It is well-known since Benjamin & Ursell (Reference Benjamin and Ursell1954) that when modelling the Faraday instability, the amplitudes, $\eta _k$, for the interface modes of wavenumber $k$ are ruled by a Mathieu equation

(3.1)\begin{equation} \ddot{\eta}_k +2\,\gamma(k)\,\dot{\eta}_k +\varOmega ^{2}(k)\,(1+ F \cos \omega t) \eta_k=0. \end{equation}

In (3.1), we define the inviscid frequency of the diffuse interface $\varOmega$ and the viscous damping term $\gamma$, both of which depend on the horizontal wavenumber $k$. Note that the decoupling of each inviscid mode is true only in the limit of small damping. However, the full analysis of this problem can be performed using the method proposed by Kumar & Tuckerman (Reference Kumar and Tuckerman1994).

The inviscid frequency $\varOmega (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, Briard et al. (Reference Briard, Gostiaux and Gréa2020) for a piecewise linear profile

(3.2)\begin{equation} \varOmega (k)=\left( \frac{\mathcal{A} G_0 k}{1+ k \delta/2}\right)^{1/2}. \end{equation}

For small wavenumbers, i.e. $k \delta \ll 1$ with $\delta$ the thickness of the interface, the classical dispersion relationship for an interface within the deep-water approximation $\varOmega =\sqrt {\mathcal {A} G_0 k}$ is recovered. In the large wavenumber limit, $k \delta \gg 1$, the interface mode reduces to $\varOmega (k)=\sqrt {2 \mathcal {A} G_0 /\delta }$, corresponding to the local buoyancy or Brunt–Väisälä frequency at the interface.

The viscous dissipation term, $\gamma (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 $\gamma _b(k)=2 \nu k^{2}$ (Lamb Reference Lamb1945; Landau & Lifshitz Reference Landau and Lifshitz2013). However, due to the velocity gradients, significant damping can also occur within the thin layer separating the two fluids. Assuming a piecewise linear vertical density profile, Briard et al. (Reference Briard, Gostiaux and Gréa2020) obtained the expression $\gamma _\delta (k)=\mathcal {A} G_0 \nu k^{2}/ \varOmega ^{2} \delta \approx \nu k/\delta$ for $k \delta \ll 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 $\delta _w=(2 \nu / \varOmega )^{1/2}$. This gives values $\delta _w \sim 1 - 2 \ \text {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, Keulegan (Reference Keulegan1959) and Miles & Benjamin (Reference Miles and Benjamin1967) 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 Thorpe (Reference Thorpe1968) as detailed in Appendix B, and leads to the following expression for the damping coefficient:

(3.3)\begin{equation} \gamma _{w}\approx \frac{\nu}{D \delta_w}= \frac{\sqrt{\nu \varOmega}}{\sqrt{2}D}.\end{equation}

Here, (3.3) thus expresses the dominant contribution of the lateral walls (in the $z$–$x$ plane) 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 the boundary layers at the walls. It can be shown that these values do not vary by more than several per cent 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 $\gamma _w$ is larger than the other contributions. The damping $\gamma _{w}$ indeed scales like $\nu ^{1/2}$ in (3.3), while it is linear in $\nu$ for the dissipation $\gamma _b$ or from the interfacial layer $\gamma _\delta$. Bechhoefer et al. (Reference Bechhoefer, Ego, Manneville and Johnson1995) have extensively discussed this aspect and they suggest using 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 favour the use of low-viscosity fluids. Note that for larger wavenumber $k$, the contributions from the interface layer regain importance and cannot be neglected.

Table 3. Values for the damping coefficients, $\gamma _w, \gamma _b, \gamma _\delta$, in $\text {s}^{-1}$ and evaluated for the largest wavelengths developing in the experiment. We assume here that the Atwood number is $\mathcal {A} =0.03$ and the thickness of the interfacial layer is $\delta =0.5 \ \text {cm}$. Here, the top boundary is taken as a wall to evaluate $\gamma_w$ (the values would be nearly the same for a free surface).

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 $\gamma _w$ from (3.3) and $\gamma _\delta$ (the latter contribution being smaller). The neutral curves of (3.1) are computed using the method proposed by Kumar & Tuckerman (Reference Kumar and Tuckerman1994) and also used in Briard et al. (Reference Briard, Gostiaux and Gréa2020) assuming different Atwood number values $\mathcal {A}$ and an initial interface thickness $\delta =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, $\varOmega (k) =\omega /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 \gamma /\omega$ (see Rajchenbach & Clamond (Reference Rajchenbach and Clamond2015), 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 $\gamma _w \sim \omega ^{1/2}$, leading to a decrease of $F_{th}$ at larger $\omega$. However, this effect is compensated at larger $k$ by the contribution from the damping at the interface scaling like $\gamma _\delta \sim \omega ^{2}$.

Figure 2. Stability diagram for (3.1) in a non-dimensional frequency $\omega /\sqrt {\mathcal {A} G_0 k_{2,0}}$ and forcing $F$ plane. The coloured 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 $\gamma =\gamma _\delta +\gamma _w$ at three different Atwood numbers $\mathcal {A}$, and considering an interface thickness $\delta =1 \ \text {cm}$. The neutral curves (thick plain lines) have a slight dependence on the Atwood number, explaining why 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 colours reveal which mode is eventually selected.

The parameters taken in the experiments with $F\geqslant 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.

3.2. Linear or nonlinear mode selection?

In this subsection, 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 favour the modes pertaining to the right unstable tongues in figure 2 (the mode reported in table 1 is also indicated by the colour of the symbol in figure 2). In addition, the selection mechanism does not apparently discriminate between the even and 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 centre 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. Moreover, 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 (${\pm }0.001$) or initial interface width (${\pm }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 slightly left-shift the instability tongues of figure 2 as the natural frequencies $\varOmega$ are decreased (the damping dominated by the viscous layer at the wall remains unchanged). In any case, 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 that 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 does not have 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 a disturbance and could not identify a source able to generate it. In any case, 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, Ciliberto & Gollub Reference Ciliberto and Gollub1984,  Reference Ciliberto and Gollub1985). Using weakly nonlinear approaches, Meron & Procaccia (Reference Meron and Procaccia1986) and Meron (Reference Meron1987) 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 $\varDelta ={\varOmega ^{2}}/{\omega ^{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 $\varDelta$ generally have opposite signs because the forcing frequency $\omega$ 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, $\varDelta < 0$, is supercritical as the nonlinear coupling damps the instability while the dominant one, $\varDelta \geqslant 0$, is subcritical as it is being reinforced by the nonlinearities. As the wave amplitude grows, the frequency of the wave tends to diminish (Thorpe Reference Thorpe1968) as well as the detuning parameter of each mode (Godrèche & Manneville Reference Godrèche and Manneville2005). Hence this left-shifts the instability bands of figure 2 and favours the subcritical modes at smaller wavelength.

3.3. 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. In any case, 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 DNSc in table 2), with $\mathcal {A} =0.03$. The frequency $\omega$ 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 $\omega =3.07$ (or equivalently $\omega / \sqrt { \mathcal {A} G_0 k_{2,0}}=2.2$). This corresponds to parameters with the two unstable modes having nearly the same exponential growth rate as $\sim {\rm e}^{\lambda \omega t}$. Indeed, the Floquet exponent $\lambda$ 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 $\omega$ 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)$ and $(3,0)$ modes prevails at this moment in figure 3. This procedure is performed both visually and by computing the Fourier modes of the interface.

Figure 3. Parameters of the DNS (symbols) in the stability diagram $\omega$–$F$. The coloured areas correspond to the subharmonic instability tongues for the modes $(2,0)$ and $(3,0)$ accounting for viscosity and diffused interface $\delta =1.8 \ \text {cm}$ (as no walls are present in the DNS, the damping coefficient $\gamma \approx \gamma _\delta$ and critical thresholds $F_{th}$ are very small). The symbols’ colours indicate which mode emerges from the simulation; mixed colours express that both modes can be observed. Two series of DNSb, DNSc (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 $\omega$. In the DNSc group, the frequency and forcing are fixed and we vary $r$. The point corresponding to DNSe is also placed.

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 $\omega$). The phenomenon can be scrutinized in more detail in the snapshots extracted from the two series in figure 4, where the mode $(3,0)$ emerges from cases with initial $r=0.25, \ 0.5$ or with the frequencies $\omega =2.6, \ 2.8$. Indeed, in the last row of figure 4, one can observe that mode $(2,0)$ prevails only for $\omega =2.4$ ($\omega /\sqrt {\mathcal {A} G_0 k_{2,0}}=1.72$), and for larger $\omega$, say $\omega =2.6$ (or $\omega /\sqrt {\mathcal {A} G_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).

Figure 4. Mode selection in six DNS of figure 3. (a) Cases corresponding to DNSc (see table 2) with $r$ varying and $\mathcal {A}=0.030$, $\omega =3.07$, $F=0.8$. The amplitude of the interface deformation is $\epsilon =3 \ \text {cm}$. (b) Cases corresponding to DNSb (see table 2) with $\omega$ varying and $\mathcal {A}=0.030$, $F=0.8$, $r=0.5$, $\epsilon =1.5 \ \text {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 \pm 0.15$ is in black.

Some specificities of the mode suppression in our Faraday experiment with miscible fluids are now addressed. We have not observed oscillations between two specific modes as in Ciliberto & Gollub (Reference Ciliberto and Gollub1984,  Reference Ciliberto and Gollub1985) or similarly Yalim, Welfert & Lopez (Reference Yalim, Welfert and Lopez2019) in the context of a stable stratification. This is notably because, at large forcing parameter, the interface grows irreversibly, allowing continuously new modes to be destabilized. The modes can change in our experiment, as already reported in Briard et al. (Reference Briard, Gostiaux and Gréa2020). But there is always a correspondence to a one-way transition from large to small wavelength for interface modes. The more complex transitions evidenced in figure 14 of Briard et al. (Reference Briard, Gostiaux and Gréa2020), for instance, refer to modes not pertaining to the same instability band or being of different natures. 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 approximation, the subcritical modes are indeed those with small wavelength; see Rajchenbach & Clamond (Reference Rajchenbach and Clamond2015) for details.

In this section, 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.

4.1. The global approach

4.1.1. 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=\bar {C} +c$ and ${\boldsymbol {U}}=\bar {{\boldsymbol {U}}} +{\boldsymbol {u}}$. A mean quantity $\bar {*}$ is obtained by averaging along the horizontal $x$ and $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 shown directly that the mean velocity is zero, $\bar {\boldsymbol {U}}=0$, due to the symmetries and the incompressibility condition. The mean vertical concentration profile $\bar {C} (z,t)$ evolves principally due to the vertical buoyancy flux $\overline {w c}$ as

(4.1)\begin{equation} \partial_t \bar{C}+ \partial_z \overline{w c}=\mathcal{D} \partial^{2} _{zz} \bar{C}.\end{equation}

In this global approach, the primary Faraday wave is thus embedded in the mean vertical concentration profile $\bar {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 developed by Hunt & Carruthers (Reference Hunt and Carruthers1990) and applied to turbulent mixing layers in Gréa (Reference Gréa2013) provides a convenient framework for expressing the dynamics of small-scale disturbances. We thus discard all the dissipative and nonlinear terms, except for the coupling between the fluctuations and the mean concentration field. In addition, the small-scale disturbance sees only a uniform mean concentration gradient at $z=0$ determined by the length $L=-1/\partial _z \bar {C}(0,t)$. Within this quasi-homogeneous approximation, the small-scale fluctuating quantities and pressure $p$ are determined by

(4.2a)$$\begin{gather} \partial _t {\boldsymbol{u}} ={-}\boldsymbol{\nabla} p -2 \mathcal{A}\,G(t)\,c {\boldsymbol{e}}_z, \end{gather}$$

(4.2b)$$\begin{gather}\partial _t c =\frac{w}{L(t)}, \end{gather}$$

(4.2c)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}}=0. \end{gather}$$

One recognizes in (4.2) the classical equations for an internal gravity wave, except that 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 $L_{int}= 6 \int _{-\infty } ^{+\infty } \bar {C} (1-\bar {C})\,{\rm d}z$ introduced by Andrews & Spalding (Reference Andrews and Spalding1990) and previously used in Briard et al. (Reference Briard, Gostiaux and Gréa2020) 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 later section different ways to evaluate it.

These waves depend 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 $\overline {w c}$, 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 Gréa & Ebo Adou Reference Gréa and Ebo Adou2018):

(4.3)\begin{equation} \ddot{c} + \frac{\dot{L}}{L}\,\dot{c} +\frac{2 \mathcal{A}\,G(t)}{L}\,c=0. \end{equation}

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, (4.3) exhibits the buoyancy frequency defined as $\varOmega _B=(2 \mathcal {A} G_0/L)^{1/2}$ and a damping term $\dot {L}/L$ expressing the variations of $\varOmega _B$ as the mixing zone width $L$ evolves. As will be seen in § 4.2, the $\varOmega _B$ frequency is relevant for the secondary instability because the shear at the nodes of the primary wave is driven directly by the wave amplitude.

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

4.1.2. The subcritical nature and the criterion for the wavebreaking

The stability of the model equation (4.3) has been extensively discussed in Gréa & Ebo Adou (Reference Gréa and Ebo Adou2018) and Briard et al. (Reference Briard, Gostiaux and Gréa2020), 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 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 $\varOmega _B$, we thus aim to find 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 develops rapidly at the interface, and sheds light on the subcritical nature of this secondary instability. This peculiarity is indeed well-known for nonlinear Mathieu systems such as (4.3) as detailed in Soliman & Thompson (Reference Soliman and Thompson1992) or Godrèche & Manneville (Reference Godrèche and Manneville2005).

In order to derive an analytic criterion for the wavebreaking, we consider the inverse mean concentration gradient $L$ having the simple form $L(t)=L_0 (1+\beta \cos (\omega t))$. The length $L$ is thus expected to be proportional to the amplitude $\eta _p$ of the primary Faraday wave, $L(t) \sim |\eta _p |$ in the laminar phase. Here, the proportionality coefficient depends on the shape of the nonlinear primary wave. Also, the parameter $\beta$ 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 frequency $\omega$ while $\eta _p$ oscillates at frequency $\omega /2$. However, the higher temporal harmonics of $L(t)$ or $\eta _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 $\beta =0$ and $\beta =0.7$. Note that the case $\beta =0$ indeed corresponds 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}$ that should be close to the one characterizing the wavebreaking $L_{wb}$ if the instability develops quickly (we still have $L_{crit}\leqslant L_{wb}$).

Figure 6. Stability analysis of (4.3) with $L(t)=L_0 (1+\beta \cos (\omega t))$ and represented in the $\varOmega _{0B}^{2}/\omega ^{2}$–$F$ plane (with $\varOmega _{0B}=(2 \mathcal {A} G_0/L_0)^{1/2}$). The two coloured areas indicate the first subharmonic tongues with $\beta =0$ and $\beta =0.7$, respectively. The dashed blue lines correspond to the approximation (4.4), while the red dashed dotted lines correspond to (B9).

An analytical approximation for the critical threshold $L_{crit}$ is also derived in Appendix B, leading to

(4.4)\begin{equation} L_{crit}\approx\frac{(4-2 F)}{1 +\beta/2}\,\frac{2 \mathcal{A} G_0 }{\omega ^{2}}, \quad \text{for} \ F \ll 1 \ \text{and} \ \beta \ll 1. \end{equation}

As shown in figure 6 (blue dashed lines), the criterion (4.4) slightly underestimates $L_\text {crit}$ at small $F$ and large $\beta$ 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 \geqslant 2$. Despite being more complicated, a better approximation can be derived corresponding to (B9) 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 $\beta$ 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 $\lambda =0.09$ for $\beta =0$ (also corresponding to the growth of the primary wave), while $\lambda =0.2$ for the case with $\beta =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.

4.2. 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 losing track of the physical mechanism responsible for the secondary instability. Also, it assumes that the secondary instability mode results only 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, $\eta$, at the node of the primary wave, are given in the inviscid limit (see the derivation in Appendix C) by

(4.5)\begin{equation} \ddot{\eta} -2 {\rm i} \mathcal{A} {k_{wb}}U \dot{\eta}+( \mathcal{A}\,G (t)\,{k_{wb}} - {k_{wb}} ^{2} U^{2}-{\rm i} \mathcal{A} {k_{wb}} \dot{U} ) \eta=0. \end{equation}

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 Kelly (Reference Kelly1965), expressing the dynamics of an oscillating sheared interface. In fact, supposing the instability amplitude and wavelength are small, respectively $k \eta \ll 1$ and $\kappa ={k_{wb}}/k \gg 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 $\omega$. In addition to the parametric resonance modes that Kelly (Reference Kelly1965) identified, Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1987) and Khenner et al. (Reference Khenner, Lyubimov, Belozerova and Roux1999) 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, which may exhibit frozen wave patterns at high forcing frequency (Wolf Reference Wolf1970; Wunenburger et al. Reference Wunenburger, Evesque, Chabot, Garrabos, Fauve and Beysens1999; Yoshikawa & Wesfreid Reference Yoshikawa and Wesfreid2011; Gaponenko et al. Reference Gaponenko, Torregrosa, Yasnou, Mialdun and Shevtsova2015; Lyubimov et al. Reference Lyubimov, Khilko, Ivantsov and Lyubimova2017; Gréa & Briard Reference Gréa and Briard2019). 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 Thorpe (Reference Thorpe1968) and Kalinichenko (Reference Kalinichenko2005) 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 Kelly (Reference Kelly1965) and Khenner et al. (Reference Khenner, Lyubimov, Belozerova and Roux1999) do not apply to our specific configuration where the acceleration oscillates at frequency $\omega$ while the shear velocity is subharmonic with frequency $\omega /2$. We thus reconsider the problem of (4.5) by taking a primary wave of the form $\eta _p(t)=a \cos (\omega t/2)$ leading to

(4.6a,b)\begin{align} U(t)={-} \frac{\omega a}{2} \sin (\omega t /2)\quad\text{and} \quad U^{2}=\frac{\omega^{2} a^{2}}{8}\,(1-\cos \omega t)= \frac{1}{2}\,\frac{\mathcal{A} G_0 k}{1+ 4 \varDelta}\,a^{2}(1-\cos \omega t), \end{align}

where in the last expression, the subharmonic resonance condition $\omega ^{2}=4 \mathcal {A} G_0 k/ (1+4 \varDelta )$ 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 $\eta = Y \exp ({\int _0^{t} {\rm i} \mathcal {A}{k_{wb}}\,U(t') \,{\rm d} t'})$ giving (at leading order in $\mathcal {A}$)

(4.7)\begin{equation} \ddot{Y}+( \mathcal{A}\,G(t)\,{k_{wb}} - {k_{wb}} ^{2} U^{2} )Y=0. \end{equation}

Due to the change of variable expression and the fact that $U$ oscillates at $\omega /2$, the response in $\eta$ will be also subharmonic. However, for small Atwood number, $\eta \approx 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 of the form

(4.8)\begin{gather} \left.\begin{gathered} \ddot{Y}+\left( P +Q \cos (\omega t) \right) Y =0, \\ \text{with} \ P=\mathcal{A} G_0 k \left( \kappa -\frac{1}{2(1+ 4 \varDelta)}\,\kappa^{2} (k a)^{2} \right) \ \text{and} \\ Q=\mathcal{A} G_0 k \left( \kappa F+ \frac{1}{2 (1+ 4 \varDelta)}\,\kappa^{2} (k a)^{2} \right). \end{gathered}\right\} \end{gather}

In figure 7, we show the stability diagram of (4.8) in a $\kappa$–$k a$ plane (the subharmonic resonance condition for the primary wave is used again). The instability tongue corresponding to the KH type modes appears for $P \leqslant 0$, which stands as the classical criterion for the inviscid KH instability, ${k_{wb}} \overline { U^{2}} \geqslant \mathcal {A} G_0$ (Chandrasekhar Reference Chandrasekhar1961). The parametric resonance zones start for $P \geqslant 0$ but also have a continuation in the opposite half-plane. Remarkably, the instability zones exhibit a very weak dependence on $\kappa$ for $\kappa \gg 1$. Therefore, at given $\kappa$ and as the primary wave amplitude grows, the perturbation passes through the successive instability zones, first the PR types then later 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, $k a$, lies in the instability KH or PR1 zones. This local theory is inviscid, which explains why the growth rates are higher at large $\kappa$. Of course, the viscosity and the thickness $\delta$ of the interface should play a role and moderate this aspect.

Figure 7. Stability curve of (4.8) for $F=0.7$ and $\varDelta =0$ in the $\kappa$–$k a$ plane or the $P$–$Q$ plane (insert) corresponding to the classical representation of the Mathieu equation. The blue coloured 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 $\kappa$–$k a$ representation. The critical wave steepness value indicated by the black dotted curve corresponds to criterion (4.9).

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 \rightarrow -\infty$, $Q \rightarrow +\infty$ detailed in Abramowitz & Stegun (Reference Abramowitz and Stegun1965), we obtain

(4.9)\begin{equation} k a _{crit}\approx \tfrac{1}{3} (1+4 \varDelta) (1+F). \end{equation}

This simple expression (4.9) corresponds indeed to the plateau (it does not depend on $\kappa$) separating the PR1 and PR2 bands in the small perturbation wavelength limit $\kappa \gg 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 the simulations.

5.1. Primary wave amplitude and inverse mean density gradient

The global and local approaches detailed in § 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 $\eta _p(t)$ (see also figure 5). It is therefore important to link these two quantities in order to compare the global and local theories in the experiments. As already remarked, we can expect the proportionality between the inverse mean gradient and the wave amplitude, $L(t) = \kappa_L |\eta _p (t)|$. This assumes that the shape of the primary wave is frozen, which also determines the coefficient $\kappa$. By considering the position $\xi _p(x,t)=\eta _p (t) \sin k x$ of a sharp sinusoidal interface, as for the standing wave in Appendix C, we can easily show that $\kappa_L = {\rm \pi}$ (the wave amplitude is supposed larger than the interface thickness $\delta$). This choice is also motivated by the expression for the finite amplitude standing wave profile given by Thorpe (Reference Thorpe1968) in deep-water approximation, 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+ \beta \cos \omega t)$, while in the local theory the primary wave amplitude evolves as $\eta _p= a \cos \omega t /2$. By expanding $|\eta _p (t)|$ in Fourier series and truncating at leading order, we simply get $L_0=2 \kappa_L a/{\rm \pi}$ and $\beta =2/3$. This allows us 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 $\kappa_L ={\rm \pi}$ (sinusoidal interface), gives

(5.1)\begin{equation} k a_{crit}=\tfrac{3}{8} (2-F)(1+4 \varDelta). \end{equation}

We now detail how to measure the inverse local gradient $L$ and the wave amplitude $\eta _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 post-processing of the images of the camera as detailed in Briard et al. (Reference Briard, Gostiaux and Gréa2020), can be a bit noisy, particularly when the secondary instabilities start at the node. The quantity $|\eta _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 \leqslant \bar {C} \leqslant 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 § 4.1.1. By assuming that the shape of the concentration profiles is frozen, we can deduce all the characteristic lengths from $L_{int}$, which for a sinusoidal interface gives $L_{int} = 2.4|\eta _p|=0.76 L$.

To ensure that these relations apply in our problem, we have measured the amplitude $|\eta _p|$ at a local maximum just before the wavebreaking using the threshold method on $\bar {C}$, the inverse mean concentration gradient $L$, and the integral length $L_{int}$. In figures 8(a) and 8(b), we see that the sinusoidal profile is a good fit for the interface position and provides a good evaluation of the mean concentration profile. If the amplitude is too small, however, 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 centre of the layer $z=0$, in accordance with the fact that the wavebreaking starts at the nodes of the primary wave. In any case, we recover the expected relation between $L$ and the wave amplitude $|\eta _p|$. Therefore the different lengths $L$ and $|\eta _p|$ can be evaluated correctly from $L_{int}$. In figure 8(c), we further check that the relation between $L_{int}$ and $|\eta _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.

Figure 8. (a) Visualization 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) $L_{int}$ plots as a function of $|\eta _p|$ at a maximum amplitude just before the wavebreaking for all the experiments of table 1. The values for $|\eta _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,b).

5.2. Wavebreaking detection

In this subsection, 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 observed directly 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 Thorpe Reference Thorpe1977). Indeed, this displacement characterizes the distance that 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 for each one-dimensional vertical transect the displacements $\delta _{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 $|\delta _{T}|$ exceeds a certain threshold (here $|\delta _{T}| > 5 - 12 \textrm {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 eye only (see figure 9). However, it is not possible to detect from this method the beginning of the ‘blurred’ region as defined in Thorpe (Reference Thorpe1968), indicating the secondary instability onset leading to the wavebreaking at the nodes.

Figure 9. Procedure for the wavebreaking detection. (a) Visual criterion from the calibrated camera image. (b) The Thorpe displacements $\delta _T$ evaluated by sorting the concentration field in each vertical transect of the calibrated image. The colourbar indicates the displacement in mm.

The Thorpe displacement can also be used to measure the vortex size at wavebreaking and thus gives the ratio $\kappa = {k_{wb}}/ k$ necessary for the local theory. To determine ${k_{wb}}$, we take the maximum displacement $\delta _{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, $\delta _{T}$ exceeds this value by more than three times. 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.

Figure 10. Time evolution of the mixing zone width and the mean concentration profiles at different instants extracted from DNSd2 (see parameters in table 2). (a) Evolution of the integral lengths $L_{int}$ and $L_{int,s}$ computed from the mean and sorted concentration profiles, respectively, as a function of $\omega t$. The green dotted line corresponds to the integral scale $L_{int,d}$ expressing the thickening of the interface by diffusion only. The star symbol at $\omega 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 illustrate the state of the interface at the same instants.

The time evolution of the integral length $L_{int}$ in a simulation revealing the subharmonic oscillations (thus $L_{int}$ oscillates at frequency $\omega$) and the growth of the primary wave are presented in figure 10(a). 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 from the sorted concentration profiles as in Briard et al. (Reference Briard, Gréa and Gostiaux2019). This quantity expresses the evolution of irreversible mixing by distinguishing the available and background potential energies (see also, for instance, Winters et al. Reference Winters, Lombard, Riley and D'Asaro1995; Peltier & Caulfield Reference Peltier and Caulfield2003; Davies Wykes & Dalziel Reference Davies Wykes and Dalziel2014). Noticeably, we observe that the growth of the irreversible mixing starts just after the wavebreaking time. In figure 10(a), 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 10(b), 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 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 10(a), 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 this point more thoroughly in the next subsection.