2.1. General equations

We consider a system of two immiscible, incompressible fluids, each with density $\rho _i$ , where the subscript $i=1$ identifies the lighter fluid and $i=2$ the denser fluid ( $\rho _2 \gt \rho _1$ ), as represented in figure 1. The fluids are confined within a cuboidal container that is oscillated at frequency $\omega$ and amplitude $a$ along one of its principal axes, specifically the horizontal $x$ -direction. In a Cartesian coordinate frame $\boldsymbol {x}=(x,y,z)$ moving with the container, the system is subjected to the ambient gravitational field $g$ , which acts in the vertical $z$ -direction, and the oscillatory acceleration in the horizontal $x$ -direction.

Figure 1. Schematic representation of the experimental container and the frozen wave instability studied in this work, accompanied by a glossary of the key non-dimensional numbers governing the system’s dynamics.

This system is governed by the incompressible Navier–Stokes equations, which describe the dynamics of the velocity field ${\boldsymbol {U}}(\boldsymbol {x},t)$ , the mass fraction of the heavier fluid $Y(\boldsymbol {x},t)$ and the pressure $P(\boldsymbol {x},t)$ . These variables, along with the spatial coordinate $\boldsymbol {x}$ and time $t$ , are non-dimensionalised using the characteristic length $a$ , oscillation frequency $\omega$ and a reference density $\rho _0=(\rho _1+\rho _2)/2$ . This yields

(2.1a) \begin{align}& \partial _t {\boldsymbol {U}}+ {\boldsymbol {U}} \boldsymbol{\cdot }(\boldsymbol{\nabla }{\boldsymbol {U}})=-\frac {1}{\rho (Y)} \boldsymbol{\nabla }P - \frac {1}{F}\boldsymbol {e}_z+ \cos ( t) \boldsymbol {e}_x+\frac {1}{{\textit{Re}}}\boldsymbol{\nabla }\boldsymbol{\cdot }\varPi _\nu +\frac {1}{\rho (Y)} \boldsymbol{F}_\gamma , \end{align}

(2.1b) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \partial _t Y+ {\boldsymbol {U}} \boldsymbol{\cdot }\boldsymbol{\nabla }Y = 0, \end{align}

(2.1c) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\, \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol {U}}=0, \end{align}

where (2.1a ) expresses the conservation of momentum, (2.1b ) the conservation of mass and (2.1c ) the incompressibility condition. The renormalised density $\rho$ is related to the mass fraction of heavy fluid $Y$ and the Atwood number $\mathcal A=(\rho _2-\rho _1)/(\rho _1+\rho _2)$ as

(2.2) \begin{equation} \frac {1}{\rho (Y)}=\frac {Y}{1+\mathcal A} +\frac {1-Y}{1-\mathcal A}. \end{equation}

In the non-dimensionalised system, the Froude number $F=a \omega ^2/g$ (also called forcing parameter) emerges in the governing equations to quantify the ratio of inertial forces to gravitational forces. This dimensionless parameter directly reflects the relative influence of the oscillatory driving force (characterised by the amplitude and frequency) compared with the stabilising effect of gravity.

In (2.1a ), the viscous effects are characterised by the Reynolds number, ${\textit{Re}}=a^2 \omega /\nu _2$ , where $\nu _2$ is the kinematic viscosity of the denser fluid. The viscous stress tensor $\varPi _\nu$ differs between the two Newtonian fluids due to their distinct viscosities: in fluid $2$ (the denser fluid), the viscous stress tensor is given by $\varPi _\nu =\boldsymbol{\nabla }{\boldsymbol {U}}+\boldsymbol{\nabla} ^T {\boldsymbol {U}}$ . In fluid $1$ (the lighter fluid), the viscous stress tensor is modified by the viscosity ratio $\varUpsilon =\nu _1/\nu _2$ , and is expressed as $\varPi _\nu =\varUpsilon (\boldsymbol{\nabla }{\boldsymbol {U}}+\boldsymbol{\nabla} ^T {\boldsymbol {U}})$ , where $\nu _1$ is the kinematic viscosity of the lighter fluid. The viscosity ratio $\varUpsilon$ introduces asymmetry in the viscous contributions between the two fluids, influencing the dynamics at the interface and throughout the flow.

The effect of surface tension $\gamma$ between the two fluids, which acts to minimise the interfacial area, is incorporated into the system through a renormalised volumetric force, $\boldsymbol {F}\gamma =2 \mathcal A \kappa \delta _s {\boldsymbol {n}}/ {\textit{We}}$ . Here, $\kappa$ denotes the curvature of the interface, $\boldsymbol {n}$ is the unit vector normal to the interface and $\delta _s$ is the surface Dirac function localised at the interface (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011; Popinet Reference Popinet2018). The formulation involves the Weber number, ${\textit{We}} = a^3 \omega ^2 (\rho _2 - \rho _1) / \gamma$ , which quantifies the relative importance of inertial forces to surface tension forces at the interface. A high Weber number indicates that inertial effects dominate, whereas a low Weber number suggests a stabilising influence of surface tension.

The system of equations ((2.1a )–(2.1c )) is supplemented with boundary conditions at the walls of the container, where the no-slip condition is imposed, i.e. ${\boldsymbol {U}}=0$ . In addition to the five non-dimensional parameters $\{\mathcal A,F,{\textit{Re}},{\textit{We}},\varUpsilon \}$ , which describe the properties of the fluids and the vibrational forcing, the system is further influenced by three geometric parameters, which are the non-dimensional lengths of the container, $\{\mathcal H \equiv H/a, \mathcal W \equiv W/a,\mathcal D \equiv D/a \}$ where $H$ , $W$ and $D$ are the height, width and depth of the container, respectively. Assuming that the dynamics of the interface is independent of initial conditions, the system is fully characterised by this set of eight non-dimensional parameters. These parameters, which are recalled in figure 1, encapsulate the fluid properties, forcing characteristics and geometric constraints of the system, providing a comprehensive description of the interface dynamics.

2.2. Inviscid theory of Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986)

