2 Equations

In this paper, we consider only two-dimensional perturbations to the background flow. This is a common assumption, by appealing to the results of Squire (Reference Squire1933) and Yih (Reference Yih1955), who showed that any three-dimensional normal mode can be associated with a two-dimensional one with smaller Richardson number and larger Reynolds number. However, this is not necessarily sufficient to show that the fastest growing mode is always a two-dimensional one (Smyth & Peltier Reference Smyth and Peltier1990). We discuss this further in § 5.

An infinitesimal normal-mode perturbation with vertical velocity $w(x,z,t)={\hat{w}}(z)\text{e}^{\text{i}k(x-ct)}$ to an inviscid Boussinesq flow with velocity profile $U(z)$ and buoyancy profile $B(z)$ must satisfy the well-known Taylor–Goldstein equation

(2.1)$$\begin{eqnarray}(U-c)(\unicode[STIX]{x2202}_{z}^{2}-k^{2}){\hat{w}}-U_{zz}{\hat{w}}=-\frac{B_{z}}{U-c}{\hat{w}}.\end{eqnarray}$$

Here, $k$ is the streamwise wavenumber of the perturbation, and $c=c_{r}+\text{i}c_{i}$ is the complex phase speed, so that the growth rate of a disturbance is given by $\unicode[STIX]{x1D70E}=kc_{i}$.

When kinematic viscosity $\unicode[STIX]{x1D708}^{\ast }$ and diffusivity of the buoyancy distribution $\unicode[STIX]{x1D705}^{\ast }$ are taken into account, (2.1) becomes the more complicated pair of equations

(2.2)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (U-c)(\unicode[STIX]{x2202}_{z}^{2}-k^{2}){\hat{w}}-U_{zz}{\hat{w}}=\text{i}k\hat{b}+\frac{1}{\text{i}k}\frac{1}{Re}(\unicode[STIX]{x2202}_{z}^{2}-k^{2})^{2}{\hat{w}},\\ \displaystyle (U-c)\hat{b}+\frac{1}{\text{i}k}B_{z}{\hat{w}}=\frac{1}{\text{i}k}\frac{1}{PrRe}(\unicode[STIX]{x2202}_{z}^{2}-k^{2})\hat{b},\quad Re\equiv \frac{U_{0}^{\ast }d_{0}^{\ast }}{\unicode[STIX]{x1D708}^{\ast }};Pr\equiv \frac{\unicode[STIX]{x1D708}^{\ast }}{\unicode[STIX]{x1D705}^{\ast }},\end{array}\right\}\end{eqnarray}$$

where length scales and time scales have been non-dimensionalised using the half-depth $d_{0}^{\ast }$ of the shear layer, and half the velocity difference $U_{0}^{\ast }$ across the shear layer, leading to a conventional definition of the Reynolds number, $Re$, and $Pr$ is the usual Prandtl number.

Following Hazel (Reference Hazel1972) and many subsequent authors, we consider the ‘Hazel’ model for the background velocity and buoyancy distributions

