On the generation of free-surface waves by instabilities in quadratic shear flows

TL;DR

This study explores how linear instabilities in quadratic shear flows generate free-surface waves, revealing the effects of flow curvature, viscosity, and forcing on stability and identifying five distinct instability mechanisms.

Contribution

It provides a comprehensive analysis of free-surface wave generation in quadratic shear flows, including new stability boundaries and the identification of five instability types.

Findings

Convex curvature suppresses long-wave instabilities.

Concave curvature suppresses short-wave instabilities.

Strong forcing increases long-wave instability at high Reynolds numbers.

Abstract

This paper investigates the generation of free-surface waves in a liquid layer driven by linear instabilities in Couette-Poiseuille (quadratic) shear flows. The base velocity profiles are characterized by a curvature parameter, and two-dimensional viscous and inviscid perturbations are analyzed across a wide parameter space of curvature, wavenumber, and Reynolds number, for fixed Froude and Bond numbers. In the inviscid limit, analytical solutions of the Rayleigh equation reveal that velocity profiles ranging from Nusselt to linear flows remain stable against the rippling instability, with long-wave growth occurring only under strong interfacial forcing, whereas weaker forcing produces well-defined stability boundaries. For the viscous problem, Orr-Sommerfeld computations and asymptotic analyses reveal that a slight convex curvature of the shear flow suppresses long-wave instabilities, while a slight concave curvature suppresses short-wave instabilities, so even small deviations from a linear profile produce qualitatively different behaviors. Furthermore, we observe that strongly forced long waves are more unstable at large $Re$ than the inviscid value they latch on to as $Re \to \infty$. Growth-rate maps highlight smooth transitions between long-wave and rippling modes and reveal an additional shear instability near the linear profile at high Reynolds numbers. Based on energy transfers and eigenfunction structures, five distinct instability types are identified: shear, rippling, long-wave interfacial, short-wave interfacial, and a composite mode that combines features of shear, rippling and long-wave interfacial instabilities at large Reynolds numbers.

Tables

Table 1. Table 1: Range of parameters reported in various experiments on air–water or air–water+glycerol two-phase flows. Key non-dimensional parameters are Reynolds number ( R ​ e = ρ 2 ​ U s ​ h 2 / μ 2 Re=\rho_{2}U_{s}h_{2}/\mu_{2} ), inverse squared Froude number ( G = g ​ h 2 / U s 2 G=gh_{2}/U_{s}^{2} ), and Bond number ( B ​ o = ρ 2 ​ g ​ h 2 2 / T Bo=\rho_{2}gh_{2}^{2}/T ). Here, U s U_{s} is the interfacial velocity, g g is the acceleration due to gravity, T T is the surface tension, ρ 2 \rho_{2} is the density, μ 2 \mu_{2} is the dynamic viscosity, and h 2 h_{2} is the depth of the bottom layer. References marked (*) estimate U s U_{s} as 2 % 2\% of free-stream velocity, while those marked (†) use surface tension values from Takamura et al. ( 2012 ) .

Paper	$U_s$ (mm/s)	$h$ (mm)	$Re$	$G$	$Bo$
Cohen & Hanratty (1965)†	40 – 220	1 – 7	$O\left(10\right)$ – $O\left({10}^3\right)$	$O\left({10}^{-1}\right)$ – $O\left(10\right)$	$\sim O\left(1\right)$
Craik (1966)	3 – 68	0.1 – 1.5	$O\left({10}^{-1}\right)-O\left({10}^2\right)$	$O\left(1\right)-O\left({10}^2\right)$	$O\left({10}^{-3}\right)-O\left({10}^{-1}\right)$
Hidy & Plate (1966)*	16 – 296	25 – 100	$O\left({10}^2\right)-O\left({10}^4\right)$	$O\left(1\right)-O\left({10}^3\right)$	$O\left({10}^2\right)-O\left({10}^3\right)$
Buckley & Veron (2016)*	10 – 330	700	$\sim O\left({10}^5\right)$	$O\left({10}^2\right)-O\left({10}^4\right)$	$\sim O\left({10}^6\right)$
Paquier et al. (2016)†*	60 – 168	35	$O\left(10\right)-O\left({10}^3\right)$	$O\left(10\right)-O\left({10}^2\right)$	$\sim O\left({10}^2\right)$
Liu et al. (2017)	44 – 387	0.195 – 1.09	$O\left(10\right)-O\left({10}^2\right)$	$O\left({10}^{-2}\right)-O\left(1\right)$	$O\left({10}^{-3}\right)-O\left({10}^{-1}\right)$

Equations

$$\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0,$$

$$\frac{\partial u}{\partial t}+U\frac{\partial u}{\partial x}+u\frac{dU}{dz}=-\frac{\partial p}{\partial x}+\frac{1}{Re}{\nabla }^{2}u,$$

$$\frac{\partial w}{\partial t}+U\frac{\partial w}{\partial x}=-\frac{\partial p}{\partial z}+\frac{1}{Re}{\nabla }^{2}w,$$

$$Re=\frac{\rho U_{s}h}{\mu },$$

$$u=0,w=0.$$

$$\frac{\partial \hat{\eta }}{\partial t}+\frac{\partial \hat{\eta }}{\partial x}=w,\hat{\eta }\frac{d^{2}U}{d{z}^2}+\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}=0,$$

$$p+\frac{2}{Re}\frac{\partial u}{\partial x}-G\left(\hat{\eta }-B{o}^{-1}\frac{{\partial }^2\hat{\eta }}{\partial x^2}\right)=0,$$

$$G=\frac{gh}{U_s^2},$$

$$Bo=\frac{\rho g{h}^2}{T}.$$

$$\{\psi ,p,\hat{\eta }\}=\{\varphi (z),f(z),\eta \}{e}^{ik(x-ct)},$$

$$(U-c)\left(\frac{d^2}{d{z}^2}-k^2\right)\varphi -U^{\prime \prime }\varphi =\frac{1}{ikRe}{\left(\frac{d^2}{d{z}^2}-k^2\right)}^2\varphi ,$$

$$\varphi (-1)=0,{\varphi }^{\prime }(-1)=0.$$

$$\eta =\frac{\varphi }{c-U},\eta U^{\prime \prime }+{\varphi }^{\prime \prime }+k^2\varphi =0,$$

$${\varphi }^{\prime \prime \prime }-3k^2{\varphi }^{\prime }+ikRe\left(U^{\prime }\varphi +(c-U){\varphi }^{\prime }\right)-ikReG\eta \left(1+\frac{k^2}{Bo}\right)=0.$$

$$U(z)=a{z}^2+(a+1)z+1,$$

$$(U(z)-c)\left({\varphi }^{\prime \prime }(z)-k^2\varphi (z)\right)-U^{\prime \prime }(z)\varphi (z)=0,$$

$$\varphi (-1)=0,$$

$${\Xi }_w(c-U{)}^2+U^{\prime }(0)(c-U)-G\left(1+k^{2}B{o}^{-1}\right)=0,\text{where}{\Xi }_w=\frac{{\varphi }^{\prime }(0)}{\varphi (0)}$$

$$\varphi (z)=e^{\kappa (\zeta (z)-{\zeta }_c)}({\zeta }_c^2-\zeta (z{)}^2)\left(C_1{\varphi }_1(z)+C_2{\varphi }_2(z)\right),$$

$${\varphi }_1(z,{\zeta }_c)=\frac{\text{Hc}\left[2+4\kappa {\zeta }_c,4\kappa {\zeta }_c,2,0,4\kappa {\zeta }_c,\frac{{\zeta }_c+\zeta (z)}{2{\zeta }_c}\right]}{{\zeta }_c-\zeta (z)},{\varphi }_2(z,{\zeta }_c)={\varphi }_1(z,-{\zeta }_c).$$

$$\zeta (z)=U^{\prime }(z)=2az+(a+1),{\zeta }_c=\zeta (z_c),\kappa =k/2a,$$

$${\Xi }_w(c-1{)}^2+(a+1)(c-1)-G(1+k^{2}B{o}^{-1})=0,$$

$${\Xi }_w=k+\frac{a+1}{1-c}+\frac{{\varphi }_1(-1){\varphi }_2^{\prime }(0)-{\varphi }_2(-1){\varphi }_1^{\prime }(0)}{{\varphi }_1(-1){\varphi }_2(0)-{\varphi }_2(-1){\varphi }_1(0)}.$$

$${\int }_{-1}^0\frac{1}{(U(z)-c{)}^2}dz-\frac{1}{G}=0.$$

$$G\left(\alpha (1+a(2c-1))+2ac(c-1)\log \left(\frac{2c-1+a+\alpha }{2c-1+a-\alpha }\right)\right)-c(c-1){\alpha }^3=0,$$

$${\left(\frac{p^{\prime }}{(U(z)-c{)}^2}\right)}^{\prime }-k^2\frac{p}{(U(z)-c{)}^2}=0,$$

$$k^2(U(0)-c{)}^{2}p(0)-G(1+B{o}^{-1}{k}^2)p^{\prime }(0)=0,p^{\prime }(-1)=0,$$

