## 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.* Reference Dauxois2021). 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 Reference Partridge, Lefauve and Dalziel2019), 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 Reference Schowalter, Van Atta and Lasheras1994; Caulfield & Peltier Reference Caulfield and Peltier2000), spanwise vortices (Salehipour, Peltier & Mashayek Reference Salehipour, Peltier and Mashayek2015) and hairpin vortices (Smyth & Winters Reference Smyth and Winters2003; Watanabe *et al.* Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019). 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.* Reference Shrestha, Nottage, Yu, Alvarez and Liu2021). To distinguish rigid-body rotation from straining motions, various decomposition based on the velocity gradient tensor have been proposed (Kolář Reference Kolář2007; Li, Zhang & He Reference Li, Zhang and He2014; Gao & Liu Reference Gao and Liu2018; Keylock Reference Keylock2018; Nagata *et al.* Reference Nagata, Watanabe, Nagata and Da Silva2020; Watanabe, Tanaka & Nagata Reference Watanabe, Tanaka and Nagata2020; Hayashi, Watanabe & Nagata Reference Hayashi, Watanabe and Nagata2021). Here, we apply to our state-of-the-art experimental datasets the new rortex–shear (RS) decomposition proposed by Liu *et al.* (Reference Liu, Gao, Tian and Dong2018) and Xu *et al.* (Reference Xu, Gao, Deng, Liu and Liu2019) 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.* Reference Shrestha, Nottage, Yu, Alvarez and Liu2021) $\boldsymbol{R}$ is referred to as the ‘liutex’ vector, however we use the term ‘rortex’ (as originally proposed by Liu *et al.* (Reference Liu, Gao, Tian and Dong2018) to emphasise its rotational nature. The 15 datasets and accompanying movies are available online (Jiang *et al.* Reference Jiang, Lefauve, Dalziel and Linden2022).

Although vortices can be produced *externally* (e.g. by artificial vortex rings impinging a density interface as a model for turbulent eddies, see Linden Reference Linden1973; Olsthoorn & Dalziel Reference Olsthoorn and Dalziel2015), they naturally develop *internally*, either from internal gravity waves (Fritts, Arendt & Andreassen Reference Fritts, Arendt and Andreassen1998) or from shear-driven instabilities leading to (usually short-lived) Kelvin–Helmholtz billows (Caulfield & Peltier Reference Caulfield and Peltier2000) or (usually long-lived) Holmboe waves (Smyth & Winters Reference Smyth and Winters2003). Lefauve *et al.* (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018*a*) 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 (Reference Smyth and Winters2003), 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 Reference Smyth and Peltier1991) and grow to large amplitude via vortex stretching’. The hypothesis of horseshoe (or hairpin) vortices was initially proposed by Theodorsen (Reference Theodorsen1952), and has proven key to the understanding of boundary-layer turbulence (Acarlar & Smith Reference Acarlar and Smith1987; Smith *et al.* Reference Smith, Walker, Haidari and Sobrun1991; Adrian Reference Adrian2007; Jiang Reference Jiang2019; Lee & Jiang Reference Lee and Jiang2019). Head & Bandyopadhyay (Reference Head and Bandyopadhyay1981) 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.* (Reference Zhou, Adrian, Balachandar and Kendall1999) and experimentally by Haidari & Smith (Reference Haidari and Smith1994). 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 Reference Lee1998; Lee & Wu Reference Lee and Wu2008; Jiang *et al.* Reference Jiang, Lee, Chen, Smith and Linden2020*a*,Reference Jiang, Lee, Smith, Chen and Linden*b*). Hairpin-like structures have also been observed in stably stratified boundary layers, experimentally by Williams (Reference Williams2014) and numerically by Atoufi, Scott & Waite (Reference Atoufi, Scott and Waite2019); 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 Reference Carpenter, Lawrence and Smyth2007). 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 Reference Smyth and Moum2000; Smyth Reference Smyth2006; Pham, Sarkar & Winters Reference Pham, Sarkar and Winters2012). Recently, using direct numerical simulations, Watanabe *et al.* (Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019) 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 (Reference Lloyd, Dorrell and Caulfield2022) 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.* (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018*a*) (see their § 3). The SID sustains a hydraulically controlled exchange flow inside a long duct (of length $L=1350$ mm) of square cross-section (of height and width $H=45$ mm). The duct is inclined at a small angle $\theta$ with respect to the horizontal and connects two large reservoirs initially filled with aqueous salt solutions of different densities $\rho _0 \pm \Delta \rho /2$. The Prandtl number is thus $Pr=\nu /\kappa \approx 700$, where $\nu$ and $\kappa$ are the kinematic viscosity and salt diffusivities, respectively.

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

