The onset and saturation of the Faraday instability in miscible fluids in a rotating environment

Abstract

We investigate the influence of rotation on the onset and saturation of the Faraday instability in a vertically oscillating two-layer miscible fluid using a theoretical model and direct numerical simulations (DNS). Our analytical approach utilizes Floquet analysis to solve a set of the Mathieu equations obtained from the linear stability analysis. The solution of the Mathieu equations comprises stable and harmonic, and subharmonic unstable regions in a three-dimensional stability diagram. We find that the Coriolis force delays the onset of the subharmonic instability responsible for the growth of the mixing zone size at lower forcing amplitudes. However, at higher forcing amplitudes, the flow is energetic enough to mitigate the instability delaying effect of rotation, and the evolution of the mixing zone size is similar in both rotating and non-rotating environments. These results are corroborated by DNS at different Coriolis frequencies and forcing amplitudes. We also observe that for $(\,f/\omega )^2<0.25$, where $f$ is the Coriolis frequency, and $\omega$ is the forcing frequency, the instability and the turbulent mixing zone size-$L$ saturates. When $(\,f/\omega )^2\geq 0.25$, the turbulent mixing zone size-$L$ never saturates and continues to grow.

1. Introduction

The interface of two immiscible fluids enclosed in a vertically vibrating container upon reaching a certain frequency and acceleration becomes unstable. This unstable phase is known as Faraday instability and results in the formation of nonlinear standing waves at the interface of the two fluids. These waves are termed as Faraday waves and were first described by Faraday (1831). The onset of Faraday instability in immiscible liquids was theoretically demonstrated by Benjamin & Ursell (1954) using linear stability analysis of the interface of an ideal fluid governed by the set of equations relevant to the system of Mathieu equations. Linear stability of the immiscible, viscous and finite depth fluid problem for one and two frequency excitation was studied by Kumar & Tuckerman (1994) and Besson, Edwards & Tuckerman (1996), respectively, using Floquet analysis. Apart from theoretical investigations, experimental (Gollub & Meyer 1983; Douady & Fauve 1988; Simonelli & Gollub 1989; Douady 1990; Müller 1993; Edwards & Fauve 1994; Binks & van de Water 1997; Kudrolli, Pier & Gollub 1998; Arbell & Fineberg 2002; Westra, Binks & Van De Water 2003; Kityk et al. 2005; Rajchenbach, Clamond & Leroux 2013; Shao et al. 2020, 2021) and numerical (Perinet, Juric & Tuckerman 2009; Takagi & Matsumoto 2011; Kahouadji et al. 2015; Takagi & Matsumoto 2015) studies were also performed by many researchers to understand the Faraday wave patterns at the interface of immiscible liquids.

The majority of the studies on Faraday waves deals with immiscible liquids. However, in the past few decades researchers have shifted their focus to understanding the dynamics of Faraday waves in miscible liquids. Zoueshtiagh, Amiroudine & Narayanan (2009) performed experiments and two-dimensional numerical simulations to investigate the Faraday instability of diffuse interfaces between pairs of miscible liquids of different densities in a rectangular cell. They observed that, above a certain forcing amplitude, the standing waves that appear on the diffuse interface are highly disorganized and interact with each other leading to the mixing of fluids followed by the disappearance of the waviness on the interface when the two fluids are mixed. They also found the mixing phase to be subharmonic, similar to the subharmonic nature of the Faraday waves in immiscible fluids. Further, this finding was supported by Diwakar et al. (2015) who used Floquet theory in conjunction with a quasisteady approximation to carry out linear analysis of Faraday instability in miscible fluids. An experimental and numerical framework similar to Zoueshtiagh et al. (2009) has been used by Amiroudine, Zoueshtiagh & Narayanan (2012) to report the exponential growth of the mixing layer thickness owing to fingering at the interface, followed by a small growth rate that demonstrates the saturation of mixing. Recently, Gréa & Adou (2018) conducted three-dimensional (3-D) numerical simulations and demonstrated that Faraday instability can generate turbulent mixing zones for strong forcing parameter and/or for sufficiently random initial condition at the interface of two miscible fluids. Further, they formulated a system of equations based on the second-order correlation spectra for turbulent quantities and used perturbation analysis to estimate the final size of the mixing zone,

(1.1)\begin{equation} L_{sat}=\frac{2\mathcal{A}g_0}{\omega^2}(2F+4), \end{equation}

where $L_{sat}$ is the final size of mixing zone (defined later in (2.5)) in the saturation state, $\mathcal {A}=(\rho _1-\rho _2)/(\rho _1+\rho _2)$ is the Atwood number expressing the density contrast between heavy fluid ($\rho _1$) and light fluid ($\rho _2$), $g_0$ is the mean acceleration, $F$ is the oscillation amplitude and $\omega$ is the forcing frequency. Gréa & Adou (2018) validated the predicted $L_{sat}$ against their simulations for a wide range of parameters within the homogeneous framework and full-inhomogeneous systems of two miscible fluids. They reported that the instability of the diffuse interface begins with a small harmonic phase, but the main instability phase is dominantly subharmonic, where the turbulent mixing zone develops. The irreversible mixing associated with the growth of the mixing zone owing to the harmonic to subharmonic transition was numerically quantified by Briard, Gréa & Gostiaux (2019) using potential energies. The total potential energy (TPE) was split into background potential energy (BPE) and available potential energy (APE). They demonstrated that the BPE, which signifies the measure of irreversible mixing, increases after the transition. The APE, which denotes the fraction of the TPE that can be converted to BPE through irreversible mixing, peaks at saturation and is partially released in the flow as BPE. This increase in the BPE causes the numerically obtained final mixing zone size $L$ to exceed the theoretically predicted $L_{sat}$ (1.1). To further elucidate the dynamics of the turbulent mixing zone driven by the Faraday instability, Briard, Gostiaux & Gréa (2020) performed experiments with fresh and salty water and supported their finding with direct numerical simulations (DNS) and theoretical predictions. They concluded that when the instability is triggered, a natural wavelength appears at the interface between the two fluids. As the amplitude of this wave increases, well-defined structures break up to produce a turbulent mixing layer. At the saturation of the instability, turbulence is inhibited and a mixing layer of final size consistent with the analytically predicted $L_{sat}$ was obtained for a wide range of parameters. Mondal & Kumar (2004), performed a linear stability analysis to investigate the effects of Coriolis force on the Faraday waves in a thin sheet of viscous fluid placed on the vibrating plate. They used Floquet theory to solve the set of equations in the form of Mathieu equations and explained the effect of rotation rate on the unstable and stable regions of the stability diagrams. They found that the subharmonic waves get suppressed with increasing rotation rates resulting in the delay of the onset of the Faraday waves. They also reported the existence of a tricritical point at the onset of the instabilities, where subharmonic, harmonic and superharmonic waves are simultaneously excited.

Many researchers explored the physics associated with the onset of Faraday instabilities and the subsequent turbulent mixing layer in two miscible fluids. However, the effect of the Coriolis force on the onset and saturation of the Faraday instabilities that drive turbulent mixing in miscible fluids is absent from the literature. In stably stratified fluids, the turbulent mixing mechanisms depend on the nature of external forcing and background stratification. The vertical time-varying displacement of the density surfaces in the ocean is mainly associated with tidal oscillations (Garrett & Kunze 2007; Sarkar & Scotti 2017). These oscillations generate internal waves that can interact with ocean currents and topography, resulting in instability and turbulent mixing in the deep ocean. In oceanic flows, the Earth's rotation significantly influences the dynamics of mixing processes in the thermocline (Noh & Long 1990; Fernando 1991; Fleury et al. 1991). For example, $N/f \approx 4.67$, where $N$ is the Brunt Väisälä frequency, in the abyssal Southern Ocean at midlatitude (Nikurashin, Vallis & Adcroft 2013; Rosenberg et al. 2015). The presence of rotation in stratified flows leads to a generation of inertial-gravity waves (Hanazaki 2002; Smith & Waleffe 2002; Liechtenstein, Godeferd & Cambon 2005; Ferrari & Wunsch 2009; Sarkar & Scotti 2017). When stratification is strong and rotation weak $N/f>1$, the motions are primarily controlled by gravity wave dynamics with maximum frequency $\approx N$, resulting in the propagation of internal gravity waves. However, when rotation is stronger than stratification $N/f<1$, inertial waves with maximum frequency $\approx f$ are generated. Both rotation and stratification produce highly anisotropic flow patterns with quasihorizontal velocity fields (Liechtenstein et al. 2005; Praud, Sommeria & Fincham 2006). However, vertically elongated column-like structures are formed in the rotation-dominant ($N/f<1$) flows, whereas horizontally layered pancake-like structures are formed in stratification-dominant flows ($N/f>1$). Praud et al. (2006) performed experiments to investigate the decaying grid turbulence in a rotating stratified fluid. They found an important effect of rotation that the increase in rotation rate ($\,f$) inhibits the turbulent kinetic energy decay. Most studies in rotating stratified flows have focused on investigating the impact of the ratio $N/f$ on the direct and inverse energy cascades depending on the forcing scales at which the energy is injected into the flow (Smith & Waleffe 2002; Deusebio, Vallgren & Lindborg 2013; Marino et al. 2013, 2014; Herbert et al. 2016; Alexakis & Biferale 2018; Pouquet et al. 2018; Li et al. 2020). Marino et al. (2013) conducted DNS and showed that the interplay between rotation and stratification plays an important role in the inverse energy cascade, where the small scales transfer energy to large scales. They found that the inverse energy cascade is most efficient in the range $1/2\leq N/f \leq 2$. These studies have demonstrated how rotation is crucial in stratified turbulent flows. Therefore, it would be interesting to explore the influence of rotation (or ratio $N/f$) on stratified turbulent flows when the external forcing is time-periodic, which is fundamentally different from that used in the literature mentioned above. The generation of turbulent kinetic and potential energy depends on the strength of the periodic forcing. This gives us an opportunity to investigate how the ratio $N/f$ affects the onset, growth and saturation/non-saturation of the turbulent mixing driven by the Faraday instability in the stably stratified miscible fluids. The Faraday instability, which arises when the interface of two fluids is subjected to time-periodic accelerations has the potential to promote mixing in stably stratified fluids. Therefore, understanding the dynamics of the Faraday instability under the influence of rotation can help to gain some insights into turbulent mixing processes in the deep ocean owing to the time-varying oscillations of density surfaces caused by tidal forcing. The present investigation of rotating miscible fluids with time-periodic accelerations will provide a framework to understand the physics associated with the mixing processes in oceanic flows. We present a linear stability analysis of the Faraday waves in miscible fluids under the influence of the Coriolis force. We aim to investigate the effects of the Coriolis force on the harmonic to subharmonic transition, the onset of the subharmonic instability phase, and the saturation phase of the instability. Our analytical model accounts for the vertical inhomogeneity in the background density profile and the Coriolis effect. We followed the classical approach of Yih (1960) to derive the linear equations of the form of the Mathieu equations and use Floquet theory to solve them. We draw the stability diagrams for the corresponding stable and unstable solutions of the Mathieu equations. We also perform DNS at different forcing amplitudes and Coriolis frequencies to check the validity of our analytical model.