$$p_0(z)=\{\begin{array}{ll}1-{\left(\frac{\zeta (z)}{1+a}\right)}^5,& 0\ge z>z_c,\\ 1,& z_c>z\ge -1,\end{array}$$

$$\frac{p^{\prime }(0)}{(1-c{)}^2}=k^2{\int }_{-1}^0\frac{p}{(U-c{)}^2}dz.$$

$$-\frac{160a^3}{(a+1{)}^5}\sim -k^2\frac{i\pi }{2\sqrt{a}(c_1{)}^{3/2}},$$

Full text

**MSC Codes ** (Optional) Please enter your MSC Codes here.

1 Introduction

Generation of waves through linear instabilities of parallel two-phase shear flows is a canonical problem with wide ranging applications. Such instabilities underline phenomena as diverse as the generation of surface waves on the ocean (Sullivan & McWilliams, 2010), spray-generation in combustion applications (Lefebvre & McDonell, 2017), oil transportation (Joseph & Renardy, 2013) and coating technology (Weinstein & Ruschak, 2004). In geophysical settings, shear driven instabilities are recognized as the fundamental mechanism of wave generation for wind-driven waves on lakes and the ocean. These waves in turn interact and modulate large-scale winds, currents, and mixing in the ocean (Young & Van Vledder, 1993; Constantin & Ivanov, 2019; Dai et al., 2010).

Depending on the application, various flow configurations have been considered in the literature, such as two-phase Couette flows, two-phase Poiseuille flows, core-annular flows and free-surface flows (Govindarajan & Sahu (2014), see Figure 2). Multiple instabilities were found in different regions of the parameter space and they were broadly classified depending on the dominant energy transfer mechanism leading to a growth in perturbation kinetic energy (Boomkamp & Miesen, 1996), and the eigenfunction structure (Charru & Hinch, 2000). Among the various mechanisms proposed for wave generation on the ocean surface (Ayet & Chapron, 2022), the Miles instability (Miles, 1957) remains the most prominent and experimentally supported model (Plant, 1982). This inviscid instability occurs due to a transfer of energy, through wave Reynolds stress, from a critical layer ($z_{c}$, defined as the vertical location where the wave phase speed ($c_{r}$) matches the background flow velocity i.e $c_{r}=U(z_{c})$) in the air phase to the interface (Janssen, 2004). The asymptotic calculation of Miles (1957) and Miles (1959) shows that this energy transfer is proportional to the air-layer velocity profile curvature ($U^{\prime\prime}(z_{c})$) implying that Miles instability will not exist in a constant shear flow (such as in the viscous sub-layer). Evidence of a critical layer in the air-layer was provided in the experiments of Buckley & Veron (2016) and Carpenter et al. (2022). A similar mechanism of instability but with the critical layer in the water phase was studied by Stern & Adam (1973) and Morland et al. (1991), indicating wave generation due to shear currents. The presence of wind stress driven currents in the ocean has been documented both in classical near equatorial wind-drift theories (Stommel, 1959) and more recent exact solutions for azimuthal equatorial flows with a free surface (Constantin & Johnson, 2016). These studies highlight the variety of velocity profile curvatures (including flow reversals; see Constantin & Johnson, 2019) that can arise in the water layer and are crucial for the onset of instability. In particular, asymptotic analysis of Shrira (1993) and Bonfils et al. (2023) show that a negative velocity profile curvature in the water layer is necessary for wave growth. A combined two-phase inviscid problem comparing the growth rates of both instabilities was studied by Young & Wolfe (2014), who found that the latter instability is important at smaller wavelengths (naming it the “rippling instability”) and can attain higher growth rates compared to Miles instability. A combined two-phase viscous problem was also studied by Zeisel et al. (2008), who identified a second instability at high friction velocities and found, counter-intuitively, that viscosity enhances Miles instability growth rates. More recently, Kadam et al. (2023) showed that this enhancement is due to a fundamental shift in the instability energy source from Reynolds stress at the critical layer to tangential stress at the interface. They also showed that viscosity reduces rippling instability growth rates, although the underlying instability mechanism is altered. The present study is motivated by the non-trivial dependence of instability mechanisms on both viscosity and background flow curvature. To this end, we focus on quadratic (Couette–Poiseuille) velocity profiles, which not only reproduce the forward and backward bulges observed in field and laboratory flows, but also remain tractable for both inviscid and viscous stability analyses. This framework enables us to disentangle the respective roles of curvature and viscosity on free surface instabilities. Some of the well-known free-surface instabilities associated with different shear flows, as discussed above, are schematically illustrated in figure 1.

Viscosity-stratified instabilities were first identified by Yih (1967) for Couette–Poiseuille base-state flows, where the velocity profile is quadratic with coefficients depending on viscosity ratio ($m=\mu_{2}/\mu_{1}$), density ratio ($r=\rho_{2}/\rho_{1}$), and height ratio ($n=h_{2}/h_{1}$). Subsequent long-wave asymptotic analyses revealed instabilities at arbitrarily small Reynolds numbers ($Re=\rho_{2}U_{s}h_{2}/\mu_{2}$). A similar instability had been noted earlier by Kapitza (1948) and studied rigorously by (Benjamin, 1957; Yih, 1963) in the context of falling films, called the Kapitza instability. The former is viscosity-induced while the latter is gravity-induced, with combined effects also studied in the literature (see Boomkamp & Miesen (1996)). Long- and short-wave instabilities were later detailed by Smith (1990), Charru & Hinch (2000), and Kalliadasis et al. (2011), with Hooper & Boyd (1983) further identifying a short-wave instability in constant-shear flows due to tangential stress at the interface, a mechanism later described by Hinch (1984). In both cases, growth rates are proportional to $|m-1|$ for $r=1$, while density stratification ($r\neq 1$) can stabilize the flow under certain conditions (Hooper & Boyd, 1983). These instabilities, dominated by interfacial mechanisms, occur in a range of natural and technological contexts from magma and glaciers to oil transport and biological flows (Govindarajan & Sahu, 2014).

Apart from inviscid and interfacial instabilities, shear-mode instabilities have also been identified. For a single-phase free-surface constant shear flow, Miles (1960) derived the critical Reynolds number ($Re_{c}\approx 203$), showing that a viscous instability occurs at high $Re$ for modes which are inviscidly neutral with a critical layer, consistent with Heisenberg’s criterion (Lin, 1946). Later, Smith & Davis (1982) considered surface-tension-gradient-driven shear flows and found much lower $Re_{c}$, while Hooper & Boyd (1987) used asymptotic theory to study shear-mode instabilities in two-phase flows. Further works such as Yiantsios & Higgins (1988), Hooper (1989), and Miesen & Boersma (1995) established that shear modes can coexist with interfacial modes while mode coalescence at high $Re$ is reported by Timoshin (1997), Özgen et al. (1998), Timoshin & Hooper (2000), and Kaffel & Riaz (2015). A recent review by Mohammadi & Smits (2016) emphasized the importance of detailed studies of shear–interfacial mode interactions and their role across parameter regimes.

In addition to theoretical studies, numerous experiments have investigated two-phase instabilities across configurations (see Govindarajan & Sahu, 2014, for details). Viscosity-stratified instabilities in oil–water systems have been studied in core-annular (Bai et al., 1992; Kouris & Tsamopoulos, 2001), circular Couette (Sangalli et al., 1995; Barthelet et al., 1995), and planar parallel flows (Charles & Lilleleht, 1965; Kao & Park, 1972). For air–water like systems, several experiments in planar channel flows are summarized in Table 1. These cover a broad range of $Re$, $G$, and $Bo$ and report distinct instabilities such as “fast” and “slow” waves (Cohen & Hanratty, 1965; Craik, 1966), viscous modifications of Miles instability (Hidy & Plate, 1966; Buckley & Veron, 2016), transitions between wrinkle and wave regimes (Paquier et al., 2016), and glaze-ice-related film instabilities (Liu et al., 2017). Interpretations vary across studies, with energy analyses attributing most of these waves to interfacial instabilities (Boomkamp & Miesen, 1996), though Reynolds stress contributions are sometimes non-negligible (Náraigh et al., 2011). These works collectively demonstrate the complex interplay of viscosity, shear, and interfacial mechanisms in air–water instabilities.

The experiments and theory above emphasize that the two-phase stability problem is governed by six non-dimensional parameters—$Re$, $G$, $Bo$, $m$, $r$, and $n$—making comprehensive characterization difficult. Fixing the phases (e.g., air–water) reduces this to three: $Re$, $G$, and $Bo$. Further simplification arises from the free-surface approximation, where the air layer is treated as passive. This approximation neglects air-side instabilities like Miles’ instability but remains justified due to the large density ratio ($\rho_{2}/\rho_{1}\gg 1$) and small viscosity ratio ($\nu_{2}/\nu_{1}\ll 1$). Young & Wolfe (2014) have shown that this approximation preserves key features of the rippling instability. In the present work, we adopt this framework and consider a family of quadratic Couette–Poiseuille velocity profiles parameterized by curvature $a$, allowing a systematic exploration of profile curvature effects without explicitly specifying $m$, $r$, and $n$.

To the best of our knowledge, a comprehensive study of how viscosity and velocity profile curvature affect rippling instability and its interaction with other modes has not been carried out. The present work addresses this gap. Section 2 describes the governing equations and base-state velocity profiles. Inviscid stability results, including analytical solutions, asymptotic expressions, stability boundaries, and complete parameter-space behavior, are discussed in §3. Viscous stability analyses, combining asymptotic and numerical approaches, are presented in §4. Section 5 characterizes instability families using energy and eigenfunction structures, and conclusions are summarized in §6.

2 Problem formulation

We consider a two-dimensional, incompressible, parallel flow in a water layer, represented by a family of quadratic base-state velocity profiles, i.e. Couette–Poiseuille flows. The domain is unbounded in the horizontal direction and bounded vertically by a free surface at the top and a rigid wall at the bottom, located at a dimensional depth $h$ below the interface. The base-state shear flow velocity at the free surface, $U_{s}$, is chosen as the characteristic velocity scale, while $h$ is chosen as the characteristic length scale, and $h/U_{s}$ as the characteristic time scale, to nondimensionalize all variables. Henceforth, all quantities in the analysis are nondimensional unless otherwise specified.

Assuming small-amplitude perturbations, the governing set of nondimensional linearized equations for the perturbation fields can be written as

[TABLE]

[TABLE]

[TABLE]

where $x$ and $z$ are the horizontal and vertical coordinates, respectively, $t$ denotes time, $u(x,z,t)$ and $w(x,z,t)$ are the perturbation velocity components, $p(x,z,t)$ is the perturbation pressure, and $U(z)$ is the base-state velocity profile. The operator $\nabla^{2}(\cdot)$ denotes the two-dimensional Laplacian in the $(x,z)$ coordinates. The Reynolds number is defined as

[TABLE]

where $\rho$ and $\mu$ are the density and dynamic viscosity of water, respectively.

At the rigid wall ($z=-1$), the no-slip and no-penetration boundary conditions are

[TABLE]

At the free surface ($z=0$), the kinematic, tangential stress, and normal stress boundary conditions, linearized about the base state, take the form

[TABLE]

[TABLE]

where $z=\hat{\eta}(x,t)$ is the free surface displacement, $G$ is the inverse squared Froude number,

[TABLE]

and $Bo$ is the Bond number,

[TABLE]

Here, $g$ is the gravitational acceleration and $T$ is the air–water surface tension. Condition (5b) expresses that the total tangential stress vanishes at the free surface, while (6) represents the balance of perturbation normal stresses with surface tension.

Introducing a perturbation streamfunction $\psi(x,z,t)$ such that $u=\partial\psi/\partial z$ and $w=-\partial\psi/\partial x$, and adopting a normal-mode decomposition

[TABLE]

equations (1)–(3) can be combined into the Orr–Sommerfeld equation (Drazin & Reid, 2004),

[TABLE]

where $k$ is the nondimensional wavenumber and $c$ is the complex phase speed of the perturbation.

At the rigid wall ($z=-1$), the boundary conditions reduce to

[TABLE]

The linearized conditions at the free surface ($z=0$) are

[TABLE]

[TABLE]

Here, primes denote differentiation with respect to $z$. For a given $k$ and specified nondimensional parameters ($Re,G,Bo$), the Orr–Sommerfeld equation (7), subject to the boundary conditions (8)–(10), yields the eigenvalue $c(k)=c_{r}(k)+ic_{i}(k)$ and eigenfunction $\phi(z;k)$. The real part $c_{r}$ represents the phase speed, while $kc_{i}$ denotes the temporal growth rate. Modes with $c_{i}>0$ are unstable, whereas those with $c_{i}<0$ are stable.

The present study focuses on the limit of strong interfacial forcing ($G=0$), where gravitational effects are negligible compared to shear. In this limit, for any finite $Bo$, equation (10) shows that surface tension also does not contribute. This asymptotic regime has been examined in several previous works (Miles, 1960; Smith & Davis, 1982; Miesen & Boersma, 1995) and is especially relevant for thin-film flows (see Table 1). For completeness, cases with nonzero $G$ are also considered; in these, gravity is finite but surface tension is assumed small, corresponding to a large Bond number $Bo\approx 10^{5}$ throughout the study.

2.1 Base-state profiles

To investigate the effect of velocity profile curvature on stability characteristics, we consider a one-parameter family of Couette–Poiseuille base-state flow profiles defined as

[TABLE]

where $a$ is a curvature parameter that may take values in the range $-\infty<a<\infty$. Any quadratic velocity distribution can be recast into this nondimensional form, making it a convenient canonical representation. The parameter $a$ thus provides a systematic way to explore the influence of curvature on instability.

Several notable flow configurations correspond to special values of $a$. For $a=0$, equation (11) reduces to a linear velocity profile, characteristic of a simple shear flow. The case $a=-1$ yields the Nusselt profile, which arises in the classical problem of a liquid film flowing down an inclined plane under gravity. For $a=3$, the profile exhibits a flow reversal near the free surface, a phenomenon observed in laboratory experiments on wind-driven water waves (Hidy & Plate, 1966; Paquier et al., 2015, 2016). Interestingly, similar flow reversals occur in large-scale geophysical settings, such as equatorial oceanic currents, though with more complex vertical structures (Stommel, 1959; Constantin & Johnson, 2016). In the asymptotic limits, $a\to-\infty$ yields a rescaled forward Poiseuille-type flow, whereas $a\to+\infty$ corresponds to a backward Poiseuille-type profile.

The parameter space of $a$ can be naturally divided into four distinct regimes, each corresponding to a qualitatively different profile shape: (i) backward bulging ($1\leq a<\infty$), which includes the flow-reversal case; (ii) monotonic convex ($0<a<1$); (iii) monotonic concave ($-1<a\leq 0$), which contains the linear shear profile; and (iv) forward bulging ($-\infty<a\leq-1$), which includes the Nusselt profile. Representative examples of these regimes are illustrated in figure 2.

3 Inviscid linear stability analysis

In the inviscid limit ($\Rey\to\infty$), the Orr–Sommerfeld equation (7) reduces to the Rayleigh equation (Drazin & Reid, 2004),

[TABLE]

which governs the eigenfunction $\phi(z)$ for a given wavenumber $k$ and complex phase speed $c$. At the rigid bottom wall ($z=-1$), only the no-penetration condition remains,

[TABLE]

while the number of independent conditions at the free surface $z=0$ reduces from three to two. Specifically, the kinematic and normal-stress conditions can be combined into a single dynamic constraint (Young & Wolfe, 2014),

[TABLE]

3.1 Analytical solution

For the family of base-state velocity profiles considered, the Rayleigh equation (12) can be solved analytically in closed form, generalizing the special case $a=3$ studied in Kadam et al. (2023). After a sequence of transformations, the solution to (12) for arbitrary $a$ may be expressed in terms of confluent Heun functions (denoted by Hc) as

[TABLE]

where

[TABLE]

[TABLE]

with $C_{1}$ and $C_{2}$ as integration constants, and $z_{c}$ denoting the critical layer depth obtained by solving $U(z_{c})-c=0$. Additional details on confluent Heun functions and their alternate representations are provided in Appendix A. Here $\zeta(z)$ is the negative of the base-state vorticity and $\zeta_{c}$ its value at the critical level.

The dispersion relation follows by enforcing boundary conditions (13) and (14), giving

[TABLE]

where

[TABLE]

For prescribed $(a,k,c,Bo)$, equation (16) yields $G$ explicitly. Conversely, for fixed $(a,G,Bo,k)$, a root-finding procedure is required to determine $c$ and the corresponding growth rate $kc_{i}$. Thus, despite its compact analytical form, the dispersion relation cannot be analyzed directly without numerical evaluation. In the sections that follow, we examine its implications: Section 3.2 addresses long-wave asymptotics for the growth rate, while Section 3.3 derives explicit expressions for the stability boundary in different regimes of $a$.

3.2 Long wave asymptotics ($k\ll 1$)

In order to obtain physical insight without resorting to a numerical evaluation of the complete dispersion relation, we perform asymptotics in the long wave ($k\ll 1$) limit. The leading-order behaviour of long waves can be obtained from the Burns integral relation (Burns, 1953), which connects the wave phase speed to the Froude number in the limit $k\to 0$:

[TABLE]

For the velocity profiles in (11), this reduces to

[TABLE]

with $\alpha^{2}=(a-1)^{2}+4ac$.

For $G>0$ and every value of curvature ($a$), equation (19) yields two neutrally stable modes ($c_{i}=0$): a prograde branch with speed exceeding the maximum background velocity, and a retrograde branch with speed below the minimum velocity (figure 3a). Thus, long waves remain neutrally stable for any $G>0$. In the limit of $G=0$, however, the phase speeds $c_{r}$ for all the curvatures, of the prograde and retrograde modes attain the maximum and minimum values of the base state velocity, respectively (black curves in figure 3a). With the maximum given by \big{(}-(a-1)^{2}/4a for $a\leq-1$ and 1 for a>-1\big{)} and minimum \big{(}0 for $a\leq 1$ and $-(a-1)^{2}/4a$ for a>1\big{)}. This coincidence with the background velocity indicates the emergence of a critical layer and thereby the possibility of a rippling-type instability at small but non-zero $k$ (Young & Wolfe, 2014). Notably, when $a<-1$, the prograde mode itself can destabilize due to the forward bulge in the profile, consistent with the observations of Bonfils et al. (2023).

Figures 3b,c illustrate this behaviour at finite $k$ for $G=0$ and representative curvature parameters $a=-3$ (blue), $3$ (red), and $0.3$ (green). The real part $c_{r}$ (panel b), computed from the dispersion relation (16), begins at the extremal velocity and approaches the surface speed $U_{s}$ as $k$ increases, confirming the existence of a critical layer for all $k$. The imaginary part $c_{i}$ (panel c) demonstrates finite growth rates across the $k$ range, with maxima occurring at $O(1)$ wavenumbers, a typical feature of rippling instability (Young & Wolfe, 2014). The logarithmic axis for $k$ highlights the distinct small- and large-$k$ asymptotic behaviour.

To capture the small-$k$ growth rates for $G=0$, we adopt the perturbative approach of Renardy & Renardy (2013), using the leading-order pressure perturbation solution as input to the pressure pertubation equation at the next order. The governing equation is

[TABLE]

subject to

[TABLE]

which reduce to $p(0)=0$ and $p^{\prime}(-1)=0$ for $G=0$.

Case I ($|a|>1$). Substituting the phase speeds from (19) into (20), the $k=0$ solution reads

[TABLE]

with the critical layer at $z_{c}=-(a+1)/2a$ where $c=U_{\min/\max}$. This piecewise form reflects the neutral continuous spectrum. For finite $k$, integration of (20) yields the relation

[TABLE]

Using $p(z)\approx p_{0}(z)$ and setting $c=U_{\min/\max}+c_{1}$ with $|c_{1}|\ll 1$ gives

[TABLE]

leading to

[TABLE]

The non-integer power law ($c_{i}\sim k^{4/3}$) reflects the singular role of the critical layer. This generalises earlier results of Kadam et al. (2023) and recovers Renardy & Renardy (2013) in the $|a|\to\infty$ limit. 2. 2.

Case II ($|a|<1$). Here the $k=0$ solutions correspond to $c=1$ (prograde) and $c=0$ (retrograde). Since the Doppler shift prevents the prograde mode phase speed from matching with the background velocity profile, only the retrograde branch can destabilize. Following Renardy & Renardy (2013), we obtain a neutral solution at $c=0$, leading to the asymptotic relation

[TABLE]

with $\bar{k}=(1-a)k/2a$, $\bar{c}=4ac/(1-a)^{2}$ and

[TABLE]

This gives $c_{i}\sim k^{4}$ for $0<a<1$, while $a<0$ yields stability. Figure 3c (black stars) confirms the agreement with the full solution.

For large $k$, the growth rate scales as $c_{i}\sim k^{-2}$ in the $G=0$ limit, consistent with Kaffel & Renardy (2011). More generally, studies by Shrira (1993) and Bonfils et al. (2023) show exponential decay for $G>0$, with the instability mechanism tied to curvature at the critical layer. A complementary small-curvature asymptotic analysis is given in Appendix B.

Thus, long waves are neutrally stable for any non-zero $G$, while for $G=0$ the presence of a critical layer leads to the rippling instability with growth rates that depend on curvature $a$. In the next section, we derive explicit expressions for the stability boundary.

3.3 Stability boundary

A closed-form analytical solution, as outlined in §3.1, enables the numerical construction of the full stability contour. However, such an approach is computationally intensive and does not provide an immediate answer to the simpler question of whether instability arises for a given set of parameters. In this regard, a stability boundary that demarcates the stable and unstable regions of the parameter space serves as a valuable diagnostic tool. We construct this boundary by invoking the extension of Howard’s semicircle theorem to free-surface shear flows, due to Yih (1972), which restricts the phase speed $c_{r}$ of any unstable mode to lie between the minimum and maximum values of the base-state velocity profile. This property has been widely used (Morland et al., 1991; Engevik, 2000; Renardy & Renardy, 2013; Young & Wolfe, 2014; Kadam et al., 2023) to identify neutral modes at $c=U_{min}$ or $c=U_{max}$, which define the boundary. Here we apply this approach to obtain explicit expressions for the stability boundary in terms of $G$, $k$, and $a$. Substituting the extremal values of the velocity profile into the governing equations provides conditions for neutral stability.

Case I ($|a|>1$): In this regime, where the flow is either forward or backward bulging, the velocity extrema are $U_{max/min}=-(a-1)^{2}/4a$, occurring at $z_{c}=-(a+1)/2a$. Setting $c=U_{min/max}$ in (12) yields the neutral eigenfunction

[TABLE]

with $\zeta(z)$ and $\kappa$ defined in (15). Enforcing the surface condition at $z=0$ gives the stability boundary

[TABLE] 2. 2.

Case II ($0<a<1$): Here only the retrograde branch can destabilize. Taking $c=U_{min}=0$ does not yield a simpler eigenfunction, but substitution into the dispersion relation (16) gives

[TABLE]

where $\textrm{H}_{c}^{\prime}$ denotes the derivative of the confluent Heun function. 3. 3.

Case III ($-1<a<0$): In this range, the prograde mode remains stable. Substituting $c=1$ into the boundary condition (14) simplifies to

[TABLE]

demonstrating that no neutral mode exists for finite $G$.

In summary, (29)–(30) can be solved to obtain the minimum unstable wavenumber $k_{min}$ as a function of $G$ and $a$. The resulting contour map is shown in figure 4a. For $G=0$, $k_{min}=0$ across all $a$, consistent with the long-wave analysis. For $G>0$ and $a>0$, $k_{min}$ generally grows with $G$ but varies non-monotonically with $a$, peaking near $a\approx 1.41$. For $a<0$, $k_{min}$ increases rapidly with $G$, reaching $O(10^{2})$ at $G=1$. Since such large cutoffs are beyond our focus, we shade $k_{min}>3$ uniformly in cyan. The white band at $0<a<1$ reflects the absence of instability, in agreement with (31). For finite $G$, the same analysis can also provide the short-wave cutoff $k_{max}$, while for $G=0$ no cutoffs exist, as seen earlier in figure 3c.

3.4 Complete solution of the dispersion relation

As outlined in § 3.1 and Appendix A, the dispersion relation (16) can be solved numerically across a wide parameter space using a root-finding routine in Mathematica (Wolfram Research, 2024), yielding the complex phase speed. Figures 3b-c show $c_{r}$ and $c_{i}$ versus $k$ for $G=0$ and $a=-3$ (blue), $3$ (red), and $0.3$ (green). Corresponding results for $G=0.02$ and $a=\pm 3$ (blue), $\pm 5$ (red) are shown in figures 4b-c. As discussed earlier, neutral modes meet the unstable branch at the long-wave cutoff $k_{min}$. At $k=k_{min}$, $c_{r}$ crosses the velocity extrema $U_{min/max}$ (gray dashed lines in figure 4b). For $G=0$, $k_{min}=0$, so $c_{r}$ coincides with $U_{min/max}$ at onset. For $k>k_{min}$, modes are unstable, and both figure sets (3, 4) show the development of a critical layer. Unlike the $G=0$ case (fig. 3b), the critical layer does not approach the surface asymptotically when $G>0$. At large $k$, $c_{r}$ departs from the surface velocity, recrossing $U_{min/max}$ and recovering neutral stability. This short-wave cutoff is due to surface tension. It is not displayed in figure 4 since the cutoff occurs at very large $k$, obscuring the long-wave features of interest. In figure 4b, the $U_{max}$ line (gray) moves downward with $a$, eventually coinciding with $c_{r}=1$ at $a=-1$. The result is a larger long-wave cutoff, smaller short-wave cutoff, and a shift in the most unstable mode.

Figures 5a-c show growth-rate contours ($kc_{i}$) in the $(a,k)$ plane for $G=0$, $0.02$, and $0.5$. For $a<0$, only the prograde branch is unstable, whereas for $a>0$ the retrograde branch is unstable. As anticipated from (31), profiles with $-1<a<0$ are always stable. For $G=0$, instability exists for all $k$ whenever $a<-1$ or $a>0$. Any nonzero $G$ produces a finite region of stability, which grows with increasing $G$. Thus, quadratic velocity profiles with $a<-1$ or $a>0$ are always unstable in the inviscid limit, and gravity and surface tension act together to stabilize disturbances when $G\neq 0$. The black curves in figures 5a-c show the stability boundaries derived in § 3.3, in good agreement with the growth-rate contours. The gray dashed lines mark the locus of maximum growth rate in the $(a,k)$ plane. For retrograde modes ($a>0$), the maximum growth rate occurs at wavenumbers $O(1)$. For prograde modes ($a<-1$), the most unstable wavenumber increases as $a$ increases toward $-1$. This reflects the diminishing forward bulge of the velocity profile, which vanishes at $a=-1$. Away from the stability boundary in figures 5a-c, growth rates increase with $a$ (see also fig. 3.3c), consistent with the dependence of instability strength on the velocity curvature at the critical layer (Shrira, 1993; Bonfils et al., 2023).

Figure 5d illustrates the streamline pattern overlaid on normalized vorticity contours ($\Omega/\Omega_{max}$) for $G=0$, $a=-3$, and $k=6.31$ (black star in fig. 5a). Two critical layers are visible (black dashed lines), since $c_{r}=U(z)$ is satisfied at two distinct depths for forward-bulging profiles ($a<-1$). Around these layers, the streamfunction tilts and vorticity contours distort strongly, consistent with the negative Reynolds stress correlation $\overline{uw}<0$ that drives the instability (Young & Wolfe, 2014).

To summarize, the inviscid stability analysis demonstrates that for $G=0$, $c_{i}\sim k^{4/3}$ when $|a|>1$ and $c_{i}\sim k^{4}$ when $|a|<1$, consistent with asymptotic results. These scalings were confirmed using solutions expressed in terms of confluent Heun functions. All quadratic profiles are unstable except those with $-1<a<0$. With increasing $G$, long-wave cutoffs emerge, stabilizing low-wavenumber disturbances. Instabilities propagate faster than the surface velocity when $a<-1$ and slower when $a>0$. In both cases, the peak growth occurs at wavenumbers $O(1)$.

The inviscid framework adopted here is the same as in earlier studies of wind-generated surface waves (Miles, 1957; Stern & Adam, 1973; Young & Wolfe, 2014), justified by the typically large Reynolds numbers of geophysical flows. In contrast, laboratory experiments on thin layers (Paquier et al., 2015, 2016) involve Reynolds numbers of $O(10^{3})$, making viscous effects non-negligible. While viscous stability analyses are common in engineering, they have rarely addressed the rippling instability directly (see Table 1 of Young & Wolfe, 2014, for a summary). In the following section, we extend the analysis to finite Reynolds numbers, to ask: How does viscosity influence the instability? Does it merely damp the inviscid mode, or does it qualitatively alter the dynamics? And how does the inviscid rippling instability compare with other viscous modes?

4 Viscous linear stability analysis

The influence of viscosity on linear stability is well established in single-phase parallel shear flows, where it generally acts as a dissipative agent (Drazin & Reid, 2004). In air–water two-phase systems, however, its role is more subtle. As discussed in the introduction, Zeisel et al. (2008) and Kadam et al. (2023) showed that viscosity can amplify Miles-type instability growth at short wavelengths, while Kadam et al. (2023) also reported viscous damping of rippling instability under conditions corresponding to the laboratory experiments of Paquier et al. (2015, 2016). Thus, it is not obvious in which regions of parameter space viscosity suppresses, and in which it enhances, the rippling instability. The picture is further complicated by the fact that viscosity stratification alone, as in air over water, can destabilize the interface even at low Reynolds numbers (see Yih, 1967). Moreover, Miles (1960) showed that a linearly sheared liquid film with a free surface is unstable to viscous perturbations, with growth rates that decay as $Re^{-1/2}$ and vanish in the inviscid limit. This is consistent with the stability of the linear profile in §3, but emphasizes that viscosity is not always simply stabilizing: depending on the configuration, it can alter the character of existing instabilities or give rise to new ones.

In this section, we present the viscous linear stability results in the $(k,Re)$ parameter plane for different values of $a$ and $G$. The Orr–Sommerfeld equation (7) with the base velocity profiles (11) does not admit closed-form solutions. Accordingly, asymptotic results for the long- and short-wave limits are first derived in §4.1. Numerical solutions of the eigenvalue problem (7)–(10) are then presented in §4.2, followed by a synthesis of the results into a complete picture of growth-rate behavior across the parameter space in §4.3.

4.1 Asymptotic evaluation of growth rates

4.1.1 Long wave asymptotics ($k\ll 1$)

In the long-wave limit, following Yih (1967), the eigenfunction $\phi$ and eigenvalue $c$ may be expanded in a regular perturbation series in $k$,

[TABLE]

Substituting these expansions into (7)–(10) and solving the resulting equations order by order in $k$, as outlined in Appendix C (see also Yih, 1967), yields explicit expressions for the first three terms in the eigenvalue,

[TABLE]

[TABLE]

The leading-order and second-order contributions, $c^{(0)}_{k\ll 1}$ and $c^{(2)}_{k\ll 1}$, are purely real and thus do not affect the growth rate. The imaginary part arises solely from the $O(k)$ term, which controls the instability. In the limit of a strongly forced interface or in the absence of gravity ($G=0$), the imaginary part of the eigenvalue reduces to

[TABLE]

Equation (34) shows that long-wave (Yih-type) instability occurs when $a<0$ or $a>1$, whereas velocity profiles between the linear case ($a=0$) and the half-parabolic case ($a=1$) remain stable, in agreement with figure 15 of Smith & Davis (1982). Moreover, $c_{i}$ scales linearly with $kRe$, the natural small parameter in this expansion (Yih, 1967). Finally, from $c^{(0)}_{k\ll 1}$ in (33), the leading order phase speed exceeds the maximum velocity in the domain for $a<0$ and falls below the minimum velocity for $a>1$, confirming the asymmetry of propagation depending on profile curvature.

4.1.2 Short wave asymptotics ($k\gg 1$)

A parallel analysis to the long-wave limit can be carried out for short waves. This approach was first developed by Hooper & Boyd (1983) for two-layer linear–linear profiles and later generalized by Yiantsios & Higgins (1988) to channel flow. In this regime the small parameter is taken as $\epsilon=k^{-1}$, and the vertical coordinate is rescaled as $\tilde{z}=kz$. The eigenfunction and eigenvalue are then expanded in regular series of $\epsilon$,

[TABLE]

Substitution of these expansions into the system (7)–(10) and collection of terms at successive orders in $\epsilon$ gives an iterative hierarchy, the details of which are outlined in Appendix D. The first four contributions to the eigenvalue are,

[TABLE]

[TABLE]

Here $Ca=Bo/(GRe)$ is the Capillary number.

From (4.1.2) it is evident that the imaginary parts up to $O(\epsilon^{3})$ are negative, and therefore act to stabilize the interface. The first destabilizing contribution appears at $O(\epsilon^{4})$. In the absence of gravity ($G=0$), the imaginary part of $c$ simplifies to

[TABLE]

Thus all lower-order imaginary terms vanish at $G=0$, leaving only the destabilizing $O(\epsilon^{4})$ contribution. This implies that short-wave instability is more widespread in parameter space when gravity is absent, and that nonzero $G$ suppresses it through the stabilizing terms in (4.1.2). According to (37), instability arises for $a>0$ and $a<-1$, while velocity profiles between the linear case ($a=0$) and the Nusselt profile ($a=-1$) remain stable to short-wave disturbances in the $G=0$ limit.

4.2 Numerical solution

Although the asymptotic solutions above reveal the existence of instabilities, their validity is confined to the limiting cases of small and large wavenumbers. To analyze the behavior at intermediate parameters, a numerical approach is required, since (7) admits no closed-form solution. We therefore employ a Chebyshev spectral collocation method following Boomkamp et al. (1997) and Trefethen (2000) to solve the eigenvalue problem (7)–(10).

4.2.1 For small curvature parameters ($|a|\ll 1$)

As established in §4.1.1 and §4.1.2, long-wave instabilities are absent for $-1<a<0$ and short-wave instabilities are absent for $0<a<1$ in the $G=0$ limit. A direct implication is that a small perturbation of the curvature parameter $a$ away from zero (corresponding to the linear profile) results in an abrupt change in both the type of instability and the range of unstable parameters. Slightly convex profiles ($a>0$) develop short-wave instabilities, while slightly concave profiles ($a<0$) develop long-wave instabilities.

This behavior is illustrated in figure 6, which shows $c_{i}$ as a function of $k$ for $a=-0.05$ and $a=0.05$ at $Re=10$. The solid lines denote the numerical solutions, while the dashed lines correspond to the asymptotic predictions (34) and (37). The asymptotic scalings $c_{i}\sim k$ for the long-wave instability and $c_{i}\sim k^{-4}$ for the short-wave instability agree well with the numerical results. Notably, for small $|a|$, a sharp transition is observed from long-wave (small negative $a$) to short-wave (small positive $a$) instability at a critical wavenumber $k_{\text{cr}}$, which remains nearly constant across a range of Reynolds numbers. As $|a|$ increases, however, the transition becomes more gradual and a gap opens between the long- and short-wave cutoffs, which broadens progressively with curvature.

4.2.2 Transition from viscous to inviscid regimes

To illustrate how viscous instabilities (§4.1) connect to the inviscid limit, figure 7 shows the variation of $c_{i}$ with $Re$ for four representative velocity profiles: $a=3$, [math], $-1$, and $-3$, at fixed wavenumber $k=1.5$ in the $G=0$ limit. Two distinct behaviors emerge. For $a=3$ (blue) and $a=-3$ (purple), $c_{i}$ increases with $Re$ and asymptotically approaches the value of $c_{i}$ in the inviscid limit (black dashed lines). By contrast, for $a=0$ (red) and $a=-1$ (yellow), $c_{i}$ decays with increasing $Re$, indicating viscous instabilities that vanish in the inviscid limit. As noted in §3 and §4.1, the linear profile ($a=0$) is stable to long-wave, short-wave, and rippling instabilities; the observed mode here corresponds to the viscous instability identified by Miles (1960), which decays as $Re^{-1/2}$, consistent with the numerical results. Similarly, the Nusselt profile ($a=-1$) is inviscidly stable, and the observed instability corresponds to the viscous continuation of the long-wave mode. In contrast, the profiles $a=\pm 3$ are unstable to all viscous mechanisms, though their small-$Re$ behavior differs. The matching slopes of the $a=-1$ and $a=-3$ curves at low $Re$ suggest a long-wave character for $a=-3$, followed by a transition to rippling instability as $Re$ increases.

The case $a=3$ is particularly noteworthy, as it corresponds to velocity profiles observed in the experiments of Paquier et al. (2015, 2016). Figures 8a–c present its stability properties at $G=0$. Figure 8a shows $c_{i}$ versus $Re$ for $k=0.1,\,0.5,$ and $1.5$. While the $k=1.5$ case (yellow) is stable at low $Re$, both $k=0.1$ (blue) and $k=0.5$ (red) exhibit finite-$Re$ growth followed by a transition to rippling instability (dotted lines) at large $Re$. The non-monotonic variation of $c_{i}$ with $Re$ becomes more pronounced as the wavelength increases, demonstrating that viscosity can amplify rather than suppress disturbances—by more than an order of magnitude in some cases. Figure 8b shows $c_{i}$ as a function of $k$ for $Re=1$, $10^{5}$, and the inviscid limit. The $Re=1$ curve (red) displays asymptotic behavior at small and large $k$, consistent with the expressions of §4.1, and also reveals an interval at $O(1)$ wavenumbers where the instability is absent, in line with §4.2.1. At $Re=10^{5}$ (blue), the spectrum overlaps with the inviscid prediction (yellow dashed) for $10^{-1}\lesssim k\lesssim 10$, while deviations at smaller $k$ reflect the delayed approach to the inviscid asymptote seen in figure 8a. This behavior is clarified in figure 8c, where $c_{i}$ is plotted against $kRe$ for several long-wave cases. The collapse of the curves onto a single master trend reveals two regimes: at small $kRe$, $c_{i}$ grows linearly with $kRe$, as predicted by the long-wave asymptotics in §4.1.1; at large $kRe$, $c_{i}$ decreases as $(kRe)^{-2/5}$ and approaches the inviscid value. The asymptote $-2/5$ is obtained through a numerical fit.

In summary, viscosity can either suppress instabilities (as for $a=0$ and $a=-1$) or enhance them (as for $a=\pm 3$), before smoothly connecting to the inviscid rippling mode. This dual role underlines that viscosity does not merely act as a dissipative effect but can also generate significant growth depending on the flow configuration and parameter regime.

4.3 Growth rate behavior: a complete picture

Given the variety of instabilities that emerge in different asymptotic regimes, as discussed in the previous subsections, it is useful to examine the contours of growth rate ($kc_{i}$) in the $(k,Re)$ plane to obtain a broader picture. These contours (together with neutral stability curves) reveal the dominant instability mechanisms, their parameter ranges, and their mutual transitions. Figures 9a–h present the growth rate contours of the most unstable mode for four representative velocity profiles—$a=-3,-1,0,$ and $3$—for $G=0$ (left panels) and $G=0.02$ (right panels).

For the forward-bulging profile ($a=-3$), two broad instability regions appear (figures 9a,b), separated by a narrow band of stability. Three distinct modes can be identified within them. At small $Re$ and $k$, the long-wave interfacial (Yih) mode dominates, with peak growth rates at $Re\sim O(10)$ and $k\sim O(1)$. At large $Re$ and moderate $k$, a bright region corresponds to the rippling instability, consistent with figure 5a, where the wavenumber of maximum growth at $Re=10^{5}$ matches that of the inviscid case. A third zone, at higher $k$, corresponds to the short-wave interfacial mode, although its growth rates decay rapidly with increasing wavenumber. Including gravity and surface tension ($G=0.02$) strongly suppresses the short-wave mode and introduces a tongue-like region at very large $Re$, reminiscent of the shear mode described by Floryan et al. (1987). While present in the $G=0$ case as well, this mode is less visible there since it is not the most unstable. The overall growth rates, however, are diminished by the presence of gravity and surface tension.

The Nusselt profile ($a=-1$) lacks both rippling and short-wave instabilities, leaving only the long-wave interfacial mode and the shear mode at large $Re$ and $k\sim O(1)$ (figures 9c,d). Interestingly, at $G=0$, the long-wave mode extends into the short-wave regime at high $Re$, though this extension is curtailed when $G=0.02$, with gravity and surface tension acting to stabilize both long and short waves. The shear mode, however, remains relatively unaffected, reflecting its $O(1)$ wavenumber character.

For the linear profile ($a=0$), the instability landscape is quite different. Figures 9e,f show that the unstable region coincides precisely with the neutral curves of earlier studies: the outer boundary from Miles (1960) (blue dots) and Smith & Davis (1982) (magenta dots), and the inner boundary from Miesen & Boersma (1995) (black dots). This agreement validates both the present formulation and numerical approach. Importantly, the instability observed for $a=0$ is not of the long-wave, short-wave, or rippling types; instead, it resembles a shear-type mode (often termed “internal mode” by Boomkamp & Miesen (1996)), with properties more akin to those of the shear instabilities in figures 9b–d. The right-hand branch of the V-shaped contour, corresponding to short waves, disappears when $G=0.02$, while a new long-wave instability region emerges. Though weaker, this mode suggests that surface tension and gravity can, counterintuitively, play a destabilizing role, as also noted by Yiantsios & Higgins (1988). At higher $G$ (e.g. $0.5$), this mode vanishes, restoring the expected stabilizing effect.

The flow-reversal profile ($a=3$) displays contours (figures 9g,h) that qualitatively resemble those of $a=-3$: a long-wave interfacial instability, a short-wave interfacial instability, and a strong rippling mode at large $Re$. In this case, however, the rippling instability is more dominant, exhibiting growth rates that exceed those of the other modes, and persisting up to $Re\sim O(10^{2})$. Interestingly, at large $Re$ and $k$, a narrow band of stability appears within the rippling region, consistent with the deviations from inviscid predictions shown earlier in figure 8b. Unlike the other cases, no shear mode is observed for $a=3$.

In summary, the viscous stability analysis maps a diverse set of instabilities across the $(k,Re)$ parameter space. Long- and short-wave asymptotics identify the ranges of curvature $a$ that support interfacial instabilities, while the full growth rate contours clarify their coexistence and transitions. A key observation is that $a=0$ marks a sharp boundary in behavior, with small positive and negative curvatures favoring short- and long-wave instabilities, respectively. More generally, viscosity is not merely dissipative: for long waves it can enhance growth compared to inviscid predictions, while in other regimes it suppresses disturbances or generates new modes such as the shear-type instability.

These findings naturally raise several questions. What is the precise energy transfer mechanism that distinguishes the viscous long-wave mode from the inviscid rippling mode? How does the transition between interfacial and rippling instabilities proceed at finite Reynolds numbers, particularly for $a<-1$ and $a>0$? Under what conditions do gravity and surface tension play a destabilizing role rather than a stabilizing one? And finally, how robust are the shear-type modes in more realistic air–water systems with viscosity and density contrasts? These questions motivate a more detailed energy budget and eigenfunction analysis, which we turn to in the following section.

5 A family of instabilities

The characterization of instabilities in two-phase parallel flows is notoriously difficult. Even with the free-surface approximation (neglecting the top layer), the problem involves five parameters: $k,Re,G,Bo$, and $a$. A widely used diagnostic is the perturbation kinetic energy equation. Boomkamp & Miesen (1996), following Hooper & Boyd (1983), classified instabilities in two-phase flows by identifying the dominant contribution to perturbation kinetic energy production. For free-surface flows, Kelly et al. (1989) derived the analogous energy balance in the context of falling films. In this framework, the perturbation kinetic energy (KE) is expressed as the sum of work done by tangential stresses (TAN), normal stresses (NOR), Reynolds stresses (REY), and viscous dissipation (DIS):

[TABLE]

The individual terms are given by

[TABLE]

[TABLE]

where overbars denote horizontal averages over a wavelength, and TAN and NOR are evaluated at the free surface ($z=0$). With the normal mode ansatz and the normal stress boundary condition, these expressions can be evaluated. For $G=0$, NOR vanishes identically. To compare the relative importance of TAN and REY, we define

[TABLE]

so that $\varepsilon>0$ implies REY dominance and $\varepsilon<0$ implies TAN dominance.

Figures 10a-d show contour plots of $\varepsilon$ for the most unstable mode in the $(k,Re)$ plane, for $G=0$ and $a=3,-3,-1,$ and [math]. Yellow regions indicate REY dominance, blue regions TAN dominance. For $a=3$ and $a=-3$ (figures 10a,b), TAN dominates for long- and short-wave interfacial modes at low $Re$, consistent with Boomkamp & Miesen (1996). At higher $Re$, REY overtakes TAN, marking a change in mechanism (Charru & Hinch, 2000). The black curve corresponds to $\textrm{REY}=\textrm{TAN}$, while the red curve denotes the appearance of a critical layer ($c_{r}=U(z_{c})$), and the arrow points towards the region where a critical layer exists. However, whether the mode in this regime is a rippling mode or a shear mode cannot be distinguished from these figures alone. The $a=-1$ (Nusselt flow) case, shown in figure 10c, illustrates a different behavior: TAN dominates across all $Re$, with no gradual transition to REY dominance. The small REY-dominated patch at large $Re$ and $O(1)$ wavenumbers corresponds to a shear mode, consistent with earlier discussions. For $a=0$ (figure 10d), the instability resembles the shear mode observed for $a=-1$ and is clearly REY-dominated. This conclusion differs from Boomkamp & Miesen (1996), who argued for TAN dominance with REY as a small but finite contribution, motivating their term “internal mode.” Our results suggest instead that the mode is primarily driven by REY, though a critical layer persists across the parameter space.

These results highlight the limitations of energy analysis alone for classifying instabilities. In particular, the transition from long-wave interfacial to REY-dominated instability at high $Re$, and the precise distinction between rippling and shear instabilities, remain ambiguous. To refine this classification, Charru & Hinch (2000), building on Hooper & Boyd (1987), introduced two useful quantities: (i) the penetration depth of disturbances, and (ii) an effective Reynolds number measuring inertial influence via the ratio of the imaginary to real parts of the vorticity eigenfunction. They classified instabilities into three regimes: (a) a shallow-viscous regime for $k\ll 1$ (disturbances penetrate the depth), (b) a deep-viscous regime for $k\gg 1$ (disturbances confined to one wavelength), and (c) an inviscid regime where penetration depth is set by the viscous scale and the effective Reynolds number is $O(1)$. Transitions between these regimes are sharp and can be demonstrated by asymptotic calculations. Extending such a framework to free-surface flows lies beyond the scope of this study, but in the following subsections we adopt this terminology to describe representative eigenfunction structures.

5.1 Long wave modes ($k\ll 1$)

In §4.2.2, it was shown that the long-wave growth rates for the flow-reversal profile ($a=3$) behave differently at small and large Reynolds numbers, especially when compared with the inviscid limit. This difference can be partly attributed to the transition of the energy source from tangential stress (TAN) to Reynolds stress (REY) for long waves (figure 10a). To explore this further, figures 11a–c present streamlines and normalized vorticity contours at $k=0.01$ for $Re=1$, $10^{5}$, and the inviscid limit, respectively.

At $Re=1$ (figure 11a), the vorticity is nearly uniform across the depth, with no appreciable horizontal displacement, and the streamfunction is localized near the surface, in-phase with both vorticity and surface displacement $\eta$. This is characteristic of the shallow-viscous regime of Charru & Hinch (2000) and agrees with the mechanism described by Smith (1990): in the long-wave, low-$Re$ limit, the leading-order solution is a constant vorticity disturbance (see Appendix C). Here the fluid depth is much smaller than both the wavelength and the viscous length scale, so vertical variations are $O(1/k)$ and horizontal shifts are $O(kRe)$. In contrast, the $Re=10^{5}$ case (figure 11b)—which lies in the parameter region where growth rates exceed the inviscid values (figure 8b) and where REY dominates over TAN without a critical layer (figure 10a)—shows markedly different eigenstructures. The vorticity no longer spans the entire depth but instead peaks in the interior, decays smoothly below this point, and exhibits a horizontal shift indicative of inertia. The streamlines remain in-phase with vorticity but tilt rightward and concentrate more strongly near the surface. The reduced penetration depth suggests correspondence with the “inviscid regime” of Charru & Hinch (2000), in which disturbances are confined to scales smaller than both the wavelength and the fluid depth. Finally, the inviscid solution (figure 11c) reveals yet another distinct structure: a critical layer is now present (black dashed lines), unlike at $Re=10^{5}$, where vorticity simply peaked in the bulk. Outside this layer, vorticity decays rapidly, while the streamfunction resembles the $Re=1$ case away from the critical layer but acquires a pronounced tilt in its vicinity, much like the behavior shown in figure 5d.

In summary, the $Re=10^{5}$ eigenstructure is not merely a midpoint between the long-wave interfacial mode ($Re=1$) and the rippling mode (Inviscid), but rather represents a qualitatively distinct regime with its own dynamical character.

5.2 Short wave modes ($k\gg 1$)

Unlike the long-wave case, short waves do not undergo a gradual transition from TAN- to REY-dominated regimes (figure 10a). Instead, a distinct region of stability separates the two, suggesting that the underlying mechanisms are fundamentally different instabilities rather than a continuous evolution. To clarify this, we examine the eigenstructures at short wavelengths. Figures 11d–f present streamlines and normalized vorticity contours for $k=20$ at $Re=1$, $10^{5}$, and the inviscid case, respectively.

At $Re=1$ (figure 11d), the streamlines remain in phase with vorticity and are concentrated near the surface, with negligible horizontal shift due to weak inertia. However, unlike the long-wave case, the streamfunction exhibits a cellular pattern of nearly circular streamlines rather than surface-attached structures. This agrees with the physical mechanism of the short-wave interfacial instability described by Hinch (1984): the disturbance horizontal velocity is positive below crests and negative below troughs, a consequence of viscosity contrast across the interface. The vorticity penetration depth is of order one wavelength, consistent with the deep-viscous regime of Charru & Hinch (2000), where the viscous length scale exceeds the wavelength. The critical layer (black dashed line) lies close to the surface, in line with the eigenvalue expression (37).

At $Re=10^{5}$ (figure 11e), the critical layer shifts slightly below the surface and exhibits a localized vorticity maximum, reminiscent of a rippling mode, but concentrated within a narrow interfacial region like a shear mode. Thus, the eigenstructure combines features of both instabilities. The streamlines, anchored at the surface, show sharp gradients across the critical layer and mark a transition to the inviscid regime of Charru & Hinch (2000), where penetration depths are smaller than a wavelength and comparable to the viscous length scale. In the inviscid case (figure 11f), the vorticity maximum resides exclusively at the critical layer, and not at the surface, identifying the instability as a purely rippling mode. Compared with the long-wave inviscid case (figure 11c), the vorticity and streamfunction structures are confined to a shallower region of the domain due to the large $k$ and the proximity of the critical layer to the interface.

Taken together, these results show that short waves highlight a sharper distinction between viscous and inviscid mechanisms than long waves. At low $Re$, they correspond to the deep-viscous interfacial mode; at moderate to high $Re$, they acquire mixed features of shear and rippling modes; and in the inviscid limit, they reduce to the classical rippling instability.

5.3 Shear mode at $a=0$

Before examining eigenfunctions in the REY-dominant region for $a=3$, it is useful to first consider the shear mode for the linear profile ($a=0$). Figures 12a–c show streamlines and normalized vorticity contours for three representative points: $(k,Re)=(1,10^{3})$, $(6,10^{3})$, and $(1.75,150)$. The first two correspond to the left and right branches of the ‘V’-shaped instability region in figure 9e, while the third lies at their intersection. The case $(k,Re)=(1,10^{3})$ (figure 12a) exhibits the hallmark features of a shear mode: vorticity concentrated near the bottom wall with a horizontal shift (Hooper & Boyd, 1987), the presence of a critical layer (Miles, 1960), and streamfunction contours filling most of the domain, consistent with the ‘A’-family of modes (Kaffel & Riaz, 2015). The Reynolds stress $(-\overline{uw})$ shows a local extremum near the critical layer (Hooper, 1989), and the perturbation kinetic energy is supplied entirely by REY, confirming that this instability is a shear mode and not the “internal mode” suggested by Boomkamp & Miesen (1996).

For $(k,Re)=(6,10^{3})$ (figure 12b), the eigenstructure shifts to the surface: vorticity is concentrated near the interface, the critical layer sits close to the surface, and the Reynolds stress maximum is also located there. In this regime, where the viscous length scale is smaller than both the wavelength and the fluid depth, Hooper & Boyd (1987) identified three possible sources of growth: boundary layers at the wall, the interface, and the critical layer. For the linear profile, the critical layer plays no role, while the interface boundary layer was expected to be stabilizing, yet here the right-branch instability clearly originates at the surface. The intersection case $(k,Re)=(1.75,150)$ (figure 12c) combines wall- and surface-driven mechanisms: the critical layer is at mid-depth, and vorticity is localized at both boundaries but with a $180^{\circ}$ phase difference. Similar eigenstructures appear in the yellow region of figure 10c, confirming that it too corresponds to a shear mode, albeit without the right branch. In general, both left and right branches occur for curvature values close to $a=0$, but their extent diminishes differently as $a$ moves away from zero, leaving only the left branch at larger $a$ (as in figure 10c).

5.4 A composite mode at large $Re$

A deviation between the growth rates at $Re=10^{5}$ and the inviscid case was attributed in figure 8b to the large-$Re$ asymptote at long waves. Hence, the large-$Re$, small-$k$ region does not correspond to a pure rippling instability. Figure 10a further shows that REY dominates the perturbation KE in this regime, which rules out a purely long-wave interfacial mechanism, while the absence of a critical layer and the lack of vorticity concentration at the wall or interface exclude a shear mode. To probe this mixed behavior, figures 13a–c present eigenfunctions at $Re=10^{4}$ for $k=0.1,1,$ and $10$. At $k=0.1$ (figure 13a), the vorticity and streamfunction resemble a long-wave interfacial mode but with a horizontal tilt at depth, similar to figure 11b despite the order-of-magnitude difference in parameters. Increasing the wavenumber to $k=1$ (figure 13b) introduces two critical layers, one near the wall and another in mid-depth. Here, the vorticity shows a wall maximum like a shear mode, while the streamlines extend across the domain like a long-wave mode but tilt as in a rippling instability. At $k=10$ (figure 13c), the eigenstructure closely resembles figure 11e: vorticity peaks at the surface and also has a local maximum near a surface-critical layer, combining surface-driven shear behavior (cf. figure 12b) with rippling features (cf. figure 11f).

A comparison with the $a=-1$ case (figures 13d–f), which lacks a rippling instability, highlights the difference. Although the eigenfunctions deviate from the classical long-wave structure of figure 11a, the variations in tilt and penetration depth are monotonic with $k$ and $Re$, consistent with inertia-modified long-wave interfacial instability (Charru & Hinch, 2000). By contrast, the $a=3$ case exhibits a composite mode: behaving like a shear–interfacial hybrid at small $k$, showing properties of all three modes at moderate $k$, and resembling a shear mode modified by rippling at large $k$. Examination of the eigenspectrum confirms that this is not due to coalescence of separate modes (Özgen, 2008), but rather a single mode with varying character across parameter space. Similar dual behavior was noted by Albert & Charru (2000) in two-phase Couette flows, where an interfacial mode at small $Re$ evolved into a shear mode at large $Re$. The richer composite nature described here, however, does not appear to have been reported previously.

6 Conclusion

This work has presented a linear stability analysis of a class of Couette–Poiseuille (quadratic) shear flows in a liquid layer with a free surface. The base-state profiles are characterized by a curvature parameter $a$, which controls the degree of profile convexity or concavity. A systematic study was carried out in the $(a,k,Re)$ parameter space, with fixed inverse squared Froude number $G$ and Bond number $Bo$. Both inviscid and viscous mechanisms were analyzed, with emphasis on how viscous modes transition into their inviscid counterparts as $Re$ increases. Growth rate contours were mapped in $(k,Re)$ space, and eigenfunction structures were examined to classify the instabilities into five types: long-wave interfacial, short-wave interfacial, shear, rippling, and composite modes. The composite mode was found to exhibit hybrid characteristics of long-wave, shear, and rippling instabilities, bridging transitional regions between them.

In the inviscid limit, closed-form solutions of the Rayleigh equation were obtained in terms of confluent Heun functions (§3.1). Long-wave asymptotics for $G=0$ showed $c_{i}\sim k^{4/3}$ for $|a|>1$ and $c_{i}\sim k^{4}$ for $0<a<1$, while $a\in[-1,0]$ were stable. Both asymptotic and small-$a$ analyses highlighted the role of the critical layer in driving rippling instability. Modes were found to be prograde for $a<0$ and retrograde for $a>0$. Numerical results showed that short waves satisfy $c_{i}\sim k^{-2}$ for all $a\notin[-1,0]$. For finite $G$, both long- and short-wave cutoffs appear, and analytical expressions for the neutral stability boundary were derived (§3.3). Growth rate maps in $(a,k)$ revealed that instability is strongest for large $|a|$ and moderate $k$.

Viscous effects were explored using asymptotic analyses for long ($k\ll 1$) and short ($k\gg 1$) waves (§4.1). For $G=0$, velocity profiles near $a=0$ exhibited a sharp change: slightly concave profiles ($a<0$) supported long-wave instabilities, while slightly convex profiles ($a>0$) supported short-wave instabilities (figure 6). Asymptotic results further showed that $a\in[0,1]$ are stable to long waves, while $a\in[-1,0]$ are stable to short waves. Numerical solutions of the Orr–Sommerfeld problem confirmed these findings. For long waves, $c_{i}$ grew linearly with $Re$ at small values, peaked at moderate $Re$, and decayed towards the inviscid asymptote at large $Re$; viscosity thus enhanced long-wave growth by nearly an order of magnitude (figure 8). In contrast, for moderate and short waves, viscous growth rates increased monotonically and then converged to inviscid values. Complete stability diagrams in $(k,Re)$ (§4.3, figure 9) revealed four main regimes: long-wave interfacial, short-wave interfacial, shear, and rippling instabilities. The rippling and long-wave regions overlapped smoothly, while shear and short-wave instabilities remained sharply distinct.

To clarify the nature of these instabilities, we utilize the energy-budget analysis (Boomkamp & Miesen, 1996), the presence of a critical layer, and eigenfunction structure (Charru & Hinch, 2000). For $a=3$, the high-$Re$, low-$k$ regime with unusually high growth rates was shown to be driven mainly by Reynolds stress, distinguishing it from both rippling and long-wave interfacial modes. Eigenfunction analysis confirmed the physical mechanisms of long- and short-wave instabilities described by Smith (1990) and Hinch (1984) respectively. For $a=0$, the unstable mode was identified as a shear mode powered solely by Reynolds stress, in contrast to the “internal mode” classification of Boomkamp & Miesen (1996). For $a=3$, the high-$Re$ unstable branch displayed a mixture of long-wave, shear, and rippling characteristics, leading us to define a new composite mode (§5.4).

Figure 14a,b present a schematic summary of the identified instability regimes in the parameter spaces $(Re,k)$ and $(Re,a)$. The colour gradients indicate the relative strength of the instabilities. In both cases, the rippling instability, which arises in the inviscid limit, forms the upper bound of the instability map. In Figure 14a ($Re$–$k$ plane), the long- and short-wave instabilities appear at small to moderate $Re$, with the long-wave mode persisting to slightly higher $Re$. The shear mode dominates at large $Re$ and $O(1)$ wavenumbers, while the composite mode serves as a transition between the long-wave instability, the shear mode, and the rippling limit. In Figure 14b ($Re$–$a$ plane), the rippling instability is prominent for $a>0$ and $a<-1$. The long-wave mode is relevant in the ranges $a>1$ and $a<0$, while the short-wave mode is active for $a>0$ and $a<-1$. The shear mode is restricted to small values of $a$, with a slight asymmetry of greater importance for negative values of $a$ with its influence diminishing as $|a|$ increases.

The present study has focused on linear instabilities in a single-layer system. While many of the same mechanisms appear in two-layer flows (Charru & Hinch, 2000), the addition of a top layer introduces further modes such as Miles’ instability (Miles, 1957) and upper-wall shear instabilities (Shapiro & Timoshin, 2005), and can alter the stability boundary and critical Reynolds number of the bottom layer modes (Miesen & Boersma, 1995). Open questions remain: How does a top layer modify bottom-layer instabilities? What new composite modes arise from coupling between layers? Beyond this, the current analysis was restricted to two-dimensional disturbances, yet in regimes where $c_{i}$ grows with $Re$, three-dimensional or oblique instabilities may dominate (Yiantsios & Higgins, 1988; Mohammadi & Smits, 2017; Guha & Frigaard, 2010; Cheng et al., 2024). Non-modal growth (Abdullah & Park, 2024; Deng et al., 2025) and the competition between linear and nonlinear effects (Francius & Kharif, 2017) should also be explored further. An initial value problem approach (Patibandla et al., 2023) could reveal the role of different unstable modes in amplifying an initial disturbance and its later time evolution. While the present work is restricted to the linear regime, it serves as a precursor by identifying parameter ranges and destabilization mechanisms that can be further explored through full direct numerical simulations (DNS) to capture their nonlinear evolution.

\backsection

[Supplementary data]

\backsection

[Acknowledgements]

\backsection

[Funding]A.R. acknowledges financial support of the Science and Engineering Research Board, Government of India, under project nos. SPR/2021/000536 and MTR/2021/000706.

\backsection

[Declaration of interests]The authors report no conflict of interest.

\backsection

[Data availability statement]

\backsection

[Author ORCIDs]

\backsection

[Author contributions]H.M. and R.P. share equal contributions to the work

Appendix A Various confluent Heun function solutions of the Rayleigh equation with a quadratic background velocity profile

As mentioned earlier, the Rayleigh equation (12) with a quadratic velocity profile can be simplified to a confluent Heun differential equation (Kadam et al., 2023)

[TABLE]

where, $\phi(z)=x(x-1)e^{\beta x}f(x)$, $\beta=k\zeta_{c}/a$, $x=(\zeta(z)+\zeta_{c})/2\zeta_{c}$, $\zeta_{c}=\zeta(z_{c})$, $\zeta(z)$ is the base-state vorticity and $z_{c}$ is defined as the location that satisfies $U(z_{c})=c$ (here $c$ is complex). The auxiliary parameters in the confluent Heun equation above are $4\beta,2\beta,2,2$ and $2\beta$, respectively. Equation (39) possesses two regular singularities at [math] and $1$ and an irregular singularity at $\infty$. One can write a Frobenius solution to (39) at $x=0$, with indicial exponents [math] and $-1$. The first linearly independent series solution can be extended for $|x|>1$ using analytic continuation (Motygin, 2018) and is defined as the confluent Heun function (‘HeunC’ in Mathematica). Note that the second series solution, with negative indicial exponent, indicates logarithmic terms in the series and is not available in popular software applications. One can instead use the series solution at the regular singularity $x=1$ as the second linearly independent solution of the differential equation (39). A Mobius transformation $x^{\prime}=x-1$, provides a confluent Heun function that is centered at $x=1$. Note that, due to the nature of the analytic continuation algorithm, evaluating the confluent Heun function at $x=1$ results in an indeterminate expression. This can be circumvented by using a confluent Heun function centered at $x=1$. Here, we present various solution forms of the Rayleigh equation (12), with the quadratic velocity profile, that can be evaluated in Mathematica.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The two terms in each of the expression sets above are linearly independent with each other. Solution forms (41) are presented in § 3.1. Expressions (42)-(43) are obtained following Ishkhanyan (2024). One can write any two forms from the seven expressions above as a solution to (12) although their linear independence should be checked. One can also change the sign of $\beta$ in the expressions above, to obtain additional solutions. There are indeed many solution forms (Slavianov & Lay, 2000) and listing them here is not the intention of the current work. Rather, it is to point out that different sets of solution forms can be chosen to ease the numerical calculation of the confluent Heun function at different regions in the parameter space. In the current work the seven solution forms above have been used.

Appendix B Small-curvature inviscid asymptotics

We assume $a$ to be a small parameter in the Rayleigh equation (12) and the boundary condition (14). A regular perturbation expansion in terms of $a$, of the eigenfunction $\phi$ and eigenvalue $c$ are given by

[TABLE]

The subscript $a\ll 1$ will be dropped henceforth for convenience. The leading order eigenfunction and eigenvalue correspond to the linear base-state flow and are given by

[TABLE]

To the next order, the Rayleigh equation can be simplified as

[TABLE]

Using the variation of parameters, a particular solution that satisfies the bottom boundary condition can be written as

[TABLE]

where,

[TABLE]

The boundary condition at $O(a)$ is,

[TABLE]

where the apastrophe indicating a derivative with $z$. Simplifying the expression above yields

[TABLE]

Please note here that, $\phi^{(0)}(0)=1$, $\phi^{\prime(0)}(0)=k\coth{k}$. For the existence of a growth rate, the RHS of (50) should be complex. This, in-turn, requires complex $g_{1}(0)$ and $g_{2}(0)$. Inspecting the integrals in (48), one finds that if the value of $c_{0}$ doesn’t match the flow velocity (somewhere in the domain) then equation (47) and hence (50) will be real. One can evaluate the imaginary parts of $g_{1}(0)$ and $g_{2}(0)$ using Plemelj theorem in the current case, as outlined by Shrira (1993), as

[TABLE]

Therefore, one obtains,

[TABLE]

[TABLE]

Finally, substituting the results above into (50), an expression for growth rate can be written as

[TABLE]

This expression is valid for all $k$ and $G$, but at small $a$. It depends on gravity and surface tension implicitly through $c^{(0)}$. Interestingly, one can find that $kc_{i}$ is positive for positive $a$ and vice versa. In otherwords, velocity profiles with $U^{\prime\prime}>0$ are unstable, however small the curvature may be (Shrira, 1993; Bonfils et al., 2023). Even a small negative curvature renders the mode stable. This explains the stark contrast across $a=0$ mentioned in § 3.2, and § 3.4. Finally, as shown in figure (15), $c_{i}$ calculated from the asymptotic expression (54) matches well with the complete numerical solution for small $a$ (= 0.1).

Appendix C Viscous long wave asymptotics

Substituting (32) into the Orr-Sommerfeld equation (7) and the corresponding viscous boundary conditions (8)-(10), one obtains the following set of equation at the leading order,

[TABLE]

[TABLE]

[TABLE]

The subscript $k\ll 1$ will be dropped henceforth for convenience. Solving equation (55) along with its boundary conditions yields the leading-order solution. An undetermined constant, which appears due to the homogeneity of the equation, is set to 1, resulting in:

[TABLE]

The order $k$ governing equations are given by,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The zeroth order solution will appear in all subsequent orders. To eliminate repetition, the coefficient multiplying the zeroth order solution can be set to zero at each order, without loss of generality. The solution at $O(k)$ can then be written as

[TABLE]

The $O(k^{2})$ governing equation is given by,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The calculation is terminated at $O(k^{2})$ because an unstable mode is already identified at $O(k)$. This provides the solution for the eigenvalue $c$ up to $O(k^{2})$. The process can, however, be extended further to determine higher-order solutions if needed. The leading order ($k=0$) solution for vorticity ($\Omega$) can be found by differentiating $\phi^{(0)}(z)$ from (56) twice with respect to $z$, which shows that the vorticity ($\Omega$) is constant.

Appendix D Viscous short wave asymptotics

As outlined in §4.1.2, the governing equations for the short wave asymptotics are first scaled by $k$, with a small parameter ($\epsilon=k^{-1}$). The solution is then expressed as a series expansion in powers of $\epsilon$. By substituting the expansion (35) into the scaled equations and the corresponding boundary conditions—where, in the scaled coordinate system, the domain extends from [math] at the interface to $-\infty$ at the wall—the leading-order governing equation and the boundary conditions can be written as:

[TABLE]

[TABLE]

[TABLE]

Solving the set of equations (60), one obtains the leading order solution as,

[TABLE]

Henceforth, the subscript $k\gg 1$ will be dropped for convenience. For the $O(\epsilon)$ terms, the governing equation and the boundary conditions become,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Similar to the long wave asymptotics, the solution must be free of a linear multiple of the $O(1)$ solution. After solving the $O(\epsilon)$ equation along with its boundary conditions, one obtains

[TABLE]

The calculation is carried out up to $O(\epsilon^{4})$, where the eigenvalue reveals an instability.