The evolution of coherent vortical structures in increasingly turbulent stratified shear layers

Abstract

We study the morphology of Eulerian vortical structures and their interaction with density interfaces in increasingly turbulent stably stratified shear layers. We analyse the three-dimensional, simultaneous velocity and density fields obtained in the stratified inclined duct laboratory experiment (SID). We track, across 15 datasets, the evolution of coherent structures from pre-turbulent Holmboe waves, through intermittent turbulence, to full turbulence and mixing. We use the rortex–shear decomposition of the local vorticity vectors into a rortex vector capturing rigid-body rotation and a shear vector. We describe the morphology of ubiquitous hairpin-like vortical structures (revealed by the rortex), similar to those commonly observed in boundary-layer turbulence. These are born as relatively weak vortices around the strong three-dimensional shearing structures of confined Holmboe waves, and gradually strengthen and deform under increasing turbulence, transforming into pairs of upward- and downward-pointing hairpins propagating in opposite directions on the top and bottom edge of the shear layer. The pair of legs for each hairpin are counter-rotating and entrain fluid laterally and vertically, whereas their arched-up ‘heads’, which are transverse vortices, entrain fluid vertically. We then elucidate how this large-scale vortex morphology stirs and mixes the density field. Essentially, vortices located at the sharp density interface on either edge of the mixing layer (mostly hairpin heads) engulf blobs of unmixed fluid into the mixing layer, whereas vortices inside the mixing layer (mostly hairpin legs) further stir it, generating strong, small-scale shear, enhancing mixing. These findings provide new insights into the role of turbulent coherent structures in shear-driven stratified mixing.

1. Introduction

Stably stratified turbulence and the enhanced mixing across density isosurfaces (isopycnals) that it accomplishes is a crucial but poorly understood component of many deep-ocean and coastal flow systems of importance under a changing climate. A ‘grand challenge’ of environmental fluid dynamics is to parameterise accurately this small-scale ‘diapycnal’ mixing in large-scale circulation models to improve predictions for the vertical transport of heat, carbon dioxide, salt and other scalars in our oceans (Dauxois et al. 2021). To complement expensive and sparse field observations, laboratory experiments have a key role to play in the effort to develop better mixing parameterisations. In this paper we use datasets obtained from such an experiment, the stratified inclined duct (SID), which sustains a two-layer exchange flow in an inclined square duct. This experiment allows us to accurately control the flow geometry and levels of interfacial turbulence by a systematic variation of two key non-dimensional flow parameters. Using newly available measurement technologies (Partridge, Lefauve & Dalziel 2019), this experiment also allows us to obtain the three-dimensional Eulerian velocity and density fields simultaneously at high spatio-temporal resolutions, and thus to study three-dimensional coherent structures like never before.