The problem formulation, which includes the governing equations and theoretical framework for linear stability analysis, is provided in § 2. The details of the numerical methodology, and simulation parameters, are given in § 3. The predictions of the linear stability analysis, and the numerical results, are explained in § 4. We conclude § 5.

2. Problem formulation

We present the details of the configuration and the governing equations in § 2.1. Then we perform the linear stability analysis and derive the Mathieu equations governing the stability of a full inhomogeneous inviscid rotating system under vertically periodic forcing in § 2.2. We determine the characteristic frequencies $\varOmega _i$ (defined later) of the Mathieu equations of the full inhomogeneous problem for the $N>f$ case in § 2.3. Based on the values of $\varOmega _i$ for the $N>f$ case, we investigate the linear theory for the limiting case of the fully developed mixing zone in § 2.3.1 and the limiting case of the interface in § 2.3.2. Finally, we compute the $\varOmega _i$ for the $N< f$ and $N=f$ cases in §§ 2.4 and 2.5, respectively.

2.1. Governing equations

We consider two miscible fluids of small contrasting densities. The lighter fluid of density $\rho _2$ is above the denser fluid of density $\rho _1$ in a cubical domain (see figure 1). This domain is subject to periodic vertical oscillations with acceleration $g(t)=g_0(1+F\cos {(\omega t))}$, where $g_0$ is the mean acceleration, $F$ is the oscillation amplitude and $\omega$ is the forcing frequency. We are assuming that the density of the mixture is linearly varying with mass concentration, with $C(\rho _1)=1$ and $C(\rho _2)=0$. The 3-D conservation equations for mass, momentum and concentration field under vertically periodic forcing $g(t)$, for unsteady incompressible flow with Boussinesq approximation, are

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

(2.1b)$$\begin{gather}\frac{\partial \boldsymbol{U}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} =-\boldsymbol{\nabla} P - f\hat{k} \times \boldsymbol{U} + \nu {\nabla}^2 \boldsymbol{U} - 2\mathcal{A}{C}g(t) \hat{k}, \end{gather}$$

(2.1c)$$\begin{gather}\frac{\partial{C}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} {C} = \kappa {\nabla}^2 {C}. \end{gather}$$

Here $\boldsymbol {U}$ represent the velocity vector with components ($U_1$, $U_2$ and $U_3$) in the $x_1$, $x_2$ and $x_3$ (vertical) directions, respectively, ${P}$ is the pressure deviation from the hydrostatic background state, $f$ is the Coriolis's frequency, $\nu$ ($\mathrm {m}^2\ \mathrm {s}^{-1}$) is the kinematic viscosity and $\kappa$ ($\mathrm {m}^2 \mathrm {s}^{-1}$) is the diffusion coefficient.

Figure 1. Schematic of the domain containing two miscible liquids subjected to vertically periodic forcing.

We decompose a variable into its mean $\langle B \rangle _H$ and fluctuating component $b$ following Reynolds decomposition,

(2.2)\begin{equation} B=\langle B \rangle _H + b, \end{equation}

where $\langle B \rangle _H$ is the average in the $x_1$ and $x_2$ horizontal directions, and $\langle b \rangle _H=0$. We consider periodic boundary conditions in both $x_1$ and $x_2$ directions and assume the quantities are statistically invariant in these directions. Therefore, we can define the horizontal average for any quantity, say $B(\boldsymbol {x},t)$, as follows:

(2.3)\begin{equation} \langle B \rangle_H({x_3},t)=\lim_{\substack{(l_{x_1},l_{x_2}) \to (\infty,\infty) }} \frac{1}{l_{x_1} l_{x_2}} \int_{{-l_{x_1}}/{2}}^{{+l_{x_1}}/{2}} \int_{{-l_{x_2}}/{2}}^{{+l_{x_2}}/{2}} B (\boldsymbol{x},t)\,\mathrm{d}\kern0.7pt x_1\,\mathrm{d}\kern0.7pt x_2. \end{equation}

Substituting Reynolds decomposition (2.2) into the governing equations (2.1a), (2.1b) and (2.1c), with mean velocity field $\langle \boldsymbol {U} \rangle _H =0$ and $\langle \boldsymbol {u} \rangle _H=\langle c \rangle _H=\langle \,p \rangle _H=0$, we obtain the equations for fluctuating velocity $\boldsymbol {u}(\boldsymbol {x},t)$ and concentration $c(\boldsymbol {x},t)$ fields as follows:

(2.4a)$$\begin{gather} \frac{\partial u_i}{\partial x_i}=0, \end{gather}$$

(2.4b)$$\begin{gather}\frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} =- \frac{\partial p}{\partial x_i} + f\epsilon_{ij3} u_j\hat{e}_3 - 2\mathcal{A}g(t)c\delta_{i3} + \frac{\partial \langle u_i u_3 \rangle_H}{\partial x_3}+ \nu {\nabla}^2 u_i, \end{gather}$$

(2.4c)$$\begin{gather}\frac{\partial c}{\partial t} + u_j \frac{\partial c}{\partial x_j} =- u_3 \frac{\partial \langle C \rangle_H}{\partial x_3} + \frac{\partial \langle u_3 c \rangle_H}{\partial x_3} + \kappa {\nabla}^2 c. \end{gather}$$

We will compute the mixing zone width $L(t)$ from the mean concentration profile (Andrews & Spalding 1990) as

(2.5)\begin{equation} L(t)=6\int_{-\infty}^{+\infty} \langle C \rangle _H (x_3,t) (1-\langle C \rangle _H (x_3,t))\,\mathrm{d}\kern0.7pt x_3. \end{equation}

We define the Brunt Väisälä (or stratification) frequency as follows:

(2.6)\begin{equation} N=\sqrt{-2\mathcal{A} g_0\frac{\partial \langle C \rangle_H}{\partial x_3}} = \sqrt{-2\mathcal{A} g_0 \varGamma}. \end{equation}

Here, $\varGamma = {\partial \langle C \rangle _H}/{\partial x_3}$ is the vertical gradient of mean concentration, which can be approximated as ${\partial \langle C \rangle _H}/{\partial x_3}=-1/L$.

2.2. Linear theory for the inhomogeneous inviscid rotating system under parametric oscillations

Using similar formulation, assumptions and analytical techniques as in Briard et al. (2020), we extend their linear stability analysis for the inviscid inhomogeneous system to rotating cases and develop a theoretical model to investigate the effects of initial stratification along with the Coriolis force on the triggering and development of instability. We define the initial stratification by Atwood number ($\mathcal {A}$) and initial mixing zone size ($L_0$). We consider a piecewise concentration profile, constant in the bottom and upper pure (unmixed) fluids and linear in the mixing layer, with a vertical gradient of concentration $\varGamma =0$, in the unmixed region and $\varGamma =\text {constant}$, in the mixed region. We also assume that the fluids are inviscid and the fluctuations in velocity and concentration are small. After neglecting the nonlinear, viscous and diffusion terms from (2.4b) and (2.4c) and then eliminating the pressure term, $u_1$ and $u_2$ from the resulting equation of fluctuating velocity $\boldsymbol {u}(\boldsymbol {x},t)$, we obtain (see details in Appendix A)

(2.7a)$$\begin{gather} \frac{\partial^2}{\partial t^2}(\nabla^2 u_3) =-f^2\frac{\partial^2 u_3}{\partial x_3^2} + 2 \mathcal{A} g(t) \varGamma \nabla^2_{H} u_3, \end{gather}$$

(2.7b)$$\begin{gather}\frac{\partial c}{\partial t} =- u_3\varGamma . \end{gather}$$

Here, $\nabla ^2_{H} \equiv \partial ^2/\partial x_1^2 + \partial ^2/\partial x_2^2$ is the horizontal Laplacian operator. To solve (2.7a) we first consider a case constant forcing, $g(t)=g_0$ and assume a solution for $u_3$ of the form

(2.8)\begin{equation} u_3(x_1,x_2,x_3,t) = \phi(x_3) \, {\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)}, \end{equation}

where $\phi (x_3)$ is the amplitude, $\varOmega$ is the temporal frequency or characteristic frequency, and $k$ and $l$ are the components of the wavevector in the horizontal directions $x_1$ and $x_2$, respectively. Substitution of (2.8) in (2.7a) yields (see details in Appendix B)

(2.9)\begin{equation} \frac{\partial^2\phi}{\partial x_3^2} + {K}^2 \frac{( N^2 -\varOmega^2)}{(\varOmega^2 - f^2 )} \phi =0. \end{equation}

Here, $K$ is the horizontal wavenumber which is defined as $K=\sqrt {k^2 + l^2}$. Now we focus on solving the case with vertically periodic forcing (see (2.7a)). The horizontal directions $x_1$ and $x_2$ are assumed to be periodic. Therefore, we apply a Fourier transform on $u_3(x_1,x_2,x_3,t)$ with respect to $x_1$ and $x_2$ such that

(2.10)\begin{equation} u_3(x_1,x_2,x_3,t) = \hat{u}_3(k,l,x_3,t)\, {\rm e}^{{\rm i} (kx_1+lx_2)}. \end{equation}

Substituting the solution (2.10) in (2.7a), we get

(2.11)\begin{equation} \frac{\partial^2}{\partial t^2} \left( \frac{\partial^2 \hat{u}_3 }{\partial x_3^2} - K^2 \hat{u}_3 \right) =-f^2 \frac{\partial^2 \hat{u}_3}{\partial x_3^2} - 2\mathcal{A} g(t) \varGamma K^2 \hat{u}_3. \end{equation}

We further decompose the amplitude $\hat {u}_3(k,l,x_3,t)$ on the basis of previously calculated solutions $\phi _i$ of the stationary system (2.9), as follows:

(2.12)\begin{equation} \hat{u}_3(k,l,x_3,t) =\sum_{i} A_i(k,l,t)\phi_i(k,l,x_3). \end{equation}

