2.1 Vorticity equation

We model the density stratified shear layer using the two-dimensional (2-D) incompressible and inviscid Navier–Stokes equations along with the mass continuity equation. The flow is in the $x{-}z$ plane with the density stratification along the vertical ( $z$ ) axis. The horizontal and vertical components of velocity are respectively $u$ and $w$ . To accurately capture the effects of the free surface we use a non-Boussinesq model. The model takes into account the effect of density inhomogeneity in the first order of the inertial terms, effects on higher-order nonlinear terms are neglected. A good way to understand the effect of density inhomogeneity is to examine the vorticity evolution equation, which for a 2-D, inviscid, density stratified fluid is given by

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}q}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}q}{\unicode[STIX]{x2202}x}+w\frac{\unicode[STIX]{x2202}q}{\unicode[STIX]{x2202}z}=\frac{1}{\unicode[STIX]{x1D70C}^{2}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}z}\right). & \displaystyle\end{eqnarray}$$

Here vorticity is given by $q\equiv \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z-\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x$ , while $\unicode[STIX]{x1D70C}$ , $p$ and $g$ denote density, pressure and acceleration due to gravity respectively. The terms on the right-hand side of (2.1) denote the baroclinic generation of vorticity. A base state that varies only along the $z$ axis is assumed, and is given by $u=\bar{u}(z)$ , $w=0$ , $\bar{q}(z)=\text{d}\bar{u}/\text{d}z$ , $p=\bar{p}(z)$ and $\unicode[STIX]{x1D70C}=\bar{\unicode[STIX]{x1D70C}}(z)$ . The base state follows hydrostatic pressure balance $\text{d}\bar{p}/\text{d}z=-\bar{\unicode[STIX]{x1D70C}}g$ . Perturbations are added to the base flow: $u=\bar{u}(z)+\tilde{u}$ , $w=\tilde{w}$ , $p=\bar{p}(z)+\tilde{p}$ , $\unicode[STIX]{x1D70C}=\bar{\unicode[STIX]{x1D70C}}(z)+\tilde{\unicode[STIX]{x1D70C}}$ and $q=\bar{q}(z)+\tilde{q}$ , where $\tilde{f}$ denotes the perturbation quantities ( $f$ is a placeholder variable). We assume the perturbations to be infinitesimal and linearize (2.1). This results in the perturbation vorticity evolution equation

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{q}}{\unicode[STIX]{x2202}t}+\bar{u}\frac{\unicode[STIX]{x2202}\tilde{q}}{\unicode[STIX]{x2202}x}=\underbrace{-\tilde{w}\frac{\text{d}\bar{q}}{\text{d}z}}_{\text{barotropic}}+\underbrace{\frac{g}{\bar{\unicode[STIX]{x1D70C}}}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x}}_{\text{gravitational baroclinic}}+\underbrace{\frac{1}{\bar{\unicode[STIX]{x1D70C}}^{2}}\left(\frac{\text{d}\bar{\unicode[STIX]{x1D70C}}}{\text{d}z}\frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}x}\right)}_{T_{3}\equiv \text{non-Boussinesq baroclinic}}. & \displaystyle\end{eqnarray}$$

Equation (2.2) can also be found in Heifetz & Mak (Reference Heifetz and Mak2015) (their equation (5)). The right-hand side of the above equation provides different sources of vorticity generation. The barotropic generation term arises due to the advection of base state vorticity gradient by the perturbation vertical velocity. The ‘gravitational baroclinic torque’ is responsible for the propagation of interfacial gravity waves at the pycnocline and it is the only baroclinic generation term present when the Boussinesq approximation is invoked. In the non-Boussinesq (large density stratification) regime, in addition to the gravitational baroclinic generation term, there is the ‘non-Boussinesq baroclinic generation term’ (hereafter referred to as $T_{3}$ ), that arises out of the density variations in the inertial terms, and is completely independent of the gravitational effects. However, even for modelling the free surface (the highly non-Boussinesq interface between air and water), $T_{3}$ has been neglected in Bakas & Ioannou (Reference Bakas and Ioannou2009) (see their equation (6)) as well as Longuet-Higgins (Reference Longuet-Higgins1998) (they do not explicitly write in terms of the vorticity equation). In § 3.2 we provide a detailed discussion so as to delineate the circumstances under which $T_{3}$ can or cannot be ignored while modelling interfaces with large density jumps (like the free surface).

2.2 The non-Boussinesq Taylor–Goldstein equation

We now present the non-Boussinesq Taylor–Goldstein equation in the inviscid and non-diffusive limit. To obtain it we assume temporal normal mode form for the perturbations given by $\tilde{f}=\hat{f}(z)\text{e}^{\text{i}\unicode[STIX]{x1D6FC}(x-ct)}$ . Here $\unicode[STIX]{x1D6FC}$ and $c$ are respectively the real wavenumber and the complex phase speed ( $c=c_{r}+\text{i}c_{i}$ ), and $f$ could represent $u$ , $w$ , $p$ , $q$ or $\unicode[STIX]{x1D70C}$ . Temporal normal mode form is substituted in the linearized Navier–Stokes equations, yielding

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \underbrace{\bar{\unicode[STIX]{x1D70C}}^{\prime }[(\bar{u}-c){\hat{w}}^{\prime }-\bar{u}^{\prime }{\hat{w}}]}_{\text{non-Boussinesq}}-\underbrace{\frac{g\bar{\unicode[STIX]{x1D70C}}^{\prime }}{\bar{u}-c}{\hat{w}}}_{\text{gravitational}}+\bar{\unicode[STIX]{x1D70C}}[(\bar{u}-c)({\hat{w}}^{\prime \prime }-\unicode[STIX]{x1D6FC}^{2}{\hat{w}})]-\underbrace{\bar{\unicode[STIX]{x1D70C}}\bar{u}^{\prime \prime }{\hat{w}}}_{\text{barotropic}}=0. & \displaystyle\end{eqnarray}$$

The total derivative with respect to $z$ is denoted by $^{\prime }$ . The non-Boussinesq Taylor–Goldstein equation obtained is same as the one obtained by Barros & Choi (Reference Barros and Choi2011), Barros & Choi (Reference Barros and Choi2014) and Carpenter, Guha & Heifetz (Reference Carpenter, Guha and Heifetz2017). Terms denoted by the braces indicate the corresponding terms in the perturbation vorticity equation (2.2). We intend to solve (2.3) for smooth profiles of base state velocity and density, and for this a numerical eigenvalue solver is necessary. The numerical solution of (2.3) poses a resolution issue. A very fine spacing in parameter space is required which raises the computational cost tremendously. An efficient way to control the resolution issue is to introduce numerical viscosity ( $\unicode[STIX]{x1D707}$ ) and diffusivity ( $\unicode[STIX]{x1D705}$ ). Therefore we numerically solve the viscous diffusive form of the non-Boussinesq Taylor–Goldstein equation, which is given by