(2.3a-d)$$\begin{eqnarray}U(z)=\tanh z,\quad B(z)=\frac{J}{R}\tanh Rz;\quad R\equiv \frac{d_{0}^{\ast }}{\unicode[STIX]{x1D6FF}_{0}^{\ast }};\quad J\equiv \frac{B_{0}^{\ast }d_{0}^{\ast }}{U_{0}^{\ast 2}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{0}^{\ast }$ is the (dimensional) half-depth of the background buoyancy layer with half-difference $B_{0}^{\ast }$ and $J$ is the bulk Richardson number. This is an extension of the Holmboe model (Holmboe Reference Holmboe1960), which has $R=1$ and is attractive because the stability boundary can be found analytically (Miles Reference Miles1961). It is close to the self-similar error function profile expected for a diffusing stratified shear layer when $Pr=R^{2}$ (Smyth, Klaassen & Peltier Reference Smyth, Klaassen and Peltier1988). It is important to note that these profiles are not steady solutions of the viscous Boussinesq equations, but we make the ‘frozen flow’ approximation (Smyth & Carpenter Reference Smyth and Carpenter2019) which is valid when $\unicode[STIX]{x1D70E}\gg 1/Re$. This inequality is not always satisfied by the instabilities we find, as discussed in § 3.4.

The gradient Richardson number $Ri_{g}$, defined as

(2.4)$$\begin{eqnarray}Ri_{g}(z)\equiv \frac{\text{d}B/\text{d}z}{(\text{d}U/\text{d}z)^{2}}=J\frac{\operatorname{sech}^{2}Rz}{\operatorname{sech}^{4}z},\end{eqnarray}$$

for the Hazel model flow, which means (for this particular flow) that at the centreline, $Ri_{g}(0)=J$. For $R\leqslant \sqrt{2}$, $J$ is the minimum of $Ri_{g}$, for $\sqrt{2}<R<2$, there are two minima of $Ri_{g}<J$ either side of local maximum $z=0$, and for $R\geqslant 2$, $Ri_{g}\rightarrow 0$ as $z\rightarrow \infty$ and $z=0$ is a global maximum (Alexakis Reference Alexakis2005). From the Miles–Howard theorem, we then deduce that inviscid HWI at arbitrarily large $J$ is only possible when $R\geqslant 2$. In fact, Alexakis (Reference Alexakis2007) showed that HWI only exists at all for $R\geqslant 2$, despite the possibility of instability at $J>1/4$ when $\sqrt{2}<R<2$.

The solution of (2.2) is performed using a MATLAB code from Smyth & Carpenter (Reference Smyth and Carpenter2019). The method is to construct a large matrix eigenvalue problem, using evenly spaced finite differences. This is a mature code, and additionally the existence of viscous Holmboe was confirmed in direct numerical simulation (DNS) of the Boussinesq equations at finite $Re$ and $Pr=1$ (Parker et al. Reference Parker, Caulfield and Kerswell2019). The boundary conditions are that $\unicode[STIX]{x2202}{\hat{w}}/\unicode[STIX]{x2202}z=\unicode[STIX]{x2202}\hat{b}/\unicode[STIX]{x2202}z=0$, i.e. frictionless, insulating boundaries, at $z=\pm L_{z}$, although all of the instabilities we discuss here are centred around the shear layer, and changing the boundary conditions would not qualitatively affect the results. All linear stability results are found using $768$ finite difference points in the vertical direction, except for the $L_{z}=20$ case which used $1024$ points, and the figures are generated from a $48\times 48$ grid of calculated growth rates.

3 Linear stability analyses

Figure 1 shows a typical example of the VHI. There is a clear distinction between those unstable modes with zero phase speed, which we identify as KHI, and the modes with non-zero phase speed, which we identify as VHI. Although the existence of unstable modes at $R=1$ with non-zero phase speed was unknown before Parker et al. (Reference Parker, Caulfield and Kerswell2019), the diagram bears a striking resemblance to the classic stability diagrams for inviscid HWI for a piecewise linear profile with a density discontinuity (Holmboe Reference Holmboe1962, figure 7) and the Holmboe model with $R>2$ (Hazel Reference Hazel1972, figure 8). Crucially, above $J=0.25$ on this diagram, the gradient Richardson number of the flow is everywhere greater than one quarter, and so we expect stability as $Re\rightarrow \infty$. In the inviscid case, as $J$ is increased the dominant KHI mode and a subdominant KHI mode converge and bifurcate into the pair of HWI modes, with opposite phase speeds. In the viscous case, the regions of the two instabilities overlap slightly and there is no clean bifurcation from one to the other. The remainder of this section will explore how the structure of stability diagrams like figure 1 change as various parameters are varied.

Figure 1. Stability diagram for the Hazel model flow profile as defined in (2.3) with $R=1$, $Re=500$, $Pr=1$, with boundaries at $z=\pm L_{z}=\pm 15$. The contours show the growth rate of two-dimensional normal-mode perturbations of wavenumber $k$, at bulk Richardson number $J$. The colours show the phase speed. The lower region, up to $J=0.25$, is KHI with zero phase speed. The upper lobe is the VHI, with non-zero phase speed. The dashed line shows the analytic stability boundary $J=k(1-k)$ for an unbounded domain in the inviscid limit (Miles Reference Miles1961). In this, and all the stability diagrams in this paper, a waviness is apparent near stability boundaries. This is a common problem in such stability diagrams (Hogg & Ivey Reference Hogg and Ivey2003; Smyth & Winters Reference Smyth and Winters2003; Carpenter et al. Reference Carpenter, Balmforth and Lawrence2010, Reference Carpenter, Tedford, Heifetz and Lawrence2013), and is associated with interpolating near sharp changes of gradient in contour plots.

Figure 2. Vorticity field for the most unstable VHI mode for $R=1$ ($J=0.2128$, $k=0.1042$, a) and $R=3$ ($J=0.8085$, $k=0.5208$, b). In the latter case, a critical layer exists at $z=0.63437$ and is marked with a dashed line.

Table 1. The various parameters used for the linear stability diagrams, as well as the maximum growth rate $\unicode[STIX]{x1D70E}^{max}$ of VHI for each set of parameters, and the phase speed $c_{r}^{max}$, wavenumber $k^{max}$ and bulk Richardson number $J^{max}$ at which they occur.

Figure 2 shows typical eigenmodes of the spanwise vorticity. With $R=1$, i.e. a density interface as wide as the shear layer, no critical layer exists. With $R=3$, the eigenmode is virtually indistinguishable from the $Re\rightarrow \infty$ case, and a critical layer is present and clearly manifests itself within the spatial structure of the mode. Both of these modes have an equivalent mode associated with the complex conjugate eigenvalue, which is identical except for a reflection in the centreline. In the $R=1$ case, we also note that the growth rate is maximised at a much lower wavenumber.

Table 1 shows the full list of parameters for which stability diagrams were produced. For each diagram, we find the maximum growth rate for VHI, i.e. the maximum of $\unicode[STIX]{x1D70E}$ such that the phase speed $c_{r}$ is non-zero, maximised over the discretised values of $k$ and $J$. Since the grids are relatively coarse, the values will not be the true maxima as no optimisation algorithm has been employed, but they give a strong indication of the trend.

3.1 Effects of domain height

The instabilities we study, KHI and HWI, were originally derived as solutions to the Taylor–Goldstein equation in an unbounded domain. There are several ways to approximate a domain of infinite height numerically, but we choose the simplest, which is to use a domain of sufficiently large, but finite, height. How large is sufficient is an important question, as a very large domain is computationally inefficient. Certainly as the height gets small compared with the wavelength of the instabilities we expect the results to change dramatically, and Hazel (Reference Hazel1972) noted how the diagrams always differ from the analytic, unbounded results at low wavenumbers. Figure 3 shows the same diagram as figure 1, but at a smaller (figure 3a) and a larger (figure 3b) domain height. Although the results are slightly different, qualitatively they are very similar, especially for $L_{z}=20$, with $L_{z}=10$ showing more instability at low wavenumbers. The maximum growth rate of the VHI region is $\unicode[STIX]{x1D70E}=1.9489\times 10^{-3}$ for $L_{z}=10$ and $\unicode[STIX]{x1D70E}=2.0869\times 10^{-3}$ for $L_{z}=20$, compared with $\unicode[STIX]{x1D70E}=2.0310\times 10^{-3}$ for $L_{z}=15$, suggesting that $L_{z}=15$ is sufficient to capture the behaviour in which we are interested.

Figure 3. As for figure 1, but with $L_{z}=10$ (a) and $L_{z}=20$ (b).

3.2 Effects of Prandtl number

Figure 4 shows the effect on the stability diagram of varying the Prandtl number. For $Pr=0.7$ (characteristic of thermally stratified air), we find a maximum growth rate of $\unicode[STIX]{x1D70E}=2.7958\times 10^{-3}$, and for $Pr=7$ (a typical value for thermally stratified water) of $\unicode[STIX]{x1D70E}=5.6834\times 10^{-4}$, compared with $\unicode[STIX]{x1D70E}=2.0310\times 10^{-3}$ for $Pr=1$. Therefore, decreasing the diffusion of buoyancy seems to have a stabilising effect on VHI. In contrast, the KHI at the bottom of the figure is virtually unchanged as $Pr$ is varied by an order of magnitude, which reinforces the idea that KHI is produced by the shear alone. Jones (Reference Jones1977) found strong instability at very low $Pr$, but we believe this to be a different effect.

Figure 4. As for figure 1, but with $Pr=0.7$ (a) and $Pr=7$ (b).

Henceforth, we give results using $Pr=R^{2}$, as proposed by Smyth et al. (Reference Smyth, Klaassen and Peltier1988). Despite the fact that the instability seems to be destabilised when $Pr$ is reduced, it is also stabilised when $R$ is decreased, as we shall see.

3.3 Effects of $R$

So far, the results we have presented have concentrated on $R=1$, the original Holmboe model. However, in the inviscid limit, HWI exists only for $R>2$ (Alexakis Reference Alexakis2007). Figure 5 shows the stability diagram at $Re=500$ over a range of $R$, with $Pr=R^{2}$. All diagrams show a region of instability with non-zero phase speed, which we identify as VHI. In the case $R=3$, the diagram is very similar to the classical diagram of an inviscid fluid (Hazel Reference Hazel1972). The unstable region above the usual band, at low wavenumbers, has $c_{r}>1$, so there is no critical layer. As $R\rightarrow 2$ from above, the inviscid results suggest that the band should narrow to a line (Alexakis Reference Alexakis2005), but instead we see a significant region of instability. In the diagrams for $R=1.5$ and $R=2$, a second band of instability is observed above the first, with reduced phase speed, and we conjecture that this may be connected with the higher Holmboe modes. This has not been investigated further, as the growth rate here is vanishingly small.

Figure 5. As for figure 1, but with (a) $R=0.5$, $Pr=0.25$, (b) $R=1.5$, $Pr=2.25$, (c) $R=2$, $Pr=4$, (d) $R=3$, $Pr=9$. Only the last of these would exhibit HWI at $Re=\infty$.

In all cases, although it is not clear from the truncated diagrams, the instability is suppressed at large $k$ by viscosity. This is in contrast to the inviscid limit, which has instability at arbitrarily large $k$ and $J$. It is only in this large $k$ limit that the wave interaction arguments can be made rigorous.

3.4 Effects of Reynolds number

The Miles–Howard theorem tells us that VHI at $R=1$ must disappear for $J>1/4$, in the inviscid limit $Re\rightarrow \infty$. This leaves many possibilities: (i) the region of instability could retreat below $J=1/4$; (ii) the region could shrink; or (iii) the growth rates could vanish but the region remain a constant size. There may or may not be some finite $Re$ above which VHI does not exist. It is also important to ask at what value of $Re$ the growth of the instability is the fastest, or indeed the relative growth rate compared with the diffusion of the background profile.

The growth rate is maximised between $Re=25$ and $Re=30$, with value $\unicode[STIX]{x1D70E}\approx 0.0105$. The relative growth rate $\unicode[STIX]{x1D70E}Re$, which is required to be large compared with unity for a physically relevant instability, was found to increase with $Re$ (at least until $Re=10\,000$), which is a curious result, since it means that despite the growth rate being maximised at a very low value of $Re$, in practice we are more likely to observe the instability at much higher $Re$.

Figure 6. (a) The growth rate $\unicode[STIX]{x1D70E}$ (left axis) and relative growth rate $\unicode[STIX]{x1D70E}Re$ (right axis), maximised over $k$ and $J$, for VHI at $R=1$, $Pr=1$, as $Re$ varies. (b) Growth rate against $J$ and $Re$, maximised over $k$, for $R=1$ and $Pr=1$. The band at the bottom of the figure is KHI, destabilised as $Re$ increases. The upper region with $J\gtrsim 1/4$ is VHI, clearly stabilised as $Re$ increases.

The critical Reynolds number $Re_{c}$ for the viscous Holmboe instability at $R=1$, $Pr=1$, below which there is no instability except KHI, was found to be $Re=4.615$. At criticality, the instability appears at $J=0.25$ and $k=0.12$. This is in contrast with KHI, which at $J=0$ was found to have $Re_{c}=0$ (Betchov & Szewczyk Reference Betchov and Szewczyk1963). In that case, viscosity has a purely stabilising effect. Figure 6 shows how the growth rate varies with $Re$.

3.4.1 Asymptotic behaviour at high $Re$

The fact that some of the VHI do not have a critical layer ($c_{r}^{max}>1$ in Table 1) suggests a regular perturbation analysis may be sufficient to capture the effect of high but finite $Re$. Defining a perturbation parameter $\unicode[STIX]{x1D716}:=1/Re\ll 1$, we may rewrite (2.2) as

(3.1)$$\begin{eqnarray}{\mathcal{L}}_{0}(c)\left(\begin{array}{@{}c@{}}{\hat{w}}\\ \hat{b}\end{array}\right)=\unicode[STIX]{x1D716}{\mathcal{L}}_{1}\left(\begin{array}{@{}c@{}}{\hat{w}}\\ \hat{b}\end{array}\right),\end{eqnarray}$$

where we have defined linear operators

(3.2)$$\begin{eqnarray}{\mathcal{L}}_{0}(c):=\left(\begin{array}{@{}cc@{}}(U-c)(\unicode[STIX]{x2202}_{z}^{2}-k^{2})-U_{zz} & -\text{i}k\\ \displaystyle \frac{1}{\text{i}k}B_{z} & (U-c)\end{array}\right)\end{eqnarray}$$

and

(3.3)$$\begin{eqnarray}{\mathcal{L}}_{1}:=\frac{1}{\text{i}k}\left(\begin{array}{@{}cc@{}}(\unicode[STIX]{x2202}_{z}^{2}-k^{2})^{2} & 0\\ 0 & \displaystyle \frac{1}{Pr}(\unicode[STIX]{x2202}_{z}^{2}-k^{2})\end{array}\right).\end{eqnarray}$$

Adopting the expansions

$$\begin{eqnarray}\displaystyle & \displaystyle \left(\begin{array}{@{}c@{}}{\hat{w}}\\ \hat{b}\end{array}\right)=\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)+\unicode[STIX]{x1D716}\left(\begin{array}{@{}c@{}}{\hat{w}}_{1}\\ \hat{b}_{1}\end{array}\right)+O(\unicode[STIX]{x1D716}^{2}), & \displaystyle \nonumber\\ \displaystyle & \displaystyle c=c_{0}+\unicode[STIX]{x1D716}c_{1}+O(\unicode[STIX]{x1D716}^{2}), & \displaystyle \nonumber\end{eqnarray}$$