### 2.2. Measurements and processing

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

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

Instead of the original datasets used and visualised in Lefauve *et al.* (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018*a*, Reference Lefauve, Partridge and Linden2019) and Partridge *et al.* (Reference Partridge, Lefauve and Dalziel2019), in this paper we use the slightly modified datasets of Lefauve & Linden (Reference Lefauve and Linden2022*a*,Reference Lefauve and Linden*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 Reference Lefauve and Linden2022*c*) .

### 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

and

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

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

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

Although subject to inherent technical limitations (summarised in Lefauve *et al.* (Reference Lefauve, Partridge and Linden2019), Appendix A, and Lefauve & Linden (Reference Lefauve and Linden2022*b*), 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.* Reference Zhou, Adrian, Balachandar and Kendall1999; Gao, Ortiz-Dueñas & Longmire Reference Gao, Ortiz-Dueñas and Longmire2011). For example, $\boldsymbol {\omega }$ is not generally aligned with the local rotation, and the maximum of $|\boldsymbol {\omega }|$ does not generally coincide with the ‘core’ of a vortex, because vorticity does not discriminate between shear and swirl (rotation).

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

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

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

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

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

### 3.2. RS decomposition

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

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

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

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

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

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

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

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

### 3.3. Rortex and shear structures

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

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

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

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

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

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

Figure 3(*b*) shows that the symmetric Holmboe flows (datasets H1 and H3) have their peak $\langle S\rangle _{xyt}$ at $z>0$; these flows feature two trains of counter-propagating waves due to the velocity interface $\langle u \rangle _{xyt}=0$ and density interface $\langle \rho \rangle _{xyt}=0$ being nearly coincident around $z\approx 0\unicode{x2013}0.15$ (see LL22a figure 3). By contrast, the asymmetric Holmboe flows (H2 and H4) have their peak $\langle S\rangle _{xyt}$ at $z<0$; these flows feature a single train of waves due to the velocity interface $\langle u \rangle _{xyt}=0$ being slightly offset from the density interface $\langle \rho \rangle _{xyt}=0$, the latter of which is around $z \approx -0.2$ (see LL22a figure 3). This peak is due to the ‘body’ of the CHW structure described in Lefauve *et al.* (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018*a*), 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)$:

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

*a*–

*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.

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

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$.

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

### 4.2. Inclination of shear structures

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

### 4.3. Inclination of rortex structures

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

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

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

In figure 7(*a*) we observe the following.

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

(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.

(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}$.

(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}$.

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.

(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).

(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).

(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).

(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 7*b*) $\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 8*b*). 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 8*a*).

### 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$).

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

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

## 5. 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}$).

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).

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 3*b*).

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

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

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

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

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

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

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

## 6. 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 Reference Haidari and Smith1994). 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.* (Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018*a*) (their paper was solely based on this dataset H4).

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

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

### 6.2. In intermittent turbulence

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

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 14*d*). 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 14*c*). The streamwise length scale of these rortices is $\approx$1–2, i.e. of the order of the shear layer depth, which is the energy-containing scale according to the turbulent kinetic energy spectra of this dataset T3 (see figure 5(*c*, *f*,*i*) of Lefauve & Linden Reference Lefauve and Linden2022*b*). However, the diameter of these rortices is $\approx$0.1–0.4, i.e. almost an order of magnitude smaller, placing them in the inertial range. This typical length/diameter ratio of 5–10 supports the schematics of figure 9, in which we had anticipated the length scale data, not given by our statistical analysis of the orientation of rortices. The non-dimensional Kolmogorov scale in these turbulent datasets was estimated as $\eta \approx 2\theta ^{-1/4}Re^{-3/4} \approx 0.02$ (Lefauve & Linden Reference Lefauve and Linden2022*b*), 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$.

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

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

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

Finally, we study the interaction between the rortex and the density fields, based on the $x$ and $y$ components of $\boldsymbol R\times \boldsymbol {\nabla }\rho$ in figure 16(*d*) and (*e*–*h*), respectively. We observe that the strongest interactions occur near the upper and lower interfaces of the stratified layer. Peaks i