2.1. Set-up and flows

The datasets analysed in this paper were collected in the SID experiment, whose set-up is described in prior publications such as Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a) (see their § 3). The SID sustains a hydraulically controlled exchange flow inside a long duct (of length $L=1350$ mm) of square cross-section (of height and width $H=45$ mm). The duct is inclined at a small angle $\theta$ with respect to the horizontal and connects two large reservoirs initially filled with aqueous salt solutions of different densities $\rho _0 \pm \Delta \rho /2$. The Prandtl number is thus $Pr=\nu /\kappa \approx 700$, where $\nu$ and $\kappa$ are the kinematic viscosity and salt diffusivities, respectively.

Increasing $\theta$ (defined to be positive when it accelerates the flow) and/or the Reynolds number $Re \propto \sqrt {g(\Delta \rho /\rho _0)H}H/\nu$ (entirely set by the density difference here) allows the experimenter to sweep through four qualitatively different flow regimes: from laminar flow with a flat interface, through finite-amplitude Holmboe waves propagating at the interface (this regime is hereafter abbreviated ‘H’), then intermittent turbulence and interfacial mixing (hereafter ‘I’) and to fully developed turbulence and mixing (hereafter ‘T’). These flow regimes have been mapped in the ($\theta,Re$) plane and their transitions have been studied extensively (Meyer & Linden Reference Meyer and Linden2014; Lefauve, Partridge & Linden Reference Lefauve, Partridge and Linden2019; Lefauve & Linden Reference Lefauve and Linden2020).

2.2. Measurements and processing

We consider 15 datasets, each corresponding to a single experiment performed at a given $\theta$ and $Re$. Four belong to the H regime (labelled H1–H4), eight to the I regime (I1–I8) and three to the T regime (T1–T3). Each dataset comprises a time-resolved series of the three-component velocity field $(u,v,w)$ and density field $\rho$ given simultaneously in 3-D volumes $(x,y,z)$, where $u$ and $x$ are the streamwise velocity and coordinate (along the duct), $v$ and $y$ are spanwise and $w$ and $z$ are ‘vertical’ (normal to both $x$ and $y$) in the frame of reference of the tilted duct. The acceleration of gravity $\boldsymbol {g}$ along the ‘true vertical’ is thus tilted with respect to this coordinate system and has components $[g\sin \theta,0,-g\cos \theta ]$ along $(x,y,z)$. See Lefauve et al. (Reference Lefauve, Partridge and Linden2019), figure 1 for a schematic and figures 3 and 4 for a snapshot of $u$ and $\rho$ in each regime.

These 3-D volumes were obtained by the novel laser-sheet-scanning technique described in Partridge et al. (Reference Partridge, Lefauve and Dalziel2019), in which simultaneous stereo particle image velocimetry (PIV) and planar laser-induced fluorescence (PLIF) are performed in successive $x$–$z$ planes. The $u,v,w,\rho$ data obtained at spanwise locations $y = y_i$ ($i= 1, 2, \ldots, n_y$) and respective times $t = t_i$ are subsequently combined in volumes containing $n_y$ planes spanning the cross-section of the duct. This makes the volumes only ‘near-instantaneous’ in the sense that each plane $(x, y_i, z, t_i)$ is separated from the previous one by a small time increment $\delta t = t_i - t_{i-1}$. The advantage of this method over earlier scanning or tomographic methods is the ability to scan relatively large volumes (here typically $200 \times 45 \times 45$ mm$^{3}$) and obtain high $x$–$z$ planar resolutions for both velocity and density. Each experiment typically captures ${\approx }300$ volumes (time snapshots), and each volume typically contains ${\approx }400\times 40\times 80$ velocity vectors in $x,y,z$, respectively.

Instead of the original datasets used and visualised in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a, Reference Lefauve, Partridge and Linden2019) and Partridge et al. (Reference Partridge, Lefauve and Dalziel2019), in this paper we use the slightly modified datasets of Lefauve & Linden (Reference Lefauve and Linden2022a,Reference Lefauve and Lindenb) (hereafter LL22a,b). These modifications are explained in LL22a (see their §§ 3.3–3.5 and figure 1) and are summarised as follows. First, they cropped early transients (in $t$) characterised by a slight net flow through the duct (sloshing between reservoirs) in order to focus on statistically steady dynamics. Second, they cropped the near-wall regions (in $y$ and $z$) in order to discard viscous boundary layers and focus on the interfacial quasi-hyperbolic-tangent ‘free shear layer’ region. Third, they non-dimensionalised the coordinates and flow variables of each individual dataset using: (i) half the density difference $\Delta \rho /2$ (after removing the mean $\rho _0$) such that $-1\le \rho \le 1$; (ii) half the actual resulting shear-layer depth, such that $-1\le z \le 1$ ; and (iii) half the actual (mean) peak-to-peak velocity magnitude, such that $-1\lesssim u \lesssim 1$ with the mean velocity extrema $\langle u \rangle _{x,y} (y=0,z=\pm 1) = \mp 1$. This cropping procedure allows for a meaningful side-by-side non-dimensional analysis of the shear-layer dynamics of all 15 datasets.

The datasets, and the associated codes and 3-D visualisation movies (with the same viewing angle and slices positions as the supplementary movies of positions as used here), can all be freely downloaded from their repository (Lefauve & Linden Reference Lefauve and Linden2022c) .

2.3. Parameters and resolution

The corresponding ‘shear-layer’ Reynolds number $Re$ and bulk Richardson number $Ri_b$ are defined as in LL22a §§ 3.2 and 3.3 as