Here $A_i(k,l,t)$ represents the time-dependent part of the amplitude of vertical velocity field modes evolving as a result of periodic forcing and $\phi _i(k,l,x_3)$ corresponds to the spatially ($x_3$) dependent amplitude of the vertical velocity field modes. Substituting the decomposition (2.12) in (2.11), we get

(2.13)\begin{equation} \sum_{i} \frac{\partial^2}{\partial t^2} \left( A_i\frac{\partial^2 \phi_i }{\partial x_3^2} - K^2 \phi_i A_i \right) =- \sum_{i}f^2 A_i \frac{\partial^2 \phi_i }{\partial {x_3}^2} - \sum_{i}2\mathcal{A} g(t) \varGamma K^2 A_i \phi_i. \end{equation}

Since we are interested in the temporally evolving amplitude ($A_i$) with periodic forcing, we eliminate $\phi _i$ and $K$ from (2.13) by substituting (2.9) to obtain (see Appendix C)

(2.14)\begin{equation} \sum_{i} \frac{\partial^2 A_i}{\partial t^2} = \sum_{i}\left( -\frac{f^2 ( N^2 -\varOmega_i^2 ) + N^2( \varOmega_i^2 - f^2 ) (1+F\cos{\omega t}) }{N^2 - f^2} \right) A_i. \end{equation}

Now, we define the amplitude of the concentration field mode $a_i$ which is related to the amplitude of vertical velocity mode $A_i$ via ${\partial c}/{\partial t} = - u_3\varGamma$ (from (2.7b)), therefore

(2.15)\begin{equation} \sum_{i} \frac{\partial a_i}{\partial t} = \sum_{i} - A_i\varGamma. \end{equation}

Substituting $A_i$ from (2.15) in (2.14) and then integrating the resulting equation with respect to $t$ and rearranging the terms, we get

(2.16)\begin{equation} \frac{\partial^2 a_i}{\partial t^2} + \left( f^2 \left(\frac{N^2 - \varOmega_i^2}{N^2 - f^2}\right) + N^2 \left(\frac{\varOmega_i^2 - f^2}{N^2 - f^2}\right) (1+F \cos(\omega t)) \right) a_i = 0. \end{equation}

This system of equations is equivalent to a set of Mathieu equations (Kovacic, Rand & Mohamed Sah 2018) of characteristic frequency $\varOmega _i$, which governs the stability of our full inhomogeneous inviscid rotating system under vertically periodic forcing. Briard et al. (2020) obtained an equation similar to (2.16) in the absence of rotation ($\,f=0$) for the full inhomogeneous inviscid system. The discrete values of characteristic frequency $\varOmega _i$ can be determined for the cases when $N>f$, $N< f$ and $N=f$, as detailed below.

2.3. Characteristic frequency $\varOmega _i$ of the full inhomogeneous problem (2.16) when $N>f$

To solve (2.9) for $N>f$, we consider the walls of the domain at $x_3 = \pm H/2$. Therefore, the boundary conditions at $x_3 = \pm H/2$ are $\phi =0$. As we are solving (2.9) for the without periodic forcing case, we assume constant $\phi ={u_3}_{top}$ at $x_3=L_0/2$ (interface between upper lighter fluid and mixed fluid) and constant $\phi ={u_3}_{bot}$ at $x_3=-L_0/2$ (interface between mixed fluid and bottom denser fluid). Here $L_0$ represents the initial mixing zone size, and $H$ represents the domain height, as shown in the figure 1. Equation (2.9) has a finite number of solutions ($\phi _i$) corresponding to the characteristic frequency $\varOmega _i$ for a given horizontal wavenumber ($K$). The solutions of (2.9) (see Appendix D for the solution steps) in the upper pure fluid region ($x_3 \geq {L_0}/{2},\varGamma =0$), bottom pure fluid region ($x_3 \leq {-L_0}/{2},\varGamma =0$) and mixed fluid region ($\lvert x_3 \rvert \leq {L_0}/{2},\varGamma =-N^2/(2\mathcal {A}g_0)$), respectively, are

(2.17a)$$\begin{gather} \phi_i \left( x_3 \geq \frac{L_0}{2} \right) =-{u_3}_{top}\frac{\displaystyle \sinh{\left( K Y_i \left(x_3-\frac{H}{2} \right) \right)}}{\displaystyle \sinh{\left( K Y_i \left(\frac{H-L_0}{2} \right) \right)} }, \end{gather}$$

(2.17b)$$\begin{gather}\phi_i \left( x_3 \leq- \frac{L_0}{2} \right) =+{u_3}_{bot}\frac{\displaystyle \sinh{\left( K Y_i \left(x_3+\frac{H}{2} \right) \right)}}{\displaystyle \sinh{\left( K Y_i \left(\frac{H-L_0}{2} \right) \right)} }, \end{gather}$$

(2.17c)$$\begin{gather}\phi_i \left( \lvert x_3 \rvert \leq \frac{L_0}{2} \right) = \frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2\sin{\left(K X_i \frac{L_0}{2} \right)} } \sin{(K X_i x_3 )} + \frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos{\left(K X_i\frac{L_0}{2} \right)} } \cos{(K X_i x_3 )}, \end{gather}$$

where,

(2.18a,b)\begin{equation} X_i=\sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}, \quad Y_i=\sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 -f^2}}.\end{equation}

Both $\phi _i$ and $\partial _{x_3}\phi_i$ are continuous at ${x_3}=\pm L_0/2$. Therefore, we differentiate equation (2.17c) with respect to $x_3$ at ${x_3}=+L_0/2$ and ${x_3}=-L_0/2$, and equate it to the derivative of (2.17a) at ${x_3}=+L_0/2$ and derivative of (2.17b) at ${x_3}=-L_0/2$, respectively. After comparing the coefficients of ${u_3}_{top}$ and ${u_3}_{bot}$ we obtain the following equations (see Appendix E):

(2.19a)$$\begin{gather} \tan\left(K \sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}\frac{L_0}{2} \right) =- \sqrt{\frac{N^2-\varOmega_i^2}{\varOmega_i^2}} \tanh{\left(K \sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 -f^2}} \left(\frac{H-{L}_0}{2} \right) \right)}, \end{gather}$$

(2.19b)$$\begin{gather}\tan\left(K \sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}\frac{L_0}{2} \right) = \sqrt{\frac{\varOmega_i^2}{N^2-\varOmega_i^2}} \frac{1}{\displaystyle \tanh{\left(K \sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 -f^2}} \left(\frac{H-{L}_0}{2} \right) \right)}}. \end{gather}$$

If we assume the vertical height of the domain to be infinite, the terms with $\tanh {()}$ can be replaced by 1. We are interested in the solutions of (2.19a) and (2.19b) in terms of $\varOmega$ for a given value of $N$, $f$ and $KL_0$. In general, the solutions of (2.19a) and (2.19b) show that for given $KL_0$, the possible values of characteristic frequency $\varOmega _i$ fall between $f$ and $N$, where $N>f$. To demonstrate this, we assume $N=3$ and $f=1$ and the solutions show that $1<\varOmega _i<3$ (as illustrated in figure 2a). All the $\varOmega _i$ become asymptotic for $KL_0 \gg 1$ and $\varOmega _0 \approx N$ represents the fully developed mixing zone case corresponding to the limiting case of homogeneous stratified turbulence as reported by Gréa & Adou (2018) and Briard et al. (2020). However, for $KL_0 \ll 1$, from (2.19b), we obtain $\varOmega _0^2 \simeq f^2 + N^2 (KL_0/2) \simeq f^2 + \mathcal {A}g_0K$ (see details in Appendix F), which corresponds to the limiting case of the interface. Notably, for $f=0$, the frequency of interface reduces to $\varOmega _0 \simeq \sqrt {N^2 (KL_0/2)} \simeq \sqrt {\mathcal {A}g_0K}$ and the limiting case of the interface presented by Benjamin & Ursell (1954) and Briard et al. (2020) is recovered.

Figure 2. Possible values of characteristic frequency $\varOmega _i$ of the full inhomogeneous problem (2.16) as a function of horizontal wavenumber $K$ and initial mixing zone width $L_0$. (a) Here $\varOmega _i$, normalized by $N$, as a function of $KL_0$ for $N>f$, where $\varOmega _i$ are the solutions of (2.19a) and (2.19b) for rotating $\,f=1$ and non-rotating $f=0$ cases with $N=3$. Solutions of (2.19a) are odd and of (2.19b) are even. For the limiting case of the interface ($KL_0\ll 1$) plot of $\varOmega _0/N$ is shown in the inset. (b) Here $\varOmega _i$, normalized by $f$, as a function of $KL_0$, where $\varOmega _i$ are solutions of (2.27a) and (2.27b) with $f=3$, $N=1$ and $KH=100$. Plot of $\varOmega _i/f$ for $f=3$, $N=1$ and $KH=1000$ is shown in the inset. Solutions of (2.27a) are odd and of (2.27b) are even.

2.3.1. Linear theory for the limiting case of the fully developed mixing zone (when $N>f$)

We now simplify the full inhomogeneous inviscid problem (general Mathieu equation (2.16)) for the homogeneous case corresponding to the fully developed mixing zone where $KL_0\gg 1$, because all the $\varOmega _i$ values are constant. First, we rewrite the (2.9) as

(2.20)\begin{equation} \frac{\partial^2\phi_i}{\partial {x_3}^2} + m^2 \phi_i = 0, \end{equation}

where

(2.21)\begin{equation} m^2=\frac{K^2(N^2-\varOmega_i^2)}{(\varOmega_i^2-f^2)}. \end{equation}

According to the approximate Wentzel–Kramers–Brillouin (WKB) solution of (2.20) (see Appendix L for WKB approximation), $m$ represents the local wavenumber in the $x_3$ (vertical) direction. So, (2.21) represents the dispersion relation for the inertial–gravity waves (Smith & Waleffe 2002), and is given by

(2.22)\begin{equation} \varOmega_i^2 = \frac{K^2N^2+m^2f^2}{K^2+m^2} = \frac{K^2}{\lvert\boldsymbol{K}\rvert^2}N^2+ \frac{m^2}{\lvert\boldsymbol{K}\rvert^2}f^2, \end{equation}

where, $\lvert \boldsymbol {K}\rvert =\sqrt {K^2+m^2}$ is the magnitude of wavevector $\boldsymbol {K}$. We introduce $\tan (\theta _i)=K/m$, where $\theta _i$ is the angle between the vertical axis and the wavevector $\boldsymbol {K}$, with components $(k,l,m)$ in the $x_1$, $x_2$ and $x_3$ (vertical) directions, respectively. After substituting the definition of angle $\theta _i$ in (2.21), we obtain the characteristic frequency as follows (see Appendix G):