In this section, we rewrite the linear inviscid stability criterion derived by Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986) in a non-dimensional form. The instability onset of frozen waves is associated with the vibrational velocity, leading to the criterion

(2.3) \begin{equation} \mathcal C=\frac {(a \omega )^2}{(a \omega )^2_c} = \frac {\mathcal A (1-\mathcal A^2){\textit{We}}^{1/2} F^{1/2}}{4} \geqslant 1. \end{equation}

The critical number $\mathcal C$ quantifies how far the system’s parameters lie from the neutral stability curve at $\mathcal C=1$ . Moreover, near neutrality the most unstable wavelength is $2 \pi \ell _\gamma$ , introducing the non-dimensional capillary length as

(2.4) \begin{equation} \ell _\gamma =\sqrt {\frac {F}{\textit{We}}}=\frac {1}{a}\sqrt {\frac {\gamma } {(\rho _2-\rho _1)g}}. \end{equation}

Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986) further derived, through a weakly nonlinear analysis, the non-dimensional amplitude of the frozen waves, $L$ , near the instability onset when (2.3) is satisfied

(2.5) \begin{equation} L^2=\frac {16 \mathcal A (1+\mathcal A)}{(1-\mathcal A)(5-16 \mathcal A^2)} \frac {F^{3/2}}{{\textit{We}}^{1/2}} \left (1-\frac {1}{\mathcal C} \right ) \end{equation}.

This criterion underscores that the instability arises only in the presence of a non-zero density contrast (or non-zero Atwood number), which is essential for generating shear at the fluid interface. Notably, the expression in (2.5) is valid for Atwood numbers less than $\sqrt {5}/4 \approx 0.559$ , corresponding to a supercritical instability. For larger Atwood numbers, the instability transitions into a subcritical regime. Moreover, the formula indicates that the dimensional amplitude of the frozen waves scales linearly with the vibrational velocity, $a \omega$ , as the system moves further away from the onset of instability.

As mentioned in the introduction, the inviscid analysis presented by Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986) has been extended to include the effects of viscosity (Lyubimov et al. Reference Lyubimov, Khilko, Ivantsov and Lyubimova2017). This extended analysis highlights, in particular, the influence of the viscosity contrast, expressed by $\varUpsilon$ , in lowering the instability onset. This extension does not alter the weakly nonlinear results, which show that, away from the instability onset, the amplitude of the frozen waves continues to scale linearly with the vibrational velocity.

In the following sections, we propose a theoretical framework to describe the behaviour of frozen waves in the strongly nonlinear inertial regime corresponding to a large critical number $\mathcal C$ and Reynolds number ${\textit{Re}}$ . This framework extends the existing analyses to capture the complex dynamics and amplitude scaling that arise far beyond the instability onset, where inertial effects dominate the system’s response.

2.3. Favre-averaged equations and vibro-equilibrium solution

To characterise the behaviour of frozen waves in the inertial regime, where ${\textit{Re}} \gg 1$ and ${\textit{We}} \gg 1$ , we focus on the equations averaged along the horizontal plane $(x,y)$ . Any quantity $Q$ can be decomposed into its Reynolds average $\overline Q$ and its fluctuation $q'$ as follows:

(2.6) \begin{equation} Q=\overline {Q}+q'. \end{equation}

However, when dealing with two fluids of significantly different densities, it is convenient to introduce the Favre average (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2013), defined as

(2.7) \begin{equation} \widetilde Q = \frac {\overline {\rho Q}}{\overline {\rho }}. \end{equation}

Similar to the Reynolds decomposition, any field $Q$ can be expressed in terms of its Favre average and fluctuation $q''$ , such that

(2.8) \begin{equation} Q=\widetilde {Q}+q''. \end{equation}

This leads the Favre-averaged equations which take the following form after discarding the viscous and surface tension terms:

(2.9a) \begin{align} \overline { \rho } \partial _t \widetilde {{\boldsymbol {U}}}+ \overline { \rho } \widetilde {{\boldsymbol {U}}} \boldsymbol{\cdot }\left(\boldsymbol{\nabla }\widetilde {{\boldsymbol {U}}}\right)+\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\overline { \rho } \widetilde {{\boldsymbol {u}}''{\boldsymbol {u}}''}\right)&=-\overline {\boldsymbol{\nabla }P} - \frac {\overline { \rho }}{F}\boldsymbol {e}_z+ \overline { \rho } \cos ( t) \boldsymbol {e}_x , \end{align}

(2.9b) \begin{align} \overline { \rho } \partial _t \widetilde {Y}+ \overline { \rho } \widetilde {{\boldsymbol {U}}} \boldsymbol{\cdot }\left(\boldsymbol{\nabla }\widetilde {Y}\right)+\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\overline { \rho } \widetilde {{\boldsymbol {u}}''Y''}\right)&= 0, \end{align}

(2.9c) \begin{align}\widetilde{\nabla \boldsymbol{\cdot} U}&=0. \end{align}

We now aim to determine a vibro-equilibrium solution for (2.9a –2.9c ), assuming homogeneity along the horizontal directions. This assumption is justified by the premise that the lateral walls of the container are sufficiently far apart relative to the scale of the frozen wave structures under consideration. This solution corresponds to ${\boldsymbol {u}}''=Y''=0$ , an oscillating Favre-averaged velocity aligned with the direction of forcing as $\widetilde {{\boldsymbol {U}}}=V(z) \sin (t) \boldsymbol {e} _x$ and a steady density $\overline {\rho }(z)$ or equivalently mass fraction profile $\widetilde Y(z)$ . We take pressure in the form