(2.1)\begin{equation} Re \equiv \underbrace{\frac{\dfrac{ \Delta U}{2}\dfrac{H}{2}}{\nu}}_{\begin{subarray}{c} \text{hydraulics}\\ \text{(input)}\end{subarray}} \ \cdot \ \underbrace{\frac{\delta u}{2}\frac{h}{2}}_{\begin{subarray}{c} \text{shear layer}\\ \text{(output)}\end{subarray}} \equiv \frac{\sqrt{g'H}H}{2\nu} \cdot \frac{\delta u h}{4} \end{equation}

and

(2.2)\begin{equation} Ri_b \equiv \underbrace{\frac{\dfrac{g}{\rho_0}\dfrac{\Delta \rho}{2} \dfrac{H}{2}}{\left(\dfrac{\Delta U}{2}\right)^{2}}}_{\begin{subarray}{c} \text{hydraulics}\\ \text{(input)}\end{subarray}} \ \cdot \ \underbrace{\frac{\dfrac{h}{2}}{\left(\dfrac{\delta u}{2}\right)^{2}}}_{\begin{subarray}{c} \text{shear layer}\\ \text{(output)}\end{subarray}} \equiv \frac{1}{4}\ \cdot\ \frac{ 2h }{(\delta u )^{2}}. \end{equation}

These parameters are consistent with the non-dimensionalisation introduced previously. These parameters consist of (i) a ‘hydraulics’ part based on input parameters, including half the peak-to-peak dimensional velocity scale $\Delta U/2 \equiv \sqrt {g'H} = \sqrt {g(\Delta \rho /\rho _0)H}$, half the duct height $H/2$ and $\nu$, and (ii) a ‘shear-layer’ rescaling based on half the non-dimensional output (measured after the hydraulic non-dimensionalisation) peak-to-peak velocity magnitude $\delta u/2$ (typically ${\approx }0.5\unicode{x2013}1.2$) and shear layer depth $h/2$ (typically ${\approx }0.5\unicode{x2013}0.7$). Note that $Re$ and $Ri_b$ were denoted as $Re^{s}$ and $Ri_b^{s}$, respectively, in LL22a,b to emphasise this specific shear-layer non-dimensionalisation.

The parameters of datasets H1-T3 are listed in table 1. For further properties, such as the mean flows, see LL22a, § 4. As a rule of thumb, increasing levels of turbulence and transitions between the H, I and T regimes, are well described by the product $\theta Re$. The historical names of datasets (H1, $\ldots$ H4, I1, $\ldots$ I8, T1, $\ldots$, T3) were based on increasing values of the product of $\theta$ with the hydraulics (input) Reynolds number $\Delta U H/(4\nu$) in Lefauve et al. (Reference Lefauve, Partridge and Linden2019) rather than on the shear-layer (output) Reynolds number $Re$ used in this paper. Nevertheless, the datasets remain approximately ordered with increasing $\theta Re$.

Table 1. List of the 15 volumetric datasets used, adapted from LL22a's table 1. Note that $\theta$ is in radians in the product $\theta Re$.

The full vector resolution of each dataset in $x,y,z,t$ are given in LL22a, table 3. As a rule of thumb, the vector resolution of the data approaches the Kolmogorov turbulent length scale (marking the end of the inertial range) in $x,z$, but it is coarser in $y$. The temporal resolution (i.e. the time it takes to reconstruct a single volume) is of order 1–4 advective time units, smaller values representing a better ‘freezing’ of the flow, see Partridge et al. (Reference Partridge, Lefauve and Dalziel2019). Fortunately, this does not mean that the structures of interest are distorted in the $x$–$y$ plane by as much as 1–4 length units. First, within the shear layer, substantial portions of the flow have speeds $|u|\ll 1$ (the maximum velocity $|u|\approx 1$ is only reached at the edges of the shear layer). Second, coherent structures tend to evolve more slowly than the background flow at these relatively low $Re$. A archetypal example is the confined Holmboe wave (CHW) in dataset H4, a slow-moving wave which was captured with very little distortion in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a). This presumably extends to turbulent flows, because the magnitude of the perturbation (turbulent) velocities are on average at most 25 % of the mean flow (see figure 1(c) in Lefauve & Linden Reference Lefauve and Linden2022b). Third, our statistics on long time series alternating forward and backward scans help cancel out these distortions, at least their asymmetry around the $x=0$ plane (where vortices concentrate).

Although subject to inherent technical limitations (summarised in Lefauve et al. (Reference Lefauve, Partridge and Linden2019), Appendix A, and Lefauve & Linden (Reference Lefauve and Linden2022b), Appendix B), we show in the following that these 15 state-of-the-art datasets deliver new insights into the time-resolved, 3-D coherent structures of shear-driven stratified turbulence.

3.1. Previous methods and $Q$-criterion

As this paper focuses on Eulerian vortical structures in shear layers, we start by addressing the delicate first step of identifying such ‘vortices’. It is known that identifying a vortex based on a specified threshold of the magnitude of the vorticity vector $\boldsymbol {\omega }=\boldsymbol {\nabla }\times \boldsymbol {u}$ is subjective and generally inappropriate (Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999; Gao, Ortiz-Dueñas & Longmire Reference Gao, Ortiz-Dueñas and Longmire2011). For example, $\boldsymbol {\omega }$ is not generally aligned with the local rotation, and the maximum of $|\boldsymbol {\omega }|$ does not generally coincide with the ‘core’ of a vortex, because vorticity does not discriminate between shear and swirl (rotation).

Several improved Eulerian vortex identification schemes based on the eigenvalues of the velocity gradient tensor $\boldsymbol \nabla \boldsymbol {u}$ have thus been developed since the 1980s, including the $Q$-criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988), $\varDelta$-criterion (Chong, Perry & Cantwell Reference Chong, Perry and Cantwell1990), $\lambda _2$-criterion (Jeong & Hussain Reference Jeong and Hussain1995) and $\lambda _{ci}$-criterion (Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999; Chakraborty, Balachandar & Adrian Reference Chakraborty, Balachandar and Adrian2005). These methods have proved effective and influential to study approximate vortex boundaries in a variety of flows.

Thus, we begin our analysis of vortical structures in figure 1 by a visualisation of vortices based on the popular $Q$-criterion, i.e. the second principal invariant of $\boldsymbol {\nabla }\boldsymbol {u}$ calculated as $Q=\frac {1}{2}(\|\boldsymbol{\mathsf{B}}\|^{2}-\|\boldsymbol{\mathsf{A}}\|^{2})$, where $\boldsymbol{\mathsf{A}} =\frac {1}{2}(\boldsymbol \nabla \boldsymbol {u} + (\boldsymbol \nabla \boldsymbol {u})^{\rm T})$ and $\boldsymbol{\mathsf{B}} =\frac {1}{2}(\boldsymbol \nabla \boldsymbol {u} - (\boldsymbol \nabla \boldsymbol {u})^{\rm T})$ are, respectively, the symmetric (or strain rate tensor) and anti-symmetric (or rotation rate tensor) components of the $\boldsymbol {\nabla }\boldsymbol {u}$. We plot a single snapshot of $Q$-criterion vortices identified by the isosurface $Q=Q_0>0$ (i.e. the local rotation exceeds strain) in a selection of 10 representative datasets having three different tilt angles: $\theta = 5^{\circ }$ in panels (a,b,d,g,j) (datasets H2, H4, I5, I8, T3, with increasing ${Re}$), $\theta = 6^{\circ }$ in panels (c,e,i) (datasets I4, I6, T2) and $\theta = 3^{\circ }$ in panels ( f,h) (datasets I7, T1). The colouring of the isosurfaces denotes the vertical location $z$ and the legend on each plot identifies the dataset, the tilt angle, the Reynolds number and the $Q$ value of the isosurface displayed.

Figure 1. Isosurfaces of the $Q$-criterion (single snapshot) in datasets (a) H2, (b) H4, (c) I4, (d) I5, (e) I6, ( f) I7, (g) I8, (h) T1, (i) T2 and (j) T3. The dataset name, tilt angle $\theta$, shear-layer Reynolds number ${Re}$ (see table 1) and the isosurface $Q$ value are listed on the right of each panel. Colours on the isosurfaces denotes the $z$ position (we show $-0.9< z<0.9$, i.e. the middle 90 % of the shear layer).

Broadly speaking, from the Holmboe to the turbulent regime (i.e. with increasing $\theta {Re}$), vortical structures evolve from individual, disconnected hairpins which start as $\Lambda$-structures without elongated trailing legs (a,b), to groups of hairpins with elongated legs (h–j). The hairpins of low-${Re}$ flows are relatively weak and less obvious (see (c,d)), especially when a corrugation appears on the isosurfaces (which we attribute to a low signal-to-noise ratio), or when the head of the hairpin is broken. However, we show in the following that these hairpin-shaped vortices are indeed robust features of Holmboe waves. In higher-${Re}$ flows, hairpins are stronger and more obvious (higher signal-to-noise ratio), with their head tending toward the edges of the shear layer ($|z|\approx 1$) and their legs stretching in the streamwise direction ($x$).

Further features are worth mentioning. Figure 1(e) shows two large hairpins in each of the upper (isosurfaces shaded in red and labelled ‘U1’ and ‘U2’) and lower (shaded in blue and labelled ‘L1’ and ‘L2’) parts of the shear layer. Figure 1( f–h) (weaker turbulence) show hairpins that are usually asymmetric in the sense that one ‘leg’ is longer than the other. In this paper we interpret ‘quasi-streamwise vortices’ (denoted by ‘SV’; e.g. see arrows SV1, SV2 in (h)) as an extreme form of asymmetric hairpin vortices. Figure 1(i,j) (the most turbulent datasets) show large-scale hairpins coexisting with small-scale, fragmentised structures distributed over the shear layer, which form hairpin packets/forests reminiscent of turbulent boundary-layers.

The supplementary movie 1 available at https://doi.org/10.1017/jfm.2022.588 shows the complete time evolution of figure 1(e–j). From these snapshots (and supplementary movie 1), we hypothesise that hairpin-like coherent vortical structures may be a common and important vortical structure in SID flows.

3.2. RS decomposition

Although the scalar $Q$-criterion identifies general rotational structures, it does not isolate rigid-body rotation from shear (Gao & Liu Reference Gao and Liu2018; Shrestha et al. Reference Shrestha, Nottage, Yu, Alvarez and Liu2021), nor does it provide the orientation of vortical structures.

Recently, a new ‘vortex’ vector called the ‘rortex’ (or ‘liutex’) was proposed by Tian et al. (Reference Tian, Gao, Dong and Liu2018) and Liu et al. (Reference Liu, Gao, Tian and Dong2018) using two successive coordinate rotations to isolate the rigid rotational part in the vorticity field, and directly point in the direction of local rigid-body rotation. The vorticity is decomposed uniquely as

(3.1)\begin{equation} \boldsymbol{\omega} = \boldsymbol{R} + \boldsymbol{S}, \end{equation}

where $\boldsymbol {R}$ is the rortex vector and $\boldsymbol {S}$ is the shear vector. This decomposition is particularly helpful in shear-driven turbulence, as in this paper. An alternative, eigenvector-based definition of the rortex was introduced in Gao & Liu (Reference Gao and Liu2018), followed by its explicit expression with a significantly improved implementation in Xu et al. (Reference Xu, Gao, Deng, Liu and Liu2019), who calculated it as

(3.2)\begin{equation} {\boldsymbol{R}=\left(1-\sqrt{1-\frac{4\lambda_{ci}^{2}}{(\boldsymbol\omega\boldsymbol{\cdot}\boldsymbol{r})^{2}}}\right) (\boldsymbol\omega\boldsymbol{\cdot}\boldsymbol{r})\boldsymbol{r}}. \end{equation}

The direction of the rortex $\boldsymbol {r}$ is the local unit real eigenvector of $\boldsymbol {\nabla } \boldsymbol {u}$, indicating the rotational axis, and $\lambda _{ci}$ is the imaginary part of the complex conjugate eigenvalues of $\boldsymbol {\nabla } \boldsymbol {u}$. This is the practical definition that we apply to our datasets in the remainder of this paper, and we deduce the shear as $\boldsymbol {S} = \boldsymbol {\omega }-\boldsymbol {R}$.

Equation (3.2) is based on the idea that $\boldsymbol \nabla \boldsymbol {u}$ has either one or three real eigenvalues. When there is only one real eigenvalue, the direction of the rortex $\boldsymbol {R}$ is aligned with the associated normalised eigenvector $\boldsymbol {r}$ selected such that $\boldsymbol \omega \boldsymbol {\cdot}\boldsymbol {r} > 0$. The conjugate pair of complex eigenvalues have imaginary parts $\pm \lambda _{ci}$ characterising the rotation about $\boldsymbol {r}$. When there are three real eigenvalues, $\lambda _{ci} = 0$ and thus $\boldsymbol {R}=\boldsymbol {0}$, meaning that all the vorticity is shear without rigid-body rotation. Equation (3.2) shows the relative importance of $\lambda _{ci}$ versus the vorticity projected onto $\boldsymbol {r}$.

The RS decomposition (3.1) can be related to the anti-symmetric rotation rate tensor $\boldsymbol{\mathsf{B}}$ from the above tensor decomposition $\boldsymbol {\nabla } \boldsymbol {u} = \boldsymbol{\mathsf{A}}+\boldsymbol{\mathsf{B}}$. As $B_{ij} = -(1/2) \epsilon _{ijk} \omega _k$ (Batchelor Reference Batchelor1967), we can write it as the sum of a ‘rortex tensor’ $-(1/2) \epsilon _{ijk} R_k$ and a ‘shear tensor’ $-(1/2) \epsilon _{ijk} S_k$, where $\epsilon _{ijk}$ is the Levi-Civita symbol.

Any flow with non-zero vorticity $\boldsymbol {\omega }$ is fundamentally ‘rotational’, in the general sense that infinitesimal fluid elements rotate. In the RS decomposition, vorticity is separated into the rortex vector $\boldsymbol {R}$, about which the fluid is in rigid-body rotation, and the remainder is identified as shear (a rotation with relative motion of fluid particles). For example, a parallel shear layer or a laminar boundary layer (Tian et al. Reference Tian, Gao, Dong and Liu2018) have no rigid-body rotation and $\boldsymbol {\omega }= \boldsymbol {S}$.

Hereafter, we denote the magnitude of the shear vector $|\boldsymbol {S}|=S$, and we refer to the magnitude of the rortex vector $|\boldsymbol {R}|=R$ as the rorticity and any coherent, solid rotational vortical structure it identifies simply as a rortex. We also use the term vortex to refer more generically to vortical structures that have not been unequivocally identified using the RS decomposition, as is the case in all the literature pre-dating 2018.

3.3. Rortex and shear structures

Figure 2 shows isosurfaces of rorticity $R$ (in red) and shear $S$ (in blue) in four datasets (H4, I6, T2 and T3), together with contour plots of $R$ in $y$–$z$ planes (cross-sectional slices, red shades) and a contour plot of $S$ in an $x$–$z$ plane (longitudinal slice, blue shades). Black contour lines highlight the value of each isosurface and its projection on the respective slices. (See supplementary movies 2, 3, 4 and 5 for the complete evolution of $R$ and $S$ as well as the velocity and density information around them for Holmboe, intermittency and turbulent regimes, respectively.)

Figure 2. Visualisation of rortex ($R$) and shear ($S$) structures in datasets (a) H4 (snapshot at time $t_n=180$), (b) I6 (at $t_n=36$); (c) T2 (at $t_n=39$) and (d) T3 (at $t_n=132$). In each panel, $y$–$z$ slices show the $R$ values, whereas the $x$–$z$ planes show the $S$ values. Red 3-D structures represent an $R$ isosurface (respectively $R=0.13, 0.5, 0.85, 0.8$ in a,b,c,d), whereas blue structures represent an $S$ isosurface (respectively $S=1.3, 1.4, 2, 1.8$ in a,b,c,d). The value of the isosurfaces is shown by a black contour on the respective slices. Note that (c,d) share the same colour scales for $R$ and $S$ as panel (b). In all isosurface panels, the spanwise edges of the shear layer $|y|>0.85$ have been blanked for clarity; the lower half $z<0$ has also been blanked in (c,d).

Figure 2(a) in the Holmboe regime (H4) shows a shearing structure (blue isosurface), which is highly reminiscent of the shape of the CHW described in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a) from isosurfaces of the spanwise component of vorticity ($\omega _y$) of the same dataset. This similarity is because the shear dominates the vorticity in the Holmboe regime, as seen by the fact that contour values for $S$ are at least five times larger than that of $R$ (see colour bars), and that the blue isosurface $S=1.3$ is 10 times larger than the red isosurface $R=0.13$. Although weaker, rortices are also observed near the ‘head’ of the Holmboe wave, in a $\Lambda$-shape similar to the $Q$-vortex from figure 1(a,b). The two ‘legs’ of the rortex flank the ‘head’ of Holmboe wave, leading to a hypothesis that the rortex may originate from the localised high-shear regions of the wave.

