3.1 Stability of step-like profiles in the $Re\rightarrow \infty$ , $R\rightarrow \infty$ limit

Before investigating the full linear problem at finite $Re$ , we perform a linear stability analysis on the simpler background density field that consists of step-like layers in the density profile, shown in figure 1(a) with a dashed line, making the further assumptions that both $Re\rightarrow \infty$ and $RePr\rightarrow \infty$ . As reviewed in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011), such profiles allow substantial analytical progress to be made as the eigenvalue problem to determine the growth rate of normal modes reduces to the solution of a polynomial equation. This simplified problem retains many of the key features of the full, finite- $Re$ problem.

We consider normal mode perturbations to the base flow of the form

(3.1a-d ) $$\begin{eqnarray}(\boldsymbol{u},\unicode[STIX]{x1D70C})=[\hat{\boldsymbol{u}}(y),\hat{\unicode[STIX]{x1D70C}}(y)]\exp (\text{i}\unicode[STIX]{x1D6FC}[x-ct]),\quad \unicode[STIX]{x1D6FC}=kh,\quad \unicode[STIX]{x1D70E}=\unicode[STIX]{x1D6FC}c_{i},\quad \boldsymbol{u}=\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D713},\end{eqnarray}$$

where $\unicode[STIX]{x1D713}$ is the streamfunction, $\unicode[STIX]{x1D6FC}$ is the non-dimensional wavenumber, assumed to be real and non-negative without loss of generality, and $c=c_{r}+ic_{i}$ is the (in general) complex phase speed, such that $\unicode[STIX]{x1D70E}$ is the growth rate of the linear normal mode. Considering such infinitesimal perturbations to the base flow defined in (2.9a,b ), the linearized form of (2.5) (with $Re\rightarrow \infty$ and $RePr\rightarrow \infty$ ) can be reduced to the appropriate form of the (inviscid) Taylor–Goldstein equation for the streamfunction eigenfunction $\hat{\unicode[STIX]{x1D713}}(y)$ ,

(3.2a,b ) $$\begin{eqnarray}\frac{\text{d}^{2}\hat{\unicode[STIX]{x1D713}}}{\text{d}y^{2}}-\unicode[STIX]{x1D6FC}^{2}\hat{\unicode[STIX]{x1D713}}+\frac{N^{2}\hat{\unicode[STIX]{x1D713}}}{(y-c)^{2}}=0,\quad N^{2}=-Ri_{B}\frac{\text{d}\bar{\unicode[STIX]{x1D70C}}}{\text{d}y}.\end{eqnarray}$$

Figure 2.

(a) Contours of growth rate $\unicode[STIX]{x1D70E}$ for the inviscid, non-diffusive step-like density profile with dispersion relation given by (3.4a,b ). The dashed line shows the large- $\unicode[STIX]{x1D6FC}$ wave interaction relationship $Ri_{B}=2\unicode[STIX]{x1D6FC}/9$ . (b) Contours of the largest modal growth rate for the viscous fluid with diffusive, smooth density profiles as defined in (3.1a–d ) for $Pr=1$ and $Re=10\,000$ . In both panels, the darkest blue outermost contour has $\unicode[STIX]{x1D70E}=0$ , there is a constant contour interval of $0.005$ . Black lines show the theoretical stability boundary defined by (3.5).

In the limiting case of interest with $R\rightarrow \infty$ , giving discontinuous steps in the density field, the buoyancy frequency reduces to two $\unicode[STIX]{x1D6FF}$ -functions localized at each density interface and $N^{2}=0$ everywhere away from these interfaces. Therefore, $\hat{\unicode[STIX]{x1D713}}$ takes the form of exponential functions in each of the three layers, subject to kinematic and dynamic matching conditions at the two interfaces:

(3.3a,b ) $$\begin{eqnarray}\displaystyle \left[\frac{\hat{\unicode[STIX]{x1D713}}}{y-c}\right]_{-}^{+}=0,\quad \left[(y-c)\frac{\text{d}\hat{\unicode[STIX]{x1D713}}}{\text{d}y}-\hat{\unicode[STIX]{x1D713}}-\frac{Ri_{B}\bar{\unicode[STIX]{x1D70C}}\hat{\unicode[STIX]{x1D713}}}{y-c}\right]_{-}^{+}=0. & & \displaystyle\end{eqnarray}$$

Applying these conditions leads straightforwardly to the dispersion relation for the phase speed