(2.4a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \bar{\unicode[STIX]{x1D70C}}^{\prime }[(\bar{u}-c){\hat{w}}^{\prime }-\bar{u}^{\prime }{\hat{w}}]+\bar{\unicode[STIX]{x1D70C}}[(\bar{u}-c)({\hat{w}}^{\prime \prime }-\unicode[STIX]{x1D6FC}^{2}{\hat{w}})-\bar{u}^{\prime \prime }{\hat{w}}]\nonumber\\ \displaystyle & & \displaystyle \quad =-\text{i}\unicode[STIX]{x1D6FC}\hat{\unicode[STIX]{x1D70C}}g-\frac{\text{i}}{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D707}[{\hat{w}}^{\prime \prime \prime \prime }-2\unicode[STIX]{x1D6FC}^{2}{\hat{w}}^{\prime \prime }+\unicode[STIX]{x1D6FC}^{4}{\hat{w}}],\end{eqnarray}$$

(2.4b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \text{i}\unicode[STIX]{x1D6FC}(\bar{u}-c)\hat{\unicode[STIX]{x1D70C}}+{\hat{w}}\bar{\unicode[STIX]{x1D70C}}^{\prime }=\unicode[STIX]{x1D705}(\hat{\unicode[STIX]{x1D70C}}^{\prime \prime }-\unicode[STIX]{x1D6FC}^{2}\hat{\unicode[STIX]{x1D70C}}).\end{eqnarray}$$

2.3 Solution of the eigenvalue problem

Equation set (2.4a,b ) can be converted into a generalized eigenvalue problem of the form $M\unicode[STIX]{x1D751}=cN\unicode[STIX]{x1D751}$ :

(2.5) $$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}cc@{}}M_{11} & M_{12}\\ M_{21} & M_{22}\end{array}\right]\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D717}_{1}\\ \unicode[STIX]{x1D717}_{2}\end{array}\right]=c\left[\begin{array}{@{}cc@{}}N_{11} & N_{12}\\ N_{21} & N_{22}\end{array}\right]\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D717}_{1}\\ \unicode[STIX]{x1D717}_{2}\end{array}\right]. & & \displaystyle\end{eqnarray}$$

Elements of the matrices are given by

(2.6) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}M_{11}=\bar{\unicode[STIX]{x1D70C}}^{\prime }[-\bar{u}D+\bar{u}^{\prime }]+\bar{\unicode[STIX]{x1D70C}}[-\bar{u}D^{2}+\unicode[STIX]{x1D6FC}^{2}\bar{u}+\bar{u}^{\prime \prime }]-{\displaystyle \frac{\text{i}}{\unicode[STIX]{x1D6FC}}}\unicode[STIX]{x1D707}[D^{4}-2\unicode[STIX]{x1D6FC}^{2}D^{2}+\unicode[STIX]{x1D6FC}^{4}],\\ M_{12}=-\text{i}\unicode[STIX]{x1D6FC}g,\quad M_{21}=\bar{\unicode[STIX]{x1D70C}}^{\prime },\quad M_{22}=\text{i}\unicode[STIX]{x1D6FC}\bar{u}-\unicode[STIX]{x1D705}[D^{2}-\unicode[STIX]{x1D6FC}^{2}],\quad \unicode[STIX]{x1D717}_{1}={\hat{w}},\unicode[STIX]{x1D717}_{2}=\hat{\unicode[STIX]{x1D70C}},\\ N_{11}=-\bar{\unicode[STIX]{x1D70C}}^{\prime }D+\bar{\unicode[STIX]{x1D70C}}[-D^{2}+\unicode[STIX]{x1D6FC}^{2}],\quad N_{12}=0,N_{21}=0,N_{22}=\text{i}\unicode[STIX]{x1D6FC}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Equation (2.5) is an eigenvalue problem for eigenvalue $c$ and eigenfunctions ${\hat{w}}$ and $\hat{\unicode[STIX]{x1D70C}}$ . The first, second and fourth derivative matrices of the total derivatives with respect to $z$ are given by $D$ , $D^{2}$ and $D^{4}$ respectively, which have been numerically evaluated using the fourth-order central difference scheme. Since central differencing is not possible at the boundaries, the boundary points have been discretized using second-order one sided finite difference scheme. The boundary conditions used for ${\hat{w}}$ are impenetrability and free slip. The impenetrable boundary condition arises due to the continuum hypothesis, and is given by

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\hat{w}}=0. & \displaystyle\end{eqnarray}$$

Furthermore, the free-slip boundary condition,

(2.8) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\text{d}\hat{u} }{\text{d}z}=0, & & \displaystyle\end{eqnarray}$$

along with the fluid incompressibility give

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}{\hat{w}}}{\text{d}z^{2}}=0. & \displaystyle\end{eqnarray}$$

We use an insulating boundary condition for $\hat{\unicode[STIX]{x1D70C}}$ , which is given by

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\hat{\unicode[STIX]{x1D70C}}}{\text{d}z}=0. & \displaystyle\end{eqnarray}$$

The insulating boundary condition is preferable in cases where there is no physical boundary, since it has minimal effect on the flow. For the upper boundary, we have assumed an ‘imaginary rigid lid’ in the air region at a certain height (chosen such that it does not adversely impact the results) above the free surface. Boundary conditions (2.7)–(2.10) have been applied to the imaginary rigid lid.

The matrix eigenvalue problem is solved using built-in functions in MATLAB. The numerical solution procedure used is similar to the one used by Smyth et al. (Reference Smyth, Moum and Nash2011). The numerical routine developed has been validated against the Boussinesq shear layer problem of Smyth et al. (Reference Smyth, Klaassen and Peltier1988) and that of surface gravity waves. Our primary interest is to capture inviscid instabilities like KH and Holmboe, hence viscosity ( $\unicode[STIX]{x1D707}$ ) and diffusivity ( $\unicode[STIX]{x1D705}$ ) introduced in (2.4a,b ) are treated as purely numerical parameters for controlling the resolution issue. The stability characteristics do not vary appreciably on changing the numerical parameters by an order of magnitude.

3.1 Base state profiles

In order to model a stratified shear layer, we use hyperbolic tangent functions in $z$ to represent the base state velocity and density profiles. The first point ensures that the derivatives of base state velocity remain finite, which is required for the numerical stability of (2.5). The second point ensures avoiding the formation of critical layers in the air region. In other words, we avoid the wind–wave instability of Miles (Reference Miles1957). The dimensional base state velocity variation is given by

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{u}(z)=U\tanh \left(\frac{z-\frac{1}{2}H}{h}\right). & \displaystyle\end{eqnarray}$$

Here $U$ (in $\text{m}~\text{s}^{-1}$ ) represents the surface current and $h$ is the half shear layer thickness. $H$ is the dimensional half-depth of the channel (i.e. the distance between free surface and bottom is $2H$ ). In this paper we have taken $h=1/2$ m. The base state velocity is depicted in figure 2(a), while the base state vorticity is depicted in figure 2(b). If $H/h$ is not large, moderately long waves inside the shear layer can feel the effect of the free surface. In the present study we have detailed the results for $H/h=10,\,4,\,3$ and $2$ . These values of $H/h$ have been chosen so that a clear transition in stability characteristics can be seen as the surface gravity waves come closer to the other waves present in the shear layer (therefore gaining an opportunity to resonantly interact with them). The dimensional base state density profile is given by