Figure 2(b) in the low-$Re$ intermittent regime (I6) shows slightly evolving $R$ and $S$ patterns, with a hairpin rortex straddling the shear, as pointed out by the arrows S1 and S2. The snapshot in figure 2(b) is for the same flow and time as that shown in figure 1(e), where two pairs of upper and lower rortices were labelled U1, U2, L1 and L2, respectively. These $S$ structures in figure 2(b) seem to originate from increasingly turbulent symmetric Holmboe waves (having upper and lower counter-propagating modes), as opposed to the asymmetric Holmboe wave of figure 2(a) (only having an upper mode).

Figure 2(c,d) in the turbulent regime (T2 and T3) shows more numerous and smaller-scale structures, which is why the lower half of the shear layer ($z<0$) was omitted for clarity. An apparently robust observation is that rortices still tend to straddle the strong shearing region. The high-shear regions tend to be pushed nearer the top and bottom edges of shear layer, and they are more aligned with the $x$ direction (less tilted) than in the Holmboe and intermittent regimes. Finally, the relative strength of rorticity versus shear also increases (see the colour bars and isosurface values in the caption), indicating an increasing correlation between rortices and turbulent mixing. Although we do not quantify mixing explicitly in this paper, we refer the reader to Lefauve & Linden (Reference Lefauve and Linden2022a) and Lefauve & Linden (Reference Lefauve and Linden2022b) who did so in these 15 datasets using a variety of measures.

Although the vortical structures identified by the threshold of rorticity or by conventional diagnostics (e.g. the $Q$-criterion) are similar in general shape, the detailed morphology and the strength of rotation are different, especially where shear strongly dominates rotation (e.g. in the Holmboe wave).

3.4. Averaged distribution of rorticity and shear

The relative strengths of $R$ and $S$ are explored quantitatively in figure 3. In figure 3(a) we plot, for all 15 experimental datasets, the averaged magnitudes $\langle R\rangle _{xyzt}$ and $\langle S\rangle _{xyzt}$ (where $\langle {\cdot }\rangle _{xyzt}$ represents the average over the whole shear layer volume and time as in LL22a,b) against the product $\theta Re$, our proxy for increasing levels of turbulence. Both $\langle R\rangle _{xyzt}$ ($\diamondsuit$) and $\langle S\rangle _{xyzt}$ ($\triangle$) are nearly constant at ${\approx }0.1\unicode{x2013}0.2$ and ${\approx }1$, respectively, when $\theta {Re}\lesssim 80$ (where $\theta$ is in radians), corresponding to the H and I regimes, but they increase with turbulence beyond this (see the dashed trend lines).

Figure 3. (a) Distribution of $\langle R\rangle _{xyzt}$ and $\langle S\rangle _{xyzt}$ for all datasets, separating H, I and T data by the solid light grey lines ($\triangle$ for $S$; $\diamondsuit$ for $R$; $\square$ for $S'$; $\bigcirc$ for $R'$; the dashed lines are trend lines the dotted line is the fitting curve for $S'$; filled colours of the symbols denote the flow regime and number shown on the right). (b) and (c) are distributions of $\langle S\rangle _{xyt}$ and $\langle R\rangle _{xyt}$ along the $z$ direction for all datasets, respectively, separating H, I and T data in different columns. Colours of the curves and shadings share the same legend in (a). The transparent shadings denote the local variability in time corresponding to one root-mean-square value of $\langle S\rangle _{xy}-\langle S\rangle _{xyt}$ (and similarly for $R$).

To understand this, we also performed the RS decomposition on the fluctuating velocity $\boldsymbol u'=\boldsymbol u - \bar {\boldsymbol u}$ where $\bar {\cdot} = \langle {\cdot } \rangle _{t}$ is the spatially varying temporal average, giving the underlying parallel shear flow $\bar {u}(x,y,z)$ (because $\bar {v},\bar {w}\approx 0$). Figure 3(a) also shows the resulting fluctuating shear $S'$ ($\square$) and rorticity $R'$ ($\bigcirc$). First, we find that $R \approx R'$ (the symbols are nearly indistinguishable) whereas $S' \ll S$, as expected in the presence of mean background shear $\partial _z \bar {u}$. Second, we find that $S'$ is only about two to three times larger than $R'$ (a weaker dominance compared with that of $S$ over $R$), and $S'$ seems to increase fairly linearly with $\theta {Re}$ even before the transition to fully turbulent flow at $\theta {Re} \approx 80$ (see the green dotted trend line). These observations suggest that the background shear plays an important role, but the details are beyond the scope of the present study, which focuses primarily on rortex structures. The remainder of this paper thus adopts the RS decomposition of the total velocity, as in the original papers of Liu et al. (Reference Liu, Gao, Tian and Dong2018) and Gao & Liu (Reference Gao and Liu2018).

To examine the strength of $S$ and $R$ along the $z$ (vertical) direction, we plot $x,y,t$ averages in figure 3(b,c), separating $S$ and $R$, respectively, as well as the datasets belonging to the H regime (top row), I regime (middle row) and T regime (bottom row) for clarity. The time variations around the averaged value measured by the local root mean square (r.m.s) are displayed as transparent shadings underlying each curve.

Figure 3(b) shows that the symmetric Holmboe flows (datasets H1 and H3) have their peak $\langle S\rangle _{xyt}$ at $z>0$; these flows feature two trains of counter-propagating waves due to the velocity interface $\langle u \rangle _{xyt}=0$ and density interface $\langle \rho \rangle _{xyt}=0$ being nearly coincident around $z\approx 0\unicode{x2013}0.15$ (see LL22a figure 3). By contrast, the asymmetric Holmboe flows (H2 and H4) have their peak $\langle S\rangle _{xyt}$ at $z<0$; these flows feature a single train of waves due to the velocity interface $\langle u \rangle _{xyt}=0$ being slightly offset from the density interface $\langle \rho \rangle _{xyt}=0$, the latter of which is around $z \approx -0.2$ (see LL22a figure 3). This peak is due to the ‘body’ of the CHW structure described in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a), and observed in figure 2(a). We also see an apparent plateau at $0< z<0.5$ in these datasets (see the ‘P’ arrow), presumably due to the ‘head’ of the CHW. With increasing levels of turbulence (I regime, middle row), this plateau in $\langle S\rangle _{xyt}$ appears to develop into another peak (see the ‘peak2’ arrow). Both the former and the newer peaks then tend to move closer to the edges ($z=\pm 1$) of the turbulent shear layer (T regime, rightmost column). Their values ${\approx }1.5\unicode{x2013}2$ is comparable to those in the H regime. A further interesting observation is that the temporal r.m.s. fluctuations of $\langle S\rangle _{xyt}$ (width of the transparent shading) increase from the H to the I regime but then decrease in the T regime. This trend reflects the high fluctuations associated with intermittency.

Figure 3(c) shows that $\langle R\rangle _{xyt}$ is nearly uniform in $z$ across the entire shear layer, with H1 and H3 having higher values, presumably due to their higher $Re$ than H2 and H4 ($Re\approx 400$ versus $200$). In the ‘late’ I regime (e.g. I7, I8) and in the T regime, a broad peak in $\langle R\rangle _{xyt}$ is centred in the middle the shear layer, with peak value that increases approximately monotonically. Finally, unlike the shear, the rortex experiences equally high or even higher temporal fluctuations in the T regime (compared with the I regime), possibly caused by the breakdown and interaction of rortices. This suggests that the emergence and increasing importance of the rortex are fundamental aspects of turbulence dynamics and mixing, justifying our greater focus on $\boldsymbol {R}$ (vortical structures) than $\boldsymbol {S}$ (shear structures) in the remainder of this paper.

4.1. WCA and orientation probability distribution functions (pdfs)

In order to remove noise and to give stronger weight to stronger rortices, we first apply a ‘conditional sampling’ method on the fields $\boldsymbol {R}(\boldsymbol {x},t)$ (containing $0.3 \times 10^{9}$ to $1 \times 10^{9}$ points per dataset, depending on the spatial resolution and length of the time series). We condition these data by the following formula at all $(\boldsymbol {x},t)$:

(4.1) \begin{equation} \boldsymbol{R} = \left\{\begin{array}{@{}ll} \boldsymbol{R}, & \text{if}\ R/R_{rms} \ge k_{th}, \\ 0, & \text{if}\ 0 < R/R_{rms} < k_{th}, \end{array}\right. \end{equation}

where $k_{th}$ is a tunable threshold level below which the data are discarded, and $R_{rms}$ is the r.m.s. of all non-zero $R$ values (before the above conditioning). Statistical processing is then performed on all non-zero $R$ after the above preconditioning.

Next, we use ‘orientation’ probability density functions (pdfs) to quantify the likelihood of the orientation of specific vectors, measured by their frequency distribution in our datasets. For any non-zero vector $\boldsymbol {V}$ (where $\boldsymbol {V}$ may represent $\boldsymbol {R}$, $\boldsymbol {S}$, etc.) we define the angles between $\boldsymbol {V}$ and the planes ${{x}}_\perp$ and ${{z'}_\perp }$ as

(4.2a–c)\begin{equation} \alpha = \angle(\boldsymbol{V},\boldsymbol{\hat{x}})-90,\quad \beta = \angle(\boldsymbol{V},\boldsymbol{\hat{z}}')-90, \quad \boldsymbol{\hat{z}}'=\frac{-\boldsymbol{g}}{g}, \end{equation}

where $\angle (\boldsymbol a,\boldsymbol b)\equiv \arccos (\boldsymbol a \boldsymbol {\cdot} \boldsymbol b/(|\boldsymbol a |\, |\boldsymbol b|))\in [0,180]$ is the angle in degrees between $\boldsymbol a$ and $\boldsymbol b$. The unit vectors are defined as follows: $\boldsymbol {\hat {x}}, \boldsymbol {\hat {z}}$ are the unit vectors along the streamwise ($x$) and wall-normal ($z$) direction of the duct; ${{x}_\perp }$ indicates the plane normal to $\boldsymbol {\hat {x}}$; $\boldsymbol {\hat {z}}'= \cos \theta \boldsymbol {\hat {z}}-\sin \theta \boldsymbol {\hat {x}}$ is the ‘true vertical’ unit vector (in the opposite direction of gravity); and ${{z'}_\perp }$ is the ‘true horizontal’ plane normal to the $\boldsymbol {\hat {z}}'$. These coordinate systems and angles (with their sign) are illustrated in figure 4 and our angle notation is summarised in table 2.

Figure 4. Definition of the angles (a) $\alpha = \angle (\boldsymbol {V},\boldsymbol {\hat {x}})-90$ and (b) $\beta = \angle (\boldsymbol {V},\boldsymbol {\hat {z}}')-90$, where $\boldsymbol {V}$ may represent $\boldsymbol S$, $\boldsymbol R$, or $\boldsymbol {\nabla }\rho$, leading to the six angles summarised in table 2. The grey cones correspond to the locus of possible $\boldsymbol {V}$ for a single value of (a) $|\alpha |$ and (b) $|\beta |$. The ‘true horizontal’ plane ${z'}_\perp$ (in green in b) is normal to the opposite direction of gravity $\boldsymbol {\hat {z}}'$, whereas the plane ${{x}_\perp }$ (in green in a) is normal to $\boldsymbol {\hat {x}}$. Finally, $\theta$ is the tilt of the duct with respect to the horizontal direction (the convention is that $\theta >0$ indicates that the flow is forced).

Table 2. Summary of the angles discussed in this paper (refer to definition (4.2a–c) and figure 4). The bottom row indicates the figures in which their distributions are shown.

Finally, to extract detailed rortex morphology from orientation pdfs, we weigh the occurrence of each value within a particular interval (histogram value) by the local value of the ‘rorstrophy’ (the squared rorticity) $R^{2}(\boldsymbol {x},t)$. This weight gives more importance to occurrences that locally coincide with high rortex values. Practically, the averaged orientation pdf $N(i,k_{th})$ of any angle $\alpha$ or $\beta$ at a value of $i \in [ -90^{\circ }, 90^{\circ }]$ and for a given conditional threshold of $k_{th}$ is calculated by

(4.3)\begin{equation} N(i,k_{th}) = \frac{\sum_{l=1}^{n_t}\sum_{j=1}^{n_i}{R(j,l)^{2}}}{n_t},\quad {\rm with}\ R(j,l) = 0 \quad \mathrm{if}\ \frac{R(j,l)}{R_{rms}}< k_{th}, \end{equation}

where $n_i$ is the count of the occurrences when the angle under consideration belongs to the interval (bin) $i\pm \delta i$, $j$ is the index for all $(x,y,z)$ data points belonging to this interval and $l$ is the time index sweeping through the $n_t$ ‘frames’ (volumes) in the dataset. Note that $\int _{-90}^{90} N \,\textrm {d}i = \langle R^{2} \rangle _{xyzt}$, i.e. the area under the curve of $N$ gives the time- and volume-averaged ‘rorstrophy’ satisfying the threshold $R/R_{rms} \ge k_{th}$. If $k_{th}=0$, the original rortex field and all existing rortices are considered, following (4.1).

The process above is a weighted conditional average (WCA): weighted by $R^{2}$ and conditioned by selecting $R/R_{rms}\ge k_{th}$. This same approach will be applied to angle frequency distributions and used to study how progressively stronger rortices are aligned with respect to ${{x}_\perp }$ and ${{z'}_\perp }$. By analogy, we also extend WCA to $\boldsymbol {S}$ (weighting by $S^{2}$ and conditioning by $S/S_{rms} \ge k_{th}$).

Before showing our results on the orientation of $\boldsymbol {R},\boldsymbol {S}$, it is worth studying the ‘volume fraction of rortices’ $f$ resulting from our conditional sampling method in (4.1) alone, without weighting. Figure 5(a) shows how the global rortex volume fraction $\langle f\rangle _{xyzt} \in [0,1]$ (the time- and volume-averaged ratio of points satisfying $R/R_{rms}\ge k_{th}$) decreases with increasing threshold level $k_{th}$. The semi-log axes and the exponential fit (dashed line) reveal that $\langle f\rangle _{xyzt}$ decreases approximately exponentially with $k_{th}$ with decay constant $\approx 1.4$. The intercept of $0.755$ at $k_{th}=0$ means that before conditioning, approximately three quarters of the shear layer volume has non-zero rortex (i.e. the velocity gradient tensor $\boldsymbol {\nabla }\boldsymbol {u}$ has a component of solid-body rotation). All 15 datasets display a similar behaviour at small $k_{th}$, but their curves spread out significantly for $k_{th}\gtrsim 2$. Increasing turbulence (curve colour transitioning from blue to red) reduces the decay rate at high $k_{th}$, confirming the intuition that turbulence leads to more frequent extreme rortex values. Though not shown here, the volume fraction of $R'$ (based on the fluctuating velocity $\boldsymbol {u}'$) is indistinguishable from that of $R$.

Figure 5. (a) Volume fraction $\langle f \rangle _{xyzt}$ of non-zero rorticity $R$ under different conditional threshold levels $k_{th}$ as defined in (4.1). Dashed line denotes the exponential fit (note the semi-log axes). (b) Variation of the volume fraction $\langle\, f \rangle _{xyt}$ of $R$ along $z$ for threshold $k_{th}=1$.

Next, figure 5(b) shows the vertical distribution of volume fraction $\langle f\rangle _{xyt}(z)$ (averaged in horizontal planes and time, but not $z$) at threshold $k_{th}=1$. The result shows that rortices mainly concentrate in the middle region of the shear layer in the intermittent and turbulent datasets, which agrees with the plots of $\langle R\rangle _{xyt}(z)$ in figure 3(c) without conditional thresholds. In particular, T2 and T3 have a robust ${\approx }20\,\%$ fraction for $|z|\lesssim 0.5$, which tapers off to only ${\approx }5\,\%$ at the edges $|z| =1$, justifying our cropping of the original datasets (see § 2.2) to restrict our attention to the $|z|\le 1$ ‘shear layer’ containing the turbulent rortices of interest. However, the H and early I regimes show slightly different tends. While datasets H1, H3 and I1 (see arrows) have a similar distribution to I6–T3, datasets H2, H4 and I4 have their minimum rortex fraction near the centre of the shear layer and their maximum at the edges. This indicates that asymmetric Holmboe waves (found at high $\theta$, low $Re$) and symmetric Holmboe waves (found at low $\theta$, high $Re$), and their respective weakly intermittent turbulent counterpart (I1–I4) have inherently different rortex distribution and dynamics along $z$. This echoes the findings of LL22a (§ 6.4) that high-$\theta$, low-$Re$ turbulence is characterised by more overturning motions and less extreme shear-dominated enstrophy than low-$\theta$, high-$Re$ turbulence.

4.2. Inclination of shear structures

Although vortical structures are the focus of this study, we start with a brief study of the inclination angles of the shear vector $\boldsymbol {S}$, remembering its background role (in fact, dominant in magnitude) in shear-driven, stratified turbulence. Figure 6 shows the pdfs of angles $\alpha _S$ and $\beta _S$ (as summarised in table 2). We applied our WCA method with threshold $k_{th}=2$ (thus excluding most of the modest shear associated with the mean flow $\partial _z \bar {u} = O(1)$), and weight $S^{2}$ (giving emphasis to large shear events). We see that large shear events are generally perpendicular to both $\boldsymbol {\hat {x}}$ and $\boldsymbol {\hat {\boldsymbol z}}'$ and thus primarily along $\boldsymbol {\hat {y}}$, indicating the dominance of the spanwise component of vorticity, which motivated the use of $\omega _y$ in our previous study (Lefauve et al. Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a). Increasing turbulence (blue to red curves) widens the pdf of both $\alpha _S$ and $\beta _S$, revealing increasingly 3-D shearing structures. These results are robust for other thresholds $k_{th}$.

Figure 6. Frequency distributions (orientation pdfs) of angles (a) $\alpha _S$ and (b) $\beta _S$ with threshold $S>k_{th}=2$, weighted by $S^{2}$. All distributions are normalised to have unit integral (like any pdf) for easier comparison between datasets.

4.3. Inclination of rortex structures

In previous studies, the orientation of vortices was obtained through the orientation of the vorticity vector (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019), two-point spatial correlations of velocity (Williams Reference Williams2014), geometrical features of some isosurfaces/contours based on conventional vortex diagnostics (Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999) and vortex vectors (Arun et al. Reference Arun, Sameen, Srinivasan and Girimaji2019). According to Arun et al. (Reference Arun, Sameen, Srinivasan and Girimaji2019), the orientation of vortical structures has been shown to be closely related to turbulence production in compressible mixing layers. We now turn to a comprehensive analysis of the orientation pdfs of the vortex vector $\boldsymbol {R}$ angles in our stratified mixing layer, which show different and more subtle behaviours than those of the shear.

In figure 7, we plot the pdfs of $\alpha _R$ (a) and $\beta _R$ (b). Each of the 15 datasets are plotted in separate subpanels, and arranged according to their positions in the $(\theta,{Re})$ plane, in order to draw connections between the observed vortex statistics and the two key flow parameters. Furthermore, the pdfs are weighted with $R^{2}$, and we use curves of increasingly dark colour (blue in (a) and red in (b)) to indicate increasingly high conditional threshold levels $k_{th}$ representative of more extreme rortices (note the semi-log scale). To facilitate comparisons of magnitudes between the 15 subpanels, we use the same axis limits in all subpanels and normalise all pdfs such that the integral of each over the interval $-90^{\circ }$ to $90^{\circ }$ gives the all-time and volume-averaged square norm of the conditioned $\langle R^{2}\rangle _{xyzt}$ (rather than 1, as in figure 6). The dashed blue diagonal lines in the background are a fit of the observed ‘overturn fractions’ in these datasets. These were calculated in LL22a (see their § 6) as the time- and volume-averaged fraction of the flow that experiences density overturnings ($ {\partial }_z \rho >0$), and the best fit in $(\theta, Re)$ space was shown to scale with $\theta ^{3.17}Re^{1.75}$.

Figure 7. WCA frequency distribution (orientation pdfs) of angles (a) $\alpha _R$ and (b) $\beta _R$ in all 15 datasets, arranged in subpanels in the $(\theta,Re)$ plane (symbols indicate precise parameters of each dataset). We use $R^{2}$ weights and plot increasingly high conditional threshold levels $k_{th}\in [0:0.5:4]$ in darker shades. All panels have the same axis limits as labelled in the H2 subpanel. The faint dashed blue lines in the background indicate the fitting of the ‘overturn fraction’ of LL22a (their figure 8c). Note the semi-log scale in all subpanels.

4.3.1. With respect to the $x_\perp$ plane

In figure 7(a) we observe the following.

1. (i) All $\alpha _R$ pdfs are statistically symmetric around $0^{\circ }$ to an excellent approximation.

2. (ii) Stronger rortices (i.e. high $k_{th}$, darker blue lines) nearly always have the same peak angle as weak rortices (i.e. low $k_{th}$, light blue lines), identified by black arrows, except in a few datasets when an additional minor peak at $\alpha _R=0^{\circ }$ appears at higher $k_{th}$, identified by blue arrows.

3. (iii) A peak angle (maximum of the weighted pdf) at $\alpha _R\approx$ $\pm 60^{\circ }$ appears in the bottom-left and top-right corners of the $(\theta,{Re})$ plane, i.e. either at low $\theta$ and low ${Re}$ or high $\theta$ and high ${Re}$.

4. (iv) In contrast to (iii), a wider and more uniform pdf across $\alpha _R\in [-60^\circ,60^\circ]$ appears in the upper-left and lower-right corners, i.e. either low $\theta$ and high ${Re}$ or at large $\theta$ and low ${Re}$.

These features are clearly visible in figure 8(a), which shows the evolution of the peak values of $\alpha _R$ for values of $k_{th}$ ranging from 1 (weak rortices) to 3 (strong rortices). Larger symbols denote stronger (and, thus, more significant) peaks, and open symbols refer to the fluctuating rortex data $\boldsymbol {R}'$ (without $\bar {\boldsymbol{u}}$), which are essentially similar to the full rortex data $\boldsymbol {R}$.

Figure 8. Peak angle(s) of the (a) $\alpha _{R}$ and (b) $\beta _{R}$ pdfs, automatically extracted from figure 7 for $k_{th} = 1$:0.5:3 ($\bigcirc$). We use $R^{2}$ weights and plot increasingly high conditional threshold levels $k_{th}\in [1:0.5:3]$ in darker colours (as for figure 7). Open symbols ($\vartriangleright$) denote the peak angle based on $\boldsymbol {R}'$ with the same $k_{th}$ range showing the influence of background shear. Symbol size indicate the strength of the peak angles (proportional to the square root of the ordinates $\sum R^{2}$ in figure 7).

We conclude that both $Re$ and $\theta$ influence the horizontal orientation of rortices in subtle ways. In turbulent flows (T regime), both weaker and stronger rortices are primarily inclined at an angle $\alpha _R\approx \pm 60^{\circ }$ to the ${{x}_\perp }$ plane, largely independent of $\theta$ (though it varies across a wider range than $Re$ does). However, at lower $Re$ values (H, I regimes), the evolution of $\alpha _R$ with $Re$ (at fixed values of $\theta$, i.e. along horizontal lines) varies depending on the value of $\theta$. For example, at $\theta =5^{\circ }$, the unimodal or uniform pdf becomes bimodal at higher $Re$ (compare the evolution H2 $\rightarrow$ H4 $\rightarrow$ I5 $\rightarrow$ I8 $\rightarrow$ T3), whereas almost the opposite happens at $\theta =2$ (evolution H3 $\rightarrow$ I1$\rightarrow$ I2 $\rightarrow$ I3). This difference again indicates that Holmboe waves at low versus high $Re$ (here coinciding respectively with the asymmetric and symmetric type of Holmboe waves) have different properties. Rortices found in high-$Re$ (symmetric) Holmboe waves H1, H3 are more akin to those found in turbulent flows. Finally, stronger rortices can exhibit a trimodal pdf (see blue arrows) in some intermittent flows (I regime). The middle peak at $\alpha _R=0^{\circ }$ (perpendicular to the $x$-axis) suggests transverse rortices or the heads of hairpin rortices, which become less dominant under stronger turbulence.

4.4. Inferred morphology

Based on the above descriptions, we now draw in figure 9 representative schematics of the morphology of rortices typical of each flow regime. Top views (in the $x$–$y$ plane) are shown in the top row (a–d), whereas side views (in the $x'$–$z'$ plane, where $x'$ is normal to $z'$) are shown in the bottom row (e–h). The strongest magnitudes $R$ are always found inside the structures and are denoted by darker shades of red (typically corresponding to $k_{th}=1$ and $2.5$).

Figure 9. Schematics of the evolution of rortex structures in increasingly turbulent flows: (a–d) top view; (e–h) side view. Left to right: ${H_{asym}}$ and ${H_{sym}}$ denote asymmetric and symmetric Holmboe regimes; I and T denote intermittently turbulent and fully turbulent regimes, respectively. Light red denotes weak rorticity $R$ (low $k_{th}$), dark red denotes strong $R$ (high $k_{th}$). Dashed lines indicate the typical angles of inclination found in figures 7 and 8.

The angle $\alpha _R$ of the rortices in (b–d) is close to that of hairpin vortices observed in unstratified turbulent boundary layers having $\alpha _R \approx 40\sim 60^{\circ }$ (Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999). For the vertical inclination, the rortices in the T regime (h) have an average $|\beta _R|\approx 11^{\circ }$ which agrees well with the hairpin rortices observed in the direct numerical simulations of stratified shear layers in Watanabe et al. (Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019) (who reported an inclination of 14$^{\circ }$) and those hairpin legs observed in unstratified wall-bounded flows in Haidari & Smith (Reference Haidari and Smith1994) and Zhou et al. (Reference Zhou, Adrian, Balachandar and Kendall1999) (who reported ${\approx }8^{\circ }$). However, these angles are lower than the average tilt of hairpin in unstratified turbulent boundary layers, which are typically ${\approx} 45^{\circ }$ (Head & Bandyopadhyay Reference Head and Bandyopadhyay1981; Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999). Interestingly, the ‘heads’ of our hairpin rortices are barely more inclined than the legs, which suggests that the ‘lift-up’ of the head may be inhibited by stratification. The legs of the symmetric Holmboe rortices ( f) are almost horizontal, showing a much clearer connection with intermittent and turbulent structures than the asymmetric Holmboe ones (e).

The schematics shown in figure 9 are ideal models of symmetric $\varLambda$ and hairpin rortices representative of the evolution of global $\alpha _R$ and $\beta _R$ statistics on either side of the turbulent transition. Instantaneous rortices in high-$Re$, turbulent, unsteady flows, are more likely to feature ‘broken’, asymmetric hairpins, such as quasi-streamwise or cane-like rortices. Nevertheless, we will show with case studies in § 6 that these ideal models are representative of actual rortices observed in instantaneous snapshots.

5.1. Distribution of the density gradient and relation to the rortex

In figure 10 we plot the WCA pdfs of the angle between $\boldsymbol {\nabla }\rho$ and the ‘true horizontal’ plane ${{z'}}_\perp$, i.e. $\beta _\rho = \angle (\boldsymbol {\nabla }\rho,\boldsymbol {\hat {z}}')-90$. We compare the pdfs under two different weights: in (a) with $|\boldsymbol {\nabla }\rho |^{2}$ to highlight the strongest density gradients (we refer to this as $\beta _{\rho 1}$); and in (b) with the squared rortex magnitude $R^{2}$ to highlight the orientation of density interfaces coinciding with strong rortices (we also impose the threshold $R>k_{th}=2$, and refer to this as $\beta _{\rho 2}$).

Figure 10. Frequency distribution (orientation pdfs) of the WCA vertical angles of the density gradient: (a) $\beta _{\rho 1}$ weighted by $|\boldsymbol {\nabla }\rho |^{2}$ without conditional threshold; (b) $\beta _{\rho 2}$ weighted by $R^{2}$ with a threshold $R>k_{th}=2$. (c) Peak angle (maximum of the weighted pdf) $\beta _{\rho 1}$ (empty $\bigcirc$) and $\beta _{\rho 2}$ (filled $\bigcirc$). We also add the horizontal angles $\alpha _{\rho 1}$ with weight $|\boldsymbol {\nabla }\rho |^{2}$ (empty $\square$) and $\alpha _{\rho 2}$ with weight $R^{2}$ (filled $\square$). Note the log vertical scale. Symbol size indicates the relative strength of the peak, i.e. the value of its ordinate in (a,b).

In (a) we find, unsurprisingly for stably stratified flows, that the strongest density gradients overwhelmingly point downwards, thus $\beta _{\rho 1}>0^{\circ }$, with a small deviation from the perfect ‘true’ vertical of $90-\beta _{\rho 1} \approx 5^{\circ }$.

By contrast, in (b) we observe a much broader pdf. The rorstrophy-weighted density gradients are much less vertical or downward-pointing, and a significant fraction point upward (see the grey shaded area ‘overturned’ for $\beta _{\rho 2}<0^{\circ }$). Density gradients that are co-located with strong rortices thus appear much more susceptible to being distorted and even overturned. The peak values of this pdf tell us about the extent of this distortion process. We observe an almost monotonic evolution from modest deviations of $90 - \beta _{\rho 2} \approx 7\unicode{x2013}10^{\circ }$ in the H datasets to more substantial deviations of ${\approx }16^{\circ }\unicode{x2013}25^{\circ }$ in the I and T datasets (see dashed trend line). In the latter datasets, and especially in T2 and T3, the left-hand tails of the pdfs decay much more slowly, with significant overturns (we recall that the overturn fraction was shown by light blue contours in figure 7, giving a typical ${\approx }3\,\%\unicode{x2013}5\,\%$ of overturned fluid in T2 and T3). These observations are robust at different $R$ conditional threshold $k_{th}$.

Figure 10(c) summarises the evolution of the peak angles of $\beta _{\rho 2}$ (colour-filled $\bigcirc$), noting in passing the difference with the evolution (or lack thereof) of $\beta _{\rho 1}$ (empty $\bigcirc$). The trend of a monotonically decreasing $\beta _{\rho 2}$ (or increasing $90 - \beta _{\rho 2}$) with $\theta Re$ is very clear. We also add the peak in horizontal angle $\alpha _\rho = \angle (\boldsymbol {\nabla }\rho,\hat {\boldsymbol {x}})-90$ ($\alpha _{\rho 1}$ weighted by $|\boldsymbol {\nabla }\rho |^{2}$ and $\alpha _{\rho 2}$ weighted by $R^{2}$) to show that they both have a consistent peak at 0$^{\circ }$.

Interpreting these results, we note that although both shear and rorticity increase with turbulence (as was shown in figure 3), the decrease in $\beta _{\rho 2}$ is presumably caused by the increasing dominance of nearly horizontal rortices having $\beta _R\approx 0^{\circ }$ (see figure 8). These strongly rotating hairpin ‘heads’ (see figure 9) lift up and overturn the flow around the $y$ axis. The symmetric Holmboe data (with high-$Re$, H1 and H3) are in this sense ‘pre-turbulent’, especially H1 (marked by a red arrow) which peaks at ${\approx} 70^{\circ }$. However, the asymmetric Holmboe data (with low-$Re$, H2 and H4, labelled with blue arrows) are again different, because the rortex remains relatively weak compared with the shear and does not visibly influence the density gradient. The density interfaces with strongest rortices are indeed inclined at the same angle as the rortices themselves, as evidenced by the fact that $\beta _{\rho 2}$ peaks in the range $40$–$50^{\circ }$ (labelled by blue arrows), which is very close to the peak $\beta _R$ in figure 8(b).

To further study the relationship between rortices and density interfaces, we define a last angle $\phi =\angle (\boldsymbol {R},\boldsymbol {\nabla }\rho )$. The evolution of its peak angle $\phi _p$, extracted from its WCA distribution, is shown in figure 11. The symbol size denotes the peak height, whereas the shading denotes the width (or spread) of the distribution at half height (see top right insert). Again we compare distributions under two different weights: $|\boldsymbol {\nabla }\rho |^{2}$ (white circles and green shading) and $R^{2}$ (red symbols and grey shading).

Figure 11. Peak angle $\phi _p$ between $\boldsymbol {R}$ and $\boldsymbol {\nabla }\rho$ based on WCA distributions. Symbol size denotes the peak height $h$ (as shown schematically in the insert). Shadings denote the width $W$ of the distribution at half peak height $h/2$ (see insert). Red-filled $\bigcirc$ and grey shading have weight $R^{2}$; white-filled $\bigcirc$ and green shading have weight $|\boldsymbol {\nabla }\rho |^{2}$. All statistics weighted by $R^{2}$ have conditional threshold $k_{th}=2$.

First, all peak angles are $\phi _{p}\approx 90^{\circ }$, revealing a new piece of information: rortices and density gradients are most frequently perpendicular, in all datasets, regardless of statistical weight. Second, distributions of $\phi$ weighted by $|\boldsymbol {\nabla }\rho |^{2}$ (green shading) are narrower than those weighted by $R^{2}$ (grey shading), having a typically spread of ${\pm }10^{\circ }$ versus ${\pm }30^{\circ }$ (in some datasets even higher, e.g. H2, H4, I4). This observation brings an important nuance to our above conclusion: while the strongest $\boldsymbol {\nabla }\rho$ are indeed frequently perpendicular to $\boldsymbol {R}$, the strongest $\boldsymbol {R}$ are less frequently perpendicular to $\boldsymbol {\nabla }\rho$. This subtle asymmetry in the relation between $\boldsymbol {R}$ and $\boldsymbol {\nabla }\rho$ is important and understandable: we expect strong density gradients to generate rortices by baroclinic torque, but we do not expect all rortices, especially the strongest ones, to be generated by this mechanism. Our stratified layers remain dominated (driven) by shear, and most rortices can be expected to be a product of instabilities that grow by extracting energy from the mean shear, especially in the centre of the shear layer where the fluid is partially mixed and density gradients are weak.

Third, although we have argued that the spread of distributions weighted by $\boldsymbol {\nabla }\rho$ (green shading) is relatively ‘narrow’, it broadens somewhat from approximately $90\pm 5^{\circ }$ (H data) to $90\pm 15^{\circ }$ (T data) as turbulence levels increase. This evolution shows that rortices become less perpendicular to even the strongest density gradients in increasingly turbulent flows, probably as a result of weaker stratification ($Ri_b \approx 0.15$ in T data) and, thus, of a weaker feedback of density in the momentum equation.

5.2. Vertical variation of the rortex–density and shear–density relationships

We now turn to the strength of $|\boldsymbol {\nabla }\rho |$ and of its relationships with $\boldsymbol {S}$ and $\boldsymbol {R}$ along the $z$ direction. Figure 12(a–c) shows $\langle |\boldsymbol {\nabla }\rho |\rangle _{xyt}(z)$, segregating the H, I and T data in different columns. We see that the initially sharp density interface broadens (from H $\rightarrow$ I $\rightarrow$ T) and ends up in (c) becoming partially mixed across $-0.5 \lesssim z \lesssim 0.5$, flanked by two weaker interfaces. This evolution is similar to that of the shear $\langle S \rangle _{xyt}(z)$ (previously shown in figure 3b).

Figure 12. Vertical distribution (along $z$) of the average (a–c) $\langle |\boldsymbol {\nabla }\rho |\rangle _{xyt}$, (d–f) $\langle |\boldsymbol S\times \boldsymbol {\nabla }\rho |\rangle _{xyt}/\langle |\boldsymbol {\nabla }\rho |\rangle _{xyt}$ and (g–i) $\langle |\boldsymbol R\times \boldsymbol {\nabla }\rho |\rangle _{xyt}/\langle |\boldsymbol {\nabla }\rho |\rangle _{xyt}$ segregating the H, I and T regimes in different columns. The dashed lines in (d–f) and (g–i) correspond to $\langle S\rangle _{xyt}$ and $\langle R\rangle _{xyt}$, respectively (repeating some information from figure 3b,c).

Figure 12(d–f) shows the average magnitude of their cross product $|\boldsymbol {S} \times \boldsymbol {\nabla }\rho |_{xyt}$ scaled by $\langle |\boldsymbol {\nabla }\rho |\rangle _{xyt}$. The profile is very close to $\langle S\rangle _{xyt}$ (shown by dashed lines), indicating that $\sin (\angle (\boldsymbol {S},\boldsymbol {\nabla }\rho ))\approx 1$, i.e. that $\boldsymbol {S}$ (primarily along $y$) is approximately perpendicular to $\boldsymbol {\nabla }\rho$ (primarily along $z$).

Figure 12(g–i) shows $|\boldsymbol {R} \times \boldsymbol {\nabla }\rho |$ scaled by $\langle |\boldsymbol {\nabla }\rho |\rangle _{xyt}$. The interaction of rortex and density gradient is distributed more evenly across $z$, with weak peaks that neither reflect the peak of $\langle R\rangle _{xyt}$ (in dashed lines) nor the peak of $\langle |\boldsymbol {\nabla }\rho |\rangle _{xyt}$ (see first row of this figure). This proves that rortices interact strongly with density gradients across the whole shear layer, rather than just at a single sharp density interface or at the two weaker interfaces on either edge of the shear layer. Although the rortex and density gradient are frequently nearly perpendicular (as we have shown in the previous figure 11), the departure from this general tendency is substantial enough, especially in turbulent flows, that we do not find any region with strongly peaked $|\boldsymbol {R} \times \boldsymbol {\nabla }\rho |$ across the shear layer.

The cross product is a particularly useful indicator because it considers not only the strength of $\boldsymbol R$ and $\boldsymbol {\nabla }\rho$ but also their angle $\phi$ (recall $|\boldsymbol R\times \boldsymbol {\nabla }\rho | \equiv |\boldsymbol R| |\boldsymbol {\nabla }\rho | \sin \phi$). Rortices perpendicular to the density gradient ($|\sin \phi |=1$) have a potential energy cost proportional to $|\boldsymbol {\nabla }\rho |$ and $Ri_b$, whereas those parallel to the density gradient ($\sin \phi =0$) are energetically neutral. However, the interpretation of this cross product is not straightforward: a ‘large’ cross product could correspond to a strong rortex barely suppressed by a weak $\boldsymbol {\nabla }\rho$, or to a weak rortex suppressed by a strong $\boldsymbol {\nabla }\rho$, which is the reason why normalised cross products were considered in figure 12.

5.3. Hypothesis for the role of rortices versus shear on mixing

To interpret the above observations on the different effects of shear versus rortex on the density field, we formulate in the following a hypothesis for their interaction and contributions to mixing. We first recall the flow visualisation in figure 2 showing strong local shear structures among the hairpin-like rortices that ‘straddle’ them. Being dominant in the vorticity field, the first contribution of the shear to mixing is to distort sharp density interfaces by shear instabilities, a process that forms vortical structures that we unequivocally identify as rortices. The role of these rortices (weaker relative to the shear) then appears to depend on their relative strength and morphology, and on those of the density gradient. When rortices are weak and density gradients are stronger, such as in the Holmboe regime, rortices tend to ‘scour’ density interfaces but do not overturn them (as pointed out by Salehipour, Caulfield & Peltier (Reference Salehipour, Caulfield and Peltier2016)), and thus result in little mixing. However, when rortices are strong and density gradients weaker, such as in turbulent regime, hairpin rortices within the shear layer (i.e. their legs) create bursting (i.e. lift-up or sweep-down events in the $z$ direction), thereby further stirring fluid within the partially-mixed layer (this is visible in the contours of velocity $v$ and $w$ in supplementary movies 2–5).

On the other hand, the most strongly rotating parts of the hairpins at the edges of the stratified layer (i.e. their heads) cause overturning and entrainment, thus broadening the mixing layer. Strong local shear then further stretches these newly created density gradients, accelerating small-scale molecular diffusion and ultimately achieving mixing.

In the next section we explore this hypothesis and consolidate our understanding of the subtle role of rortices.

6.1. In Holmboe waves

Starting with a snapshot of asymmetric Holmboe flow H4 in figure 13(a,b), we first observe in panel (a) that the streamwise vorticity ($\omega _x$, see colours) is mainly concentrated either under the two sides of the wave ‘head’ or on the two sides of the wave ‘body’ (shown by an isosurface of $S$ in grey). The ‘rortex line’ (black lines labelled A, B and C), equivalent to a streamline but based on the rortex vector $\boldsymbol {R}$, connects the regions of high opposite-signed $\omega _x$ in a $\Lambda$ shape. This indicates that the rortex we observed in § 3 likely originates from the ‘confined Homboe wave’ of Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a) (their paper was solely based on this dataset H4).

Figure 13. Rortex–density interaction in asymmetric Holmboe dataset H4. Comparison between experimental (a,b) and linear stability (c,d) results. (a) Streamwise vorticity ($\omega _x$, in colour) in three $y$–$z$ planes I, II and III, at time $t_n = 180$, superimposed on an isosurface of shear structure ($S=1.3$, in grey). (b) Spanwise component of $\boldsymbol R\times \boldsymbol {\nabla }\rho$ (in colour), superimposed with the contours of $\partial _z\rho$ (grey contour) in the $x$–$z$ mid-plane $y = 0$. Blue lines indicate the contour of $S=1.3$, and vectors indicate velocity in the plane. (c) Same visualisation as (a) but for the fastest-growing mode of the 3-D ‘confined Holmboe instability’ computed in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a), superimposed on the isosurface $S=2.5$. (d) Same visualisation as (b). Only the structures within the regions of $-0.8< y,z<0.8$ are shown, and the $z$ axis is stretched by a factor of 3 as in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a).

The $y$ component of $\boldsymbol R\times \boldsymbol {\nabla }\rho$ in $y$–$z$ planes (in colours) is plotted in (b). The strongest interaction is located near the density interface where $|\boldsymbol {\nabla }\rho |$ (in grey) is largest (see the regions labelled D, E and F), recalling that here $|\boldsymbol {\nabla }\rho |$ is an order of magnitude larger than $R$ and that the two vectors are nearly perpendicular. High values of $(\boldsymbol R\times \boldsymbol {\nabla }\rho )_y$ are also found on either side of the wave head (see the region labelled G), due to the high rorticity $R$. The velocity profiles within the $x$–$z$ plane reveal a more inflectional, and thus potentially unstable, region above the density interface, especially near the wave head (labelled k in the figure). We believe all these characteristics are important in asymmetric Holmboe waves.

In figure 13(c,d) we show similar visualisations but for the numerical solution corresponding to the fastest growing (or most unstable) linear mode computed on the two-dimensional experimental base flow in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a) (these data are available on the repository Lefauve et al. Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018b). The agreement between the observed ‘CHW’ (top row) and the numerically predicted ‘confined Holmboe instability’ (bottom row) is excellent. The only discrepancy lies in the absence of the strong interaction region that appears near the wave head in the linear solution (labelled G in c). We conclude that this particular feature around the head (observed in the experimental data but absent from the linear solution) is likely caused by nonlinearities. Conversely, most other details of the rortex/density dynamics discussed previously can be attributed to linear instability dynamics, which are significantly modulated by the spanwise confinement in the square duct geometry (an effect studied in Ducimetière et al. Reference Ducimetière, Gallaire, Lefauve and Caulfield2021).