(keeping $k$ fixed), then

(3.4)$$\begin{eqnarray}{\mathcal{L}}_{0}(c)={\mathcal{L}}_{0}(c_{0})+\unicode[STIX]{x1D716}c_{1}{\mathcal{L}}_{0}^{\prime }+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

where

(3.5)$$\begin{eqnarray}{\mathcal{L}}_{0}^{\prime }:=\frac{\unicode[STIX]{x2202}{\mathcal{L}}_{0}}{\unicode[STIX]{x2202}c}=\left(\begin{array}{@{}cc@{}}-(\unicode[STIX]{x2202}_{z}^{2}-k^{2}) & 0\\ 0 & -1\end{array}\right).\end{eqnarray}$$

Inserting these expansions into (3.1), we have

(3.6)$$\begin{eqnarray}{\mathcal{L}}_{0}(c_{0})\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)+\unicode[STIX]{x1D716}{\mathcal{L}}_{0}(c_{0})\left(\begin{array}{@{}c@{}}{\hat{w}}_{1}\\ \hat{b}_{1}\end{array}\right)+\unicode[STIX]{x1D716}c_{1}{\mathcal{L}}_{0}^{\prime }\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)=\unicode[STIX]{x1D716}{\mathcal{L}}_{1}\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)+O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$