(3.4a,b ) $$\begin{eqnarray}\displaystyle c^{2}=\frac{1}{9}+\frac{Ri_{B}}{2\unicode[STIX]{x1D6FC}_{D}}\pm \sqrt{\frac{2Ri_{B}}{9\unicode[STIX]{x1D6FC}_{D}}+\frac{Ri_{B}^{2}\sinh ^{2}2\unicode[STIX]{x1D6FC}/3}{4\unicode[STIX]{x1D6FC}_{D}^{2}\sinh ^{2}4\unicode[STIX]{x1D6FC}/3}},\quad \unicode[STIX]{x1D6FC}_{D}=\frac{\unicode[STIX]{x1D6FC}}{2}(\coth 4\unicode[STIX]{x1D6FC}/3+\coth 2\unicode[STIX]{x1D6FC}/3). & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The term under the square root is positive for all positive wavenumbers, and so $c^{2}$ is real and any unstable modes must have zero (real) phase speed. The flow is unstable when $c^{2}<0$ , which occurs in the wavenumber interval defined implicitly by

(3.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FC}\sinh 2\unicode[STIX]{x1D6FC}}{18\sinh \unicode[STIX]{x1D6FC}\cosh \unicode[STIX]{x1D6FC}/3\sinh 2\unicode[STIX]{x1D6FC}/3}<Ri_{B}<\frac{\unicode[STIX]{x1D6FC}\sinh 2\unicode[STIX]{x1D6FC}}{18\cosh \unicode[STIX]{x1D6FC}\sinh \unicode[STIX]{x1D6FC}/3\sinh 2\unicode[STIX]{x1D6FC}/3}.\end{eqnarray}$$

We plot the associated contours of growth rate in figure 2(a). Following the original argument of Taylor (Reference Taylor1931), we also plot with a dashed line the predicted long-wave interaction condition for this flow, which is $Ri_{B}=2\unicode[STIX]{x1D6FC}/9$ . This condition corresponds to a Doppler-shifted interaction between the internal gravity waves on the two interfaces, and so clearly shows that this is a TCI. As discussed in detail in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011), such wave interactions can be understood and quantified at all wavenumbers, but here we are principally interested in the classification of the different types of instability, and the large-wavenumber limit is adequate for such classification. The flow is only unstable for $Ri_{B}>1/6$ . In an unbounded flow with constant shear (see for example Sutherland (Reference Sutherland2010) for a further discussion), long waves are linearly unstable for $0<Ri_{B}<1$ , $0<\unicode[STIX]{x1D6FC}\ll 1$ , which is consistent with the linear stability of unstratified inviscid plane Couette flow as demonstrated by Case (Reference Case1960). Due to the presence of the horizontal boundaries in this idealized flow, this long-wave instability is shifted to higher values $1/6<Ri_{B}<1/2$ . The global maximum growth rate across all $\unicode[STIX]{x1D6FC}$ and $Ri_{B}$ is $\unicode[STIX]{x1D70E}=0.085$ , and occurs for $Ri_{B}=0.46$ and $\unicode[STIX]{x1D6FC}=1.75$ .

3.2 Stability of smooth density profiles at finite $Re$ and $Pe$

As discussed in the Introduction, we wish to investigate the nonlinear evolution of the TCI numerically, and so it is necessary to consider initially the linear stability of a flow with smooth density profiles at finite Reynolds number $Re$ and Péclet number $Pe=RePr$ . We therefore consider the linear stability of the base flow defined in (2.9a,b ) at finite Reynolds number with $Pr=1$ and $R=30$ . Linearizing the full (finite $Re$ ) governing equations (2.5) about the base flow (2.9a,b ), and assuming normal mode perturbations of the form (3.1a–d ) we obtain a stratified version of an Orr–Sommerfeld equation coupling the vertical velocity eigenfunction $\hat{v}$ and the density eigenfunction $\hat{\unicode[STIX]{x1D70C}}$ ,