6.2. In intermittent turbulence

A snapshot of I6 is shown in figure 14. In this intermittently turbulent flow, the life cycle (appearance and disappearance) of rortices is chaotic. Following the approach of ‘kernel’ studies, this snapshot was selected at a time when rortices are in a relatively complete form (before their breakdown). Figure 14(a) shows two hairpin rortices detected by the isosurface $R=0.6$ (with its head pointing up and the other pointing down, respectively labelled Ru and Rd) travelling in opposite directions. The rortex Ru is stronger than Rd, with the strength of their heads being $R\approx 2.2$ and $R\approx 0.65$, respectively. As time evolves (see supplementary movie 4), Ru is lifted higher up, its head finally reaches outside the shear layer while its legs are stretched through the upper layer. The overall inclination angle of Ru is ${\approx }30^{\circ }$, with its head inclined more steeply at ${\approx }60^{\circ }$ and its leg inclined less steeply at ${\approx }20^{\circ }$. The inclination of this particular rortex is larger than the peak value (the maximum of the WCA pdf), see the I6 data in figures 7 and 8. Meanwhile, the weaker rortex Rd is inclined as a whole at $13^{\circ }$ and its head does not lift up. The life of Rd is relatively shorter as it quickly breaks down into smaller eddies.

Figure 14. Rortex–density interaction in intermittent dataset I6 at time $t_n=36$. (a) Isosurface of $R = 0.6$ with colour denoting the $z$ position. (b) Spanwise component of $\boldsymbol R\times \boldsymbol {\nabla }\rho$ in a $y$–$z$ (labelled I in a). (c,d) Density and vertical component of the density gradient in the same $y$–$z$ plane and in (e, f) a $x$–$z$ plane (labelled II in a). The black contour in (c) is for $R = 0.6$; the vectors in (c,e) show the (subsampled) velocity field; green contours in (d, f) are overturned regions where $\partial _z \rho \ge 0$; the blue lines in ( f) show $S=1.5$. Only a subvolume (in $x$) is shown here for better visualisation.