(2.10) \begin{equation} P(z,x,t)=P_0-\frac {1}{F}\int _0 ^{z} \overline { \rho }(z') {\textrm d}z'+x \mathcal{P} \cos t, \end{equation}

where $\mathcal P$ is a constant expressing the longitudinal pressure gradient to be defined. We obtain from the momentum equation, (2.9a )

(2.11) \begin{equation} V(z)=1- \frac {\mathcal P}{\overline { \rho }}=1- \frac {\mathcal P}{1-\mathcal A}+\frac {2 \mathcal A \widetilde Y \mathcal P}{1-\mathcal A^2} . \end{equation}

The corresponding vertical shear is then given by

(2.12) \begin{equation} \mathcal S (z,t) =\partial _z V \sin t=2 \mathcal A \frac {\mathcal P}{1-\mathcal A^2}\frac {{\textrm d}\widetilde Y}{{\textrm d}z} \sin t . \end{equation}

The vibro-equilibrium solution is therefore determined by the mean vertical profile of density or mass fraction.

We further consider the following structure for the solution:

1. (i) the upper region, $z\geqslant L/2$ , contains pure light fluid, characterised by $\overline \rho =1-\mathcal A$ or $\widetilde Y=0$ ;

2. (ii) the lower region, $z\leqslant -L/2$ , contains pure heavy fluid, similarly characterised by $\overline \rho =1+\mathcal A$ or $\widetilde Y=1$ ;

3. (iii) the intermediate region $ -L/2 \leqslant z \leqslant L/2$ represents the mixing zone, where the transition between the two fluids occurs. The parameter $L$ defines the width of this mixing zone, corresponding to the frozen wave’s amplitude.

However, to fully characterise the vibro-equilibrium solution, an additional equation must be derived to complement the resolution of (2.11). This supplementary equation is obtained from the conservation of volume, expressed as

(2.13) \begin{equation} \int _{-\infty }^{+\infty } V(z) {\textrm d}z=0. \end{equation}

Notably, the volume-conservation condition, in the case where $L$ is small compared with the height of the container and the mixing zone is centred in the middle of the container, leads to the particular solution

(2.14) \begin{equation} \mathcal P=1-\mathcal A^2 , \end{equation}

with $V=\mathcal A$ in the upper region and $V=-\mathcal A$ in the bottom region. Conversely, if (2.14) is not satisfied, the velocity in the upper and lower regions may not be opposite.

As demonstrated by (2.12), the mean shear responsible for the frozen wave instability is proportional to the gradient of the Favre-averaged mass fraction. Consequently, as the amplitude of the frozen wave increases, both the shear and the mean mass fraction gradient decrease, eventually stabilising the instability and leading to the saturation of the frozen wave’s amplitude, as detailed in Gréa & Briard (Reference Gréa and Briard2019).

In the following section, we analyse the stability of the vibro-equilibrium solution to derive a criterion for determining the saturated amplitude of the frozen wave.

2.4. Stability analysis

The equations governing the Favre fluctuations in the inertial regime can be derived by subtracting the Favre-averaged equations, (2.9a –2.9c ), from the system equations (2.1a –2.1c ). This process yields

(2.15a) \begin{align}&\quad \partial _t {\boldsymbol {u}}''+\widetilde {{\boldsymbol {U}}} \boldsymbol{\cdot }(\boldsymbol{\nabla }{\boldsymbol {u}}'')+{\boldsymbol {u}}''\boldsymbol{\cdot }\big(\boldsymbol{\nabla }\widetilde {{\boldsymbol {U}}}\big)+{\boldsymbol {u}}''\boldsymbol{\cdot }(\boldsymbol{\nabla }{\boldsymbol {u}}'')-\frac {1}{\overline {\rho }}\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\overline { \rho } \widetilde {{\boldsymbol {u}}''{\boldsymbol {u}}''}\big)\nonumber \\&\qquad =\frac {2 \mathcal A}{1-\mathcal A ^2} Y''(\overline {\boldsymbol{\nabla }P}+(\boldsymbol{\nabla }P)') -\frac {1}{\overline {\rho }} (\boldsymbol{\nabla }P)', \end{align}

(2.15b) \begin{align}& \partial _t Y''+\widetilde {{\boldsymbol {U}}} \boldsymbol{\cdot }(\boldsymbol{\nabla }Y'')+{\boldsymbol {u}}''\boldsymbol{\cdot }\left(\boldsymbol{\nabla }\widetilde {Y}\right)+{\boldsymbol {u}}''\boldsymbol{\cdot }(\boldsymbol{\nabla }Y'')-\frac {1}{\overline {\rho }}\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\overline { \rho } \widetilde {{\boldsymbol {u}}''Y''}\big)= 0, \end{align}

(2.15c) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol {u}}''=0. \end{align}

These equations describe the dynamics of the fluctuating components relative to the Favre averages, ${\boldsymbol {u}}''$ and $Y''$ , capturing the contributions to the overall behaviour of the frozen waves in the inertial regime.

By linearising the system (2.15a –2.15c ) around the vibro-equilibrium solution, we obtain

(2.16a) \begin{align} \partial _t {\boldsymbol {u}}''+\widetilde {{\boldsymbol {U}}} \boldsymbol{\cdot }(\boldsymbol{\nabla }{\boldsymbol {u}}'')+{\boldsymbol {u}}''\boldsymbol{\cdot }(\boldsymbol{\nabla }\widetilde {{\boldsymbol {U}}})&= \frac {2 \mathcal A}{1-\mathcal A ^2} Y'' \overline {\boldsymbol{\nabla }P} -\frac {1}{\overline {\rho }} (\boldsymbol{\nabla }P)', \end{align}

(2.16b) \begin{align} \partial _t Y''+\widetilde {{\boldsymbol {U}}} \boldsymbol{\cdot }(\boldsymbol{\nabla }Y'')+{\boldsymbol {u}}''\boldsymbol{\cdot }(\boldsymbol{\nabla }\widetilde {Y})&= 0, \end{align}

(2.16c) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol {u}}''&=0. \end{align}

This system can be further expressed for each components of the velocity ${\boldsymbol {u}}''$ and replacing the mean solution by the vibro-equilibrium solution expressed in § 2.3 as

(2.17a) \begin{align} \partial _t u_x^{\prime\prime}+V \sin (t) (\partial _x u_x^{\prime\prime})+\mathcal S (t,z) u_z^{\prime\prime}&= 2 \mathcal A \frac {\mathcal P}{1-\mathcal A^2}\cos (t) Y'' -\frac {1}{\overline {\rho }} (\partial _x P)', \end{align}