At leading order, the inviscid Taylor–Goldstein equation is recovered,

(3.7)$$\begin{eqnarray}{\mathcal{L}}_{0}(c_{0})\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)=0,\end{eqnarray}$$

which by the Miles–Howard theorem has only wavelike solutions $c_{0}\in \mathbb{R}$ when $R=1$ and $J>1/4$. In this case, the diagonal elements of ${\mathcal{L}}_{0}(c)$ are real and the off-diagonal elements are purely imaginary, so that for a solution we must have $\arg \hat{b}_{0}=\arg {\hat{w}}_{0}\pm \unicode[STIX]{x03C0}/2$ in the absence of a critical layer. Without loss of generality, we may choose the phase so that ${\hat{w}}_{0}$ is real and $\hat{b}_{0}$ is purely imaginary.

With an inner product on the space of vectors

(3.8)$$\begin{eqnarray}\left\langle \left(\begin{array}{@{}c@{}}w_{1}\\ b_{1}\end{array}\right)\!,\left(\begin{array}{@{}c@{}}w_{2}\\ b_{2}\end{array}\right)\right\rangle :=\int _{L_{z}}^{L_{z}}(w_{1}^{\ast }w_{2}+b_{1}^{\ast }b_{2})\,\text{d}z,\end{eqnarray}$$