Moving on to figure 14(b), we see four regions of strong interaction (labelled 1 to 4) evidenced by the $y$-component of $\boldsymbol {R}\times \boldsymbol {\nabla }\rho$ in the $y$–$z$ plane labelled I in (a). Again, the weaker rortex Rd seems to have a stronger interaction with the density gradient, just as in dataset H4 in the above section. This is due to a stronger density gradient collocated with Rd, as is clear in panels (c,d) (showing $\rho$ and $\partial _z \rho$ in the same $y$–$z$ plane).

In figure 14(c) we see that dense fluid (in red) under Rd is lifted up after being driven laterally outward by the rortex. On the contrary, light fluid is driven downward by the legs of Ru. The combined effect of these two rortices stirs fluid of different density (blue and red) around and make them meet in the middle region (see the white arrows in c) where fluid is well mixed (cyan). This region between the two pairs of legs is where density overturns are observed (labelled A and B in figure 14d). Note that the density field near Ru is better mixed than near Rd, i.e. the asymmetry between Ru and Rd carries over to the density field.

We now turn our attention to figure 14(e, f) showing the velocity vectors and density within the $x$–$z$ plane ‘II’ in (a). The strong head of Ru is visible and causes lift-up and sweep-down, just like the legs (see the white arrows). The overturns seen in (d) are visible between the two rortices (labelled C and D in f). We also note that the density interface (white dashed line) tilts towards Rd. Strong shearing structures (blue lines) are found either at the edge of the stratified shear layer (aligned with high density gradients, see the region labelled F) or near the centre of a strong rortex (see the region labelled E). This upper region of high $S$ allows us to infer a relationship between the tilted rortex and the tilted region of high shear. In the H4 snapshot (figure 13), the rortex was weak, likely created by localised shear around the Holmboe wave crest. By contrast, in this I6 snapshot, the rortex is further strengthened and stretched, and the lift up of its head also lifts up the high shear region between its two legs.

