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

(2.15) \begin{equation} U(z) = a z^2 + (a+1)z + 1, \end{equation}

where $a$ is a curvature parameter that may take values in the range $-\infty \lt a \lt \infty$ . Any velocity profile of the bottom layer in a two-phase Couette–Poiseuille flow can be recast into this non-dimensional form, making it a convenient canonical representation (see Appendix A for a derivation). 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$ , (2.15) reduces to a linear velocity profile, characteristic of a simple shear flow. The case $a=-1$ yields the half-parabolic profile (also called Nusselt profile), which typically arises in the classical problem of a liquid film flowing down an inclined plane under gravity and can also be seen in a two-phase Couette–Poiseuille flow for certain viscosity and height ratios. For $a=3$ , the profile exhibits a flow reversal near wall, a phenomenon observed in laboratory experiments on wind-driven water waves (Hidy & Plate Reference Hidy and Plate1966; Paquier et al. Reference Paquier, Moisy and Rabaud2015, Reference Paquier, Moisy and Rabaud2016) on water layers confined in the streamwise direction leading to a net zero flow rate. Interestingly, similar flow reversals occur in large-scale geophysical settings, such as equatorial oceanic currents, though with more complex vertical structures (Stommel Reference Stommel1959; Constantin & Johnson Reference Constantin and Johnson2016). Velocity profiles with $a\lt -1$ can be seen, for example, in falling films with counter-current gas flows in the top layer (Vellingiri, Tseluiko & Kalliadasis Reference Vellingiri, Tseluiko and Kalliadasis2015; Schmidt et al. Reference Schmidt, Ó Náraigh, Lucquiaud and Valluri2016; Ishimura et al. Reference Ishimura, Mergui, Ruyer-Quil and Dietze2023). 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 divided into four distinct regimes, each corresponding to a qualitatively different profile shape: (i) backward bulging ( $1 \leqslant a \lt \infty$ ), which includes the flow-reversal case; (ii) monotonic convex ( $0 \lt a \lt 1$ ); (iii) monotonic concave ( $-1 \lt a \leqslant 0$ ), which contains the linear shear profile; (iv) forward bulging ( $-\infty \lt a \leqslant -1$ ), which includes the Nusselt profile. Representative examples of these regimes are illustrated in figure 2.

Figure 2.

Base-state velocity profiles for different values of the curvature parameter $a$ , highlighting the four regimes: I – backward bulging ( $a \in [1,\infty )$ ); II – monotonic convex ( $a \in (0,1)$ ); III – monotonic concave ( $a \in [-1,0]$ ); IV – forward bulging ( $a \in (-\infty ,-1]$ ).

Figure 2. Long description

A line graph showing base-state velocity profiles for different values of the curvature parameter, highlighting the four regimes: backward bulging, monotonic convex, monotonic concave, and forward bulging. The x-axis represents the variable U(z), and the y-axis represents the variable z. The graph includes multiple lines representing different values of the curvature parameter: a equals 3, a equals 0.7, a equals 0, a equals negative 0.7, and a equals negative 3. Each line corresponds to a different regime, with the regions labeled as I, II, III, and IV. The graph illustrates how the velocity profiles change across these regimes. All values are approximated.

3.1. Analytical solution

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

(3.4) \begin{equation} \phi (z) = e^{\kappa (\zeta (z)-\zeta _c)}\big(\zeta _c^2-\zeta (z)^2\big)\left (C_1 \phi _1(z) + C_2\phi _2(z)\right ), \end{equation}

where

(3.5) \begin{align} \phi _1(z, \zeta _c) &= \dfrac {{\textrm{H}}_{\textrm{c}}\left [2+4\kappa \zeta _c,4\kappa \zeta _c,2,0,4\kappa \zeta _c,\dfrac {\zeta _c+\zeta (z)}{2\zeta _c}\right ]}{\zeta _c-\zeta (z)},\hspace {.5cm}\phi _2(z, \zeta _c) = \phi _1(z, -\zeta _c). \nonumber \\ \zeta (z) &= U^{\prime}(z) = 2az + (a+1),\hspace {1cm} \zeta _c = \zeta (z_c),\hspace {1cm} \kappa = k/2a, \end{align}

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 B. 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 (3.2) and (3.3), giving

(3.6) \begin{equation} \varXi _w(c-1)^2 + (a+1)(c-1) - G\big(1+k^2Bo^{-1}\big) = 0, \end{equation}

where

(3.7) \begin{equation} \varXi _w = k + \dfrac {a+1}{1-c} + \dfrac {\phi _1(-1)\phi^{\prime}_2(0)-\phi _2(-1)\phi^{\prime}_1(0)}{\phi _1(-1)\phi _2(0)-\phi _2(-1)\phi _1(0)}. \end{equation}

For prescribed $(a,k,c,Bo)$ , (3.6) 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 analysed directly without numerical evaluation. In the sections that follow, we examine its implications: § 3.2 addresses long-wave asymptotics for the growth rate, while § 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 Reference Burns1953), which connects the wave phase speed to the Froude number in the limit $k\to 0$ :

(3.8) \begin{equation} \int _{-1}^{0} \dfrac {1}{(U(z) - c)^2} \, {\textrm d}z - \dfrac {1}{G} = 0. \end{equation}

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

(3.9) \begin{equation} G\left (\alpha (1+a(2c-1)) + 2ac(c-1)\log \left (\dfrac {2c-1+a+\alpha }{2c-1+a-\alpha }\right )\right )-c(c-1)\alpha ^3 = 0, \end{equation}