we may define the adjoint operator ${\mathcal{L}}_{0}^{\dagger }(c)$ to ${\mathcal{L}}_{0}(c)$ via

(3.9)$$\begin{eqnarray}\left\langle \left(\begin{array}{@{}c@{}}w_{1}\\ b_{1}\end{array}\right),{\mathcal{L}}_{0}(c)\left(\begin{array}{@{}c@{}}w_{2}\\ b_{2}\end{array}\right)\right\rangle =\left\langle {\mathcal{L}}_{0}^{\dagger }(c)\left(\begin{array}{@{}c@{}}w_{1}\\ b_{1}\end{array}\right),\left(\begin{array}{@{}c@{}}w_{2}\\ b_{2}\end{array}\right)\right\rangle ,\end{eqnarray}$$

which gives

(3.10)$$\begin{eqnarray}{\mathcal{L}}_{0}^{\dagger }(c):=\left(\begin{array}{@{}cc@{}}(\unicode[STIX]{x2202}_{z}^{2}-k^{2})(U-c)U_{zz} & \displaystyle -\frac{1}{\text{i}k}B_{z}\\ \text{i}k & (U-c)\end{array}\right).\end{eqnarray}$$

The existence of a non-trivial solution to (3.7) implies a non-trivial solution to

(3.11)$$\begin{eqnarray}{\mathcal{L}}_{0}^{\dagger }(c_{0})\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}^{\dagger }\\ \hat{b}_{0}^{\dagger }\end{array}\right)=0,\end{eqnarray}$$

which is the adjoint eigenfunction. Again, we can choose ${\hat{w}}_{0}^{\dagger }$ to be real and $\hat{b}_{0}^{\dagger }$ to be imaginary.

The $O(\unicode[STIX]{x1D716})$ terms in (3.6) give

(3.12)$$\begin{eqnarray}{\mathcal{L}}_{0}(c_{0})\left(\begin{array}{@{}c@{}}{\hat{w}}_{1}\\ \hat{b}_{1}\end{array}\right)={\mathcal{L}}_{1}\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)-c_{1}{\mathcal{L}}_{0}^{\prime }\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right),\end{eqnarray}$$

so taking the inner product with the adjoint eigenfunction, we have

(3.13)$$\begin{eqnarray}\left\langle \left(\begin{array}{@{}c@{}}{\hat{w}}_{0}^{\dagger }\\ \hat{b}_{0}^{\dagger }\end{array}\right),{\mathcal{L}}_{0}(c_{0})\left(\begin{array}{@{}c@{}}{\hat{w}}_{1}\\ \hat{b}_{1}\end{array}\right)\right\rangle =\left\langle \left(\begin{array}{@{}c@{}}{\hat{w}}_{0}^{\dagger }\\ \hat{b}_{0}^{\dagger }\end{array}\right),{\mathcal{L}}_{1}\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)\right\rangle -c_{1}\left\langle \left(\begin{array}{@{}c@{}}{\hat{w}}_{0}^{\dagger }\\ \hat{b}_{0}^{\dagger }\end{array}\right),{\mathcal{L}}_{0}^{\prime }\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)\right\rangle .\end{eqnarray}$$