(2.17b) \begin{align} \partial _t u_y^{\prime\prime}+V \sin (t)(\partial _x u_y^{\prime\prime})&= -\frac {1}{\overline {\rho }} (\partial _y P)', \end{align}

(2.17c) \begin{align} \partial _t u_z^{\prime\prime}+V \sin (t) (\partial _x u_z^{\prime\prime})&= -\frac {1}{F}\frac {2 \mathcal A}{1-\mathcal A ^2} Y'' -\frac {1}{\overline {\rho }} (\partial _z P)', \end{align}

(2.17d) \begin{align} \partial _t Y''+V \sin (t)(\partial _x Y'')+u_z^{\prime\prime}\frac {{\textrm d}\widetilde Y}{{\textrm d}z} &= 0, \end{align}

(2.17e) \begin{align} \partial _x u_x^{\prime\prime}+\partial _y u_y^{\prime\prime}+\partial _z u_z^{\prime\prime}&=0. \end{align}

We now analyse the stability of the system, (2.17a –2.17e ), locally at a given position $z$ by considering perturbations with wavelengths much smaller than the characteristic length scale of the mean mass fraction profile $\widetilde Y$ . This approach corresponds to the rapid distortion theory (Hunt & Carruthers Reference Hunt and Carruthers1990; Pope Reference Pope2000), which allows for solutions expressed in terms of a single time-evolving Fourier wavevector $\boldsymbol {k}(t)$

(2.18) \begin{equation} \left \{\!\! \begin{array}{c} {\boldsymbol {u}}''=\hat {\boldsymbol {u}} ''(t,\boldsymbol {k}(0)) e^{i \boldsymbol {k}(t) \,\boldsymbol{\cdot }\,\boldsymbol {x}}, \\ Y''=\hat Y ''(t,\boldsymbol {k}(0)) e^{i \boldsymbol {k}(t) \,\boldsymbol{\cdot }\,\boldsymbol {x}},\\ P'=\hat P '(t,\boldsymbol {k}(0)) e^{i \boldsymbol {k}(t) \,\boldsymbol{\cdot }\,\boldsymbol {x}} .\end{array} \right . \end{equation}

After eliminating the pressure, we obtain a close system for $\hat u''_z$ , $\hat Y''$ and $k_z$

(2.19a) \begin{align}& \partial _t \hat u_z^{\prime\prime}-2 \frac {k_x k_z}{k^2} \mathcal S \hat u_z^{\prime\prime} = \left ( \cos (t) \frac {\mathcal P}{1-\mathcal A^2}\mathbb P_{xz}-\frac {\mathbb P_{zz}}{F(1-\mathcal A ^2)}\right ) 2 \mathcal A \hat Y'', \end{align}

(2.19b) \begin{align}&\qquad\qquad\qquad\qquad\qquad \partial _t \hat Y''+\hat u_z^{\prime\prime}\frac {{\textrm d}\widetilde Y}{{\textrm d}z} = 0, \end{align}

(2.19c) \begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \partial _t k_z= -\mathcal S k_x, \end{align}

where, in (2.19a ), the projector is defined as $\mathbb P_{ij}(\boldsymbol {k})=\delta _{ij}-k_i k_j/k^2$ . This system is very close to the one derived in Gréa & Briard (Reference Gréa and Briard2019), corresponding to (4a–4c) in its dimensional form. The key difference here is that the forcing and the shear are slightly modified to account for non-Boussinesq effects. The stability of a linear system with time-periodic coefficients can be analysed using Floquet theory. Following the approach in Gréa & Briard (Reference Gréa and Briard2019), and neglecting the distortion of the vertical wavenumber $k_z$ (justified if $F \gg 1$ ), this gives

(2.20) \begin{align} &\ddot Y''+2 A B \sin (t) \dot Y''+A(1+B \cos (t)) Y''=0,\nonumber \\&\quad \text{with} \ A=-\frac {{\textrm d}\widetilde Y}{{\textrm d}z} \frac {2 \mathcal A \sin ^2 (\theta )}{F (1-\mathcal A^2)}, \ \text{and} \ B=F \mathcal P \cot (\theta ) \cos (\phi ), \end{align}

where $\theta$ and $\phi$ are the angles corresponding to the spherical coordinates of the wavevector $\boldsymbol {k}$ . It can be demonstrated from Floquet analysis that the stability condition for (2.20) corresponds to

(2.21) \begin{equation} A\leqslant \frac {2}{2+B^2} \ \text{or} \ \left | \frac {{\textrm d}\widetilde Y}{{\textrm d}z}\right | \leqslant \frac {1-\mathcal A^2}{\mathcal A F \mathcal P^2}. \end{equation}

Assuming that the maximum of the Favre-averaged mean mass fraction gradient serves as a good measure of the frozen wave amplitude $L$ , we can derive from criterion (2.21) a simple estimation for the non-dimensional frozen wave amplitude in the inertial regime

(2.22) \begin{equation} L=\mathcal A \frac {\mathcal P^2}{(1-\mathcal A^2)} F. \end{equation}

Equation (2.22) shows that the amplitude of the frozen waves (in dimensional form) scales with the square of the vibrational velocity in the inertial regime. This expression also generalises the formula derived in Gréa & Briard (Reference Gréa and Briard2019) by incorporating a non-Boussinesq correction, thereby extending its validity to larger Atwood numbers. Importantly, the frozen wave amplitude in the inertial regime is independent of the wavelength selected during the growth phase. As a result, for a given Froude and Atwood number, it is possible to obtain frozen waves of identical amplitude but with different wavelengths, as will be shown in simulations.

To derive the amplitude of the frozen waves, we assumed that surface tension and viscous effects are negligible. To evaluate the consistency of this hypothesis, it is useful to consider the non-dimensional capillary length, given by $\ell _\gamma$ , and verify that it remains significantly smaller than the observed amplitude of the frozen waves. Taking $\mathcal P=1-\mathcal A^2$ corresponding to symmetric shear profiles for simplicity, we obtain using (2.22)

(2.23) \begin{equation} \frac {L}{\ell _\gamma }= 4 \mathcal C. \end{equation}

This expression offers a new interpretation of the critical parameter $\mathcal C$ as the ratio between the frozen wave amplitude and the capillary length. This confirms that $\mathcal C$ must be sufficiently large to ensure that the system operates within the inertial regime.

Similarly, viscous effects can be assessed through the non-dimensional Stokes length for the heavy fluid, given by $\delta _\nu =\sqrt {2 \nu _2/\omega }/a=\sqrt {2 /{\textit{Re}}}$ . This expression shows that, at large values of the Froude number $F$ and Reynolds number ${\textit{Re}}$ , the ratio between the frozen wave amplitude and the viscous length becomes large.

As recalled earlier in § 2.2, the most amplified wavelength near the neutral curve is $2 \pi \ell _\gamma$ (Lyubimov & Cherepanov Reference Lyubimov and Cherepanov1986). If this linearly selected wavelength is ‘frozen’ into the final pattern of amplitude $L$ , the frozen wave’s angle $\alpha$ is given by

(2.24) \begin{equation} \tan (\alpha /2)=\frac {\pi \ell _\gamma }{L}=\frac {\pi }{4 \mathcal C}. \end{equation}

Equation (2.24) shows that, as the system departs from neutrality into the inertial regime (i.e. $\mathcal C$ large), the frozen wave angles sharpen – up to the point where secondary instabilities intervene.

In the next section, we present experimental results on frozen waves, designed to evaluate the validity of this theoretical framework.

3.1. Experimental set-up

The frozen wave experiments were performed in a cuboidal container with an inner height $H=40$ cm, width $W=12$ cm and depth $D=3$ cm (see figure 2). The container is constructed from Plexiglas, with its lateral sides reinforced by an aluminium frame. It was carefully designed to withstand the high accelerations generated by the shaker, while preventing resonance modes induced by the forcing. Additionally, the reduced width $W$ of the container was chosen to minimise cavitation at the lateral walls, a phenomenon caused by the large horizontal pressure gradient resulting from the forcing.

Figure 2. (a) Front view of the cuboidal container, highlighting the interface between the silicone oil and the LST heavy fluid. (b) Overview of the experimental set-up, featuring the shaker, two high-speed cameras and the container.

Figure 3. Parameters used for the frozen wave experiments are presented in a Weber–Froude $(\textit{We},F)$ representation, with marker colours corresponding to the Reynolds number ${\textit{Re}}$ values. Note that multiple experiments may be conducted for a single parameter. The dashed line represents the stability criterion proposed by Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986), as given in (2.3). Additionally, visualisations from the high-speed camera showcasing the frozen wave interface in the saturation regime are included in the figure.