(3.6) $$\begin{eqnarray}\displaystyle -\text{i}\unicode[STIX]{x1D6FC}c\left(\frac{\text{d}^{2}}{\text{d}y^{2}}-\unicode[STIX]{x1D6FC}^{2}\right)\hat{v} & = & \displaystyle -\text{i}\unicode[STIX]{x1D6FC}y\left(\frac{\text{d}^{2}}{\text{d}y^{2}}-\unicode[STIX]{x1D6FC}^{2}\right)\hat{v}+Ri_{B}\unicode[STIX]{x1D6FC}^{2}\hat{\unicode[STIX]{x1D70C}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{Re}\left(\frac{\text{d}^{2}}{\text{d}y^{2}}-\unicode[STIX]{x1D6FC}^{2}\right)^{2}\hat{v},\end{eqnarray}$$

(3.7) $$\begin{eqnarray}-\text{i}\unicode[STIX]{x1D6FC}c\hat{\unicode[STIX]{x1D70C}}=-\text{i}\unicode[STIX]{x1D6FC}y\hat{\unicode[STIX]{x1D70C}}-\hat{v}\frac{\text{d}}{\text{d}y}\bar{\unicode[STIX]{x1D70C}}+\frac{1}{RePr}\left(\frac{\text{d}^{2}}{\text{d}y^{2}}-\unicode[STIX]{x1D6FC}^{2}\right)\hat{\unicode[STIX]{x1D70C}}.\end{eqnarray}$$

Substituting the base flow profile for $\bar{\unicode[STIX]{x1D70C}}$ from (2.9a,b ), these equations define an eigenvalue problem for the complex phase speed $c$ , which can be solved numerically (see Smyth, Moum & Nash (Reference Smyth, Moum and Nash2011) for more details of the numerical procedure used).

In figure 2(b), we plot contours for the growth rate $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D6FC}c_{i}$ as $\unicode[STIX]{x1D6FC}$ and $Ri_{B}$ is varied for a flow with $Re=10\,000$ and $Pr=1$ . We expect the growth rates at this Reynolds number to resemble the inviscid step-like limit to some extent, and so we also plot the stability boundary from figure 2(a). It is apparent that the structure of the stability boundary has changed qualitatively. In particular, analogously to the results considered by Caulfield (Reference Caulfield1994), there are now multiple instability branches. The central branch, which still has the highest growth rates for a given value of $Ri_{B}$ (at least for this choice of middle layer thickness), is the viscous analogue of the TCI discussed in the previous section. However, there are now two other apparent branches of instability, each with non-zero (real) phase speed.

These secondary branches always involve two instabilities with equal growth rate and equal in magnitude but oppositely signed real phase speed. However, instabilities in these branches always have substantially smaller growth rates than the instabilities in the branch which we identify as corresponding to instabilities of TCI type, and so we first discuss in detail the instabilities in this dominant branch. There are several ways in which the finite values of the Reynolds number $Re$ and the Péclet number $Pe=RePr$ appear to have affected the properties of the TCI. There is no longer instability predicted to occur for $\unicode[STIX]{x1D6FC}\rightarrow 0$ , reflecting the singular perturbation associated with the introduction of the no-slip boundary condition at the channel walls and the associated fact that long-wavelength modes are inevitably affected by the boundaries. For small non-zero wavenumbers there appears to be a close connection between the two stability boundaries, although at large wavenumbers, due to the finite thickness of the density interfaces for the flow with $Re=$ 10 000, there is a non-trivial difference.

Figure 3.

Variation of maximum growth rate $\unicode[STIX]{x1D70E}_{m}(Re)$ with $Re$ . Dotted lines mark the zero growth rate and the maximum growth rate $\unicode[STIX]{x1D70E}_{\infty }=0.0907$ for the inviscid limit with smooth density profiles.

We have solved the linear stability eigenvalue problem for a wide range of Reynolds numbers at this Prandtl number $Pr=1$ . In figure 3, we plot against Reynolds number the maximum growth rate over all $Ri_{B}$ and $\unicode[STIX]{x1D6FC}$ at each Reynolds number, defined as

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{m}(Re)=\max _{Ri_{B},\unicode[STIX]{x1D6FC}}[\unicode[STIX]{x1D70E}(Re,Ri_{B},\unicode[STIX]{x1D6FC})].\end{eqnarray}$$

We find that the flow is linearly unstable to two-dimensional normal modes for $Re\geqslant 181$ (i.e. $\unicode[STIX]{x1D70E}_{m}>0$ for $Re\geqslant 181$ ), and that $\unicode[STIX]{x1D70E}_{m}\rightarrow \unicode[STIX]{x1D70E}_{\infty }=0.0907$ as $Re\rightarrow \infty$ . The $Re=10\,000$ flow has a very similar maximum growth rate of $\unicode[STIX]{x1D70E}_{m}=0.0902$ . However, between these two extremes, $\unicode[STIX]{x1D70E}_{m}(Re)$ varies non-monotonically with $Re$ , with a global (across $Re$ at $Pr=1$ ) maximum of $\unicode[STIX]{x1D70E}_{m}(700)=0.137$ .

Figure 4.

