2.1 Holmboe instabilities

Holmboe (Reference Holmboe1962) considered the stability of an idealised two-layer stratification embedded within a shear layer of finite depth, and demonstrated that, for sufficiently strong stratification, such inviscid flows are always linearly unstable over a finite band of wavenumbers. He considered a symmetric configuration, with the density interface located at the mid-point of the shear layer, so there are two normal modes (with equal and opposite phase speeds) with equal growth rates in this instability band of wavenumbers.

This instability, which we refer to as the ‘traditional’ Holmboe instability, can be interpreted physically as being due to a resonance between infinitesimal internal wave-like perturbations localised at the density interface and infinitesimal Rayleigh (or ‘vorticity’) wave-like perturbations localised at the edges of the shear layers (Baines & Mitsudera Reference Baines and Mitsudera1994; Caulfield Reference Caulfield1994). This interpretation is a particular example of the ‘wave interaction theory’ for the identification and classification of instabilities (see Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2013) for a review). The symmetric problem considered by Holmboe (Reference Holmboe1962) has two such resonances simultaneously, leading to two instabilities: one localised in the upper half of the shear layer, the other in the lower half of the shear layer. However, this physical internal wave–vorticity wave interaction instability mechanism is generic and gives rise to instabilities sharing the same essential features as the traditional Holmboe instability in more general stratified shear flows with finite depth shear layers with embedded ‘sharp’ density interfaces. It is this generic class of instabilities that we, for simplicity, refer to as Holmboe instabilities (HIs).

Lawrence, Browand & Redekopp (Reference Lawrence, Browand and Redekopp1991) generalised the original stability calculations of Holmboe to consider asymmetric flows where the sharp density interface and the mid-point of the shear layer are not coincident. This breaks the symmetry of the traditional stability problem, giving rise to two instabilities with different growth rates, phase speeds and wavenumber bands for a given stratification. For flows in domains that are not bounded vertically, the resonance between the infinitesimal internal wave-like perturbation localised at the density interface and the vorticity wave-like perturbation localised at the further edge of shear layer gives rise to the assumed dominant instability band, with larger maximum predicted growth rate, smaller magnitude phase speed and smaller wavenumbers. Haigh & Lawrence (Reference Haigh and Lawrence1999) further demonstrated that vertical confinement by sufficiently close top and bottom boundaries could switch the dominant branch of instability, consistent with the seminal work of Hazel (Reference Hazel1972).

2.2 Holmboe waves: laboratory observations and numerical simulations

Finite-amplitude manifestations of HIs, which we refer to as Holmboe waves (HWs), have been observed experimentally in a wide range of shear flows containing relatively sharp density interfaces (see e.g. Macagno & Rouse Reference Macagno and Rouse1961; Thorpe Reference Thorpe1968; Browand & Winant Reference Browand and Winant1973; Maxworthy & Browand Reference Maxworthy and Browand1975; Caulfield et al. Reference Caulfield, Peltier, Yoshida and Ohtani1995; Zhu & Lawrence Reference Zhu and Lawrence2001; Hogg & Ivey Reference Hogg and Ivey2003; Tedford, Pieters & Lawrence Reference Tedford, Pieters and Lawrence2009; Carpenter et al. Reference Carpenter, Tedford, Rahmani and Lawrence2010; Meyer & Linden Reference Meyer and Linden2014). Generically, cusped waves are observed at the interface, propagating at a phase speed intermediate between the propagation speed of the density interface and the maximum speed of the sheared fluid, in qualitative agreement with linear stability theory.

Zhu & Lawrence (Reference Zhu and Lawrence2001) observed both symmetric (up–down) and asymmetric (only up or only down) HWs in a laboratory exchange flow. They compared their observations with two different normal mode stability calculations. One calculation used the piecewise-linear profiles originally considered by Holmboe (Reference Holmboe1962), while the other used smooth parallel velocity and density profiles, calculated as solutions to the Taylor–Goldstein equation by Nishida & Yoshida (Reference Nishida and Yoshida1990). They found only qualitative agreement with the predicted phase speeds and wavelengths, perhaps unsurprisingly due to the highly idealised nature of the assumed velocity and density profiles compared to the actual laboratory flows. Tedford et al. (Reference Tedford, Pieters and Lawrence2009) and Carpenter et al. (Reference Carpenter, Tedford, Rahmani and Lawrence2010) conducted further experiments, and the latter made the first comparison with three-dimensional direct numerical simulations (DNS). They attributed the discrepancy between experiments and linear theory to wave stretching (increasing wavelength by the streamwise convective acceleration of the mean flow), and the discrepancy between DNS and linear theory to wave merging (increasing wavelength by vortex merging events). They also found that the wave amplitude was larger in the DNS than in the experiments, although the simulations did not correspond exactly to the experimental situation, in particular in the choice of the Schmidt number.

In the first two-dimensional DNS of HWs, Smyth et al. (Reference Smyth, Klaassen and Peltier1988) established that the counter-propagating trains of HWs above and below the (symmetrically located) density interface accelerate when the crests of the characteristic cusped waves are aligned vertically, and slow down when they are laterally far apart. As is apparent in their figure 7, this behaviour is associated with the relative location of the elliptical spanwise vortices that develop in front of the density crests, with the slowing occurring when the two vortices are close to being aligned vertically (corresponding to density crests being far apart). Indeed, for sufficiently weak stratification, these two vortices ‘lock’ and then roll up the interface into an array or ‘train’ of elliptical billows, i.e. the finite-amplitude manifestation of the Rayleigh instability of a finite-depth shear layer (see Rayleigh Reference Rayleigh1879, and for further discussion of this phenomenon for ‘marginal’ HWs, Smyth & Peltier Reference Smyth and Peltier1991). Such structures are commonly referred to as ‘Kelvin–Helmholtz billows’, even though Kelvin (Thomson Reference Thomson1871) and Helmholtz (Reference Helmholtz1878) only considered a discontinuous velocity profile, which does not exhibit the scale selection characteristic of billow trains.

The elliptical vortices of HWs develop close to the critical layers of the linear instability (i.e. where the predicted phase speed of the linear instability matches the flow speed). Due to the relatively strong baroclinic vorticity generation associated with the cusped waves, these vortices do not typically correspond to the maximum vorticity magnitude within the flow, although they typically induce the entrainment of ‘wisps’ of fluid from the crest of the propagating cusped waves.

Smyth & Peltier (Reference Smyth and Peltier1991) used secondary linear Floquet stability analysis, (i.e. the linear stability analysis of non-parallel flow fields extracted from two-dimensional numerical simulations under the assumption that these flow fields were ‘frozen’ in time) in an attempt to understand how a two-dimensional saturated HW becomes three-dimensional, breaks down and hence has the potential to undergo the transition to turbulence. Simulations at moderate Reynolds numbers suggested that HWs could be robust and identifiable over relatively long times (see for example Smyth & Winters Reference Smyth and Winters2003, Carpenter, Lawrence & Smyth Reference Carpenter, Lawrence and Smyth2007 and Smyth, Carpenter & Lawrence Reference Smyth, Carpenter and Lawrence2007). However, through simulations at significantly higher Reynolds number, Salehipour et al. (Reference Salehipour, Caulfield and Peltier2016) demonstrated that finite-amplitude symmetric HWs are prone to a wide range of secondary instabilities of both convective and shear-driven types, which destroy the coherence of the primary instability and trigger a relatively long period (at least in comparison to KH billow trains) of turbulent motions characterised by enhanced turbulent dissipation and irreversible mixing. Nevertheless, these behaviours at higher Reynolds numbers may well be inherently linked to the transient run-down nature of such simulations, with no forcing acting to re-energise the initial shear flow, and the initial profiles and computational domain geometry being chosen to be susceptible to specific (and typically monochromatic) linear instabilities.

2.3 Limitations of previous research

In this paper, we attempt to address some of the limitations of the previous studies of HWs, which we classify in three categories: steadiness; dissipation; and three-dimensionality.

2.3.2 Dissipation

The majority of the experiments and numerical calculations of which we are aware considered flows at relatively low Reynolds numbers of $O(10{-}100)$ (based on the definition (3.4) adopted later in this paper), with the recent exception of Salehipour et al. (Reference Salehipour, Caulfield and Peltier2016). Moreover, the ‘run-down’ character means that a strictly restricted amount of energy, set at the beginning of the flow evolution, is available to dissipate, and that this dissipation is inherently limited in time. Indeed, the stratified turbulence literature (see e.g. the review of Ivey et al. (Reference Ivey, Winters and Koseff2008)) increasingly highlights the importance of the ‘buoyancy Reynolds number’ $Re_{b}$ , defined as

(2.1) $$\begin{eqnarray}\displaystyle Re_{b}\equiv \frac{\unicode[STIX]{x1D716}_{0}}{\unicode[STIX]{x1D708}N_{0}^{2}}. & & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D716}_{0}$ is a characteristic (volume-averaged) value of the rate of dissipation of total kinetic energy $\unicode[STIX]{x1D716}\equiv 2\unicode[STIX]{x1D708}\unicode[STIX]{x1D634}_{ij}\unicode[STIX]{x1D634}_{ij}$ (where $\unicode[STIX]{x1D634}_{ij}\equiv (\unicode[STIX]{x2202}_{x_{i}}u_{j}+\unicode[STIX]{x2202}_{x_{j}}u_{i})/2$ is the symmetric strain rate tensor) and $N_{0}$ is a characteristic value of the buoyancy frequency $N\equiv \sqrt{-(g/\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}}$ . This non-dimensional number, sometimes referred to as the ‘activity parameter’ or $Gn$ (Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), is the (square) ratio of the frequencies associated with kinetic energy dissipation $\sqrt{\unicode[STIX]{x1D716}/\unicode[STIX]{x1D708}}$ and buoyancy $N$ . An alternative and illuminating interpretation is to express $Re_{b}$ as a ratio of length scales: $Re_{b}=(L_{o}/L_{k})^{4/3}$ where $L_{o}=(\unicode[STIX]{x1D716}/N^{3})^{1/2}$ is the Ozmidov scale, the vertical scale below which turbulent eddies are not significantly affected by stratification, and $L_{k}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D716})^{1/4}$ is the Kolmogorov scale, below which viscous dissipation dominates. The ratio $Re_{b}$ therefore measures the range of scales that are not significantly affected either by stratification (suppressing vertical scales ${\gtrsim}L_{o}$ ) or by viscosity (damping scales ${\lesssim}10L_{k}$ ). Although caution needs to be exercised when comparing different studies and flow geometries where $Re_{b}$ may be defined in different ways, consensus is developing, not least motivated by the arguments of Gibson (Reference Gibson and Nihoul1980) and Gibson (Reference Gibson1999), that $Re_{b}\gtrsim 20{-}30$ is required for the flow to exhibit any of the key characteristics of stratified turbulence (see Bartello & Tobias Reference Bartello and Tobias2013; Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016 for further discussion).