(2.23)\begin{equation} \varOmega_i^2=N^2\sin^2{(\theta_i)}+f^2\cos^2{(\theta_i)}. \end{equation}

The above relation shows that $\varOmega _i$ depends only on the direction of wavevector $\boldsymbol {K}$ and not on the magnitude of $\boldsymbol {K}$, which is consistent with the observation that all $\varOmega _i$ are constant for $KL_0\gg 1$. Now, we substitute the characteristic frequency $\varOmega _i$ from (2.23) in (2.16) to get (see Appendix H)

(2.24)\begin{equation} \frac{\partial^2 a_i}{\partial t^2} + (\,f^2 \cos^2{(\theta_i)} + N^2 \sin^2{(\theta_i)} ( 1 + F \cos{(\omega t)} )) a_i= 0 . \end{equation}

We define non-dimensional time as $\tau =\omega t$ and replace $\theta =\theta _i$ and $a=a_i$ such that (2.24) becomes

(2.25)\begin{equation} \frac{\partial^2 a}{\partial \tau^2} + \left( \frac{f^2 \cos^2{(\theta)}}{\omega^2} + \frac{N^2 \sin^2{(\theta)}}{\omega^2} (1+F\cos{(\tau)} ) \right) a = 0. \end{equation}

We solve this Mathieu equation (2.25) using Floquet theory, which provides the solutions to linear differential equations with periodic coefficients. Details of solving (2.25) using Floquet theory are in Appendix K. The resulting 3-D stability diagram for (2.25) is illustrated in figure 3. The stability diagram of the general equation (2.16) for the full inhomogeneous inviscid rotating system is similar to (2.25) except for the eigenvalues.

Figure 3. A 3-D stability diagram obtained by solving the Mathieu equation (2.25) in the parameter space $(N^2 \sin ^2(\theta )/\omega ^2,f^2 \cos ^2(\theta )/\omega ^2,F )$. All solutions inside the red and cyan coloured regions are unstable in nature with subharmonic and harmonic responses, respectively. Solutions are stable outside the unstable regions.

2.3.2. Linear theory for the limiting case of the interface (when $N>f$)

In the limiting case $KL_0 \ll 1$, where the frequency of the well-defined interface is $\varOmega _0^2 \simeq f^2 + N^2 (KL_0/2) \simeq f^2 + \mathcal {A}g_0K$, the full inhomogeneous inviscid problem (see (2.16)) reduces to the interface problem (see derivation in Appendix I), which is written as

(2.26)\begin{equation} \frac{\partial^2 a}{\partial \tau^2} + \left( \frac{f^2}{\omega^2} + \frac{\mathcal{A}g_0K}{\omega^2} (1+F \cos(\tau) ) \right)a = 0. \end{equation}

The stability diagram for the interface problem is obtained by solving (2.26) using Floquet theory (see Appendix K). We show this stability diagram in the $(\mathcal{A}g_0K/{\omega ^2}-F)$ plane at different ${f^2}/{\omega ^2}$ in a later figure and will discuss it in the results section (§ 4.3).

2.4. Characteristic frequency $\varOmega _i$ of the full inhomogeneous problem (2.16) when $N< f$

To solve (2.9) for $N< f$, we use the same boundary conditions and piecewise background concentration profile (for simplicity) as used for the previous case $N>f$. The final solutions ($\phi _i$) (see Appendix J for solution steps) are used to derive an equation for $\varOmega _i$ by ensuring that $\phi _i$ and $\partial _{x_3}\phi_i$ are continuous at ${x_3}=\pm L_0/2$ and using the same procedure as in the earlier case of $N>f$ (Appendix E). The resulting equations are written as

(2.27a)$$\begin{gather} \tan\left(K \sqrt{\frac{\varOmega_i^2 - N^2}{f^2 - \varOmega_i^2}}\frac{L_0}{2} \right) =- \sqrt{\frac{\varOmega_i^2 - N^2}{\varOmega_i^2}} \tan{\left(K \sqrt{\frac{\varOmega_i^2}{f^2 - \varOmega_i^2 }} \left(\frac{H-{L}_0}{2} \right) \right)}, \end{gather}$$

(2.27b)$$\begin{gather}\tan\left(K \sqrt{\frac{\varOmega_i^2 - N^2}{f^2 - \varOmega_i^2 }}\frac{L_0}{2} \right) = \sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 - N^2}} \frac{1}{\displaystyle \tan{\left(K \sqrt{\frac{\varOmega_i^2}{f^2 - \varOmega_i^2 }} \left(\frac{H-{L}_0}{2} \right) \right)}}. \end{gather}$$

We solve (2.27a) and (2.27b) to obtain the characteristic frequency $\varOmega _i$. Figure 2(b) illustrates $\varOmega _i$, normalized by $f$, as a function of $KL_0$ for $f=3$, $N=1$ and $KH=100$. The solutions $\varOmega _i$ at higher $KH$ ($=1000$) are shown in the inset of the figure 2(b). Note that here $KL_0< KH$. For all $KL_0$ each $\varOmega _i$ is constant with possible values between $N$ and $f$, i.e. $N<{\varOmega _i}< f$ or $0.333<{\varOmega _i}/{f}<1$. Moreover, for $KH=100$ and $1000$ each corresponding $\varOmega _i$ are approximately equal. The maximum characteristic frequency $\varOmega _0\sim f$ corresponds to the frequency of inertial waves in a stratified rotating fluid, when rotation is stronger and stratification is weak ($N< f$) (Ferrari & Wunsch 2009 and Davidson 2013, chapters 3 and 4). Therefore, the full inhomogeneous inviscid problem (2.16) can be simplified for the $N< f$ case because all the $\varOmega _i$ are constant for all $KL_0$ and $KH$. We rewrite the expression (2.21) for $N< f$ as follows:

(2.28)\begin{equation} m^2=\frac{K^2(\varOmega_i^2-N^2)}{(\,f^2-\varOmega_i^2)}, \end{equation}

where $m$ denotes the local vertical ($x_3$) wavenumber from WKB approximation (see Appendix L). Notably, the dispersion relation becomes the same as that for $N>f$ which is given by (2.22) and (2.23). The observation of the constant values of all $\varOmega _i$ for each $KL_0$ and $KH$ (see figure 2b) agrees with the independence of $\varOmega _i$ on the magnitude of wavevector and dependence on the direction ($\theta _i$) of the wavevector shown by the dispersion relation (2.23). Since the expression of $\varOmega _i$ is same for the $N< f$ and $N>f$ cases, (2.16) for the full inhomogeneous inviscid problem simplifies to (2.25). As a result, the 3-D stability diagram for $N< f$ is the same as that for the $N>f$ case.

2.5. Characteristic frequency $\varOmega _i$ of the full inhomogeneous problem (2.16) when $N=f$

When the stratification and rotation are of equal strength ($N=f$), the dispersion relation for the inertial-gravity waves (2.22) gives $\varOmega = N = f$. Here, we use the concept of group velocity $\boldsymbol {c}_g$, which is defined as a gradient of $\varOmega$ in the wavenumber space ($k,l,m$) and shows the velocity at which energy propagates away from a disturbance in the form of wave packets (Davidson 2013, chapters 3 and 4). Since $\varOmega = N = f$, $\boldsymbol {c}_g=\partial \varOmega /\partial \boldsymbol {K}_i=0$, indicating that the wave energy does not propagate or there is no dispersal of energy by waves, and there is no anisotropy in the large-scale eddies.

3. Numerical methodology and simulation parameters

The governing equations (2.1a), (2.1b) and (2.1c) are discretized using the finite-difference method on a staggered grid arrangement. We store the velocity fields at the cell faces and the pressure and concentration fields at the cell centres. A second-order central finite-difference scheme is employed to compute all the spatial derivatives. For the advancement in time, we use an explicit third-order Runge–Kutta method except for the diffusion terms, which are solved implicitly using the Crank–Nicolson method (Pal 2020; Pal & Chalamalla 2020). The pressure Poisson equation is employed to project a velocity field into a divergence-free space and is solved using a parallel multigrid iterative solver to obtain the dynamic pressure. This numerical solver has been extensively validated and used for several DNS of stratified free-shear and wall-bounded turbulent flows (Brucker & Sarkar 2010; Pal, de Stadler & Sarkar 2013; Pal & Sarkar 2015; Pal 2020; Pal & Chalamalla 2020; Naskar & Pal 2022a,b). To control the spurious reflections from the disturbances propagating out of the domain, we use sponge regions near the top and bottom boundaries, where damping functions gradually relax the values of the variables to their corresponding values at the boundary. These damping functions are added on the right-hand side of the momentum equation (2.1b) as explained (Brucker & Sarkar 2010). The sponge region is far away from the mixing zone such that it does not affect the dynamics of the mixing of fluids. Periodic boundary conditions are used in the $x_1$ and $x_2$ (horizontal) directions for all the variables. In the $x_3$ (vertical) direction, the top and bottom boundaries are walls, with Dirichlet and Neumann boundary conditions as follows:

(3.1)\begin{equation} U_1=U_2=U_3=0, \quad \frac{\partial P}{\partial x_3}=\frac{\partial C}{\partial x_3}=0, \text{at} \ x_3={\pm} H/2. \end{equation}

The initial concentration profile (Briard et al. 2019, 2020) in the present simulations is given as

(3.2)\begin{equation} C(\boldsymbol{x},t=0 )=\frac{1}{2} \left(1-\tanh{\left(\frac{2x_3}{\delta}\right)}\right), \end{equation}

where the parameter $\delta$ is used to change the initial mixing zone size-$L_0$. All of the cases simulated here have the same $\delta =0.035$ and corresponding $L_0=0.096$. Broadband fluctuations (Pal & Sarkar 2015) are imposed on the initial concentration profile (3.2), with an initial energy spectrum as follows:

(3.3)\begin{equation} E(\mathcal{K})=(\mathcal{K}/\mathcal{K}_0)^4\, \mathrm{e}^{-2(\mathcal{K}/\mathcal{K}_0)^2}, \end{equation}

where, $\mathcal {K}_0=30$. The fluctuations are localized at the centreline and are damped exponentially away from the centreline by multiplying the fluctuating concentration $c'$ with a damping function given by

(3.4)\begin{equation} g(x_3)=\mathcal{B}\,\mathrm{e}^{-\frac{1}{2}({x_3}/{\delta})^2}, \end{equation}