Structure of the vertical velocity for the instability with maximum growth rate $\unicode[STIX]{x1D70E}_{m}(Re)$ at $Pr=1$ : (a) $Re=700$ (with $\unicode[STIX]{x1D6FC}=2.65$ ); (b) $Re=10\,000$ (with $\unicode[STIX]{x1D6FC}=1.99$ ). Dashed lines show the location of the density interfaces.

Figure 4 shows the structure of the vertical velocity of two wavelengths of the instability with maximum growth rate at both $Re=700$ and $Re=10\,000$ . In each case the eigenfunction is concentrated about the density interfaces. Unsurprisingly, the eigenfunction is more diffuse for the flow with $Re=700$ , apparently allowing a stronger interaction between the two density interfaces, since the eigenfunctions are not so strongly localized at each interface, even though the associated wavenumber $\unicode[STIX]{x1D6FC}=2.65$ is somewhat larger than that associated with the flow at $Re=10\,000$ , $\unicode[STIX]{x1D6FC}=1.99$ . In addition, for a given amplitude, the eigenfunction associated with the $Re=700$ flow has stronger gradients, yielding an associated vorticity field with larger maxima and minima.

There is an additional consideration in the linear problem at intermediate $Re$ which modifies the spectral properties of the eigenvalue problem at wavenumbers and bulk Richardson numbers close to those associated with the instability with the maximum growth rate $\unicode[STIX]{x1D70E}_{m}(Re)$ . In figure 5, we plot contours of growth rate as a function of $\unicode[STIX]{x1D6FC}$ and $Ri_{B}$ for various $Re$ with $Pr=1$ . It is clear that the low-wavenumber instability branch apparent in figure 2(b) shrinks in size as $Re$ decreases, and indeed completely disappears for sufficiently small $Re$ . As already noted, the growth rates of the globally most unstable TCI vary non-monotonically with Reynolds number, reaching a maximum at $Re=700$ . The position, size and shape of the region of the $\unicode[STIX]{x1D6FC}{-}Ri_{B}$ plane subject to TCI varies with decreasing Reynolds number. The principal effect of this is that the minimum wavenumber below which the system is stable increases, consistent with the increasingly stabilizing effect of the walls on the largest wavelengths as the Reynolds number is decreased. Nevertheless, in broad terms the unstable TCI branch does not change its location substantially in the $\unicode[STIX]{x1D6FC}{-}Ri_{B}$ plane as $Re$ varies.

Figure 5.

Contours of growth rate of the most unstable instability at $Pr=1$ and various values of $Re$ : (a) $300$ ; (b) $400$ ; (c) $700$ ; (d) $1000$ ; (e) $2000$ ; (f) $3000$ ; (g) $5000$ ; (h) 10 000. The central, high-growth-rate branches have zero phase speed and represent TCI. The secondary branches of instabilities have non-zero phase speeds. Contours begin at $\unicode[STIX]{x1D70E}=0$ and have 0.008 increments.

On the other hand, as the Reynolds number is varied, the location of the other branch of instability (apparent at low wavenumbers in figure 2 b) does vary significantly in the $\unicode[STIX]{x1D6FC}{-}Ri_{B}$ plane. This branch clearly moves towards the region of the unstable TCI branch as the Reynolds number is reduced from 10 000, and the two branches begin to overlap. Indeed, at $Re=700$ , which has the global maximum growth rate at $Pr=1$ , this unstable branch straddles the unstable TCI branch. As the Reynolds number is reduced further, this other unstable branch emerges at smaller $\unicode[STIX]{x1D6FC}$ , or equivalently larger $Ri_{B}$ relative to the unstable TCI branch, and the two branches of instability separate once more. At $Re=700$ , the most unstable mode over all $Ri_{B}$ and $\unicode[STIX]{x1D6FC}$ is a TCI occurring at $Ri_{B}=0.61$ and $\unicode[STIX]{x1D6FC}=2.65$ , with growth rate $\unicode[STIX]{x1D70E}_{m}(700)=0.137$ . At these values of $Ri_{B}$ and $\unicode[STIX]{x1D6FC}$ there is also a pair of instabilities with equal and opposite non-zero real phase speed, $c=\pm 0.674$ . It is therefore possible for a primary TCI and another pair of instabilities to be predicted in parallel to be unstable in the same periodic flow domain.

3.3 Dimension of the primary instabilities

Since we observe non-monotonic variation of growth rate of the various instabilities with both $Re$ and $Ri_{B}$ , it is important to investigate whether three-dimensional modes, with a wavenumber vector oriented at an angle to the mean flow, are actually more unstable than the two-dimensional modes considered above. Squire’s theorem (Squire Reference Squire1933) establishes that the primary linear instability of an unstratified viscous fluid is expected to be two-dimensional at a given Reynolds number $Re_{0}$ if