In the stratified inclined duct experiment, we estimate the volume-averaged $Re_{b}\approx Re\sin \unicode[STIX]{x1D703}$ (see (3.4) and (3.6)) and CHWs are found to exist for $Re_{b}$ up to ${\approx}50$ , controlled by increasing either the Reynolds number $Re$ or the tilt angle of the duct $\unicode[STIX]{x1D703}$ (see § 3.4.1 for more details). We are therefore able to generate waves that maintain relatively high dissipation for hundreds of advective times units, making this experiment a viable attempt to reproduce inertia-dominated flows of geophysical relevance in the laboratory.

2.3.3 Three-dimensionality

Temporal linear stability analyses of parallel shear flows with a background streamwise velocity $U(z)\hat{\boldsymbol{x}}$ focus on wave-like disturbances, i.e. normal modes, of the form

(2.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F9}^{\prime }(\boldsymbol{x}_{\bot },z,t)\equiv \hat{\unicode[STIX]{x1D6F9}}(z)\exp (\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}_{\bot }+\unicode[STIX]{x1D70E}t), & & \displaystyle\end{eqnarray}$$

where $\hat{\unicode[STIX]{x1D6F9}}(z)$ is the one-dimensional complex eigenfunction of any perturbation variable; $z$ is the vertical coordinate; $\boldsymbol{x}_{\bot }\equiv x\hat{\boldsymbol{x}}+y\hat{\boldsymbol{y}}$ is the horizontal coordinate vector; $\boldsymbol{k}\equiv |\boldsymbol{k}|(\cos \unicode[STIX]{x1D6FD}\hat{\boldsymbol{x}}+\sin \unicode[STIX]{x1D6FD}\hat{\boldsymbol{y}})$ is the real horizontal wave vector; and $\unicode[STIX]{x1D70E}$ is the complex growth rate. Conventionally, $\boldsymbol{k}$ is taken to be aligned with the background flow along $x$ , i.e. $\unicode[STIX]{x1D6FD}=0$ . In this paper, we refer to such an analysis as ‘2P-1B’, for two-dimensional perturbations on a one-dimensional base flow and to the more general analysis with $\unicode[STIX]{x1D6FD}\neq 0$ as ‘2.5P-1B’, for its intermediate character between two- and three-dimensional perturbations. The focus on 2P-1B analyses can generically be justified by appealing to a corollary of ‘Squire’s theorem’ (Squire Reference Squire1933) in unstratified flows at finite Reynolds number, and the analogous ‘Yih’s theorem’ (Yih Reference Yih1955; Smyth et al. Reference Smyth, Klaassen and Peltier1988) in (inviscid) stratified flow.

Both theorems rely on the observation that ‘three-dimensional’ normal modes aligned at some angle $\unicode[STIX]{x1D6FD}\neq 0$ to the background flow, having growth rate $\unicode[STIX]{x1D70E}$ , wavenumber $|\boldsymbol{k}|$ in a flow with Reynolds number $Re$ and bulk Richardson number $Ri$ are equivalent to ‘two-dimensional’ normal modes having $\unicode[STIX]{x1D6FD}=0$ with lower growth rate $(\unicode[STIX]{x1D70E}\cos \unicode[STIX]{x1D6FD})$ and larger wavenumber $(|\boldsymbol{k}|/\text{cos}\,\unicode[STIX]{x1D6FD})$ in an equivalent flow with lower Reynolds number $(Re\cos \unicode[STIX]{x1D6FD})$ and higher bulk Richardson number $(Ri/\text{cos}^{2}\,\unicode[STIX]{x1D6FD})$ . This statement is summarised by the relation (see Smyth & Peltier Reference Smyth and Peltier1990)

(2.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}(|\boldsymbol{k}|,\unicode[STIX]{x1D6FD},Re,Ri)=\cos \unicode[STIX]{x1D6FD}\;\unicode[STIX]{x1D70E}\left(\frac{|\boldsymbol{k}|}{\cos ~\unicode[STIX]{x1D6FD}},0,Re\cos \unicode[STIX]{x1D6FD},\frac{Ri}{\cos ^{2}\unicode[STIX]{x1D6FD}}\right). & & \displaystyle\end{eqnarray}$$

It is then natural to consider the locus of the two-dimensional modes ( $\unicode[STIX]{x1D6FD}=0$ ) with maximum growth rate $\unicode[STIX]{x1D70E}_{m}$ across all wavenumbers for given values of the external flow parameters $Re$ and $Ri$ . For an unstratified flow, if $\unicode[STIX]{x1D70E}_{m}(Re)$ decreases more slowly than $1/Re$ (i.e. $\unicode[STIX]{x1D70E}_{m}(Re)\times Re$ increases with $Re$ ), then the most unstable two-dimensional modes have higher growth rates than any three-dimensional mode, a circumstance which is generic for shear flow instabilities. In particular, if instability occurs only above a critical Reynolds number (i.e. for $Re\geqslant Re_{c}$ ), this condition implies that the marginal instability mode must be two-dimensional (as any three-dimensional mode in a flow with $Re=Re_{c}$ experiences an equivalent flow with $Re<Re_{c}$ , and hence is stable).

The equivalent statement for (inviscid) stratified flows is, if $\unicode[STIX]{x1D70E}_{m}(Ri)$ increases more slowly than $\sqrt{Ri}$ (i.e. $\unicode[STIX]{x1D70E}_{m}(Ri)/\sqrt{Ri}$ decreases with increasing $Ri$ ), then the most unstable two-dimensional modes have higher growth rates than any three-dimensional mode. For stratified flows, the dominance of two-dimensional modes is not so clear cut, as HIs (for example) only occur in stratified flows, and occur for all bulk Richardson numbers, thus making it at least plausible that $\unicode[STIX]{x1D70E}_{m}(Ri)$ increases more rapidly than $\sqrt{Ri}$ . Indeed, in a stratified flow at finite Reynolds number, Smyth & Peltier (Reference Smyth and Peltier1990) identified a small parameter range in $Re{-}Ri$ space where $\unicode[STIX]{x1D70E}_{m}(Ri)$ does increase more rapidly than $\sqrt{Ri}$ and HIs at an angle $\unicode[STIX]{x1D6FD}\neq 0$ were predicted to be dominant, but it has proved difficult to test these predictions either experimentally or numerically.

It is important to recognise that such ‘three-dimensional’ modes considered by the 2.5P-1B analysis (in this precise sense of travelling at angle $\unicode[STIX]{x1D6FD}\neq 0$ to the background flow) are strictly only well defined in infinite domains with translational invariance in $x$ and $y$ . As another example of a ‘three-dimensional’ stability analysis, Drazin (Reference Drazin1974) considered the linear instability of an inviscid flow with two layers of different density, where the upper layer has a non-zero streamwise velocity which varies slowly in the spanwise direction. This flow is prone to a (truly) three-dimensional variant of the conventional Kelvin–Helmholtz instability, yet the background flow was still not confined or localised in any sense. Such flow instabilities, just like the flow instabilities ‘propagating at an angle’, are therefore unlikely to be observed experimentally. Furthermore, the modal structure (2.2) assumed for such perturbations inevitably imposes the same spanwise periodicity on the perturbation density and velocity components, thereby strongly restricting their possible spatial structure. As we discuss in more detail in this paper, perturbations with a more complex spatial structure are, in fact, consistent with our experimental observations.

Indeed, experimental observations of HWs are typically compared to HI ‘2P-1B’ predictions (often representing $U(z)$ by hyperbolic-tangent or error-function vertical profiles in unbounded domains). Analogously, three-dimensional numerical simulations typically consider periodic boundary conditions in both of the horizontal directions, avoiding potential issues with horizontal boundary conditions interacting with HWs. However, viscous effects and the presence of confining side walls in an experimental tank make any experimental mean flow inherently dependent on the spanwise ( $y$ ) and vertical ( $z$ ) coordinates and affect any waves which may develop by requiring the perturbations to decay to zero at these walls. Many of the geophysical flows mentioned in § 2.3.1 also exhibit two- or three-dimensional base flows and significant spatial confinement. It is therefore desirable to understand the effects of this inherent three-dimensionality and confinement and model them appropriately. Our experimental apparatus is designed to obtain data to achieve this.

3.1 The stratified inclined duct

We consider the flow in a duct connecting two reservoirs containing aqueous salt solutions of densities $\unicode[STIX]{x1D70C}_{0}\pm \unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$ (figure 1). The seal between the duct and the barrier separating the two reservoirs is made of flexible rubber, allowing the duct to tilt at a relatively small angle $\unicode[STIX]{x1D703}$ from the horizontal. The duct has length $L=1350$  mm and has a square cross-section of height and width $H=45$  mm (aspect ratio $L/H=30$ ). At the start of the experiment, the duct is opened, initiating an exchange flow (with zero net volume flux) between the reservoirs.

Figure 1.

Experimental flow geometry: an inclined duct (1) connects a reservoir of dense fluid (2) and a reservoir of light fluid (3). All fluid is seeded with particles for stereo particle image velocimetry (sPIV), and fluid initially in (3) is dyed for planar laser induced fluorescence (PLIF). The laser beam emitted from (4) is directed to the pair of oscillating mirrors in (5), before entering the light-sheet-producing optics mounted on the linear traverse (6). The scanning light sheet sweeps a measurement volume (7) whose successive planes are imaged for sPIV by cameras (8) and (9) and for PLIF by camera (10). The measurement volume inset shows the coordinate system and notation used, as discussed in § 3.2.