(3.2) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}(z)=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{01}\left(1-A_{t,\,0}\tanh \left[5\left({\displaystyle \frac{z-H}{h}}\right)\right]\right) & 1.2H\geqslant z\geqslant 0.75H,\\[12.0pt] \unicode[STIX]{x1D70C}_{12}\left(1-A_{t}\tanh \left[5\left({\displaystyle \frac{z-\frac{1}{2}H}{h}}\right)\right]\right) & 0.75H\geqslant z\geqslant 0.\end{array}\right. & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D70C}_{01}=(\unicode[STIX]{x1D70C}_{ww}+\unicode[STIX]{x1D70C}_{air})/2=491~\text{kg}~\text{m}^{-3}$ is the mean density of air and warm water, while the Atwood number corresponding to the air–warm water interface (i.e. free surface) is given by $A_{t,\,0}=(\unicode[STIX]{x1D70C}_{ww}-\unicode[STIX]{x1D70C}_{air})/(\unicode[STIX]{x1D70C}_{ww}+\unicode[STIX]{x1D70C}_{air})=9.97\times 10^{-1}\approx 1$ . Similarly the mean density for cold water and warm water is given by $\unicode[STIX]{x1D70C}_{12}=(\unicode[STIX]{x1D70C}_{cw}+\unicode[STIX]{x1D70C}_{ww})/2=988.35~\text{kg}~\text{m}^{-3}$ and the Atwood number corresponding to the pycnocline is given by

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle A_{t}=\frac{\unicode[STIX]{x1D70C}_{cw}-\unicode[STIX]{x1D70C}_{ww}}{\unicode[STIX]{x1D70C}_{cw}+\unicode[STIX]{x1D70C}_{ww}}=\frac{10.3~\text{kg}~\text{m}^{-3}}{(993.5+983.2)\,\text{kg}~\text{m}^{-3}}\approx 0.005. & \displaystyle\end{eqnarray}$$

The densities of cold water and warm water are respectively given by $\unicode[STIX]{x1D70C}_{cw}$ and $\unicode[STIX]{x1D70C}_{ww}$ .

Figure 2.

The base state plots for (a) horizontal velocity, $\bar{u}$ in $\text{m}~\text{s}^{-1}$ , (b) corresponding vorticity, $\bar{q}(z)=\text{d}\bar{u}/\text{d}z$ in $\text{s}^{-1}$ , and (c) density, $\bar{\unicode[STIX]{x1D70C}}$ in $\text{kg}~\text{m}^{-3}$ corresponding to $H/h=10$ (which indicates the location of the free surface). The horizontal lines indicate the domain extent of the water body for each $H/h$ . For example, when $H/h=3$ , the vertical domain of the water body is confined between the horizontal lines indicated by ‘ $3$ ’ (which means that the free surface and the bottom boundary will respectively occur at the upper and lower horizontal lines indicated by ‘ $3$ ’).

Base state density has the units of $\text{kg}~\text{m}^{-3}$ , and is plotted in figure 2(c).

A key parameter of importance is the bulk Richardson number, $J$ , given by

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle J=2A_{t}\frac{gh}{U^{2}}. & \displaystyle\end{eqnarray}$$

For reporting stability characteristics we have non-dimensionalized the wavenumber, $\unicode[STIX]{x1D6FC}$ , the phase speed, $c$ , and the growth rate $\text{Im}\{\unicode[STIX]{x1D6FC}c\}$ . The wavenumber is non-dimensionalized by the shear layer half-thickness, $h$ , while the phase speed is non-dimensionalized by the surface velocity, $U$ . Growth rates are non-dimensionalized by the shear scale $U/h$ .

3.2 Is ignoring the non-Boussinesq baroclinic term legitimate for modelling the air–water free surface?

As noted in § 2.1, a couple of previous studies (Longuet-Higgins Reference Longuet-Higgins1998; Bakas & Ioannou Reference Bakas and Ioannou2009) that considered the effect of the free surface on the shear instability underneath have ignored the term $T_{3}\equiv (\bar{\unicode[STIX]{x1D70C}})^{-2}(\text{d}\bar{\unicode[STIX]{x1D70C}}/\text{d}z)(\unicode[STIX]{x2202}\tilde{p}/\unicode[STIX]{x2202}x)$ appearing in (2.2). For the purpose of comparison we have analysed the free surface effect using two different models, namely the ‘Boussinesq free surface model’ (which ignores $T_{3}$ ) and the ‘non-Boussinesq free surface model’ (which fully solves (2.2)). Here we try to understand the implications of assuming $T_{3}=0$ . There are two possible interpretations of $T_{3}=0$ .

(i) $\text{d}\bar{\unicode[STIX]{x1D70C}}/\text{d}z=0$ . This is a valid approximation in the Boussinesq limit (where $\bar{\unicode[STIX]{x1D70C}}\approx \text{const.}$ ). In our setting, $T_{3}$ will only matter at the free surface, where obviously $\text{d}\bar{\unicode[STIX]{x1D70C}}/\text{d}z\neq 0$ , in fact approximately a delta function. Hence assuming $\text{d}\bar{\unicode[STIX]{x1D70C}}/\text{d}z=0$ for the air–water density jump is not correct.

(ii) $\unicode[STIX]{x2202}\tilde{p}/\unicode[STIX]{x2202}x=0$ . In the inviscid limit

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}x}=-\bar{\unicode[STIX]{x1D70C}}\left[\frac{\mathbb{D}\tilde{u} }{\mathbb{D}t}+\tilde{w}\bar{q}\right], & \displaystyle\end{eqnarray}$$

where $\bar{q}=\text{d}\bar{u}/\text{d}z$ and $\mathbb{D}\tilde{u} /\mathbb{D}t\equiv \unicode[STIX]{x2202}\tilde{u} /\unicode[STIX]{x2202}t+\bar{u}\unicode[STIX]{x2202}\tilde{u} /\unicode[STIX]{x2202}x$ is the linearized material derivative of $\tilde{u}$ ; see (A 8). For surface gravity waves, the free surface is a vortex sheet ( $\tilde{u}$ changes sign above and below the free surface). To elaborate this point, consider a two-fluid system, an upper fluid with density $\unicode[STIX]{x1D70C}_{1}$ and lower fluid with density $\unicode[STIX]{x1D70C}_{2}$ . For analytical simplicity, both fluids are assumed to be infinitely deep. In the absence of any background flow, the velocity just below (above) the density interface is

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{\pm }=\pm \unicode[STIX]{x1D714}\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x-\unicode[STIX]{x1D714}t)}, & \displaystyle\end{eqnarray}$$

where $+(-)$ indicates below (above), and $\unicode[STIX]{x1D714}=\sqrt{g\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})/(\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2})}$ is the frequency; see Kundu, Cohen & Dowling (Reference Kundu, Cohen and Dowling2012, chap. 7.7). The velocity at the interface is the average of the two: $\tilde{u} =(\tilde{u} _{+}+\tilde{u} _{-})/2=0$ . Substitution of $\tilde{u} =0$ in (3.5) yields $\unicode[STIX]{x2202}\tilde{p}/\unicode[STIX]{x2202}x=0$ , provided $\bar{q}=0$ . Hence $T_{3}=0$ . Proceeding further and using the normal mode ansatz, we obtain the dispersion relation for surface gravity waves. We note here in passing that such intricacies do not appear in the standard technique for deriving the surface gravity wave dispersion relation, which simplifies the problem from the beginning by using potential flow approximation in each layer, and then makes use of kinematic and dynamic boundary conditions; see Kundu et al. (Reference Kundu, Cohen and Dowling2012, chap. 7.2). However in a generalized scenario where background shear is present at the density interface (i.e.  $\bar{q}\neq 0$ ), equation (3.5) directly reveals that $\unicode[STIX]{x2202}\tilde{p}/\unicode[STIX]{x2202}x=-\bar{\unicode[STIX]{x1D70C}}\bar{q}\tilde{w}\neq 0$ , implying $T_{3}\neq 0$ . In summary, $T_{3}=0$ , as implicitly demanded by the Boussinesq free surface model, is not correct when background shear is present at the free surface.

Figure 3.