These coherent flow structures, and especially vortical structures, exist across a wide spectrum of spatio-temporal scales and play an important role in the processes of turbulent bursting and mixing. Previous studies identified a range of vortical structures in stratified shear layers (i.e. in a nearly parallel layer of vorticity not caused by a solid wall, and that embeds a density interface), such as streamwise or quasi-streamwise vortices (Schowalter, Van Atta & Lasheras 1994; Caulfield & Peltier 2000), spanwise vortices (Salehipour, Peltier & Mashayek 2015) and hairpin vortices (Smyth & Winters 2003; Watanabe et al. 2019). In the shear-driven flows of interest in this paper, it is meaningful to define vortical structures after an appropriate treatment of the ‘contaminating’ shear (Shrestha et al. 2021). To distinguish rigid-body rotation from straining motions, various decomposition based on the velocity gradient tensor have been proposed (Kolář 2007; Li, Zhang & He 2014; Gao & Liu 2018; Keylock 2018; Nagata et al. 2020; Watanabe, Tanaka & Nagata 2020; Hayashi, Watanabe & Nagata 2021). Here, we apply to our state-of-the-art experimental datasets the new rortex–shear (RS) decomposition proposed by Liu et al. (2018) and Xu et al. (2019) to decompose the three-dimensional (3-D) vorticity vector into a rigid-body rotation vector, the ‘rortex’ vector $\boldsymbol{R}$, and a shear vector $\boldsymbol{S}$. In some literature (e.g. Shrestha et al. 2021) $\boldsymbol{R}$ is referred to as the ‘liutex’ vector, however we use the term ‘rortex’ (as originally proposed by Liu et al. (2018) to emphasise its rotational nature. The 15 datasets and accompanying movies are available online (Jiang et al. 2022).

Although vortices can be produced externally (e.g. by artificial vortex rings impinging a density interface as a model for turbulent eddies, see Linden 1973; Olsthoorn & Dalziel 2015), they naturally develop internally, either from internal gravity waves (Fritts, Arendt & Andreassen 1998) or from shear-driven instabilities leading to (usually short-lived) Kelvin–Helmholtz billows (Caulfield & Peltier 2000) or (usually long-lived) Holmboe waves (Smyth & Winters 2003). Lefauve et al. (2018a) described and explained the origin of ‘confined’ Holmboe waves in the SID experiment, a typical example of long-lived coherent vortical structures, which they visualised using a simple vorticity threshold. The 3-D development of the Holmboe-wave instability was studied numerically by Smyth & Winters (2003), who noted that ‘Loop structures in the density field associated with hairpin-like vortices are a conspicuous feature of turbulent Holmboe waves. These structures are initiated by secondary instabilities (in one case this resembled the localised convective instability described by Smyth & Peltier 1991) and grow to large amplitude via vortex stretching’. The hypothesis of horseshoe (or hairpin) vortices was initially proposed by Theodorsen (1952), and has proven key to the understanding of boundary-layer turbulence (Acarlar & Smith 1987; Smith et al. 1991; Adrian 2007; Jiang 2019; Lee & Jiang 2019). Head & Bandyopadhyay (1981) used smoke visualisation to investigate the evolution of hairpins in a turbulent boundary layer with increasing Reynolds number, and they found that elongated hairpin vortices were inclined at a characteristic angle of approximately $40\unicode{x2013}50^{\circ }$ to the wall. These hairpins were observed to be less elongated and more isolated at low Reynolds number, and to agglomerate and become very elongated at high Reynolds number. The inclination and evolution (generation and regeneration) of hairpins was subsequently studied in more detail, numerically by Zhou et al. (1999) and experimentally by Haidari & Smith (1994). A hypothesis based on soliton-like coherent structures has been put forward to explain the bursting process and the generation of hairpins in wall-bounded flows (Lee 1998; Lee & Wu 2008; Jiang et al. 2020a,b). Hairpin-like structures have also been observed in stably stratified boundary layers, experimentally by Williams (2014) and numerically by Atoufi, Scott & Waite (2019); they are apparently similar to those found in unstratified boundary layers.

In stably stratified shear layers (not visibly influenced by top and/or bottom walls), such as deep ocean overflows, exchange flows through straits or saltwater intrusions in estuaries, Kelvin–Helmholtz or Holmboe instabilities (found in weakly and strongly stratified flows, respectively) can grow in a symmetric or asymmetric fashion depending on the vertical offset between the centres of the velocity profile and the density profile (Carpenter, Lawrence & Smyth 2007). Hairpin vortices have been associated with these wave trains, especially after they succumb to secondary instabilities (i.e. further symmetry breaking in the third dimension), whose breakdown creates fully 3-D turbulence (Smyth & Moum 2000; Smyth 2006; Pham, Sarkar & Winters 2012). Recently, using direct numerical simulations, Watanabe et al. (2019) found inclined hairpin vortices throughout the stratified shear layer, and argued that turbulent mixing was very active at the length scales close the streamwise extent of the hairpins. In stratified plane Poiseuille flow, Lloyd, Dorrell & Caulfield (2022) found numerically that hairpin vortices arise far from the wall and interact with strong buoyancy gradient, inducing internal wave breaking.

However, despite tantalising numerical evidence of their existence and role in stratified shear-driven mixing, hairpin vortices have until now not been described in comparable laboratory flows. It also remains unclear: (i) how they develop from pre-turbulent flows (especially Holmboe waves) and evolve in increasingly turbulent flows; and (ii) how they interact with density interfaces and participate in density overturning, stirring and ultimately mixing. These are the two questions that we address in this paper.

In § 2 we introduce our experimental datasets and explain their relevance to our objectives. In § 3 we visualise vortical structures first by a traditional method, and then by our new method based on the RS decomposition of vorticity, to build intuition for the subsequent statistical analyses. In § 4 we reveal the detailed morphology of the ‘rortices’ identified by the rortex (and, to a lesser extent, of the shear) by a ‘weighted conditional averaging’ (WCA) method. In § 5 we study the interaction between rortices and density gradients, and in § 6 we confirm and augment these statistical results with ‘case studies’ on instantaneous snapshots. In § 7 we synthesise these results and propose a tentative model for the origin and role of hairpin vortices in stratified shear layers, and we conclude in § 8.

2. Experimental datasets

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. (2018a) (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 2014; Lefauve, Partridge & Linden 2019; Lefauve & Linden 2020).

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. (2019), 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. (2019), 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. (2018a, 2019) and Partridge et al. (2019), in this paper we use the slightly modified datasets of Lefauve & Linden (2022a,b) (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 2022c) .

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. (2019) 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. (2019). 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. (2018a). 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 2022b). 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. (2019), Appendix A, and Lefauve & Linden (2022b), 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. Identification of vortical structures

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. 1999; Gao, Ortiz-Dueñas & Longmire 2011). 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 1988), $\varDelta$-criterion (Chong, Perry & Cantwell 1990), $\lambda _2$-criterion (Jeong & Hussain 1995) and $\lambda _{ci}$-criterion (Zhou et al. 1999; Chakraborty, Balachandar & Adrian 2005). 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 2018; Shrestha et al. 2021), 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. (2018) and Liu et al. (2018) 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 (2018), followed by its explicit expression with a significantly improved implementation in Xu et al. (2019), 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 1967), 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. 2018) 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. (2018a) 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 (2022a) and Lefauve & Linden (2022b) 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. (2018) and Gao & Liu (2018).

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. (2018a), 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. Detailed morphology of rortices