The container is securely attached laterally to the shaker, which delivers the periodic forcing (see again figure 2). The forcing can typically reach frequencies as high as $200$ Hz, velocities up to $1$ m s⁻¹ and accelerations as large as $62$ g. The specific parameters for each experimental run are detailed later in figure 3. Notably, the different experimental runs were primarily characterised by the forcing frequency and velocity, ensuring that all conditions remained within the operational limits of the shaker. The sinusoidal forcing was gradually introduced through a slow ramp spread over several hundred vibration periods to maintain system stability. To ensure the specified forcing was accurately delivered, the output from the shaker was carefully monitored using two accelerometers positioned at different locations on the container (see § 3.3).

The light fluid used in the experiments is silicone oil (BLUESIL 604V50), with a density of $\rho _1 = 957 \pm 3$ kg m⁻³ and a kinematic viscosity of $\nu _1 = (55 \pm 2.5) \times 10^{-6}$ m $^2$ s⁻¹. The uncertainties are primarily attributed to room temperature variations within the range of 20–25 °C. The heavy fluid is LST heavy liquid (a solution of lithium heteropolytungstates in water), characterised by a density of $\rho _2 = 2850 \pm 2.5$ kg m⁻³ and a kinematic viscosity of $\nu _2 = (4.1 \pm 0.3) \times 10^{-6}$ m $^2$ s⁻¹. This fluid has been previously employed in studies of turbulent mixing in Rayleigh–Taylor experiments (Roberts & Jacobs Reference Roberts and Jacobs2016).

Consequently, the two-fluid system is defined by an Atwood number of $\mathcal A = 0.497 \pm 0.001$ and a viscosity contrast of $\varUpsilon = 13.4 \pm 0.1$ . The interface between the silicone oil and the LST heavy liquid is positioned at the vertical midpoint of the container. Additionally, the surface tension between the two fluids, measured using the pendant drop method, is found to be $\gamma = 22.5 \pm 2.5$ mN m⁻¹.

To provide a reference point, we evaluate the critical vibrational velocity for our experimental set-up using the inviscid criterion given in (2.3), yielding $(a \omega ) _c=0.34$ m s⁻¹. Associated with this threshold is the most unstable wavelength, which is governed by the capillary length. This critical wavelength is given by $2 \pi \ell _\gamma a=6.92$ mm, indicating that the container width is sufficient to allow the instability to fully develop. These considerations underscore the need for a high-performance shaker in order to explore the frozen wave dynamics well beyond the onset, within the inertial regime.

Two Phantom high-speed cameras were positioned in front of the container (see figure 2) to capture the dynamics of the two fluids and their interface at a frame rate of 400 frames per second and a resolution of $1024 \times 768$ pixels. The cameras were arranged at different viewing angles, and all recordings were conducted in the laboratory frame.

3.2. Phenomenology

During the Vibmix test campaign, 32 experimental runs were performed, covering a wide range of Froude and Weber numbers, as summarised in table 1 and illustrated in figure 3 (see also the supplemental video). The forcing parameters were chosen to ensure that the container’s lateral displacement remained small compared with its width, maintaining a ratio $\mathcal{W} \gt 20$ . Several parameter sets – such as runs f50v0.500a-c, corresponding to a forcing frequency of 50 Hz and a vibrational velocity of $0.5$ m s⁻¹ as shown in figure 4 – were repeated multiple times, confirming the repeatability of the results.

Table 1. Parameters of the experiments conducted during the VIBMIX campaign. The last column reports the measured amplitudes of the mixing zone in the stationary regime.