Growth rate contours in the $\unicode[STIX]{x1D6FC}{-}J$ plane for the numerical stability of smooth profiles. Each row represents growth rates for a given $H/h$ . (a–c) $H/h=10$ , (d–f) $H/h=4$ , (g–i) $H/h=3$ and (j–l) $H/h=2$ . (a,d,g,j), (b,e,h,k) and (c,f,i,l) respectively correspond to rigid lid, Boussinesq free surface and non-Boussinesq free surface.

Table 1.

Growth rates and phase speeds for the maximum growth rate mode of the smooth profiles. ‘H’ stands for Holmboe. To keep parity in nomenclature, we have referred to even the remnants of KH (in Boussinesq free surface as well as non-Boussinesq free surface cases) as ‘KH’. These modes are strongly affected by the surface gravity waves.

3.3 Variation of stability characteristics with $H/h$

Here we have studied the effect of varying the distance between the free surface and the pycnocline (i.e.  $H/h$ ) on the submerged stratified shear layer. As mentioned previously, the proximity between the surface gravity waves and the vorticity and interfacial gravity waves existing at the shear layer can affect the stability characteristics. A comparison of stability characteristics of a stratified shear layer in the presence of a rigid lid, a Boussinesq free surface and an actual/non-Boussinesq free surface for $H/h=10$ are shown in figure 3(a–c). In this case we see a very good match between the rigid lid, the Boussinesq free surface and the non-Boussinesq free surface, which is expected since the free surface is located quite far away. The maximum growth rate and the corresponding phase speed also agree very well, as can be seen from table 1.

If we look at the first column of figure 3, which corresponds to the rigid-lid case, it becomes clear that decreasing $H/h$ significantly decreases the maximum growth rate. The KH instabilities (the downward pointing closed curves having higher growth rates) are more significantly affected than the Holmboe modes (the upward pointing open curves with lesser growth rates). Moreover, short waves are stabilized by this process, as is evident from the leftward shift of the right stability boundary. The results are in good agreement with those of Hazel (Reference Hazel1972) and Haigh & Lawrence (Reference Haigh and Lawrence1999).

The results obtained on replacing the rigid lid by a non-Boussinesq free surface are significantly different; see the last column of figure 3. Decreasing $H/h$ has little effect on the maximum growth rate; see table 1. However, the stability boundaries change quite dramatically, and furthermore, new modes appear. The distinct identity of the KH mode observed for $H/h=10$ is lost in $H/h=3$ ; a part of it separates out as a ‘bubble’ and becomes part of the Holmboe branch. The other part near $J=0$ axis leads to a new mode, the ‘SG-LV mode’, which is distinctly visible for $H/h=3$ and $2$ . The SG-LV mode always has the highest growth rate. Comparing with the rigid-lid case (first column), it can be straightforwardly argued that the SG-LV mode would not exist if there were no surface gravity waves. We will show later that the SG-LV mode is a result of the interaction between waves $2$ and $6$ in figure 1.

3.4 Eigenfunction analysis of a few important modes

Eigenfunctions are useful in understanding the physical structure of the modes. The eigenfunction of the perturbation streamfunction, $\hat{\unicode[STIX]{x1D713}}$ , can be obtained from the vertical velocity ${\hat{w}}$ via the relation $\hat{\unicode[STIX]{x1D713}}=-\text{i}{\hat{w}}/\unicode[STIX]{x1D6FC}$ ( ${\hat{w}}$ is obtained by solving the eigenvalue problem (2.5)). Since eigenfunctions are non-unique, we normalize $\hat{\unicode[STIX]{x1D713}}$ by its maximum value for reporting purposes. Apart from $\hat{\unicode[STIX]{x1D713}}$ , we have also plotted the contours of perturbation streamfunction $\tilde{\unicode[STIX]{x1D713}}=\text{Re}\{\hat{\unicode[STIX]{x1D713}}\text{e}^{\text{i}\unicode[STIX]{x1D6FC}x}\}$ for one wavelength of the disturbance. The direction of the perturbed velocity field would be tangential to the contour lines. The eigenfunction of the density perturbations $\hat{\unicode[STIX]{x1D70C}}$ is also of significant interest. For plotting purposes, we have normalized it by the maximum value of the corresponding $\hat{\unicode[STIX]{x1D713}}$ . The contours of perturbation density $\tilde{\unicode[STIX]{x1D70C}}$ are plotted for one wavelength. This spatial variation reveals information about the gravity waves present in the system, their relative magnitudes and phases. Here we have chosen to analyse two specific unstable modes, the Holmboe mode for $H/h=10$ and the SG-LV mode for $H/h=3$ .

3.4.1 Unstable Holmboe mode for $H/h=10$

Holmboe instability results due to the interaction between counter-propagating vorticity waves and interfacial gravity waves. The maximum growth rate mode for Holmboe instability corresponds to $(\unicode[STIX]{x1D6FC},J)=(0.6775,0.3498)$ in figure 3(c). The corresponding streamfunction and density perturbation contours and eigenfunctions have been plotted in figure 4. The peaks of $|\hat{\unicode[STIX]{x1D713}}|$ occur at the pycnocline (upper peak) and the lower vorticity gradient extremum (lower peak), verifying the fact that Holmboe instability arises due to the interaction between an interfacial gravity wave and a vorticity wave. We should note here that the free surface introduces a small asymmetry, favouring the ‘leftward Holmboe mode’ over the ‘rightward Holmboe mode’. This is because the leftward propagating vorticity wave (which exists at the upper extrema of base vorticity gradient), being closer to the free surface, is more affected. Therefore, the leftward propagating Holmboe mode, arising due to the interaction between the rightward propagating vorticity wave and the leftward propagating interfacial gravity wave, becomes the dominant mode of the two Holmboe modes. This leftward propagating mode, with phase speed $c_{r}=-0.4155$ (see table 1), curves into the lower fluid region, as clearly evidenced in figure 4(f).

Figure 4.

Perturbation streamfunction and density for the non-Boussinesq free surface case for $H/h=10$ . The plot corresponds to the maximum growth rate mode of the Holmboe instability, given by $(\unicode[STIX]{x1D6FC},J)=(0.6775,0.3498)$ . (a) Norm of the perturbation streamfunction eigenfunction $\hat{\unicode[STIX]{x1D713}}$ versus $z$ , and (b) contours of the perturbation streamfunction $\tilde{\unicode[STIX]{x1D713}}$ . The perturbation density characteristics near the free surface ( $z=10$ ): (c) norm of the perturbation density eigenfunction $\hat{\unicode[STIX]{x1D70C}}$ versus $z$ , and (d) contours of the perturbation density $\tilde{\unicode[STIX]{x1D70C}}$ . The perturbation density characteristics near the pycnocline ( $z=5$ ): (e) norm of the perturbation density eigenfunction $\hat{\unicode[STIX]{x1D70C}}$ versus $z$ , and (f) contours of the perturbation density  $\tilde{\unicode[STIX]{x1D70C}}$ .

A few salient features of Holmboe instability in the presence of a free surface are worth noticing. The value of $\hat{\unicode[STIX]{x1D70C}}$ is a few orders of magnitude smaller at the free surface than that at the interface; see figure 4(c–f). The same figures also show that the density contours at the free surface and the interface are $\unicode[STIX]{x03C0}$ phase-shifted. These two features are the hallmark of the classic ‘baroclinic/internal/varicose mode’ in layered flows in the absence of background velocity shear; see Sutherland (Reference Sutherland2010, chap. 2). The analogy of the Holmboe instability with the baroclinic mode seems appropriate because, like the baroclinic mode, the Holmboe instability is driven by the baroclinic torque produced at the interface.