(3.9) $$\begin{eqnarray}\left.\frac{\text{d}}{\text{d}s}[Re\unicode[STIX]{x1D70E}_{max}(Re)]\right|_{Re=Re_{0}}>0,\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{max}$ is the maximum growth rate over $\unicode[STIX]{x1D6FC}$ at fixed Reynolds number and $\text{d}/\text{d}s$ is the derivative in the direction of increasing Reynolds number along the curve of maximum growth rates. An analogous result due to Yih (Reference Yih1955) establishes that the primary linear instability of an inviscid stratified flow is expected to be two-dimensional at a given bulk Richardson number $Ri_{0}$ if

(3.10) $$\begin{eqnarray}\left.\frac{\text{d}}{\text{d}s}\left[\frac{\unicode[STIX]{x1D70E}_{max}(Ri_{B})}{Ri_{B}^{1/2}}\right]\right|_{Ri_{B}=Ri_{0}}<0,\end{eqnarray}$$

where now $\unicode[STIX]{x1D70E}_{max}(Ri_{0})$ is the maximum growth rate over $\unicode[STIX]{x1D6FC}$ at fixed bulk Richardson number, and $\text{d}/\text{d}s$ is the derivative in the direction of increasing bulk Richardson number along the curve of maximum growth rates. For the modes considered here, both of the above conditions are satisfied, and so the primary instabilities are expected to be two-dimensional.

4.1 Simulation P600

Figure 8.

Total energy density (black), kinetic energy density (blue) and potential energy density (red) as a function of time for simulation P600 with $Re=600$ , $Pr=300$ , $Ri_{B}=0.52$ and $L_{x}=6.31$ . Red circles mark the times shown in figure 9. Inset: Time evolution of the total energy (red line) compared to the predicted energy growth of the primary TCI (black line).

Figure 8 shows the evolution of the perturbation energy densities from simulation P600 of the flow with $Re=600$ , $Pr=300$ , $Ri_{B}=0.52$ and $L_{x}=6.31$ . As is apparent in the inset, the initial energy growth is clearly associated with the linear growth rate of the most unstable TCI. The perturbation energy associated with the TCI reaches a peak amplitude at $t\approx 97$ , and subsequently decays whilst slowly oscillating, suggesting strongly that it is not appropriate to describe the TCI as having a ‘saturated’ finite-amplitude manifestation. This behaviour is qualitatively different from the behaviour of flows susceptible to primary KHI, which are known experimentally, numerically and observationally to saturate at finite amplitude as periodic arrays of the classical Kelvin–Helmholtz elliptical ‘billows’.

In this flow, soon after the decay of perturbation energy from its peak value, rapid oscillations appear, as $E$ once again increases, reaching a local maximum at $t\approx 160$ . This local maximum is followed by a sequence of rapid oscillations of period $t_{fast}=4.8\pm 0.1$ . This is very close to the predicted period of oscillation for the primary VHWI of $t_{HWI}\approx 4.7$ . During this time the flow is modulated by a slow oscillation of period $t_{slow}\approx 56$ until the energy reaches its global maximum around $t\approx 320$ , when it then decays at a substantially slower rate. Although the correspondence is not so clear as with the primary TCI, the oscillatory growth of the perturbation energy may be associated with the growth and eventual appearance of the primary VHWI predicted at these parameter values. This suggests that it is indeed possible for qualitatively different primary instabilities to develop in parallel in a stratified shear flow. However, it is important to appreciate that there are undoubtedly profoundly nonlinear dynamics occurring within this flow as well, as for $t\gtrsim 270$ the rapid oscillation transitions to another oscillation with a higher frequency, giving a period $t_{new}=2.8\pm 0.1$ . The system ultimately settles into this higher-frequency, higher-amplitude oscillation for $t\gtrsim 320$ .

Figure 9.

Snapshots of the P600 flow evolution at the red circles shown in figure 8 for $t=77.6$ , 103.6, 119.6, 233.6, 399.6, and 401.6 (a–f, respectively). (Left) Total density field, (middle) perturbation vorticity field and (right) mean horizontal velocity (red), mean density (blue) and notional gradient Richardson number (black). Videos of the flow evolution are available as supplementary material.