where $\mathcal {B}=0.7$ is the maximum intensity of fluctuation at the centreline. The size of the computational domain is $l_{x_1}=2{\rm \pi} \ \mathrm {m}$ and $l_{x_2}=2{\rm \pi} \ \mathrm {m}$ in the horizontal directions. The height of the domain is $l_{x_3}=2H=3.5{\rm \pi} \ \mathrm {m}$ which includes sponge region of thickness $H_s=0.64{\rm \pi} \ \mathrm {m}$ at top and bottom boundaries. We have used uniform grids $N_{x_1}=N_{x_2}=512$ in the horizontal direction, while in the vertical direction $x_3$, non-uniform grids $N_{x_3}=512$ are used. The grid is clustered at the vertical centre region of thickness $1.53{\rm \pi} \ \mathrm {m}$ with $\Delta x_{3_{min}}=\Delta x_1=\Delta x_2$, which is fine enough to resolve all the length scales in the mixing zone. We perform simulations with different parameters listed in table 1 to verify our theoretical solutions. The forcing frequency $\omega$ is chosen according to the saturation criterion (1.1) for all simulation cases with Coriolis frequency $f=0$, such that $L_{sat}\simeq 2.45$ at all forcing amplitudes $F$. We use the same $\omega$ for rotation cases at the corresponding $F$ to compare non-rotation and rotation cases.

Table 1. Parameters for the simulation: here $F$ is the forcing amplitude, $\omega$ is the forcing frequency, $f$ is the Coriolis frequency and $L_{end}$ is the final mixing zone size. The Atwood number $\mathcal {A}=0.01$, initial mixing zone width $L_0=0.096 \ \mathrm {m}$, gravitational acceleration $g_0=10\ \mathrm {m\ s}^{-2}$, kinematic viscosity $\nu =1\times 10^{-4}\ \mathrm {m}^2\ {\rm s}^{-1}$ and diffusion coefficient $\kappa =1\times 10^{-4}\ \mathrm {m}^2\ {\rm s}^{-1}$ is used for all of the cases (Briard et al. 2019).

4. Results

The 3-D stability diagram of our Mathieu equation (2.25) in the homogeneous limit corresponding to the fully developed mixing zone case for $N>f$ is shown in figure 3. The stability diagram (figure 3) is also valid for the $N< f$ case since the Mathieu equation (2.25) is the same for both $N>f$ and $N< f$ cases. The stability diagram shows the effect of the Coriolis force on the triggering and saturation/non-saturation of the instability through representative examples (line segments of different colours). The coloured curves correspond to the neutral stability curves, which separate the $(N^2 \sin ^2(\theta )/\omega ^2,f^2 \cos ^2(\theta )/\omega ^2 )$-plane into regions of stable and unstable (subharmonic or harmonic) solutions. Inside the red-coloured subharmonic instability tongues, the possible frequencies of the growing solutions are only odd multiple of the $\omega /2$, which is half of the external forcing frequency. The solutions inside the cyan-coloured harmonic instability tongues can have frequencies that are integer multiples of $\omega$. Other solutions that lie outside the instability tongues are stable.

4.1. Exploration of the stability diagram in the absence of the Coriolis force

Figure 4 shows the stability curves only for initial stratification conditions in the absence of the Coriolis force ( $f=0$) in the $(N/\omega )^2-F$ parameter plane. We obtained this stability curve similar to that of Gréa & Adou (2018) and Briard et al. (2020). The stable regions are separated by the unstable subharmonic (inside the red stability curve) and harmonic (inside the cyan stability curve) regions. Multiple frequencies are excited for a given initial stratification condition $N_0$ at a fixed $F$ and $\omega$, corresponding to an angle $\theta$ between $0$–${\rm \pi} /2$, as demonstrated by the horizontal brown line segment delineated by two crosses ($\times )$ in figure 4. This segment ranges from $0$ (left-hand end; $\theta =0$) to $(N_0/\omega )^2$ (right-hand end; $\theta ={\rm \pi} /2$), where $(N_0/\omega )^2$ represents the initial stratification condition corresponding to the initial mixing zone width-$L_0$ (because $(N/\omega )^2=2\mathcal {A}g_0/(L \omega ^2) \propto 1/L$, from (2.6)). When the initial condition $(N_0/\omega )^2$ is in the first stable region, the entire segment (shown with a diamond in figure 4) lies in the stable region. The corresponding $\theta$-modes are stable, and the instability does not trigger. However, if the initial conditions $(N_0/\omega )^2$ are in the first subharmonic tongue (segment with circles) the instability is triggered immediately by the unstable subharmonic $\theta$-modes, resulting in the rapid growth of $L$ and decrease in $(N/\omega )^2$, as indicated by the arrows pointing to the left. Therefore, in a non-rotating system $f=0$, Gréa & Adou (2018) deduced the following sufficient condition for the appearance of the instabilities:

(4.1)\begin{equation} \left(\frac{N}{\omega} \right)^2>\mathcal{G}(0,F). \end{equation}

Here $\mathcal {G}(0, F)=1/(2F+4)$ represents the marginal stability curve, shown by the black dashed curve in figure 4. This curve lies on the left-most boundary of the first subharmonic tongue. Finally, the mixing zone size saturates ($L_{sat}$) owing to the saturation of Faraday instability when the instability condition, (4.1), is no longer valid. The estimation of the saturated mixing zone size in this situation will be as follows (Gréa & Adou 2018):

(4.2)\begin{equation} \left(\frac{N}{\omega}\right)^2=\mathcal{G}(0,F)\ \mathrm{or}\ L_{sat}=\frac{2\mathcal{A}g_0}{\omega^2}\frac{1}{\mathcal{G}(0,F)}. \end{equation}

The black squares indicate the $L_{sat}$ on the marginal stability curve. When the initial condition is in the first harmonic tongue (segment with squares) or far away ($(N_0/\omega )^2 > \mathcal {G}(0, F)$), the mixing zone size grows due to both subharmonic and harmonic instabilities, and $(N/\omega )^2$ decreases. The transition from harmonic to subharmonic occurs when no more harmonic $\theta$-modes are excited, such that $(N/\omega )^2$ will cross the left stability curve of the first harmonic tongue. The corresponding mixing zone size is $L_{tr}$ and is demonstrated by the black circle (see figure 4). We will also show this transition using numerical simulations in section § 4.4. Afterwards, the mixing zone size saturates at $L_{sat}$ (as shown by the black square on the marginal stability curve) when the instability condition 4.1 is no longer valid. Therefore, the mixing zone size-$L$ saturates for all cases without rotation.

Figure 4. Stability diagram of the Mathieu equation (2.25) for stratification ($N$) initial condition in the absence of rotation (Coriolis frequency, $f=0$). The red and cyan-coloured stability curves separate the stable and unstable (subharmonic and harmonic) regions. The brown horizontal line segment represent the frequencies excited, corresponding to an angle $\theta$ between $0$–${\rm \pi} /2$, for the initial condition $(N_0/\omega )^2$ (right-hand end ($\times$) at $\theta ={\rm \pi} /2$) in the different tongue-like zones. Arrows indicate the evolution of $(N/\omega )^2$ as the mixing zone size-$L$ increases.

4.2. Exploration of the stability diagram in the presence of the Coriolis force

Figure 5 shows the top view of the 3-D stability diagram (figure 3) at forcing amplitude $F=1$. Various frequencies are excited for given $f$ and $(N_0/\omega )^2$, corresponding to an angle $\theta$ between 0–${\rm \pi} /2$, as represented by the inclined segments in figure 5. In contrast to the horizontal segment in the ‘without rotation’ cases, the inclined segment in the presence of rotation ranges from $[(N/\omega )^2,0]$ to $[(N/\omega )^2,(\,f/\omega )^2]$ for a given value of $f$. The right-hand end ($\times$) of the inclined segment (at $\theta ={\rm \pi} /2$) corresponds to the initial stratification condition $(N_0/\omega )^2$, while the left-hand end ($\times$) (at $\theta =0$) corresponds to the given $f$.

Figure 5. Top view of the 3-D stability curve (figure 3) at forcing amplitude $F = 1$. The inclined segments represent the frequencies excited, corresponding to an angle $\theta$ between $0$ to ${\rm \pi} /2$, for $(N/\omega )^2$ (right end ($\times$) at $\theta ={\rm \pi} /2$) and $(\,f/\omega )^2$ (left end ($\times$) at $\theta =0$). The initial conditions $(N_0/\omega )^2$ for given $(\,f/\omega )^2$ are in the first left-most (a) stable region and subharmonic tongue, and (b) harmonic tongue. Arrows indicate the evolution of $(N/\omega )^2$ as the mixing zone size-$L$ increases.

Similar to the ‘without rotation’ case, we define the condition for the occurrence of instabilities in the presence of rotation as follows:

(4.3)\begin{equation} \left(\frac{N}{\omega} \right)^2>\mathcal{G}(\,f^2/\omega^2, F), \end{equation}

where $\mathcal {G}(\,f^2/\omega ^2, F)$ represents the extreme left-hand boundary of the first subharmonic tongue. Beyond $\mathcal {G}$, Faraday instability is triggered by the unstable subharmonic or harmonic $\theta$-modes, resulting in the growth of $L$ and decrease of $(N/\omega )^2$. To understand the stability diagram in the presence of rotation, we study the following possible configurations based on the initial condition $(N_0/\omega )^2$ and a given $(\,f/\omega )^2$.

1. (i) If the initial condition is in the first left-hand most stable region such that $(N_0/\omega )^2 < \mathcal {G}(\,f^2/\omega ^2, F)$, then two scenarios, namely $(\,f/\omega )^2<0.25$ and $(\,f/\omega )^2 \geq 0.25$, are possible. For $(\,f/\omega )^2<0.25$, all the excited frequencies lie in the stable region, and the corresponding $\theta$-modes are stable, as depicted by the inclined pink line segment in figure 5(a). Therefore, the instability is not triggered, and the mixing zone size will remain at its initial value of $L(t=0)$. Note that we have defined the instability condition, (4.3), based on the fact that if at least one $\theta$-modes fall in any unstable tongues solution will grow, resulting in instability. Therefore, we define the saturation state of the Faraday instability when the instability condition (4.3) is no longer valid and $(\,f/\omega )^2<0.25$ as follows:

   (4.4) \begin{gather} \left(\frac{N}{\omega}\right)^2=\mathcal{G}(\,f^2/\omega^2, F)\ \mathrm{or}\ L_{sat}=\frac{2\mathcal{A}g_0}{\omega^2}\frac{1}{\mathcal{G}(\,f^2/\omega^2, F)},\quad (\mathrm{when}\ (\,f/\omega)^2<0.25). \end{gather}
   For $(\,f/\omega )^2 \geq 0.25$, some of the excited frequencies are in the first unstable subharmonic tongue (red), and the corresponding $\theta$-modes are unstable as demonstrated by the dark blue line segment (figure 5a). Therefore, instabilities trigger, $L$ grows, and $(N/\omega )^2$ decreases as indicated by the arrows directing to the left. The mixing zone size-$L$ will never reach a saturation state because, at every instant, some frequencies will always be excited in the subharmonic region, and the corresponding $\theta$-modes will be unstable, as depicted by the dashed blue segment in figure 5(a). Therefore, the instability condition 4.3 is always satisfied, and $L$ will continue to grow. However, the growth rate of $L$ might be slow at later stages because the $\theta$-modes that fall in the first unstable tongue decrease as $(N/\omega )^2$ decreases.