The left-hand side of this equation is zero by construction, and so

(3.14)$$\begin{eqnarray}c_{1}=\frac{\displaystyle \left\langle \left(\begin{array}{@{}c@{}}{\hat{w}}_{0}^{\dagger }\\ \hat{b}_{0}^{\dagger }\end{array}\right),{\mathcal{L}}_{1}\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)\right\rangle }{\displaystyle \left\langle \left(\begin{array}{@{}c@{}}{\hat{w}}_{0}^{\dagger }\\ \hat{b}_{0}^{\dagger }\end{array}\right),{\mathcal{L}}_{0}^{\prime }\left(\begin{array}{@{}c@{}}{\hat{w}}_{0}\\ \hat{b}_{0}\end{array}\right)\right\rangle }.\end{eqnarray}$$

Observe from (3.3) and (3.5) that ${\mathcal{L}}_{1}$ and ${\mathcal{L}}_{0}^{\prime }$ are purely imaginary and real respectively. By our choice of phase it is clear that the numerator is therefore imaginary and the denominator real, so $c_{1}$ is purely imaginary. It can be similarly shown that $c_{2}$ is real so the next contribution to the growth rate is at $O(1/Re^{3})$. We further observe that both ${\mathcal{L}}_{1}$ and $c_{1}{\mathcal{L}}_{0}^{\prime }$ are purely imaginary, so from (3.12) we deduce that ${\hat{w}}_{1}$ and $\hat{b}_{1}$ are purely imaginary and real respectively, the opposite situation to ${\hat{w}}_{0}$ and $\hat{b}_{0}$.

We compute $c_{1}$, and thus the growth rate of the instability as $Re$ becomes large $\unicode[STIX]{x1D70E}=-\text{i}kc_{1}/Re+O(1/Re^{3})$, using (3.14) as follows. First, we make an initial guess of $c_{0}$ from the real part of $c$ from a numerical linear stability analysis at $Re=10\,000$. Secondly, we use this approximate $c_{0}$ in the inverse iteration eigenvalue algorithm to solve both (3.7) and (3.11). Finally, we directly evaluate (3.14) using a trapezoidal quadrature rule for the inner products.

The results for two sample parameter values are shown in figure 7(a). Figure 7(b) shows calculated zeroth- and first-order modes. We thus see that the viscous Holmboe modes are a destabilisation of a stable, propagating mode in the inviscid limit. Viscosity acts to break the exact $\unicode[STIX]{x03C0}/2$ phase difference between the vertical velocity and the buoyancy modes. The fact that VHI varies smoothly between modes with and without critical layers, in particular in figure 5(c), suggests this regular perturbation analysis will extend to the case where a critical layer exists, and that the critical layer is not important to the dynamics.

Figure 7. (a) Asymptotic (dashed) and numerical (solid) values of growth rate for VHI at $J=1$, $k=0.5$ (lower) and $J=0.5$, $k=0.25$ (upper). (b) Modes for the $J=0.5$, $k=0.25$ case. The zeroth order modes have been scaled for clarity. In both panels, $R=1$, $Pr=1$ and $L_{z}=15$, corresponding to figure 1.

4 Nonlinear evolution

Smyth & Peltier (Reference Smyth and Peltier1990) showed that at low Reynolds numbers, the linear evolution of HWI is insufficiently fast to overcome the diffusion of the background flow. This leads to the possibility that VHI, for which the growth rates are always small, never physically manifests when the background flow is allowed to diffuse. We consider the nonlinear evolution, which allows us to see whether the viscous Holmboe instability develops the classic counter-propagating vortices of HWI. We use the same DNS code as Parker et al. (Reference Parker, Caulfield and Kerswell2019) to solve the full Boussinesq equations, which is pseudospectral in the streamwise direction and utilises finite differences in the vertical. In the present case, the background flow is allowed to diffuse.

Here we present the results of two direct numerical simulations (restricted to two dimensions) with $R=1.5$, a case for which no HWI is predicted in the inviscid limit. We take $Re=4000$, a compromise between maximising the relative growth rate (see § 3.4) and minimising the spatial resolution. We chose a domain width of $L_{x}=20$, which permits multiple unstable modes. Figure 8 shows the results of a calculation with $J=0.1$, for which we expect a Kelvin–Helmholtz instability to develop to finite amplitude. A linear stability analysis predicts exponential growth rates of $\unicode[STIX]{x1D70E}=0.1244$ and $\unicode[STIX]{x1D70E}=0.0981$ for mode 1 and mode 2 disturbances ($k=\unicode[STIX]{x03C0}/10$ and $k=\unicode[STIX]{x03C0}/5$) respectively, in both cases with zero phase speed. We use a relatively large initial perturbation of random noise in low wavenumber Fourier and Hermite modes, which, along with the comparable growth rates for the two unstable modes, leads to an incoherent, but nevertheless recognisable, Kelvin–Helmholtz billow. At the large $Re$ studied, this rapidly breaks down into turbulence, and significant mixing is achieved, although it is important to remember that this DNS is restricted to two dimensions, and so the specific characteristics of the mixing are likely to be unphysical.