Figure 9 shows snapshots of the simulation P600 flow at different points during the evolution of the system, enabling us to interpret the observed dynamics of the energy in terms of flow structures. At each time, we plot the total density distribution $\unicode[STIX]{x1D70C}_{tot}(x,y,t)$ , (capturing the underlying three-layer structure) and the perturbation vorticity distribution, $\unicode[STIX]{x1D714}(x,y,t)$ as defined in (2.8), which is the vorticity in the $z$ -direction after the subtraction of the constant (unit in this non-dimensionalization) associated with the laminar background velocity $\boldsymbol{U}=y\hat{\boldsymbol{x}}$ . We also plot the streamwise averages of the density distribution $\overline{\unicode[STIX]{x1D70C}}(y,t)$ and the streamwise velocity $\overline{u}(y,t)$ , where

(4.4a,b ) $$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D70C}}(y,t)=\frac{1}{L_{x}}\int _{0}^{L_{x}}\unicode[STIX]{x1D70C}_{tot}\,\text{d}x,\quad \overline{u}(y,t)=\frac{1}{L_{x}}\int _{0}^{L_{x}}u\,\text{d}x. & & \displaystyle\end{eqnarray}$$

In addition, we plot a notional gradient Richardson number $Ri_{n}(y,t)$ from these horizontally averaged fields defined as

(4.5) $$\begin{eqnarray}\displaystyle Ri_{n}(y,t)=-Ri_{B}\frac{{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}y}}}{\left({\displaystyle \frac{\unicode[STIX]{x2202}\overline{u}}{\unicode[STIX]{x2202}y}}\right)^{2}}, & & \displaystyle\end{eqnarray}$$

which proves to be a useful measure for the extent to which nonlinear processes modify the linear stability properties of the flow. In particular, snapshots of $Ri_{n}(y,t)$ allow us to identify whether the fundamental three-layer structure of the flow subject to a larger scale, close to linear shear, persists or is substantially modified by nonlinear processes. Videos of the flow evolution are available as supplementary material at https://doi.org/10.1017/jfm.2016.686.

At early times, as is apparent in figure 9(a), the two density interfaces form inward-facing (i.e. towards the midpoint of the flow) cusps and are drawn towards each other at these points. This is typical of the nonlinear manifestation of the TCI observed experimentally by  Caulfield et al. (Reference Caulfield, Peltier, Yoshida and Ohtani1995) and the previous two-dimensional numerical simulations of Lee & Caulfield (Reference Lee and Caulfield2001). These cusps are naturally induced by the (positive) vorticity redistributing into elliptical vortices, spanning the entire depth of the middle layer, a characteristic finite-amplitude manifestation of the TCI. However, the induced cusps also lead to baroclinic vorticity generation, which is strongly negative at the interfaces’ cusps. Unlike the KHI at finite amplitude, it is clear from figure 9(b,c) that the TCI consists of a coupled vortical structure, with a central (positive) vortical core flanked by two negative vortical wings. We also see in figure 9(b) that four vortical cores have formed from the mode-two primary TCI, and in 9(c) that two of these cores have grown to have larger streamwise extent than the other two cores. The formation of an extra pair of vortical cores is due to the attempted pinching together of the cusps on opposite density interfaces becoming misaligned due to the background shear, and so inducing two extra cusps on each interface.

The vortical cores appear to be tilted into the base shear, and thus have the potential to extract energy from the base flow transiently through the Orr mechanism (Orr Reference Orr1907). However, this tilting also induces modifications of the density field, leading to annihilation of the different-signed vorticity distributions, and hence the primary TCI instability is almost entirely sheared out by the background flow, as is apparent in figure 9(d). The shearing out process has also lessened the extent of the asymmetry between the alternating large and small streamwise extent of the vortical cores. Indeed, the spatial periodicity of the primary TCI is typically rapidly destroyed by the secondary instabilities, and the mode-two primary TCI is only briefly apparent within the flow.

It is apparent in figure 9(d), showing the horizontally averaged density and velocity profiles and the associated notional gradient Richardson number $Ri_{n}(y,t)$ , that the growth (and decay) of the primary TCI does not modify the base flow qualitatively, in that there remains a clear bimodal peak in the Richardson number, associated with a smooth variation of velocity across the entire flow, with two relatively ‘sharp’ interfaces persisting to separate three close-to-constant density layers.