2. (ii) The instability condition, (4.3), is always satisfied when $(N_0/\omega )^2$ is chosen to be in the first subharmonic (red) region, and therefore, instabilities are triggered by the unstable $\theta$-modes. For $(\,f/\omega )^2<0.25$, the most unstable $\theta$-modes are excited in the first subharmonic region, as depicted by the solid orange segment in the figure 5(a). Therefore, $L$ grows and $(N/\omega )^2$ decreases. The instability will achieve a state of saturation when no subharmonic unstable $\theta$-modes are excited, i.e. the instability condition, (4.3), is no longer satisfied. The corresponding mixing zone size, $L_{sat}$, is given by the saturation condition, (4.4). Now $L_{sat}$ may lie anywhere on the left stability curve of the first subharmonic tongue, as shown by the orange ellipse (see figure 5a). For $(\,f/\omega )^2 \geq 0.25$, the solid green segment (see figure 5a) represents the initially excited frequencies. For this case, the mixing zone size-$L$ will never reach a saturation state because some unstable $\theta$-modes are always excited in the subharmonic region, as demonstrated by the dashed green segment.

3. (iii) Now we consider the case when the initial condition $(N_0/\omega )^2$ is beyond the first subharmonic tongue. For the initial condition in the first harmonic tongue (cyan), the grey segment (see figure 5b) represents the case for $(\,f/\omega )^2<0.25$. Owing to the validity of (4.3), both harmonic and subharmonic instabilities are triggered by the unstable $\theta$ -modes. For some initial time, the most unstable $\theta$-modes ($\theta \simeq {\rm \pi}/2$) are excited in the harmonic region, so the mixing zone size-$L$ grows due to harmonic instability and thus $(N/\omega )^2$ decreases. The growth of $L$ is followed by the harmonic to subharmonic transition, and the corresponding mixing zone size $L_{tr}$ is shown by the solid grey ellipse on the left boundary of the first harmonic tongue (see figure 5b). Afterwards, $L$ grows due to subharmonic instability if $\theta$-modes are excited in the first unstable subharmonic region. We will also demonstrate the harmonic to subharmonic transition, followed by a subharmonic instability, using numerical simulations in section § 4.4. Finally, the mixing zone size saturates at $L_{sat}$ (shown by a dashed grey ellipse) when the saturation condition, (4.4), is satisfied. For $(\,f/\omega )^2 \geq 0.25$, only the harmonic instability is triggered, as shown by the solid black segment in figure 5(b). After some time, during the growth of $L$, the $(N/\omega )^2$ can be in the stable region between the first subharmonic and harmonic tongues, as shown by the dash–dotted black segment (figure 5b). No unstable $\theta$-modes are excited in this stable region, hence $L$ will not grow. Therefore, subharmonic instability will not be triggered in the present case. This suggests that the instability condition 4.3 is not valid if the entire inclined segment lies in any of the stable regions beyond the first subharmonic tongue for a given $(\,f/\omega )^2$ and $(N/\omega )^2$. Now, we recall that the present theoretical model does not include the diffusion term. However, in realistic configurations, $L$ can diffuse and grow owing to the molecular diffusion resulting in a decrease in $(N/\omega )^2$. As a result, the right-hand end of the dash–dotted black segment $(N/\omega )^2$ enters the first subharmonic tongue. The subharmonic instability then leads to the growth of $L$. However, the mixing zone size-$L$ never achieves a state of saturation because, at all times, some frequencies are excited in the subharmonic region as depicted by the dashed black segment in figure 5(b).

Figures 6(a), 6(b) and 6(c) illustrate the stability diagram of the Mathieu equation (2.25) at $F=0.75, 2$ and $3$, respectively. We can observe that the first subharmonic (red) and harmonic (cyan) tongues become wide with increasing forcing amplitude $F$ (see comparison among figure 6a–c). Therefore, the right-hand boundary of the first subharmonic tongue (red) shifts towards the right-hand side, and the left-hand boundary of the first harmonic tongue (cyan) moves towards the left-hand side. As a result, the stable region between these tongues shrinks. As $F$ increases more subharmonic frequencies are excited and the possibility that the frequencies fall in the stable region between the first subharmonic and harmonic tongues decreases. Consequently, at higher $F$, the instabilities are always triggered for all values of $(\,f/\omega )^2$ and $(N/\omega )^2$ satisfying the instability condition, (4.3), unlike the lower $F$ cases where the entire segment can lie in the large stable region between the first subharmonic and harmonic tongues. Additionally, for all $(\,f/\omega )^2$, the excitation of more subharmonic frequencies triggers the instability quickly, resulting in a rapid increase in $L$. This signifies that for higher $F$, the effect of rotation on triggering the instabilities is negligible. The numerical simulations presented in section § 4.4 will further support these findings.

Figure 6. Top view of the 3-D stability curve (figure 3) at three different values of the forcing amplitude (a) $F=0.75$, (b) $F=2$ and (c) $F=3$. The green segment shown in (a) represents the frequencies excited for the ‘without rotation’ case $(\,f/\omega )^2=0$ (left end ($\times$) at $\theta =0$) with the initial condition $(N_0/\omega )^2$ (right end ($\times$) at $\theta ={\rm \pi} /2$) in the first harmonic tongue.

To summarize, the mixing zone saturates at $L_{sat}$ for $(\,f/\omega )^2<0.25$, whereas for $(\,f/\omega )^2 \geq 0.25$, the mixing zone never saturates and continues to grow due to the excitation of unstable subharmonic $\theta$-modes at all times.

4.3. Effect of Coriolis force on the onset of Faraday instability from a sharp interface

Figure 7(a) demonstrates the stability diagram of the interface problem given by (2.26) in the ($\mathcal {A}g_0K/\omega ^2,F$)-plane. The marginal stability curves are shown for different values of rotation rates $f^2$, normalized by $\omega ^2$, ranging from 0 to 0.3. Corresponding to $f^2/\omega ^2$, the region bounded by the solid curves represents the subharmonic ($SH$) instability tongue, whereas the region bounded by the dashed curves represents the harmonic ($H$) instability tongue. The stable regions are outside the instability tongues, which are separated by the solid and dashed neutral stability curves. The growth rates ($\sigma \neq 0$) of unstable modes for $f^2/\omega ^2=0,0.1$ and $0.2$ inside the first subharmonic tongue are depicted in the figure 7(b). These growth rates ($\sigma$) are determined using Floquet analysis (see Appendix K). Corresponding to $f^2/\omega ^2$ the solid lines represent $\sigma = 0$, whereas dash–dotted lines represent $\sigma \neq 0$ where $\sigma$ ranges from 0.02 to 0.14 with an interval of 0.02. Inside the instability tongues, the well-defined interface between two fluids is unstable to any perturbation resulting in the generation of standing waves with wavelength $\lambda$ ($=2{\rm \pi} /K$). For example, at $F=0.75$ for $f^2/\omega ^2=0$, the subharmonic waves can be excited with wavelength $\lambda ^{sh}$ between $\lambda ^{sh}_{min}$ and $\lambda ^{sh}_{max}$ as indicated by the black circles on the horizontal orange dashed line in figure 7(a). We can observe that the first subharmonic tongue moves upwards at higher rotation rates $f^2/\omega ^2 \geq 0.25$. Consequently, the subharmonic Faraday waves are not triggered at lower forcing amplitudes $F<2$ and $f^2/\omega ^2 \geq 0.25$, signifying that the Coriolis force suppresses the formation of subharmonic waves at relatively higher rotation rates. However, the subharmonic instability can always be excited by increasing the forcing amplitude above a critical value $F_c$ which is defined as the minima of the first subharmonic tongue for a given $f^2/\omega ^2$. In contrast, harmonic waves can always be triggered because harmonic tongues do not move upwards in the presence of Coriolis force, as shown in figure 7(a).

Figure 7. Panel (a) shows the stability diagram of the interface problem (Mathieu equation) (2.26) in the parameter space $(\mathcal {A}g_0K/\omega ^2,F )$ for different rotation rates $f^2/\omega ^2$. The tongue-like region bounded by the solid curves of respective $f^2/\omega ^2$ correspond to the subharmonic ($SH$) instability, whereas the region bounded by the dashed curves correspond to the harmonic ($H$) instability. The stable solutions lies outside the $SH$ and $H$ instability tongues. Panel (b) is an enlargement of the first left-most subharmonic tongues for $f^2/\omega ^2=0, 0.1$ and $0.2$ in (a). For the corresponding $f^2/\omega ^2=0$, solid curves denote the neutral stability curve for the growth rate $\sigma =0$, whereas dash–dotted curves denote the contours of different $\sigma$ at intervals of 0.02 ranging from 0.02 to 0.14. Arrow lines indicate the increasing order of $\sigma$. Panel (c) demonstrates the shrinkage of the first subharmonic tongue with increasing $f^2/\omega ^2$ in the ( $f^2/\omega ^2 + \mathcal {A}g_0K/\omega ^2,F$)-plane.