The magnitudes of $\boldsymbol {R}$ and $\boldsymbol {S}$ in the previous section provided us with the general trends of their spatial structures. This section tackles their more detailed morphology, and in particular the 3-D orientation of hairpin rortices, revealed by a comprehensive statistical analysis of our $\boldsymbol {R}$ datasets.

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. 2018a). 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. 2019), two-point spatial correlations of velocity (Williams 2014), geometrical features of some isosurfaces/contours based on conventional vortex diagnostics (Zhou et al. 1999) and vortex vectors (Arun et al. 2019). According to Arun et al. (2019), 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.3.2. With respect to the $z'_\perp$ plane (true horizontal)

In figure 7(b), we observe the following.

1. (i) All $\beta _R$ pdfs are also nearly statistically symmetric around $0^{\circ }$, even at the highest tilt angles $\theta =5^{\circ }$ and $6^{\circ }$, confirming a certain symmetry of rortices with respect to ${{z'}_\perp }$ (the true horizontal plane) rather than ${{z}_\perp }$ (the $x$-$y$ plane based on the duct coordinate system).

2. (ii) Increasingly strong rortices are inclined at increasingly steep angles to the true horizontal plane, i.e. the peak $\beta _R$ moves away from $0^{\circ }$ (see the dashed trend lines in the figure), especially at low ${Re}$ and high $\theta$ (e.g. H2, H4); but this tendency diminishes somewhat (i.e. the trend line becomes more vertical) in turbulent flows (e.g. T1–T3) or at low $\theta$ and low ${Re}$ (e.g. H1, H3).

3. (iii) Both weak and strong rortices have a relatively narrow peak $\beta _R=0^{\circ }$ (of width ${\approx }\pm 10^{\circ }$) in the bottom-left corner of the $(\theta,Re)$ plane (H1, H3, I1). This peak widens slightly to ${\approx }\pm 20^{\circ }$ for the strong rortices in the top-right corner (I8, T1, T2, T3).

4. (iv) The other pdfs (top-left and bottom-right corners of the $(\theta,Re)$ plane) tend to be more uniform or bimodal, in particular H2, H4 which have two clear peaks at $\beta _R\approx \pm (35\sim 55)^{\circ }$.

The evolution of the peak values for $\beta _R$ (and $\beta _{R'}$ in open $\vartriangleright$) are plotted in figure 8(b). Holmboe flows are again split into two categories. In asymmetric H2 and H4 (and to some extent I1), $\Lambda$-rortices (without elongated tails) are inclined to the horizontal plane at an angle about ${\pm }35\sim 55^{\circ }$, with their head inclined more steeply. In symmetric H1 and H3, rortices lie close to the horizontal plane at $0\sim \pm 10^{\circ }$. Supplementary movies 2–5 show the complete time evolution of these structures. In turbulent datasets, the inclination angle becomes smaller with increasing $\theta {Re}$ (see inclined dashed line). However, there is a slight difference between pdfs of $\beta _{R}$ and $\beta _{R'}$ (open and closed symbols) in that the mean shear seems to suppress the lift-up of rortices (i.e. $\beta _R<\beta _R'$). In H2 and H4, the central peak (the maximum of the weighted pdf in figure 7b) $\beta _R=0^{\circ }$ observed for weak rortices ($k_{th}<1$) seem less due to Holmboe waves than to the mean background shear $\partial _z \bar {u}, \partial _y \bar {u}$, because this peak disappears entirely when considering the pdf $\beta _{R'}$ (based on the fluctuation $\boldsymbol {R}'$, see figure 8b). However, the mean shear does not appear to ‘contaminate’ the spanwise inclination of rortices ($\alpha _R\approx \alpha _R'$) in the more turbulent datasets (see figure 8a).

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. 1999). 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. (2019) (who reported an inclination of 14$^{\circ }$) and those hairpin legs observed in unstratified wall-bounded flows in Haidari & Smith (1994) and Zhou et al. (1999) (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 1981; Zhou et al. 1999). 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. Interaction with density gradients

The statistical study in § 4 showed that the strength and orientation of the rortices change in the different flow regimes (from Holmboe to turbulence). The morphology and behaviour of rortices are expected to be influenced by the density (or buoyancy) field in stably stratified flows having relatively large bulk Richardson numbers between $Ri\approx 0.15\unicode{x2013}0.55$ (in H flows) and $Ri\approx 0.1\unicode{x2013}0.2$ (in I and T flows). In this section we examine the interaction between the rortex $\boldsymbol {R}$ and density gradient $\boldsymbol {\nabla } \rho$. We first apply our WCA method to the 3-D density gradient, before studying the averaged strength of interaction between rortex, shear and density gradients along the vertical. Finally, we discuss the complex relationship between vortical structures and mixing.

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 (2016)), 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. Case studies: instantaneous snapshots

The above hypothesis is based on the statistics of 15 datasets using spatial and temporal averaging, which reveals the general characteristics of structures. To verify that these characteristics are indeed representative of the actual flow phenomenology, we now study instantaneous volumetric snapshots of the rortex–density dynamics. We focus on relatively isolated rortices, inspired by the approach taken in ‘kernel’ studies of turbulent boundary layers for the interaction of vortical structures and the generation of near-wall bursts (Haidari & Smith 1994). We study three datasets in turn: H4, I6 and T3, which cover the three key flow regimes of interest.

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. (2018a) (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. (2018a), 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. (2018a).

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. (2018a) (these data are available on the repository Lefauve et al. 2018b). 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. 2021).

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 2022b). 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 2022b), 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. (2020) and Boukharfane et al. (2021) 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 2007). 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. Synthesis and discussion