Based on the above observation, we can describe rortex Ru as a strong ‘stirrer’ of weakly stratified fluids, and rortex Rd as a weaker ‘revolving door’ entraining denser, more strongly stratified fluid into the mixing zone, and pulling pre-mixed fluid away from it. However this ‘revolving door’ remains weak compared to the stratification at the interface and cannot destroy it entirely (in addition, the nature of the exchange flow in SID ensures that unmixed fluid continually replaces mixed fluid, thereby sustaining such interfaces).

6.3. In fully developed turbulence

In figures 15 and 16 we select a representative snapshot in T3. Here, rortices are more complex, making it more difficult to inspect isolated structures. Figure 15(a) shows the isosurface of $R=0.6$ (grey region), $\rho$ (in colour) and $\boldsymbol {u}$ vectors in two $y$–$z$ planes (labelled as p1 and p2). The regions where rortices intersect these two planes are numbered 1 to 9. Rortices 1 and 2 are the two legs of a large hairpin, whose head has been partly truncated in this figure because it protrudes outside the shear layer $|z| \le 1$ within which our analysis is restricted. The strong rortices 1, 2 and 6 near the edges of the upper density interface move fluid laterally (around their rortex axis), thereby entraining lighter fluid (in blue) downward and neutral fluid (in green) upward. However, due to the buoyancy restoring force, the vertical flow is less vigorous than the spanwise flow (note the arrow length). This entrainment pattern agrees with our previous observation in I6 (figure 14c). The streamwise length scale of these rortices is $\approx$1–2, i.e. of the order of the shear layer depth, which is the energy-containing scale according to the turbulent kinetic energy spectra of this dataset T3 (see figure 5(c, f,i) of Lefauve & Linden Reference Lefauve and Linden2022b). However, the diameter of these rortices is $\approx$0.1–0.4, i.e. almost an order of magnitude smaller, placing them in the inertial range. This typical length/diameter ratio of 5–10 supports the schematics of figure 9, in which we had anticipated the length scale data, not given by our statistical analysis of the orientation of rortices. The non-dimensional Kolmogorov scale in these turbulent datasets was estimated as $\eta \approx 2\theta ^{-1/4}Re^{-3/4} \approx 0.02$ (Lefauve & Linden Reference Lefauve and Linden2022b), although it is subject to large uncertainties (see the discussion at the end of their § 5.1.3, which also considers a smaller Kolmogorov scale of order 0.006). The most likely estimate $\eta \approx 0.02$ would make these rortices of length ${\approx }50\eta \unicode{x2013}100\eta$ and diameter ${\approx }5\eta \unicode{x2013}20\eta$.

Figure 15. Rortex–density interaction in turbulent dataset T3 at time $t_n=41$. (a) Isosurface of $R$ = 0.6 (in grey), with two planes (p1 and p2) showing $\rho (y,z)$. (b) Vortex lines (based on $\boldsymbol {R}$) with colour denoting the strength of rortex, with the same two planes showing $\partial _z\rho (y,z)$. Green contours show density overturns; blue lines show $S=1.5$. Only a subvolume (in $x$) is shown here for clearer visualisation. Blue numbers show the length scale.

Figure 16. Further rortex–density interaction in turbulent dataset T3 at time $t_n=41$ (as in figure 15). (a–d) Contours in the $y = 0$ plane of $\rho,\partial _z\rho,S$ and the $x$-component $(\boldsymbol R\times \boldsymbol {\nabla }\rho )_x$, respectively. (e–h) Contours in $y$–$z$ planes of the $y$-component $(\boldsymbol R\times \boldsymbol {\nabla }\rho )_y$. The solid lines in (a–d) indicate the position of the planes (I, II, III and IV) in (e–h), whereas the dashed lines indicate the position of planes p1 and p2 in figure 15.

Inspecting now the density gradient in figure 15(b), we observe that rortex 1 acts again as the typical ‘revolving door’ described in I6, which allows for density transport across the relatively sharp upper density interface. Rortex lines, whose colour indicate the strength $R$, intersect the $y$–$z$ planes in regions 1, 2, 5 and 6. The upper density interface and the corresponding high shear region (with $S=1.5$, see blue contours) between rortices 1 and 2 is visibly distorted. Between the upper and lower density interfaces (in the mixing region), rortices 3, 4, 5 and 7 stir the fluid, but the overturned region ($ {\partial }_z\rho \ge 0$, green contours) is spotty due to the weak stratification $ {\partial }_z\rho \approx 0$.

Figure 16 offers visualisation of the same snapshot in complementary planes (the $x$–$z$ plane with $y=0$ in (a–d) and four different $y$–$z$ planes in e–h). Based on $\rho$ (a) and $\partial _z\rho$ (b) the upper density interface is more irregular and slightly less stratified compared with the lower density interface. Velocity vectors also show that vertical motions are somewhat suppressed near the lower interface. A strong rortex can deform and even break a nearby strong density gradient in a vertical ‘eruption’ process across the interface (see region D where the hairpin head of the rortices 1 and 2 appears in figure 15). Large overturning regions (filled in green, see regions A, B and C) are usually close to regions of high density gradients, sometimes forming a ‘sandwich’ configuration.

Comparing $\partial _z \rho$ and $S$ contours in figure 16(b,c) reveals a good correlation between them (see also supplementary movie 5). Although this paper mainly discusses rotational structures, it should be kept in mind that shear-driven instabilities are clearly important, and probably even dominant, in the process of turbulent production and mixing. The dominant contribution of shear to the vorticity was also reported for isotropic turbulence in Nagata et al. (Reference Nagata, Watanabe, Nagata and Da Silva2020) and Boukharfane et al. (Reference Boukharfane, Er-raiy, Alzaben and Parsani2021) and may be a generic feature of turbulence.

