1. Introduction
Shear generated surface or internal waves in density stratified flows, particularly at high Reynolds numbers, have occupied a central space in hydrodynamic stability research since the pioneers (Helmholtz 1868; Rayleigh 1895; Thomson Reference Thomson1871; Taylor Reference Taylor1931; Goldstein Reference Goldstein1931). That important outstanding problems persist, attest to the complexity of the subject and highlight its relevance to wide ranging situations of geophysical, astrophysical or industrial interest. Adverse effects of climate change seemingly call for renewed attention into geophysical problems involving surface and internal waves; this is so because processes accompanying these often generate complicated, multiscale phenomena, both in the ocean and the atmosphere with long-term climate consequences. Examples include wave breaking (surface or internal gravity waves in the ocean and atmosphere) and associated mixing (Peltier & Caulfield Reference Peltier and Caulfield2003, section 6.2.1 in Fritts & Alexander Reference Fritts and Alexander2003 etc.) or ocean spray generation (Deike Reference Deike2022). Due to the multiscale nature, parameterisation of some of these processes in large-scale numerical models of wave or climate is necessary (Staquet & Sommeria Reference Staquet and Sommeria2002; Cavaleri et al. Reference Cavaleri2007; Janssen Reference Janssen2008; Perlin, Choi & Tian Reference Perlin, Choi and Tian2013) and in turn calls for accurate knowledge of the micro-physics. We refer the reader to figure 9 in Garrett & Munk (Reference Garrett and Munk1979) (figure
$5$
in Thorpe Reference Thorpe1975) for a caricature of some of these processes in the ocean.
The mechanism through which an overlying, moving gas or a liquid raises waves at the interface of an underlying, immiscible liquid layer (the common example being wind blowing over the deep ocean) has been of near constant interest since (Reference Jeffreys1925); see Pizzo, Deike & Ayet (Reference Pizzo, Deike and Ayet2021) for a recent summary and Mitsuyasu (Reference Mitsuyasu2002) for a historical overview of the wind-wave problem. The configuration also applies to industrial situations involving gas–liquid (Matas, Marty & Cartellier Reference Matas, Marty and Cartellier2011) or liquid–liquid combinations involving thick or thin films (Cohen & Hanratty Reference Cohen and Hanratty1965; Boomkamp & Miesen Reference Boomkamp and Miesen1996) subject to shear. Our interest here is restricted to the regime of geo/astro-physical interest involving two fluid layers (gas–liquid or liquid–liquid depending on the density ratio) with shear in the upper layer and a sharp density jump. Such a configuration is of oceanic interest at low density ratio (air–water) (Miles Reference Miles1957) and of astrophysical as well as geophysical interest at higher values (Pouliquen, Chomaz & Huerre Reference Pouliquen, Chomaz and Huerre1994; Alexakis, Young & Rosner Reference Alexakis, Young and Rosner2002). In the base state (also referred to as the ‘background state’) depicted in figure 1, a shear profile is imposed in the upper layer (the exponential velocity profile studied first in Morland & Saffman (Reference Morland and Saffman1993) with a sharp density interface at
$z=0$
, representative of a quasi-laminar wind profile), with quiescent lower fluid (stable density stratification). We study the linear stability of this configuration analytically and the time evolution of the most unstable mode computationally, within the inviscid framework, accounting for gravity and surface tension. Our analysis is restricted to unstable waves sufficiently long compared with the shear layer thickness (i.e. the background vorticity layer) as well as relevant capillary and viscous length scales (we return to viscous effects at the end of the study), but sufficiently short such that the deep layer approximation remains applicable. Within these simplifying approximations, we study transition in the nature of the fastest growing mode in the non-dimensional space comprising Froude number (
$\textit{Fr}$
) (ratio of inertia to gravity) and density ratio (
$\delta$
) (upper to lower fluid), at a large Bond (
$Bo\gg 1$
) number (ratio of gravity to surface tension). The precise definitions of these non-dimensional numbers are introduced in § 2.
Exponential velocity profile with quiescent fluid below, in the base (background) state with upper fluid velocity varying as
$U_u(z)=U_{\infty }[1-\exp (-z/\varDelta )]$
. A sharp interface separates two fluids of differing densities in a statically stable configuration:
$\rho (z)=\rho _l,\; z\lt 0$
and
$\rho (z)=\rho _u,\; z \gt 0$
with
$\rho _l \gt \rho _u$
.
The key result of our study is the demonstration, at high Bond number and moderately low Froude number, of a novel transition (for the fastest growing mode) from the Kelvin–Helmholtz (KH) (Thomson Reference Thomson1871; Helmholtz 1868) to the Holmboe (H) instability (Holmboe Reference Holmboe1962) and eventually to the Miles (Miles Reference Miles1957) critical layer instability, as the density ratio is reduced from
$0.9$
through
$0.5$
to the air–water value (
$0.001$
). While transition from the KH to the H instability (Smyth & Peltier Reference Smyth and Peltier1989; Lawrence, Browand & Redekopp Reference Lawrence, Browand and Redekopp1991, Reference Lawrence, Haigh and Zhu1998; Smyth & Winters Reference Smyth and Winters2003) for other background profiles or from the Miles to the KH instability (Morland & Saffman Reference Morland and Saffman1993; Alexakis et al. Reference Alexakis, Calder, Dursi, Rosner, Truran, Fryxell, Zingale, Timmes, Olson and Ricker2004a
) for the current profile have been reported, we demonstrate here the occurrence of all three aforementioned, inviscid, shear instabilities within the single base-state model of figure 1; to our knowledge, this is the first such demonstration.
Our analytical results receive support from comparison of the temporal evolution of the interface shape, obtained through solving the nonlinear, inviscid equations of motion. Simulations conducted using the open-source code Basilisk (Popinet & Collaborators Reference Popinet2013–Reference Popinet2025) show very good agreement with theoretical predictions in the linear regime. In the nonlinear regime, we note qualitative differences between the KH spirals (at
$\delta =0.9$
) and contrast these to sheared, cusp-like interface deformation emitting spume droplets for H-like (
$\delta =0.5$
) waves; these numerical observations are consistent with the numerical and experimental observations of Lawrence, Haigh & Zhu (Reference Lawrence, Haigh and Zhu1998) distinguishing the KH from the H instability for asymmetric velocity profiles. At
$\delta =0.1$
, simulations display features in the nonlinear regime, reminiscent of ripples on finite-amplitude, capillary-gravity Stokes waves. Their appearance typically accompanies an abrupt reduction of the parent wave amplitude. In the rest of this study, we refer to waves appearing for
$\delta \gg 0.1$
as internal waves consistent with experimental literature (Mercier, Vasseur & Dauxois Reference Mercier, Vasseur and Dauxois2011). The analysis of Yih (Reference Yih1960) showed that waves on a surface of density discontinuity are a limiting case of internal (gravity) waves appearing in continuously, stratified flows, thus supporting this nomenclature. In contrast, at very low density ratios (
$\ll 0.1$
), we describe surface waves.
1.1. Literature review
In the low density ratio regime (air–water,
$\delta =0.001$
), linear stability of the base state of figure 1 has been studied by several authors for a range of Froude and Bond numbers. For example, refer to Morland & Saffman (Reference Morland and Saffman1993), Leblanc (Reference Leblanc2007), Carpenter, Guha & Heifetz (Reference Carpenter, Guha and Heifetz2017), Bonfils et al. (Reference Bonfils, Mitra, Moon and Wettlaufer2022), Kadam, Patibandla & Roy (Reference Kadam, Patibandla and Roy2023) for linear stability studies with distinct perspectives and to Leblanc (Reference Leblanc2007), Leblanc (Reference Leblanc2007) for weakly nonlinear investigations; these studies are relevant to wind-driven, water waves with
$\delta =0.001$
. Apparently, situations involving immiscible fluids subject to shear at density ratios significantly higher than the air–water value also occur commonly. Fresh and salt water, for example, have a density ratio
${\gt } 0.95$
; a sharp interface is maintained on experimental time scales between the two; see Mercier et al. (Reference Mercier, Vasseur and Dauxois2011). Pouliquen et al. (Reference Pouliquen, Chomaz and Huerre1994) studied the KH and H instability experimentally using fluids of a high density ratio
$\sim 0.92$
(water and equal mixture of silicone oil with 1–2–3–4-tetrahydronaphthalene; see their table 1), comparing their experimental results against inviscid stability theory. That an inviscid theory suffices to describe their experimental results involving somewhat viscous fluids is justified from the estimate of the viscous diffusion time scale on p.
$282$
of Pouliquen, Chomaz & Huerre (Reference Pouliquen, Chomaz and Huerre1994) along with the Reynolds number (see p.
$293$
). Beyond geophysical scales, Alexakis et al. (Reference Alexakis, Young and Rosner2002) report that a density ratio of
$0.1$
in shear flows occurs in astrophysical situations involving accretion from a normal star to a white dwarf star; see the discussion on pp.
$1$
–
$2$
of Alexakis et al. (2002).
To our knowledge, the only analytical study of the stability of the exponential profile of figure 1 at a larger density ratio (compared with air–water) has been carried out by Alexakis et al. (Reference Alexakis, Young and Rosner2002, Reference Alexakis, Calder, Dursi, Rosner, Truran, Fryxell, Zingale, Timmes, Olson and Ricker2004a
). These authors studied the moderately low Froude number and large wavenumber limit theoretically, obtaining an analytical expression for the growth rate (in the linear regime) at an arbitrary density ratio; no analysis of the nature of the instability or any transition therein with varying density ratio was reported, as is done here. Alexakis, Young & Rosner (Reference Alexakis, Young and Rosner2004b
) studied the weakly nonlinear regime for a small Froude number (the weak wind limit) and small density ratio, obtaining amplitude equations that were then solved numerically. These solutions validated the prediction (Haberman Reference Haberman1972) that, for a growing mode (Miles Reference Miles1957), the phase change across the critical layer expected from linear theory saturates to zero as the wave amplitude increases with time. The nonlinear simulations by Alexakis et al. (Reference Alexakis, Calder, Dursi, Rosner, Truran, Fryxell, Zingale, Timmes, Olson and Ricker2004a
) were however in the moderate-Froude-number
$ \textit{Fr}=2.3{-}9.5$
and low-Mach-number regime, at a fixed density ratio of
$0.1$
. Qualitative agreement with linear theory was reported for the growth of gravitational potential energy of unstable Fourier mode(s), although significant scatter is apparent in their comparison; see Alexakis et al. (Reference Alexakis, Calder, Dursi, Rosner, Truran, Fryxell, Zingale, Timmes, Olson and Ricker2004a
, Figure 4). At these Froude numbers, the authors reported regimes separating the Miles from the KH instability, although their KH model is similar to that in figure 2(b) with a jump at the interface, and not the more general continuous and piecewise linear (PL) model of figure 2(a). We will demonstrate later on that at
$\delta =0.1$
(as is the case for Alexakis et al. Reference Alexakis, Calder, Dursi, Rosner, Truran, Fryxell, Zingale, Timmes, Olson and Ricker2004a
) and low Froude number, the instability in the exponential model of figure 1 begins to lose its signature at the critical location, indicative of a potential change in the instability mechanism.
(a) The PL profile with a shear layer of thickness
$\varDelta$
. The upper fluid background velocity is given by
$U_u(z)=U_\infty ({z}/{\varDelta }),\; 0\leqslant z \leqslant \varDelta$
and
$U_u(z)=U_{\infty },\; z \geqslant \varDelta$
. (b) The classic KH discontinuous profile, obtained when
$\varDelta \rightarrow 0$
in the PL profile in panel (a). In both panels, the lower fluid is quiescent in the base state, i.e.
$U_l(z)=0,\; z\lt 0$
. The statically stable, base-state density profile is
$\rho (z)=\rho _l,\; z\lt 0$
and
$\rho (z)=\rho _u,\; z \gt 0$
with
$\rho _l \gt \rho _u$
.
We note the important study by Morland & Saffman (Reference Morland and Saffman1993), which explored transitions from the Miles instability to the discontinuous KH instability in the exponential model, at a low (air–water) density ratio and for a large
$\textit{Fr}$
and infinite Bond number. Importantly, they proved the continuation of the Howard semi-circle theorem to shear flow with a sharp density interface; this had been done by Yih (Reference Yih1972) albeit for surface waves (
$\delta =0$
) and also subsequently by Alexakis et al. (Reference Alexakis, Young and Rosner2002) for arbitrary
$\delta$
. In § 3.1 we will see that even at
$ \textit{Fr} \sim 3000 \gg 1$
, the growth rate of the fastest growing mode in the exponential model agrees much better with that of the PL model (figure 2
a) rather than the discontinuous KH model of figure 2(b) considered by Morland, Saffman & Yuen (Reference Morland, Saffman and Yuen1991) and Alexakis et al. (Reference Alexakis, Young and Rosner2002). With this brief survey above, it becomes clear that the nature of the linear instability for the exponential model of figure 1 as a function of the density ratio and transitions, if any, for the fastest growing mode are not known. This is demonstrated in the present study by employing theory and numerical simulations.
1.2. Inviscid, shear instabilities and the exponential model
For the exponential profile in figure 1, the background shear layer thickness is approximately
$\approx \varDelta$
. With a sharp density interface placed within this shear layer, three inviscid instabilities can potentially manifest, depending on the parameter regime. Two of these (the KH and H), are approximately understood using PL velocity profiles (their continuations to smooth profiles also being well known) while for the third instability (Miles Reference Miles1957), curvature of the background velocity profile is crucial. For completeness, we recap some salient features.
Recall that the KH instability occurs through an inflexional instability mechanism (see discussion on pp.
$277{-}278$
of Pouliquen et al. (Reference Pouliquen, Chomaz and Huerre1994) applicable to a homogeneous fluid under shear (i.e. no stratification). The textbook KH profile is usually the base state of figure 2(b) (a vortex sheet), representing the limiting case of an inflectional velocity profile (Michalke Reference Michalke1964) such as a
$\tanh$
profile; the vortex sheet of figure 2(b) is unstable to all wavenumbers when the density ratio (
$\delta \equiv \rho _u/\rho _l$
) is unity. Density stratification plays a stabilising role for the sheet – only sufficiently large wavenumbers are unstable when
$0 \lt \delta \lt 1$
. The KH instability can also manifest in the PL model of figure 2(a) for perturbations whose wavelength
$\lambda$
is large compared with the shear layer thickness (
$\lambda \gg \varDelta$
). By replacing the smooth, background vorticity gradient (
${\rm d}^2U_u/{\rm d}z^2$
) of figure 1 with the apparently simpler PL profile of figure 2(a) with a jump in background vorticity at
$z=\varDelta$
, the KH instability can also be interpreted as resonance between waves on the two vorticity interfaces at
$z=0$
and
$z=\varDelta$
; see discussion around figure
$1$
in Baines & Mitsudera (Reference Baines and Mitsudera1994) and around equation
$33$
in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011) for smooth profiles.
In contrast to the KH instability, the H instability (Holmboe Reference Holmboe1962) occurs due to the destabilising influence of stable stratification under shear, manifesting in situations where the thickness of the shear layer is much larger than that of the density layer. Employing the same PL profile of figure 2(a), one can interpret the H instability as resonance between a wave on the density and that on the vorticity interface(s); see figure 16(d) in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011). The difference between the H and KH instability is seen most readily when the centre of the shear layer coincides with the density interface (symmetric profiles) leading to stationary (KH) and travelling (H) modes. The profiles of figures 1 and 2(a) however are asymmetric, the latter represents a limiting case of the profile in Lawrence, Browand & Redekopp (Reference Lawrence, Browand and Redekopp1991) and we will return to this asymmetry later in our study.
The third instability admissible for the profile of figure 1 is the critical layer instability of Miles (Reference Miles1957). Although the Rayleigh equation is strictly not singular at any real value of the vertical coordinate
$z$
for a growing mode with complex
$c$
, for our present study, we define the critical location as the (real) location
$z_c$
where
$U_u(z=z_c)=c_r(k)$
(
$c_r$
denotes the real part of the phase speed
$c$
), this being done for unstable modes. Unlike the H and KH instabilities, non-zero curvature of the background profile at the critical location is crucial to the instability (see below). This instability, more broadly, belongs to the ‘class B’ family of instabilities of Benjamin (Reference Benjamin1960, Reference Benjamin1963). As explained in Benjamin (Reference Benjamin1960) last paragraph on p.
$514$
), the instability ‘takes the form of waves progressing in the flow direction at very nearly the same speed as free surface waves in the flexible medium, the waves being amplified by the action of the flow which supplies sufficient energy to counterbalance internal dissipation’. At low density ratio, Miles (Reference Miles1957) demonstrated that, for this instability, the growth rate of a Fourier mode is proportional to
$(-{\rm d}^2U_u/{\rm d}z^2)_{z_c}$
,
$U_u(z)$
being the background shear flow in the upper fluid (subscript ‘u’ for upper). In regimes where this critical location dominates the instability mechanism, we will verify that for the fastest growing Fourier mode, a near-discontinuous jump in the (inviscid) Reynolds stress appears at the critical location and upto the interface (also see figure
$2$
in Boomkamp & Miesen (Reference Boomkamp and Miesen1996). (We recall that the critical layer plays a crucial role in the viscous instability of non-inflectional, velocity profiles in unstratified, parallel shear flows, e.g. the plane Poiseuille flow (Lighthill & Rosenhead (Reference Lighthill and Rosenhead1963) figure IX.8.) or TS instability of the Blasius boundary layer (Ruban, Gajjar & Walton (Reference Ruban, Gajjar and Walton2023)
$2$
). The associated work done by this shear stress leads to extraction of energy from the background shear primarily in the region between the critical location and the interface. This will be contrasted against situations where the energy extraction occurs in the entire shear layer with no distinct signature of the Reynolds stress at the critical location; this is despite the mathematical existence of a critical location for every growing Fourier mode, as ensured by the semi-circle theorem.
While we have discussed the three instabilities above in an apparent mutually exclusive manner, critical locations appear within the inviscid framework (Taylor–Goldstein equation) for the KH as well as H instabilities with smooth velocity (and density) profiles. Refer to figures
$5$
and
$7$
in Smyth & Peltier (Reference Smyth and Peltier1989), where phase jumps and peaks are seen at the critical location for the eigenfunction and the Reynolds stress for unstable KH and H modes, respectively. We bear this aspect in mind when distinguishing between critical layer Miles instabilities at low density ratios and H and KH instabilities at higher values. In addition, comparison of instabilities obtained from curved background profiles (figure 1) with their straight line counterparts (figure 2), as is done throughout our study, despite being common, requires caution; see, e.g. Morland et al. (Reference Morland, Saffman and Yuen1991) for the distinction in the transition mechanism between instability in the PL profile and that of the exponential profile for a sheared lower fluid (
$\delta =0$
). For background profiles with a smooth second derivative, the logarithmic singularity in Tollmien’s inviscid solution (applicable to a neutral mode, provided such a mode exists) gets regularised and spread over a finite thickness for a growing mode. Note that the neutral mode is singular at a (real)
$z$
where the second derivative of the background profile is non-zero and the modal phase speed matches the background speed. Interestingly, it is possible to describe this critical layer region via kinks present in PL profiles; a large number of kinks being necessary for an accurate representation; see Carpenter et al. (Reference Carpenter, Guha and Heifetz2017). Closely related to our study, we also note the interesting study by Eaves & Balmforth (Reference Eaves and Balmforth2019), where several diagnostics, both analytical and numerical, in order to distinguish the Taylor-Caulfield instability from the H and KH instability are discussed.
This study is organised as follows. In § 2 we discuss the temporal, linear stability of the exponential profile in figure 1 obtaining the dispersion relation. The same is also done for the PL profile in figure 2(a). Section 3 presents a detailed exploration of the parameter space for high and low Froude numbers with increasing density ratio, for high Bond numbers. Section 4 presents numerical simulations validating linearised predictions but also presenting several novel results in the nonlinear regime. In § 5 we revisit the validity of the background state of figure 1, particularly for higher density ratios, close to unity. We conclude with a summary in § 6.
2. Temporal, linear stability analysis
The base state (figure 1) comprises a stably stratified (
$\rho _l \gt \rho _u$
) configuration with a sharp density jump at
$z=0$
and shear (only) in the upper fluid:
\begin{align} & U(z) = \left \{ \begin{array}{ll} U_u(z)=U_{\infty }\left (1-\exp \left (-\dfrac {z}{\varDelta }\right )\right ), & \qquad z \gt 0, \\[6pt] U_l(z)=0, & \qquad z \lt 0 ,\end{array}\right . \\[-12pt] \nonumber \end{align}
\begin{align} & \rho (z) = \left \{ \begin{array}{ll} \rho _u, & \qquad z \gt 0, \\[6pt] \rho _l, & \qquad z \lt 0. \end{array}\right .\\[10pt]\nonumber \end{align}
In the base state the interface is flat at
$z=0$
. The perturbed state is represented by the streamfunctions
$\psi _u(x,z,t)$
and
$\psi _l(x,z,t)$
(for the upper and lower fluid, respectively), while the deviation of the interface from its flat state is measured by
$\eta (x,t)$
(see figure 1). The vertical variation of the perturbation streamfunction (
$\varPsi _{u/l}$
) is governed by the Rayleigh equations below (subscripts
$u$
or
$l$
refer to upper and lower layers, respectively):
\begin{align} & \big \lbrace U_{u/l}(z)-c(k)\big \rbrace \big(D^2 - k^2 \big)\varPsi _{u/l}(z) - D^2U_{u/l}(z)\varPsi _{u/l}(z) = 0. \end{align}
Here
$D\equiv {\text{d}}/{\text{d}z}$
and
$\varPsi _{u/l}(z)$
is related to the streamfunction
$\psi _{u/l}(x,z,t)$
via a travelling wave, Fourier mode representation, viz. (we omit the
$u$
or
$l$
subscripts below)
\begin{eqnarray} &&\psi (x,z,t) = \dfrac {1}{2}\varPsi (z)\exp \left [i k(x-ct)\right ] + \text{c.c.},\quad u \equiv \dfrac {\partial \psi }{\partial z},\quad w\equiv -\dfrac {\partial \psi }{\partial x}. \end{eqnarray}
Here
$i\equiv \sqrt {-1}$
,
$u$
and
$w$
are the horizontal and vertical, perturbed velocity components, respectively, while c.c. denotes the complex conjugate. The perturbed interface (vertical displacement) is similarly represented using complex notation as
\begin{eqnarray} \eta (x,t)=\dfrac {1}{2}\eta _0\exp \left [i k(x-ct)\right ] + \text{c.c.}, \end{eqnarray}
where
$\eta _0$
is a complex constant. Equation (2.2) represents eigenvalue problems for determining complex functions
$\varPsi _{u}(z),\varPsi _{l}(z)$
and
$c(k)$
subject to boundary conditions (see below) for the chosen base-state profile in (2.1) and for real
$k$
. The boundary conditions are the (linearised) kinematic boundary condition
\begin{align} c\eta _0 = \varPsi _u(0)= \varPsi _l(0), \end{align}
and the jump in normal stress (pressure in the inviscid framework here) due to surface tension (
$T$
) is
\begin{eqnarray} P_u(0) - P_l(0) = \left (\rho _u-\rho _l\right )g\eta _0- Tk^2\eta _0, \end{eqnarray}
where
$p(x,z,t)=({1}/{2})P(z)\exp [i (kx-ct) ]+\text{c.c.}$
is the perturbation pressure (with appropriate subscripts). The expressions relating the perturbation streamfunction
$\varPsi (z)$
to the vertically varying part of the perturbation pressure (
$P(z)$
) are
\begin{align} P_u(0)=-\rho _u\left [(U_u(z)-c)D\varPsi _{u} - \left (DU_u\right )\varPsi _{u}\right ]_{z=0} ,\quad P_l(0)=c\rho _l\left (D\varPsi _{l}\right )_{z=0}. \end{align}
In addition to the boundary conditions above, it is also necessary to ensure that the solution to (2.2) does not diverge as
$z\rightarrow \pm \infty$
and this determines two constants of integration in expressions for
$\varPsi _{l}(z)$
and
$\varPsi _{u}(z)$
. As shown by J. Miles in Morland & Saffman (Reference Morland and Saffman1993), see Russell (Reference Russell1994, § 3.3) for details, the Rayleigh equation for the upper fluid may be converted into the Gauss hypergeometric equation whose analytical properties are well known. The Rayleigh equation for the lower fluid has constant coefficients and is thus straightforward to solve. We thus obtain expressions for the perturbation streamfunction in both fluids as
\begin{align} \varPsi _u(z) = M\,\, \exp \left [-k\left (z-z_c\right )\right ]\, _2F_1\left (\alpha ,\beta ,\gamma ;\dfrac {\exp \left (-\frac{z}{\varDelta }\right )}{1 - \dfrac{c}{U_{\infty }}}\right ), \,\, \varPsi _{l}(z) = N\exp \left (kz\right ), \end{align}
with
$\alpha \equiv k\Delta + \sqrt {1 + k^2\varDelta ^2}, \beta \equiv k\Delta - \sqrt {1 + k^2\varDelta ^2}$
and
$\gamma \equiv 1 + 2k\Delta$
, where
$M$
and
$N$
are (complex) constants of integration. In (2.8a
),
$\;_2F_1(\boldsymbol{\cdot },\boldsymbol{\cdot },\boldsymbol{\cdot };\boldsymbol{\cdot })$
is the Gauss hypergeometric function defined by the series
\begin{eqnarray} _2F_1(\alpha ,\beta ,\gamma ;t) = 1 + \dfrac {\alpha \beta }{\gamma }\dfrac {t}{1!} + \dfrac {\alpha \left (\alpha +1\right )\beta \left (\beta +1\right )}{\gamma \left (\gamma +1\right )}\dfrac {t^2}{2!} + \cdots ,\quad \quad |t| \lt 1 .\end{eqnarray}
Note that with
$t \equiv {\exp (-z/\varDelta )}/{(1-(c/U_{\infty }))}$
, as
$z$
varies in the range
$\infty \lt z \leqslant z_c$
, we have the variation
$0\lt t \leqslant 1$
assuming that
$z_c$
is real (for a neutral mode) and determined using
$U_u(z_c)=c(k)$
. For the range
$0 \leqslant z \leqslant z_c\Rightarrow ({U_{\infty }}/({U_{\infty }-c_r}) ) \geqslant t \geqslant 1$
assuming that
$c \lt U_{\infty }$
. This necessitates the continuation of
$\;_2F_1(\boldsymbol{\cdot },\boldsymbol{\cdot },\boldsymbol{\cdot };t)$
for
$|t|\gt 1$
for which, for example, the Euler’s integral representation of
$_2F_1$
may be used (see further discussion in Andrews et al. Reference Andrews, Askey and Roy1999, p.
$65$
):
\begin{align} & _2F_1(\alpha ,\beta ,\gamma ;t) \equiv \dfrac {\varGamma (\gamma )}{\varGamma (\beta ) \varGamma (\gamma -\beta )}\int _{0}^{1} \dfrac {r^{\beta -1}(1-r)^{\gamma -\beta -1}}{\left ( 1- rt\right )^{\alpha }}\,\text{d}r. \end{align}
Using the kinematic boundary condition (2.5a,b ) in (2.8), we obtain
\begin{align} \varPsi _{u}(z) = c\eta _0 \exp \left (-kz\right )\kern-1pt\left [\kern-1pt\dfrac {_2F_1\left (\alpha ,\beta ,\gamma ; \dfrac{\exp \left (- \dfrac{z}{\varDelta }\right )}{1 - \dfrac{c}{U_{\infty }}}\right )}{_2F_1\left (\alpha ,\beta ,\gamma ;\dfrac{1}{1-\dfrac{c}{U_{\infty }}}\right )}\kern-1pt\right ]\!, \, \varPsi _l(z) = c\eta _0\exp \left (kz\right )\!.\end{align}
We also note that
$t=1$
is a branch point of
$_2F_1(\alpha ,\beta ,\gamma ;t)$
; this generates a (near) jump in phase of the streamfunction
$\varPsi _u(z)$
in (2.11a
) for a Miles growing mode (i.e. a
$k$
for which
${\rm Im}(c)= c_i\gt 0$
), as the critical location
$z_c$
is crossed (
$t=1^{+}$
). This in turn also generates the phase shift in perturbation pressure at
$z=0^{+}$
compared with the interface
$\eta$
. This phase shift produces the Jeffreys (Reference Jeffreys1925) mechanism of generating wind waves, extensively referred to in Miles (Reference Miles1957).
Using expressions for
$\varPsi _{u}(z)$
and
$\varPsi _{l}(z)$
from (2.11) in (2.7) and the normal stress condition in (2.6), we obtain the following dispersion relation. In non-dimensional form, this may be expressed as
\begin{align} & \left [\delta \left \lbrace \kappa - \left (\dfrac {1}{1+2\kappa }\right )\left (\dfrac {\textit{Fr}}{\textit{Fr}-\tilde {c}}\right )\dfrac {_2F_1\left (\alpha +1,\beta +1,\gamma +1;\dfrac{\textit{Fr}}{\textit{Fr}-\tilde {c}}\right )}{_2F_1\left (\alpha ,\beta ,\gamma ;\dfrac{\textit{Fr}}{\textit{Fr}-\tilde {c}}\right )}\right \rbrace + \kappa \right ]\tilde {c}^2 \\\nonumber & \qquad \qquad -\delta \textit{Fr}\tilde {c} - \left (1-\delta \right ) - \dfrac {1}{Bo}\kappa ^2 = 0,\end{align}
where
\begin{eqnarray} \tilde {c} \equiv \dfrac {c}{\sqrt {g\Delta }},\; \kappa \equiv k\Delta ,\; \textit{Fr} \equiv \dfrac {U_{\infty }}{\sqrt {g\Delta }},\; \delta \equiv \dfrac {\rho _u}{\rho _l},\;Bo \equiv \dfrac {\rho _lg\varDelta ^2}{T}. \end{eqnarray}
Here
$\textit{Fr}$
,
$Bo$
and
$\delta$
represent the Froude number, Bond number and the density ratio, respectively. For a given
$\kappa$
(wavenumber), the complex roots (
$\tilde {c}$
) to (2.12) are obtained numerically using the findroot function from the mpmath library (Johansson Reference Johansson2024) in Python, which provides arbitrary-precision root finding based on Muller’s method (Müller Reference Müller1956). In the supplementary material (§ 1) we provide benchmarking with the results of Morland & Saffman (Reference Morland and Saffman1993) and Young & Wolfe (Reference Young and Wolfe2014) for validation of the numerical procedure of solving (2.12). In the following sections, we extensively use the symbols
$c_g,c_{gc}$
and
$c_{_{\textit{KH}}}$
whose formulae are provided below in the deep limit.
\begin{align} &\text{Unsheared surface-gravity wave speed:} && c_g^2 \equiv {\left ( \dfrac {1-\delta }{1+\delta } \right ) \dfrac {g}{k}}. \\[-12pt] \nonumber \end{align}
\begin{align} &\text{Unsheared capillary-gravity wave speed :} && c_{gc}^{2} \equiv \left (\dfrac {1-\delta }{1+\delta }\right )\dfrac {g}{k} + \dfrac {Tk}{\rho _l\left (1+\delta \right )} .\\[-12pt] \nonumber \end{align}
\begin{align} & \text{The KH wave speed:} && c_{_{\textit{KH}}} \equiv \left (\dfrac {\delta }{1+\delta }\right )U_\infty. \\[10pt] \nonumber \end{align}
3. Exploring the parameter space
In the following analysis, we explore the nature of the instability in the non-dimensional space of Froude number (
$\textit{Fr}$
) and density ratio (
$\delta$
) for given value(s) of the Bond number (
$Bo$
). Wherever relevant, we compare the results of the PL profile of figure 2(a) with that of the exponential profile of figure 1, both profiles chosen to have the same vertical scale
$\varDelta$
. We investigate four representative cases by varying
$\delta$
and
$\textit{Fr}$
at fixed
$Bo$
: (i) high
$ \textit{Fr} = 3194.38$
, low
$\delta = 0.001$
; (ii) high
$ \textit{Fr} = 3194.38$
, high
$\delta = 0.5$
; (iii) moderately low
$ \textit{Fr} = 4.79$
, low
$\delta = 0.001$
; and (iv) moderately low
$ \textit{Fr} = 4.79$
, high
$\delta = 0.5$
. For the benefit of the reader, table 1 summarises key observations, as will be inferred from following analysis. We commence our parametric exploration with the high-Froude-number case with zero surface tension, i.e. the infinite Bond case.
Summary of observations based on fastest growing modes.
3.1. Infinite bond (
$Bo\rightarrow \infty$
) and high Froude number (
$ \textit{Fr} \gg 1$
)
3.1.1. Low density ratio (
$\delta \ll 1$
):
To our knowledge the first study that investigated transition in the nature of the instability for the exponential profile of (2.1) with a sharp density jump at
$z=0$
was by Morland & Saffman (Reference Morland and Saffman1993). These authors solved the Rayleigh equation numerically for the exponential profile (quiescent water in base state) for air–water (
$\delta =0.001$
) taking into account gravity but setting surface tension to zero (
$Bo\rightarrow \infty$
). Among other conclusions, Morland & Saffman (Reference Morland and Saffman1993) noted a change in the nature of the instability at
$ \textit{Fr} \gg 1$
(for constant
$\delta$
), transitioning from ‘Miles unstable’ (critical layer induced instability) gravity waves with approximate phase speed
$c_g^2$
(in deep water, (2.14a
)) to KH unstable waves (i.e. waves appearing on the base state of figure 2
b) with a phase speed given by
$c_{_{\textit{KH}}}$
(2.14c
). Such KH-like behaviour is anticipated for unstable modes of wavelength
$\lambda$
satisfying
$\varDelta \ll \lambda \lt ({2\pi U_{\infty }^2\delta }/{g (1-\delta ^2 ))}$
(see discussion below (3.2) in Morland & Saffman (Reference Morland and Saffman1993)). It is also intuitively clear from examining the exponential profile in (2.1) that taking the limit
$\varDelta \rightarrow 0$
(
$ \textit{Fr}\rightarrow \infty$
) for fixed
$U_{\infty },\;g$
and density ratio
$\delta$
, deforms the exponential profile in figure 1 towards a uniform velocity profile (
$U_{\infty }$
) in the upper fluid; in other words, the classic KH base-state profile of figure 2(b). (For the exponential profile, the
$ \textit{Fr}\rightarrow \infty$
limit for fixed
$\delta \ll 1$
as discussed here is qualitatively different from
$\delta \rightarrow\kern-1pt 0,\kern-1pt\textit{Fr}\rightarrow\kern-1pt\infty$
, with the limits taken in that order; see discussion below equation (
$7.17$
) in Young & Wolfe (Reference Young and Wolfe2014).) In dimensional variables, the dispersion relation for this KH base state is (Drazin & Reid (Reference Drazin and Reid2004), Equation (4.20))
\begin{eqnarray} \dfrac {c(k)}{U_{\infty }} = \dfrac {\delta }{1+\delta } \pm \left [\left (\dfrac {1-\delta }{1+ \delta }\right )\dfrac {g}{kU_{\infty }^2} - \dfrac {\delta }{(1+\delta )^2}\right ]^{1/2}. \end{eqnarray}
The relation (3.1) predicts that unstable KH modes travel with the speed
$c_{_{\textit{KH}}}$
, independent of wavenumber. At any
$U_{\infty }$
, sufficiently long waves are stable (due to gravity), while shorter waves are destabilised by the shear and grow exponentially in time.
We also compare predictions from the dispersion relation (2.12) for the exponential profile in the high-Froude-number limit with that of the PL velocity profile indicated in figure 2(a). In terms of the non-dimensional variables and parameters introduced earlier, the dispersion relation (cubic polynomial in
$\tilde {c}$
) for the PL profile of figure 2(a) can be shown to be
\begin{align} \begin{split} & \big ( \textit{Fr}(1 - 2 \kappa ) + 2\tilde {c}\kappa \big ) \boldsymbol{\cdot }\left ( \dfrac {1 - \delta }{1 + \delta }(-2 + \tilde {c}\textit{Fr}) - \tilde {c} \textit{Fr} + 2\tilde {c}^2 \kappa \right ) \\ & \qquad + {\rm e}^{-2\kappa } \boldsymbol{\cdot }\textit{Fr} \boldsymbol{\cdot }\left \{ \tilde {c} \textit{Fr} \left (1 - \dfrac {1 - \delta }{1 + \delta }\right ) + 2\left (\dfrac {1 - \delta }{1 + \delta }\right ) \big(1 - \tilde {c}^2 \kappa \big) \right \}=0, \quad \tilde {c}\neq 0. \end{split} \end{align}
Plotted in figures 3(a) and 3(b) are the growth rate and phase speed of the exponential profile (labelled as ‘Exp’) as obtained from (2.12), for the KH profile from (3.1) and for the PL profile from (3.2), at a high Froude number, viz.
$ \textit{Fr}=3194.39$
. A few qualitative features are discernible readily. On account of having an additional length scale
$\varDelta$
compared with the classic KH model in figure 2(b), the stability characteristics of the PL profile (‘
$\boldsymbol{\cdot }\times \boldsymbol{\cdot}$
’ symbols) differ qualitatively at large wavenumbers from the KH model (dashed pink line), as seen in figure 3. Wavenumbers with
$\kappa \geqslant O_m(0.1)$
(
$\mathcal{O}_m$
denotes order of magnitude) are stable in the PL model – this is the reason for the abrupt termination of the ‘
$\boldsymbol{\cdot }\times \boldsymbol{\cdot}$
’ symbols (PL model) around
$\kappa \sim 10^{-1}$
in figure 3(a), above this wavenumber the growth rate goes to zero rapidly. For the exponential profile, also recall that a lower cutoff in wavenumber exists at infinite depth (Young & Wolfe (Reference Young and Wolfe2014), Kadam et al. (Reference Kadam, Patibandla and Roy2023)); wavenumbers less than this cutoff are stable (Young & Wolfe Reference Young and Wolfe2014) as seen from the abrupt decay in growth rate in figure 3(a) at low
$\kappa$
.
The infinite Bond (
$Bo\rightarrow \infty$
), high Froude (
$ \textit{Fr}=3194.38$
) number and low (air–water) density ratio (
$\delta =0.001$
) case. (a) Growth rate (
$\kappa \tilde {c}_i$
) and (b) phase speed (
$c_r/c_{_{\textit{KH}}}$
) as a function of the non-dimensional wavenumber
$\kappa$
for the KH dispersion relation from (3.1) (pink dashed curve) for the exponential profile (solid blue curve) from (2.12) and for the PL profile from the dispersion relation in (3.2) (‘
$\boldsymbol{\cdot }\times \boldsymbol{\cdot}$
’ symbols) for the base states of figures 1 and 2. We see that the growth rate and phase speed are very similar for
$\kappa \rightarrow 0$
. The KH dispersion relation of (3.1) does not contain the Froude number as a parameter and represents the
$ \textit{Fr}\rightarrow \infty$
limit. The red and green dots in both panels depict a small wavenumber and the fastest growing wavenumber, respectively.
Figures 3(a) and 3(b) also show that as
$\kappa \rightarrow 0$
(long wave indicated by the red dot), the growth rate and phase speed of the exponential, as well as the PL modes, asymptote to the KH values provided by expression (3.1). The fastest growing mode for the exponential profile is indicated by the green dot in both figures. The red dot corresponds to a Fourier mode of wavelength
$\approx 230$
times the wavelength of the green dot. It is clear that both are long wavelength modes compared with the shear layer thickness (
$\lambda /\varDelta \approx 3.1\times 10^4$
for the red dot and
$\approx 1.4\times 10^2$
for the green dot). A qualitative difference is apparent between these two, long wavelength modes: the fastest growing mode (green dot) of the exponential profile grows somewhat slower than the KH growth estimate of (3.1), travelling at a speed that is nearly twenty times the KH prediction; cf. figure 3(b). Interestingly, its growth rate and phase speed both agree quite well with that of the PL model (‘
$\boldsymbol{\cdot }\times \boldsymbol{\cdot}$
’ symbol).
Having established that the fastest growing mode of the exponential profile in the high-Froude-number (and infinite Bond number and low density ratio) limit has a growth rate and phase speed agreeing with the corresponding mode of the PL profile, it becomes necessary to inquire about the nature of this mode in the PL profile. For this, we recall a few essential stability properties of the PL profile. This base state of figure 2(a), by virtue of having a density interface (zero thickness) much smaller in thickness than the shear layer, admits the H instability (interaction between a vorticity and a density interface) as well as the KH instability (interaction between two vorticity interfaces; see Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011) and references therein). Under the Boussinesq approximation, depending on the choice of parameters on the Richardson number
$J\equiv ({1-\delta }){g\Delta }/({U_{\infty }^2}({1+\delta }) )$
versus wavenumber (
$\kappa$
) space, the instability can be dominantly one of the two types or even of a mixed nature. As noted by Lawrence et al. (Reference Lawrence, Browand and Redekopp1991), the PL profile that we employ is an asymmetric profile, implying that the centre of its shear layer (
$z = \varDelta /2 \gt 0$
; see figure 2
a) does not coincide with the location of the density interface (
$z=0$
). Indeed, the velocity profile of figure 2(a) may be recovered from that of Lawrence et al. (Reference Lawrence, Browand and Redekopp1991) by setting
$u_2=0$
and their asymmetry parameter
$\epsilon$
to unity (refer to figure
$1$
of Lawrence et al. Reference Lawrence, Browand and Redekopp1991). Unlike symmetric profiles, for asymmetric profiles, the transition between the KH and H instability is not readily apparent on the
$J-\kappa$
plane, as both modes have non-zero phase speeds; see Carpenter et al. (Reference Carpenter, Balmforth and Lawrence2010, §
$B$
, p. 9) and the related study by Barros & Choi (Reference Barros and Choi2011) for non-Boussinesq effects. Using wave resonance ideas reviewed in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011), one can use the phase speed versus wavenumber plane to conclude if the dominant contribution comes from the KH (vorticity–vorticity resonance) or the H instability (vorticity–density resonance) for the fastest growing mode of figure 3(a). We defer this exercise to § 3.2.1 (the low-Froude-number limit), where we show that at low Froude number and low as well as moderately high density ratio, growing modes in the PL model arise due to resonance between a surface-gravity mode and a vorticity wave (i.e. the H instability dominates over the KH instability).
Figures 4(a) and 4(b) present the vertical variation of the (inviscid) Reynolds stress
$-\overline {\rho uw}(z)$
and (twice) the perturbation kinetic energy
$\overline {u^2+w^2}(z)$
for both wavenumbers (green and red dots) of figure 3. In these (and subsequent figures) the Reynolds stress and perturbation kinetic energy have been normalised by their respective values at
$z=0$
. The overbar indicates averaging along the horizontal coordinate over one wavelength, i.e.
\begin{eqnarray} \overline {uw}(z, t) \equiv \dfrac {1}{\lambda }\int _{x}^{x+\lambda }u(x,z,t)\;w(x,z,t) \,\text{d}x. \end{eqnarray}
Recall that the base-state vorticity (
$DU_u(z)$
from (2.1)) is non-zero and positive in the approximate range
$0 \lt z \lt \varDelta$
in the upper fluid. This corresponds to
$0\lt z/z_c \lt \varDelta /z_c = 51.4,(\text{wavelength}\; \lambda /z_c\approx 7.2\times 10^3)$
for the green curves and
$0\lt z/z_c \lt \varDelta /z_c =832.2, (\lambda /z_c\approx 2\times 10^7)$
for the red curves, in both panels of figure 4. It is clear that the perturbation kinetic energy
$ (\overline {u^2+w^2} )(z)$
as well as the Reynolds stress (
$ (-\overline {\rho uw} )(z)$
) becomes non-zero close to the onset of the shear layer for both modes. While the fastest growing mode (green curves) reaches its plateau very close to its maximum value at the critical layer (
$z/z_c = 1$
) in both panels, the KH mode (red) plateaus significantly before the critical layer is reached, i.e.
$z/z_c \sim 100$
. One may derive an equation for the rate of growth of perturbation kinetic energy, viz.
$E_k(z) \equiv ({1}/{2})\rho (\overline {u^2+w^2} )(z)$
, which is (see Smyth & Peltier (Reference Smyth and Peltier1989), Equation (
$4$
))
\begin{eqnarray} \dfrac {\partial E_k(z)}{\partial t} = -(DU(z))\overline {\rho uw}(z) - D\overline {wp}(z). \end{eqnarray}
Vertical (
$z$
) profiles of (a) the Reynolds stress and (b) the perturbation kinetic energy in the upper fluid, normalised by their values at the location of the undisturbed interface at
$z=0$
, for the exponential profile. Non-dimensional parameters:
$(Bo\rightarrow \infty ,\textit{Fr}=3194.8)$
. In panels (a,b) the red and green curves correspond to the low wavenumber and fastest growing modes from figure 3, respectively. The vertical coordinate
$z$
is scaled here by the coordinate
$z_c$
of the critical layer where the base velocity
$U_{u}(z_{c})$
matches the real part of the phase speed of an unstable mode, i.e.
$c_{r}(k)$
. The blue curve is for the fastest growing mode with
$\delta =0.5$
.
In (3.4) the first term on the right-hand side is the rate of work done by the Reynolds stress and has the physical interpretation of a source term involving energy (per unit volume) extraction from the base state, when it is positive. The second term however, represents the (
$z$
) divergence of the flux of perturbation kinetic energy/volume and acts to redistribute the available perturbation kinetic energy (see discussion in Smyth & Peltier (Reference Smyth and Peltier1989), p.
$3705$
). The vertical variation of Reynolds stress in figure 4(a) shows that, for the fastest growing mode (green curve), the dominant energy extraction from the base state occurs in the region between the critical layer and the interface, i.e.
$0\lt z/z_c\lt 1$
. In contrast, for the low wavenumber KH mode (red curve), a significant proportion of the energy extracted from the base-state shear also occurs in the region where the shear commences
$z\approx \varDelta$
and up to the critical location, i.e.
$z/z_c = 1$
.
3.1.2. Higher density ratio (
$\delta \lt 1$
)
We now move towards a higher density ratio and the high-Froude-number regime. As discussed, the PL profile admits instabilities even for the unstratified case
$(\delta =1)$
due to resonance between the vorticity interfaces at
$z=0$
and
$z=\varDelta$
; see the Boussinesq approximation in Lawrence et al. (Reference Lawrence, Browand and Redekopp1991, Figure 3d) in the Boussinesq approximation). While the growth rates of the fastest growing mode in the exponential model agree reasonably well with the corresponding PL mode (figure 5
a), the phase speed differs noticeably (figure 5
b). The corresponding Reynolds stress and perturbation kinetic energy for the fastest growing mode are plotted in blue in both panels of figure 4. Note the appearance of a peak in the perturbation kinetic energy at the location where the shear layer commences, i.e.
$\varDelta /z_c\approx 3.29$
, we have checked that no such peak is detected for the fastest growing mode at the air–water density ratio (green curves corresponding to
$\delta =0.001$
in figure 4). This change in the nature of the variation of
$c_r/c_{_{\textit{KH}}}$
with wavenumber
$\kappa$
for unstable wavenumbers as the density ratio is varied is seen both in the PL model and in the exponential model in figure 6, where a systematic increase of density ratio up to unity is reported.
The infinite Bond (
$Bo\rightarrow \infty$
), high Froude (
$ \textit{Fr}=3194.38$
) number and high density ratio (
$\delta =0.5$
) case. (a) Growth rate (
$\kappa \tilde {c}_i$
) and (b) phase speed (
$c_r/c_{_{\textit{KH}}}$
) as a function of the non-dimensional wavenumber
$\kappa$
for the KH dispersion relation from (3.1) (pink dashed curve) for the exponential profile (solid blue curve) from (2.12) and for the PL profile from the dispersion relation in (3.2) (‘+’ symbols) for the base states of figures 1 and 2. The green dot indicates the fastest growing mode in both panels. The phase speeds of the PL and the exponential profile are markedly different.
Variation of phase speed with wavenumber for varying density ratios. Here
$Bo\rightarrow \infty ,\;\textit{Fr}=3194.38$
for all cases. All dashed curves are for the exponential model of (2.1).
3.2. Infinite bond (
$Bo\rightarrow \infty$
) and low Froude number
$(\textit{Fr} \sim 1)$
We next turn to the case of infinite Bond number, and moderately low Froude number, examining the effect of the variation of density ratio on the instability mechanism in the exponential profile. Before doing this, it is useful to digress and highlight the instability mechanisms in the PL profile of figure 2(a) within this parameter regime of interest as a function of density ratio; these results will serve as reference for the exponential profile emphasising the role of background curvature. In this context, we note the study by Morland et al. (Reference Morland, Saffman and Yuen1991) who compared free surface waves subject to shear (only in lower layer with
$\delta =0$
) generated by an exponentially decaying velocity profile (infinite depth) against those from a PL profile. The instability for their PL profile (Morland, Saffman & Yuen Reference Morland, Saffman and Yuen1991) was shown to be a H-like instability, arising due to resonance (referred to as ‘collision’ in Morland et al. (Reference Morland, Saffman and Yuen1991)) between a retrograde capillary-gravity mode and a prograde vorticity mode. Here ‘prograde’ and ‘retrograde’ modes are defined (following Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011)) as those that travel leftward (negative velocity) or rightward (positive velocity), respectively, in a frame of reference that moves at the local background speed at the vorticity or density interface; see the explanation around equation (
$25$
) in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011). We see in the next subsection that an analogous H instability occurs in our PL profile of figure 2(a) at a low density ratio and low Froude number.
3.2.1. The instabilities in the PL profile at low
$\textit{Fr}$
: density ratio (
$\delta$
) variation
Here, we interpret the instabilities admitted by the PL profile of figure 2(a) at a low Froude (infinite Bond) number and for a varying density ratio, employing wave resonance ideas reviewed in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011). Figure 7(a) shows the phase speed versus wavenumber (refer to the inset for the narrow range of wavenumbers, which are unstable and shaded) for the PL profile at
$\delta =0.001$
and
$ \textit{Fr}\approx 10$
. The three solid blue lines in the figure are obtained from the numerical solution to the dispersion relation given by (3.2). Recall that the density interface at
$z=0$
(see figure 2
a) supports stable prograde and retrograde surface-gravity waves in the absence of shear (labelled as
$c_{g}^{+}$
(dashed red line) and
$c_{g}^{-}$
(dashed cyan line) in figure 7
a). The vorticity interface when present alone at
$z=\varDelta$
and in the absence of stratification supports a stable retrograde vorticity wave (
$c_v(\varDelta )$
, green dashed line in figure 7
a). Similarly, the vorticity interface at
$z=0$
when present alone supports a stable prograde vorticity wave (
$c_v(0)$
, yellow dashed line); note that one of the blue curves computed from the numerical solution to (3.2) asymptotes to the green dashed line in figure 7(a). However, the other two blue curves from the numerical solution to (3.2) are indistinguishable from the two
$c_{g}$
modes, indicating the negligible effect of shear on the phase speed of these.
(
$Bo\rightarrow \infty ,\;\textit{Fr}\approx 10$
) The normalised phase speed (
$c_r/U_{\infty }$
) variation with wavenumber
$\kappa$
for the isolated prograde (
$c_v(0)$
) and retrograde (
$c_v(\varDelta )$
) vorticity modes as estimated using Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011, Equation 25) and isolated surface-gravity modes with
$c_{g}^{+/-}$
(2.14a
). The three solid blue curves are the roots of the dispersion relation given by (3.2). The grey shaded regions represent unstable wavenumbers. In the inset of panel (a), unstable wavenumbers are indicated as black solid curve. (a) At low
$\delta =0.001$
, the instability is the result of resonance between the prograde gravity-capillary wave and the retrograde vorticity wave. (b) At higher
$\delta =0.9$
, the instability is the result of resonance between the prograde and the retrograde vorticity waves. Consistent with the Howard semi-circle theorem (Morland & Saffman Reference Morland and Saffman1993), unstable wavenumbers have their phase speeds satisfying
$0\lt c_r/U_{\infty } \lt 1$
in both panels.
The inset in figure 7(a) makes it clear that the instability in the PL model (in this regime) arises due to resonance between the prograde gravity wave (
$c_{g}^{+}$
) at
$z=0$
and the slightly modified version of the retrograde vorticity wave at
$z=\varDelta$
(blue curve asymptotic to the green dashed line). Refer to the inset where that the coalescence of these aforementioned curves produces the band of unstable wavenumbers indicated by the shaded region; we conclude that these wavenumbers in figure 7(a) are H unstable following the definition of the H instability (Carpenter et al. Reference Carpenter, Tedford, Heifetz and Lawrence2011). In figure 7(b) a qualitatively different behaviour is seen at a much higher density ratio of
$\delta =0.9$
(
$ \textit{Fr}\approx 10$
). Now the band of unstable wavenumbers (hashed region) arise from resonance between the two vorticity waves at
$z=0$
and
$z=\varDelta$
, as might be intuitively expected from a consideration of the PL configuration, excluding stratification, i.e. at
$\delta =1$
. In figure 7(b), these wavenumbers are KH unstable, propagating nearly at the mean background speed
$U_{\infty }/2$
, consistent with (3.1). It is thus established that for a low Froude number, the PL model of figure 2(a) admits the H instability at a low density ratio, transitioning to the KH instability at a higher density ratio closer to unity. How does inclusion of background profile curvature affect this H instability? We answer this next.
3.2.2. Low density ratio (
$\delta \ll 1$
) – instability in the exponential profile
Having discussed the two instabilities in the PL profile, we now turn to studying the effect of background curvature at the low Froude number and low density ratio regime where the PL profile was studied in the previous subsection. Because of the presence of a critical location where the background profile curvature is non-zero, instead of the H instability we now expect to recover the classic Miles (Reference Miles1957), wind-wave instability, dominated by critical location dynamics. In this context, it is relevant to mention the recent study by Carpenter et al. (Reference Carpenter, Guha and Heifetz2017), which extends the wave interaction mechanism discussed in § 1, in the context of the KH and H instability, to the Miles instability as well. Interestingly, they show that the Miles critical layer instability (at the air–water density ratio) can be interpreted as resonance between a vortex sheet at the critical layer and a surface wave at the density interface.
We note that our choice of low
$ \textit{Fr} \approx 4.8$
in this subsection is in the same ballpark range as Young & Wolfe (Reference Young and Wolfe2014); see their figure
$11$
where their smallest Froude number of
${\sim} 1.4$
corresponds to the air–water oceanic situation. In this parameter regime the fastest growing mode has features of a classic ‘Miles mode’ of instability (Miles Reference Miles1957) dominated by the critical layer (Miles Reference Miles1957; Young & Wolfe Reference Young and Wolfe2014; Kadam, Patibandla & Roy Reference Kadam, Patibandla and Roy2023); figure 8 presents data to support this. Unlike the high Froude cases studied previously, the fastest growing mode in the moderate Froude case (figure 8
a) travels nearly at the surface-gravity wave speed, viz.
$c_{g}$
. Note that
$c_{g}/c_{_{\textit{KH}}}= (({1-\delta ^2})/{\delta ^2} )^{1/2}\kappa ^{-1/2}\textit{Fr}^{-1} \gg 1$
for the air–water density ratio, as for the fast growing mode
$\kappa = O_m(1)$
with
$ \textit{Fr} = O_m(1)$
; see figure 8(a). The fact that at sufficiently low
$\delta$
, unstable, prograde modes travel nearly at the corresponding surface-gravity wave speed in still air/water is originally due to Miles (Reference Miles1957). This has been derived formally by several subsequent studies following Miles’ procedure of expanding the eigenfunction and eigenvalue in the Rayleigh equation, employing a density ratio as a small perturbation parameter, e.g. see below equation (
$3.10$
) in Alexakis et al. (Reference Alexakis, Young and Rosner2002) or figure
$3$
in Young & Wolfe (Reference Young and Wolfe2014) where capillarity has also been taken into account. In what follows, we illustrate several attributes of the fastest growing Miles mode and use these as reference to contrast the change observed, as the density ratio is increased significantly beyond that of air–water.
The infinite Bond (
$Bo\rightarrow \infty$
), low Froude (
$ \textit{Fr} = 4.79$
) number and low (air–water) density ratio (
$\delta =0.001$
) case. (a) Growth rate (
$\kappa \tilde {c}_i$
) and phase speed (
$c_r/c_{\text{g}}$
) as a function of the wavenumber
$\kappa$
. The shaded region is that beyond which it is numerically challenging to detect growth, due to a very small growth rate (
$\lt 10^{-12}$
). (b) Vertical variation of Reynolds stress and perturbation kinetic energy for the fastest growing mode represented by dots in (a) with
$\lambda /z_c\approx 32$
.
Figure 8(b) presents the vertical variation of perturbation kinetic energy and the Reynolds stress for the fastest growing ‘Miles mode’ (dots in figure 8
a). Note the classic near-discontinuous behaviour of the Reynolds stress at the critical location
$z/z_c=1$
. This is consistent with Lin (Reference Lin1954) (their equation (
$3.10$
)) who showed, following Tollmien (Reference Tollmien1935), that the Reynolds stress
$\tau (z) \equiv \overline {-\rho u w}$
across the critical location
$z_c$
(at which
$D^2U_u \lt 0$
and
$DU_u\neq 0$
) jumps by the amount
\begin{eqnarray} \tau (z=z_c^{+}) - \tau (z=z_c^{-}) = \dfrac {\pi }{k}\left (\dfrac {D^2U_u}{|DU_u|}\right )_{z=z_c}\rho \overline {w_c^2} \end{eqnarray}
for a singular neutral mode, i.e. a neutral solution to the Rayleigh equation for which there is a critical layer with
$D^2U_u(z_c), DU_u(z_c)\neq 0$
; refer to Tollmien’s inviscid solution in equations (
$22.25$
) and (
$22.26$
) in Drazin & Reid (Reference Drazin and Reid2004) or around equation (
$2.7$
) in Stewartson (Reference Stewartson1981). In (3.5),
$w_c$
is the vertical perturbation velocity at the critical location
$z_c$
. Although figure 8(b) has been plotted for the fastest growing Miles mode (and not the neutral mode to which (3.5) applies), even for the fastest growing Miles mode
$c_i/c_r \approx 0.001 \ll 1$
, leading to the near-discontinuous behaviour seen in the Reynolds stress (see further discussion on this around (3.6)). In figure 8(b) the shear layer is between
$0 \lt z/z_c \lt \varDelta /z_c\approx 3.4$
and it is clear that the energy extraction (work done by the Reynolds stress term in (3.4)) from the base state occurs only within the thin region between the critical layer and the interface. We also contrast figure 8(b) (where
$ \textit{Fr} \sim O(1)$
) with its counterparts presented earlier in figure 4
$(\textit{Fr}\gg 1)$
. Unlike in figure 8(b), sizable contributions to both perturbation kinetic energy and Reynolds stress in figure 4 were noted in the region between the start of the shear layer and the critical location (
$1 \lt z/z_c \lt \varDelta /z_c$
).
Figure 9 shows the (real part) of the perturbation streamfunction for the fastest growing mode in figure 8(a). A clear bend in the streamfunction contours is visible at
$z=z_c$
(figure
$13$
(a) in Young & Wolfe (Reference Young and Wolfe2014) notes the same behaviour), indicating the phase change in the perturbation streamfunction, characteristic of the critical layer. This in turn generates the phase change in perturbation pressure across
$z/z_c=1$
, causing the perturbation pressure (see (2.7) relating the streamfunction
$\varPsi (z)$
to perturbation pressure) to have a component in phase with
${\partial \eta }/{\partial x}$
. As lucidly explained in Miles (Reference Miles1959), at the air–water density ratio the pressure component in-phase with
${\partial \eta }/{\partial x}$
is more efficient than the one in-phase with
$-\eta$
in generating instability; this leads to a much lower threshold of instability for the critical wind speed when gravity and capillarity are both accounted for in the exponential model as compared with the corresponding discontinuous KH model (see figure
$1$
in Kadam et al. (Reference Kadam, Patibandla and Roy2023)).
Figure 10(a) notes the variation of the perturbation pressure
$p_u(x,z=0^{+},t=0)$
. If the perturbation pressure just above the interface (air side) were to be entirely in phase with
${\partial \eta }/{\partial x}$
, the phase difference with
$\eta$
(black dashed curve in figure 10
a) would be
$\pi /2$
, as for the Jeffrey model discussed earlier (see discussion around their equation (
$1.1$
) in Bonfils et al. (Reference Bonfils, Mitra, Moon and Wettlaufer2022)). For the fastest growing Miles mode in figure 10(a), the phase angle (
$\approx 72^{\circ }$
) is much closer to
$\pi /2$
compared with the angle
$\pi$
of the KH model (red dashed curve) – we note in passing that the angle
$\pi$
for the KH model changes with increasing
$\delta$
. Figure 10(b) notes the vertical variation of the kinetic energy flux (second term in (3.4)) for the fastest growing Miles mode, note the sign change in the flux around
$z_c$
, indicating that a thin layer centred around the critical layer redistributes the bulk of the generated perturbation kinetic energy towards the interface and a much smaller fraction away from it. In the budget for wave energy integrated vertically over the entire domain, the integration of this flux-divergence term in (3.4) equals the potential energy of the interface (purely gravitational in the current case); see Young & Wolfe (Reference Young and Wolfe2014, Equation (3.2)).
(a) Variation of perturbation pressure for the fastest growing mode at
$z=0^{+}$
with horizontal coordinate
$x$
in figure 8(a). Here
$ \textit{Fr} = 4.79$
,
$\delta =0.001$
,
$Bo\rightarrow \infty$
. The dotted black line represents the real part of the perturbed interface given by
$\eta$
in (2.4). The blue curve corresponds to perturbation pressure for the fastest growing Miles mode in figure 8(a). The red dashed line represents the perturbation pressure
$p(z=0^{+})$
obtained from the KH model of figure 2(b) for
$\delta =0.001$
and (in CGS units)
$U_{\infty }=1000$
. (b) The vertical variation of the second term on the right-hand side of (3.4).
We next turn to examining the effect of increase of density ratio (
$\delta$
) on the fastest growing mode in the low-Froude-number regime of current interest. In order to rationalise the behaviour of the dispersion relation of (2.12), it is useful to compare the numerical solution to this against the approximation obtained in the low density ratio limit by expanding
$\tilde {c}(\kappa ;\delta )$
as
$\tilde {c}=\tilde {c}_0(\kappa ) + \delta \tilde {c}_1(\kappa ) + \mathcal{O}(\delta ^2)$
for fixed Froude and Bond numbers.
For the zero surface-tension limit relevant here, we use the expression upto
$\mathcal{O}(\delta )$
obtained by Carpenter et al. (Reference Carpenter, Guha and Heifetz2017) (see equation (
$15$
) in their study). This, expressed in our notation is
\begin{align} \tilde {c} \approx \tilde {c}_0 + \delta \tilde {c}_1 = \tilde {c}_0 + \dfrac {\delta }{2\kappa }\big [\textit{Fr} + \tilde {c}_0 \big (\varXi (\tilde {c}_0, \kappa ;\textit{Fr} ) - \kappa \big )\big ], \end{align}
with (
$\tilde {c}_0 \neq \textit{Fr}$
)
\begin{align} \tilde {c}_0 \equiv \kappa ^{-1/2},\; \varXi (\tilde {c}_0, \kappa ;\textit{Fr}) &\equiv -\kappa + \left (\dfrac {1}{1+2\kappa }\right )\left (\dfrac {\textit{Fr}}{\textit{Fr}-\tilde {c}_0}\right ) \nonumber\\[5pt]&\quad\dfrac {_2F_1\left (\alpha +1,\beta +1,\gamma +1;\dfrac{\textit{Fr}}{\textit{Fr}-\tilde {c}_0}\right )}{_2F_1\left (\alpha ,\beta ,\gamma ;\dfrac{\textit{Fr}}{\textit{Fr}-\tilde {c}_0}\right )}, \end{align}
the definitions of
$\alpha ,\beta$
and
$\gamma$
being available below (2.8). We remind the reader that the perturbative expansion (3.6) assumes that
$\delta \rightarrow 0,\;\textit{Fr}=\mathcal{O}(1)$
, thus excluding the high Froude (
$ \textit{Fr}\rightarrow \infty$
) limit studied earlier.
Figures 11(a) and 11(b) compare the numerical solutions to the complete dispersion relation from (2.12) (solid lines) with the approximation in (3.6) (big and small dots). One notes the increase in growth rate beyond the linear dependence (for
$\delta =0.5$
) on
$\delta$
predicted by (3.6). Note from figure 11(b) that as
$\delta$
is increased, the correction to the phase speed (Doppler shifting) compared with
$\tilde {c}_{g}$
(in still air–water
$\tilde {c}_0 \approx \tilde {c}_g$
) for the fastest growing mode becomes significant. We assert that this increase in
$\tilde {c}_r$
in figure 11(b) indicates a change in the nature of the instability for the fastest growing mode as
$\delta \rightarrow 0.5$
. To support this argument, we compare the growth rate (
$\kappa \tilde {c}_i$
) at
$\delta =0.5$
for the exponential model (solid red curve in figure 11
a) with the PL model (dashed black line in figure 11
a), noting the good agreement between the two around a region of wavenumbers close to the fastest growing mode; note also the match in the phase speed in figure 11(b), although there is evident a mismatch for wavenumbers centred around the fastest growing mode (dashed black curve and solid red curve). We noted in the last section that at a low Froude numbers and low density ratio, the instability in the PL model is the H instability. It is further shown in the supplementary material (§ 2) that the H instability persists until
$\delta \approx 0.5$
at low
$\textit{Fr}$
for the PL model (although KH contributions are also apparent; see figure 2(a) in the supplementary material). Given the good agreement seen in the growth rates between the PL and the exponential model in figure 11(a) at
$\delta =0.5$
, we thus intuitively anticipate that the exponential model also admits the H instability around
$\delta =0.5$
.
(a) Growth rate and (b) phase speed variation with
$\kappa$
for the exponential profile. Non-dimensional parameters:
$ \textit{Fr} = 4.79$
,
$Bo\rightarrow \infty$
. In both panels, solid lines are from the full dispersion relation in (2.12) while large and small dots are from the approximation (3.6). The dashed black line in both panels is from the PL model in figure 2(a) for the same non-dimensional parameters as above. The curve with stars and dots in panel (a) represents a comparison with formulae A31 in Appendix A of Alexakis et al. (Reference Alexakis, Young and Rosner2004b
), valid for small Froude numbers, large wavenumbers and assuming that
$c_i \ll c_r$
. The black dot in panel (b) is for the fastest growing mode for the exponential model (solid red curve).
Our aforementioned assertion of a change from the Miles instability with increasing
$\delta$
is further supported by the vertical variation in Reynolds stress (figure 12
a) and perturbation kinetic energy (figure 12
b) plotted for the fastest growing mode at every
$\delta$
. Note in figure 12(a) the near-discontinuous jump for
$\delta =0.001$
at the critical layer noted earlier, giving way at
$\delta =0.5$
, to a smooth transition through
$z=z_c$
bereft of any distinguishing signature at this location; notice also that around
$\delta =0.5$
, significant contributions appear to the Reynolds stress from the range
$z/z_c\gt 1$
unlike the Miles case discussed earlier in figure 8(b). Why does the vertical variation of Reynolds stress in figure 12(a) for the fastest growing mode show a qualitative change as
$\delta$
is increased? We address this question in § 3.4; a noteworthy observation from this section is as
$\delta$
increases from
$0.001$
to
$0.5$
in figure 12(a), the ratio
$c_i/c_r$
(for the fastest mode) increases from
$\approx 0.001$
to
$\approx 0.42$
.
Vertical variation of (a) Reynolds stress and (b) perturbation kinetic energy for the fastest growing mode for every
$\delta$
. The non-dimensional parameters are
$ \textit{Fr} = 4.79$
and
$Bo\rightarrow \infty$
. In panel (b), note the transition from having a dominant peak at the critical layer to one only at the interface.
3.3. Effect of capillarity
We discuss the effect of capillarity in this section, restricting ourselves to the low-Froude-number regime where interesting transitions were recorded. As remarked earlier, the fastest growing mode studied so far are long waves (compared with the shear layer thickness as well as the air–water capillary length) and, thus, have subdominant contribution from surface tension. For consistency, we only present high-Bond-number results here (
$Bo=\mathcal{O}_m(88)$
) relevant to geophysical and astrophysical situations. A discussion of low-Bond-number instabilities is eschewed because these can involve shorter waves where the viscous contribution to Reynolds stresses as well as the Tollmien–Schlichting instability can become important (Miles Reference Miles1962). In addition, for shorter wavelengths, one also expects the inviscid rippling instability (Stern Reference Stern1973; Morland, Saffman & Yuen Reference Morland, Saffman and Yuen1991; Shrira Reference Shrira1993; Young & Wolfe Reference Young and Wolfe2014; Kadam, Patibandla & Roy Reference Kadam, Patibandla and Roy2023) when shear is also included in the water layer, as discussed earlier for Morland et al. (Reference Morland, Saffman and Yuen1991); the growth rates for this can be significantly more than that of the longer wavelength modes, especially at a low density ratio. An important difference here, compared with the infinite-Bond-number case is the existence of a critical Froude number
$ \textit{Fr}_c\gt 0$
for instability. This is easily rationalised at
$\delta \ll 1$
. Recall that prograde, capillary-gravity waves propagating at the interface of still water/air (which are destabilised by shear in the upper fluid) have a non-zero, minimum phase speed (unlike the case of surface-gravity modes), thereby leading to a critical Froude number based on this minimum speed.
Figure 13(a) presents the growth-rate curves for varying Froude numbers at a low density ratio (
$\delta =0.001$
). We note that
$ \textit{Fr}_c\approx 0.5$
and the case with
$ \textit{Fr}=0.6$
in the figure 13(a), is very close to the instability onset. The
$ \textit{Fr}=4$
case (figure 13
b) shows that the fastest growing mode in this case travels at almost
$c_{gc}$
. Examination of the Reynolds stress and perturbation kinetic energy for
$ \textit{Fr}=4$
in figures 14(a) and 14(b) establish that these correspond to a classic Miles mode with distinct (near) discontinuities and peaks, respectively, at
$z_c$
. Note the interesting behaviour of the fastest growing mode for
$ \textit{Fr}=0.6$
in these figures. There is a sharp near discontinuity at
$z_c$
(cf. figure 14
a) for
$\tau \equiv -\overline {\rho u w}$
but no corresponding signature at
$z_c$
in figure 14(b) for the kinetic energy. To conclude this section, figures 15(a) and 15(b) note the
$z$
variation of the Reynolds stress and perturbation kinetic energy (for the fastest growing mode in each case) with increasing density ratio; similar to figure 12 a clear change in the nature of variation of both metrics at
$z/z_c=1$
is observed.
(a) Growth rate versus wavenumber for varying
$ \textit{Fr}\gt \textit{Fr}_c\approx 0.5$
. For both panels,
$\delta =0.001,$
$Bo = 87.813$
. The shaded region indicates the range of wavenumbers, beyond which the growth rate is numerically very small. (b) Growth rate and phase speed for
$ \textit{Fr} = 4$
. The Froude number (
$ \textit{Fr} = 4$
) here has approximately the same value as in figure 8(a), which is for an infinite Bond number compared with
$Bo=87.813$
here.
(a) Reynolds stress and (b) kinetic energy variation for fastest growing mode. Here
$\delta =0.001,$
$Bo = 87.81$
. In panel (b), note the non-monotonic behaviour of the peak at the critical layer with increasing
$\textit{Fr}$
.
(a) Vertical variation of Reynolds stress and (b) kinetic energy for the fastest growing mode. Here
$ \textit{Fr}=4,\;Bo= 87.81$
. In panel (b), note the transition from having a dominant peak at the critical layer to one at the interface, with increasing
$\delta$
.
3.4. Why does the Reynolds stress become continuous at
$z=z_c$
with increasing
$\delta$
?
We return to the question: Why does the Reynolds stress for the fastest growing mode show qualitative change with increasing
$\delta$
? This behaviour is observed both in figure 12(a) (infinite Bond number) and in figure 15(a) (
${Bo}= 87.81$
). It turns out that this can be understood from Lin’s relation for the vertical variation of Reynolds stress
$\tau (z)$
for an unstable mode with phase speed
$c = c_r + ic_i$
. Using background velocity as the independent variable (in place of
$z$
) this is
\begin{eqnarray} \dfrac {\text{d}\tau }{\text{d}U_u} = F(U_u)\left [\dfrac {c_i}{\big (U_u(z)-c_r\big )^2+c_i^2}\right ], \end{eqnarray}
where
$F(U_u(z))\equiv ({1}/{k}) ({\text{d}Z}/{\text{d}U_u} )\overline {\rho _uw^2(z)},\; Z(z)\equiv DU_u$
. For our exponential profile (2.1),
$DU_u$
can be rewritten in terms of
$U_u$
as
\begin{eqnarray} Z \equiv DU_u= \dfrac {U_{\infty }}{\varDelta }\left [1-\dfrac {U_u}{U_{\infty }}\right ]. \end{eqnarray}
Using (3.9) in (3.8) simplifies to the following non-dimensional form:
\begin{eqnarray} -\left (\dfrac {1}{\rho _u\overline {w^2(z)}}\right )\left (\dfrac {\text{d}\tau }{\text{d}\tilde {U}_u}\right ) = \dfrac {1}{\kappa }\left [\dfrac {\tilde {c}_i}{\big (\tilde {U}_u(z)-\tilde {c}_r\big )^2+\tilde {c}_i^2}\right ]. \end{eqnarray}
Here
$\tilde {U}_u \equiv U_u/\sqrt {g\Delta }$
. Importantly, even for the fastest growing Miles mode in figure 8,
$c_i/c_r \ll 1$
. Now using the observation that
$\displaystyle \lim\nolimits_{y\rightarrow 0}({1}/{\pi }) [ {y}/({(x-x_0)^2+y^2} )]\rightarrow \delta (x-x_0)$
(Lin (Reference Lin1954)),
$\delta (x)$
being the Dirac delta function, we see that the quantity inside the square bracket in expression (3.10) for the fastest growing Miles mode has a near delta function behaviour and this generates the jump in the Reynolds stress in figure 8(b).
To visualise this, in figure 16(a) we plot the right-hand side of (3.10) as a function of the vertical coordinate (
$z$
) for the fastest growing mode in each case. As expected, the low Froude number (
$ \textit{Fr}=4$
) and low density ratio (
$\delta =0.001$
) has the signature of a Dirac delta function (blue dashed curve) leading to the jump behaviour of
$\tau$
in figure 15(a). For the dashed green and red curves in figure 16(a), notice the widening as the ratio
$c_i/c_r$
increases significantly. For
$\delta =0.5$
(red curve with
$c_i/c_r=0.37$
), a large fraction of the area under the curve comes from
$z/z_c\gt 1$
with a smaller contribution from
$z/z_c \lt 1$
. This generates the observed smooth behaviour of
$\tau$
through
$z_c$
in figure 15(a). Figure 16(b) justifies the transition seen with increasing density ratio from a slightly different viewpoint. Recall that Tollmien’s inviscid solutions have a logarithmic singularity when the curvature at the critical location is finite (with
$ (DU_u )_c\neq 0$
). Figure 16(b) computes the phase shift in such a term in a small region around
$z/\text{Re}{ (z_c )}=1$
(
$z_c$
is a complex number with
${\rm Im}(z_c)\gt 0$
, and is computed for the fastest growing mode using
$U_u(z_c)=c(k)$
) with
$c(k)$
being complex for the fastest growing mode. It is seen that a qualitative change occurs with increasing
$\delta$
, viz. the abrupt jump in phase seen at
$z/\text{Re} (z_c )$
for the argument of the logarithmic term gives way to a smooth variation in figure 16(b). These aforementioned arguments pinpoint the source of the differences seen in the vertical variation of
$\tau (z)$
between the low and high density ratio cases.
(a) Vertical variation of the right-hand side of (3.10) for the exponential profile for different density ratios. In each case, the fastest growing mode has been chosen. (b) Change in the phase of
$(z-z_c)\log (z-z_c)$
. This term appears in Tollmien’s singular inviscid solution; see expressions (
$22.25$
) and (
$22.26$
) in Drazin & Reid (Reference Drazin and Reid2004). For (b),
$ \textit{Fr}=4,Bo= 87.81$
.
3.5. Distinguishing between the H and KH instability for asymmetric profiles
Given the good comparison between the growth rates at
$\delta =0.5$
of the PL profile and the exponential profile in figure 11(a), we hypothesise that at low Froude and high Bond numbers, the nature of the fastest growing mode in the exponential profile undergoes a transition – it changes from being dominated by the critical location for
$\delta \ll 1$
to the H instability for
$\delta \approx 0.5$
and eventually the KH instability at even larger
$\delta$
. We emphasise that the KH to H transition seen clearly in the PL model occurs smoothly, in particular at
$\delta =0.5$
there seems to be contributions from both (see § 2 of the supplementary material). For the asymmetric, exponential background profile of figure 1, the vertical variation of the perturbation kinetic energy and the Reynolds stress do not show any qualitative change as
$\delta$
is changed from
$0.5$
to
$0.9$
, this may be ascertained from figures 15(a) and 15(b). This lack of distinction however is not indicative; even for symmetric curved background profiles, where the KH and H instability can be distinguished clearly based on phase speed, the corresponding Reynolds stress and the perturbation kinetic energy do not show qualitative differences in their vertical distribution. Smyth & Peltier (Reference Smyth and Peltier1989) comment on this in the first paragraph on p.
$9$
(note that the word inside the bracket has been added by us), ‘the two (Holmboe) modes illustrated in figures
$\ldots$
are qualitatively very similar to KH modes
$\ldots$
’, emphasising the challenge of distinction between the two based on these metrics. One may, for example, use the diagnostic proposed in Carpenter et al. (Reference Carpenter, Balmforth and Lawrence2010) to partition the growth rate from the exponential profile into a H and a KH contribution; see their figure 9. We do not do this here, instead taking up numerical simulations of the incompressible, Euler’s equations with surface tension and gravity; in the nonlinear regime captured through these simulations, clear differences are seen to emerge with increasing density ratio. These simulations clearly distinguish between the H and KH instability, as described next.
4. Numerical simulations – linear and nonlinear regime
In this section we present results from numerical simulations obtained by solving the incompressible, Euler’s equations using the open-source code Basilisk (Popinet & Collaborators Reference Popinet2013–Reference Popinet2025). We restrict ourselves to the low-Froude-number and high-Bond-number regime discussed earlier. To be consistent with the inviscid stability model discussed so far, we solve the incompressible Euler’s equations with gravity and surface tension, in Basilisk, with an interface separating two fluids. The equations that are solved read as
\begin{align} \boldsymbol{\nabla .u} & =0, \\[-12pt] \nonumber \end{align}
\begin{align} \dfrac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{\nabla }.(\boldsymbol{u\otimes u}) & = - \dfrac {\boldsymbol{\nabla }p}{\rho } + g + \dfrac {T}{\rho }\kappa \delta _s \boldsymbol{n}, \\[-12pt] \nonumber \end{align}
\begin{align} \dfrac {\partial f}{\partial t}+\boldsymbol{\nabla .}(f\boldsymbol{u}) & =0,\\[10pt]\nonumber \end{align}
where
$\boldsymbol{u},\,p,\,\kappa$
,
$T$
and
$f$
are the velocity field, pressure field, interface curvature, surface tension and the volume-fraction field, respectively. The velocity field is
$\boldsymbol{u} = (u,w)$
, where
$u$
and
$w$
are horizontal and vertical components, respectively. Basilisk implements the volume-of-fluid algorithm in a one fluid formulation for reconstruction and advection of the interface separating the two fluids. This is done along with time marching of the accompanying velocity and pressure fields. A volume-fraction field
$f(x,z,t)$
is defined that takes values
$\in [0,1]$
. Specifically, in all computational cells containing only the upper or lower fluid,
$f$
takes the value
$0$
or
$1$
. Thus, by definition the interface is present in those cells where
$0 \lt f \lt 1$
, the value representing the fraction of the cell that is occupied by one of the two fluids. Given an
$f$
field at any time instant, a reconstruction algorithm computes the interface using a PL approximation for each computational cell where
$0 \lt f \lt 1$
. The
$f$
field is advanced in time employing a discrete version of (4.1c
); volume fluxes through cell walls being computed geometrically ensuring high accuracy of volume conservation; see Tryggvason et al. (Reference Tryggvason, Scardovelli and Zaleski2011) and references therein. The density and viscosity fields are computed as a weighted average of their respective values in the two phases, employing the local value of
$f$
as the weight.
It is important to note that the volume-of-fluid algorithm constructs a sharp representation of the interface. This representation is not just for visualisation, but is also explicitly employed to estimate the volume fluxes for each computational cell, while time stepping (4.1c
). Surface tension is incorporated as a body force into the momentum equations and curvature is computed using a generalised, height-function algorithm (see §
$5$
in Popinet & Collaborators Reference Popinet2013–Reference Popinet2025). The code has been benchmarked extensively, particularly at high Reynolds numbers including inviscid test cases (Popinet Reference Popinet2003; Kayal et al. Reference Kayal, Sanjay, Yewale, Kumar and Dasgupta2025). We have, in particular, extensively employed the adaptive grid feature of this code in our current simulations. The criteria for adaptivity is based on changes to the volume-fraction field
$f$
, curvature
$\kappa$
and the velocity field, when subject to local grid refinement. The computational domain length has been chosen to match the wavelength
$\lambda$
of the fastest growing mode, obtained from theory. The smallest local grid size is
$\text{(domain size)}/2^{13}=\lambda /2^{13}$
. This level of adaptivity is referred to as
$2^{13}$
, for example, in figure 19(a) and similarly in other figures. In some of the simulation results that follow, qualitative features are observed to remain unchanged under grid refinement, although quantitative changes are apparent in certain cases. A sample simulation file for the simulation reported in figure 20 is also provided (Kumar Reference Kumar2025).
Figure 17(a) depicts the simulation domain along with boundary conditions while figure 17(c) shows a sample adaptive grid employed. The initial conditions in our simulations are obtained from linear theory presented in § 2 (figure 17
b). In both upper and lower fluids, we initialise the velocity field using a velocity field represented as base state (2.1) plus perturbation (2.3). The perturbation field is obtained from the analytically computed eigenfunctions in (2.11). Note that there is only one free parameter in our linear theory, viz. the value of
$\eta _0$
. The value of this is chosen such that the steepness
$k\eta _0=0.02$
for all simulations, consistent with the linear theory assumed (small steepness). The interface is initialised using (2.4), employing the fastest growing wavenumber in each case, as computed from the dispersion relation in (2.12) for a Fourier mode travelling rightwards.
(a) A schematic of the computational domain along with boundary conditions implemented. (b) Vertical variation of fastest growing eigenmodes and background shear. A superposition of background shear with the fastest eigenmode with a small amplitude is used as the initial velocity field in all simulations. Here
$(\hat {U},\hat {u},\hat {w}) \equiv (U/U_{\infty }, u/U_{\infty },w/U_{\infty })$
. (c) Sample adaptive grid.
Figure 18(a) shows simulations for the case for
$ \textit{Fr}=4, Bo=87.81$
and
$\delta =0.01$
. The dimensional parameters for figure 18(a) are
$U_{\infty } = 199.57$
cm s−1,
$\varDelta =2.54$
cm,
$\rho _l=1$
gm cm
$^{-3}$
and
$\rho _l=0.01$
gm cm
$^{-3}$
. For these parameters, the fastest growing mode in the exponential profile has wavelength
$\approx 20.9$
cm. We have seen earlier in figure 15(a) that the Reynolds stress vertical distribution with a sharp jump at the critical location is nearly identical for
$\delta =0.001$
(air–water) as well as
$\delta =0.01$
, indicative of the critical location dominating the instability. In figure 18(a) we present the evolution of the interface from theory and simulations for this Miles regime. Note that time is measured in units of the duration (
$\tau$
) taken for the wave to travel its wavelength; notice the slow growth in figure 18(b). As the wave amplitude significantly exceeds its initial value, a mismatch with linear theory is apparent by
$t/\tau =5$
(five wave periods). We note that simulations for the air–water density ratio (
$\delta =0.001$
) were also attempted; however, due to the extremely small growth rate, long simulations are necessary. Computational inaccuracies in such long runs prevented detection of any exponential growth and, thus, are not reported here. Figure 18 is the smallest density ratio (
$\delta =0.01$
) where we were able to detect growth and compare against linear theory, while being in the critical layer dominated regime. Figure 18(b) reports the interface displacement at the centre of the domain and a slow exponential growth along with a reasonably good match with theory. Note that simulations were carried out at three adaptive resolutions:
$2^{11}$
, for example, represents the upper limit of the number of grid points along each direction. A reasonably good grid convergence is apparent in figure 18(a).
Panels (a i–a iv) compare the interface shape from theory and simulations. Non-dimensional parameters:
$Bo=87.81,\; \textit{Fr}=4,\;\delta =0.01$
and
$\varDelta /\lambda =0.12$
(rounded to two decimal places). (b) The time signal obtained by tracking interface displacement at
$x=0$
(horizontal midpoint of the domain). Here
$\tau$
is the time taken for a Fourier mode to travel one wavelength.
Figure 19 is at the same Froude number as figure 18 but for a density ratio one hundred times that of air–water. We recall that even at these elevated values of
$\delta$
, the Reynolds stress for the fastest growing mode in figure 15(a) retains a strong signature at the critical layer similar to the Miles instability for air–water. An interesting phenomena is apparent at
$t/\tau = 5.6$
; see the inset in figure 19(a v). A series of ripples developed on the backward face of the parent wave – some grid sensitivity is noted in their wavelengths. The inset to this subpanel describes the interface at a slightly later time
$t/\tau =6.5$
where the ripples are more clearly seen. These are reminiscent of ripples on capillary-gravity Stokes waves; see figure
$2$
in Hung & Tsai (Reference Hung and Tsai2009) for the vortical structure of such ripples. For comparison, we have provided in panel (c), an irrotational, left to right propagating Stokes wave (
$Bo=384.615,\; \textit{Fr}=0.4389$
), computed in the limit of
$\delta \rightarrow 0$
with gravity and surface tension. This is a nonlinear, travelling wave solution to the potential flow equations and qualitatively similar structures are seen on the Stokes wave with a marked asymmetry between the forward and the backward face. Qualitative similarities may be seen between panel (c) and the inset in panel (a v) of figure 19. We have also compared numerical simulations for
$\delta =0.01$
(figure 18) and for
$\delta =0.1$
(figure 19) in the nonlinear regime upto the same value of
$t/\tau _{_G}\approx 3.5$
(data not shown here). Here
$\tau _{_G}\equiv (kc_i )^{-1}$
is the time scale of the growth of the instability (which depends on the density ratio) and is distinct from
$\tau = 2\pi /(kc_r)$
(as is done in figures 18 and 19). We observe that the finite-amplitude wave for
$\delta =0.01$
has a far smaller steepness and exhibits tiny, ripple-like structures on its surface, compared with that for
$\delta =0.1$
where the finite-amplitude wave is significantly steeper (slope) and displays larger-amplitude ripples (upto grid precision).
Panels (a i)–(a v) compare the interface shape from theory and simulations. Non-dimensional parameters:
$Bo=87.81,\; \textit{Fr}=4,\;\delta =0.1$
and
$\varDelta /\lambda =0.12$
(rounded to two decimal places). (b) Growth with time at
$x=0$
. Notice the mismatch at larger time. (c) A finite-amplitude, capillary-gravity Stokes wave containing finer-scale ripples comoving with the wave. The shape is computed by solving the nonlinear, potential flow equations including gravity and surface tension following Shelton et al. (Reference Shelton, Milewski and Trinh2025). For this nonlinear wave which is obtained by evolving the gravity-capillary Stokes wave in time, the initial steepness is
$0.24$
,
$ \textit{Fr} \equiv {c_p}/{\sqrt {g\lambda }}=0.4389$
and
$Bo \equiv {\rho g \lambda ^2}/{T} = 384.615$
, with
$\delta =0$
(no upper fluid). The shape of the wave as depicted here is obtained by propagating the initial Stokes wave and taking a snapshot at
$t/\tau _S=2$
. Here
$\lambda ,c_p$
and
$\tau _S$
represent the wavelength, nonlinear phase speed and the time taken to travel one wavelength respectively for the capillary-gravity Stokes wave.
As anticipated, a clear signature of differences in the nonlinear regime is seen between figures 20 and 21, as density is increased from
$\delta =0.5$
to
$\delta =0.9$
, respectively. These differences are further seen readily in the collage of images provided in figure 22, all of them approximately at the same value of
$t/\tau$
. One notes transitions in the nonlinear regime around
$\delta \approx 0.3$
and
$\delta \approx 0.7$
. In figure 20 we particularly note the formation of a cusp-like structure at the crest that emits spume droplets (panel d) (Deike Reference Deike2022); we treat this as a signature of the H instability. Compare figure 20(d) with the shape of the interface extracted from Lawrence et al. (Reference Lawrence, Haigh and Zhu1998) in figure 23(a) that is an asymmetric H instability. With a further increase in density ratio
$\delta =0.9$
, we recover the classic KH spirals in figure 21(d), comparing this to figure 23(b), symptomatic of the asymmetric KH instability. Figures 20 and 21, as well as the collage of images in figure 22, thus show clear evidence that the nonlinear deformation regime is quite different for
$\delta =0.5$
compared with
$\delta =0.9$
. While the former bears similarity to the asymmetric H instability of Lawrence et al. (Reference Lawrence, Haigh and Zhu1998) in figure 23, the latter is very much the asymmetric KH instability in the same figure. We thus validate our hypothesis that the fastest growing mode in the exponential profile indeed undergoes a transition from the Miles to the H instability and further to the KH instability, as the density ratio is increased.
(a,b) Comparison of the interface shape with linear theory. (c,d) Finite-amplitude regime with H-like waves emitting spume-like droplets (Deike Reference Deike2022) from the tip. Here
$Bo=87.81,\; \textit{Fr}=4,$
$\delta =0.5$
.
(a,b) Comparison of the interface shape with linear theory. (c,d) Finite-amplitude regime, note the formation of classic KH spirals in the last two panels. Here
$Bo=87.81,\; \textit{Fr}=4,\;\delta =0.9$
.
Results are shown for (a)
$t/\tau =2.863$
,
$\varDelta /\lambda =0.12$
(rounded to two decimal places); (b)
$t/\tau =2.897$
,
$\varDelta /\lambda =0.12$
; (c)
$t/\tau =2.434$
,
$\varDelta /\lambda =0.12$
; (d)
$t/\tau =2.278$
,
$\varDelta /\lambda =0.12$
; (e)
$t/\tau =2.237$
,
$\varDelta /\lambda =0.12$
; (f)
$t/\tau =2.270$
,
$\varDelta /\lambda =0.12$
; (g)
$t/\tau =2.530$
,
$\varDelta /\lambda =0.13$
; (h)
$t/\tau =2.489$
,
$\varDelta /\lambda =0.13$
; (i)
$t/\tau =2.573$
,
$\varDelta /\lambda =0.13$
. Non-dimensional parameters:
$Bo=87.81,\; \textit{Fr}=4$
. The final set of images at
$\delta$
close to unity bear remarkable similarity to the experimental images in figure
$14.3$
of Cushman-Roisin & Beckers (Reference Cushman-Roisin and Beckers2011). We thank an anonymous referee for pointing us to these images.
Data extracted from figure 9 of Lawrence et al. (Reference Lawrence, Haigh and Zhu1998). Compare panels (a) (asymmetric H) and (b) (asymmetric KH) of this figure with panel (d) of figures 20 and 21, respectively.
In order to establish the qualitative validity of these inviscid simulations, we have also compared these against viscous simulations employing a realistic value for the viscosity of the fluids, at the density ratios of interest. To account for viscosity, the viscous term
$\nu \boldsymbol{\nabla} ^2\boldsymbol{u}$
is added to the right-hand side of (4.1b
). This addition converts the system described by (4.1) into the two-dimensional, Navier–Stokes equations, where
$\nu$
is the kinematic viscosity. Two fluid combinations (water and perfluorodecalin) with a density ratio of
$0.5$
and fresh water and brine with a density ratio of
$0.9$
are considered. We emphasise that these viscous simulation results are preliminary and are reported for consistency – a careful comparison of the solution to the Rayleigh equation against the numerical solution to the Orr–Sommerfeld equations for the exponential profile is necessary and is currently underway. The viscous simulations have been carried out employing free-slip boundaries, thus ensuring that boundary layers do not develop at the computational boundaries. These simulations are reported in figure 24(a,b), comparing these against their inviscid counterparts; a good comparison in general is apparent.
Comparison of inviscid and viscous simulations at the same
$Bo$
and
$\textit{Fr}$
. (a) Inviscid simulation (as in figure 20) versus a viscous case with perfluorodecalin (
$\mu =5.01$
mPa s, heavier phase) and water (
$\mu =1$
mPa s, lighter phase), with density ratio
$\delta = 0.5$
. (b) Inviscid simulation (as in figure 21) versus a viscous case with brine (
$\mu =1.05$
mPa s, heavier phase) and fresh water (
$\mu =1$
mPa
$\boldsymbol{\cdot}$
s, lighter phase), with
$\delta = 0.9$
. For the viscous results, we solve the incompressible Navier–Stokes equations with kinematic viscosity
$\nu \gt 0$
in both fluids. Viscosity has a relatively minor effect within the time window considered.
5. Effect of change of
$\delta$
on the background profile: the double exponential profile
So far, we have considered inviscid linear stability of the background shear and density profiles in (2.1), with shear only in the upper fluid. Recall that within the inviscid approximation, any shear profile in either fluid (upper or lower) along with a flat density interface (
$\eta =0$
) with arbitrary
$\rho _u,\rho _l \gt 0$
and a hydrostatic pressure field is an exact solution to the inviscid momentum and mass equations, as well as the boundary conditions (see § 3 in the supplementary material). In experiments or naturally occurring situations involving viscous, immiscible fluids, one however expects that as the density ratio
$\delta \rightarrow 1^{-}$
(such as in this study), the viscous shear stress exerted by the upper fluid on the lower one generates significant momentum transfer to the lower fluid. (We thank an anonymous referee for pointing this out). For the air–water configuration, the importance of viscous shear in either fluids was studied by Valenzuela (Reference Valenzuela1976) in their solution to the Orr–Sommerfeld equation (see their figure
$3b$
for the background state). Viscosity, interestingly, can also generate a viscous H instability distinct from the inviscid H instability discussed so far; this was discovered and analysed in Parker et al. (Reference Parker, Caulfield and Kerswell2020) for smooth shear and density profiles of the tan hyperbolic form. These results in turn imply that, as
$\delta \rightarrow 1^{-}$
, for our inviscid linear stability results to be physically relevant and realisable, we need to consider background states with shear in both fluids. This is in contrast to (2.1), which is expected to be an accurate physical representation in the limit of
$\delta \rightarrow 0$
only.
We demonstrate below that the transition in the nature of instability observed with varying
$\delta$
for the background state of (2.1) also occurs for profiles with background shear in both fluids. Towards this, we first compute the self-similar velocity profile at the interface of two viscous, immiscible fluids with a flat density interface; this background state seems to have been first reported by Lock (Reference Lock1951). Figure 25 depicts the velocity profile as a function of the similarity coordinate
$\eta$
(see figure caption for a definition) obtained by solving the equations of Lock (Reference Lock1951) (see figure 25 caption for details). Expectedly, as
$\delta$
is increased (for a fixed viscosity ratio), the background shear in the lower fluid becomes comparable to that in the upper one. In the same figure, for comparison, the ‘double exponential’ approximation (expression (5.1a
); Young & Wolfe Reference Young and Wolfe2014) to the viscous, self-similar velocity profile of Lock (Reference Lock1951) are also provided. This double exponential profile is given by
\begin{align} & U(z) = \left \{ \begin{array}{ll} U_u(z)=U_{\infty } - \left (U_{\infty } - U_{0}\right )\exp \left (-\dfrac {z}{\Delta _u}\right ), & \qquad z \gt 0 ,\\ U_l(z)=U_0\exp \left (\dfrac {z}{\Delta _l}\right ), & \qquad z \lt 0, \end{array}\right . \end{align}
\begin{align} & \rho (z) = \left \{ \begin{array}{ll} \rho _u, & \qquad z \gt 0, \\ \rho _l, & \qquad z \lt 0, \end{array}\right.\\[10pt]\nonumber \end{align}
where
$U_0$
is the background interface velocity at
$z=0$
.
Vertical variation of the horizontal velocity profile with
$\delta$
, i.e.
${u_{u}}/{U_{\infty }}= f_{1}(\eta _{1}),\; {u_{l}}/{U_{\infty }}= f_{2}(\eta _{2})$
, for fixed viscosity ratio (
$\mu ^{u}/\mu ^{l}=0.02$
) corresponding to air–water and
$Bo \rightarrow \infty$
. Note that
$\eta =\eta _{1}$
and
$\eta _{2}$
for upper and lower fluids, respectively. The profiles (dashed) are obtained by numerically solving equation (
$9$
) (upper layer) and equation (
$15$
) (lower layer) in Lock (Reference Lock1951) governing
$f_1,f_2$
, subject to boundary conditions provided therein. We choose the far-field velocity in the lower fluid
$U_2=0$
as
$\eta _2 \rightarrow -\infty$
(Lock Reference Lock1951), employing the shooting method. The variable
$\eta _{1} \equiv y\sqrt {{U_{\infty }}/{\nu _u x}}$
with
$\eta _2 = \sqrt {( {\nu _u}/{\nu _l})}\eta _{1}$
. In each case, we fit a double exponential velocity profile of the form provided in expressions (5.1) with the (dimensional) coordinate
$z$
in the upper and lower fluids being replaced by (non-dimensional)
$\eta _{1}$
and
$\eta _{2}$
, respectively, and the dimensional variables
$\Delta _u$
and
$\Delta _l$
replaced by (non-dimensional)
$\mathcal{A}$
and
$\mathcal{B}$
, respectively. The fitted parameters are (a)
$U_0\approx 0.01386$
,
$\mathcal{A}=1.1777$
,
$\mathcal{B}=5.766$
; (b)
$U_0\approx 0.1685$
,
$\mathcal{A}=1.607$
,
$\mathcal{B}=3.659$
; (c)
$U_0\approx 0.2046$
,
$\mathcal{A}=3.408$
,
$\mathcal{B}=3.408$
. For all three panels and the double exponential profiles in these,
$ \textit{Fr}=4$
.
Expressions (5.1) generalise the earlier background state of (2.1); setting
$U_0 = 0$
in the former recovers the latter. A good qualitative match between the double exponential and the viscous profiles is apparent in figure 25. Note, in particular, that for the air–water density ratio
$\delta = 0.001 \ll 1$
, the shear in water is negligibly small (
$U_0 \approx 0$
, leftmost panel of figure 25) and the viscous as well as the double exponential profile are modelled well, by the single exponential profile of (2.1). This is not so at higher
$\delta \rightarrow 1^{-}$
(middle and right panels of figure 25) where the background shear in the lower fluid becomes substantial. Our choice of the double exponential profile as an approximation to the self-similar viscous profile (Lock Reference Lock1951) is partly motivated by the observation that the (inviscid) linear stability problem can also be analytically solved for this profile (Young & Wolfe Reference Young and Wolfe2014). This is in contrast to the profile of Lock (Reference Lock1951), for which only numerical solutions to the Rayleigh equation seem feasible. An additional and important advantage is that the change in the sign of
$D^2U(z)$
across
$z=0$
for expression (5.1a
) implies that the KH instability can be expected for this profile in the unstratified limit
$\delta =1$
, as a consequence of Rayleigh’s inflection point theorem.
Due to the good qualitative match seen between the self-similar, viscous profiles (Lock Reference Lock1951) and the double exponential profiles in figure 25, we proceed next with the intuitive expectation that linear, inviscid stability results, including transitions in the instability mechanism obtained for the background state (5.1), also apply to the viscous profiles of Lock (Reference Lock1951). We thus examine the inviscid, linear stability of the background state (5.1) as a function of the density ratio
$\delta$
.
The stability analysis of (5.1) is carried out in an exactly similar manner as § 2 earlier, the only difference being that both fluids now have shear in their background state. The Rayleigh equations in both fluids may be solved for this background state in terms of the Gauss hypergeometric function; we omit details as these can be found in Young & Wolfe (Reference Young and Wolfe2014). Note that two additional non-dimensional numbers appear for the background state of (5.1a
), (5.1b
), viz.
$\Delta _l/\Delta _u$
and
$U_s/U_{\infty }$
: both these ratios were zero in our analysis so far for the background state of (2.1).
The vertical variation of Reynolds stress and perturbation kinetic energy as a function of
$\delta$
are depicted in figures 26(a) and 26(b), respectively. Note that the parameters of the double exponential profile were chosen such that the rippling instability (Young & Wolfe Reference Young and Wolfe2014) does not manifest. In figure 26, one notes the same transition in the vertical variation of either metrics as was noted in figure 15 for the single exponential profile. Furthermore, results from nonlinear simulations initialised with eigenfunctions obtained from linear stability analysis of the double exponential profile of (5.1) are shown in figure 27. It thus becomes clear that our earlier results, notably the transition from Miles to the H and eventually to the KH instability, with increasing
$\delta$
remain applicable also to the double exponential background state of (5.1a
) and (5.1b
). The applicability of our results, spanning a range of density ratios from very low (
$\delta \ll 1$
) to close to unity (
$\delta \sim O_m(1)$
), is thus apparent.
(a) Vertical variation of Reynolds stress. (b) Perturbation kinetic energy for the double exponential background profile of expressions (5.1). Here
$ \textit{Fr}=4,\; Bo=87.813,\;\Delta _l/\Delta _u = 0.1,\; U_s/U_{\infty }= 0.1$
.
Evolution of the interface initialised with eigenmodes obtained from linear stability analysis of the background state in (5.1). Results are shown for (a)
$t/\tau =8.317$
, (b)
$t/\tau =3.0$
and (c)
$t/\tau =2.55$
, where
$\tau$
is the time taken to travel one wavelength. Here
$ \textit{Fr}=4,\; Bo=87.813,\;\Delta _l/\Delta _u = 0.212,\; U_s/U_{\infty }= 0.1$
.
6. Conclusions
In this study we have identified a novel transition in the nature of the fastest growing mode for the exponential profile with a sharp density interface. This transition has been studied at a low Froude number and high Bond number, and is observed with varying density ratio. At
$\delta \ll 1$
, the fastest growing mode has features typical of the Miles instability (Miles Reference Miles1957) with a near-discontinuous phase change in the eigenfunction across
$z_c$
; the extremely low value of
$c_i/c_r$
even for the fastest growing mode in this regime makes it behave nearly like a singular neutral mode (Maslowe & Spiteri Reference Maslowe and Spiteri2013), leading to energy extraction only below the critical location and upto the density interface, i.e. a subset of the shear layer. As the density ratio is increased beyond
$0.1$
, the value of
$c_i/c_r$
increases proportionally. For a growing mode, while the semi-circle theorem always ensures that
$U_u(z_c)=c_r(k)$
, the sharp jump in Reynolds stress at
$z_c$
gives way to a smooth variation through
$z_c$
, implying energy extraction throughout the shear layer. We rationalise this transition employing an analytical expression for the Reynolds stress, originally due to Lin (Reference Lin1954). The aforementioned smoothing is also shown to occur in the phase of the leading-order term bearing the logarithmic singularity in Tollmien’s inviscid, series solutions. At
$\delta =0.5$
, we show that the growth rate for the fastest mode compares well against that of the PL profile. This suggests that the transition from the H to KH instability seen in the PL profile with increasing
$\delta$
can also manifest for the exponential model, albeit modified by the presence of the critical location in the latter. We emphasise that, for the PL profile, due to a lack of profile curvature (except at the corners), the unstable eigenmodes are not singular even though they have a critical layer (see similar observations for stable eigenfunctions in Kadam et al. (Reference Kadam, Patibandla and Roy2023) for the case of a water layer at
$\delta =0$
with a linear shear profile). We have further validated that, for the asymmetric profile under consideration and with increasing
$\delta$
, the Reynolds stress and perturbation kinetic energy do not bear any signature of a potential transition from the H to KH instability that we anticipate. To establish this, we turn to nonlinear simulations. A pictorial summary of key results is shown in figure 28.
A pictorial summary of key results. The colour contours are informed by results obtained from theory and simulations. These have been generated synthetically, chosen so as to highlight the transition in instability mechanism with increasing
$\delta$
.
In addition to the excellent validation of linear predictions at least for a few wave periods, our simulations clearly show qualitative differences in the nonlinear regime between
$\delta =0.5$
and
$\delta =0.9$
. The former is associated with a finite-amplitude wave with a sheared cusp emitting spume-like droplets from its tip. In marked contrast, at
$\delta =0.9$
, nonlinear simulations recover the classic KH spirals. We show that these differences are consistent with the classification scheme proposed in Lawrence et al. (Reference Lawrence, Haigh and Zhu1998), distinguishing the H from KH instability for asymmetric profiles.
An important takeaway from our study is the observation that the Miles instability persists for
$\delta =0.01$
(at low
$\textit{Fr}$
and high
$Bo$
), i.e. even at ten times the air–water ratio. Given that the growth rate of this instability scales with density ratio, this plausibly makes it easier to experimentally measure growth rates of the Miles instability by exploring fluid combinations other than air–water. Our study also reports several comparisons against the PL profile in high- and low-
$\textit{Fr}$
regimes, establishing the qualitative role of background profile curvature. The role of viscosity requires careful analysis involving the Orr–Sommerfeld equation. This is underway and will be reported subsequently.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2026.11588.
Acknowledgements
R.D. acknowledges Professor Anubhab Roy for suggesting the Holmboe instability at the early stage of this study. Figure 19(c) was generated by Nikhil Yewale and is acknowledged. We thank the organisation of the International Research Network (IRN) under the auspice of CNRS for the workshop held at Univ. Paris, Saclay, France in June 2025 where this work was presented. We thank Professor Herbert E. Huppert for helpful comments. We thank the anonymous referees for several helpful remarks concerning interpretations in terms of wave interactions.
Funding
We gratefully acknowledge financial support from DST-SERB (Govt. of India) grants MTR/2019/001240, CRG/2020/003707, SPR/2021/000536, MoE-STARS/STARS-2/2023-0595 on various problems concerning nonlinear ocean waves, wave breaking and spray formation. S.R. is supported through IIT Bombay Seed Grant RD/0522-IRCCSH0-020. The Ph.D. tenure of A.K. is supported by the SERB grant (CRG/2020/003707) and is gratefully acknowledged. Computational time on the IIT Bombay High Performance Computing (HPC) facility is gratefully acknowledged.
Declaration of interests
The authors report no conflict of interest.