Synthesising our previous statistics and case studies on rotational structures and density interfaces, we now propose a simplified model for the evolution of their morphology under increasing turbulence.

7.1. Origin of hairpin rortices

Lefauve et al. (2018a) 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. (2021). This vorticity intensification is reminiscent of the wave warping process (Hama & Nutant 1963) and the vortex concentration in wall-bounded flows (Smith 1984). 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. 2018a, 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 (1991) and Smyth (2006), 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 2000; Jiang et al. 2020a), 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 (2003). 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 (2018) and Zhou et al. (2017), 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. 1999). 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 1973) 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 2015).

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 2011; Pham et al. 2012), or from spanwise vortices hypothesised in stratified plane Poiseuille flow to be due to local quasi-linear processes and internal waves (Lloyd et al. 2022).

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 1975; Haidari & Smith 1994). 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$.

8. Conclusions

In this paper, we have investigated the morphology of coherent vortical structures, and more specifically ‘rortices’, and their relation to density gradients. We have adopted an empirical (data-driven) approach based on the analysis of 15 experimental datasets of increasingly turbulent stratified shear layers obtained by exchange flow in a long inclined square duct.

Using the standard $Q$-criterion, we first observed (figure 1) that coherent vortical structures are mainly hairpin-like (e.g. $\Lambda$-structures, hairpin structures, cane structures, quasi-streamwise structures) irrespective of the flow regime. After splitting vorticity into rigid-body rotation (rortex vector $\boldsymbol {R}$) and anti-symmetric shearing motion (shear vector $\boldsymbol {S}$) in figure 2, we examined their averaged magnitude $R, S$ and vertical distribution in all datasets (figure 3). We have found that the shear $S$ always dominates, although the rorticity $R$ increases significantly in the intermittent and turbulent regimes, and shows coherent structures reminiscent of other turbulent shear flows.