Finally, we study the interaction between the rortex and the density fields, based on the $x$ and $y$ components of $\boldsymbol R\times \boldsymbol {\nabla }\rho$ in figure 16(d) and (e–h), respectively. We observe that the strongest interactions occur near the upper and lower interfaces of the stratified layer. Peaks in the $x$ component of $\boldsymbol {R}\times \boldsymbol {\nabla } \rho$, denoted by $(\boldsymbol {R}\times \boldsymbol {\nabla } \rho )_x$, are usually centred at the concentration of a pair of opposite values of the $y$-component $(\boldsymbol {R}\times \boldsymbol {\nabla }\rho )_y$ which again suggests a hairpin structure (see the pairs of D1–D2, E1–E2, G1–G2, F1–F2 and H1–H2 in (e–h) and compare with the structures of (d) at nearby $x$ locations). However, in some datasets, the opposite pair have unequal magnitude (such as G1 and G2 in (g)); these rortices are then not complete hairpins, but rather cane-like or quasi-streamwise rortices, such as a rortex with a single leg connecting to a spanwise-oriented head (Adrian Reference Adrian2007). Further comparison between (b,d) and ( f) shows that the ‘sandwich’ region A is consistent with the regions of K and E1–E2 pair, which indicates that the overturning is closely related to adjacent rortices.

7.1. Origin of hairpin rortices

Lefauve et al. (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a) described the (asymmetric) CHW in dataset H4 (here visualised in figures 2(a) and 13). This long-lived 3-D wave is dominated by a shear structure having a wide ‘body’ and a narrower ‘head’, which are well predicted by the most unstable mode of a suitable linear stability analysis. The same study pointed out that the shear in the CHW is closely related to the coherent intensification of spanwise vorticity due to confinement by the lateral walls of the duct, which was studied systematically in Ducimetière et al. (Reference Ducimetière, Gallaire, Lefauve and Caulfield2021). This vorticity intensification is reminiscent of the wave warping process (Hama & Nutant Reference Hama and Nutant1963) and the vortex concentration in wall-bounded flows (Smith Reference Smith1984). Therefore, we speculate that the formation of rortices in SID flows is similar to the scenario of wave-induced vortex generation in a transitional boundary layer, and that $\Lambda$-rortices originate from 3-D CHWs.

A cartoon of the development of the CHW and its nearby rortices is shown in figure 17. First, the confined Holmboe instability (Lefauve et al. Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018a, CHI) arises due to a resonant interaction between a vorticity wave (on the edge of the shear layer) and a gravity wave (at the density interface). It then grows to a finite amplitude, intensifying and warping the background spanwise vorticity ($\omega _y$) until it nonlinearly saturates to the observed CHW characteristic shear structure (indicated by the dashed line in figure 17a). Subsequently, a 3-D state (with strong gradients of vorticity) likely develops according to Smyth & Peltier (Reference Smyth and Peltier1991) and Smyth (Reference Smyth2006), evidenced by the nucleation of streamwise vorticity ($\omega _x$) on either side of the wave (figure 13a,b). These sites are located under the wave ‘head’ (or below the cusp) and above the wave body (or around the ‘neck’). Finally, rorticity further concentrates at these two locations, forming a $\Lambda$-shaped rortex (see the embryonic rortex vectors in both the experimental and linear stability results of figure 13. The newly generated rortices are labelled A and B in figure 17(a), corresponding to the two nucleation sites labelled A and B in figure 13(a,c).

Figure 17. Schematic view of (a) the origin of a hairpin rortex from Holmboe wave and its evolution in the (b) intermittent regime and (c) turbulent regime. Dashed lines indicate shear structure, and green segmented tubes indicate rortices (having direction $\boldsymbol {R}$). The background colour indicates density stratification. Abbreviations are: LU, lift up; SD, sweep down; LM, lateral movement; R, rortex vector, HS, high shear.

This above process is analogous to the wave-induced $\Lambda$-vortex scenario described in boundary layers (Lee Reference Lee2000; Jiang et al. Reference Jiang, Lee, Chen, Smith and Linden2020a), in that the amplification of a 3-D wave (soliton-like coherent structure) is a key initiator of a vortex, indicating that such structures may be generic to shear-driven turbulence, with and without walls. The phenomenon of Holmboe-wave-induced rortices was also suggested by the numerical simulations of Smyth & Winters (Reference Smyth and Winters2003). Because of the restoring force in the stratified layer, the induced embryonic $\Lambda$-rortex does not significantly distort the sharp density interface. The strength of the rortex is about $10\,\%$ that of the shear, thus the flow remains shear-dominated, although we recall that shear dominates over rotation even in isotropic turbulence. As pointed out by Salehipour, Peltier & Caulfield (Reference Salehipour, Peltier and Caulfield2018) and Zhou et al. (Reference Zhou, Taylor, Caulfield and Linden2017), in strongly stratified flows such as in the H datasets ($Ri\approx 0.1\unicode{x2013}0.6$), self-organisation occurs through ‘scouring’ motions that keep the density interface sharp and robust.

With increasing $\theta$ and $Re$, the $\Lambda$-structure increases in amplitude as the rorticity is amplified by stretching, thus a more readily identifiable hairpin rortex appears (labelled C in figure 17b). Its head inclines more steeply, due to the competition between the shear-induced stretching of the mean profile and the self-induced velocity of the $\Lambda$-vortex (Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999). The rortex ejects stratified fluid away from the interface, and its elongated legs stretch into the body of the former Holmboe wave, creating transient bursts consisting of lift-up (LU) and sweep-down (SD) events, as shown in figure 17(b) by the green and blue arrows, respectively. This stirs fluid above or below their original positions, creating a net buoyancy flux and production of turbulent kinetic energy from the mean shear. At this stage, the interface becomes more unstable, overturns more frequently, becomes thicker (see region labelled E in the figure), and the Holmboe wave becomes shorter-lived. Owing to the presence of a strong rortex and of its induced bursting behaviours, the shear corresponding to the CHW is intensified and becomes more unstable. This localised high shear further stimulates a new rortex (e.g. the embryonic rortex labelled D in the figure), which jointly contributes to create more intense intermittent fluctuations.

Although the model shown in figure 17 is inspired from asymmetric Holmboe data, we believe the rortex generation mechanism is similar in symmetric Holmboe data. In symmetric Holmboe waves, the typical $\Lambda$-rortex is more streamwise (as depicted in figure 9b), and the streamwise stretching of the rortex legs plays a significant role during the transition to turbulence. The more horizontal the hairpin rortex is, the easier it is for it to stir fluid up and down between its legs, generating vertical motions.

The source of the mixing in these flows is unequivocally ‘internal’ (following the classification of Turner Reference Turner1973) in the sense that the rortex develops internally within the shear layer, which contrasts with externally introduced vortices, such as the vortex rings produced by actuating a pump (Olsthoorn & Dalziel Reference Olsthoorn and Dalziel2015).

The above hypothesis on hairpin rortices originating from pre-turbulent CHWs differs from other mechanisms proposed in slightly different flows. For example, hairpin-like vortices may develop from shear-layer-stretched spanwise tubes formed during the breakdown of Kelvin–Helmholtz billows (Pham & Sarkar Reference Pham and Sarkar2011; Pham et al. Reference Pham, Sarkar and Winters2012), or from spanwise vortices hypothesised in stratified plane Poiseuille flow to be due to local quasi-linear processes and internal waves (Lloyd et al. Reference Lloyd, Dorrell and Caulfield2022).

7.2. Role of hairpin heads and legs on stirring

At higher $\theta$ and $Re$, there is a stronger rortex–density interaction, more bursting and overturning because: (i) the rortices are stronger and more horizontal (i.e. they have a smaller inclination angle to the $x$–$y$ plane, see figure 9); (ii) a weaker stratification ($Ri_b \approx 0.1\unicode{x2013}0.2$) is less able to suppress these vigorous vertical motions.

Figure 17(c) shows a schematic model of rortex–density interface interaction in such conditions representative of the late intermittent regime and the turbulent regime. Both quasi-streamwise and hairpin-like rortices are observed in the thicker, partially mixed layer bounded by two interfaces (see, e.g., figure 2c,d). Here, for simplicity, we only sketch one hairpin rortex on the upper interface.

The hairpin rortex consisting of a head and two legs is inclined to the true horizontal plane at a small angle (5$^{\circ} \sim 15^{\circ }$). The legs with opposite rotation lift up denser fluid in the middle of the rortex (see green arrow LU1 in figure 17c), whereas on the outer side of the rortex legs, lighter fluid is swept down (see blue arrow SD1). This bursting cycle also enhances a lateral exchange of momentum along the legs (see the black arrows labelled LM1 and LM2), which produces sufficiently large spanwise stress for a burst regeneration (Landahl Reference Landahl1975; Haidari & Smith Reference Haidari and Smith1994). Compared to vertical motions, such lateral motions in stratified turbulence are more energetic due to the absence of buoyancy restoring force in the spanwise direction.

The appearance of a distinct rortex head (or transverse rortex, along $y$), as a manifestation of the strong concentration of spanwise rorticity, provides another bursting cycle, labelled LU2 and SD2. Dense fluid is lifted up (LU2) by the head, and overturned by the corresponding sweep-down (SD2). Importantly, the lift-ups caused by the rortex legs and head creates a localised, inclined high-shear region slightly above the rortex (dashed line, labelled HS). Secondary instabilities of this localised high-shear produce further rortices which enhance the mixing, in a fractal-like fashion.

In figure 18, we further represent rortices according to their position: (i) on the interface of the stratified layer or (ii) within the partially mixed layer, away from the interfaces. The green longitudinal cylinder indicates a longitudinal rortex (labelled RL, e.g. streamwise vortex, legs of hairpin vortex or cane-like vortex), whereas the magenta spiral indicates a transverse rortex (labelled RT, e.g. spanwise vortex, hairpin head). As explained in § 5, the rortex across the interface acts as a ‘revolving door’ lifting inner (pre-mixed) denser fluid away from the upper interface and entraining outer lighter fluid into the shear layer. Both RL and RT can be candidates for this role, as shown in figure 18(a). The localised breakup of the interface is closely related to this bursting cycle (see rortex 1, 2 and 3 for example). However, the rortices within the stratified shear layer (and away from the interface, e.g. rortex 4 and 5) contribute little to this revolving door, being instead primarily responsible for further lateral stirring of the pre-mixed fluid, causing the spanwise density gradient to be further smoothed, making the interior better mixed. For rortices across the interfaces, this lateral rotation deforms the shear layer and promotes the inner-outer exchange of fluid (see rortex 6 in the cross-sectional view of figure 18b).

Figure 18. Cartoon of turbulent rortices near a density interface and within the partially mixed layer. (a) Side view ($x$–$z$ plane). (b) Cross-sectional view ($y$–$z$ plane). RL: longitudinal rortex (green) with $\boldsymbol {R}$ primarily along $x$; RT: transverse rortex (magenta) with $\boldsymbol {R}$ primarily along $y$.