Figure 7(a) demonstrates that the first subharmonic tongue shifts leftwards as $f^2/\omega ^2$ increases from $0$ to $\leq 0.25$. In addition, all inside U-tongues for growth rates ranging from $\sigma = 0.02$ to $0.14$ shifts upwards as $f^2/\omega ^2$ increases from 0 to 0.2 (see figure 7b). For example, the minima of U-tongue for $\sigma = 0.02$ moves from $F\approx 0.162$ to $0.27$ and $0.76$ when $f^2/\omega ^2$ increases from $f^2/\omega ^2=0$ to 0.1 and 0.2, respectively, as indicated by the orange squares in figure 7(b). At constant forcing amplitude $F$, the growth rate of the most unstable mode decreases with an increase in rotation rate from $f^2/\omega ^2=0$ to $<0.25$. For example, at $F=0.75$ indicated by the horizontal dashed orange line in figure 7(b), the growth rate of the most unstable mode is $\sigma \approx 0.095$, $0.06$ and $0.02$ for $f^2/\omega ^2=0$, 0.1 and 0.2, respectively. These results suggest that at a fixed $F$, the Coriolis force suppresses the excitation of the subharmonic waves. As a result, the subharmonic waves will take time to develop and grow, thereby delaying the onset of the subharmonic instability. We conclude that at lower forcing amplitudes, the subharmonic instability suppressing and delaying effect of the Coriolis force increases with an increase in rotation rates and eventually becomes so strong that no subharmonic instability is triggered at higher rotation rates, i.e. when $F<2$ and $f^2/\omega ^2 \geq 0.25$. However, the subharmonic instability can always be excited by increasing $F$ above a critical value $F_c$. Here we can recall that if molecular diffusion exists, the subharmonic instability can also be triggered from the diffuse interface when $F<2$ and $f^2/\omega ^2 \geq 0.25$, but not from the sharp interface, as discussed in the previous section (§ 4.2).

Additionally, the first subharmonic tongue and the inside U-tongues shrink when $f^2/\omega ^2$ increases from 0 to 0.2 (see figure 7b). This shrinkage can be clearly observed in the $(\,f^2/\omega ^2 + \mathcal {A}g_0K/\omega ^2,F)$-plane shown in figure 7(c). Notably, the shrinkage with increasing $f^2/\omega ^2$ is significant at lower $F$, and becomes insignificant at higher $F$. This is illustrated by the following example: in comparison with $f^2/\omega ^2=0$ the shrinkage of subharmonic tongue for $f^2/\omega ^2=0.2$ is $\approx 80\,\%$ at $F=0.75$ and $\approx 20\,\%$ at $F=3$. It suggests that fewer most unstable modes are available at lower $F$ than at higher $F$, which can excite and grow the subharmonic waves. According to the previous discussion, at fixed $F$, the growth rate of the most unstable mode is less for $f^2/\omega ^2=0.2$ than for $f^2/\omega ^2=0$, resulting in the suppression and delay of the subharmonic instability by the Coriolis force. This observation is only true at lower $F$ owing to the additional fact that shrinkage of the subharmonic tongue results in the availability of fewer unstable modes for $f^2/\omega ^2=0.2$ than that for $f^2/\omega ^2=0$ which are also responsible for the excitation and growth of the subharmonic waves. In contrast to the lower $F$, a large number of most unstable modes are available at higher $F$, nearly equivalent to that for $f^2/\omega ^2=0$, that can rapidly grow the subharmonic waves. Therefore, we can deduce that at higher forcing amplitudes, the flow is energetic enough to overcome the efficacy of rotation in suppressing and delaying the onset of subharmonic instability.

For each $f^2/\omega ^2$, the U-tongues for different $\sigma$ move upwards as $\sigma$ increases, as indicated by the arrows in figure 7(b). For $f^2/\omega ^2=0$ and at lower forcing amplitudes approximately $F\leq 1$, the minima of all U-tongues deviate slightly from $\mathcal {A}g_0K/\omega ^2=0.25$ as shown by a vertical dashed black line in figure 7(b). As a result, for the ‘without rotation’ case $f^2/\omega ^2=0$, at a given $F$ the wavenumber $K^{sh}$ corresponding to the growth rate of the most unstable mode at the onset of subharmonic instability can be approximated as $K^{sh}=\omega ^2/(4\mathcal {A}g_0)$ (Briard et al. 2020). In the presence of rotation, the first subharmonic tongue ($\sigma =0$) and the inside U-tongues ($\sigma \neq 0$) move leftwards with an increase in $f^2/\omega ^2$ until $\leq 0.25$. Similar to the case $f^2/\omega ^2=0$, for approximately $F\leq 1$ the minima of U-tongues depart slightly from $\mathcal {A}g_0K/\omega ^2=0.15$ and $0.05$ for $f^2/\omega ^2=0.1$ and $0.2$ cases, respectively. Consequently, for $f^2/\omega ^2\leq 0.25$, the approximation of $K^{sh}$ corresponding to maximum growth rate can be written as

(4.5)\begin{equation} K^{sh}=\frac{\omega^2}{\mathcal{A}g_0}\left(\frac{1}{4}- \left(\frac{f}{\omega}\right)^2\right). \end{equation}

Substituting $\lambda ^{sh}=2{\rm \pi} /K^{sh}$ in (4.5), we get

(4.6)\begin{equation} \lambda^{sh}=\frac{2{\rm \pi}\mathcal{A}g_0}{\omega^2} \left(\frac{1}{{1}/{4}-\left(\,{f}/{\omega}\right)^2}\right). \end{equation}

Therefore, this approximation (4.6) of the wavelength of the initial most amplified mode starting from the well-defined sharp interface signifies that $\lambda ^{sh}$ increases (or $K^{sh}$ decreases) with an increase in $f^2/\omega ^2(\leq 0.25)$. With further increase in $f^2/\omega ^2(>0.25)$, $\lambda ^{sh}$ starts to decrease (or $K^{sh}$ increase) since the first subharmonic tongue begins to shift rightwards (see figure 7a). Similar effects of the Coriolis force on the onset of subharmonic waves in a thin sheet of viscous fluid were reported by Mondal & Kumar (2004), who performed a linear stability analysis by solving the set of Mathieu equations using Floquet theory.

4.4. Numerical results

In the analytical model, the instability of the interface under the influence of rotation is presented in terms of the mixing zone size-$L(\propto \omega ^2/N^2)$. We present the time evolution of the mixing zone size-$L$ (see (2.5)) for forcing amplitude $F=0.75$ in figure 8(a) obtained from our simulations. The initial stratification condition corresponding to the initial mixing zone size-$L_0=0.096$ is $N_0^2/\omega ^2=4.64$ that satisfies the instability condition (4.3) and ensures the onset of the Faraday instability. For the case without rotation (F075f/$\omega$0) at $\omega t\simeq 58$, $L\simeq 0.48$ and begins to oscillate with a small amplitude that rapidly grows (shown in the inset of figure 8a indicated by black arrow) owing to the subharmonic instability resulting in turbulent mixing. According to the stability diagram in figure 6(a) (green line segment), the $N^2/\omega ^2\simeq 0.93$ corresponding to $L\simeq 0.48$ lies in the stable region between the first subharmonic and harmonic tongues. When $L$ increases, $N^2/\omega ^2$ decreases and enters the first subharmonic tongue resulting in the rapid growth of $L$ and turbulent mixing. After $\omega t\simeq 130$, an asymptotic state is achieved, where the mixing zone size oscillates without any further increase, signifying the saturation of instability. The final size of the mixing zone is $L_{end}\simeq 3.1$ for F075f/$\omega$0. In the stability diagram (figure 6a, green line segment), $N^2/\omega ^2\simeq 0.144$ corresponding to $L_{end}\simeq 3.1$ lies in the first stable region and signifies the saturation state of the instability.

Figure 8. Evolution of the mixing zone size-$L$ with non-dimensional time $\omega t$ for without rotation and rotation cases (a) at $F=0.75$, (b) at $F=1$, (c) at $F=2$ and (d) at $F=3$. Horizontal dashed black lines represent $L_{sat}$ calculated from the prediction (1.1) for $f=0$ cases. Inset labelled by transition region is the enlarged view of small oscillation before the harmonic to subharmonic transition $L_{tr}$. Insets indicated by arrows correspond to the onset of subharmonic instability (shown by square box).

The presence of rotation at $F = 0.75$ impedes the growth of the instabilities. For the case with $f/\omega =0.48$ the subharmonic instability is triggered at $\omega t\simeq 340\unicode{x2013}355$, $L\simeq 1.08\unicode{x2013}1.12$, $N^2/\omega ^2\simeq 0.412\unicode{x2013}0.398$, as illustrated by the case F075f/$\omega$48 in figure 8(a) (see the inset indicated by a green arrow). When the Coriolis frequency is increased to $f/\omega =0.59$, the subharmonic instability begins at $\omega t\simeq 635\unicode{x2013}650$, $L\simeq 1.48\unicode{x2013}1.52$, $N^2/\omega ^2\simeq 0.3\unicode{x2013}0.293$, as depicted by the case F075f/$\omega$59 in figure 8(a) (see the inset indicated by a red arrow). Therefore, the onset of the subharmonic instability causing turbulent mixing is significantly delayed owing to rotation. The large Coriolis frequency ( $f/\omega =0.59$) also suppresses the amplitude of oscillations in $L$ in the subharmonic instability phase. For $f/\omega =0.48$ ( $f^2/\omega ^2=0.23$) the instability saturates ($L_{end}\simeq 3.25$) and small oscillations in the asymptotic state are damped out by diffusion, whereas for $f/\omega =0.59$ ( $f^2/\omega ^2=0.35$) the instability never saturates and hence the mixing zone size-$L$ continues to evolve with oscillations. This result is consistent with our theoretical prediction that the instability saturates for $f^2/\omega ^2<0.25$ but never saturates for $f^2/\omega ^2\geq 0.25$.