This flow is a two-layer stratified shear flow, forced by gravity in two distinct ways. The first contribution is a hydrostatic pressure gradient in each layer, resulting from the end of each duct sitting in reservoirs containing fluids of different densities, and is sufficient to drive an exchange flow at $\unicode[STIX]{x1D703}=0^{\circ }$ . The second contribution depends on the tilt angle: the resulting non-zero projection of gravity along the duct provides additional buoyancy forces. A positive angle $\unicode[STIX]{x1D703}>0^{\circ }$ (defined here by the duct being raised in the denser reservoir) reinforces the exchange flow by accelerating the light layer up (to the left) and the heavy layer down (to the right). Each individual reservoir measures $1\times 0.2\times 0.5$  m and holds approximately 100 l, allowing us to maintain a statistically steady exchange flow for several minutes. This steady flow continues until the controls at the ends of the duct (see below) are flooded by the accumulation of fluid of a different density coming from the other reservoir.

As outlined in § 1, Meyer & Linden (Reference Meyer and Linden2014) described and mapped four qualitatively different flow states obtained by varying the main parameters $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}$ and $\unicode[STIX]{x1D703}$ , using shadowgraph observations and a larger square duct with length $L=3000$  mm and height $H=100$  mm. At low $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}$ and $\unicode[STIX]{x1D703}$ , the flow is laminar with an undistorted interface. At larger values of either or both parameters, persistent coherent HWs distort the interface. Higher values of the forcing give rise to spatio-temporally intermittent turbulence at the interface, and eventually, fully developed turbulence with a thick mixing layer. Their measurements of layer-averaged velocities through mass flux balances support the hypothesis that the flow is hydraulically controlled at both ends of the duct for $\unicode[STIX]{x1D703}\gtrsim 1^{\circ }$ . The flow is subcritical (with respect to long, small-amplitude, interfacial disturbances) inside the duct, and undergoes two transitions at each end: a smooth one to a supercritical state shortly after exiting the duct, and an abrupt one further downstream (in the form of a jump) from supercritical to subcritical, to match the reservoir conditions. These hydraulic controls ensure that the exchange flow taking place through the duct is maximal, i.e. that its layer-averaged velocities obey the Froude number condition (3.1) (see Armi Reference Armi1986).

The practical consequence of this hydraulic control is that positive angles cannot increase the exchange flow rate, or mean flow velocities, beyond a certain value. The additional kinetic energy input provided by gravity through a positive tilt angle must be dissipated by the flow within the duct, either by increased wave activity or turbulence. We are therefore able to control the dissipation rate by the tilt angle, and to maintain target dissipation rates in a statistically steady sense for extended periods of time. The dissipation achieved is expressed non-dimensionally as a buoyancy Reynolds number in § 3.2, and depends on two free experimental parameters: the Reynolds number (when appropriately defined as in (3.4) and the tilt angle.

The inclined duct experiment is thus ideally suited to study HWs for the following reasons.

1. (i) It naturally sets up a flow with large shear layer thickness/density interface thickness ratio. The low diffusivity of salt compared to momentum (or high Schmidt number $Sc\approx 700$ , as defined in § 1), and the continuous supply of unmixed fluid from both reservoirs in the open duct mean that the interface thickness ratio resulting from interfacial diffusion of vorticity compared to interfacial diffusion of density is of the order of $\sqrt{Sc}=O(10)$ .

2. (ii) It naturally sets up a flow with large bulk Richardson number. The hydraulic control sets an upper bound for the speeds achieved inside the duct and, as we shall see in § 3.2 and from (3.5), effectively sets the overall or bulk Richardson number to an $O(0.1{-}1)$ constant, regardless of the forcing parameters.

3. (iii) Both symmetric and asymmetric HWs can be observed. Experiments at low angles $\unicode[STIX]{x1D703}\lesssim 2^{\circ }$ typically exhibit flows where the density interface is close to the centre of the shear layer and so is prone to symmetric HWs (i.e. two waves propagating in opposite directions with similar amplitude). In contrast, experiments at higher angles $\unicode[STIX]{x1D703}\gtrsim 2^{\circ }$ exhibit an offset between the mid-point of the shear layer and the density interface sufficient to give rise to asymmetric (i.e. one-sided) HWs, where the instability with smaller growth rate is apparently being sufficiently suppressed such that it does not reach an observable ‘wave’ with finite amplitude. The mechanism leading to the offset of the interface is thought to be related to the complex pressure adjustment necessary to match the different hydrostatic pressures in each reservoirs when the duct is tilted. We will use this empirical knowledge to our advantage in this paper, deliberately selecting parameters for which asymmetric HWs are observed.

3.2 Notation and non-dimensional parameters

The notation we use in the paper is shown in the measurement volume inset in figure 1. The streamwise $x$ axis is aligned along the duct and the spanwise $y$ axis across the duct, making the $z$ axis tilted at an angle $\unicode[STIX]{x1D703}$ from the vertical (resulting in a non-zero streamwise projection of gravity). All coordinates are centred in the middle of the duct, such that $-L/2\leqslant x\leqslant L/2$ and $-H/2\leqslant y,z\leqslant H/2$ . The velocity vector field has components $\unicode[STIX]{x1D66A}(x,y,z,t)=(u,v,w)$ along $x,y,z$ , respectively, and we denote the density field by $\unicode[STIX]{x1D70C}(x,y,z,t)$ .

Only a limited number of parameters are believed to play important roles in this experiment. The geometrical parameters are $L$ , $H$ , $\unicode[STIX]{x1D703}$ and the dynamical parameters are the acceleration due to gravity $g$ , the non-dimensional density difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}$ , the kinematic viscosity of water $\unicode[STIX]{x1D708}$ and the molecular diffusivity of salt $\unicode[STIX]{x1D705}_{s}$ . The relative density difference ( ${<}0.1\,\%$ in the experiment reported here) is sufficiently small so that the Boussinesq approximation requiring that $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}\ll 1$ is valid, and density differences only play a dynamically relevant role through the reduced gravity $g^{\prime }\equiv g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}$ . Since the six free parameters $L,H,\unicode[STIX]{x1D703},g^{\prime },\unicode[STIX]{x1D708},\unicode[STIX]{x1D705}_{s}$ have two independent dimensions (length and time), it is possible to construct four independent non-dimensional parameters, as we do below.

It is important to note at this point that the experiment does not have any imposed velocity scale. However, to a good approximation, the hydraulic controls at the ends of the duct require that the composite Froude number

(3.1) $$\begin{eqnarray}\displaystyle G^{2}\equiv F_{1}^{2}+F_{2}^{2}=1, & & \displaystyle\end{eqnarray}$$

where $F_{i}^{2}\equiv \langle u_{i}^{2}\rangle _{y,z}/(g^{\prime }h_{i})$ is the Froude number of layer $i$ and $\langle \cdot \rangle _{y,z}$ denotes spanwise and vertical averaging over the depth $h_{i}$ of each layer. The symmetry of this Boussinesq exchange flow implies that both layers have equal depth and Froude number in the central section of the duct

(3.2) $$\begin{eqnarray}\displaystyle F^{2}\equiv F_{1}^{2}(x=0)=F_{2}^{2}(x=0), & & \displaystyle\end{eqnarray}$$

such that for each layer

(3.3a,b ) $$\begin{eqnarray}\displaystyle F=\frac{\langle u\rangle _{y,z}}{\sqrt{g^{\prime }(H/2)}}=\pm \frac{1}{\sqrt{2}},\quad \text{or}\quad \langle u\rangle _{y,z}=\pm \frac{\sqrt{g^{\prime }H}}{2}. & & \displaystyle\end{eqnarray}$$

Due to the moderate values of $Re$ , the mean velocity profile is affected by viscosity and we find that the peak velocities in each layer are approximately twice the layer-averaged values, i.e. $\max |u|\approx \sqrt{g^{\prime }H}$ . We choose to non-dimensionalise velocities by this characteristic peak value, or half the total (peak-to-peak) velocity jump $\unicode[STIX]{x0394}U/2\equiv \sqrt{g^{\prime }H}$ . We define $(\tilde{u} ,\tilde{v},\tilde{w})\equiv (u,v,w)/(\unicode[STIX]{x0394}U/2)$ such that $-1\lesssim \tilde{u} \lesssim 1$ (note that in general $|\tilde{v}|,|\tilde{w}|\ll |\tilde{u} |$ ; see § 6.1.2). For consistency, we choose $H/2$ as the length scale, such that $-1\leqslant {\tilde{y}},\tilde{z}\leqslant 1$ , and $-L/H\leqslant \tilde{x}\leqslant L/H$ . Consequently, the natural time scale for non-dimensionalisation is the advective time $H/\unicode[STIX]{x0394}U=H/(2\sqrt{g^{\prime }H})$ . The resulting non-dimensional time units will be referred to as advective time units (ATU). The relevant Boussinesq density field is $\tilde{\unicode[STIX]{x1D70C}}\equiv (\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0})/(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2)$ , where we use $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$ to non-dimensionalise density such that $-1\leqslant \tilde{\unicode[STIX]{x1D70C}}\leqslant 1$ .

Using the previously defined scales, it is natural to construct the Reynolds number

(3.4) $$\begin{eqnarray}\displaystyle Re\equiv \frac{\displaystyle \frac{\unicode[STIX]{x0394}U}{2}\frac{H}{2}}{\unicode[STIX]{x1D708}}=\frac{\sqrt{g^{\prime }H}H}{2\unicode[STIX]{x1D708}}, & & \displaystyle\end{eqnarray}$$

which, for a given duct geometry and fluid, is a function of the driving density difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ alone.

The Froude condition (3.1) artificially adds another dimensional parameter, the velocity scale, to our previous set of six parameters. It will prove useful in the rest of the paper (e.g. for the governing equations (5.1)) to define the overall Richardson number of the flow, expressed as the non-dimensional product of the density, length and inverse square velocity scales

(3.5) $$\begin{eqnarray}\displaystyle Ri\equiv \frac{\displaystyle \frac{g}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\frac{H}{2}}{\left(\displaystyle \frac{\unicode[STIX]{x0394}U}{2}\right)^{2}}=\frac{1}{4}, & & \displaystyle\end{eqnarray}$$