Figure 9 shows the same calculation with $J=0.67$, which maximises the growth rate for VHI at this wavenumber. Again, both modes 1 and 2 are unstable, with growth rates $\unicode[STIX]{x1D70E}=4.1121\times 10^{-4}$ and $\unicode[STIX]{x1D70E}=1.7012\times 10^{-4}$ respectively, and phase velocities $c_{r}=\pm 1.0211$ and $c_{r}=\pm 1.0056$. Since the phase speeds are greater than $1$, no critical layer exists for these instabilities. In this case, the relative growth rate clearly does not satisfy $\unicode[STIX]{x1D70E}Re\gg 1$, so we require a large initial perturbation to trigger significant instability. The strong asymmetry of this random perturbation means that a Holmboe ‘wave’ is apparent only on one side of the interface. Despite the lack of a critical layer, a ‘cusped wave’ very reminiscent of classic HWI (Alexakis Reference Alexakis2009; Salehipour et al. Reference Salehipour, Caulfield and Peltier2016) is apparent, and grows large enough for a clear vortex to be apparent. This vortex is responsible for some mixing, which can be observed when comparing the long time vorticity distribution above the interface, where the vortex exists, to below, where no strong VHI was triggered. However, this mixing is relatively weak compared with the diffusion of the background profile. It is difficult to define a speed precisely for the nonlinear wave, but it appears to be close to $1$. Both the background flow velocity at the level of the vortices and the phase velocity of the linear instability are also approximately equal to $1$. Animations of both evolving flows are available as supplementary materials.

Figure 8. The total vorticity field of a two-dimensional nonlinear simulation of the Boussinesq equations at $Re=4000$, $Pr=2.25$, $L_{x}=20$, $L_{z}=10$ and $J=0.1$. The initial state is a background field with $R=1.5$, plus a perturbation of random noise in the first sixth of the horizontal Fourier modes, and the first five Hermite polynomials in the vertical. Two domain widths are shown horizontally: (a) $t=0$, showing the random initial conditions; (b) $t=20$, showing the Kelvin–Helmholtz billow that has begun to develop; (c) $t=40$, showing that the billow has saturated and is starting to break down; (d) $t=60$, showing that the KHI has led to (two-dimensional) turbulence. An animation of the evolving flow is available as supplementary movies 1 and 2 at https://doi.org/10.1017/jfm.2020.340.

Figure 9. As for figure 8, but with $J=0.4355$: (a) $t=0$, showing the random initial conditions; (b) $t=20$, showing that a ‘cusped wave’ is apparent, characteristic of HWI at finite amplitude; (c) $t=35$, showing that a leftwards-propagating vortex is now visible above the shear layer; (d) $t=110$, showing that the vortex has weakened as the mixing layer diffuses away. An animation of the evolving flow is available as supplementary movies 3 and 4.

5 Discussion and conclusions

In this paper, we have described a new, inherently viscous instability and have demonstrated that it shares many of the characteristic features of the classic, inviscid Holmboe wave instability, namely manifesting as a propagating vortex on either side of the mixing layer and appearing to be caused by the interaction of internal gravity waves on a shear interface. Since it exists in regions of parameter space where no instability is predicted in the inviscid limit, we term it the viscous Holmboe instability, or VHI. The instability we have described is distinct from the ‘viscous Holmboe wave instability’ found by Eaves & Caulfield (Reference Eaves and Caulfield2017) in plane Couette flow, which required non-slip and non-penetration effects in the presence of a rigid boundary, whereas we have shown that boundaries only weakly affect the instability, and the VHI discussed here is truly an instability of a stratified shear layer. Despite the similarities to inviscid HWI, it has significant differences from the classical case: it exists when the density interface is not sharp compared with the shear layer; it can have a phase speed greater than the maximum fluid velocity; and it is destabilised by viscosity. When there is no critical layer, a simple perturbation analysis shows that the VHI arises by viscosity disrupting the perfectly out-of-phase velocity and buoyancy components of the neutrally stable inviscid limit. Our work has made the ‘frozen flow’ approximation that requires $\unicode[STIX]{x1D70E}$ to be large compared with $1/Re$ for the instability to grow quickly compared with the diffusion of the background profile, but we did not find this to be the case. Indeed, our perturbation analysis shows that $\unicode[STIX]{x1D70E}\sim 1/Re$ as $Re\rightarrow \infty$. Numerically, we find that $\unicode[STIX]{x1D70E}$ is small compared with $1/Re$ when $Re\lesssim 10^{2}$, despite the fact that $\unicode[STIX]{x1D70E}$ is maximised for $Re\approx 25$ and only just rises above $1/Re$ for $Re\gtrsim 10^{3}$. This leads to the curious situation that although this is an instability which requires viscosity to exist, the effect of the instability relative to the diffusion of the background flow appears to be greater as $Re$ is increased.