Figure 3 also shows snapshot insets of the interface in the saturated regime for representative parameter sets. In the red and orange insets, well-formed trochoidal frozen waves appear at moderate Weber and Froude numbers. As these numbers increase (green inset), secondary instabilities gradually distort the waves, leading to atomisation and turbulence (blue and purple insets). Crucially, the ultimate morphology – trochoidal versus atomised – can be predicted by the dimensionless capillary length, $\ell _\gamma =\sqrt {F/{\textit{We}}}$ , with larger $\ell _\gamma$ suppressing secondary breakup (yielding clear trochoidal shapes) and smaller $\ell _\gamma$ allowing instabilities to dominate (producing atomisation).

When intact, the frozen waves exhibit a wavelength close to the most amplified linear mode – around 7 mm – but small variations occur both between runs and within a single experiment, reflecting nonlinear interactions. Crucially, as the critical parameter $\mathcal C$ grows, the frozen wave’s angle sharpens (equivalently, the steepness increases), consistent with our theoretical scaling provided by (2.24).

Figure 4. Visualisation of the experiments runs f50v0.500a-c in the stationary regime corresponding to a frequency $\omega /2 \pi =50$ Hz and a vibrational velocity $0.5$ m s⁻¹. The measures for the frozen waves amplitudes, performed by the methodology described in the manuscript, are also reported.

The effects of large density contrasts are evident not only in the breaking of the up–down symmetry but also in the influence of the lateral walls, which induce an inclination of the interface. These effects become more pronounced under strong forcing.

All the runs presented in figure 3 exhibit a destabilisation of the interface, even though some are predicted to be stable according to the inviscid theory of Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986), (2.3). These discrepancies can likely be attributed to viscosity contrasts destabilising the interface, as reported in numerous studies (Lyubimov et al. Reference Lyubimov, Khilko, Ivantsov and Lyubimova2017).

No drift instability of the frozen waves, as described in García Bennett et al. (Reference García Bennett, Salgado Sánchez, Gligor and Porter2021), was observed in our experiments. This absence can likely be attributed to the geometry of the container, whose width is much smaller than its height ( $W/H=0.3$ ), thereby inhibiting the development of such instabilities.

For runs with large Froude parameters, cavitation phenomena can develop at the lateral walls, making it challenging to further investigate frozen waves in the inertial regime.

3.3. Measures of the frozen wave amplitudes

In this section, we analyse the amplitude of the frozen waves observed in the experiments in greater detail, with the goal of evaluating the accuracy of the predictions provided by (2.22).

Figure 5. (a) Typical renormalised acceleration signals as a function of time $t$ for experiments f150v0.500a (b) and f50v0.500a (c) corresponding to parameters described in table 1. The grey area represents the acceleration measured by the accelerometers fixed to the tank, while the red dashed line shows the acceleration estimated from image displacement. The insert shows the temporal frequencies computed from the acceleration profile. (d) Images of the frozen waves from the same experiment, with red circles indicating the tracked points used by the post-processing algorithm to reconstruct the acceleration profile from the image sequence.

To analyse the images from the cameras, horizontal motion was removed during post-processing using a point-tracking algorithm based on the Lucas–Kanade method, as implemented in OpenCV. This algorithm estimates optical flow for a sparse set of features – typically corners – to determine horizontal displacement. The measured motion aligned well with data from two accelerometers, as shown in figure 5, confirming the consistency and accuracy of the correction procedure.

Figure 6. (a) Time evolution of horizontally averaged images for two experiments with (b) $F=48$ , ${\textit{We}}=11.2$ (f150v0.500a) and (c) $F=62$ , ${\textit{We}}=134.4$ (f100v1.000). The averaging process is performed in the central region of the container to eliminate wall effects. The dashed blue lines, obtained using Otsu’s method, are used to evaluate the frozen wave amplitude, which is also represented by the grey lines in the figures. (d) Visualisation of the interface for the same experiments in the stationary regime, taken at the moment marked by the vertical dashed grey lines in the left-side images.

A spatio-temporal representation of the experiments is obtained by horizontally integrating the images within the central region of $6$ cm width of the container to minimise wall effects. The upper and lower boundaries of the mixing zone are then identified using Otsu’s thresholding method. To determine the amplitude $L$ at saturation, measurements are taken and averaged during the stationary regime, as illustrated in figure 6.

Notably, a comparison between the tank’s acceleration profile shown in figure 5 and the time evolution of the interface in figure 6 reveals that the development of the frozen wave instability is well synchronised with the applied acceleration.

Figure 7. Rescaled square root of the frozen wave amplitudes as a function of the square of the renormalised vibrational velocity $\mathcal C=(a \omega ) ^2/ (a \omega )^2_c$ . The error bars for the experimental data (grey circles) mainly reflect uncertainties in the surface tension and the pixel size of the images. The hexagon markers represent the two-dimensional and three-dimensional simulation results described in § 4.

The measurements of the square root of the frozen waves amplitudes are presented in figure 7 as a function of the renormalised square of the vibrational velocity, defined by (2.3). The results are compensated, as outlined by (2.22), to recover the constant $\mathcal P$ in the inertial regime. Remarkably, for $\mathcal C \geqslant 2$ , the data converge well toward a constant value, providing strong support for the stability analysis developed in § 2.4. Furthermore, the constant corresponds to $\mathcal P=1-\mathcal A^2$ , which aligns with the symmetric shear profile solution derived from the volume-conservation condition.

Since the experimental measurements do not provide access to the velocity field, we present a numerical study in the next section to further validate the prediction given by (2.22).

4.1. Numerical set-up