by definition of $\unicode[STIX]{x0394}U$ . Note that the value of $Ri$ is set and not a free parameter; it is simply an equivalent formulation of the Froude condition (3.1): $Ri=1/(8F^{2})=1/(4G^{2})$ and $G^{2}=1$ . We choose the Schmidt number $Sc$ (defined in § 1), $L/H$ and $\unicode[STIX]{x1D703}$ as the last parameters. In summary, we have a total of four free independent non-dimensional parameters: $Sc$ , $L/H$ , $\unicode[STIX]{x1D703}$ , $Re$ and one imposed parameter $Ri$ . For the apparatus considered, we have $Sc=700$ , $L/H=30$ , $Ri=1/4$ , and we have the freedom to vary $\unicode[STIX]{x1D703}$ and $Re$ (through varying $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ ).

A characteristic buoyancy Reynolds number $Re_{b}$ (see § 1 and (2.1)) can be estimated for $\unicode[STIX]{x1D703}>0^{\circ }$ , assuming complete dissipation of the instantaneous power gained by the fluid moving at a layer-averaged velocity $\unicode[STIX]{x0394}U/4$ in the streamwise gravity field $g^{\prime }\sin \unicode[STIX]{x1D703}$ , i.e. the characteristic dissipation scale is $\unicode[STIX]{x1D716}_{0}\approx g^{\prime }(\unicode[STIX]{x0394}U/4)\sin \unicode[STIX]{x1D703}$ . Using the characteristic buoyancy frequency $N_{0}^{2}=g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/(2\unicode[STIX]{x1D70C}_{0})/(H/2)=g^{\prime }/H$ , we find that

(3.6) $$\begin{eqnarray}\displaystyle Re_{b}\approx \frac{g^{\prime }\displaystyle \frac{\unicode[STIX]{x0394}U}{4}\sin \unicode[STIX]{x1D703}}{\unicode[STIX]{x1D708}\displaystyle \frac{g^{\prime }}{H}}=\frac{\unicode[STIX]{x0394}UH}{4\unicode[STIX]{x1D708}}\sin \unicode[STIX]{x1D703}=Re\sin \unicode[STIX]{x1D703}. & & \displaystyle\end{eqnarray}$$

Note that this is a volume-averaged estimate that includes the dissipation due to viscous stresses on the duct walls, and that the actual dissipation can be highly heterogeneous in space. Therefore, it is necessary to be cautious when comparing the specific numerical values of $Re_{b}$ reported here with other studies where other definitions have been used. For positive angles $\unicode[STIX]{x1D703}>0^{\circ }$ , the three main qualitative bifurcations in flow regimes have been found to scale with $Re\sin \unicode[STIX]{x1D703}$ over a wide range of $Re$ and $\unicode[STIX]{x1D703}$ . Based on hundreds of shadowgraph observations, HWs are observed for $20\lesssim Re_{b}\lesssim 50$ , intermittent turbulence is observed for $50\lesssim Re_{b}\lesssim 100$ and steady turbulence for $Re_{b}\gtrsim 100$ (Lefauve Reference Lefauve2018).

Henceforth, we drop the tildes and, unless explicitly stated otherwise, use non-dimensional variables throughout.

3.3 Three-dimensional volumetric sPIV/PLIF measurements

We obtained time-resolved, three-dimensional, volumetric measurements of the three components of velocity using a novel method that combines successive two-dimensional planar stereo particle image velocimetry (sPIV) measurements to construct a volume. The flow was illuminated by a thin, pulsed, vertical laser sheet, which was rapidly scanned back and forth in the spanwise $y$ direction (figure 1), sweeping the volume of interest. The fundamental technical challenge with this approach is to obtain pairs of images, separated by a small time interval appropriate for PIV, at the same spanwise locations without relying on the overlap of successive laser sheets achieved either from the continuous traverse of an excessively thick laser sheet (which would compromise the spanwise resolution) or from a very slow scanning speed (which would compromise the near-instantaneous character of the measurements). The novelty of the system employed here is the conversion of the continuous motion of a traverse into a discontinuous motion of the light sheet. The laser beam was directed onto a pair of oscillating mirrors controlled by galvanometers (for accurate and fast positioning) that deflect its position before entering the optics mounted on a fast-moving linear traverse. This allowed us to control the successive positions of the pulsed laser sheet independently of the traverse, and thus obtain pairs of images at the same spanwise location while continuously scanning the thin laser sheet across the measurement volume. A detailed description of the scanning system and other aspects of the sPIV/planar laser induced fluorescence (PLIF) measurements can be found in Partridge, Lefauve & Dalziel (Reference Partridge, Lefauve and Dalziel2018).

The disadvantage of this approach is that the measurements are not instantaneous as successive planes have a time lag. However, this time lag can be made very small by increasing the scanning speed while maintaining the co-location of laser sheet pairs. As discussed in more detail in Partridge et al. (Reference Partridge, Lefauve and Dalziel2018), this method has three main advantages over existing methods for obtaining three-dimensional experimental diagnostics of a stratified flow across a volume: it is less sensitive to unavoidable residual refractive index variations; it yields a higher spatial resolution; and crucially, it allows simultaneous measurements of the density field.

Simultaneous measurements of the density field were achieved by combining successive PLIF measurements in the same fashion. PLIF requires an additional camera imaging the concentration of a fluorescent dye marking one of the fluids (in our case, the less dense solution, initially in region (3) of figure 1). The dye emits light at a slightly longer wavelength than the wavelength of the absorbed laser light. Narrow bandpass filters in front of the PIV cameras and a suitable long-pass filter for the PLIF camera allow separation of the two signals.

3.4 Experimental details

4.1 Instantaneous snapshots

Figure 2 shows four snapshots of the flow (in non-dimensional units), in the vertical mid-plane of the duct $y=0$ . The spanwise component of vorticity $\unicode[STIX]{x1D714}_{y}=\unicode[STIX]{x2202}_{z}u-\unicode[STIX]{x2202}_{x}w$ (a,c,e,g) is shown together with the density field $\unicode[STIX]{x1D70C}$ (b,d,f,h) at four different times $t=0,250,300,335$ (our choice of times is motivated by the analysis in § 4.3). A movie showing the full dynamics is available as supplementary movie 1 (https://doi.org/10.1017/jfm.2018.324). Both fields exhibit a quasi-periodic wave pattern corrugating the interface by upward-pointing cusps that propagate leftwards (with negative phase speed), as will be shown in § 4.2. The vorticity field shows concentration of negative vorticity that spans the interface downward to the right in tilde-shaped bands. The wavelength of these patterns shows some temporal variation, but both fields are coupled in the sense that their wavelength changes together. It is this wave that we refer to as a confined Holmboe wave (CHW).

A small amount of mixed fluid is present above the sharp interface, consistent with thin wisp-like ejections characteristic of asymmetric HWs, as discussed by Carpenter et al. (Reference Carpenter, Lawrence and Smyth2007) (see their figure 19). No wisp ejection was detected within the measurement volume; mixed fluid was either advected from other regions inside the duct where ejection took place and/or was a residual from the start of the experiment that had not been flushed out due to small velocities at the density interface.

Figure 2.

Snapshots of spanwise vorticity $\unicode[STIX]{x1D714}_{y}$ (a,c,e,g) and density $\unicode[STIX]{x1D70C}$ (b,d,f,h) in the mid-plane $y=0$ at $t=0$ (a,b), $t=250$ (c,d), $t=300$ (e,f) and $t=335$ (g,h). Axes are to scale. Throughout the paper, all data are shown in non-dimensional units.

To gain insight into the three-dimensional structure of the vorticity field, figure 3 shows, for the same times, the isosurface $\unicode[STIX]{x1D714}_{y}=-\unicode[STIX]{x1D714}_{max}/2=-2.5$ , where $\unicode[STIX]{x1D714}_{max}=5.0$ is the maximum value of $|\unicode[STIX]{x1D714}_{y}|$ found in the domain. This value represents a good intermediate value that shows the structure of interest. (Values $\unicode[STIX]{x1D714}_{y}\gtrsim -2$ mostly show the $x$ -independent mean shear while values $\unicode[STIX]{x1D714}_{y}\lesssim -3$ show only a small part of the wave structure.) To make observation easier, the $z$ axis is stretched by a factor of three, so it is important to remember that the actual structure is closer to horizontal than displayed.

The vorticity structure is confined both in the spanwise ( $y$ ) and vertical ( $z$ ) directions and takes the form of inclined, distorted, prolate spheroids with a wide (in $y$ ) ‘body’ and a progressively (in the negative $x$ direction) narrower ‘head’ (as indicated by the arrows in figure 3 a). A chevron shape at the posterior side of the body is partially visible in figure 3(c,d). The structures overlap one another and are connected by their bodies at the edges (the extent of the connection between structures naturally depends on the chosen isosurface). It is unknown whether the observed regular annular banding of the isosurface is a real, coherent signal or the result of the isosurface rendering of noise in the data. A movie showing the full dynamics is available as supplementary movie 2.

Figure 3.

Isosurfaces of spanwise vorticity $\unicode[STIX]{x1D714}_{y}=-2.5$ at the same times as in figure 2. The $z$ axis is stretched by a factor of 3 and only $z\in [-0.5,0.3]$ is shown.

4.2 Spatio-temporal behaviour

Figure 4 shows the spatio-temporal behaviour of the CHW. Slices of the vorticity field (figure 4 a) and density field (figure 4 b) taken at the mean positions of the respective interfaces ( $z=-0.11$ for the velocity interface, $z=-0.22$ for the density interface) are stacked in time to form a spatio-temporal (Hovmöller) diagram.

The wave propagates leftwards in both fields, i.e. with phase speed $c<0$ . The coupling of the fields is confirmed, as both travel at the same instantaneous phase speed and undergo the same gradual change in wavelength. The range of wavelengths observed from this figure is $\unicode[STIX]{x1D706}\in [3.5,6.0]$ , corresponding to a range in wavenumber $k\in [1.05,1.80]$ . The range of phase speeds, as determined by the slope of the characteristics, is $c\in [-0.18,0.05]$ .

There is evidence of a preferred wavelength, as splitting events occur when the wavelength becomes ‘too large’, presumably due to wave stretching (e.g. at $t\approx 180$ and $t\approx 270$ ), consistent with some development in $x$ of the mean flow (Tedford et al. Reference Tedford, Pieters and Lawrence2009; Carpenter et al. Reference Carpenter, Tedford, Rahmani and Lawrence2010). This development of the mean flow along $x$ is an inevitable consequence of the convective acceleration of the flow as it travels along the duct. In experiments for which gravity forcing dominates over the longitudinal pressure gradient ( $\unicode[STIX]{x1D703}\gtrsim \tan ^{-1}(H/L)\approx 2^{\circ }$ , as in the present case) this effect is small and the flow is nearly parallel. Splitting events require a difference in phase speed between the parent and daughter structures (stretching–splitting) and we indeed observe for a short time a daughter wave propagating with slightly positive phase speed while the parent continues at unchanged, negative phase speed (at $x\approx -10$ , $t\approx 180$ and $x\approx -13$ , $t\approx 270$ in figure 4).

Figure 4.

(a,b) Spatio-temporal diagrams taken in the vertical mid-plane and horizontal interfacial plane: (a) $\unicode[STIX]{x1D714}_{y}(y=0,z=-0.11)$ and (b) $\unicode[STIX]{x1D70C}(y=0,z=-0.22)$ . Dashed sloping lines are characteristics with phase speed $c=-0.078$ and wavelength $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}/1.46=4.3$ . Black solid horizontal lines indicate the start ( $t=300$ ) and end ( $t=335$ ) of the phase-averaging window, represented by thick sloping lines. (c,d) Two views of the isosurface of phase-averaged spanwise vorticity $\langle \unicode[STIX]{x1D714}_{y}\rangle =-2.5$ . Two wavelengths are shown by repeating the structure along $x$ . Note the slight lack of periodicity evidenced by the discontinuity at $x=4.30$ . The $z$ axis is stretched by a factor of 3 and only $z\in [-0.5,0.3]$ is shown.