with $\alpha ^2 = (a-1)^2 + 4 a c$ .

For $G\gt 0$ and every value of curvature ( $a$ ), (3.9) yields two neutrally stable long-wave modes ( ${k=} \,c_i=0$ ): a prograde branch with speed exceeding the maximum background velocity, and a retrograde branch with speed below the minimum velocity (figure 3 a). Thus, long waves remain neutrally stable for any $G\gt 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 3 a). With the maximum given by $ (-(a-1)^2/4a$ for $a\leqslant -1$ and 1 for $a\gt -1 )$ and minimum $ (0$ for $a\leqslant 1$ and $-(a-1)^2/4a$ for $a\gt 1 )$ . 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 Reference Young and Wolfe2014). Notably, when $a\lt -1$ , the prograde mode itself can destabilize due to the forward bulge in the profile, consistent with the observations of Bonfils et al. (Reference Bonfils, Mitra, Moon and Wettlaufer2023).

Figure 3.

(a) Long-wave phase speeds obtained from (3.9) for different inverse squared Froude numbers ( $G$ ), as a function of curvature parameter $a$ . The shaded blue region corresponds the background velocity values $U(z)$ for each $a$ . For $G=0$ , (b) the real part $c_r$ , and (c) the imaginary part $c_i$ of the complex phase speed computed by numerically solving (3.6) are plotted against wavenumber $k$ for various $a$ . Black curves with markers in (c) indicate the asymptotes calculated in § 3.2, while the grey dashed lines mark the location of maximum $c_i$ for each $a$ .

Figure 3. Long description

The image contains three graphs. The first graph (a) shows long-wave phase speeds as a function of curvature parameter for different inverse squared Froude numbers. The shaded blue region indicates background velocity values. The second graph (b) plots the real part of the complex phase speed against wavenumber for various values. The third graph (c) shows the imaginary part of the complex phase speed against wavenumber, with black curves indicating asymptotes and grey dashed lines marking the location of maximum growth rate for each value. All values are approximated.