The direct numerical simulations were conducted using a finite-volume solver implemented in the open-source code Basilisk (Popinet Reference Popinet2009, Reference Popinet2015). The solver is well suited for fluid dynamics problems involving interfacial instabilities (Fuster & Rossi Reference Fuster and Rossi2021). It employs a coupled volume-of-fluid/embedded boundary method to track and reconstruct interfaces between immiscible fluids and solid walls (Tavares et al. Reference Tavares, Josserand, Limare, Lopez-Herrera and Popinet2024). Viscous terms are treated implicitly, while a projection method is used to handle velocity–pressure coupling. In this study, the governing equations are discretised using a tree-based adaptive Cartesian grid with a second-order scheme in both time and space. The minimum grid spacing is approximately ten times smaller than the capillary length. For simulations representative of the experiments, we typically use 1024 grid points along the width of the domain, resulting in a spatial resolution of $0.12$ mm. In comparison, the capillary length is estimated as $\ell _\gamma a=(\gamma / (\rho _2-\rho _1) g)^{1/2}=1.1$ mm. We conducted a convergence study to demonstrate that the main features of the frozen waves are accurately captured at this resolution.

The initial conditions assume the fluid is at rest, and the external acceleration is applied through a ramp function, as in the experiments. Small random initial perturbations are introduced at the interface, characterised by an annular power spectrum, following the approach described in Thévenin et al. (Reference Thévenin, Gréa, Kluth and Nadiga2025). Boundary conditions include a no-slip condition, and a contact angle of $90^o$ is imposed on all boundaries. The resulting system is solved using an iterative cycle with a small relative tolerance to limit numerical instabilities.

The numerical investigation of frozen waves comprises three series of simulations, each designed with a specific objective:

1. (i) series A: these simulations were designed to closely replicate the experimental configuration, enabling a direct comparison between numerical and experimental results;

2. (ii) series B: this series investigated the effect of lateral confinement on the frozen wave amplitudes by varying the width $W$ (or equivalently $\mathcal W$ ) of the domain;

3. (iii) series C: in these simulations, the wall boundary conditions were replaced with periodic boundary conditions in the $(x,y)$ plane. Additionally, the forcing acceleration was replaced by a time-periodic pressure gradient applied along the $x$ -direction, analogous to the approach of Bunner & Tryggvason (Reference Bunner and Tryggvason2002) in bubbly flows.

The simulation series includes both three-dimensional (3-D) simulations, designed to closely replicate the phenomenology observed in experiments, and more computationally efficient 2-D simulations, which enable a broader exploration of the parameter space. Notably, in the 2-D configurations, the resolution can be increased up to $8192$ grid points across the width of the domain, ensuring a minimum of 10 grid points per capillary length for accurate resolution of the interfacial dynamics.

4.2. Results

Figure 8. Visualisation of the interfaces of the frozen wave instability in the stationary regime from 3-D Basilisk simulations, configured to replicate the experimental set-up (series A). The simulations (a–e) are sampled for different values of the critical parameter $\mathcal C$ , achieved by varying $a$ while keeping the forcing frequency fixed at $\omega /2 \pi =200$ Hz. The inset shows a zoomed view of the interface highlighting the development of secondary Faraday instabilities.

Figure 9. (a) Spatio-temporal diagram of the mean interface height $\eta (x,t)$ and (b) temporal evolution of $L(t)$ obtained from 2-D Direct Numerical Simulation (DNS) using the nominal parameters of the experiments with walls and a forcing condition of $0.5$ m s⁻¹ and $200$ Hz (run f200v0.500). The ramp profiles of accelerations, converging towards the terminal value of $62$ g, are also indicated in the figures corresponding to the red curves. The initial perturbations have wavenumbers $k_0$ confined within the bands: from top to bottom, $[10,20]$ , $[7,14]$ and $[6,12]$ , respectively. The insets show a snapshot in the steady regime.

A representative example of the 3-D simulations from series A, based on the experimental configuration, is shown in figure 8. These simulations use parameters close to those of the experiments, with the light fluid characterised by $\rho _1 = 950$ kg m⁻³ and $\nu _1 = 52 \times 10^{-6}$ m $^2$ s⁻¹, and the heavy fluid by $\rho _2 = 2850$ kg m⁻³ and $\nu _2 = 4.5 \times 10^{-6}$ m $^2$ s⁻¹. Additionally, the surface tension is set to 20 mN m⁻¹. The simulations were carried out for various values of the control parameter, $\mathcal{C} \in \{0.58,\ 1.3,\ 2.3,\ 3.6,\ 5.19\}$ , allowing the exploration of different regimes of the frozen wave instability, from its onset to the inertial-dominated regime.

The simulations successfully reproduce the amplitude of the frozen waves observed in the experiments and validate at large $\mathcal C$ the predictions given by (2.22) (see figure 7). They also accurately capture the experimental wavelength and the vibro-equilibria effects near the walls. Furthermore, the simulations shed light on secondary Faraday instabilities (see also supplemental movie 2), as reported by Shevtsova et al. (Reference Shevtsova, Gaponenko, Yasnou, Mialdun and Nepomnyashchy2016), Lyubimova et al. (Reference Lyubimova, Ivantsov, Garrabos, Lecoutre and Beysens2019), which converge toward a homogeneous layer in the inertial regime.

While the initial perturbation of the interface is not known in the experiments, it can be fully controlled in the simulations. In figure 9, we present results from 2-D simulations designed to reproduce run f200v0.500 (see table 1). The initial conditions correspond to a periodic perturbation of the interface $\eta$ , with Fourier modes $k_0$ (where $k_0=1$ corresponds to the width of the tank) confined within specified bands. The final stationary regime exhibits a relatively weak dependence on the initial conditions, with the number of frozen waves settling between 8 and 10, in good agreement with the experimental observations. The formation of these frozen waves – particularly their wavelength – can be monitored using spatio-temporal diagrams of the interface. The instability requires some time to develop, as the acceleration ramp gradually reaches its final value. The dominant wavelength that emerges is likely governed by capillary effects as expected by the inviscid theory of Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986), characterised by a critical wavelength of $2 \pi \ell _\gamma a=6.9$ mm, corresponding to a wavenumber $ \simeq 17$ . Interestingly, the initial wavelength continues to grow until the system reaches the final stationary regime, while the orientation of the interface and the angle of the frozen waves remain approximately constant.

Figure 10. (a) Favre-averaged mean mass fraction profiles, $\widetilde {Y}$ , and (b) amplitude of the oscillating Favre-averaged velocity, $V$ , for the different 3-D simulations of series A. These simulations, corresponding to varying critical parameters, are the same as those presented in figure 8. The averaging process is conducted in the central region of the container to eliminate the influence of the wall effects.