4.3 Phase-averaged properties

To characterise the ‘typical’ structure of the CHW and reduce its spatio-temporal complexity, we phase-average the flow variables over a single wavelength by following the wave along a characteristic, an operation denoted by $\langle \cdot \rangle$ . To make this calculation, we determine the phase speed $c$ and wavelength $\unicode[STIX]{x1D706}$ by considering the interval $t\in [300,335]$ towards the end of the experiment. This period, indicated by solid horizontal lines in figure 4(a,b), was selected as it is during this period that the wave has reached a very nearly constant negative phase speed and a steady wavelength, making the analysis easier and more meaningful. We determine the phase speed $c=-0.078$ and wavelength $\unicode[STIX]{x1D706}=4.30$ (wavenumber $k=1.46$ ) by fitting in time the location of the vorticity minimum contained within the wavelength denoted by two thick solid sloping lines in figure 4(a,b). The dashed sloping lines having slope $c$ and spacing $\unicode[STIX]{x1D706}$ , drawn half a wavelength from the solid lines, extend over the entire spatio-temporal plot and demonstrate that the phase-averaged properties are representative of the flow at earlier times.

The three-dimensional structure of the resulting phase-averaged vorticity in the region $x\in [0,4.30]$ is shown in figure 4(c,d), using the same isosurface level $\langle \unicode[STIX]{x1D714}_{y}\rangle =-2.5$ as in figure 3. For better visualisation, this region has been replicated to $x\in [4.30,8.60]$ ) to show two wavelengths for the structure. The two panels show views from different angles. Naturally, we recover the features identified in figure 3, but are able to discern further details. The structure is remarkably symmetrical about $y=0$ and very nearly periodic (see the small discontinuity at $x=\unicode[STIX]{x1D706}=4.30$ where the replicated structures join), despite the slow stretching in $x$ discussed above. It is also triply connected, and its tilt with respect to $x$ (the duct axis) is $10^{\circ }$ (maximum angle at the inflection point).

We discuss further details of the phase-averaged structure of the CHW (including the other flow variables $\unicode[STIX]{x1D70C},u,v,w$ ) in § 6 when comparing it with the linear confined Holmboe instability (CHI) predicted by the stability analysis introduced in the next section § 5.

5.1 Three-dimensional perturbations, two-dimensional base flow (3P-2B)

We assume a steady parallel two-dimensional base flow $\boldsymbol{U}\equiv [U(y,z),0,0]$ and one-dimensional density distribution $R(z)$ (note that (5.1) do not support steady states with this base flow and $R(y,z)$ ). We superimpose three-dimensional infinitesimal perturbations, so that the full velocity field is

(5.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D66A}\equiv \left[\begin{array}{@{}c@{}}U(y,z)\\ 0\\ 0\end{array}\right]+\unicode[STIX]{x1D700}\left[\begin{array}{@{}c@{}}u^{\prime }(x,y,z,t)\\ v^{\prime }(x,y,z,t)\\ w^{\prime }(x,y,z,t)\end{array}\right], & & \displaystyle\end{eqnarray}$$

and similarly $\unicode[STIX]{x1D70C}\equiv R(z)+\unicode[STIX]{x1D700}\unicode[STIX]{x1D70C}^{\prime }(x,y,z,t),p\equiv P(z)+\unicode[STIX]{x1D700}p^{\prime }(x,y,z,t)$ , ( $|\unicode[STIX]{x1D700}|\ll 1$ ). We assume that each perturbation variable $\unicode[STIX]{x1D713}^{\prime }$ has periodic wavelike behaviour in $x$ and $t$ ,

(5.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}^{\prime }(x,y,z,t)\equiv \hat{\unicode[STIX]{x1D713}}(y,z)\exp (\text{i}kx+\unicode[STIX]{x1D70E}t), & & \displaystyle\end{eqnarray}$$

where the real part is implied. Note the fundamental distinction between this 3P-2B form and the 2.5P-1B form in (2.2). We tackle the temporal stability problem by specifying the perturbation wavenumber $k>0\in \mathbb{R}$ and solve the following eigenvalue problem for the (temporal) growth rate $\unicode[STIX]{x1D70E}\in \mathbb{C}$ and eigenfunctions $\hat{\unicode[STIX]{x1D713}}(y,z)\in \mathbb{C}$ :

(5.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D6FB}^{2} & & \\ & \unicode[STIX]{x1D6FB}^{2} & \\ & & {\mathcal{I}}\end{array}\right]\left[\begin{array}{@{}c@{}}\hat{v}\\ {\hat{w}}\\ \hat{\unicode[STIX]{x1D70C}}\end{array}\right]=\left[\begin{array}{@{}ccc@{}}{\mathcal{L}}_{v} & {\mathcal{L}}_{vw} & {\mathcal{L}}_{v\unicode[STIX]{x1D70C}}\\ {\mathcal{L}}_{wv} & {\mathcal{L}}_{w} & {\mathcal{L}}_{w\unicode[STIX]{x1D70C}}\\ & {\mathcal{L}}_{\unicode[STIX]{x1D70C}w} & {\mathcal{L}}_{\unicode[STIX]{x1D70C}}\end{array}\right]\left[\begin{array}{@{}c@{}}\hat{v}\\ {\hat{w}}\\ \hat{\unicode[STIX]{x1D70C}}\end{array}\right], & & \displaystyle\end{eqnarray}$$

where

