Effect of free surface on submerged stratified shear instabilities

TL;DR

This study investigates how shallow water domains and free surface effects influence stratified shear instabilities, revealing new instability modes and the limitations of the rigid-lid approximation in modeling such flows.

Contribution

It introduces new instability modes caused by free surface interactions and compares shallow and deep domain behaviors, advancing understanding of stratified shear flow dynamics.

Findings

Discovery of two new instabilities: SG-IG and SG-LV modes.

Surface effects significantly alter instability characteristics in shallow domains.

Rigid-lid approximation is valid only for large vertical domains.

Abstract

Here we have considered the effects of shallowness of the domain as well as the air-water free surface on the stratified shear instabilities of the fluid underneath. First, we numerically solve the non-Boussinesq Taylor-Goldstein equation for smooth velocity and density profiles of a model shear layer with a free surface. When the depth of the fluid is relatively shallow compared to the shear layer thickness, the surface gravity waves existing at the free surface come closer to the waves existing in the shear layer. This can lead to resonant wave interactions, making the flow unstable to more varieties of modal instabilities. To obtain a deeper understanding of the instability mechanisms, we have performed analytical studies with broken-line profiles. Furthermore, reduced order broken-line profiles have also been developed, based on which dispersion diagrams are constructed. Through these diagrams we have underpinned the resonantly interacting waves leading to each type of instability. Two new instabilities have been found; one of them, referred to as the surface gravity - interfacial gravity (SG-IG) mode, arises due to the interaction between a surface gravity wave and an interfacial gravity wave, and would be absent if there is no internal density stratification. The other one - the surface gravity - lower vorticity (SG-LV) mode, which arises due to the interaction between a surface gravity wave and the lower vorticity wave, surpasses Kelvin-Helmholtz instability to become the most unstable mode. Except for Holmboe instability, remarkable differences are observed in all other instabilities occurring in shallow domains when the air-water interface is replaced by a rigid lid. We infer that the rigid-lid approximation is valid for large vertical domains. We have also shown that if shear is absent at the free surface, our problem can be modeled using a Boussinesq type approximation.

Tables

Table 1. 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 refereed 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.

Case	$H/h$	Mode	$\alpha$	$J$	$\gamma$	$c_r$
Rigid lid	$10$	KH	$0.4474$	$0$	$0.1893$	$0$
	$10$	H	$0.6773$	$0.3498$	$0.0781$	$\pm 0.4140$
	$4$	KH	$0.4724$	$0$	$0.1769$	$0$
	$4$	H	$0.6815$	$0.3498$	$0.0768$	$\pm 0.4155$
	$3$	KH	$0.4892$	$0$	$0.1566$	$0$
	$3$	H	$0.6898$	$0.3498$	$0.0737$	$\pm 0.4191$
	$2$	KH	$0.4849$	$0$	$0.0966$	$0$
	$2$	H	$0.8152$	$0.5106$	$0.0618$	$\pm 0.4864$
Boussinesq free surface	$10$	KH	$0.4431$	$0$	$0.1893$	$-0.0001$
	$10$	H	$0.6775$	$\mathrm{0.0.3498}$	$0.078$	$-0.4146$
	$4$	SG-LV/KH	$0.4641$	$0$	$0.1816$	$-0.0153$
	$4$	H	$0.6856$	$0.3498$	$0.0767$	$-0.4158$
	$3$	SG-LV	$0.4850$	$0$	$0.1677$	$-0.0270$
	$3$	KH	$0.4139$	$0.03004$	$0.1035$	$-0.0034$
	$3$	H	$0.6982$	$0.3498$	$0.0734$	$-0.4198$
	$2$	SG-LV	$0.4933$	$0$	$0.1224$	$-0.0375$
	$2$	H	$0.8738$	$0.5106$	$0.05931$	$-0.4897$
Non-Boussinesq free surface	$10$	KH	$0.4431$	$0$	$0.1894$	$-0.0051$
	$10$	H	$0.6775$	$0.3498$	$0.0788$	$-0.4155$
	$4$	SG-LV	$0.4557$	$0$	$0.1884$	$-0.0513$
	$4$	KH	$0.46826$	$0.0152$	$0.15445$	$-0.00695$
	$4$	H	$0.6898$	$0.3498$	$0.07708$	$-0.41638$
	$3$	SG-LV	$0.4808$	$0$	$0.1844$	$-0.1010$
	$3$	KH	$0.4599$	$0.019$	$0.1209$	$-0.0073$
	$3$	H	$0.715$	$0.3498$	$0.0729$	$-0.4197$
	$2$	SG-LV	$0.5477$	$0$	$0.1658$	$-0.1789$
	$2$	KH	$0.3219$	$0.0227$	$0.0374$	$-0.0150$
	$2$	H	$0.8695$	$0.5106$	$0.0569$	$-0.4827$

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

Case	$H/h$	Mode	$\alpha$	$J_2$	$\gamma$	$c_r$
Boussinesq/Non-Boussinesq free surface	$10$	KH	$0.3973$	$0$	$0.20$	$-0.001$
	$10$	H	$0.7975$	$0.3973$	$0.1422$	$\pm 0.4667$
	$4$	SG-LV	$0.4163$	$0$	$0.1954$	$-0.0942$
	$4$	SG-IG	$0.7004$	$0.0045$	$0.0435$	$0.1106$
	$4$	KH	$0.3713$	$0.0185$	$0.1480$	$0$
	$4$	H	$0.8235$	$0.4388$	$0.1415$	$\pm 0.4844$
	$3$	SG-LV	$0.4313$	$0$	$0.1880$	$-0.1945$
	$3$	SG-IG	$0.7064$	$0.004$	$0.0507$	$0.0730$
	$3$	KH	$0.2662$	$0.023$	$0.0814$	$0.0128$
	$3$	H	$0.8565$	$0.4947$	$0.1389$	$\pm 0.5091$
	$2$	SG-LV	$0.4573$	$0$	$0.1597$	$-0.3980$
	$2$	SG-IG	$0.8425$	$0.005$	$0.0466$	$0.0532$
	$2$	H	$0.9655$	$0.7616$	$0.1251$	$\pm 0.6003$