The simulations provide access to the mean mass fraction profiles, $\widetilde {Y}$ , and the oscillating Favre-averaged velocity, $V$ , as shown in figure 10. The quasi-linear mean mass fraction profiles suggest that their local gradient serves as a good measure of the frozen wave amplitudes, i.e. $\left | ({{\textrm d}\widetilde Y}/{{\textrm d}z})\right |=1/L$ . Additionally, the $V$ profiles reach values close to $\pm \mathcal A$ just outside the mixing zone in the pure fluids and exhibit a symmetric shear structure. Although the velocity profiles decay away from the interface – reflecting the influence of the top and bottom walls – this decay does not break the overall symmetry. Consequently, the volume-conservation condition can still be applied unchanged to derive the theoretical predictions, recovering the constant $\mathcal P=1-\mathcal A^2$ .

Figure 11. Visualisation of the frozen waves in the stationary regime for 2-D simulations (series B) at (a) $\mathcal C=2.3$ , (b) $5.19$ , (c) $9.22$ and for different domain widths $\mathcal W$ . (d) Rescaled square root of the frozen wave amplitudes as a function of the renormalised width $\mathcal W$ for the series B simulations and Vibmix experiments.

However, the influence of the sidewalls can still be anticipated. In our simulations, this can be systematically investigated by varying the channel width $W$ or equivalently the non-dimensional number $\mathcal W=W/a$ in the simulations. The results of 2-D simulations from series B between $\mathcal{C} = 2.3$ and $\mathcal{C} = 9.22$ at different width parameters $\mathcal{W}$ are presented in figure 11. The effect of lateral confinement, which diminishes the amplitudes of frozen waves, is clearly visible in the simulation snapshots.

To better quantify this effect, figure 11 also shows the rescaled square root of the frozen wave amplitudes, $L$ , as a function of $\mathcal{W}$ . The simulation results collapse onto a single curve, suggesting that confinement effects are not strongly dependent on the parameter $\mathcal{C}$ when $\mathcal{C}$ is large. Additionally, the results indicate that lateral confinement has only marginal effects when $\mathcal{W}$ exceeds 200. This criterion proves valuable for the design of future experiments, as increasing the length of the container is desirable but also increases the likelihood of cavitation effects. Furthermore, it highlights that confinement may influence the frozen wave amplitudes observed in the Vibmix experiments, underscoring the importance of accounting for geometric constraints when interpreting experimental results.

Figure 12. Examples of 2-D homogeneous simulations from series C in the stationary regime, illustrating (a) unbroken, (b) partially fragmented and (c) fully fragmented frozen waves. A resolution of at least ten grid points per capillary length ensures faithful capture of interface dynamics. Panels (a–c, left) Instantaneous frozen wave patterns (heavy fluid in black). Panels (a–c, centre) Favre-averaged mass fraction profiles $\widetilde Y$ at several time instants within one oscillation period. Panels (a–c, right) Corresponding Favre-averaged velocity profiles $\widetilde U$ .

Figure 13. (a) Parameters corresponding to the series C homogeneous simulations are represented in the Weber–Froude plane, with markers coloured according to the Reynolds number. The dashed line represents the neutral stability curve proposed by Lyubimov & Cherepanov (Reference Lyubimov and Cherepanov1986), as described in (2.3). (b) Rescaled square root measures of the frozen wave amplitudes from the same series C simulations are plotted as a function of the renormalised square of the vibrational velocity, $\mathcal{C}$ . The triangle markers indicate the parameters and measures from the simulation in figures 12 and 14.

Figure 14. (a) Spatio-temporal diagram of the mean interface height $\eta (x,t)$ and (b) temporal evolution of the renormalised $L(t)$ obtained from 2-D homogeneous DNS with parameters indicated in the figure. The insets show a snapshot in the steady regime.

To eliminate the effects of lateral confinement, we present in figures 12 and 13 the results from the homogeneous simulations of series C, which cover a wide range of parameters in the ${\textit{Re}},F,We$ space, including the parameters used in the experimental set-up. These simulations employ horizontal periodic boundary conditions and an imposed oscillating pressure gradient in the $x$ -direction to drive the flow. For instance, the first example in figure 12 shares the same parameters as the case shown in figure 8 with $\mathcal C=2.3$ . By eliminating the influence of the walls, the periodic simulation recovers the amplitude of a simulation with $\mathcal{W}\gg 1$ and the sawtooth pattern is no longer tilted toward the sides.

The results confirm that, when the critical parameter $\mathcal C$ and the Reynolds number ${\textit{Re}}$ are sufficiently large, the predicted frozen wave amplitudes, as given by (2.22), are accurately recovered. This validates the theoretical scaling with the square of the vibrational velocity in the inertial regime.

Notably, the simulations also highlight the influence of the Reynolds number: at low ${\textit{Re}}$ , viscous effects can counteract inertia and inhibit the formation of frozen waves, whereas at high ${\textit{Re}}$ , the flow tends to become turbulent and homogenised. Configurations with high Reynolds numbers, which align more closely with the assumptions underlying the theory, are therefore better captured by the model. In contrast, at lower ${\textit{Re}}$ , the frozen wave amplitudes are slightly underestimated.

At the same time, the simulations confirm that the dominant wavelength is selected during the early stages of the instability and is strongly correlated with the capillary length. In figure 14, we present simulations performed at fixed Atwood, Froude and Reynolds numbers, while varying the capillary length by adjusting the surface tension. This variation leads to different Weber and critical wavenumbers, which in turn influences the selected wavelength. As the amplitude of the frozen wave increases, the wavelength also evolves slightly until the system reaches the saturation regime. As a result, configurations with identical frozen wave amplitudes but different wavelengths – and therefore different angles – can be observed. It is worth noting that this observation refines the original interpretation of frozen waves proposed by Wolf (Reference Wolf1970), who drew an analogy with the Kapitza pendulum. Our results indicate that it is not the angle that is universal, but rather the amplitude of the frozen waves.