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

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. 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 vorticity field into a pure rotational part (the rortex vector), and a non-rotational part (the 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. Each hairpin's pair of legs are counter-rotating and entrain fluid laterally and vertically, and 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, while 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.


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 will 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 , 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). Among them are: streamwise or quasistreamwise vortices (Schowalter et al. 1994;Caulfield & Peltier 2000), hairpin vortices (Smyth & Winters 2003;Watanabe et al. 2019), spanwise vortices (Salehipour et al. 2015), columnar vortices (Billant & Chomaz 2000;Waite & Smolarkiewicz 2008), and pancake vortices (Fincham et al. 1996;Riley & Lelong 2000;Riley & DeBruynKops 2003). Most of these were identified by Eulerian criteria, either based on a threshold of vorticity (Basak & Sarkar 2006;Mashayek & Peltier 2012) or through the eigenvalues of the velocity gradient tensor (Hunt et al. 1988;Chong et al. 1990;Jeong & Hussain 1995;Zhou et al. 1999;Chakraborty et al. 2005). In the shear-driven flows of interest in this paper, meaningful turbulent vortical structures must be defined after an appropriate treatment of the 'contaminating' mean shear (Shrestha et al. 2021). Recent efforts have been devoted to decompose the velocity gradient tensor into a rotational and shearing part (Li et al. 2014;Keylock 2018;Nagata et al. 2020;Watanabe et al. 2020;Hayashi et al. 2021). Here, we apply to our state-of-the-art experimental datasets the new Rortex-Shear (RS) decomposition proposed by Wang et al. (2019) and Xu et al. (2019) to decompose the three-dimensional (3-D) vorticity field into a 'pure rotation' field and a 'pure shear' field.
Although vortices (i.e. coherent regions of pure rotation) 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 et al. 1998) or, more typically, from sheardriven 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 turbulent boundary layer with increasing Reynolds number, and they found that elongated hairpin vortices were inclined at a characteristic angle of approximately 40 − 50 • 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 details, 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 et al. (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 stronglystratified 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 et al. 2007). Hairpin vortices have been associated with those 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 et al. 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 et al. (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 their important 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 preturbulent 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 will 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 Rortex-Shear 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' method. In § 5 we study the interaction between rortices and density gradients. In § 6 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 § 7.

Setup and flows
The datasets analysed in this paper were collected in the Stratified Inclined Duct (SID) experiment, whose setup 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 θ with respect to the horizontal and connects two large reservoirs initially filled with aqueous salt solutions of different densities ρ 0 ± ∆ρ/2. The Prandtl number is thus P r = ν/κ ≈ 700, where ν and κ are the average kinematic viscosity and salt diffusivities, respectively.
Increasing θ (defined to be positive when it accelerates the flow) and/or the Reynolds number Re ∝ g(∆ρ/ρ 0 )HH/ν (primarily set by the density difference) 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 (θ, Re) plane and their transitions have been studied extensively (Meyer & Linden 2014;Lefauve et al. 2019;Lefauve & Linden 2020).

Measurements and processing
We consider 15 datasets, each corresponding to a single experiment performed at a given θ 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 ρ 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 g along the 'true vertical' is thus tilted with respect to this coordinate system and has components [g sin θ, 0, −g cos θ] along (x, y, z). See Lefauve et al. (2019), figure 1 for a schematic and figures 3-4 for a snapshot of u and ρ 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, ρ data obtained at spanwise locations y = y i (i = 1, 2, . . . , 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 δ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 × 45 × 45 mm 3 ) and obtain high x − z planar resolutions for both velocity and density. Each experiment typically captures ≈ 300 volumes (time snapshots), and each volume typically contains ≈ 400 × 40 × 80 velocity vectors in x, y, z, respectively.
Instead of the original datasets used and visualised in Lefauve et al. (2018aLefauve et al. ( , 2019; 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 nondimensionalised the coordinates and flow variables of each individual dataset using (i) half the density difference ∆ρ/2 (after removing the mean ρ 0 ) such that −1 ρ 1; (ii) half the actual resulting shear-layer depth, such that −1 z 1; and (iii) half the actual (mean) peak-to-peak velocity magnitude, such that −1 u 1 with the mean velocity extrema u x,y (y = 0, z = ±1) = ∓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 in this paper), can all be freely downloaded from their repository Lefauve & Linden (2022c) .

Parameters and resolution
The corresponding 'shear-layer' Reynolds number Re and bulk Richardson number Ri b are defined as in LL22a § 3.2-3.3 as and (2.2) These parameters are consistent with the non-dimensionalisation introduced above. These parameters consist of (i) a 'hydraulics' part based on input parameters, including half the peak-to-peak dimensional velocity scale ∆U/2 ≡ √ g H = g(∆ρ/ρ 0 )H, half the duct height H/2 and ν, and (ii) a 'shear-layer' rescaling based on half the nondimensional output (measured after the hydraulic non-dimensionalisation) peak-to-peak velocity magnitude δu/2 (typically ≈ 0.5 − 1.2) and shear layer depth h/2 (typically ≈ 0.5−0.7). Note that Re and Ri b were denoted as Re s and Ri s b , respectively, in LL22a,b to emphasise this specific shear-layer non-dimensionalisation.
The parameters of datasets H1-T3 are shown 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 θRe. The historical names of datasets (H1, . . . H4, I1, . . . I8, T1, . . ., T3) were based on increasing values of the product of θ with the hydraulics (input) Reynolds number ∆U H/(4ν) in Lefauve et al. (2019) rather than on the shear-layer (output) Reynolds number Re used in this paper. Nevertheless, datasets remain approximately ordered with increasing θRe.
As a rule of thumb, the vector resolution of the data approaches the Kolmogorov turbulent lengthscale (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. The full vector resolution of each dataset in x, y, z, t are given in LL22a, Appendix A, and LL22b Appendix B), we will show below that these 15 state-of-the-art datasets deliver new insights into the time-resolved, 3-D coherent structures of sheardriven stratified turbulence.

Previous methods and Q-criterion
Since 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 ω = ∇×v is subjective and generally inappropriate (Zhou et al. 1999;Gao et al. 2011). For example, ω is not generally aligned with the local rotation, and the maximum of |ω| does not generally coincide with the 'core' of a vortex, since vorticity does not discriminate between shear and swirl (rotation).
Several improved Eulerian vortex identification schemes based on the eigenvalues of the velocity gradient tensor ∇v have thus been developed since the 1980s, including the Q-criterion (Hunt et al. 1988), ∆-criterion (Chong et al. 1990), λ 2 -criterion (Jeong & Hussain 1995), λ ci -criterion (Zhou et al. 1999;Chakraborty et al. 2005), etc. 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 ∇v calculated as Q = ||B|| 2 − ||A|| 2 , where A = 1 2 (∇v + (∇v) T ) and B = 1 2 (∇v − (∇v) T ) are, respectively, the symmetric (or strain rate tensor) and anti-symmetric (or rotation rate tensor) components of the ∇v. We plot a single snapshot of Q-criterion vortices identified by the isosurface Q > 0 (i.e. the local rotation exceeds strain) in a selection of 10 representative datasets having three different tilt angles: θ = 5 • in panels (a,b,d,g,j ) (datasets H2, H4, I5, I8, T3, with increasing Re), θ = 6 • in panels (c,e,i ) (datasets I4, I6, T2) and θ = 3 • in panels (f,h) (datasets I7, T1). The colouring of the iso-surfaces denotes the vertical location z and the legend on each plot identify the dataset, the tilt angle, the Reynolds number and the Q value of the iso-surface displayed.
Broadly speaking, from the Holmboe to the turbulent regime (i.e. with increasing θRe), vortical structures evolve from individual, disconnected hairpins which start as Λ-structures without elongated trailing legs (panels a,b), to groups of hairpins with elongated legs (panels h,i,j ). The hairpins of low-Re flows are relatively weak and less obvious (see panels c,d ), especially when a corrugation appears on the iso-surfaces (which we attribute to a low signal-to-noise ratio), or when the head of the hairpin is broken. However, we will show below that these hairpin-shaped vortices are indeed robust features of Holmboe waves. In higher-Re flows, hairpins are stronger and more obvious (higher  (see table 1) and the iso-surface Q value are listed on the right of each panel. Colours on the iso-surfaces denotes the z position (we show −0.9 < z < 0.9, i.e. the middle 90% of the shear layer). signal-to-noise ratio), with their head tending toward the edges of the shear layer (|z| ≈ 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. Panels (f,g,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 panels h) as an extreme form of asymmetric hairpin vortices. Panels (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 shows the complete time evolution of panels (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.

Rortex-Shear decomposition
The scalar Q-criterion (i.e. the isosurface Q = Q 0 , where Q 0 is the threshold) is inevitably subjective to the (somewhat arbitrary) value of the threshold. Furthermore, Q-vortices do not represent rigid rotation since they are contaminated by shear .
Recently, a new 'vortex' vector called the 'rortex' (or 'liutex') was proposed by Liu et al. (2018) and  to isolate the rigid rotational part from the shear, and directly point in the direction of local fluid rotation. The vorticity is decomposed uniquely as ω = R+S, where R is the rortex vector, and S is the shear vector (indicating a non-rotational, anti-symmetric part of ∇v). This Rortex-Shear (RS) decomposition thus separates pure rotational structures from shearing structures. This decomposition is particularly helpful in shear-driven turbulence, as in this paper. Furthermore, Xu et al. (2019) showed that the RS decomposition was related to a special (transposed) Schur form of the velocity gradient tensor (rotating the ∇v into a suitable block-triangular matrix). This relationship allowed the use of efficient library functions to speed up the RS decomposition by avoiding the tedious original algorithms requiring coordinate transformations.
According to Xu et al. (2019) and Wang et al. (2019), the rortex could then be explicitly calculated as The direction of the rortex r is the local unit real eigenvector of ∇v, indicating the rotational axis, and λ ci is the imaginary part of the complex conjugate part of eigenvalues of ∇v. This is the practical definition that we apply to our datasets in the remainder of this paper. It is based on the idea that ∇v can have either one or three real eigenvalues.
When there is only one real eigenvalue, the direction of the rortex R is aligned with the associated normalised eigenvector r selected such that ω · r > 0. The conjugate pair of complex eigenvalues have imaginary parts ±λ ci characterising the rotation about r.
When there are three real eigenvalues, λ ci = 0 and thus R = 0, meaning that all the vorticity is shear without rotation. Equation 3.1 shows the relative importance of λ ci vs the vorticity projected onto r.
Hereafter, we denote the magnitude of the shear vector |S| = S, and we shall refer to the magnitude of the rortex vector |R| = R as the rorticity and the vortical structure it identifies simply as a rortex. We shall also use the term vortex to refer more The value of the iso-surfaces are 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, so has been the lower half z < 0 in (c-d ).
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. Figure 2 shows an iso-surface of both rortex 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 (crosssectional 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 iso-surface and its projection on the respective slices. (See the 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 pure shearing structure (blue isosurface), which is highly reminiscent of the shape of the confined Holmboe wave described in Lefauve et al. (2018a) from iso-surfaces of the spanwise component of vorticity (ω y ) of the same dataset. This similarity is because the shear is very strong 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 ten times larger than the red isosurface R = 0.13. This indicates that shear accounts for a substantial part of the vorticity according to ω = R + S. Although weaker, rortices are also observed near the 'head' of the Holmboe wave, in a Λ-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 smallerscale structures, such that 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 vs shear also increases (see the colour bars and isosurface values in the caption), indicating an increasing importance of rortices in turbulent mixing.

Averaged distribution of rorticity and shear
The relative strengths of R and S are explored quantitatively in figure 3. In panel (a) we plot, for all 15 experimental datasets, the averaged magnitudes R xyzt and S xyzt (where · xyzt represents the average over the whole shear layer volume and time as in LL22a,b) against the product θRe (our proxy for increasing levels of turbulence). Both R xyzt (♦) and S xyzt ( ) are nearly constant at ≈ 0.1 − 0.2 and ≈ 1, respectively, when θRe 80 (where θ 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 v = v −v where· = · t is the spatially varying temporal average, giving the underlying parallel shear flowū(x, y, z) (sincev,w ≈ 0). Figure 3(a) also shows the resulting fluctuating shear S ( ) and rortex R ( ). First, we find that R ≈ R (the symbols are nearly indistinguishable) whereas S S, as expected in the presence of background shear ∂ zū . Second, we find that S is only about two to three times larger than R (a weaker dominance compared to that of S over R), and S seems to increase fairly linearly with θRe even before the transition to fully turbulent flow at θRe ≈ 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 .
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 (panel b) and R (panel c) 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 are displayed as transparent shadings underlying each curve. Figure 3(b) shows that the symmetric Holmboe flows (datasets H1 and H3) have their peak S xyt at z > 0; these flows feature two trains of counter-propagating waves due to the velocity interface u xyt = 0 and density interface ρ xyt = 0 being nearly coincident around z ≈ 0 − 0.15 (see LL22a figure 3). By contrast, the asymmetric Holmboe flows Figure 3. (a) Distribution of R xyzt and S xyzt for all datasets, separating H, I, and T data by the solid grey lines ( for S; ♦ for R; for S ; for R ; the dashed lines are trend lines of the distribution; 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 distribution of S xyt and R xyt along z direction for all datasets, respectively, separating H, I, and T data in different rows. 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 the S and R (i.e. S xy − S xyt and similarly for R).
(H2 and H4) have their peak S xyt at z < 0; these flows feature a single train of waves due to the velocity interface u xyt = 0 being slightly offset from the density interface ρ xyt = 0, the latter of which is around z ≈ −0.2 (see LL22a figure 3). This peak is due to the 'body' of the confined Holmboe wave 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 'plateau' arrow), presumably due to the 'head' of the confined Holmboe wave. With increasing levels of turbulence (I regime, middle row), this plateau in S 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 = ±1) of the turbulent shear layer (T regime, bottom row). Their values ≈ 1.5 − 2 is comparable to those in the H regime. A further interesting observation is that the temporal root mean square fluctuations of S 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 R 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 ≈ 400 vs 200). In the 'late' I regime (e.g. I7, I8) and in the T regime, a broad peak in R xyt is centered 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 to the I regime), possibly caused by the breakdown and interaction of rortices. This suggests that the emergence and increasing importance of the rortex is a fundamental aspect of turbulence dynamics and mixing, justifying our greater focus on R (vortical structures) than S (shear structures) in the remainder of this paper.

Detailed morphology of rortices
The magnitudes of R and 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 R datasets. 4.1. Weighted conditional averaging (WCA) and orientation probability distribution functions (pdf 's) In order to remove noise and to give stronger weight to stronger rortices, we first apply a 'conditional sampling' method on the fields R(x, t) (containing between 0.3 − 1 × 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 (x, t): where k th is a tunable threshold level below which the data are discarded, and R rms is the root mean square 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 'orientiation probability density functions' (pdf's) to quantify the likelihood of the orientation of specific vectors, measured by their frequency distribution in our datasets. For any non-zero vector V (where V may represent R, S, etc) we define the angles between V and the planes x ⊥ and z ⊥ as: 180] is the angle in degrees between a and b.
The unit vectors are defined as follows:x,ẑ are the unit vectors along the streamwise (x) and wall-normal (z) direction of the duct; x ⊥ indicates the plane normal tox;ẑ = cos θẑ − sin θx is the 'true vertical' unit vector (in the opposite direction of gravity); and z ⊥ is the 'true horizontal' plane normal to theẑ . 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 pdf's, 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 (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 α or β at a value of i ∈ [−90 • , 90 • ] and for a given conditional threshold of k th is calculated by The 'true horizontal' plane z ⊥ (in green in b) is normal to the opposite direction of gravityẑ , while the plane x ⊥ (in green in a) is normal tox. Finally, θ is the tilt of the duct with respect to the horizontal direction (the convention is that θ > 0 indicates that the flow is forced).   where n i is count of the occurrences when the angle under consideration belongs to the interval (bin) i ± δ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 90 −90 N di = R 2 xyzt , i.e. the area under the curve of N gives the time-and volumeaveraged 'rorstrophy' satisfying the threshold R/R rms k th . If k th = 0, the original rortex field and all existing rortices are considered, following (4.1). The above-described weighted (by R 2 ) and conditional (by selecting only R/R rms k th ) averaging will be applied to angle frequency distributions and used to study how progressively stronger rortices are aligned with respect to x ⊥ and z ⊥ . By analogy, we will also extend our weighted conditional averaging (hereafter WCA) to S (weighing by S 2 and conditioning by S/S rms k th ).
Before showing our results on the orientation of R, S, it is worthwhile 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 f xyzt ∈ [0, 1] (the time and volume-averaged ratio of points satisfying R > k th ) decreases with increasing threshold level k th . The semi-log axes and the exponential fit (dashed line) reveal that f t decreases approximately exponentially with k th with decay constant ≈ 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 ∇v has a pure rotation component). All 15 datasets display a similar behaviour at small k th , but their curves spread out significantly for k th 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 v ) is indistinguishable from that of R.
Next, figure 5(b) shows the vertical distribution of volume fraction f 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 R xyt (z) in figure 3(c) without conditional thresholds. In particular T2 and T3 have a robust ≈ 20 % fraction for |z| 0.5, which tapers off to only ≈ 5 % at the edges |z| = 1, justifying our cropping of the original datasets (see § 2.2) to restrict our attention to the |z| 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 broadly 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 θ, low Re) and symmetric Holmboe waves (found at low θ, 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-θ, low-Re turbulence is characterised by more overturning motions and less extreme shear-dominated enstrophy than low-θ, high-Re turbulence.

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 S, remembering its background role (in fact, dominant in magnitude) in shear-driven, stratified turbulence. Figure 6 shows the pdf's of angles α S and β 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 ∂ zū = O(1)), and weight S 2 (giving emphasis to large shear events). We see that large shear events are generally perpendicular to bothx andẑ and thus primarily alonĝ y, indicating the dominance of the spanwise component of vorticity, which motivated the use of ω y in our previous study (Lefauve et al. 2018a). Increasing turbulence (blue to red curves) widens the pdf of both α S and β S , revealing increasingly 3-D shearing structures. These results are robust for other thresholds k th .

Inclination of rortex structures
We now turn to a comprehensive analysis of the orientation pdf's of the vortex vector R angles, which show different and more subtle behaviours than those of the shear. We plot in figure 7 the pdf's of α R (panel a) and β R (panel b). Each of the 15 datasets are plotted in separate subpanels, and arranged according to their position in the (θ, Re) plane, in order to draw connections between the observed vortex statistics and the two key flow parameters. Furthermore, the pdf's are weighted with R 2 , and we use curves of increasingly dark colour (blue in panel a and red in panel 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 sub-panels, we use the same axis limits in all sub-panels and normalise all pdf's such that the integral of each over the interval −90 • to 90 • gives the all-time and volume-averaged square norm of the conditioned R 2 xyzt (rather than 1, as in figure 6). The dashed grey 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 (∂ z ρ > 0), and the best fit in (θ, Re) space was shown to scale with θ 3.17 Re 1.75 .

With respect to the x ⊥ plane
In figure 7(a) we observe the following: (i) All α R pdf's are statistically symmetric around 0 • 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 α R = 0 • appears at higher k th , identified by blue arrows.
(iii) A peak angle (maximum of the weighted pdf) at α R ≈ ±60 • appears in the bottomleft and top-right corners of the (θ, Re) plane, i.e. either at low θ and low Re or high θ and high Re.
(vi) In contrast to (iii), a wider and more uniform pdf across α R ∈ [-60∼60] appears in the upper-left and lower-right corners, i.e. either low θ and high Re s or at large θ and low Re s .
These features are clearly captured in figure 8(a), which shows the evolution of the peak values of α 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 R (withoutū), which are essentially similar to the full rortex data R.
We conclude that both Re and θ 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 α R ≈ ±60 • to the x ⊥ plane, largely independent of θ (though it varies across a wider range than Re does). However, at lower Re values (H, I regimes), the evolution of α R with Re (at fixed values of θ, i.e. along horizontal lines) varies depending on the value of θ. For example, at θ = 5 • , the unimodal or uniform pdf becomes bimodal at higher Re (compare the evolution H2 → H4 → I5 → I8 → T3), whereas almost the opposite happens at θ = 2 (evolution H3 → I1→ I2 → I3). This difference again indicates that Holmboe waves at low vs 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 α R = 0 • (perpendicular to the x-axis) suggests transverse rortices or the heads of hairpin rortices, which become less dominant under stronger turbulence.

With respect to the z ⊥ plane (true horizontal)
In figure 7(b), we observe the following: (i) All β R pdf's are also nearly statistically symmetric around 0 • , even at the highest tilt angles θ = 5 • and 6 • , confirming a certain symmetry of rortices with respect to z ⊥ (the true horizontal plane) rather than z ⊥ (the x-y plane based on the duct coordinate system); (ii) Increasingly strong rortices are inclined at increasingly steep angles to the x − y plane, i.e. the peak β R moves away from 0 • (see the dashed trend lines in the figure), especially at low Re and high θ (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 θ and low Re (e.g. H1, H3).
(iv) The other pdf's (top-left and bottom-right corners of the (θ, Re) plane) tend to be more uniform or bimodal, in particular H2, H4 which have two clear peaks at β R ≈ ±(35 ∼ 55) • . The evolution of the peak values for β R (and β R in open ) are plotted in figure 8(b). Holmboe flows are again split in two categories. In asymmetric H2 and H4 (and to some extent I1), Λ-rortices (without elongated tails) are inclined to the horizontal plane at an angle about ±35∼55 • , with their head inclined more steeply. In symmetric H1 and H3, rortices lie close to the horizontal plane at 0 ∼ ±10 • . Supplementary movie 2 shows the complete time evolution of these structures. In turbulent datasets, the inclination angle becomes smaller with increasing θRe (see inclined dashed line). However, there is a slight difference between pdf's of β R and β R (open and closed symbols) in that the mean shear seems to suppress the lift-up of rortices (i.e. β R < β R ). In H2 and H4, the central peak (the maximum of the weighted pdf in figure 7b) β R = 0 • observed for weak rortices (k th < 1) seem less due to Holmboe waves than to the mean background shear ∂ zū , ∂ yū (which are not 'pure shear'), since this peak disappears entirely when considering the pdf β R (based on the fluctuation R , see figure 8b). However, the mean shear does not appear to 'contaminate' the spanwise inclination of rortices (α R ≈ α R ) in the more turbulent datasets (see figure 8a).

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 (panels a-d ), while side views (in the x − z plane, where x is normal to z ) are shown in the bottom row (panels 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 α R of the rortices in panels (b-d ) is close to that of hairpin vortices observed in unstratified turbulent boundary layers having α R ≈ 40 ∼ 60 • (Zhou et al. 1999). For the vertical inclination, the rortices in the T regime (panel h) have an average |β R | ≈ 11 • 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 • ) and those hairpin legs observed in unstratified wall-bounded flows in Haidari & Smith (1994); Zhou et al. (1999) (who reported ≈ 8 • ). However, these angles are lower than the average tilt of hairpin head in unstratified turbulent boundary layers, which are typically around 45 • (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 'liftup' of the head is inhibited by stratification. The legs of the symmetric Holmboe rortices (panel f ) are almost horizontal, showing a much clearer connection with intermittent and turbulent structures than the asymmetric Holmboe ones (panel e).
These schematics are ideal models of symmetric Λ and hairpin rortices representative of the evolution of global α R and β 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, as suggested by figure 1 (though based on the Q criterion).

Interaction with density gradients
Rortices are inevitably influenced by the density (or buoyancy) field in stably-stratified flows having bulk Richardson numbers between Ri ≈ 0.15 − 0.55 (in H flows) and Ri ≈ 0.1 − 0.2 (in I and T flows). In this section we examine the interaction between the rortex R and density gradients ∇ρ vectors. We first apply our weighted conditional averaging (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, three examples, or case studies, are discussed to illustrate aspects of the complex relationship between vortical structures and mixing.

Distribution of the density gradient and relation to the rortex
In figure 10 we plot the WCA pdf's of the angle between ∇ρ and the 'true horizontal' plane z ⊥ , i.e. β ρ = ∠(∇ρ,ẑ ) − 90. We compare the pdf's under two different weights: in panel (a) with |∇ρ| 2 to highlight the strongest density gradients (we refer to this as β ρ1 ); and in panel (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 pdf as β ρ2 ).
In panel (a) we find, unsurprisingly for stably-stratified flows, that the strongest density Figure 10. Frequency distribution (orientation pdf's) of the WCA vertical angles of the density gradient: (a) βρ1 weighted by |∇ρ| 2 without conditional threshold; (b) βρ2 weighted by R 2 and using a threshold of R > k th = 2. (c) Peak angle (maximum of the weighted pdf) βρ1 (empty ) and βρ2 (filled ). We also add the horizontal angles αρ1 with weight |∇ρ| 2 (empty ) and αρ2 with weight R 2 (filled ). Note the log vertical scale. Symbol size indicates the relative strength of the peak, i.e. the value of its ordinate in (a,b).
gradients overwhelmingly point downwards, thus β ρ1 > 0 • , with a small deviation from the perfect 'true' vertical of 90 − β ρ1 ≈ 5 • . By contrast, in panel (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 β ρ2 < 0 • ). 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 − β ρ2 ≈ 7 − 10 • in the H datasets to more substantial deviations of ≈ 16 − 25 • 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 pdf's decay much more slowly, with significant overturns (we recall that the overturn fraction was shown by light grey contours in figure 7, giving a typical ≈ 3 − 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 β ρ2 (colour-filled ), noting in passing the difference with the evolution (or lack thereof) of β ρ1 (empty ). The trend of a monotonically decreasing β ρ2 (or increasing 90 − β ρ2 ) with θRe is very clear. We also add the peak in horizontal angle α ρ = ∠(∇ρ,x) − 90 (α ρ1 weighted by |∇ρ| 2 and α ρ2 weighted by R 2 ) to show that they both have a consistent peak at 0 • .
Interpreting these results, we note that although both shear and rorticity increase with turbulence (as was shown in figure 3), the decrease in β ρ2 is presumably caused by the increasing dominance of nearly horizontal rortices having β R ≈ 90 • (see figure 8). These Figure 11. Peak angle φp between R and ∇ρ 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 and grey shading have weight R 2 ; white-filled and green shading have weight |∇ρ| 2 . All statistics weighted by R 2 have conditional threshold k th = 2.
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 'preturbulent', especially H1 (marked by a red arrow) which peaks at a 70 • . However, the asymmetric Holmboe data (with low-Re, H2 and H4, labelled with blue arrows) are again different, since the rortex remains relatively weak compared to 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 β ρ2 peaks at 30 ∼ 40 • (labelled by blue arrows), which is very close to the peak β R in figure 8(b).
To further study the relationship between rortices and density interfaces, we define a last angle φ = ∠(R, ∇ρ). The evolution of its peak angle φ p , extracted from its WCA distribution, is shown in figure 11. Symbol sizes denote the peak height, while the shading denote the width (or spread) of the distribution at half height (see top right insert). Again we compare distributions under two different weights: |∇ρ| 2 (white circles and green shading) and R 2 (red symbols and grey shading).
First, all peak angles φ p ≈ 90 • , revealing a new piece of information: rortices and density gradients are most frequently perpendicular, in all datasets, regardless of statistical weight.
Second, distributions of φ weighted by |∇ρ| 2 (green shading) are narrower than those weighted by R 2 (grey shading), having a typically spread of ±10 • vs ±30 • (in some datasets even higher, e.g. H2, H4, I4). This observation brings an important nuance to our above conclusion: while the strongest ∇ρ are indeed frequently perpendicular to R, the strongest R are less frequently perpendicular to ∇ρ. This subtle asymmetry in the relation between R and ∇ρ 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 ∇ρ (green shading) is relatively 'narrow', it broadens somewhat from approximately 90 ± 5 • (H data) to 90 ± 15 • (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 ≈ 0.15 in T data) and thus of a weaker feedback of density in the momentum equation.

Vertical profiles of their interaction
We now turn to the strength of |∇ρ| and of its interaction with S and R along the z direction. In figure 12(a-c) (top row) we show |∇ρ| xyt (z), segregating the H, I, and T data in different columns. We see that the initially sharp density interface broadens (from H → I → T) and ends up in panel (c) becoming partially mixed across −0.5 z 0.5, flanked by two weaker interfaces. This evolution is similar to that of the shear S xyt (z) (previously shown in figure 3b).
In figure 12(d-f ) (middle row) we plot the averaged magnitude of their cross product |S × ∇ρ| xyt scaled by |∇ρ| xyt . The profile is very close to S xyt (shown by dashed lines), indicating that sin (∠(S, ∇ρ)) ≈ 1, i.e. that S (primarily along y) is approximately perpendicular to ∇ρ (primarily along z). Figure 12(g-i ) (bottom row) shows |R × ∇ρ| scaled by |∇ρ| xyt . The interaction of rortex and density gradient is distributed more evenly across z, with weak peaks that neither reflect the peak of R xyt (in dashed lines) nor the peak of |∇ρ| 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 find any region with strongly peaked |R × ∇ρ| across the shear layer.

Hypothesis for the role of rortices and shear on mixing
To interpret the above observations on the different roles of shear vs rortex on the density field, we formulate below 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 depends 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 can hardly destroy them, resulting 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, 3, 4 and 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 will explore this hypothesis and consolidate our understanding of the subtle role of rortices on mixing.

Case studies: instantaneous snapshots of the rortex-density interactions
The above hypothesis is based on the statistics of 15 data sets using spatial and temporal averaging, which reveals 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 shall study three datasets in turn: H4, I6 and T3, which cover the three key flow regimes of interest.

H4 case
Starting with a snapshot of asymmetric Holmboe flow H4 in figure 13(a-b), we first observe in panel (a) that the streamwise vorticity (ω 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 R, connects the regions of high opposite ω x in a Λ 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).
To show the interaction between the coherent rortex and the density gradient, we plot the y component of R × ∇ρ in y-z planes (in colours) in panel (b). The strongest interaction is located near the density interface where |∇ρ| (in grey) is largest (see the regions labelled D, E and F), recalling that here |∇ρ| is an order of magnitude larger than R and that the two vectors are nearly perpendicular. High values of (R × ∇ρ) y are also found on either sides 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 twodimensional experimental base flow in Lefauve et al. (2018a) (these data are available on the repository Lefauve et al. (2018b)). The agreement between the observed 'confined Holmboe wave' (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 panel 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 above can be attributed to purely linear instability dynamics, which are significantly modulated by the spanwise confinement in the square duct geometry (for more details, see Ducimetière et al. (2021)).

I6 case
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). Panel (a) shows two hairpin rortices detected by the iso-surface R = 0.6 (with its head pointing up and the other 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 ≈ 2.2 and R ≈ 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 ≈ 30 • , with its head inclined more steeply at ≈ 60 • and its leg inclined less steeply at ≈ 20 • . The inclination of this particular rortex is larger than the peak value (the maximum of the WCA pdf), see the I6 data in figure 7 and 8. Meanwhile, the weaker rortex Rd is inclined as a whole at 13 • and its head does not lift up. The life of Rd is relatively shorter as it quickly breaks down into small eddies.
Moving on to 14(b), we see four regions of strong interaction (labelled 1 to 4) evidenced by the y-component of R × ∇ρ in the y − z plane labelled I in panel (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 ρ and ∂ z ρ 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 panel 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 panel 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 panel (d ) are visible between the two rortices (labelled C and D in panel f ). We also note that the density interface (white dashed line) tilts towards Rd. Strong shearing structures (blue lines) are found either at 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 (besides, the nature of the exchange flow in the SID experiment ensures that unmixed fluid continually replaces mixed fluid, thereby sustaining such interfaces).

T3 case
In figures 15-16 we select a representative snapshot in T3. In this turbulent flow, rortices are more complex, making it more difficult to inspect isolated structures. Figure 15(a) shows the iso-surface of R = 0.6 (grey region), ρ (in colour), and 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 since it protrudes outsides the shear layer |z| 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 (notice the arrow length). This entrainment pattern agrees with our previous observation in I6 (figure 14c).
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, interact and combine, thereby increasing mixing, but the overturned region (∂ z ρ 0, green contours) is rather spotty due to the weak stratification ∂ z ρ ≈ 0. Figure 16 offers visualisation of the same snapshot in complementary planes (the x − z plane with y = 0 in panels a-d and four different y − z planes in panels e-h). Based on ρ (panel a) and ∂ z ρ (panel b) the upper density interface is more irregular and slightly less stratified compared to the lower one. 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 ∂ z ρ and S contours in figure 16(b,c) reveals a good correlation between them (see also the 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.
Finally, we study the interaction between the rortex and the density fields, based on the x and y components of R × ∇ρ 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 R×∇ρ, denoted by (R×∇ρ) x , are usually centred at the concentration of a pair of opposite values of the y-component (R × ∇ρ) y which again suggests a hairpin structure (see the pairs of D1-D2, E1-E2, G1-G2, F1-F2 and H1-H2 in panels e-h and compare to the structures of panel d at nearby x locations). However, in some datasets, the opposite pair have unequal magnitude (such as G1 and G2 in panel 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 panels (b,d ) and panel (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.

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.

Origin of hairpin rortices
Lefauve et al. (2018a) described the (asymmetric) confined Holmboe wave (CHW) in dataset H4 (here visualised in figures 2a 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 Λ-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 (ω 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 (ω 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 Λ-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).
This above process is analogous to the wave-induced Λ-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 Λ-rortex does not significantly distort the sharp density interface. The strength of the rortex is about 10 % that of the shear, so the flow is still clearly shear-dominated. As pointed out by Salehipour et al. (2018) and Zhou et al. (2017), in strongly-stratified flows such as in the H datasets (Ri ≈ 0.1 − 0.6), selforganisation occurs through "scouring" motions to keep the density interface sharp and robust.
With increasing θ and Re, the Λ-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 Λ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. Due 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 Λ-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 (1979)) 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 confined Holmboe waves differs from other mechanisms proposed in slightly different flows, such as the spanwise instabilities of Kelvin-Helmholtz rollers in unstable shear layers (Pham & Sarkar 2011;Pham et al. 2012) or the internal waves and quasi-linear processes in stratified plane Poiseuille flow (Lloyd et al. 2022).

Role of hairpin rortices in turbulent mixing
At higher θ 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 ≈ 0.1 − 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 • to 15 • ). The legs with opposite rotation lift up denser fluid in the middle of the rortex (see green arrow LU1 in figure 17c), while 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), while the magenta spiral indicates transverse rortex (labelled RT, e.g. spanwise vortex, hairpin head). As explained in § 5.4, 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 LV and TV 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).

Conclusions
In this paper, we investigated the morphology of coherent vortical structures, and more specifically 'rortices', and their relation to density gradients. We adopted an empirical (data-driven) approach based on the analysis of 15 state-of-the-art 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. Λ-structures, hairpin structures, cane structures, quasi-streamwise structures) irrespective of the flow regime. After splitting vorticity into pure rotational part (rortex vector R) and non-rotational part (shear vector S) in figure 2, we examined their averaged magnitude R, S and vertical distribution in all datasets (figure 3). We 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 shear flows.
We then studied the morphology of these coherent R structures, or 'rortices', using detailed statistics (weighted conditional averaging) on the inclination angles of R with the longitudinal and vertical axes (figures 4-8). This allowed us to draw in a key schematic (figure 9) the evolution of typical rortices from asymmetric (high-θ low-Re) and symmetric (low-θ 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, while 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 R or S with ∇ρ (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; on the contrary, 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-average 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 instability 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-16).
To synthesise these statistical and structural results, we proposed a cartoon model for the evolution of vortical structures and density interfaces in figure 17. First, we hypothesise that Λ-rortices originate from the 3-D confined Holmboe waves (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. Vortices present across the upper or lower density interface act as evolving 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 much smaller scales (not currently accessible to experimental measurements), rortices and intense shear cause turbulent kinetic energy dissipation and irreversible mixing, which are two further key energy fluxes. However, our analysis in this paper has been largely kinematic; the dynamical role of coherent rortices vs shear in shaping the turbulent energetics -and in particular the efficiency of mixing -remains an open question.