Figure 3(b,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$ (figure 3 b), computed from the dispersion relation (3.6), 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$ (figure 3 c) demonstrates finite growth rates for any finite non-zero $k$ , with maxima occurring at $O(1)$ wavenumbers, a typical feature of rippling instability (Young & Wolfe Reference Young and Wolfe2014). 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 (Reference Renardy and Renardy2013), using the leading-order pressure perturbation solution as input to the pressure perturbation equation at the next order. The governing equation is

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

subject to

(3.12a,b) \begin{align} k^2 (U(0)-c)^2 p(0) - G\big(1 + Bo^{-1} k^2\big)p^{\prime}(0) = 0, \quad p^{\prime}(-1) = 0, \end{align}

which reduce to $p(0)=0$ and $p^{\prime}(-1)=0$ for $G=0$ . One can obtain (3.10) by substituting $\phi = p^{\prime}(z)/(U-c)$ in (3.1) and writing in the self-adjoint form.

1. (i) Case I ( $|a|\gt 1$ ). Substituting the phase speeds from (3.9) into (3.10), the $k=0$ solution reads

   (3.13) \begin{equation} p_0(z)= \begin{cases} 1 - \left (\dfrac {\zeta (z)}{1+a}\right )^5,& 0 \geqslant z \gt z_c,\\ 1, & z_c \gt z \geqslant -1, \end{cases} \end{equation}
   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 (3.10) yields the relation
   (3.14) \begin{equation} \dfrac {p^{\prime}(0)}{(1-c)^2} = k^2 \int _{-1}^{0} \dfrac {p}{(U -c)^2} \, {\textrm d}z. \end{equation}
   Using $p(z)\approx p_0(z)$ and setting $c=U_{min /max }+c_1$ with $|c_1|\ll 1$ gives
   (3.15) \begin{equation} -\frac {160 a^3}{(a+1)^5} \sim -k^2 \frac {i \pi }{2 \sqrt {a} (c_1)^{3/2}}, \end{equation}
   leading to
   (3.16) \begin{equation} c_1 \sim \left (\dfrac {\pi (a+1)^5}{320 a^3\sqrt {|a|}}\right )^{2/3}k^{4/3}e^{i\pi /3}. \end{equation}
   The non-integer power law ( $c_i \sim k^{4/3}$ ) reflects the singular role of the critical layer. This generalizes earlier results of Kadam et al. (Reference Kadam, Patibandla and Roy2023) and recovers Renardy & Renardy (Reference Renardy and Renardy2013) in the $|a|\to \infty$ limit.

2. (ii) Case II ( $|a|\lt 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 (Reference Renardy and Renardy2013), we obtain a neutral solution at $c=0$ , leading to the asymptotic relation

   (3.17) \begin{equation} \bar {c} \sim \dfrac {\bar {k}^2p_0(-1)}{2} + \dfrac {\bar {k}^4p_0^2(-1)}{8}\left (\pi i - \ln \left (\dfrac {\bar {k}^2p_0(-1)}{4}\right )\right ) + O(\bar {k}^6), \end{equation}
   with $\bar {k}=(1-a)k/2a$ , $\bar {c}=4ac/(1-a)^2$ and
   (3.18) \begin{equation} p_0(-1) = \left (\dfrac {\zeta (0)^5-\zeta (-1)^5}{5(1-a)^5}-\dfrac {2\zeta (0)^3-\zeta (-1)^3}{3(1-a)^3} +\dfrac {\zeta (0)-\zeta (-1)}{1-a}\right ). \end{equation}
   This gives $c_i\sim k^4$ for $0\lt a\lt 1$ , while $a\lt 0$ yields stability. Figure 3(c) (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 (Reference Kaffel and Renardy2011). More generally, studies by Shrira (Reference Shrira1993) and Bonfils et al. (Reference Bonfils, Mitra, Moon and Wettlaufer2023) show exponential decay for $G\gt 0$ , with the instability mechanism tied to curvature at the critical layer. A complementary small-curvature asymptotic analysis is given in Appendix C.

Thus, long waves are neutrally stable for any non-zero $G$ (i.e. $c_i=0$ at $k=0$ ), 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 (Reference Yih1972), 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. Reference Morland, Saffman and Yuen1991; Engevik Reference Engevik2000; Renardy & Renardy Reference Renardy and Renardy2013; Young & Wolfe Reference Young and Wolfe2014; Kadam et al. Reference Kadam, Patibandla and Roy2023) 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.

1. (i) Case I ( $|a|\gt 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 (3.1) yields the neutral eigenfunction

   (3.19) \begin{equation} \phi (z)= \begin{cases} \dfrac {(1+a)\left (\kappa \zeta (z) \cosh [\kappa \zeta (z)] - \sinh [\kappa \zeta (z)]\right )}{\zeta (z) \left (\kappa (1+a) \cosh [\kappa (1+a)] - \sinh [\kappa (1+a)]\right )}, & 0 \geqslant z \gt z_c,\\ 0, & z_c \gt z \geqslant -1, \end{cases} \end{equation}
   with $\zeta (z)$ and $\kappa$ defined in (3.4). Enforcing the surface condition at $z=0$ gives the stability boundary
   (3.20) \begin{equation} G = \dfrac {(1+a)^3}{16a^2(1+k^2Bo^{-1})}\left (-6a + \dfrac {(1+a)^2 k^2}{k(1+a)\coth [k(1+a)/2a]-2a}\right ). \end{equation}

2. (ii) Case II ( $0\lt a\lt 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 (3.6) gives

   (3.21) \begin{equation} G = \dfrac {1}{1+k^2 Bo^{-1}} \left (\dfrac {a \textrm{H}^{\prime}_c\big [2+4\kappa (a-1), 4\kappa (a-1), 2, 0, 4\kappa (a-1), \frac {a}{a-1}\big ]}{(a-1)\textrm{H}_c\big [2+4\kappa (a-1), 4\kappa (a-1), 2, 0, 4\kappa (a-1), \frac {a}{a-1}\big ]}+k-a\right ), \end{equation}
   where $\textrm{H}^{\prime}_c$ denotes the derivative of the confluent Heun function.

3. (iii) Case III ( $-1\lt a\lt 0$ ). In this range, the prograde mode remains stable. Substituting $c=1$ into the boundary condition (3.3) simplifies to

   (3.22) \begin{equation} G \left (1 + \dfrac {k^2}{Bo}\right )= 0 \quad \Rightarrow \quad G = 0, \end{equation}
   demonstrating that no neutral mode exists for finite $G$ .

In summary, (3.20)–(3.21) 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 4(a). For $G=0$ , $k_{min}=0$ across all $a$ , consistent with the long-wave analysis. For $G\gt 0$ and $a\gt 0$ , $k_{min}$ generally grows with $G$ but varies non-monotonically with $a$ , peaking near $a\approx 1.41$ . For $a\lt 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}\gt 3$ uniformly in cyan. The white band at $0\lt a\lt 1$ reflects the absence of instability, in agreement with (3.22). 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 3(c).

Figure 4.

(a) Contour plot of the long-wave cutoff $k_{min}$ in the $(G,a)$ plane. The region with $k_{min}\gt 3$ is shown in cyan to highlight the $k_{min}\lt 3$ region. For $G=0.02$ , (b) the real part $c_r$ and (c) the imaginary part $c_i$ of the complex phase speed, plotted against wavenumber $k$ for different $a$ . Grey dashed lines in (b) indicate the extrema of the base-state velocity for each $a$ , and the crossings with solid lines mark the emergence of critical layers.

Figure 4. Long description

The image contains three graphs analyzing wave instabilities in two-phase shear flows. The first graph is a contour plot showing the long-wave cutoff in the plane, with the region where kmin is greater than 3 highlighted in cyan. The second graph plots the real part of the complex phase speed against wavenumber for different values of a, with grey dashed lines indicating the extrema of the base-state velocity for each a. The third graph plots the imaginary part of the complex phase speed against wavenumber for different values of a. The graphs illustrate how different parameters affect the stability and phase speed of waves in shear flows.

3.4. Complete solution of the dispersion relation

As outlined in § 3.1 and Appendix B, the dispersion relation (3.6) can be solved numerically across a wide parameter space using a root-finding routine in Mathematica (Wolfram Research 2024), yielding the complex phase speed. Figure 3(b,c) shows $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 figure 4(b,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}$ (grey dashed lines in figure 4 b). For $G=0$ , $k_{min}=0$ , so $c_r$ coincides with $U_{min/max}$ at onset. For $k\gt k_{min}$ , modes are unstable, and both figure sets (figures 3 and 4) show the development of a critical layer. Unlike the $G=0$ case (figure 3 b), the critical layer does not approach the surface asymptotically when $G\gt 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 4(b), the $U_{max}$ line (grey) 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.

Figure 5.

Growth-rate contours ( $kc_i$ ) in the $(a,k)$ plane for $G=$ (a) 0, (b) 0.02 and (c) 0.5. Grey dashed lines in (a–c) trace the locus of maximum growth rates. Black curves in (b) and (c) mark the stability boundaries from § 3.3. Growth rates below $10^{-10}$ are omitted. Panel (d) shows streamlines overlaid on normalized vorticity contours ( $\varOmega /\varOmega _{max}$ ) for $G=0$ , $a=-3$ and $k=6.31$ – black star in panel (a). Black dashed lines indicate critical-layer locations.

Figure 5. Long description

The image contains four panels showing growth-rate contours in the plane for different values of G. Panel (a) shows the contours for G equals 0, with a grey dashed line tracing the locus of maximum growth rates. Panel (b) shows the contours for G equals 0.02, with black curves marking the stability boundaries. Panel (c) shows the contours for G equals 0.5, also with black curves marking the stability boundaries. Panel (d) shows streamlines overlaid on normalized vorticity contours for specific parameters, with black dashed lines indicating critical-layer locations. The grey dashed lines in panels (a) to (c) trace the locus of maximum growth rates. All values are approximated.

Figure 5(a–c) show growth-rate contours ( $kc_i$ ) in the $(a,k)$ plane for $G=0$ , $0.02$ and $0.5$ . For $a\lt 0$ , only the prograde branch is unstable, whereas for $a\gt 0$ the retrograde branch is unstable. As anticipated from (3.22), profiles with $-1\lt a\lt 0$ are always stable. For $G=0$ , instability exists for all $k$ whenever $a\lt -1$ or $a\gt 0$ . Any non-zero $G$ produces a finite region of stability, which grows with increasing $G$ . Thus, quadratic velocity profiles with $a\lt -1$ or $a\gt 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 figure 5(a–c) show the stability boundaries derived in § 3.3, in good agreement with the growth-rate contours. The grey dashed lines mark the locus of maximum growth rate in the $(a,k)$ plane. For retrograde modes ( $a\gt 0$ ), the maximum growth rate occurs at wavenumbers $O(1)$ . For prograde modes ( $a\lt -1$ ), the most unstable wavenumber increases as $a$ increases towards $-1$ . This reflects the diminishing forward bulge of the velocity profile, which vanishes at $a=-1$ . Away from the stability boundary in figure 5(a–c), growth rates increase with $a$ (see also figure 3 c), consistent with the dependence of instability strength on the velocity curvature at the critical layer (Shrira Reference Shrira1993; Bonfils et al. Reference Bonfils, Mitra, Moon and Wettlaufer2023).

Figure 5(d) illustrates the streamline pattern overlaid on normalized vorticity contours ( $\varOmega /\varOmega _{max}$ ) for $G=0$ , $a=-3$ , and $k=6.31$ (black star in figure 5 a). Two critical layers are visible (black dashed lines), since $c_r=U(z)$ is satisfied at two distinct depths for forward-bulging profiles ( $a\lt -1$ ). Around these layers, the stream function tilts and vorticity contours distort strongly, consistent with the negative Reynolds stress correlation $\overline {uw}\lt 0$ that drives the instability (Young & Wolfe Reference Young and Wolfe2014). Note that the overbar denotes spatial average over a wavelength.

To summarize, the inviscid stability analysis demonstrates that for $G=0$ , $c_i\sim k^{4/3}$ when $|a|\gt 1$ and $c_i\sim k^4$ when $|a|\lt 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\lt a\lt 0$ . With increasing $G$ , long-wave cutoffs emerge, stabilizing low-wavenumber disturbances. Instabilities propagate faster than the surface velocity when $a\lt -1$ and slower when $a\gt 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 Reference Miles1957; Stern & Adam Reference Stern and Adam1973; Young & Wolfe Reference Young and Wolfe2014), justified by the typically large Reynolds numbers of geophysical flows. In contrast, laboratory experiments on thin layers (Paquier et al. Reference Paquier, Moisy and Rabaud2015, Reference Paquier, Moisy and Rabaud2016) involve Reynolds numbers of $O(10^3)$ , making viscous effects non-negligible. In addition, the layer near the bottom wall can become unstable resulting in new instability modes. While viscous stability analyses are common in engineering, they have rarely addressed the rippling instability directly (see table 1 of Young & Wolfe (Reference Young and Wolfe2014) 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.1. Asymptotic evaluation of growth rates

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

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

(4.1) \begin{equation} \phi (z) = \sum _{n = 0}^{\infty } k^n \phi ^{(n)}_{k \ll 1}(z), \qquad c = \sum _{n = 0}^{\infty } k^n c^{(n)}_{k \ll 1}. \end{equation}

Substituting these expansions into (2.11)–(2.14) and solving the resulting equations order by order in $k$ , as outlined in Appendix D (see also Yih Reference Yih1967), yields explicit expressions for the first three terms in the eigenvalue,

(4.2) \begin{align} c^{(0)}_{k \ll 1} &= 1-a, \qquad c^{(1)}_{k \ll 1} = \frac {iRe}{15}\left (4a(a-1) - 5G\right ), \nonumber \\ c^{(2)}_{k \ll 1} &= \frac {1}{1008} \big (179 a^3 Re^2 - 256 a^2 Re^2 - 222 a G Re^2 + 77 a Re^2 + 2016 a + 98 G Re^2\big ). \end{align}

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 and surface tension ( $G=0$ ), the imaginary part of the eigenvalue reduces to

(4.3) \begin{equation} c_i = \frac {4}{15}\,kRe\,a(a-1) + O(k^3). \end{equation}

Equation (4.3) shows that long-wave (Yih-type) instability occurs when $a\lt 0$ or $a\gt 1$ , whereas velocity profiles between the linear case ( $a=0$ ) and the half-parabolic case ( $a=1$ ) remain stable. Moreover, $c_i$ scales linearly with $kRe$ , the natural small parameter in this expansion (Yih Reference Yih1967). Notably, from $c^{(0)}_{k \ll 1}$ in (4.2), the leading-order phase speed exceeds the maximum velocity in the domain for $a\lt 0$ and falls below the minimum velocity for $a\gt 1$ , revealing an asymmetry of propagation depending on profile curvature. For non-zero $G$ , expression (4.2) indicates a stabilization due to gravity for $G\gt 4a(a-1)/5$ . For a given $G$ , the range of curvatures ( $a$ ) bounded by $(1-\sqrt {1+5G})/2$ and $(1+\sqrt {1+5G})/2$ are stable to the long-wave instability.

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

A similar analysis to the long-wave limit can be carried out for short waves. This approach was first developed by Hooper & Boyd (Reference Hooper and Boyd1983) for two-layer linear–linear profiles and later generalized by Yiantsios & Higgins (Reference Yiantsios and Higgins1988) 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$ :

(4.4) \begin{equation} \phi (\tilde {z}) = \sum _{n=0}^{\infty } \epsilon ^n \phi _{k \gg 1}^{(n)}(\tilde {z}), \qquad c = \sum _{n=0}^{\infty } \epsilon ^n c_{k \gg 1}^{(n)}. \end{equation}

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

(4.5) \begin{align} c_{k \gg 1}^{(0)} &= 1 - \frac {i}{2Ca}, \!\quad\! c_{k \gg 1}^{(1)} = - \frac {3 i Re}{16 Ca^2}, \!\quad\! c_{k \gg 1}^{(2)} = \frac {1}{32 Ca}\left [16(a+1)Re - 5 i \frac {Re^2}{Ca^2}\right ] {- \dfrac {i}{2}GRe}, \end{align}

(4.6) \begin{align} c_{k \gg 1}^{(3)} &= \frac {Re}{512 Ca}\left [-416 a + 348 (a+1)\frac {Re}{Ca} - 192 i G Re - 85 i \frac {Re^2}{Ca^3}\right ], \nonumber \\ c_{k \gg 1}^{(4)} &= \frac {(a+1)Re^3}{Ca^3} + \frac {23 i (a+1)^2 Re^2}{32 Ca} + \frac {(a+1)}{2}GRe^2 + \frac {i Re a (a+1)}{2} - \frac {79 a Re^2}{64 Ca^2} \nonumber \\ &\quad - \frac {407 i Re^4}{2048 Ca^5} - \frac {15 i G Re^3}{32 Ca^2}. \end{align}

Here $\textit{Ca} = \textit{Bo}/(\textit{GRe})$ is the capillary number.

From (4.5) 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 and surface tension ( $G=0$ ), the imaginary part of $c$ simplifies to

(4.7) \begin{equation} c_i = \dfrac {1}{2}\,a(a+1)\textit{Re}\,\epsilon ^4 + O(\epsilon ^5). \end{equation}

Thus, all lower-order imaginary terms vanish at $G=0$ , leaving only the destabilizing $O(\epsilon ^4)$ contribution. Instability arises for $a\gt 0$ and $a\lt -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. One can find the effects of surface tension by taking the distinguished limit of $G \to 0$ while keeping $G/Bo$ finite. Such a scenario may be possible, for example in Liu et al. (Reference Liu, Chen, Bond and Hu2017) (see table 1). This limit can be achieved by setting $G = 0$ and keeping $Ca$ finite in (4.5). One finds, straightaway, that a non-zero $Ca^{-1}$ results in a stabilizing term in the leading order itself. This implies that short-wave instability is more widespread in parameter space when gravity and surface tension are absent, and that any non-zero $G{/Bo}$ suppresses it through the stabilizing terms in (4.5). Finally, collecting the real parts of (4.5), one finds that the short-wave mode travels slightly faster than the surface velocity for $a\gt -1$ and slightly slower than the surface velocity for $a\lt -1$ .

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 analyse the behaviour at intermediate parameters, a numerical approach is required, since (2.11) admits no closed-form solution. We therefore employ a Chebyshev spectral collocation method following Boomkamp et al. (Reference Boomkamp, Boersma, Miesen and Beijnon1997) and Trefethen (Reference Trefethen2000) to solve the eigenvalue problem (2.11)–(2.14). A comparison of current numerical results with calculations of Miesen & Boersma (Reference Miesen and Boersma1995) are given in Appendix F showing a good match.

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 $0\lt a\lt 1$ and short-wave instabilities are absent for $-1\lt a\lt 0$ 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\gt 0$ ) develop short-wave instabilities, while slightly concave profiles ( $a\lt 0$ ) develop long-wave instabilities. This behaviour is illustrated in figure 6, which shows $c_i$ as a function of $k$ for $a=-0.05$ and $a=0.05$ at $\textit{Re}=10$ . The solid lines denote the numerical solutions, while the dashed lines correspond to the asymptotic predictions (4.3) and (4.7). 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_{{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. The origin of this stark asymmetry across $a = 0$ can be understood physically. For a horizontal liquid layer driven by wind, a slight inclination of the container introduces a wall-parallel component of gravity: tilting in one direction yields $a \gt 0$ (short-wave unstable), while the opposite tilt yields $a \lt 0$ (long-wave unstable). This correspondence is reflected in the abrupt sign change of $c_i$ across $a = 0$ in figure 6 and in the contrasting stability boundary structure visible in figure 6.

Figure 6.

The variation of $c_i$ with wavenumber $k$ shown for velocity profiles slightly perturbed from a linear profile: $a = 0.05$ corresponds to a mildly convex profile, and $a = -0.05$ corresponds to a mildly concave profile. Solid lines represent numerical results, and dashed lines indicate the asymptotic predictions from §§ 4.1.1 and 4.1.2.

Figure 6. Long description

The line graph illustrates the variation of with wavenumber for velocity profiles slightly perturbed from a linear profile. The x-axis represents the wavenumber (k) on a logarithmic scale ranging from 10^−2 to 10^1, while the y-axis represents the growth rate (c_i) on a logarithmic scale ranging from 10^−6 to 10^0. Two sets of data are presented: one for a mildly convex profile (a_0.05) and another for a mildly concave profile (a_−0.05). The solid blue line represents the numerical calculation for a = 0.05, and the dashed blue line shows the short wave asymptotics for a = 0.05. Similarly, the solid yellow line represents the numerical calculation for a = −0.05, and the dashed yellow line indicates the long wave asymptotics for a = −0.05. The graph shows that for the convex profile, the growth rate increases with wavenumber up to a certain point and then decreases, following a c_i proportional to k^−4 trend. For the concave profile, the growth rate initially increases with wavenumber and then decreases, following a c_i proportional to k trend. All values are approximated.

Figure 7.

The variation of $c_i$ as a function of ${\textit{Re}}$ for four curvature parameters: $a = 3$ (blue curve), $a = 0$ (red curve), $a=-1$ (yellow curve) and $a = -3$ (purple curve). The black dashed lines overlapping the blue and purple solid lines represent the inviscid limits for $a = 3$ and $a = -3$ , respectively. The grey dashed curve that scales as $\textit{Re}^{-1/2}$ corresponds to the asymptotic solution for $a = 0$ , as derived by Miles (Reference Miles1960).

Figure 7. Long description

The line graph displays the variation of c_i as a function of Re for four different curvature parameters. The blue curve represents a equals 3, the red curve represents a equals 0, the yellow curve represents a equals negative 1, and the purple curve represents a equals negative 3. The black dashed lines overlapping the blue and purple solid lines indicate the inviscid limits for a equals 3 and a equals negative 3, respectively. The grey dashed curve, which scales as Re to the power of negative one half, corresponds to the asymptotic solution for a equals 0, as derived by Miles in 1960. All values are approximated.

Figure 8.

For $G = 0$ and $a = 3$ , the variation of $c_i$ as function of (a) Reynolds number ${\textit{Re}}$ , (b) wavenumber $k$ and (c) the scaled variable $kRe$ . In (a), three cases are considered corresponding to wavenumbers $0.1$ (blue), $0.5$ (red) and $1.5$ (yellow), respectively. The continuous curves are from numerical calculation and dotted lines indicate the corresponding inviscid $c_i$ values. In (b), three Reynolds numbers are considered: $\textit{Re} = 1$ (red continuous curve); $\textit{Re} = 10^5$ (blue continuous curve); inviscid (yellow dashed curve). In (c), four wavenumbers are considered: $k = 0.01$ (blue continuous curve), $k = 0.02$ (red dashed curve), $k = 0.05$ (yellow dotted curve) and $k = 0.1$ (purple dot dashed curve). The thin black lines depict the asymptotes at small and large $kRe$ .

Figure 8. Long description

The image contains three graphs. The first graph shows the variation of a parameter as a function of Reynolds number for three different wavenumbers: 0.1, 0.5, and 1.5. The continuous curves are from numerical calculations, and the dotted lines indicate the corresponding inviscid values. The second graph illustrates the variation of the parameter as a function of wavenumber for three different Reynolds numbers: 1, 100000, and inviscid. The third graph presents the variation of the parameter as a function of the scaled variable for four different wavenumbers: 0.01, 0.02, 0.05, and 0.1. The thin black lines depict the asymptotes at small and large scaled variables.

4.3. Growth rate behaviour: 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. Figure 9(a–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$ (figure 9 a,c,e,g) and $G=0.02$ (figure 9 b,d,f,h).

Figure 9.

Contours of the growth rate ( $k c_i$ ) in the $k$ – ${\textit{Re}}$ plane for four representative velocity profile curvatures: (a,b) $a = -3$ , (c,d) $a = -1$ , (e,f) $a = 0$ and (g,h) $a = 3$ . Panels (a), (c), (e) and (g) show the results for $G = 0$ limit (i.e. no gravity or surface tension), while panels (b), (d), (f) and (h) correspond to $G = 0.02$ . Panel (e) includes the neutral stability curves from Miles (Reference Miles1960), Smith & Davis (Reference Smith and Davis1982) and Miesen & Boersma (Reference Miesen and Boersma1995) for comparison. Stream function plots for points denoted by a cross ( $\times$ ) in (e) are shown in figure 11.

Figure 9. Long description

The image contains eight contour plots arranged in a 4_2 grid. Each plot shows the growth rate in the plane for different velocity profile curvatures. The x-axis represents the wavenumber (k) and the y-axis represents the Reynolds number (Re). The color bar on the right indicates the growth rate values, ranging from 10^−6 to 10^2. Panels (a), (c), (e), and (g) show results for the limit with no gravity or surface tension, while panels (b), (d), (f), and (h) correspond to a specific value of G. Panel (e) includes neutral stability curves from various studies for comparison. The plots illustrate how the growth rate varies with different parameters and conditions.

For the forward-bulging profile ( $a=-3$ ), two broad instability regions appear (figure 9 a,b), separated by a narrow band of stability. Three distinct modes can be identified within them. At small ${\textit{Re}}$ and $k$ , the long-wave interfacial (Yih) mode dominates, with peak growth rates at $\textit{Re}\sim O(10)$ and $k\sim O(1)$ . At large ${\textit{Re}}$ and moderate $k$ , a bright region corresponds to the rippling instability, consistent with figure 5(a), where the wavenumber of maximum growth at $\textit{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 ${\textit{Re}}$ , reminiscent of the shear mode described by Floryan, Davis & Kelly (Reference Floryan, Davis and Kelly1987). 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 ${\textit{Re}}$ and $k\sim O(1)$ (figure 9 c,d). Interestingly, at $G=0$ , the long-wave mode extends into the short-wave regime at high ${\textit{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. Figure 9(e,f) show that the unstable region coincides precisely with the neutral curves of earlier studies: the outer boundary from Miles (Reference Miles1960) (blue dots) and Smith & Davis (Reference Smith and Davis1982) (magenta dots), and the inner boundary from Miesen & Boersma (Reference Miesen and Boersma1995) (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 (Reference Boomkamp and Miesen1996)), with properties more akin to those of the shear instabilities in figure 9(b–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 (Reference Yiantsios and Higgins1988). At higher $G$ (e.g. $0.5$ ), this mode vanishes, restoring the expected stabilizing effect.

The flow-reversal profile ( $a=3$ ) displays contours (figure 9 g,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 ${\textit{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 $\textit{Re}\sim O(10^2)$ . Interestingly, at large ${\textit{Re}}$ and $k$ , a narrow band of stability appears within the rippling region, consistent with the deviations from inviscid predictions shown earlier in figure 8(b). 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 behaviour, with small positive and negative curvatures favouring short- and long-wave instabilities, respectively. More generally, viscosity is not merely dissipative: for long waves it can enhance growth compared with 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\lt -1$ and $a\gt 0$ ? These questions motivate a more detailed energy budget and eigenfunction analysis, which we turn to in the following section.

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 TAN to REY for long waves (figure 10 a). To explore this further, figure 11(a–c) present streamlines and normalized vorticity contours at $k=0.01$ for $\textit{Re}=1$ , $10^5$ and the inviscid limit, respectively.

At $\textit{Re}=1$ (figure 11 a), the vorticity is nearly uniform across the depth, with no appreciable horizontal displacement, and the value of the stream function decreases rapidly away from the surface, while being in-phase with both vorticity and surface displacement $\eta$ . This is characteristic of the shallow-viscous regime of Charru & Hinch (Reference Charru and Hinch2000) and agrees with the mechanism described by Smith (Reference Smith1990): in the long-wave, low- ${\textit{Re}}$ limit, the leading-order solution is a constant vorticity disturbance (see Appendix D). 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 $\textit{Re}=10^5$ case (figure 11 b) – which lies in the parameter region where growth rates exceed the inviscid values (figure 8 b) and where REY dominates over TAN without a critical layer (figure 10 a) – 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 (Reference Charru and Hinch2000), in which disturbances are confined to scales smaller than both the wavelength and the fluid depth. Finally, the inviscid solution (figure 11 c) reveals yet another distinct structure: a critical layer is now present (black dashed lines), unlike at $\textit{Re}=10^5$ , where vorticity simply peaked in the bulk. Outside this layer, vorticity decays rapidly, while the stream function resembles the $\textit{Re}=1$ case away from the critical layer but acquires a pronounced tilt in its vicinity, much like the behaviour shown in figure 5(d).

In summary, the $\textit{Re}=10^5$ eigenstructure is not merely a midpoint between the long-wave interfacial mode ( $\textit{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 10 a). 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. Figure 11 d–f) present streamlines and normalized vorticity contours for $k=20$ at $\textit{Re}=1$ , $10^5$ and the inviscid case, respectively.

At $\textit{Re}=1$ (figure 11 d), 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 stream function 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 (Reference Hinch1984): 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 (Reference Charru and Hinch2000), 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 (4.7).

At $\textit{Re}=10^5$ (figure 11 e), 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 (Reference Charru and Hinch2000), where penetration depths are smaller than a wavelength and comparable to the viscous length scale. In the inviscid case (figure 11 f), 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 11 c), the vorticity and stream function 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 ${\textit{Re}}$ , they correspond to the deep-viscous interfacial mode; at moderate to high ${\textit{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$ ). Figure 12(a–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-hand branches of the ‘V’-shaped instability region in figure 9(e), while the third lies at their intersection. The case $(k,Re)=(1,10^3)$ (figure 12 a) exhibits the hallmark features of a shear mode: vorticity concentrated near the bottom wall with a horizontal shift (Hooper & Boyd Reference Hooper and Boyd1987), the presence of a critical layer (Miles Reference Miles1960) and stream function contours filling most of the domain, consistent with the ‘A’-family of modes (Kaffel & Riaz Reference Kaffel and Riaz2015). The Reynolds stress $(-\overline {uw})$ shows a local extremum near the critical layer (Hooper Reference Hooper1989), and the perturbation KE is supplied entirely by REY, confirming that this instability is a shear mode and not the ‘internal mode’ suggested by Boomkamp & Miesen (Reference Boomkamp and Miesen1996).

Figure 12.

The streamlines of the disturbance flow field superimposed on contours of the normalized vorticity ( $\varOmega /\varOmega _{max}$ ) for $G = 0, a = 0$ , and for representative parameters given by $(k,Re) =$ (a) ( $1,10^3$ ), (b) ( $6, 10^3$ ) and (c) ( $1.75, 150$ ). These values are denoted by a cross ( $\times$ ) in figure 9(e). The black dashed lines indicate the critical layer.

Figure 12. Long description

The image displays streamlines of the disturbance flow field superimposed on contours of the normalized vorticity for different parameters. The parameters are represented by three subfigures: (a) with Reynolds number 1000 and wavenumber 1, (b) with Reynolds number 1000 and wavenumber 6, and (c) with Reynolds number 150 and wavenumber 1.75. Each subfigure shows the flow patterns and vorticity distribution. The black dashed lines indicate the critical layer in each subfigure.

Figure 13.

The streamlines of the disturbance flow field superimposed on the contours of the normalized disturbance vorticity ( $\varOmega /\varOmega _{max}$ ) for $ G = 0$ and $ Re = 10^{4}$ . Panels (a), (b) and (c) correspond to $k = 0.1$ , $k=1$ and $k=10$ , respectively, with $a=3$ . Panels (d), (e) and (f) correspond to $k = 0.1$ , $k=1$ and $k=10$ , respectively, with $a=-1$ .

Figure 13. Long description

The image contains six panels of contour plots showing streamlines and disturbance vorticity for different values of a and k. Each panel represents a different combination of these parameters. Panels (a), (b), and (c) correspond to a = 3 with k values of 0.1, 1, and 10, respectively. Panels (d), (e), and (f) correspond to a = −1 with the same k values. The color scale on the right indicates the normalized disturbance vorticity, ranging from −1 to 1. The streamlines are depicted with arrows indicating the direction of flow. The plots illustrate how the disturbance flow field and vorticity change with varying values of a and k.

For $(k,Re)=(6,10^3)$ (figure 12 b), 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 (Reference Hooper and Boyd1987) 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 12 c) 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 10(c), confirming that it too corresponds to a shear mode, albeit without the right branch. In general, both left- and right-hand branches occur for curvature values close to $a=0$ , but their extent diminishes differently as $a$ moves away from zero, leaving only the left-hand branch at larger $a$ (as in figure 10 c).

5.4. A composite mode at large ${\textit{Re}}$

A deviation between the growth rates at $\textit{Re} = 10^5$ and the inviscid case was attributed in figure 8(b) to the large- ${\textit{Re}}$ asymptote at long waves. Hence, the large- ${\textit{Re}}$ , small- $k$ region does not correspond to a pure rippling instability. Figure 10(a) 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 behaviour, figure 13(a–c) present eigenfunctions at $\textit{Re}=10^4$ for $k=0.1,1$ and $10$ . At $k=0.1$ (figure 13 a), the vorticity and stream function resemble a long-wave interfacial mode but with a horizontal tilt at depth, similar to figure 11(b) despite the order-of-magnitude difference in parameters. Increasing the wavenumber to $k=1$ (figure 13 b) 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 13 c), the eigenstructure closely resembles figure 11(e): vorticity peaks at the surface and also has a local maximum near a surface-critical layer, combining surface-driven shear behaviour (cf. figure 12 b) with rippling features (cf. figure 11 f).

A comparison with the $a=-1$ case (figure 13 d–f), which lacks a rippling instability, highlights the difference. Although the eigenfunctions deviate from the classical long-wave structure of figure 11(a), the variations in tilt and penetration depth are monotonic with $k$ and ${\textit{Re}}$ , consistent with inertia-modified long-wave interfacial instability (Charru & Hinch Reference Charru and Hinch2000). 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 attributable to a coalescence of interfacial and shear modes (Özgen Reference Özgen2008), nor to a coalescence of two shear modes (Tilley, Davis & Bankoff Reference Tilley, Davis and Bankoff1994). It also differs from the mode merging reported by Ishimura et al. (Reference Ishimura, Mergui, Ruyer-Quil and Dietze2023, Reference Ishimura, Mergui, Ruyer-Quil and Dietze2025) between a downstream-propagating Kapitza mode and an upstream-propagating short-wave mode, since both the long-wave and rippling instabilities in the present study propagate in the same direction and are never found to coexist as separate branches at any point in the parameter space. Furthermore, the composite mode is distinct from the internal mode of Miesen & Boersma (Reference Miesen and Boersma1995) and Schmidt et al. (Reference Schmidt, Ó Náraigh, Lucquiaud and Valluri2016), where tangential stress dominates the perturbation KE with only a small Reynolds stress contribution. Rather, the composite mode is a single branch whose character varies continuously across the parameter space, with the dominant energy source shifting from tangential stress to Reynolds stress as ${\textit{Re}}$ increases. Similar dual behaviour was noted only by Albert & Charru (Reference Albert and Charru2000) in two-phase Couette flows, where an interfacial mode at small ${\textit{Re}}$ behaved like a shear mode at large ${\textit{Re}}$ . The richer composite nature described here, however, does not appear to have been reported previously.