3.4.2 Unstable SG-LV mode for $H/h=3$

The maximum SG-LV mode growth rate for $H/h=3$ is given by $(\unicode[STIX]{x1D6FC},J)=(0.4808,0)$ (refer to figure 3 i). Eigenfunctions of the perturbation streamfunction and density are plotted in figure 5. Figure 5(c–f) shows that $|\hat{\unicode[STIX]{x1D70C}}|$ has two peaks, one at the free surface and the other at the pycnocline, the former being two orders of magnitude greater than the latter. Furthermore, the surface and the interfacial gravity waves are nearly in phase. This configuration resembles the classic ‘barotropic/external/sinuous mode’ in layered flows; see Sutherland (Reference Sutherland2010, chap. 2). Although from the eigenfunction plot the mechanism behind the SG-LV instability is not very clear, yet a very important conclusion can be drawn: the interfacial gravity wave plays nearly no role in the instability process. In other words, the instability has to be due to the interaction between the vorticity waves present in the shear layer and the surface gravity waves.

As is evident from table 1, the highest growth rate of the SG-LV mode for $H/h=3$ is comparable to the KH mode in a nearly unbounded flow (compare with the growth rate of the KH mode for $H/h=10$ ). The phase speed of this mode is ${\approx}-0.1$ , and is intermediate between the KH and Holmboe modes for the same $H/h$ .

Figure 5.

As in figure 4 but for $H/h=3$ . The plot corresponds to the maximum growth rate mode of the SG-LV instability, given by $(\unicode[STIX]{x1D6FC},J)=(0.4808,0)$ .

4.1 Broken-line profiles and reduced-order models

To complement the stability analyses of smooth profiles outlines in § 3, we undertake stability analyses of broken-line profiles in this section. In these profiles, vorticity (vertical gradient of velocity) and density are piecewise constant while velocity is piecewise linear. Broken-line profiles are useful in identifying the mechanisms behind different instabilities occurring in the flow since they allow the waves in the system to become localized. Furthermore, using such profiles, (2.3) can be solved analytically. An equivalent form of (2.3) is given by

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \{\bar{\unicode[STIX]{x1D70C}}[(\bar{u}-c){\hat{w}}^{\prime }-\bar{u}^{\prime }{\hat{w}}]\}^{\prime }-\frac{\bar{\unicode[STIX]{x1D70C}}^{\prime }g}{\bar{u}-c}{\hat{w}}-\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D6FC}^{2}(\bar{u}-c){\hat{w}}=0. & \displaystyle\end{eqnarray}$$

One or more waves can exist at a material interface located arbitrarily at $z=z_{i}$ . A material interface arises due to the discontinuity in the base state vorticity and/or base state density. In each layer (fluid between two consecutive interfaces), both (2.3) and (4.1) are reduced to $\bar{\unicode[STIX]{x1D70C}}(\bar{u}-c)[{\hat{w}}^{\prime \prime }-\unicode[STIX]{x1D6FC}^{2}{\hat{w}}]=0$ . The discrete eigenspectrum can be obtained by assuming $(\bar{u}-c)\neq 0$ , yielding

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\hat{w}}^{\prime \prime }-\unicode[STIX]{x1D6FC}^{2}{\hat{w}}=0. & \displaystyle\end{eqnarray}$$

The continuity of vertical velocity across an interface is imposed:

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x27E6}\,{\hat{w}}\,\unicode[STIX]{x27E7}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x27E6}....\unicode[STIX]{x27E7}$ denotes the difference across the interface. Continuity of pressure across the interface is ensured by the dynamic condition, obtained by integrating (4.1) across the interface from $z_{i}-\unicode[STIX]{x0394}z$ to $z_{i}+\unicode[STIX]{x0394}z$ and letting $\unicode[STIX]{x0394}z\rightarrow 0$ ,

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\unicode[STIX]{x27E6}\bar{\unicode[STIX]{x1D70C}}(\bar{u}-c){\hat{w}}^{\prime }-\bar{\unicode[STIX]{x1D70C}}\bar{u}^{\prime }{\hat{w}}-\bar{\unicode[STIX]{x1D70C}}\frac{g{\hat{w}}}{\bar{u}-c}\right.\unicode[STIX]{x27E7}=0. & \displaystyle\end{eqnarray}$$

Figure 6.

Base state velocity and density plots for the broken-line profiles. (a) Boussinesq/non-Boussinesq free surface case. (b) Rayleigh/Kelvin–Helmholtz free surface model. (c) Reduced-order SG-LV model. (d) Reduced-order SG-IG model.

Figure 7.

Growth rate contours in the $\unicode[STIX]{x1D6FC}{-}J_{2}$ plane for the broken-line profiles. Each row represents growth rates for a given $H/h$ . (a–c) $H/h=10$ , (d–f) $H/h=4$ , (g–i) $H/h=3$ and (j–l) $H/h=2$ . (a,d,g,j), (b,e,h,k) and (c,f,i,l) respectively correspond to Boussinesq free surface, non-Boussinesq free surface and Rayleigh/Kelvin–Helmholtz free surface.

To capture the essence of the base state profiles given in (3.1) and (3.2), the corresponding broken-line base state velocity and density profiles are given in the following equations, and are schematically shown in figure 6(a):

(4.5a ) $$\begin{eqnarray}\displaystyle \bar{u}(z)=\left\{\begin{array}{@{}ll@{}}U & \phantom{-}z\geqslant h,\\ U{\displaystyle \frac{z}{h}} & \phantom{-}h\geqslant z\geqslant -h,\\ -U & -h\geqslant z\geqslant -H,\end{array}\right. & & \displaystyle\end{eqnarray}$$

(4.5b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}(z)=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{1} & z>H,\\ \unicode[STIX]{x1D70C}_{2} & H>z>0,\\ \unicode[STIX]{x1D70C}_{3} & 0>z>-H.\end{array}\right. & & \displaystyle\end{eqnarray}$$

$T_{3}$

${\mathcal{D}}_{1}(c,\unicode[STIX]{x1D6FC};J_{1};J_{2};H/h)=0$

${\mathcal{D}}_{1}$

(4.6) $$\begin{eqnarray}\displaystyle {\mathcal{D}}_{1}=\left|\begin{array}{@{}cccccc@{}}\text{e}^{\unicode[STIX]{x1D6FC}H/h}(-P_{1}+P_{2}) & \text{e}^{-\unicode[STIX]{x1D6FC}H/h}(P_{1}+P_{2}) & 0 & 0 & 0 & 0\\ \text{e}^{\unicode[STIX]{x1D6FC}} & \text{e}^{-\unicode[STIX]{x1D6FC}}(-2P_{1}+1) & 0 & 2\text{e}^{-\unicode[STIX]{x1D6FC}}P_{1} & 0 & 0\\ 0 & 0 & 0 & 2\unicode[STIX]{x1D6FC}c & -{\displaystyle \frac{J_{2}}{c}} & -2\unicode[STIX]{x1D6FC}c-{\displaystyle \frac{J_{2}}{c}}\\ 0 & 0 & 0 & 0 & P_{3} & P_{4}\\ \text{e}^{2\unicode[STIX]{x1D6FC}} & 1 & -\text{e}^{2\unicode[STIX]{x1D6FC}} & -1 & 0 & 0\\ 0 & 0 & 1 & 1 & -1 & -1\end{array}\right|,\quad & & \displaystyle\end{eqnarray}$$

where

(4.7) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}P_{1}=\unicode[STIX]{x1D6FC}(1-c),\quad P_{2}={\displaystyle \frac{J_{1}}{1-c}},\quad P_{3}=-\text{e}^{-\unicode[STIX]{x1D6FC}}\left[{\displaystyle \frac{2\unicode[STIX]{x1D6FC}(1+c)\text{e}^{\unicode[STIX]{x1D6FC}}}{\text{e}^{\unicode[STIX]{x1D6FC}}-\text{e}^{\unicode[STIX]{x1D6FC}(2H/h\,-1)}}}+1\right],\\[12.0pt] P_{4}=-\text{e}^{\unicode[STIX]{x1D6FC}}\left[{\displaystyle \frac{2\unicode[STIX]{x1D6FC}(1+c)\text{e}^{\unicode[STIX]{x1D6FC}(2H/h\,-1)}}{\text{e}^{\unicode[STIX]{x1D6FC}}-\text{e}^{\unicode[STIX]{x1D6FC}(2H/h\,-1)}}}+1\right].\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Since $\unicode[STIX]{x1D70C}_{1}\ll \unicode[STIX]{x1D70C}_{2}$ and $\unicode[STIX]{x1D70C}_{2}\approx \unicode[STIX]{x1D70C}_{3}$ we have