At forcing amplitude $F=1$, the time evolution of $L$ for ‘without rotation’ (F1f/$\omega$0) and ‘with rotation’ (F1f/$\omega$48 and F1f/$\omega$59) cases are illustrated in figure 8(b). For all the cases, the instability condition, (4.3), is satisfied by the initial stratification condition $N_0^2/\omega ^2=4.25$ corresponding to the initial mixing zone size-$L_0=0.096$ that ensures the onset of the Faraday instability. Before the rapid growth of $L$ small oscillations of varying periods (shown in the inset labelled by a transition region of figure 8b) is observed for both the ‘without rotation’ ( $f/\omega =0$) and ‘with rotation’ ( $f/\omega =0.48$ and $f/\omega =0.59$) cases. These small oscillations represent the harmonic phase following which the harmonic to subharmonic transition (Briard et al. 2019) occurs. For $f/\omega =0.48$ and $0.59$, this transition occurs at $\omega t\simeq 60\unicode{x2013}70$ and $65\unicode{x2013}75$ (see inset of figure 8b labelled by transition region) followed by the stable region until $\omega t\simeq 150\unicode{x2013}160$ and $290\unicode{x2013}300$, respectively. No harmonic or subharmonic modes are excited in the stable region, and $L$ grows due to diffusion. We also observe similar stable regions for F075f/$\omega$48 and F075f/$\omega$59 in figure 8(a) (see the Supplementary movies 1 and 2 are available at https://doi.org/10.1017/jfm.2023.744 of the concentration field for $F = 0.75$ and $1$). The presence of this stable region delays the onset of the subharmonic instability. This verifies our theoretical prediction in § 4.2 that no instability is triggered if the excited frequencies lie in the stable region between the first subharmonic and harmonic tongues for $f^2/\omega ^2\geq 0.25$. However, the presence of diffusion causes $L$ to grow and eventually subharmonic instability triggers. Based on the numerical results we can also conclude that the same observation is also true for $f^2/\omega ^2<0.25$ cases if all the modes are excited in the stable region. Here $L$ rapidly grows after the triggering of the subharmonic instability resulting in turbulent mixing. Similar to $F = 0.75$, we observe that $L$ saturates for F1f/$\omega$48, but continues to grow for F1f/$\omega$59.

The effect of rotation on the evolution of $L$ at higher forcing amplitudes $F=2$ and $F=3$ is demonstrated in figures 8(c) and 8(d), respectively. The initial stratification condition is $N_0^2/\omega ^2=3.26$ for $F=2$ and $N_0^2/\omega ^2=2.6$ for $F=3$ corresponding to the initial mixing zone size-$L_0=0.096$ which ensures the onset of Faraday instability. At $F=2$, the harmonic to subharmonic transition ($L_{tr}\simeq 0.5$) occurs at $\omega t\simeq 13$ for F2f/$\omega$0, F2f/$\omega$48 and F2f/$\omega$59 as indicated by the black circle in the inset of figure 8(c). Further increasing the forcing amplitude to $F=3$, the transition ($L_{tr}\simeq 0.375$) occurs at $\omega t\simeq 6.5$ for the F3f/$\omega$0, F3f/$\omega$48 and F3f/$\omega$59 cases, as shown by the black circle in the inset of figure 8(d). The instability delaying effect of rotation is subdued as the forcing amplitude increases. This observation is consistent with our analytical model that predicts the shrinkage of the stable region between the unstable harmonic and subharmonic tongues with increasing forcing amplitudes (see stability diagrams at $F = 2, 3$ in figures 6b and 6c, respectively). For cases with $f^2/\omega ^2=0$ and $f^2/\omega ^2=0.23$, the instability saturates, and the oscillations in the asymptotic state are damped. However, the instability never saturates for cases with $f^2/\omega ^2=0.35$ and $L$ continues to evolve with oscillations. Note that at each corresponding forcing amplitude $F=0.75,1,2$ and $3$, the values of the final mixing zone size $L_{end}$ for the corresponding rotation cases $f/\omega =0.48$ ( $f^2/\omega ^2<0.25$) are very close to the non-rotation cases ( $f/\omega =0$). We report $L_{end}$ values in table 1. Therefore, for the parameters used in the present simulations, the main difference between the ‘with rotation’ and ‘without rotation’ cases at each $F$ comes principally from the growth of the instability.

Fourier transforms of the derivatives of $L$ for different time intervals, each starting from $t=0$, are performed to capture the harmonic to subharmonic transition in our simulations. For case F1f/$\omega$0, the normalized Fourier transform of ${\rm d}L/{\rm d}t$ as a function of the normalized frequency $f/(\omega/2{\rm \pi} )$ is depicted in figure 9(a–d). Two peaks of different frequencies are present in $L$. The first and the second peaks represent the subharmonic and harmonic regimes, respectively. We can observe that for $\omega t \leq 45$, the second peak is higher (see figure 9a), denoting the dominance of the harmonic regime. At $\omega t = 45.92$ (figure 9b) the amplitude of both the peaks are similar (see figure 9b) indicating the transition from harmonic to subharmonic regime. Note that the small oscillations in $L$ also rapidly grow from $\omega t\simeq 45$ (see figure 8b) owing to the onset of the subharmonic instability. When $\omega t > 45$ the first peak becomes the dominating peak as illustrated in figures 9(c) and 9(d), signifying the primary role of the subharmonic instabilities in turbulent mixing. The Fourier transform of ${\rm d}L/{\rm d}t$ for F1f/$\omega$48 demonstrates the harmonic to subharmonic transition, stable region and the onset of the subharmonic instability. The dominating second peak at $\omega t = 65.4$ in figure 9(e) shows small oscillations of $L$ (see inset of figure 8b indicated by a green arrow) are in the harmonic region. The amplitude of the first and second peak remains approximately unchanged until $\omega t\simeq 150$ (figures 9f and 9g) which represents the stable region. At $\omega t\simeq 150\unicode{x2013}160$, the first peak starts increasing, and the second peak starts decreasing. Eventually the first peak becomes the dominant (figure 9h) marking the onset of the subharmonic instability at $\omega t\simeq 150\unicode{x2013}160$. An increase in the Coriolis frequency to $f/\omega =0.59$ impedes the harmonic to subharmonic transition, expands the stable region, and delays the onset of the subharmonic instability (figures not shown).

Figure 9. Normalized Fourier transform of ${\dot L}$ as a function of normalized frequency, (a–d) for F1f/$\omega 0$ case, (e–h) for F1f/$\omega 48$ case, for different time intervals each starting from $\omega t=0$ to $\omega t$ mentioned on the top of each figure. The first peak at $f /(\omega /2 {\rm \pi})\simeq 1$ corresponds to the subharmonic regime and second peak at $f /(\omega /2 {\rm \pi})\simeq 2$ corresponds to the harmonic regime (Briard et al. 2019).

Figure 10 depicts the 3-D concentration field in the mixing zone for ‘without rotation’ F1f/$\omega$0 and ‘with rotation’ F1f/$\omega$48 cases. The random wave patterns on the interface at $\omega t=41.5$ for F1f/$\omega$0 indicate the inception of the Faraday instability as shown in figure 10(a). The turbulent mixing caused by the subharmonic instability is demonstrated at $\omega t=70.6$ for F1f/$\omega$0 in figure 10(b) followed by the beginning of instability saturation and fully saturated instability in figures 10(c) and 10(d) at $\omega t=120$ and $\omega t\sim 152$, respectively. For F1f/$\omega$48, the stable region after the harmonic to subharmonic transition is depicted by figure 10(f) at $\omega t=120$. The subsequent figures 10(g) and 10(h) for F1f/$\omega$48 at $\omega t \sim 212$ and $329$ show a significant delay in turbulent mixing and saturation of the instability with respect to F1f/$\omega$0. Movies (1)–(4) of the concentration field at the $x_2 = {\rm \pi}$ plane for different forcing amplitudes are included as Supplementary movies, to demonstrate the influence of rotation on the evolution of the mixing zone.

Figure 10. Visualization of 3-D concentration field, (a–d) for F1f/$\omega$0 case, (e–h) for F1f/$\omega$48 case, at different time instants. Colour map: blue for $C=0$ (lighter fluid) and red for $C=1$ (denser fluid). For better visualization of the mixing of fluids, the pure fluids are made transparent.

5. Conclusions

The effect of rotation on the onset and saturation of Faraday instability in two miscible fluids subjected to vertical oscillations is investigated by means of theoretical analysis and DNS. A linear stability analysis is performed in the presence of rotation, yielding a set of Mathieu equations for the full inhomogeneous inviscid rotating problem. These equations are investigated further for the limiting case of the interface and homogeneous case when $N>f$. The $N< f$ case analysis yields the same Mathieu equations as the homogeneous case for $N>f$. The analytical model utilizes Floquet theory to solve a set of Mathieu equations derived from linear stability analysis. On solving the set of Mathieu equations, a 3-D stability diagram is obtained that consists of stable and unstable harmonic and subharmonic regions. In the limiting case of interface, we found that increasing rotation rates suppress and delay the onset of subharmonic instability at lower forcing amplitudes ($F$). Eventually, no subharmonic instability is triggered at relatively strong rotation rates. However, the subharmonic instability can always be excited by increasing $F$ above a critical value $F_c$. In contrast, the instability suppressing and delaying effect of rotation diminishes with an increase in $F$. The reason for this occurrence is the increase in the area of the first subharmonic unstable tongue resulting in the majority of excited modes being most unstable. We also predicted a wavelength ($\lambda ^{sh}$) of the subharmonic waves corresponding to the initially most amplified mode starting from the sharp interface, revealing that $\lambda ^{sh}$ increases with an increase in $f^2/\omega ^2(\leq 0.25)$.

Further, we have performed simulations at different forcing amplitudes and rotation rates to verify the findings of our theoretical model. The harmonic to subharmonic transition followed by the onset of the subharmonic instability, and its saturation, is shown in terms of the time evolution of the mixing zone size $L$. For $F = 0.75$ and $1$, the evolution of the mixing zone size manifests a significantly slow growth under the influence of rotation. This slow growth of $L$ represents the stable region between the harmonic to subharmonic transition and before the onset of the subharmonic instabilities. After this stable region, the subharmonic instability triggers, and $L$ rapidly grows, resulting in turbulent mixing. At forcing amplitudes of $2$ and $3$, the vertical vibrations are strong enough to diminish the instability delaying effect of rotation. Therefore, the mixing zone size for the cases with $F = 2$ and $3$ with rotation grows similar to the non-rotating cases. We also perform the Fourier transform of the rate of change of $L$ to capture the harmonic to subharmonic transition and demonstrate the delay in this transition in terms of the peak frequencies.

Finally, we have compared the simulations for $f/\omega =0.48$ and $f/\omega =0.59$, which correspond to $(\,f/\omega )^2 < 0.25$ and $(\,f/\omega )^2 \geq 0.25$, respectively, on the 3-D stability diagram, at different forcing amplitudes. We find that for $f/\omega =0.48$, the instability saturates, and the mixing zone size asymptotes. However, for $f/\omega =0.59$, the instability does not saturate, and the mixing zone size keeps growing. This confirms our theoretical observation that for $(\,f/\omega )^2 < 0.25$ the instability saturates, whereas it never saturates for $(\,f/\omega )^2 \geq 0.25$ because some $\theta$-modes are always excited in the first subharmonic tongue at all times.

The turbulent mixing, attributed to the onset of the subharmonic instability, is associated with exchanges among the TPE, the BPE and the APE. We aim to extend the present investigation to explore the influence of rotation on these energy exchanges through reversible buoyancy flux, irreversible diapycnal mixing and irreversible kinetic energy dissipation to quantify the irreversible mixing efficiency.