During shearing out of the remnants of the nonlinear evolution of the TCI, mode-two Holmboe waves appear above the upper density interface and below the lower density interface. We believe that these waves are the eventual nonlinear appearance of the primary VHWI present in the original flow profile, which has a substantially smaller growth rate $\unicode[STIX]{x1D70E}_{HWI}(600)=0.020$ than the primary TCI which has growth rate $\unicode[STIX]{x1D70E}_{m}(600)=0.101$ . As already noted, the observed period of oscillation of the perturbation energy is close to the predicted period of oscillation of this primary VHWI. Furthermore, as is apparent in the perturbation vorticity field, particularly in the video available as supplementary material, the appearance of this instability is closely correlated with the appearance of propagating vorticity perturbations localized near the channel walls. It is important to appreciate that we do not initialize the flow with the modal form of this instability, but rather only with the modal form of the primary TCI and small-amplitude noise. Due to the combination of this different initialization, and the relative growth rates of the two primary instabilities, it is unsurprising that this primary instability reaches finite amplitude substantially later than the TCI.

Before the vortices are sheared out fully, the primary VHWI has grown sufficiently to produce large-amplitude nonlinear Holmboe waves which travel along the density interfaces. These nonlinear Holmboe waves prevent the vortices from developing any subsequent streamwise asymmetry, and they remain approximately equally spaced throughout the remaining evolution. The nonlinear Holmboe waves reach a peak total energy at $t\approx 320$ , before slowly decaying. The characteristic feature of a (one-sided) nonlinear Holmboe wave is a patch of vorticity (in this case with positive sign) above the density interface that pulls wisps of fluid from the layer on the other side of the interface into the layer where the vorticity patch is. In this system, the nonlinear Holmboe waves that propagate consist of two wavelengths (i.e. they are ‘mode-two’ instabilities), and there are two patches of vorticity at each interface, as can be seen in figure 9(e,f). These patches grow in intensity and size over time, interacting with the perturbation vorticity localized near to each of the channel walls. Indeed this nonlinear interaction appears to mediate the continued growth of the nonlinear Holmboe waves and initiates the slow decay for $t\gtrsim 320$ , and is specific to PCF at finite $Re$ due to the non-slip boundary condition at the channel walls.

It is clear that a primary (albeit inherently viscous) HWI of the initial base flow has grown in ‘parallel’ to the primary TCI, even though its growth rate is appreciably smaller. Therefore, as hypothesized in Caulfield et al. (Reference Caulfield, Peltier, Yoshida and Ohtani1995), it is possible for multiple primary instabilities to grow to finite amplitude in layered stratified shear flows, with different growth rates and phase speeds (and hence critical layer locations). Furthermore, here the HWI is apparently more robust than the TCI, calling into question the conventional assumption that the most unstable mode of linear theory will continue to dominate at finite amplitude as $t\rightarrow \infty$ .

4.2 Simulation S5000

Figure 10.

Total energy density (black), kinetic energy density (blue) and potential energy density (red) as a function of time for the simulation with $Re=5000$ and $Pr=70$ . Red circles mark the times shown in figure 11. Inset: Time evolution of the total energy (red line) compared to the predicted energy growth of the primary TCI (black line).

Figure 11.

Snapshots of the S5000 flow evolution at the red circles shown in figure 10 for $t=70$ , 85, 95, 100, 125, 140, 200, 380 and 400 (a–i, respectively). (Left) Total density, (middle) perturbation vorticity field and (right) mean horizontal velocity (red), mean density (blue) and notional gradient Richardson number (black). Videos of the flow evolution are available as supplementary material.

Turning our attention to the higher-Reynolds-number flow, figure 10 shows the evolution of the perturbation energy densities from simulation S5000 for the flow with $Re=5000$ , $Pr=70$ , $Ri_{B}=0.51$ and $L_{x}=6.72$ . There is a clear first peak in energy at $t\approx 80$ that results from the growth of the most unstable normal mode, and it is followed by another peak in energy at $t\approx 95$ . After this second maximum, the energy subsequently decays rapidly. There is then an onset of ‘fast’ oscillations, and the energy recovers until a maximum at $t\approx 275$ , setting into a slower oscillation for $t\gtrsim 150$ . After $t\approx 275$ , the flow transitions into a yet slower oscillation coupled with an appreciably longer time-scale nonlinear modulation. This behaviour continues for $t\gtrsim 300$ and is reminiscent of the energetics of the nonlinear Holmboe wave state found by Balmforth et al. (Reference Balmforth, Roy and Caulfield2012). Figure 11 shows snapshots of the flow at different points during the evolution of the system, which we now discuss in detail. Videos of the flow evolution are available as supplementary material.

Figure 12.

(a) Detail of the mean horizontal velocity (red) and mean density (blue) at $t=125$ , taken from figure 11(e). (b) Vertical velocity $v$ for a most unstable mode of these mean profiles.