We then studied the morphology of vector-based coherent structures, rortices, about which the flow is in rigid-body rotation, using detailed statistics (WCA) of the orientation of the rortex $\boldsymbol {R}$ (figures 4–8). This allowed us to draw in a key schematic (figure 9) the evolution of typical rortices from asymmetric (high-$\theta$ low-$Re$) and symmetric (low-$\theta$ high-$Re$) Holmboe waves (H regime), to intermittent (or transitional) flow (I regime) and eventually to turbulent flows (T regime). The two types of Holmboe waves have different rortices, the high-$Re$ symmetric ones being more similar to those found in the turbulent regime. In the intermittent (transitional) regime, some datasets show similarities to the asymmetric H regime, whereas others show similarities to the symmetric H or T regimes. Strong transverse rortices (which we attribute to ‘wide’ hairpin heads) are most clearly observed in the I regime.

Applying a similar statistical analysis to the density gradients and to the cross products of $\boldsymbol {R}$ or $\boldsymbol {S}$ with $\boldsymbol {\nabla }\rho$ (figures 10–12), we found that increasingly turbulent density interfaces were increasingly steeply inclined (with respect to their true horizontal equilibrium) as a result of weaker stratification and increasing interaction with rotational structures (rortices) across the entire shear layer. By contrast, strong interaction between density gradients and shear only occurs on the edges of the partially mixed region. We also found that while rortices can be generated baroclinically by strong density gradients, the strongest rortices were not associated with strong density gradients; in contrast, they appear in the partially mixed region presumably as a result of shear-driven instabilities.

To complement and validate our insights based on time- and volume-averaged statistics, we examined in figure 13–16 the instantaneous rortex and density interface morphologies from a representative snapshot in three datasets representative of each regime. In the H regime (figure 13), the region above the density interface and near the wave head is the most unstable. The position of initial streamwise vorticity concentration agrees with the position of embryonic rortices due to a supposed secondary, nonlinear instability (which are not present in the corresponding structure predicted by a linear stability analysis). In the I regime (figure 14) the region between a pair of hairpin rortices (pointing up and down) is the most overturned and mixed. Strong shearing structures are located either at edge of the shear layer and aligned with region of high density gradient, or near the centre of a strong rortex. In the T regime, overturnings are frequently sandwiched by two high-density-gradient regions near the upper and lower density interfaces and flanked by nearby rortices (figures 15 and 16).

To synthesise these statistical and structural results, we proposed a model for the evolution of vortical structures and density interfaces (figure 17). First, we hypothesise that $\Lambda$-rortices originate from the 3-D CHW, described in Lefauve et al. (2018a), through the formation of a highly localised shear region, the nucleation of secondary instabilities and longitudinal roll-up. In turbulent flows, this rortex is strengthened and arches up, creating characteristic hairpin rortices that stir the fluid around multiple axes. Their effect on stirring (ultimately leading to mixing) is explained in figure 18. Rortices present across the upper or lower density interface act as revolving doors that drive (mixed) fluid away from the interface and entrain outer (unmixed) fluid into the mixing region. Both longitudinal rortices (e.g. hairpin legs) and transverse rortices (e.g. hairpin head) are candidates for this role. However, rortices present within the mixing region are mainly responsible for further stirring the pre-mixed fluids by lift-up/sweep-down events as well as by strong lateral movement.

These large-scale rotational stirring motions explain the generation of a vertical buoyancy flux and the production of turbulent kinetic energy from the mean shear, which are both key energy fluxes in stratified turbulence. At smaller scales, high strain rates and density gradients cause turbulent kinetic energy dissipation and irreversible mixing, respectively. However, our analysis in this paper has been purely kinematic; the dynamical role of coherent rortices versus shear in shaping the turbulent energetics, and in particular the efficiency of mixing, remains an open question.