(4.8a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle J_{1}=\frac{\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x1D70C}_{2}+\unicode[STIX]{x1D70C}_{1}}\frac{gh}{U^{2}}\approx \frac{gh}{U^{2}},\quad J_{2}=\frac{2(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})}{\unicode[STIX]{x1D70C}_{3}+\unicode[STIX]{x1D70C}_{2}}\frac{gh}{U^{2}}\approx \frac{\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2}}{\unicode[STIX]{x1D70C}_{2}}\frac{gh}{U^{2}}. & \displaystyle\end{eqnarray}$$

The dispersion relation ${\mathcal{D}}_{1}=0$ is sixth-order due to the presence of six waves in the system, namely two vorticity waves (each present at a vorticity jump), two surface gravity waves (present at the free surface) and two interfacial gravity waves (present at the pycnocline). The dispersion relation is solved using the MATLAB routine ‘roots’. Growth rate contours for the non-Boussinesq free surface case are shown in the second column of figure 7. If $T_{3}=0$ (which is our inference), the first and second columns of figure 7 should be exactly the same. This can be confirmed by observing that the first and second columns of figure 7 are indeed identical.

One important objective behind performing broken-line analysis is to find qualitative as well as quantitative agreement with the corresponding smooth profile. Comparison between the third column of figure 3 and the second column of figure 7 clearly shows that the broken-line profile thoroughly captures the essence of different types of instabilities existing in the system. The quantitative agreement between smooth and broken-line profiles can be observed by comparing table 1 with table 2. We also note that the thin unstable protuberances (which extend to infinity) visible in figure 7(d–l) are not observed in figure 3, since the growth rate values in these protuberances are very small (nearly zero); even a tiny amount of viscosity is enough to dampen these instabilities. In fact, even in the purely inviscid case depicted in figure 7(d–l), these narrow unstable regions are obtained only after very careful computations.

Table 2.

Growth rates and phase speeds for the maximum growth rate mode of the broken-line profiles

The presence of six waves in the system makes it difficult to pinpoint the mechanisms leading to each type of instability. This motivated us to construct simplified or reduced-order models having a lesser number of waves. Since SG-LV is a new kind of instability and the mechanism leading to its formation is not yet fully understood, a reduced-order set-up is constructed that has a Rayleigh/Kelvin–Helmholtz type velocity profile and a free surface. This implies that the pycnocline, and therefore the two interfacial gravity waves, have been removed from the system. This profile has been studied by both Longuet-Higgins (Reference Longuet-Higgins1998) and Bakas & Ioannou (Reference Bakas and Ioannou2009), and is given in figure 6(b). We have termed this model as the ‘Rayleigh/Kelvin–Helmholtz free surface’ model. The base velocity and density profiles are as follows:

(4.9a ) $$\begin{eqnarray}\displaystyle \bar{u}(z)=\left\{\begin{array}{@{}ll@{}}U & \phantom{-}z\geqslant h,\\ U{\displaystyle \frac{z}{h}} & \phantom{-}h\geqslant z\geqslant -h,\\ -U & -h\geqslant z\geqslant -H,\end{array}\right. & & \displaystyle\end{eqnarray}$$

(4.9b ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D70C}}(z)=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{1} & z>H,\\ \unicode[STIX]{x1D70C}_{2} & H>z>-H.\end{array}\right. & & \displaystyle\end{eqnarray}$$

Following a procedure similar to the non-Boussinesq free surface case, a dispersion relation ${\mathcal{D}}_{2}(c,\unicode[STIX]{x1D6FC};J_{1};H/h)=0$ can be obtained for the Rayleigh/Kelvin–Helmholtz free surface model. This dispersion relation is fourth-order (due to the presence of four waves in the system), and is obtained by evaluating the following determinant and equating it to zero:

(4.10) $$\begin{eqnarray}\displaystyle {\mathcal{D}}_{2}=\left|\begin{array}{@{}ccccc@{}}\text{e}^{\unicode[STIX]{x1D6FC}H/h}(-P_{1}+P_{2}) & \text{e}^{-\unicode[STIX]{x1D6FC}H/h}(P_{1}+P_{2}) & 0 & 0 & 0\\ \text{e}^{\unicode[STIX]{x1D6FC}} & \text{e}^{-\unicode[STIX]{x1D6FC}}(-2P_{1}+1) & 0 & 2\text{e}^{-\unicode[STIX]{x1D6FC}}P_{1} & 0\\ 0 & 0 & 0 & -P_{5} & P_{5}+P_{6}\\ \text{e}^{2\unicode[STIX]{x1D6FC}} & 1 & -\text{e}^{2\unicode[STIX]{x1D6FC}} & -1 & 0\\ 0 & 0 & -\text{e}^{-\unicode[STIX]{x1D6FC}} & -\text{e}^{\unicode[STIX]{x1D6FC}} & -P_{6}\end{array}\right|. & & \displaystyle\end{eqnarray}$$

Here $P_{5}=-2\text{e}^{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D6FC}(1+c)$ , $P_{6}=\text{e}^{\unicode[STIX]{x1D6FC}(2H/h\,-1)}-\text{e}^{\unicode[STIX]{x1D6FC}}$ . Contours of constant growth rate for varying $H/h$ are plotted in the third column of figure 7. To have parity between the different growth rate contours of figure 7, we have plotted the Rayleigh/Kelvin–Helmholtz free surface model in the $\unicode[STIX]{x1D6FC}{-}J_{2}$ plane. The Rayleigh/Kelvin–Helmholtz free surface model has no $J_{2}$ inherently (since there is no pycnocline), but we can scale $J_{1}$ by $(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})/\unicode[STIX]{x1D70C}_{2}$ to yield $J_{2}$ . For large values of $J_{2}$ , the third column of figure 7 reveals that the Rayleigh/Kelvin–Helmholtz free surface model very closely resembles Rayleigh/KH instability. The SG-LV branch of instability arising from the non-Boussinesq free surface case as well as the Rayleigh/Kelvin–Helmholtz free surface model are qualitatively and quantitatively similar. This implies that the interfacial gravity waves do not play any role in the SG-LV mode. We obtained the same conclusion from the eigenfunction analysis of smooth profiles in § 3.4.2.