Rayleigh/Kelvin-Helmholtz with free surface	$10$	R	$0.4003$	all	$0.2011$	$-0.00088$
	$4$	R	$0.4043$	$0.3228$	$0.1793$	$0$
	$4$	SG-LV	$0.4163$	$0$	$0.1961$	$-0.0937$
	$3$	R	$0.3813$	$0.9790$	$0.1436$	$0$
	$3$	SG-LV	$0.4303$	$0$	$0.1886$	$-0.1946$
	$2$	SG-LV	$0.4553$	$0$	$0.1600$	$-0.40$
SG-LV model	$2$	SG-LV	$0.4161$	$0$	$0.1542$	$-0.5751$
SG-IG model	$2$	SG-IG	$0.7716$	$0.0026$	$0.0505$	$0.0394$

Equations

$$\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}+w\frac{\partial q}{\partial z}=\frac{1}{\rho^{2}}\Big{(}\frac{\partial\rho}{\partial z}\frac{\partial p}{\partial x}-\frac{\partial\rho}{\partial x}\frac{\partial p}{\partial z}\Big{)}.$$

$$\frac{\partial\tilde{q}}{\partial t}+\bar{u}\frac{\partial\tilde{q}}{\partial x}=\,\,\,\,\underbrace{-\tilde{w}\frac{d\bar{q}}{dz}}_{\text{\clap{Barotropic}}}\,\,\,\,\,\,\,\,\,\,\,\,+\,\,\,\,\,\,\,\,\,\,\,\,\underbrace{\frac{g}{\bar{\rho}}\frac{\partial\tilde{\rho}}{\partial x}}_{\text{\clap{Gravitational baroclinic}}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\underbrace{\frac{1}{\bar{\rho}^{2}}\Big{(}\frac{d\bar{\rho}}{dz}\frac{\partial\tilde{p}}{\partial x}\Big{)}}_{\text{\clap{$T_{3}\equiv$Non-Boussinesq baroclinic}}}.$$

$$\underbrace{\bar{\rho}^{\prime}\big{[}(\bar{u}-c)\hat{w}^{\prime}-\bar{u}^{\prime}\hat{w}\big{]}}_{\text{\clap{Non-Boussinesq}}}-\underbrace{\frac{g\bar{\rho}^{\prime}}{\bar{u}-c}\hat{w}}_{\text{\clap{Gravitational}}}\,\,\,+\,\,\,\bar{\rho}\big{[}(\bar{u}-c)(\hat{w}^{\prime\prime}-\alpha^{2}\hat{w})\big{]}-\underbrace{\bar{\rho}\bar{u}^{\prime\prime}\hat{w}}_{\text{\clap{Barotropic}}}=0.$$

$${\bar{\rho }}^{\prime }[(\bar{u}-c){\hat{w}}^{\prime }-{\bar{u}}^{\prime }\hat{w}]+\bar{\rho }[(\bar{u}-c)({\hat{w}}^{\prime \prime }-{\alpha }^2\hat{w})-{\bar{u}}^{\prime \prime }\hat{w}]=-\mathrm{i}\alpha \hat{\rho }g-\frac{\mathrm{i}}{\alpha }\mu [{\hat{w}}^{\prime \prime \prime \prime }-2{\alpha }^2{\hat{w}}^{\prime \prime }+{\alpha }^4\hat{w}],$$

$$\mathrm{i}\alpha (\bar{u}-c)\hat{\rho }+\hat{w}{\bar{\rho }}^{\prime }=\kappa ({\hat{\rho }}^{\prime \prime }-{\alpha }^2\hat{\rho }).$$

$$\left[\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right]\left[\begin{array}{c}{\vartheta }_1\\ {\vartheta }_2\end{array}\right]=c\left[\begin{array}{cc}{N}_{11}& N_{12}\\ N_{21}& N_{22}\end{array}\right]\left[\begin{array}{c}{\vartheta }_1\\ {\vartheta }_2\end{array}\right].$$

$$M_{11}={\bar{\rho }}^{\prime }[-\bar{u}D+{\bar{u}}^{\prime }]+\bar{\rho }[-\bar{u}{D}^2+{\alpha }^2\bar{u}+{\bar{u}}^{\prime \prime }]-\frac{\mathrm{i}}{\alpha }\mu [D^4-2{\alpha }^2{D}^2+{\alpha }^4],$$

$$M_{12}=-\mathrm{i}\alpha g,M_{21}={\bar{\rho }}^{\prime },M_{22}=\mathrm{i}\alpha \bar{u}-\kappa [D^2-{\alpha }^2],{\vartheta }_1=\hat{w},{\vartheta }_2=\hat{\rho },$$

$$N_{11}=-{\bar{\rho }}^{\prime }D+\bar{\rho }[-D^2+{\alpha }^2],N_{12}=0,N_{21}=0,N_{22}=\mathrm{i}\alpha .$$

$$\hat{w}=0.$$

$$\frac{d\hat{u}}{dz}=0,$$

$$\frac{d^2\hat{w}}{d{z}^2}=0.$$

$$\frac{d\hat{\rho }}{dz}=0.$$

$$\bar{u}(z)=U\tanh\bigg{(}\frac{z-\frac{1}{2}H}{h}\bigg{)}.$$

$$\bar{\rho}(z)=\left\{\begin{array}[]{cc}\rho_{01}\Bigg{(}1-A_{t,\,0}\tanh\bigg{[}5\bigg{(}\dfrac{z-H}{h}\bigg{)}\bigg{]}\Bigg{)}&\quad 1.2H\geq z\geq 0.75H,\\ \\ \rho_{12}\Bigg{(}1-A_{t}\tanh\bigg{[}5\bigg{(}\dfrac{z-\frac{1}{2}H}{h}\bigg{)}\bigg{]}\Bigg{)}&\quad 0.75H\geq z\geq 0.\end{array}\right.$$

$$A_t=\frac{{\rho }_{cw}-{\rho }_{ww}}{{\rho }_{cw}+{\rho }_{ww}}=\frac{10.3\text{kg}{\text{m}}^{-3}}{(993.5+983.2)\text{kg}{\text{m}}^{-3}}\approx \mathrm{0.005.}$$

$$J=2A_t\frac{gh}{U^2}.$$

$$\frac{\partial \tilde{p}}{\partial x}=-\bar{\rho }\left[\frac{\mathbb{D}\tilde{u}}{\mathbb{D}t}+\tilde{w}\bar{q}\right],$$

$${\tilde{u}}_{\pm }=\pm \omega {\mathrm{e}}^{\mathrm{i}(\alpha x-\omega t)},$$

$$\{\bar{\rho }[(\bar{u}-c){\hat{w}}^{\prime }-{\bar{u}}^{\prime }\hat{w}]{\}}^{\prime }-\frac{{\bar{\rho }}^{\prime }g}{\bar{u}-c}\hat{w}-\bar{\rho }{\alpha }^2(\bar{u}-c)\hat{w}=0.$$

$${\hat{w}}^{\prime \prime }-{\alpha }^2\hat{w}=0.$$

$$[[\hat{w}]]=0,$$

$$[[\bar{\rho }(\bar{u}-c){\hat{w}}^{\prime }-\bar{\rho }{\bar{u}}^{\prime }\hat{w}-\bar{\rho }\frac{g\hat{w}}{\bar{u}-c}]]=0.$$

$$\bar{u}(z)=\left\{\begin{array}[]{cc}U&\quad z\geq h,\\ \\ U\dfrac{z}{h}&\quad h\geq z\geq-h,\\ \\ -U&\quad-h\geq z\geq-H,\end{array}\right.$$

$$\bar{\rho}(z)=\left\{\begin{array}[]{cc}\rho_{1}&\quad z>H,\\ \\ \rho_{2}&H>z>0,\\ \\ \rho_{3}&0>z>-H.\end{array}\right.$$

$${\mathcal{D}}_1=\left|\begin{array}{cccccc}{\mathrm{e}}^{\alpha H/h}(-P_1+P_2)& {\mathrm{e}}^{-\alpha H/h}(P_1+P_2)& 0& 0& 0& 0\\ {\mathrm{e}}^{\alpha }& {\mathrm{e}}^{-\alpha }(-2P_1+1)& 0& 2{\mathrm{e}}^{-\alpha }{P}_1& 0& 0\\ 0& 0& 0& 2\alpha c& -\frac{J_2}{c}& -2\alpha c-\frac{J_2}{c}\\ 0& 0& 0& 0& P_3& P_4\\ {\mathrm{e}}^{2\alpha }& 1& -{\mathrm{e}}^{2\alpha }& -1& 0& 0\\ 0& 0& 1& 1& -1& -1\end{array}\right|,$$

$$P_{1}=\alpha(1-c),\,P_{2}=\frac{J_{1}}{1-c},\,P_{3}=-{\rm e}^{-\alpha}\Bigg{[}\frac{2\alpha(1+c){\rm e}^{\alpha}}{{\rm e}^{\alpha}-{\rm e}^{\alpha\big{(}2H/h\,-1\big{)}}}+1\Bigg{]},\\ P_{4}=-{\rm e}^{\alpha}\Bigg{[}\frac{2\alpha(1+c){\rm e}^{\alpha\big{(}2H/h\,-1\big{)}}}{{\rm e}^{\alpha}-{\rm e}^{\alpha\big{(}2H/h\,-1\big{)}}}+1\Bigg{]}.$$

$$J_1=\frac{{\rho }_2-{\rho }_1}{{\rho }_2+{\rho }_1}\frac{gh}{U^2}\approx \frac{gh}{U^2},J_2=\frac{2({\rho }_3-{\rho }_2)}{{\rho }_3+{\rho }_2}\frac{gh}{U^2}\approx \frac{{\rho }_3-{\rho }_2}{{\rho }_2}\frac{gh}{U^2}.$$

$$\bar{u}(z)=\left\{\begin{array}[]{cc}U&\quad z\geq h,\\ \\ U\dfrac{z}{h}&\quad h\geq z\geq-h,\\ \\ -U&\quad-h\geq z\geq-H,\end{array}\right.$$

$$\bar{\rho}(z)=\left\{\begin{array}[]{cc}\rho_{1}&\quad z>H,\\ \\ \rho_{2}&H>z>-H.\end{array}\right.$$

$${\mathcal{D}}_2=\left|\begin{array}{ccccc}{\mathrm{e}}^{\alpha H/h}(-P_1+P_2)& {\mathrm{e}}^{-\alpha H/h}(P_1+P_2)& 0& 0& 0\\ {\mathrm{e}}^{\alpha }& {\mathrm{e}}^{-\alpha }(-2P_1+1)& 0& 2{\mathrm{e}}^{-\alpha }{P}_1& 0\\ 0& 0& 0& -P_5& P_5+P_6\\ {\mathrm{e}}^{2\alpha }& 1& -{\mathrm{e}}^{2\alpha }& -1& 0\\ 0& 0& -{\mathrm{e}}^{-\alpha }& -{\mathrm{e}}^{\alpha }& -P_6\end{array}\right|.$$

$$\bar{u}(z)=\left\{\begin{array}[]{cc}U\dfrac{2\,h}{H+h}\left(\dfrac{z}{h}\right)-U\dfrac{H-h}{H+h}&\quad z\geq-h,\\ \\ U\dfrac{z}{h}&\quad-h\geq z\geq-H.\\ \end{array}\right.$$

Full text

††volume: xxx

1 Introduction

Flows in the natural environment like lakes, estuaries and oceans are shallow, density stratified and exhibit a free surface. In presence of a background velocity shear, which may arise due to wind forcing, and/or exchange flow, it is possible for a stable density stratified flow to become unstable. The resulting instabilities are known as stratified shear instabilities, and a few well known examples are Kelvin-Helmholtz (KH) instability, Holmboe instability and Taylor-Caulfield (TC) instability. These instabilities are often observed in the pycnocline region (sharp density changes due to salinity, temperature, or both) in natural flows. Instabilities cause the interfacial gravity waves present at the pycnocline to grow and break, which often leads to intense turbulent mixing, and thereby profoundly affect the aquatic environment.

Shear instabilities arising in the environment have been conventionally modeled in a fluid of infinite vertical extent. The effect of the free surface (interface between air and water) is neglected, hence the Boussinesq approximation can be effectively applied. Boussinesq approximation neglects the density variation effect in the inertial terms but considers it in the buoyancy term Turner (1979). This approximation gives accurate results when the density differences are small compared to the mean background density (for example, slight density differences between warm (fresh) and cold (salty) waters). Therefore, Boussinesq approximation may not hold if the effect of the free surface (i.e. large density difference between air and water) is taken into account.

There have been numerous studies on linear stability analyses and direct numerical simulations of Boussinesq stratified shear instabilities, e.g. Smyth et al. (1988), Smyth & Peltier (1989), Smyth & Peltier (1991), Lawrence et al. (1991), Sutherland & Peltier (1992), Alexakis (2005), Carpenter et al. (2007), Smyth et al. (2007), Carpenter et al. (2010), Guha et al. (2013) and Rahmani et al. (2014). These studies have considered domains that are much larger vertically than the shear layer thickness (so as to emulate an infinite vertical domain). Hence the free surface is too far away to play any significant role in the instability processes occurring in the pycnocline. For numerical implementation, a “rigid lid” boundary condition at the air-water interface ensures that it has no dynamics. An important step towards understanding flows in natural environment is to take into account the finite vertical extent of the domain, while still using Boussinesq approximation. This is achieved by considering only the water body (of finite vertical extent) and neglecting the air above. The implementation is similar to that of unbounded (infinite extent) flows mentioned above - rigid lid is used as the upper boundary. Hazel (1972) and Haigh & Lawrence (1999) have shown that the presence of a rigid lid close to the shear layer significantly affects the stability characteristics. There have been, however, only a few studies which have considered the effect of the free surface on submerged shear instabilities Longuet-Higgins (1998); Bakas & Ioannou (2009). Here the fluid below is homogeneous, and similar to the case previously mentioned, the domain extends up to the free surface (air above is neglected). While these studies have definitely made a significant advancement over the rigid lid approximation in capturing the non-trivial effects of the free surface on submerged shear instabilities, the free surface modeling is still an approximate one - the non-Boussinesq effects have not been considered. Since the free surface is indeed non-Boussinesq, involving huge density jump between air and water, it needs to be modeled carefully. The correct equations for a non-Boussinesq interface have been outlined in Barros & Choi (2011) and Heifetz & Mak (2015). Barros & Choi (2011) have analyzed the non-Boussinesq effects for Holmboe instability while Heifetz & Mak (2015) have analyzed it for Taylor-Caulfield instability.

A key aspect of the free surface is that it supports surface gravity waves. These waves can interact with the different vorticity and interfacial gravity waves that are supported in the stratified shear layer. A schematic provided in figure 1 shows that there are six waves in the system, two surface gravity waves (marked by $1$ and $2$), two interfacial gravity waves at the pycnocline (marked by $4$ and $5$), and two vorticity waves, one at each vorticity jump (marked by $3$ and $6$). Mean flow profile of the shear layer Doppler shifts four out of these six waves (waves $4$ and $5$ are not Doppler shifted since $\bar{u}=0$ there). Shear instabilities can be conceptually understood in terms of resonant interaction at a distance between counter-propagating waves Holmboe (1962); Sakai (1989); Baines & Mitsudera (1994); Caulfield (1994); Heifetz et al. (1999); Heifetz & Methven (2005); Carpenter et al. (2013); Guha & Lawrence (2014). In a counter-propagating system of two waves (each present at its own interface), the intrinsic phase speed of the waves should be opposite to each other. Furthermore, each wave’s intrinsic phase speed should be opposite to the local mean flow (unless the local mean velocity is zero). For example, classic KH (or Rayleigh) instability will result via an interaction between waves $3$ and $6$ of figure 1, while classic Holmboe will be due to waves $3$ and $4$, as well as $5$ and $6$. The surface gravity waves can resonantly interact with (at least) one of the waves existing in the submerged shear layer. Intuitively, we can expect two additional interactions simply by searching for counter-propagating configurations: (i) waves $2$ and $4$ - surface gravity wave interacting with interfacial gravity wave (SG-IG), and (ii) waves $2$ and $6$ - surface gravity wave interacting with lower vorticity wave (SG-LV). In this paper we will explore whether these instabilities are actually possible. Intuitively, we can also expect the interactions to be more prominent if the free surface is not very far from the pycnocline. We note here in passing that waves $1$ and $3$ give a false impression of counter-propagating configuration. Careful observation reveals that the condition for counter-propagation is violated since the phase speed of wave $1$ is not opposite to the local mean flow.

The objective of the present study is to consider the effects of shallowness and the free surface on the stratified shear instabilities of the fluid below. The pycnocline is assumed to be sharp, so as to capture both KH and Holmboe modes. In natural settings, the free surface in many occasions is not far from the pycnocline as compared to the shear layer thickness. Furthermore, natural flows being mostly shallow compared to the shear layer thickness, one can expect a non-trivial effect of the free surface on the instabilities occurring inside the fluid. For example, while studying stratified shear instabilities in Ishikari River estuary Yoshida et al. (1998) observed KH and Holmboe at the interface (pycnocline) when there is no wind. However in presence of wind the situation is very different - “When the wind blows, gravity waves often arise to destroy the interface. The details of the complicated mechanisms of this transformation remain undescribed, however, as do those responsible for transport of salt to the surface once the interface is gone”.

The paper is organized as follows. In §2 we outline the fundamental equations describing the flow physics, specifically the non-Boussinesq Taylor-Goldstein equation. We also briefly state the numerical procedure to solve this eigenvalue problem. In §3 we perform numerical stability analysis on smooth base state profiles similar to figure 1. These profiles provide a good approximation of realistic density stratified shear flows in presence of a free surface. The broken-line profile corresponding to the continuous profile, which is required for providing a mechanistic understanding of the instabilities identified in §3, is discussed in §4. Reduced order broken-line profiles are devised to single out the key waves that generate different instability mechanisms due to the presence of the free surface. Discussions and conclusions are stated in §5.

2 Governing equations and eigenvalue problem

2.1 Vorticity equation

We model the density stratified shear layer using the $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

[TABLE]

Here vorticity is given by $q\equiv\partial u/\partial z-\partial w/\partial x$, while $\rho$, $p$ and $g$ denote density, pressure and acceleration due to gravity respectively. The terms on the right hand side of (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)=d\bar{u}/dz$, $p=\bar{p}(z)$ and $\rho=\bar{\rho}(z)$. The base state follows hydrostatic pressure balance $d\bar{p}/dz=-\bar{\rho}g$. Perturbations are added to the base flow: $u=\bar{u}(z)+\tilde{u}$, $w=\tilde{w}$, $p=\bar{p}(z)+\tilde{p}$, $\rho=\bar{\rho}(z)+\tilde{\rho}$, 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 (1). This results in the perturbation vorticity evolution equation

[TABLE]

Equation (2) can also be found in Heifetz & Mak (2015) (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 modeling the free surface (the highly non-Boussinesq interface between air and water), $T_{3}$ has been neglected in Bakas & Ioannou (2009) (see their equation (6)) as well as Longuet-Higgins (1998) (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 modeling 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){\rm e}^{{\rm i}\alpha(x-ct)}$. Here $\alpha$ and $c$ are respectively the real wavenumber and the complex phase speed ($c=c_{r}+{\rm i}c_{i}$), and $f$ could represent $u$, $w$, $p$, $q$ or $\rho$. Temporal normal mode form is substituted in the linearized Navier-Stokes equations, yielding

[TABLE]

Total derivative with respect to $z$ is denoted by ′. The non-Boussinesq Taylor-Goldstein equation obtained is same as the one obtained by Barros & Choi (2011), Barros & Choi (2014), and Carpenter et al. (2017). Terms denoted by the braces indicate the corresponding terms in the perturbation vorticity equation (2). We intend to solve (3) for smooth profiles of base state velocity and density, and for this numerical eigenvalue solver is necessary. The numerical solution of (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 ($\mu$) and diffusivity ($\kappa$). Therefore we numerically solve the viscous diffusive form of the non-Boussinesq Taylor-Goldstein equation, which is given by

[TABLE]

A complete derivation of (4a)-(4b) is provided in Appendix A. The Boussinesq limit of the above equation set can be found in the appendix of Smyth et al. (2011). In the next subsection we briefly describe the numerical strategy to solve the eigenvalue problem (4a)-(4b).

2.3 Solution of the eigenvalue problem

Equation set (4a)-(4b) can be converted into a generalized eigenvalue problem of the form $M\boldsymbol{\vartheta}=cN\boldsymbol{\vartheta}$ :

[TABLE]

Elements of the matrices are given by

[TABLE]

[TABLE]

[TABLE]

Equation (5) is an eigenvalue problem for eigenvalue $c$ and eigenfunctions $\hat{w}$ and $\hat{\rho}$. 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

[TABLE]

Furthermore, the free-slip boundary condition,

[TABLE]

along with the fluid being incompressibility gives,

[TABLE]

We use insulating boundary condition for $\hat{\rho}$, which is given by

[TABLE]

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 (6)-(8) have been applied to that 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. (2011). The numerical routine developed has been validated against the Boussinesq shear layer problem of Smyth et al. (1988) and that of surface gravity waves. Our primary interest is to capture inviscid instabilities like KH and Holmboe, hence viscosity ($\mu$) and diffusivity ($\kappa$) introduced in (4a)-(4b) 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 Numerical stability analysis of the smooth profiles

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. Since we are interested in analyzing the effect of the free surface on the shear layer below, we need to extend the velocity profile in the air region. This extension has been made such that (i) the velocity profile is continuously differentiable at the free surface, and (ii) the velocity remains nearly constant with $z$ away from the free surface. The first point ensures that the derivatives of base state velocity remain finite, which is required for the numerical stability of (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 (1957). The dimensional base state velocity variation is given by:

[TABLE]

Here $U$ (in m/s) represents the surface current and $h$ is the half shear layer thickness. $H$ is the dimensional half-depth of the channel (i.e.  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 getting an opportunity to resonantly interact with them). The dimensional base state density profile is given by

[TABLE]

Here $\rho_{01}=(\rho_{ww}+\rho_{air})/2=491\,\,\textrm{kg}/\textrm{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}=(\rho_{ww}-\rho_{air})/(\rho_{ww}+\rho_{air})=9.97\times 10^{-1}\approx 1$. Similarly the mean density for cold water and warm water is given by $\rho_{12}=(\rho_{cw}+\rho_{ww})/2=988.35\,\,\textrm{kg}/\textrm{m}^{3}$ and the Atwood number corresponding to the pycnocline is given by

[TABLE]

The densities of cold water and warm water are respectively given by $\rho_{cw}$ and $\rho_{ww}$.

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

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

[TABLE]

For reporting stability characteristics we have non-dimensionalized the wavenumber, $\alpha$, the phase speed, $c$, and the growth rate $\Im\{\alpha c\}$. The wavenumber is non-dimensionalized by 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 modeling the air-water free surface?

As noted in §2.1, a couple of previous studies Longuet-Higgins (1998); Bakas & Ioannou (2009) that considered the effect of the free surface on the shear instability underneath have ignored the term $T_{3}\equiv(\bar{\rho})^{-2}(d\bar{\rho}/dz)(\partial\tilde{p}/\partial x)$ appearing in (2). For the purpose of comparison we have analyzed 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)). Here we try to understand the implications of assuming $T_{3}=0$. There are two possible interpretations of $T_{3}=0$:

(i) $d\bar{\rho}/dz=0$: This is a valid approximation in the Boussinesq limit (where $\bar{\rho}\approx\textrm{constant}$). In our setting, $T_{3}$ will only matter at the free surface, where obviously $d\bar{\rho}/dz\neq 0$, rather is approximately a delta function. Hence assuming $d\bar{\rho}/dz=0$ for the air-water density jump is not correct.

(ii) $\partial\tilde{p}/\partial x=0$: In the inviscid limit

[TABLE]

where $\bar{q}=d\bar{u}/dz$ and $\mathbb{D}\tilde{u}/\mathbb{D}t\equiv\partial\tilde{u}/\partial t+\bar{u}\partial\tilde{u}/\partial x$ is the linearized material derivative of $\tilde{u}$; see (37). 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, upper fluid with density $\rho_{1}$ and lower fluid with density $\rho_{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

[TABLE]

where $+(-)$ indicates below (above), and $\omega=\sqrt{g\alpha(\rho_{2}-\rho_{1})/(\rho_{1}+\rho_{2})}$ is the frequency; see Kundu et al. (2012, Chapter 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 (13) yields $\partial\tilde{p}/\partial 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. (2012, Chapter 7.2). However in a generalized scenario where background shear is present at the density interface (i.e. $\bar{q}\neq 0$), (13) directly reveals that $\partial\tilde{p}/\partial x=-\bar{\rho}\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.

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 presence of a rigid lid, a Boussinesq free surface, and an actual/non-Boussinesq free surface for $H/h=10$ are shown in figures 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. 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 that of Hazel (1972) and Haigh & Lawrence (1999).

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 “bubble” and becomes a 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 wouldn’t 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{\psi}$, can be obtained from the vertical velocity, $\hat{w}$ via the relation $\hat{\psi}=-{\rm i}\hat{w}/\alpha$ ($\hat{w}$ is obtained by solving the eigenvalue problem (5)). Since eigenfunctions are non-unique, we normalize $\hat{\psi}$ by its maximum value for reporting purposes. Apart from $\hat{\psi}$, we have also plotted the contours of perturbation streamfunction $\tilde{\psi}=\Re\{\hat{\psi}{\rm e}^{{\rm i}\alpha 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{\rho}$ is also of significant interest. For plotting purposes, we have normalized it by the maximum value of the corresponding $\hat{\psi}$. The contours of perturbation density, $\tilde{\rho}$ 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 analyze 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 $(\alpha,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{\psi}|$ 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 gets 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).

A few salient features of Holmboe instability in presence of a free surface are worth noticing. The value of $\hat{\rho}$ is a few orders of magnitude smaller at the free surface than that at the interface; see figures 4(c)-(f). The same figures also show that the density contours at the free surface and the interface are $\pi$ phase shifted. These two features are the hallmark of the classic “baroclinic/internal/varicose mode” in layered flows in absence of background velocity shear; see Sutherland (2010, Chapter 2). The analogy of Holmboe instability with baroclinic mode seems appropriate because, like baroclinic mode, 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 $(\alpha,J)=(0.4808,0)$ (refer to figure 3(i)). Eigenfunctions of perturbation streamfunction and density are plotted in figure 5. Figures 5(c)-(f) show that $|\hat{\rho}|$ 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 (2010, Chapter 2). Although from the eigenfunction plot the mechanism behind 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 evident from table 1, the highest growth rate of the SG-LV mode for $H/h=3$ is comparable to KH in a nearly unbounded flow (compare with the growth rate of KH mode for $H/h=10$). The phase speed of this mode is $\approx-0.1$, and is intermediate between KH and Holmboe modes for the same $H/h$.

4 Broken-line profiles and mechanistic picture of instability

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 it allows the waves in the system to become localized. Furthermore, using such profiles, (3) can be solved analytically. An equivalent form of (3) is given below:

[TABLE]

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 (3) and (14) are reduced to $\bar{\rho}(\bar{u}-c)[\hat{w}^{\prime\prime}-\alpha^{2}\hat{w}]=0$. The discrete eigenspectrum can be obtained by assuming $(\bar{u}-c)\neq 0$, yielding

[TABLE]

The continuity of vertical velocity across an interface is imposed:

[TABLE]

where $\llbracket....\rrbracket$ denotes the difference across the interface. Continuity of pressure across the interface is ensured by the dynamic condition, obtained by integrating (14) across the interface from $z_{i}-\Updelta z$ to $z_{i}+\Updelta z$ and letting $\Updelta z\rightarrow 0$,

[TABLE]

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

[TABLE]

Base state profiles (18a)-(18b) are first solved using the Boussinesq version of (3), which ignores density variation in the inertia terms. Growth rate contours corresponding to the Boussinesq free surface case are shown in the first column of figure 7. From the discussion in §3.2 it was inferred that if shear is absent at the free surface (which is the case here), then the non-Boussinesq and Boussinesq modes would yield the same result. In other words, to demonstrate that the term $T_{3}$ in (2) is zero, we solve the base state profiles given in (18a)-(18b) using the complete equation (3). In this regard we solve (15) between two consecutive interfaces and use (16) and (17) as boundary conditions. Two additional boundary conditions that are used are impervious bottom boundary and wave evanescence very far away from the free surface.

Thereby we obtain a $6$th order dispersion relation for the non-Boussinesq free surface case, which is given by $\mathcal{D}_{1}(c,\alpha;J_{1};J_{2};H/h)=0$. In determinant form $\mathcal{D}_{1}$ is written as

[TABLE]

where

[TABLE]

Since $\rho_{1}\ll\rho_{2}$ and $\rho_{2}\approx\rho_{3}$ we have

[TABLE]

The dispersion relation $\mathcal{D}_{1}=0$ is $6$th 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 with 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 extends to infinity) visible in figures 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 figures 7(d)-(l), these narrow unstable regions are obtained only after very careful computations.

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 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 (1998) and Bakas & Ioannou (2009), 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:

[TABLE]

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

[TABLE]

Here $P_{5}=-2{\rm e}^{\alpha}\alpha(1+c)$, P_{6}={\rm e}^{\alpha\big{(}2H/h\,-1\big{)}}-{\rm e}^{\alpha}. Contours of constant growth rate for varying $H/h$ are plotted in 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 $\alpha-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 $(\rho_{3}-\rho_{2})/\rho_{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 (1998) and Bakas & Ioannou (2009), 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 behind the formation of 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). Base state profile for this model is given below:

[TABLE]

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

[TABLE]

where $P_{7}=2h/H+1$. It is a $3$rd order equation due to 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 to be $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 $(\rho_{3}-\rho_{2})/\rho_{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 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 the 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, 2008). 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)).

Apart from 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 $\alpha=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 below (and also shown in figure 6(d)):

[TABLE]

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

[TABLE]

and $P_{8}=2\alpha c/(1-{\rm e}^{2\alpha 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).

4.2 Dispersion diagrams

A standard approach to understanding shear instabilities is via dispersion diagrams (Craik, 1988), 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 ($\omega_{i}$, which signifies the growth rate) is shown in figures 10(a) and 11(a), while the real part ($\omega_{r}$) is shown in figures 10(b) and 11(b). For a given $\alpha$, 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 $3$ distinct $\omega_{r}$ values corresponding to each $\alpha$, except for the range $0.16<\alpha<0.56$. In this range, two roots coalesce into one. In fact, as $\alpha$ increases from [math], two constituent waves come close together and coalesce at $\alpha=0.16$, and then bifurcate into two waves at $\alpha=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 similar dispersion diagram (number of curves would be different because SG-LV has lesser roots/waves), as can be found in Longuet-Higgins (1998) and Bakas & Ioannou (2009).

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 figures 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 $\alpha-\omega_{r}$ plane signifies the resonant condition (i.e. waves have the same $\omega_{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 SG-LV, the dispersion relation of the isolated lower vorticity wave after appropriate Doppler shift is given by (Sutherland, 2010)

[TABLE]

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

[TABLE]

The two dispersion curves, $\mathpzc{V}^{+}$ (rightward vorticity wave) and $\mathpzc{SG_{1}}^{-}$ (leftward surface gravity wave) cross near $\alpha=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 $\alpha$, 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 $\omega_{r}$ of interacting waves. Since eigenfunction of a wave decays exponentially, interaction is nearly zero when $\alpha$ is large. Hence for large $\alpha$, 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, $\mathpzc{SG_{1}}^{+}$, 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 can not 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<\alpha<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 below (Sutherland, 2010):

[TABLE]

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

[TABLE]

The isolated waves $\mathpzc{SG_{2}}^{-}$ (leftward surface gravity wave) and $\mathpzc{IG}^{+}$ (rightward interfacial gravity wave) intersect at $\alpha\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, $\mathpzc{SG_{2}}^{+}$, which does not interact with any other wave in the system.

From figure 11(b) it can be seen that the dispersion relations of 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 opposing direction (Heifetz & Methven, 2005; Carpenter et al., 2013; Guha & Lawrence, 2014).

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 are 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 (29). The same is true for all the isolated waves yielding SG-LV and SG-IG modes; see (26)-(29). 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.

5 Discussions and Conclusions

In this paper we have considered the effects of the vertical extent of the domain and the air-water free surface on the stratified shear instabilities of the fluid below. In the existing literature it is often customary to replace the air-water interface by a rigid lid, thereby simplifying the problem under the Boussinesq approximation. However the air-water interface is a free surface - a dynamically evolving boundary that supports surface gravity waves. Since shear instabilities arise from resonant wave interactions, it is possible for the surface gravity waves to resonantly interact with the other waves present in the stratified shear layer (e.g. vorticity waves and interfacial gravity waves). This may modify the existing instabilities, and furthermore, can lead to newer instabilities. The coupling can be made more feasible by bringing the free surface closer to the shear layer. By considering the vertical domain to extend only up to the free surface, and furthermore, by approximating the free surface as a “Boussinesq” interface, Longuet-Higgins (1998) and Bakas & Ioannou (2009) used broken-line profiles to understand how surface waves interact with the Rayleigh instability in the shear layer (they considered the fluid underneath to be homogeneous, i.e. an unstratified water body). In §3.2 we have shown that the air-water free surface can be treated as a Boussinesq interface (by neglecting the non-Boussinesq baroclinic torque $T_{3}$) only when the background shear is absent at the free surface. Since the major objective of this paper is to understand how the free surface impacts the submerged stratified shear instabilities, accurate modeling of the free surface is essential. For this reason, our vertical domain does not end at the air-water interface (free surface) but extends up into the air region. Including the air region may apparently seem redundant but a critical analysis reveals that when shear is present at the free surface, the boundary condition no more remains the well known dynamic boundary condition (i.e. the unsteady Bernoulli’s equation). Thus we treat the free surface as an internal interface between air and water, and implement a free-slip, no penetration boundary up in the air region.

In order to capture the effect of the free surface on the stratified shear instabilities occurring at the pycnocline, we have developed a code that numerically solves the non-Boussinesq Taylor-Goldstein equation. Numerical linear stability analysis is performed on smooth base state velocity and density profiles, respectively given in (9) and (10). First we have studied the simpler case when the air-water interface is replaced by a rigid lid, but the non-dimensional vertical height of the domain, $H/h$, is varied. Our results are in agreement with similar previous studies by Hazel (1972) and Haigh & Lawrence (1999). The maximum growth rate of the KH instability decreases significantly from $0.19$ to $0.096$ on reducing $H/h$ from $10$ to $2$. However, as predicted by Haigh & Lawrence (1999), Holmboe instability is little affected by changes in $H/h$. If rigid-lid approximation is not implemented, i.e. the air-water interface is treated as a non-Boussinesq free surface, the stability boundaries are found to change drastically on decreasing $H/h$. Only for large values of $H/h$, rigid lid becomes a valid approximation since, in this scenario, the free surface has practically no effect on the shear layer. Hence we infer that rigid-lid approximation is somewhat misleading and therefore should be applied with caution.

For higher values of $H/h$, the most unstable mode (for all bulk Richardson numbers) is due to the Rayleigh/KH instability, which arises from the coupling between two vorticity waves, each existing at an extrema of the base vorticity gradient (for broken-line profiles, this would translate to the jump in the base vorticity profile). However, as $H/h$ is decreased, a new mode of instability appears. This mode, which we refer to as the SG-LV mode, becomes the most unstable one, even surpassing KH. In fact, for very low values of $H/h$, KH is nearly non-existent, and the dominant instability is due to SG-LV; see table 1. We note here that similar instability has been observed by Longuet-Higgins (1998) and Bakas & Ioannou (2009) while using broken-line profiles. Holmboe instability is found to be quite resilient to variations in $H/h$, similar to what observed in the rigid-lid case. Hence it can be concluded that rigid lid is a valid approximation even for shallow domains if one is only interested in studying Holmboe instability.

The eigenfunction analysis performed in §3.4 reveals that Holmboe instability is analogous to the baroclinic mode in two-layered flows - the surface and interface are $\pi$ shifted in phase, and furthermore, the surface elevation is insignificant in comparison to that of the interface. This is probably the reason why Holmboe instability is relatively unaffected by variations in $H/h$. The SG-LV mode, on the other hand, is analogous to the barotropic mode in two-layered flows - the surface and interface (pycnocline) are nearly in phase, and the interface elevation is insignificant in comparison to that of the surface. This also implies that the pycnocline plays an insignificant role in this instability.

To complement the numerical stability analysis of the smooth profiles, and furthermore, to obtain a simplified understanding of the instabilities in the system, we have also performed stability analyses of broken-line profiles. A qualitative as well as quantitative agreement between the stability analysis of smooth and broken-line profiles is observed; compare figure 3 with figure 7 and table 1 with table 2. The SG-LV mode is clearly observed as a distinct branch of instability which arises as $H/h$ is decreased (i.e. the depth is made shallower). An additional type of instability is observed in the broken-line profile that is hidden in the analysis of smooth profile (we refer to it as the SG-IG mode). Like SG-LV, the SG-IG instability is also observed when the depth is quite shallow (e.g. $H/h\lesssim 4$). The broken-line profile paved the idea to construct a few reduced order broken-line models so as to underpin the wave interactions that lead to different instabilities. We have constructed three reduced order models, viz. Rayleigh/Kelvin-Helmholtz with free surface, SG-LV and SG-IG; see figure 6.

Through the reduced order SG-LV model we have conclusively shown that the SG-LV mode arises primarily because of the interaction between the leftward propagating surface gravity wave and the rightward propagating vorticity wave, i.e. interaction of waves $2$ and $6$ of figure 1. The understanding is augmented by studying the dispersion diagrams of the full system, as well as the stable isolated waves that are suspected to be behind these instabilities. The crossing of the dispersion curves of the isolated stable waves in the $\alpha-\omega_{r}$ plane corresponds to the unstable region in the $\alpha-\omega_{i}$ plane, confirming our initial assumption. We emphasize here that, although this mode of instability has been observed previously by Longuet-Higgins (1998) and Bakas & Ioannou (2009) (they refer to it as the branch II instability), the fundamental reason behind it was not clearly known. Neither was it known whether such instabilities exist in smooth profiles, and if yes, then how well does it compare with the broken-line counterpart. In this paper we have been able to address all these points conclusively.

An approach similar to SG-LV is taken to understand the SG-IG mode. A reduced order SG-IG model is constructed to underpin the instability mechanism. Finally, with the help of dispersion diagrams, we confirmed that the SG-IG mode is indeed an interaction of the leftward surface gravity wave (wave-$2$) and rightward interfacial gravity wave (wave-$4$). This mode, to the best of our knowledge, has not been reported previously in the literature.

In summary, we have performed a comprehensive study on the effect of free surface on the stratified shear instabilities underneath. For shallow flows, i.e., when the free surface is relatively closer to the shear layer, the free surface significantly affects the ensuing shear instabilities. The surface gravity waves resonate with the different waves present at the shear layer, thereby modifying the “well-known” instabilities, and more importantly, giving rise to two new instabilities. These important dynamics won’t be captured if the air-water interface is modeled as a rigid lid. When shear is present at the free surface, the non-Boussinesq baroclinic torque may become significant and therefore can strongly affect the stability characteristics. Moreover, the non-Boussinesq baroclinic torque is absent when shear is absent at the free surface. Thus, although there is an $\mathcal{O}(1)$ density variation at the free surface, a Boussinesq like approximation is sufficient (that is, one can only consider the gravitational part of the baroclinic torque) when free surface has no background shear. Finally, we point out here that for analytical simplicity, we have considered the pycnocline to be at the mid-depth of the water body. However in real aquatic environments, the pycnocline is usually closer to the free surface. Thus the surface waves can have a more stronger influence on the submerged shear instabilities than that considered in this paper. Future experimental studies and/or Direct numerical simulation (DNS) may be able to shed more light into SG-LV and SG-IG instabilities, especially their non-linear evolution and three dimensional structures.

Appendix A Derivation of Non-Boussinesq Taylor-Goldstein Equation

A full derivation of the non-Boussinesq viscous diffusive Taylor-Goldstein equation is given here. A 2D flow in the $x$-$z$ plane is considered. The horizontal, $x$ component of velocity is given by $u$ while the vertical, $z$ component of velocity is given by $w$. The dynamic viscosity (which is assumed constant) and mass density of fluid are given by $\mu$ and $\rho$ respectively. Acceleration due to gravity is given by $g$, while $t$ denotes time. The governing equations of the problem are given below.

Incompressible continuity equation:

[TABLE]

Navier-Stokes equation for the $x$-momentum:

[TABLE]

Navier-Stokes equation for the $z$-momentum:

[TABLE]

Advection diffusion of the stratifying agent:

[TABLE]

Equation of state relating stratifying agent to density:

[TABLE]

In equations (33) and (34), $\theta$ is the stratifying agent like temperature or salinity, $\eta$ is the (constant) molecular diffusivity of the stratifying agent, while $\beta$ is the linear coefficient relating density to changes in stratifying agent. The quantities $\rho_{0}$ and $\theta_{0}$ are the base state density and some reference value of the stratifying agent concentration. Combining equations (33) and (34) we obtain an advection diffusion equation for mass density given by

[TABLE]

We assume a base state that depends only on $z$, and is given by $u=\bar{u}(z)$, $w=0$, $p=\bar{p}(z)$ and $\rho=\bar{\rho}(z)$. The base state is also assumed to be under hydrostatic balance: $d\bar{p}/dz=-\bar{\rho}g$. Infinitesimal perturbations, denoted by $\tilde{f}$ (where $f$ is a placeholder variable), are added to the base state and then substituted in (30)-(32), and (35). These equations after linearization yields

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We have assumed perturbations of the normal mode form $\tilde{f}(x,z,t)=\hat{f}(z){\rm e}^{{\rm i}\alpha(x-ct)}$, where $\alpha$ is the wavenumber and $c=c_{r}+{\rm i}c_{i}$ is the complex phase speed. Such form on substitution in (36)-(39) yields

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here ′ denotes $d/dz$. Making the substitution $\hat{u}={\rm i}\hat{w}^{\prime}/\alpha$ we get

[TABLE]

Taking the total derivative of the above equation with respect to $z$ we obtain,

[TABLE]

Expressing $\hat{p}^{\prime}$ in terms of other variables from (42) we get

[TABLE]

$\hat{p}^{\prime}$ from the previous equation can be substituted in (44) to give

[TABLE]

[TABLE]

Equations (45)-(46) form the non-Boussinesq viscous diffusive Taylor-Goldstein equations. Inviscid limit gives rise to $\mu\rightarrow 0$, while the non-diffusive limit yields $\kappa\rightarrow 0$. In the limiting condition of inviscid and non-diffusive flow, we obtain

[TABLE]

[TABLE]

Combining (47) and (48) we get

[TABLE]

The above equation is the non-Boussinesq Taylor-Goldstein equation (inviscid and non-diffusive limit of (45)-(46)), and is same as the one obtained by Barros & Choi (2011), Barros & Choi (2014) and Carpenter et al. (2017).