(5.5a ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{v} & \equiv & \displaystyle -\text{i}kU\unicode[STIX]{x1D6FB}^{2}+\text{i}k(\unicode[STIX]{x2202}_{yy}U-\unicode[STIX]{x2202}_{zz}U)-2\text{i}k\unicode[STIX]{x2202}_{z}U\unicode[STIX]{x2202}_{z}+Re^{-1}\unicode[STIX]{x1D6FB}^{4},\end{eqnarray}$$

(5.5b ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{w} & \equiv & \displaystyle -\text{i}kU\unicode[STIX]{x1D6FB}^{2}+\text{i}k(\unicode[STIX]{x2202}_{zz}U-\unicode[STIX]{x2202}_{yy}U)-2\text{i}k\unicode[STIX]{x2202}_{y}U\unicode[STIX]{x2202}_{y}+Re^{-1}\unicode[STIX]{x1D6FB}^{4},\end{eqnarray}$$

(5.5c ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{vw} & \equiv & \displaystyle 2\text{i}k(\unicode[STIX]{x2202}_{yz}U+\unicode[STIX]{x2202}_{z}U\unicode[STIX]{x2202}_{y}),\end{eqnarray}$$

(5.5d ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{wv} & \equiv & \displaystyle 2\text{i}k(\unicode[STIX]{x2202}_{yz}U+\unicode[STIX]{x2202}_{y}U\unicode[STIX]{x2202}_{z}),\end{eqnarray}$$

(5.5e ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{\unicode[STIX]{x1D70C}} & \equiv & \displaystyle -\text{i}kU+(Re\,Sc)^{-1}\unicode[STIX]{x1D6FB}^{2},\end{eqnarray}$$

(5.5f ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{v\unicode[STIX]{x1D70C}} & \equiv & \displaystyle Ri(\cos \unicode[STIX]{x1D703}\,\unicode[STIX]{x2202}_{yz}-\text{i}k\sin \unicode[STIX]{x1D703}\,\unicode[STIX]{x2202}_{y}),\end{eqnarray}$$

(5.5g ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{w\unicode[STIX]{x1D70C}} & \equiv & \displaystyle Ri\{\cos \unicode[STIX]{x1D703}\,(\unicode[STIX]{x2202}_{zz}-\unicode[STIX]{x1D6FB}^{2})-\text{i}k\sin \unicode[STIX]{x1D703}\,\unicode[STIX]{x2202}_{z}\},\end{eqnarray}$$

(5.5h ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{\unicode[STIX]{x1D70C}w} & \equiv & \displaystyle -\unicode[STIX]{x2202}_{z}R.\end{eqnarray}$$

${\mathcal{I}}$

$\unicode[STIX]{x1D6FB}^{2}\equiv -k^{2}+\unicode[STIX]{x2202}_{yy}+\unicode[STIX]{x2202}_{zz}$

$\unicode[STIX]{x1D6FB}^{4}\equiv k^{4}+\unicode[STIX]{x2202}_{yyyy}+\unicode[STIX]{x2202}_{zzzz}+2\unicode[STIX]{x2202}_{yy}\unicode[STIX]{x2202}_{zz}-2k^{2}(\unicode[STIX]{x2202}_{yy}+\unicode[STIX]{x2202}_{zz})$

(5.6) $$\begin{eqnarray}\displaystyle \hat{u} =-\frac{\text{i}}{k}(\unicode[STIX]{x2202}_{y}\hat{v}+\unicode[STIX]{x2202}_{z}{\hat{w}}). & & \displaystyle\end{eqnarray}$$

The derivation of (5.4)–(5.6) is given in appendix A.

We represent the rigid, impermeable walls of the duct by imposing no-slip boundary conditions for velocity perturbations: $\hat{v}={\hat{w}}=0$ along $y,z=\pm 1$ . To impose $\hat{u} =0$ along $y,z=\pm 1$ , (5.6) requires that we further impose $\unicode[STIX]{x2202}_{y}\hat{v}=0$ along $y=\pm 1$ and $\unicode[STIX]{x2202}_{z}{\hat{w}}=0$ along $z=\pm 1$ . Physically, no salt diffuses through the walls of the duct, hence the boundary conditions for the density perturbation are $\unicode[STIX]{x2202}_{y}\hat{\unicode[STIX]{x1D70C}}=0$ along $y=\pm 1$ and $\unicode[STIX]{x2202}_{z}\hat{\unicode[STIX]{x1D70C}}=0$ along $z=\pm 1$ .

We solve the eigenvalue problem (5.4) numerically, using second-order finite-difference discretisation with uniform spacing in both $y$ and $z$ . A resolution of $120\times 120$ proved to be sufficient for convergence of the results, and required solving a $(3\times 118^{2})\times (3\times 118^{2})$ eigenvalue problem. The results are presented in § 6.1.

5.2 Two-dimensional perturbations, one-dimensional base flow (2P-1B)

To contextualise the results of the formulation introduced above, we will also discuss the results of the more commonly used 2P-1B analysis. This represents the special case for which the spanwise velocity perturbation $v$ and all variations in the spanwise direction $y$ are ignored, so that the full velocity field is

(5.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D66A}_{\dagger }\equiv \left[\begin{array}{@{}c@{}}U(z)\\ 0\\ 0\end{array}\right]+\unicode[STIX]{x1D700}\left[\begin{array}{@{}c@{}}u_{\dagger }^{\prime }(x,z,t)\\ 0\\ w_{\dagger }^{\prime }(x,z,t)\end{array}\right], & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D70C}_{\dagger }\equiv R(z)+\unicode[STIX]{x1D700}\unicode[STIX]{x1D70C}_{\dagger }^{\prime }(x,z,t),p_{\dagger }\equiv P(z)+\unicode[STIX]{x1D700}p_{\dagger }^{\prime }(x,z,t),$ with all perturbations

(5.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{\dagger }^{\prime }(x,z,t)\equiv \hat{\unicode[STIX]{x1D713}}_{\dagger }(z)\exp (\text{i}kx+\unicode[STIX]{x1D70E}t). & & \displaystyle\end{eqnarray}$$

Note that this 2P-1B form corresponds to the special case $\unicode[STIX]{x1D6FD}=0$ of the 2.5-1B form in (2.2). Here and below, the subscript/superscript $\dagger$ distinguishes the variables and operators of the 2P-1B analysis from those of the 3P-2B analysis.

The governing equations for this problem are obtained from (5.4)–(5.6) by removing $\hat{v}$ and setting $\unicode[STIX]{x2202}_{y}=0$ , i.e.

(5.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D6FB}_{\dagger }^{2} & \\ & {\mathcal{I}}\end{array}\right]\left[\begin{array}{@{}c@{}}{\hat{w}}_{\dagger }\\ \hat{\unicode[STIX]{x1D70C}}_{\dagger }\end{array}\right]=\left[\begin{array}{@{}cc@{}}{\mathcal{L}}_{w}^{\dagger } & {\mathcal{L}}_{w\unicode[STIX]{x1D70C}}^{\dagger }\\ {\mathcal{L}}_{\unicode[STIX]{x1D70C}w}^{\dagger } & {\mathcal{L}}_{\unicode[STIX]{x1D70C}}^{\dagger }\end{array}\right]\left[\begin{array}{@{}c@{}}{\hat{w}}_{\dagger }\\ \hat{\unicode[STIX]{x1D70C}}_{\dagger }\end{array}\right], & & \displaystyle\end{eqnarray}$$

where

(5.10a ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{w}^{\dagger } & \equiv & \displaystyle -\text{i}kU\unicode[STIX]{x1D6FB}_{\dagger }^{2}+\text{i}k\unicode[STIX]{x2202}_{zz}U+Re^{-1}\unicode[STIX]{x1D6FB}_{\dagger }^{4},\end{eqnarray}$$

(5.10b ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{\unicode[STIX]{x1D70C}}^{\dagger } & \equiv & \displaystyle -\text{i}kU+(Re\,Sc)^{-1}\unicode[STIX]{x1D6FB}_{\dagger }^{2},\end{eqnarray}$$

(5.10c ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{w\unicode[STIX]{x1D70C}}^{\dagger } & \equiv & \displaystyle Ri\{\cos \unicode[STIX]{x1D703}(\unicode[STIX]{x2202}_{zz}-\unicode[STIX]{x1D6FB}_{\dagger }^{2})-\text{i}k\sin \unicode[STIX]{x1D703}\,\unicode[STIX]{x2202}_{z}\},\end{eqnarray}$$

(5.10d ) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{\unicode[STIX]{x1D70C}w}^{\dagger } & \equiv & \displaystyle -\unicode[STIX]{x2202}_{z}R,\end{eqnarray}$$

$\unicode[STIX]{x1D6FB}_{\dagger }^{2}\equiv -k^{2}+\unicode[STIX]{x2202}_{zz}$

$\unicode[STIX]{x1D6FB}_{\dagger }^{4}\equiv k^{4}+\unicode[STIX]{x2202}_{zzzz}-2k^{2}+\unicode[STIX]{x2202}_{zz}$

$\hat{u} _{\dagger }=-(i/k)\unicode[STIX]{x2202}_{z}{\hat{w}}_{\dagger }$

$z$

$z=\pm 1$

${\hat{w}}_{\dagger },\hat{u} _{\dagger }$

${\hat{w}}_{\dagger }=\unicode[STIX]{x2202}_{z}{\hat{w}}_{\dagger }=0$

$\hat{\unicode[STIX]{x1D70C}}_{\dagger }$

$\unicode[STIX]{x2202}_{z}\hat{\unicode[STIX]{x1D70C}}_{\dagger }=0$

5.3 Experimental base flow

In order to obtain linear stability predictions which are most relevant for comparison with the CHW observed in the experiment, we use averages of the experimental flow. To obtain a representative, canonical $U(y,z)$ for the 3P-2B analysis, the streamwise velocity $u$ was first averaged along $x$ and over 50 ATU before the start of the phase-averaging window, i.e. for $t\in [250,300]$ , and symmetrised about $y=0$ : $U(y,z)\equiv [\langle u\rangle _{x,t}(y,z)+\langle u\rangle _{x,t}(-y,z)]/2$ (figure 5 a). The profile used for the 2P-1B analysis is $U(z)\equiv U(y=0,z)$ (figure 5 b). Note that since the flow is close to steady, the duration of $50$  ATU for the phase-averaging window is arbitrary and the results were only weakly sensitive to this choice.

The measured density field was averaged over the same $x$ and $t$ windows, with an additional averaging over $y\in [-1,1]$ since only the one-dimensional base density profile $R(z)$ enters the stability problem. Being distorted by a finite-amplitude HW, the interface does not lie on $z=\text{const}$ . Since the linear stability properties of HIs are known to be sensitive to the thickness of the interface, it is important to ensure that $R(z)$ has a thickness representative of undistorted local profiles. We achieve this by conditionally averaging over $x$ and $y$ , i.e. shifting all profiles vertically in order to align them at the interface before averaging them (shown in dashed green in figure 5 c). Finally, our aim is to elucidate the origin of CHWs (i.e. to identify and understand the underlying physical mechanisms leading to their observed behaviour) and not to capture either the initial transient of the experiment or the nonlinear wisp ejection dynamics, which presumably produces the slightly mixed region in the region $-0.2\lesssim z\lesssim 0.1$ described in § 4.1. We therefore eliminate this mixed region by fitting a hyperbolic tangent function (shown in solid black) to the lower (unmixed) region ( $z\leqslant -0.2$ ) of the previously obtained conditional average. The best-fit base profile we use for the 3P-2B and 2P-1B analyses is $R(z)\equiv -\tanh \{(z-z_{0})/\unicode[STIX]{x1D6FF}\}$ with $z_{0}=-0.22$ and $\unicode[STIX]{x1D6FF}=0.047$ .

Figure 5.

Base flow for stability analysis using experimental flow averaged over $t\in [250,300]$ before the phase-averaging window. (a) Streamwise velocity $U(y,z)$ for the 3P-2B analysis. (b) $U(z)=U(y=0,z)$ for the 2P-1B analysis. (c) Density: in dashed green, $\langle \unicode[STIX]{x1D70C}\rangle _{x,y,t}$ conditioned at the interface, and in solid black, $R(z)$ the $\tanh$ fit with the mixed layer removed as used for the 3P-2B and 2P-1P analyses (see text for details).

6.1 3P-2B results

6.1.1 Dispersion relation

The dispersion relation $\unicode[STIX]{x1D70E}(k)$ , where $\unicode[STIX]{x1D70E}\equiv \unicode[STIX]{x1D70E}_{r}+\text{i}\unicode[STIX]{x1D70E}_{i}$ (growth rate $\unicode[STIX]{x1D70E}_{r}$ and phase speed $c\equiv -\unicode[STIX]{x1D70E}_{i}/k$ ), is shown for growing modes $\unicode[STIX]{x1D70E}_{r}>0$ only in figure 6 (thick blue curve). We observe a lower wavenumber band of positive growth rate signalling an unstable mode for the range of wavenumbers $k\in [0.8,2.2]$ , as well as another band with lower growth rate for $k\in [2.3,5.0]$ . The fastest-growing mode (with largest $\unicode[STIX]{x1D70E}_{r}$ ) is predicted to occur at $k=1.32$ with growth rate $\unicode[STIX]{x1D70E}_{r}=0.0824$ and negative phase speed $c=-0.213$ . The other mode of instability has maximum growth rate of $\unicode[STIX]{x1D70E}_{r}=0.0192$ and positive phase speed $c=0.55$ . As discussed in § 2, the existence of a faster and a slower-growing Holmboe modes with different phase speed magnitudes, growth rates and wavenumbers is typical of asymmetric profiles (where the relatively sharp density interface and mid-point of the shear layer are not coincident). In the following we will focus on the faster-growing, leftward-propagating band of instability, which we identify as a confined Holmboe instability (CHI).

Figure 6.

Dispersion relation: (a) positive growth rate $\unicode[STIX]{x1D70E}_{r}$ and (b) phase speed $c$ . 3P-2B analysis (thick blue curve) and 2P-1B analysis (green dotted curve). (The blue star and green square mark the respective theoretically predicted fastest-growing modes.) Experimental range observed and values used for the phase average are respectively indicated by grey shading and black dashed lines.

The ranges of wavenumbers and phase speeds observed in the experimental CHW (figure 4 a,b) are shown as grey shading: $k\in [1.05,1.80]$ , $c\in [-0.18,0.05]$ , and the phase-average values ( $k=1.46$ , $c=-0.078$ ) as black dashed lines. The wavenumber of the fastest-growing mode compares well with the experimental values, and with the phase-averaged structure ( $k=1.32$ versus $k=1.46$ ). The agreement in phase speed is qualitatively good overall, despite the fastest-growing mode travelling significantly faster (in magnitude) than the representative phase-average value ( $c=-0.213$ versus $c=-0.078$ ). There is overlap between the unstable branch and the intersection of the two perpendicular grey boxes in figure 6(b), i.e. the shortest observed CHWs having $k\in [1.5,1.8]$ and moving at the highest speeds $c\in [-0.18,-0.15]$ are consistent with the slower-growing modes of the CHI.

A variety of causes can be invoked for the discrepancy in phase speed.

1. (i) The base velocity (i.e. with $\unicode[STIX]{x2202}_{x}U(y,z)=0$ ) is only an approximation to the slowly developing, non-parallel experimental flow, and such streamwise variation naturally affects the wavelength and speed of HWs.

2. (ii) The base flow used in the stability calculation is a spatio-temporal average from the unsteady experimental flow that contains the wave. The most obvious effects of this wave on the density profile have been partially eliminated in determining $R(z)$ , but not in determining $U(y,z)$ . It is not a priori obvious that this intuitive choice of base flow should correctly predict the linear instability that grows to finite amplitude and hence is observed in the experimental flow.

3. (iii) The experimental flow is not the result of an initial value problem. The nonlinear wave state selected by the flow does not have to result from the saturation of a monochromatic instability (the fastest growing mode), and there is therefore no reason to think that the nonlinear phase speed should match closely the phase speed of the fastest-growing mode.

To summarise, the analysis of the dispersion relation reveals reasonably good agreement between the predicted CHI and the observed CHW, despite the challenge of inferring the nonlinear properties of our flow from a linear analysis. In the next section we focus on comparing the three-dimensional structure of the fastest-growing CHI with the CHW.

6.1.2 Eigenfunctions

We turn to the three-dimensional structure of the eigenfunction of the fastest-growing CHI. To provide meaningful comparison with experimental data, the vorticity eigenfunction is added to the base flow vorticity

(6.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{y}^{\ast }\equiv \unicode[STIX]{x2202}_{z}U+\unicode[STIX]{x1D6FC}\,\text{Re}\{\hat{\unicode[STIX]{x1D714}}_{y}(y,z)\exp (\text{i}kx+\unicode[STIX]{x1D719})\}, & & \displaystyle\end{eqnarray}$$

where the superscript $^{\ast }$ denotes this reconstruction and $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D719}$ are constants to be determined. The phase parameter $\unicode[STIX]{x1D719}$ was simply chosen to compare directly with the phase of figure 4. The scaling parameter $\unicode[STIX]{x1D6FC}$ contains the information relative to the finite amplitude of the wave observed in the experiment. It was determined by imposing the condition that the volume-averaged root mean square (r.m.s.) of the wave perturbation in the reconstructed vorticity $\unicode[STIX]{x1D714}_{y}^{\ast }-\unicode[STIX]{x2202}_{z}U$ is equal to that of the wave perturbation in the phase-averaged vorticity $\langle \unicode[STIX]{x1D714}_{y}\rangle -\unicode[STIX]{x2202}_{z}U$ , such that

(6.2) $$\begin{eqnarray}\displaystyle \langle (\unicode[STIX]{x1D714}_{y}^{\ast }-\unicode[STIX]{x2202}_{z}U)^{rms}\rangle _{x,y,z}=\frac{\unicode[STIX]{x1D6FC}}{2}\sqrt{\langle \hat{\unicode[STIX]{x1D714}}_{y}^{2}\rangle _{y,z}}=\langle (\langle \unicode[STIX]{x1D714}_{y}\rangle -\unicode[STIX]{x2202}_{z}U)^{rms}\rangle _{x,y,z}=0.31, & & \displaystyle\end{eqnarray}$$

where r.m.s. denotes the operator $\unicode[STIX]{x1D713}^{rms}(x,y,z)\equiv \sqrt{(\unicode[STIX]{x1D713}-\langle \unicode[STIX]{x1D713}\rangle _{x,y,z})^{2}}$ . The first equality is a simple consequence of the definition in (6.1), the second equality is our condition and the third is the experimentally measured value.

The isosurface $\unicode[STIX]{x1D714}_{y}^{\ast }=-2.5$ of the resulting reconstruction for two wavelengths is shown in figure 7. A direct comparison with the phase average $\langle \unicode[STIX]{x1D714}_{y}\rangle =-2.5$ of figure 4(c,d) reveals excellent agreement between the three-dimensional structure of the CHI and the CHW. The tilt of the structure with respect to the duct axis is $8^{\circ }$ (maximum angle at the inflection point), slightly smaller than the $10^{\circ }$ obtained from the experimental phase average. The general geometry, including the detail of the triple connectedness of the head-to-tail connection, is faithfully reproduced. Supplementary movie 3 offers a panoramic visualisation of the vorticity isosurface of the CHW and CHI, allowing for more detailed comparison.

Figure 7.

Fastest-growing mode of the CHI: isosurface $\unicode[STIX]{x1D714}_{y}^{\ast }=-2.5$ (see (6.1) for definition of $\unicode[STIX]{x1D714}_{y}^{\ast }$ ). We show a double wavelength under two different views (a) and (b) for direct comparison with figures 4(c) and (d), respectively. As in those figures, the $z$ axis is stretched by a factor of 3 and only $z\in [-0.5,0.3]$ is shown for clearer visualisation.

To reveal the structure of the wave further and to allow for more detailed comparison, we build, by analogy with (6.1), a more complete picture of the wave field by considering the following variables associated with the fastest-growing mode,

(6.3a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\ast } & \equiv & \displaystyle R(z)+\unicode[STIX]{x1D6FC}\,\text{Re}\{\hat{\unicode[STIX]{x1D70C}}(y,z)\exp (\text{i}kx+\unicode[STIX]{x1D719})\},\end{eqnarray}$$

(6.3b ) $$\begin{eqnarray}\displaystyle u^{\ast } & \equiv & \displaystyle U(y,z)+\unicode[STIX]{x1D6FC}\,\text{Re}\{\hat{u} (y,z)\exp (\text{i}kx+\unicode[STIX]{x1D719})\},\end{eqnarray}$$

(6.3c ) $$\begin{eqnarray}\displaystyle v^{\ast } & \equiv & \displaystyle \unicode[STIX]{x1D6FC}\,\text{Re}\{\hat{v}(y,z)\exp (\text{i}kx+\unicode[STIX]{x1D719})\},\end{eqnarray}$$

(6.3d ) $$\begin{eqnarray}\displaystyle w^{\ast } & \equiv & \displaystyle \unicode[STIX]{x1D6FC}\,\text{Re}\{{\hat{w}}(y,z)\exp (\text{i}kx+\unicode[STIX]{x1D719})\},\end{eqnarray}$$

using the previously determined value of $\unicode[STIX]{x1D6FC}$ . We plot these wave variables in planar cuts in figure 8 (experimentally observed CHW) and figure 9 (theoretically predicted fastest-growing CHI). These figures combine views in the vertical mid-plane $y=0$ , horizontal interfacial planes ( $z=-0.22$ and $z=-0.11$ for density and velocity, respectively) and one cross-sectional plane for $\unicode[STIX]{x1D714}_{y},\unicode[STIX]{x1D70C},u,v,w$ for a single wavelength.

There is good agreement between the observed and predicted structure in each of these variables. Even though only the magnitude of the perturbation vorticity field was matched with the experimental value, the other variables in figure 9 have amplitudes very close to those in figure 8 (colour bars have the same limits in both figures and black contours are the same). The spanwise and vertical velocities $v$ and $w$ are small relative to the streamwise velocity, with a maximum of approximately 6 % of the maximum of $u$ . The quantitative and qualitative agreement extends to a number of interesting features as follows.

1. (i) The sheared vertical pattern in $w$ , typical of Holmboe modes, is clearly identifiable in m of figures 8 and 9.

2. (ii) The horizontal chevron pattern identified earlier in the three-dimensional visualisation of the vorticity $\unicode[STIX]{x1D714}_{y}$ is again visible in b (highlighted by the contours). It is also found in the streamwise velocity $u$ (h), where it contributes the main signal in this interfacial plane where the mean velocity is zero.

3. (iii) The relative phase of the density $\unicode[STIX]{x1D70C}$ and vorticity $\unicode[STIX]{x1D714}_{y}$ reveals that the ‘head’ of the vorticity structure is shifted slightly left of the top of the upward-pointing density cusp, while its ‘body’ appears to rest on the sloping interface of the cusp for which $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D70C}<0$ (a,b,d,e). This gradient is responsible for the baroclinic production of negative vorticity (reinforcing the negative mean shear). The density crests and troughs have a distinct convex structure (their vertical displacement peaks in the middle of the duct $y=0$ and decays near the boundaries $y=\pm 1$ ), making the gradients $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D70C}$ larger in the middle of the duct, presumably responsible for confinement of the vorticity and the nonlinear maintenance of the CHW at finite amplitude.

4. (iv) The variables clearly have different symmetries about the mid-plane $y=0$ , as is apparent in e,h,k,n. The structure of the CHI correctly predicts that the density $\hat{\unicode[STIX]{x1D70C}}$ , streamwise velocity perturbation $\hat{u}$ and vertical velocity perturbation ${\hat{w}}$ are even functions about this mid-plane, while the spanwise velocity perturbation $\hat{v}$ is an odd function, implying that the CHI is of varicose rather than sinuous form. It is important to stress that we did not impose these symmetries; they arose naturally as solutions to the 3P-2B stability eigenvalue problem. Furthermore, and perhaps unsurprisingly, the 2.5P-1B description of a normal mode travelling at an angle as described in § 2.3.3 is not appropriate to describe this structure, as the spanwise periodicity of $\hat{v}$ is clearly different from the other perturbation quantities $\hat{u}$ , ${\hat{w}}$ and $\hat{\unicode[STIX]{x1D70C}}$ . The odd symmetry of $\hat{v}$ seems to be inherently related (for reasons which are not yet fully understood) to the spanwise boundary conditions $\hat{v}=\unicode[STIX]{x2202}_{y}\hat{v}=0$ imposed by no-slip and incompressibility at the duct walls at $y=\pm 1$ . Indeed, a very similar eigenstructure was observed when we solved the 3P-2B problem with these boundary conditions for the spanwise uniform base flow $U(z)$ . (For reasons of brevity we do not discuss this problem further here.)

5. (v) A consequence of the specific symmetry of the CHW/CHI is the horizontal, in-plane divergence and convergence of $v$ (k,l) around the upward density crest (e,f). Due to the relatively large stratification ( $Ri=1/4$ ), the horizontal flow around the density crest is energetically preferred to the vertical flow above it. The convex structure of the density crest allows the flow to move around it horizontally. An important consequence is that it gives rise to relatively large spanwise gradients $|\unicode[STIX]{x2202}_{y}v|$ , positive in the centre of the duct $y\in [-0.5,0.5]$ and negative near the boundaries $|y|\in [0.5,1]$ . These gradients have a vortex stretching effect on $\unicode[STIX]{x1D714}_{y}$ through the term $\unicode[STIX]{x1D714}_{y}\unicode[STIX]{x2202}_{y}v$ , producing negative vorticity in the centre (reinforcing the mean shear), and positive vorticity near the boundaries (weakening the mean shear). They could play an important role in explaining the spanwise confinement (largest values of $|\unicode[STIX]{x1D714}_{y}|$ in the centre) and maintenance of the CHW at finite amplitude.

Figure 8.

Experimental phase-averaged spanwise vorticity $\langle \unicode[STIX]{x1D714}_{y}\rangle$ (a–c), density $\langle \unicode[STIX]{x1D70C}\rangle$ (d–f), and velocity components $\langle u\rangle$ (g–i), $\langle v\rangle$ (j–l), $\langle w\rangle$ (m–o). The panels show the vertical mid-plane $y=0$ (a,d,g,j,m), horizontal interfacial plane ( $z=-0.11$ for velocity and $z=-0.22$ for density) (b,e,h,k,n) and a cross-sectional plane $x=0.72$ (c,f,i,l,o). The colour bars are the same for a given variable. Ten black contours equally spaced across the range of the colour bar are also shown, except for (b) and (h) where the ten contours are across $[-4,-1]$ and $[-0.25,0.25]$ respectively, to highlight features of interest. Dashed lines represent the planes in the other panels.

Figure 9.

Fastest-growing CHI of the 3P-2B linear stability analysis superimposed on the base flow (see (6.1), (6.2) and (6.3) for definitions of the variables). All panels plotted as in figure 8 for direct comparison.

6.2 2P-1B results

6.2.2 Eigenfunctions

By analogy with the approach followed in § 6.1.2, the one-dimensional vorticity eigenfunction of the fastest-growing mode was added to the one-dimensional base flow such that

(6.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{y\dagger }^{\ast }\equiv \unicode[STIX]{x2202}_{z}U+\unicode[STIX]{x1D6FC}_{\dagger }\,\text{Re}\{\hat{\unicode[STIX]{x1D714}}_{y\dagger }(z)\exp (\text{i}kx+\unicode[STIX]{x1D719}_{\dagger })\}, & & \displaystyle\end{eqnarray}$$

where the amplitude $\unicode[STIX]{x1D6FC}_{\dagger }$ was determined in the same fashion as for the 3P-2B analysis, i.e. by matching the r.m.s. of the vorticity eigenfunction to the experimental perturbation in the mid-plane $y=0$ :

(6.5) $$\begin{eqnarray}\displaystyle \langle (\unicode[STIX]{x1D714}_{y\dagger }^{\ast }-\unicode[STIX]{x2202}_{z}U)^{rms}\rangle _{x,z}=\langle (\langle \unicode[STIX]{x1D714}_{y}\rangle -\unicode[STIX]{x2202}_{z}U)^{rms}(y=0)\rangle _{x,z}=0.58. & & \displaystyle\end{eqnarray}$$

The other flow variables $\unicode[STIX]{x1D70C}_{\dagger }^{\ast },u_{\dagger }^{\ast },w_{\dagger }^{\ast }$ were then obtained using $\unicode[STIX]{x1D6FC}_{\dagger }$ similarly to (6.3a ), (6.3b ), (6.3d ) and are shown with $\unicode[STIX]{x1D714}_{y\dagger }^{\ast }$ in figure 10.

From figure 10(a), it appears that the ‘head’ of the distinctive vorticity structure observed in the experimental CHW (figure 8 a) and predicted by the CHI (figure 9 a) is lacking. The vertical velocity pattern in figure 10(d) has a similar appearance to figures 8(m) and 9(m) but has a larger amplitude (approximately 60 % larger, requiring different colour bar limits of $\pm 0.1$ instead of the previously used $\pm 0.06$ ). The relative magnitude of the vertical velocity $w$ versus the spanwise vorticity $\unicode[STIX]{x1D714}_{y}$ (and hence the streamwise velocity $u$ ) in the 2P-1B analysis is therefore different from the equivalent relative magnitude in the 3P-2B analysis and compares poorly with the experimental observations. Moreover, this analysis is, by construction, incapable of reproducing the three-dimensional structure of the CHW discussed in § 6.1.2, which we have shown depends on the spanwise coordinate $y$ and spanwise velocity $v$ .

Figure 10.

Fastest-growing mode of the 2P-1B analysis (see text for definition of the variables). Note the higher colour bar limits of ${\hat{w}}$ when comparing with figures 8 and 9. As in those previous figures, ten equally spaced black contours span the range of the colour bar for each panel.

6.3 Possible relevance to geophysical field observations

The distinctive tilde shape of the CHW and CHI is reminiscent of the field observations of Geyer et al. (Reference Geyer, Lavery, Scully and Trowbridge2010) in the continuously forced stratified shear flow of the Connecticut River estuary. Through the combined use of acoustic backscattering and high frequency conductivity sensors around the sharp density interface (pycnocline), they observed a near-periodic signal of large sheared tilde-shaped regions of high density gradients with average wavelength around 10 m and height around 1 m (see their figures 2b and 3b, where the $z$ axis has been stretched by factors of 12 and 6, respectively). Those structures are oriented in the same direction as our vorticity structures (the mean shear has the same direction in both our work and their work) and we estimate their maximum angle at the inflection point around $5^{\circ }{-}6^{\circ }$ , somewhat smaller than that of the CHW ( $10^{\circ }$ ) and CHI ( $8^{\circ }$ ). They identified these regions as the braids connecting Kelvin–Helmholtz billows, and argued convincingly, in particular from the conductivity measurements, that these inclined ‘braids’ were the locations of the most vigorous turbulent motions, and that the large density gradients resulted from turbulent mixing within the braids, an inherently high- $Re$ phenomenon which was unlikely to be observed experimentally or numerically. Mashayek & Peltier (Reference Mashayek and Peltier2012a ,Reference Mashayek and Peltier b ) identified secondary instabilities localised in the braid regions connecting Kelvin–Helmholtz billows as a possible explanation for the observations of Geyer et al. (Reference Geyer, Lavery, Scully and Trowbridge2010). Here we have identified a primary instability, the CHI, whose finite-amplitude manifestation is also consistent with these observations. Since the background profiles of the estuary flow are not known with sufficient accuracy to distinguish between flows susceptible to the Holmboe instability or flows susceptible to the Rayleigh instability, we cannot eliminate either of these explanations at this stage. However, unlike the simulated secondary braid instabilities of Mashayek & Peltier (Reference Mashayek and Peltier2012a ,Reference Mashayek and Peltier b ), which occur as a transient event in an initial value problem, the observed CHWs occur in a flow in which the forcing is maintained in time as is the case in the estuary flow.

Preliminary experimental observations of flows in the stratified inclined duct in the intermittently turbulent and fully turbulent regimes revealed a variety of smaller-scale, shorter-lived vorticity structures whose appearance resembles that of the CHW. In particular, two-dimensional ( $y=0$ ) particle image velocimetry measurements at high temporal resolution in the intermittent regime ( $\unicode[STIX]{x1D703}=4^{\circ },Re=940$ – not reported here) strongly suggest that those small-scale structures appear to come from the cascade break-down of larger-scale structures akin to the CHW through successive stretching and splitting.

We, therefore, argue that it is at least possible that the angled structures we observed experimentally (CHW) and predicted theoretically (CHI), through their local coherent intensification of spanwise vorticity due to the spanwise confinement, also have the potential to leave a ‘signature’ at higher $Re$ , even when incoherent, turbulent small-scale motions have been triggered. In particular, in the spirit of the dynamical systems approach discussed in the introduction, we conjecture that the CHW might be sufficiently robust, even at substantially higher $Re$ , to constitute an alternative nucleation site for secondary instabilities and hence turbulent break-down, distinct from the strained braid-like region that occurs between finite-amplitude Kelvin–Helmholtz billows proposed by Corcos & Sherman (Reference Corcos and Sherman1976) and invoked by Geyer et al. (Reference Geyer, Lavery, Scully and Trowbridge2010).