As in the lower-Reynolds-number flow simulation P600 discussed above, the evolution of the unstable linear TCI results in cusps and streaks of vorticity forming on the density interfaces, pulling them towards each other, as can be seen in figure 11(a). These vorticity streaks then begin to roll up to form two coherent vortical cores, i.e. a primary mode-two instability develops. Similarly to the lower-Reynolds-number flow, the vortical cores begin to lose energy as they are sheared out. However, as the vortices are sheared out, there is once more a misalignment of the cusps on the density interfaces, resulting in the creation of two more vortical cores located between the primary pair of vortical cores, as is evident in figure 11(b). The two newly created vortical cores begin to grow in size, and appear to squeeze the two original vortical cores. This causes complex spatio-temporal vorticity dynamics within the two original vortical cores, as is evident in figure 11(c,d). The squeezing of the two original vortical cores causes a cascade of secondary braid-like instabilities that results in small-scale disorder in the intermediate layer, and a destruction of the cusps on the two density interfaces due to vigorous mixing in the intermediate layer. The disordered vorticity field is evident in figure 11(e).

The mean velocity profile during this disordered state exhibits increased shear over the laminar profile at each of the two density interfaces, as can be clearly seen in figure 12(a), where we show a detailed version of the horizontal mean profiles $\overline{u}$ and $\overline{\unicode[STIX]{x1D70C}}$ found in figure 11(e). This increased shear spontaneously creates eight regions of negative vorticity, four above the upper density interface and four below the lower density interface, as seen in figure 11(f). These vortices also pull wisps of fluid from the intermediate density layer into both the upper and lower density layers, and are clearly associated with the parasitic appearance of a nonlinear Holmboe wave. This HWI is categorically not a primary instability for this flow at this wavelength, and so the appearance of a parasitic nonlinear Holmboe wave has occurred in ‘series’ to the nonlinear development and subsequent decay of the primary TCI. Furthermore, this secondary HWI is of ‘classical’ type, associated with the localized intensified shear in the vicinity of each density interface clearly seen in figure 12, and not the perturbation vorticity localized near the channel walls. The vigorous mixing associated with the braid-instability-driven decay of the evolution of the primary TCI has rearranged the mean fields in just such a way as to cause the flow to become susceptible to the appearance of HWI. Indeed, a linear stability analysis of a nominal flow corresponding to the horizontally averaged mean profiles shown in figure 12(a) predicts unstable HWI modes with both positive and negative phase speed, and wavelength corresponding to 25 % of the horizontal extent of the computational domain. One of these mode-four modes (with positive phase speed) is plotted in figure 12(b). This prediction is consistent with the nonlinear structures which are observed to develop subsequently, although there is of course no a priori reason why the stability properties of the horizontally averaged total flow should correspond precisely to observations when the actual real flow exhibits non-trivial spatial and temporal variability. Soon after the appearance of mode-four nonlinear Holmboe waves there is a coarsening to mode-two nonlinear Holmboe waves, as is apparent from the two vortices above the upper density interface and two vortices below the lower density interface in figure 11(g).

The mode-two nonlinear Holmboe wave state gradually increases in energy during $150\lesssim t\lesssim 275$ . As for the P600 simulation discussed above, the peak in energy at $t\approx 275$ appears to correspond with the mode-two state interacting with the flow boundaries, and this mediates the end of continued energy growth. As the flow begins to decay, there is yet another coarsening event, and the flow settles onto a large-amplitude mode-one (i.e. with the same wavelength as the streamwise length of the entire computational domain) nonlinear Holmboe wave for $t\gtrsim 350$ . This final state, which is visible in figure 11(h,i), is reminiscent of the large-amplitude nonlinear Holmboe wave found by Balmforth et al. (Reference Balmforth, Roy and Caulfield2012), and its appearance is associated with parasitic development on the remnants of the evolution of the primary TCI.

Although such two-dimensional simulations at high Reynolds number should be treated with caution, two key conclusions consistent with the results of Balmforth et al. (Reference Balmforth, Roy and Caulfield2012) can be drawn. First, the TCI is clearly ‘fragile’ and prone to strong secondary two-dimensional instabilities at high Reynolds number. Second, those secondary instabilities can modify the base flow in a way that then leads to the ‘parasitic’ development of HWI in ‘series’ at later times in layered stratified shear flows. This is apparently due to the local intensification of shear at each of the relatively sharp density interfaces, as is apparent in figure 12.