This work is a study of how viscosity affects the Holmboe wave instability as certain parameters are varied. There are many possible extensions which have been examined for (classic) HWI, including considering the effects of compressibility (Witzke, Silvers & Favier Reference Witzke, Silvers and Favier2015), surface tension (Pouliquen, Chomaz & Huerre Reference Pouliquen, Chomaz and Huerre1994), and relaxing the Boussinesq approximation (Umurhan & Heifetz Reference Umurhan and Heifetz2007; Churilov Reference Churilov2019). We briefly investigated the possibility that the higher Holmboe modes described by Alexakis (Reference Alexakis2005, Reference Alexakis2007, Reference Alexakis2009) are also destabilised by viscosity at low $R$, and did indeed find a further band of instability with very small growth rates. Our work has been entirely restricted to two dimensions. Although this is a common assumption when studying linear instabilities of shear flows, there is no physical basis for this, and indeed we would fully expect to see the fastest growing mode being three-dimensional in some regions of parameter space, based on the results of Smyth & Peltier (Reference Smyth and Peltier1990). A third dimension would also significantly affect the nonlinear evolution of the instability at high $Re$.

Despite the lack of a sharp density interface relative to the shear layer for the parameters for which we have found instability, we would certainly still expect internal gravity waves to be present on the interface. There is no reason we are aware of, a priori, to think that these could not resonate with the vorticity waves to cause instability. The wave-resonance descriptions of stratified shear instabilities have been mainly qualitative, except in the cases of piecewise constant density and vorticity profiles, which would be physically inconsistent at finite viscosity. Recent attempts to analyse the components of resonances (Carpenter et al. Reference Carpenter, Balmforth and Lawrence2010; Eaves & Balmforth Reference Eaves and Balmforth2019) and to understand better the dynamics of the resonant system (Heifetz & Guha Reference Heifetz and Guha2018, Reference Heifetz and Guha2019) have relied on analysis which requires perturbations to be inviscid, and these certainly would not apply in the low $Re$ regimes we have described. Though the theory of wave resonance has given useful insight in many situations, it is clearly not the full picture. One major outstanding question is how the Miles–Howard criterion may relate to the wave resonance picture. Baines & Mitsudera (Reference Baines and Mitsudera1994) give an argument from critical layer theory, although the authors themselves admit that this gives neither a necessary nor sufficient criterion for stability.

Most of the unstable regions of the viscous Holmboe instability for $R\leqslant 2$ have $|c_{r}|>1$, so there is no critical layer. Therefore, Lindzen’s wave over-reflection hypothesis for the mechanism of stratified shear instabilities, as well as other interpretations based on the existence of a critical layer, such as the wave–particle interaction described by Churilov (Reference Churilov2019), cannot apply. This is in contrast with the viscous instability described by Miller & Lindzen (Reference Miller and Lindzen1988), in which the viscosity was thought to enable over-reflection at the critical layer. As discussed by Smyth & Peltier (Reference Smyth and Peltier1989), it could be possible that the instability is associated with over-reflection of a wave with a different phase speed, which therefore could itself have a critical layer, but this makes an intuitive explanation much harder. Since the wave over-reflection theory is not a predictive explanation of the instability in this case, it does not seem useful here, although it has certainly proven important in many other circumstances.

Under carefully controlled parameters, we have been able to show significant nonlinear growth of the viscous Holmboe instability at $R=1.5$ and $Re=4000$, from initial noise, leading to secondary instabilities and transition to disorder. This primary instability has no critical layer. Nevertheless, most of the regions of instability we have studied, with $R<2$, have much lower growth rates. We conclude that the viscous Holmboe instability is unlikely to be particularly significant in physical processes. In addition to this, for typical values of Prandtl number in the atmosphere ($Pr\approx 0.7$) we see very small growth rates and for typical values of $Pr\approx 7$ in the oceans, we see the full classical HWI, since in this case $R$ is usually large.

Despite these caveats, we have demonstrated the definite existence of an instability which bears a striking resemblance to HWI, but violates many of the supposed prerequisite conditions. We therefore suggest that any instability in a stratified shear layer be considered Holmboe instability if it manifests as propagating vortices on either side of the shear layer, regardless of the relative width of the density interface, the presence of critical layers or the minimum value of the gradient Richardson number.