We emphasize here that SG-LV mode has also been observed by Longuet-Higgins (Reference Longuet-Higgins1998) and Bakas & Ioannou (Reference Bakas and Ioannou2009), and is referred to as the ‘branch II instability’. From the previous literature it is not clear exactly what leads to the formation of this mode, i.e. which waves play the central role in this instability. Based on the intuitive understanding of resonant interactions between counter-propagating waves, we argue that only two waves are essential for the SG-LV instability, the leftward moving surface gravity wave (wave $2$ in figure 1) and the rightward moving vorticity wave existing at the lower vorticity interface (wave $6$ in figure 1).

To prove the above-mentioned hypothesis as well as to gain a deeper understanding into the formation of the SG-LV mode, we have constructed a further reduced-order model that consists only of two oppositely propagating surface gravity waves and a rightward propagating vorticity wave. We call this the SG-LV model; see figure 6(c). The base state profile for this model is given by

(4.11a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \bar{u}(z)=\left\{\begin{array}{@{}ll@{}}U{\displaystyle \frac{2\,h}{H+h}}\left({\displaystyle \frac{z}{h}}\right)-U{\displaystyle \frac{H-h}{H+h}} & \phantom{-}z\geqslant -h,\\ U{\displaystyle \frac{z}{h}} & -h\geqslant z\geqslant -H,\end{array}\right.\end{eqnarray}$$

(4.11b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \bar{\unicode[STIX]{x1D70C}}(z)=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{1} & z>H,\\ \unicode[STIX]{x1D70C}_{2} & H>z>-H.\end{array}\right.\end{eqnarray}$$

The dispersion relation ${\mathcal{D}}_{3}(c,\unicode[STIX]{x1D6FC};J_{1};H/h)=0$ is obtained by equating the following determinant to zero:

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{D}}_{3}=\left|\begin{array}{@{}cc@{}}\text{e}^{\unicode[STIX]{x1D6FC}H/h}(-P_{1}+P_{2}+P_{7}) & \text{e}^{-\unicode[STIX]{x1D6FC}H/h}(P_{1}+P_{2}+P_{7})\\ P_{3}+\text{e}^{-\unicode[STIX]{x1D6FC}}(1-P_{7}) & P_{4}+\text{e}^{\unicode[STIX]{x1D6FC}}(1-P_{7})\end{array}\right|, & \displaystyle\end{eqnarray}$$

where $P_{7}=2h/H+1$ . It is a third-order equation due to the presence of three waves in the system. To have the same Doppler shift as in the Rayleigh/Kelvin–Helmholtz free surface model, we kept the velocity at the free surface as $U$ . A consequence of keeping the same Doppler shift is that the shear in the system changes, hence the (non-dimensional) growth rates also change since they depend on the shear scale. The shear scale for the SG-LV model is chosen such that the maximum growth rate of this model is equal to that of the Rayleigh/Kelvin–Helmholtz free surface model. Like the Rayleigh/Kelvin–Helmholtz free surface model, $J_{1}$ has been scaled by $(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})/\unicode[STIX]{x1D70C}_{2}$ in the reduced-order SG-LV model in order to get an equivalent $J_{2}$ . This would maintain parity between different growth rate contours. A comparison between the non-Boussinesq free surface, Rayleigh/Kelvin–Helmholtz free surface and SG-LV models for $H/h=2$ can be made from figure 8. We have chosen $H/h=2$ specifically because the SG-LV mode is more prominent for small values of $H/h$ . Growth rates and stability boundaries of the SG-LV model qualitatively as well as quantitatively agree with both the Rayleigh/Kelvin–Helmholtz free surface and the non-Boussinesq free surface cases. This clearly establishes that SG-LV results because of the interaction between waves $2$ and $6$ .

A small but important point worth mentioning is that there is a non-zero shear at the free surface in the SG-LV model. This shear term would modify the characteristics of the surface gravity waves (Ehrnström & Villari Reference Ehrnström and Villari2008). This is the main reason behind the minor qualitative differences between the Rayleigh/Kelvin–Helmholtz free surface model (figure 8 b) and the SG-LV model (figure 8 c).

Figure 8.

Growth rate contours for $H/h=2$ case for (a) the non-Boussinesq free surface, (b) the Rayleigh/Kelvin–Helmholtz profile with free surface and (c) the reduced-order SG-LV model containing a free surface and a lower vorticity jump.

Apart from the SG-LV mode, there is yet another instability mode that is only apparent in the broken-line profile. This is the second narrow branch observed near $\unicode[STIX]{x1D6FC}=1$ and $J_{2}=0.01$ in the figure 7(k). We have termed it as the ‘SG-IG mode’. To understand and isolate the essential features of this mode, we have constructed the reduced-order SG-IG model, the base state profiles of which are given by (and also shown in figure 6 d)

(4.13a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \bar{u}(z)=U\frac{h}{H}\left(\frac{z}{h}\right),\end{eqnarray}$$

(4.13b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \bar{\unicode[STIX]{x1D70C}}(z)=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}_{1} & z>H,\\ \unicode[STIX]{x1D70C}_{2} & H>z>0,\\ \unicode[STIX]{x1D70C}_{3} & 0>z>-H.\end{array}\right.\end{eqnarray}$$

The dispersion relation for the SG-IG model is given by ${\mathcal{D}}_{4}(c,\unicode[STIX]{x1D6FC};J_{1};J_{2};H/h)=0$ , where

(4.14) $$\begin{eqnarray}\displaystyle {\mathcal{D}}_{4}=\left|\begin{array}{@{}cc@{}}\text{e}^{\unicode[STIX]{x1D6FC}H/h}\left(-P_{1}+P_{2}+{\displaystyle \frac{h}{H}}\right) & \text{e}^{-\unicode[STIX]{x1D6FC}H/h}\left(P_{1}+P_{2}+{\displaystyle \frac{h}{H}}\right)\\ -P_{8}-{\displaystyle \frac{J_{2}}{c}} & -P_{8}\text{e}^{2\unicode[STIX]{x1D6FC}H/h}-{\displaystyle \frac{J_{2}}{c}}\end{array}\right|, & & \displaystyle\end{eqnarray}$$

and $P_{8}=2\unicode[STIX]{x1D6FC}c/(1-\text{e}^{2\unicode[STIX]{x1D6FC}H/h})$ . Four waves are present in the SG-IG model – two surface gravity waves and two interfacial gravity waves – hence ${\mathcal{D}}_{4}=0$ is a fourth-order relation. Base state shear in the non-Boussinesq free surface case is different from the SG-IG model. Like SG-LV, here also we have chosen an appropriate shear scale for the SG-IG model so that the maximum growth rate matches that obtained for the non-Boussinesq free surface case. Figure 9 shows the growth rate contours for the case of non-Boussinesq free surface and the SG-IG model for $H/h=2$ . The SG-IG branch of the non-Boussinesq free surface case and the reduced-order SG-IG model are in qualitative and quantitative agreement. Some small differences arise because shear is present at the free surface in the reduced-order SG-IG model (just like the SG-LV model). Our analysis strongly points to the fact that the lower branch of instability appearing in the non-Boussinesq free surface case (figure 7 k) is due to the interaction between counter-propagating surface gravity and interfacial gravity waves (i.e. waves $2$ and $4$ of figure 1).

Figure 9.

Contours of constant growth rate in the $\unicode[STIX]{x1D6FC}{-}J_{2}$ plane for $H/h=2$ , depicting the SG-IG mode for (a) the non-Boussinesq free surface and (b) the reduced-order SG-IG model that can only support surface gravity waves and interfacial gravity waves.

Figure 10.

Dispersion diagram for the reduced-order SG-LV model for $J_{2}=0$ (which corresponds to the maximum growth rate case). (a) Growth rate of the unstable mode. (b) Real frequency, indicated by dark lines, obtained by solving ${\mathcal{D}}_{3}=0$ . The instability region is marked by two small parallel lines. The lighter lines represent the real frequencies associated with the isolated stable waves. The real frequencies of isolated waves cross near the point of maximum growth rate.

4.2 Dispersion diagrams

A standard approach to understanding shear instabilities is via dispersion diagrams (Craik Reference Craik1988), which is followed here to focus on the two new types of instability, viz. SG-LV and SG-IG. The dispersion diagrams of SG-LV and SG-IG models are respectively plotted in figures 10 and 11. The dark lines correspond to the solutions of each dispersion relation ( ${\mathcal{D}}_{3}=0$ for SG-LV and ${\mathcal{D}}_{4}=0$ for SG-IG); the imaginary part of the frequency ( $\unicode[STIX]{x1D714}_{i}$ , which signifies the growth rate) is shown in figures 10(a) and 11(a), while the real part ( $\unicode[STIX]{x1D714}_{r}$ ) is shown in figures 10(b) and 11(b). For a given $\unicode[STIX]{x1D6FC}$ , there should be $n$ roots (signifying $n$ waves) as demanded by the dispersion relation. For the SG-LV model $n=3$ while for SG-IG $n=4$ . Figure 10(b) reveals that there are three distinct $\unicode[STIX]{x1D714}_{r}$ values corresponding to each $\unicode[STIX]{x1D6FC}$ , except for the range $0.16<\unicode[STIX]{x1D6FC}<0.56$ . In this range, two roots coalesce into one. In fact, as $\unicode[STIX]{x1D6FC}$ increases from $0$ , two constituent waves come close together and coalesce at $\unicode[STIX]{x1D6FC}=0.16$ , and then bifurcate into two waves at $\unicode[STIX]{x1D6FC}=0.56$ . A pair of complex conjugate roots (growing and decaying normal modes) then become a possible solution, and its presence is verified from figure 10(a). The dispersion diagram obtained by solving the Rayleigh/Kelvin–Helmholtz free surface model would produce a similar dispersion diagram (the number of curves would be different because the SG-LV has fewer roots/waves), as can be found in Longuet-Higgins (Reference Longuet-Higgins1998) and Bakas & Ioannou (Reference Bakas and Ioannou2009).

Figure 11.

The same as figure 10 but for the reduced-order SG-IG model (here $J_{2}=0.0026$ ).

The entire procedure applied to the SG-LV model can also be applied to understand the dispersion diagrams of the SG-IG model given in figure 11(a,b). Furthermore, the understanding obtained from the dispersion diagrams (both SG-LV and SG-IG) can be augmented by adding the dispersion relations of individual waves in isolation. The location where two isolated waves cross each other in the $\unicode[STIX]{x1D6FC}{-}\unicode[STIX]{x1D714}_{r}$ plane signifies the resonant condition (i.e. waves have the same $\unicode[STIX]{x1D714}_{r}$ , which is the Doppler-shifted frequency). The dispersion relations of the isolated waves are plotted as lighter lines in figures 10(b) and 11(b). Focusing on the SG-LV, the dispersion relation of the isolated lower vorticity wave after appropriate Doppler shift is given by (Sutherland Reference Sutherland2010)

(4.15)

Similarly, the dispersion relation for the isolated surface gravity waves affected by linear shear after appropriate Doppler shift is given by (Ehrnström & Villari Reference Ehrnström and Villari2008)

(4.16)

The two dispersion curves, (rightward vorticity wave) and (leftward surface gravity wave) cross near $\unicode[STIX]{x1D6FC}=0.41$ in figure 10(b). This wavenumber corresponds to the most unstable mode, as can be clearly seen from figure 10(a). For large values of $\unicode[STIX]{x1D6FC}$ , the isolated dispersion curves asymptote to the dispersion curves obtained from solving ${\mathcal{D}}_{3}=0$ . We note here that solutions of ${\mathcal{D}}_{3}=0$ give $\unicode[STIX]{x1D714}_{r}$ of interacting waves. Since the eigenfunction of a wave decays exponentially, interaction is nearly zero when $\unicode[STIX]{x1D6FC}$ is large. Hence for large $\unicode[STIX]{x1D6FC}$ , each curve in the dispersion diagram basically represents an isolated wave. The fast neutral mode obtained from the dispersion relation matches very closely with the isolated rightward surface gravity wave, , which does not interact with any other wave in the system. This is because the intrinsic phase speed of the rightward surface gravity wave is not opposed to the local mean flow. Hence it cannot form a counter-propagating pair with any of the other waves present in the system and thus does not take part in any of the unstable modes.

Similar to the SG-LV model, the instability arising in the SG-IG model for $0.64<\unicode[STIX]{x1D6FC}<0.84$ can also be understood in terms of interacting waves in isolation. To conclusively demonstrate that the SG-IG mode arises due to the interaction between a surface gravity wave and an interfacial gravity wave, we consider them in isolation and plot their dispersion relations with lighter lines in figure 11(b). The dispersion relation for the stable isolated interfacial gravity waves is given by (Sutherland Reference Sutherland2010)

(4.17)

For stable isolated surface gravity waves in the presence of shear, the dispersion relation is given by

(4.18)

The isolated waves (leftward surface gravity wave) and (rightward interfacial gravity wave) intersect at $\unicode[STIX]{x1D6FC}\approx 0.77$ in figure 11(b). This wavenumber corresponds to the most unstable mode in figure 11(a). The fast neutral mode obtained from the dispersion relation matches very closely with the isolated rightward surface gravity wave, , which does not interact with any other wave in the system.

From figure 11(b) it can be seen that the dispersion relations of the isolated leftward surface gravity wave and isolated leftward interfacial gravity wave also intersect. This intersection does not lead to an exponential instability because the intrinsic phase speeds of the leftward surface gravity wave and the leftward interfacial gravity wave are not opposite to each other. Wave interaction theory requires that the intrinsic phase speeds of exponentially growing waves have to be in the opposing direction (Heifetz & Methven Reference Heifetz and Methven2005; Carpenter et al. Reference Carpenter, Tedford, Heifetz and Lawrence2013; Guha & Lawrence Reference Guha and Lawrence2014).

Finally, we emphasize here that the identification of isolated waves which resonate to produce the observed instabilities is non-trivial. For Rayleigh, Kelvin–Helmholtz, Holmboe or Taylor–Caulfield instabilities, identification of the isolated waves is comparatively easier, since in those systems deep water internal gravity waves and/or deep water vorticity waves are present. Dispersion relations of these waves are well established. In our case, the component waves of both SG-LV and SG-IG are not that simple, and their dispersion relations are not very well known. For example, in the formation of the SG-IG mode, the ‘SG wave’ is not a deep water surface gravity wave. It turns out to be an intermediate surface gravity wave modified by linear shear, whose dispersion relation is given by (4.18). The same is true for all the isolated waves yielding SG-LV and SG-IG modes; see (4.15)–(4.18). Without the construction of minimal models like SG-LV (see figure 6 c) and SG-IG (see figure 6 d), identification of the key waves leading to these instabilities would not have been possible.