Geometry of stratified turbulent mixing: local alignment of the density gradient with rotation, shear and viscous dissipation

Abstract

We introduce a geometric analysis of turbulent mixing in density-stratified flows based on the alignment of the density gradient in two orthogonal bases that are locally constructed from the velocity gradient tensor. The first basis connects diapycnal mixing to rotation and shearing motions, building on the recent ‘rortex–shear decomposition’ in stratified shear layers (Jiang et al., J. Fluid Mech., vol. 947, 2022, A30), while the second basis connects mixing to the principal axes of the viscous dissipation tensor. Applying this framework to datasets taken in the stratified inclined duct laboratory experiment reveals that density gradients in locations of high shear tend to align preferentially (i) along the direction of minimum dissipation and (ii) normal to the plane spanned by the rortex and shear vectors. The analysis of the local alignment across increasingly turbulent flows offers new insights into the intricate relationship between the density gradient and dissipation, and thus diapycnal mixing.

1. Introduction

The physics-based parametrisation of turbulent mixing in density-stratified flows is a fundamental challenge in geophysical and environmental fluid dynamics (Caulfield 2020; Dauxois et al. 2021). This challenge requires an understanding of the small-scale mechanisms that drive transport and irreversible mixing across isopycnals (Gregg et al. 2018). These mechanisms rely on the interaction between the velocity and scalar fields. In the past, the strain-rate tensor has been used to analyse the geometry of dissipation in shear flows (see e.g. Ashurst et al. 1987). The alignment of the principal directions of the strain-rate tensor with the scalar gradient plays a crucial role in the small-scale mixing and cascade processes (Garcia, Gonzalez & Paranthoën 2005). In the context of sheared stratified turbulence, Smyth (1999) investigated how the direction of the scalar gradient aligns with the principal directions of the time-varying strain-rate tensor during the development of Kelvin–Helmholtz instability. Unlike previous numerical studies relying on the strain-rate tensor, the present work uses three-dimensional (3-D) experimental datasets and the pseudo-dissipation tensor (to be defined in § 2.4) to study the alignment statistics with the scalar gradient.

Using simultaneous 3-D velocity and density data in the ‘stratified inclined duct’ (SID) laboratory experiment, Jiang et al. (2022) analysed the morphology of coherent (Eulerian) vortical structures and their evolution from preturbulent Holmboe waves to fully developed sheared turbulence. They explained how turbulent hairpin vortices across the density interface engulf unmixed fluid into the stratified interface, while vortices within the shear layer generate further stirring and small-scale shear, enhancing irreversible mixing. Although significant emphasis has been placed on understanding the interaction between vortices and stratification in previous research (see also Watanabe et al. 2019), the alignment of velocity and density gradients under different shear strengths, and the connection with viscous dissipation, remain unclear.

The objective of this paper is to tackle these issues by using the same datasets as Jiang et al. (2022) and quantifying the relation between kinematic rotational and shearing flow structures and the dynamics and energetics of the flow, an open question (Riley 2022) essential to improve parametrisations. Specifically, we will address the two following questions. First, how is viscous dissipation connected to the homogenisation of the density field that indicates mixing? Second, how do shear and rotation, the building blocks of vorticity, contribute distinctively to the mixing of the density field? To answer these questions, we investigate the local geometric alignment of the density gradient with characteristic directions associated with the local viscous dissipation tensor and with the rotation and shear vectors.

In § 2 we introduce our experimental datasets and new geometric framework consisting of two local orthogonal bases constructed from the velocity gradient tensor. We quantify the alignment of the density gradient in the first basis (spanned by viscous dissipation eigenvectors) in § 3.1, and in the second basis (spanned by shear and rotation vectors) in § 3.2. We then study in § 3.3 how this alignment varies under increasing turbulent intensity, and in § 3.4 examine the link between alignment and a standard mixing coefficient in the most turbulent flow, emphasising the physical insights gained. Finally, we conclude in § 4.

2. Methodology

2.1. Experiment and data processing

We collected data in the SID experiment (sketched in figure 1) in which a salt-stratified, shear-driven flow is sustained in a square duct of length $1350$ mm and cross-section $H=45$ mm (aspect ratio $30$) connecting two reservoirs of different densities $\rho _0\pm \Delta \rho /2$, giving an Atwood number $At \equiv \Delta \rho /(2\rho _0)$. The duct can be tilted at a small angle $\theta$ with respect to the horizontal to add an along-duct gravitational forcing and increase turbulence levels. This experiment has been described in more detail in Meyer & Linden (2014) and Lefauve & Linden (2020). The data acquisition and processing pipeline, originally introduced in Lefauve & Linden (2022a) (§ 3.3), comprises four steps illustrated in figure 1.

Figure 1. Set-up with the four-step data acquisition and processing pipeline, transforming volumetric velocity and density measurements into non-dimensional ‘shear-layer’ datasets with a peak-to-peak velocity jump, density jump, and shear-layer height of 2.

In step 1, we simultaneously measure the time-resolved, 3-D density field $\rho ^d(\boldsymbol {x}^d,t^d)$ and three-component velocity field $\boldsymbol {u}^d=(u^d,v^d,w^d)(\boldsymbol {x}^d,t^d)$ using scanning laser-induced fluorescence and stereo particle image velocity (Partridge, Lefauve & Dalziel 2019). The superscript $d$ highlights that the data are still in the original dimensional units. The $x^d$-axis is aligned locally with the duct, making the $z^d$-axis tilted at an angle $\theta$ with respect to the true vertical (defined by the direction of gravity $\boldsymbol {g}$), while the $y^d$-axis is the spanwise direction.

In step 2, the velocity $u$ and density $\rho -\rho _0$ are then averaged in $x^d$ and in $t^d$, yielding three key profiles: the vertical profile of the streamwise velocity, $u^d(z^d)$ (in red), sampled at $y^d=0$, the vertical profile of density, $\rho ^d(z^d)$ (in grey), also sampled at $y^d=0$, and the spanwise profile of the streamwise velocity, $u^d(y^d)$ (in blue), sampled at the height corresponding to maximum flow speed, $z^d_{max}$.

In step 3, all data are then non-dimensionalised. For length we use half the shear-layer height $h_u/2$ (measured from the locations of minimum to maximum $u^d(z^d)$), for velocity we use half the peak-to-peak velocity jump $\delta u /2$, as shown in the figure, and for density we use half the fixed maximum density jump $\Delta \rho /2$. This restricts all non-dimensional shear-layer variables $z^s,u^s,\rho ^s$ between $-1$ and $1$, and is consistent with the following Reynolds and bulk Richardson numbers, respectively,

(2.1a,b)\begin{equation} Re \equiv \left( \left.\frac{h_u}{2} \frac{\delta_u}{2}\right)\right/\nu, \quad Ri_b \equiv \left.\frac{g}{\rho_0} \frac{\Delta\rho}{2}\frac{h_u}{2}\right/\left(\frac{\delta_u}{2}\right)^2. \end{equation}

Note that salt stratification yields a large Prandtl number $Pr\equiv \nu /\kappa \approx 700$, where $\nu$ is the kinematic viscosity of water and $\kappa$ the diffusivity of salt. The thickness of the density interface, defined as the spacing between the values $\rho ^d=\pm \tanh (1) = \pm 0.76$, becomes $2/R_h$, where $R_h$ is the thickness ratio between the shear layer and the density interface. This step allows datasets obtained in flows at different values of $\Delta \rho$ (or Atwood number $At$) and tilt angle $\theta$ to be compared side-by-side with accurate non-dimensional parameters based on the actual measured velocity profiles.

In step 4, we define the shear layer as the region $|z^s|\le 1$ and $|y^s|\le L_y$ (where $2L_y$ is the central portion where $|u^s|>0.7$) and discard the regions outside this. This step avoids including wall-bounded, unstratified flow (inherent to this experiment) in our statistical analysis of the interfacial, shear-driven stratified flow data of interest to this paper. In the remainder of the paper, we drop the $s$ superscripts and use shear-layer variables throughout.

2.2. Dataset

The dataset, described in Lefauve & Linden (2022a) and freely available in Lefauve & Linden (2022c), includes 15 experiments with increasing levels of turbulent intensity, controlled by the product of the tilt angle $\theta$ and Reynolds number $\mbox { {Re}}$. The data comprise four flows in the Holmboe regime (H1–H4), showing small-scale, travelling waves whose scouring motion maintain a relatively sharp density interface $R_h\in [7.2, 11.3]$; eight flows in the intermittently turbulent regime (I1–I8) with bursts of increasingly longer-lived turbulence whose overturning mixing results in a thicker interface ($R_h \in [2.2, 5.8]$); and three flows in the turbulent regime (T1–T3), where the flow never relaminarises and the density interface is thickest ($R_h\in [1.8, 2.1]$).

Our analysis focuses on two typical cases: H1, a Holmboe wave flow at $(\theta,\mbox { {Re}},Ri_b,R_h)=(1^\circ, 381, 0.567,8.9)$; and T3, a vigorously turbulent flow at $(\theta,\mbox { {Re}},Ri_b,R_h)=(5^\circ, 1145, 0.147, 1.9)$. For detailed parameters of these and the remaining datasets used in this paper, see Lefauve & Linden (2022a, table 1).

2.3. Vorticity decomposition and structural coordinates

We build on the kinematic analysis of Jiang et al. (2022) using the local and instantaneous rortex–shear decomposition of vorticity $\boldsymbol {\nabla } \times \boldsymbol {u} \equiv \boldsymbol \omega = \boldsymbol {R}+\boldsymbol {S}$, appropriate when $\boldsymbol {\nabla }\boldsymbol {u}$ has only a single real eigenvalue. The ‘rortex vector’, capturing rigid-body rotation, is defined as (Gao & Liu 2019)

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

where $\hat {\boldsymbol {r}}$ and $\lambda _{ci}$ are, respectively, the local unit real eigenvector and imaginary part of the complex-conjugate pair of eigenvalues of the velocity gradient tensor $\boldsymbol {\nabla }\boldsymbol {u}$. Points where such a complex-conjugate pair does not exist have no local rigid-body rotation and are thus excluded from our analysis. The remaining vorticity $\boldsymbol {S}\equiv \boldsymbol {\omega } - \boldsymbol {R}$ is the called ‘shear vector’ (with magnitude $S$) as it is dominated by shearing motions (especially in SID where turbulence is fed by a sustained background shear). Since the shear $\boldsymbol {S}$ is not necessarily perpendicular to $\boldsymbol {R}$, it can be further decomposed into ‘rotational shear’ $\boldsymbol {S}_r=(\boldsymbol {S} \boldsymbol {\cdot } \hat {\boldsymbol {r}})\hat {\boldsymbol {r}}$ (aligned with $\boldsymbol {R}$) and a residual ‘non-rotational shear’ $\boldsymbol {S}_{n}=\boldsymbol {S} - \boldsymbol {S}_r$, as shown in figure 2(a). We define the total rotation as $\boldsymbol {R}_t \equiv \boldsymbol {R} + \boldsymbol {S}_r$ and construct the right-handed orthonormal basis $(\hat {\boldsymbol {n}},\hat {\boldsymbol {r}},\hat {\boldsymbol {f}})$ based on the unit non-rotational shear and rotation vectors

(2.3a–c)\begin{equation} \hat{\boldsymbol{n}}= \frac{\boldsymbol{S}_{n}}{|\boldsymbol{S}_{n}|}, \quad \hat{\boldsymbol{r}} = \frac{\boldsymbol{R}_{t}}{|\boldsymbol{R}_{t}|}, \quad \hat{\boldsymbol{f}} \equiv \hat{\boldsymbol{n}} \times \hat{\boldsymbol{r}}. \end{equation}

Figure 2. Sketch of the orthogonal bases (a) $(\hat {\boldsymbol {n}}, \hat {\boldsymbol {r}}, \hat {\boldsymbol {f}})$ based on rigid-body rotation $\hat {\boldsymbol {r}}$ in the velocity gradient tensor and non-rotational vorticity $\hat {\boldsymbol {n}}$; (b) $(\hat {\boldsymbol {d}_1},\hat {\boldsymbol {d}_2},\hat {\boldsymbol {d}_3})$ based on the eigendirections of the pseudo-dissipation tensor. In (a), the red dashed arrow labelled $\boldsymbol {u}_r$ represents velocity induced by rigid-body rotation and the green band labelled $\boldsymbol {u}_n$ represents the likely velocity distribution induced by the non-rotational shear. In (b), the shaded cones represent the typical range of alignments found in our results.

2.4. Directions of dissipation

To gain dynamical information, we define the local and instantaneous pseudo-dissipation tensor (analogous to the left Cauchy–Green tensor in solid mechanics):

(2.4)\begin{equation} \boldsymbol{\mathsf{D}} \equiv \frac{2}{\mbox{{Re}}} \boldsymbol{\nabla}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}^T \quad\text{or} \quad D_{ij} \equiv \frac{2}{\mbox{{Re}}} \frac{\partial u_i}{\partial x_k}\frac{\partial u_j}{\partial x_k}. \end{equation}

There is a close relation between the rate of kinetic energy dissipation typically used in the stratified turbulence literature $\epsilon \equiv (2/\mbox { {Re}}) \boldsymbol{\mathsf{E}} : \boldsymbol{\mathsf{E}}$ (where $\boldsymbol{\mathsf{E}} \equiv (\boldsymbol {\nabla }\boldsymbol {u} + \boldsymbol {\nabla }\boldsymbol {u}^T )/2$ is the strain-rate tensor) and the pseudo-dissipation $\tilde {\epsilon }\equiv (1/\mbox { {Re}}) \boldsymbol {\nabla } \boldsymbol {u} : \boldsymbol {\nabla } \boldsymbol {u}$ (Pope 2000, § 5.3), namely

(2.5)\begin{equation} \tilde{\epsilon} = \frac{1}{2} \textrm{tr}({\boldsymbol{\mathsf{D}}}) = \frac{1}{\mbox{{Re}}} \dfrac{\partial u_i}{\partial x_j} \dfrac{\partial u_i}{\partial x_j} = \epsilon - \frac{1}{\mbox{{Re}}} \dfrac{\partial u_i}{\partial x_j} \dfrac{\partial u_j}{\partial x_i} \approx \epsilon. \end{equation}

Our motivation for using $\boldsymbol{\mathsf{D}}$ (rather than $\boldsymbol{\mathsf{E}}$) to study dissipation in this paper is that $\boldsymbol{\mathsf{D}}$ is closely related to the eigenbasis of $\boldsymbol {\nabla }\boldsymbol {u}$, upon which structural identification schemes are often built, like the above rortex–shear decomposition. The advantage of $\boldsymbol{\mathsf{D}}$ is that (unlike $\boldsymbol {\nabla }\boldsymbol {u}$) it is symmetric positive semi-definite and therefore has three real non-negative eigenvalues $\lambda _1 \ge \lambda _2 \ge \lambda _3 \ge 0$, such that $\tilde {\epsilon }=(\lambda _1+ \lambda _2+\lambda _3)/2$, and an associated basis of real, orthogonal unit eigenvectors $(\hat {\boldsymbol {d}_1},\hat {\boldsymbol {d}_2}, \hat {\boldsymbol {d}_3})$, sketched in pink in figure 2(b).

The interpretation of this basis is provided by the unique polar decomposition $\boldsymbol {\nabla }\boldsymbol {u}=\boldsymbol{\mathsf{V}} \boldsymbol {\cdot } \boldsymbol{\mathsf{R}}$, where $\boldsymbol{\mathsf{R}}$ is the rotation tensor (a proper-orthogonal tensor) and $\boldsymbol{\mathsf{V}}$ is the left stretch tensor (a real, symmetric, positive semi-definite tensor; Spencer 1980, § 9.2). Since $\boldsymbol{\mathsf{R}}$ is orthogonal then $\boldsymbol{\mathsf{D}}=\frac {2}{\mbox { {Re}}}\boldsymbol{\mathsf{V}}^2$ and $\boldsymbol{\mathsf{V}}$ has the same eigenvectors as $\boldsymbol{\mathsf{D}}$ with associated eigenvalues $\sqrt {\mbox { {Re}}\lambda _i/2}$. The linear transformation $\boldsymbol {\delta u}=\boldsymbol {\nabla }\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\delta x}$ (the first order change of velocity a small $\boldsymbol {\delta x}$ away) is decomposed into a rotation followed by a local stretch (if $\sqrt {\mbox { {Re}}\lambda _i/2}>1$) or compression (if $\sqrt {\mbox { {Re}}\lambda _i/2}<1$) along the basis vectors $\hat {\boldsymbol {d}_i}$. As the dissipation equals (half) the sum of these principal stretches squared, we interpret $\hat {\boldsymbol {d}_1},\hat {\boldsymbol {d}_2}, \hat {\boldsymbol {d}_3}$ as representing the directions of maximum, intermediate and minimum dissipation, respectively.

According to Wu, Ma & Zhou (2016), the dissipation is related to the stress resulting from the surface deformation rate, which constitutes a viscous resistance to changes in the direction and area of surface elements. Therefore, the maximum (minimum) dissipation direction can be seen as the axis along which fluid elements can efficiently (inefficiently) dissipate kinetic energy through cumulative effects of stretching (or compression) and rotation.

2.5. Density-gradient ratios

We define the ‘dissipation angles’ as $\zeta _i \equiv \angle (\boldsymbol {\nabla }\rho,\hat {\boldsymbol {d}_i})$, with $i=1, 2, 3$, and the ‘structure angles’ as $\zeta ^\phi \equiv \angle (\boldsymbol {\nabla }\rho,\hat {\boldsymbol {\phi }})$, with $\hat {\boldsymbol {\phi }}=\hat {\boldsymbol {n}},\hat {\boldsymbol {r}}, \hat {\boldsymbol {f}}$. The density gradient can then be decomposed either along the directions of dissipation defined by the unit vectors $(\hat {\boldsymbol {d}_1},\hat {\boldsymbol {d}_2},\hat {\boldsymbol {d}_3})$,

(2.6)\begin{equation} \boldsymbol{\nabla}\rho = \sum_{i=1,2,3} \boldsymbol{\nabla}{\rho}_i, \quad \text{where} \ \boldsymbol{\nabla}{\rho}_i = |\boldsymbol{\nabla}{\rho}|\cos(\zeta_i)\hat{\boldsymbol{d}_i}, \end{equation}

or along the structural direction as

(2.7)\begin{equation} \boldsymbol{\nabla}\rho = \sum_{\phi=n,r,f} \boldsymbol{\nabla}{\rho}^\phi, \quad \text{where} \ \boldsymbol{\nabla}{\rho}^\phi = |\boldsymbol{\nabla}{\rho}|\cos(\zeta^\phi)\hat{\boldsymbol{\phi}}. \end{equation}

In (2.6), each component $\boldsymbol {\nabla }{\rho }_i$ can be further decomposed along $(\hat {\boldsymbol {n}},\hat {\boldsymbol {r}},\hat {\boldsymbol {f}})$:

(2.8)\begin{equation} \boldsymbol{\nabla} \rho_{i}= \sum_{\phi=n,r,f} \boldsymbol{\nabla} \rho_{i}^\phi, \quad \text{where} \ \boldsymbol{\nabla} \rho_{i}^\phi = (\boldsymbol{\nabla} \rho_{i} \boldsymbol{\cdot} \hat{\boldsymbol{\phi}})\hat{\boldsymbol{\phi}}. \end{equation}

The 12 components in (2.6) and (2.8) encode the alignment of the density gradient with 6 physically meaningful directions corresponding to dissipation and rotation/shear. We propose the following $3+3+9=15$ squared-density-gradient ratios (SDGRs):

(2.9a–c)\begin{equation} \mathscr{M}_{i} \equiv \frac{|\boldsymbol{\nabla} \rho_i|^2}{|\boldsymbol{\nabla} \rho|^2} = \cos^2\zeta_i, \quad \mathscr{M}^\phi \equiv \frac{|\boldsymbol{\nabla} \rho^\phi|^2}{|\boldsymbol{\nabla} \rho|^2} = \cos^2\zeta^\phi, \quad \mathscr{M}_{i}^\phi \equiv \frac{|\boldsymbol{\nabla} \rho^\phi_{i}|^2}{|\boldsymbol{\nabla} \rho_{i}|^2}. \end{equation}

Here, repeated indices do not imply summation; where summation is required, we use explicitly $\sum _j$. The last nine SDGRs encode the projection of $(\hat {\boldsymbol {d}_1},\hat {\boldsymbol {d}_2},\hat {\boldsymbol {d}_3})$ in $(\hat {\boldsymbol {n}},\hat {\boldsymbol {r}},\hat {\boldsymbol {f}})$ and $\mathscr{M}_{i}^\phi$ is related to the rotation matrix $\boldsymbol{\mathsf{T}}_{i\phi }$ between these bases. The use of orthonormal bases implies that each group of SDGR sum to one, namely $\sum _{i} \mathscr {M}_i = \sum _{\phi } \mathscr {M}^\phi = \sum _{\phi } \mathscr {M}^\phi _1= \sum _{\phi } \mathscr {M}^\phi _2 = \sum _{\phi } \mathscr {M}^\phi _3 = 1$.

3. Results

3.1. Alignment of density gradient with directions of dissipation

We start by examining the local alignment between $\boldsymbol {\nabla }\rho$ and $\hat {\boldsymbol {d}_1},\hat {\boldsymbol {d}_2},\hat {\boldsymbol {d}_3}$ quantified by the three angles $\zeta _i$. Since $\cos ^2{\zeta _2} =1-\cos ^2{\zeta _1} - \cos ^2{\zeta _3}$, we focus our analysis on the two independent angles $\zeta _1$ and $\zeta _3$.

Figure 3 compares the joint probability density function (p.d.f.) of $\cos \zeta _1$ and $\cos \zeta _3$ in the Holmboe wave flow (dataset H1, panel a) and turbulent flow (dataset T3, panel b). The blue lines represent conditional p.d.f.s for the two angles separately at increasingly large values of $k=S/S_{rms}$, where $S_{rms}$ is the root-mean-square shear averaged in space and time, starting with the interval $(0,0.5]$ (lighter shade) and ending with $(2.5,3]$ (darkest shade). Note that the high-shear region often corresponds to a higher dissipation (dominated by $\lambda _1$). We find that, in both flows, $\boldsymbol {\nabla } \rho$ is most often nearly anti-parallel to $\hat {\boldsymbol {d}_3}$ ($\cos \zeta _3 \approx -1$) and nearly perpendicular to $\hat {\boldsymbol {d}_1}$ ($\cos \zeta _1 \approx 0$) and, therefore, nearly perpendicular to $\hat {\boldsymbol {d}_2}$. This trend is clear even in the turbulent flow, despite more spread around the mean values than in the Holmboe flow.

Figure 3. Joint p.d.f.s of the alignment between $\boldsymbol {\nabla }\rho$ and the direction of maximum and minimum dissipation $\hat {\boldsymbol {d}_1}$, $\hat {\boldsymbol {d}_3}$ (angles $\zeta _1,\zeta _3$ respectively) in H1 (a) and T3 (b). Conditional p.d.f.s are shown in blue, with darker lines indicating higher shear thresholds $k=S/S_{rms}\in (0,0.5], \ldots, (2.5,3]$. The vertical scale of blue lines has been amplified by 3 (a) and 1.5 (b) for better visualisation. The top-right insets show the joint p.d.f. conditioned at $k\in (2, 2.5]$. The top-left insets compare the alignment of $\boldsymbol {\nabla }\rho$ with the intermediate and minimum strain directions ($\hat {\boldsymbol {e}_2}$, $\hat {\boldsymbol {e}_3}$) to the alignment of $\hat {\boldsymbol {d}_2}$ with $\boldsymbol {\omega }$.

The conditional p.d.f.s and the top-right insets, restricting the joint p.d.f.s to high-shear regions $k\in (2,2.5]$, show that this trend increases in regions of high shear. Figure 2 illustrates (with shaded cones) this preferential alignment of the density gradient $\boldsymbol {\nabla }\rho$ along the (negative) direction of least dissipation $\hat {\boldsymbol {d}_3}$ and thus perpendicular to $\hat {\boldsymbol {d}_1},\hat {\boldsymbol {d}_2}$.

The top-left insets of figure 3 shows the p.d.f. of the alignment of $\boldsymbol {\nabla }\rho$ with the intermediate and minimum eigendirections of the strain-rate tensor ($\hat {\boldsymbol {e}_2}$ and $\hat {\boldsymbol {e}_3}$). It shows a preferential $\boldsymbol {\nabla }\rho \perp \hat {\boldsymbol {e}_2}$ (similarly to $\hat {\boldsymbol {d}_2}$) but a deviation from alignment to $\hat {\boldsymbol {e}_3}$ by $\approx 40^\circ$ (see black arrows). This observation aligns with prior findings in sheared, stratified flow, which have shown imperfect alignment between $\boldsymbol {\nabla } \rho$ and the compressive strain direction (Ashurst et al. 1987; Smyth 1999). This finding implies that both maximum stretching and compressive strain influence the orientation of $\boldsymbol {\nabla }\rho$ in the continuously shear-driven, high-$Pr$ SID flow, promoting alignment with the direction of minimal dissipation.

Isopycnals oriented in the directions of maximum and intermediate dissipation (i.e. $\boldsymbol {\nabla }\rho \parallel \hat {\boldsymbol {d}_3}$) are particularly susceptible to stretching and diffusion. The reduced alignment between $\boldsymbol {\nabla }\rho$ and $\hat {\boldsymbol {d}_3}$ in T3, in contrast with H1, suggests increased turbulence-induced energy dissipation, resulting in a more effectively mixed shear layer. Therefore, the deflection of $\boldsymbol {\nabla }\rho$ from the $\hat {\boldsymbol {d}_3}$ direction to the $\hat {\boldsymbol {d}_1}$–$\hat {\boldsymbol {d}_2}$ plane is an outcome of this mixing process.

We also find that the vorticity $\boldsymbol {\omega }$ has a strong and robust preferential alignment along $\hat {\boldsymbol {d}_2}$ in both flows (see top-left insets of figure 3). This echoes earlier studies that observed vorticity aligning with the intermediate principal direction $\hat {\boldsymbol {e}_2}$ of the strain-rate tensor in both shear and isotropic flows (Ashurst et al. 1987). It implies that the intermediate directions of the strain and dissipation tensors may share comparable functions related to the net production of enstrophy.

3.2. Alignment of density gradient with rotation and shear

We now examine the local alignment between $\boldsymbol {\nabla }\rho$ and the rotation/shear basis using the two independent angles $\zeta ^r,\zeta ^n$ in figure 4, similar to figure 3. In both flows H1 and T3, two regions have high probability density (see arrows I and II). Region I corresponds to the rim where $\cos ^2{\zeta ^n} + \cos ^2{\zeta ^r} \approx 1$, i.e. nearly within the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane and thus $\boldsymbol {\nabla }\rho \perp \hat {\boldsymbol {f}}$. Note that $\hat {\boldsymbol {f}} = \hat {\boldsymbol {n}} \times \hat {\boldsymbol {r}}$ and corresponds to the direction of non-rotational straining (shearless) motions. Region II corresponds to the centre where $\cos {\zeta ^n} = \cos {\zeta ^r} \approx 0$, i.e. $\boldsymbol {\nabla } \rho \parallel \hat {\boldsymbol {f}}$. These two alignment properties between the $\boldsymbol {\nabla } \rho$ and $\hat {\boldsymbol {f}}$ indicate two distinct states of mixing, which we will elaborate on later in the paper.

Figure 4. Joint p.d.f.s of the alignment between $\boldsymbol {\nabla }\rho$ and non-rotational shear $\hat {\boldsymbol {n}}$ (angle $\zeta ^n$) and rotation $\hat {\boldsymbol {r}}$ (angle $\zeta ^r$) for H1 (a) and T3 (b), similar to figure 3. The top-left and top-right insets show the joint p.d.f. conditioned at $k\in (0,0.5]$ and $k\in (2, 2.5]$, respectively, with the dashed line in (b) indicating the asymmetry of the distribution.

We can gain further insight by conditioning the statistics of $\zeta ^n$ and $\zeta ^r$ on the shear $S$, as in the previous section. The blue curves in figure 4(a,b) show that $\cos \zeta ^n$ and $\cos \zeta ^r$ both become more confined to 0 with high shear (darker blue curves), thus belonging increasingly to region II. This is also evidenced by the joint p.d.f.s conditioned at $k\in (2,2.5]$ (top-right insets), showing a strong peak at the origin. By contrast, the joint p.d.f.s conditioned at $k\in (0,0.5]$ (top-left insets) show that $\zeta ^n$ and $\zeta ^r$ belong to region I. This suggests that $\boldsymbol {\nabla } \rho \parallel \hat {\boldsymbol {f}}$ in regions of high shear and $\boldsymbol {\nabla } \rho \perp \hat {\boldsymbol {f}}$ in regions of low shear. The asymmetry seen in the insets for the turbulent case (see dashed lines in I and II in figure 4b) also suggest that $\boldsymbol {\nabla } \rho$ is slightly more aligned with $\hat {\boldsymbol {n}}$ than with $\hat {\boldsymbol {r}}$.

Although Jiang et al. (2022) discovered that $\boldsymbol {\nabla }\rho$ vectors frequently exhibit perpendicular orientation relative to $\hat {\boldsymbol {r}}$ within regions of pronounced stratification, the present findings highlight the importance of the alignment between $\boldsymbol {\nabla } \rho$ and the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane, and the dependence on shear strength. We recall that in SID flows, the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane is (statistically) preferentially inclined by $0^\circ$–$15^\circ$ to the ‘true horizontal’ plane (Jiang et al. 2022), i.e. approximately $\hat {\boldsymbol {f}} \parallel \hat {\boldsymbol {g}}$ (gravity). In a stably stratified shear layer, we expect preferentially $\boldsymbol {\nabla } \rho \parallel \hat {\boldsymbol {g}}$ (hydrostatic equilibrium), and therefore $\boldsymbol {\nabla } \rho \parallel \hat {\boldsymbol {f}}$, which is what is found at high-shear regions where the stratification exerts sufficient strength to mitigate the stirring caused by neighbouring vortical structures. However, in regions with weak shear, such as the middle layer of turbulent flow, $\boldsymbol {\nabla } \rho$ reorients from $\boldsymbol {\nabla } \rho \parallel \hat {\boldsymbol {f}}$ to $\boldsymbol {\nabla } \rho \parallel \hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane (see figure 2), which arises as a result of nearby flow motions that distort local fluid parcels and disrupt isopycnals. In these unstable regions, figure 4(b) suggests that rotation is more effective than non-rotational shear in distorting the isopycnals. We believe that the intricate interplay between these flow motions and buoyancy may affect the local instability due to their different time scales.

3.3. Squared-density-gradient ratios

We now examine in figures 5 and 6 the alignment of $\boldsymbol {\nabla }\rho$ in the two orthogonal bases using the SDGRs defined in (2.9a–c). The overbar $\overline{\mathscr {M}}$ averages an SDGR over the entire available shear-layer volume ($x,y,z$) and time ($t$), as in figure 5, whereas brackets $\langle \mathscr {M}\rangle (z)$ average in $x$, $y$ and $t$, retaining the $z$ dependence, as in figure 6(b).

Figure 5. Averaged SDGRs defined in (2.9a–c) for the 15 experimental datasets. (a–d) Variation with $\theta \mbox { {Re}}$ of $\overline {\mathscr {M}_{i}}$ (a), $\overline {\mathscr {M}^\phi }$ (b), $\overline {\mathscr {M}_{1}^\phi }$ (c), $\overline {\mathscr {M}}_{3}^\phi$ (d).

Figure 6. (a) Schematic diagram of overall values of the SDGRs represented by the thickness of arrows and text. (b) Vertical profiles of the seven key SDGRs for the turbulent flow T3, together with the buoyancy frequency $N^2$. The blue shaded region in (b) is re-used in the analysis of figure 7.

Figure 7. Correlation between the turbulent flux coefficient $\varGamma$ and the six SDGRs (a) $\langle \mathscr {M}_{1}\rangle, \langle \mathscr {M}_{1}^r\rangle, \langle \mathscr {M}_{1}^f\rangle$, (b) $\langle \mathscr {M}_{3}\rangle, \langle \mathscr {M}_{3}^r\rangle, \langle \mathscr {M}_{3}^f\rangle$ in the T3 case within the shear layer ($z\in [-0.5,0]$). The size of symbols indicates $|z|$, and the darker colours indicate higher stratification $N^2$.

Figure 5 compares the averaged SDGRs in the 15 flows ranging from H1 to T3, with increasing turbulent intensity controlled by the product $\theta \mbox { {Re}}$ (where $\theta$ is in radians). In panel (a), we find that typically $\overline {\mathscr {M}_3}>\overline {\mathscr {M}_2}>\overline {\mathscr {M}_1}$, confirming the results of § 3.1 that $\boldsymbol {\nabla }\rho \parallel \hat {\boldsymbol {d}_3}$ is preferential, although this trend of dominant $\overline {\mathscr {M}_3}$ decreases slightly as $\theta \mbox { {Re}}$ increases. For $\overline {\mathscr {M}^\phi }$ in panel (b), the distribution is similar to $\overline {\mathscr {M}_i}$ in (a). Both $\overline {\mathscr {M}^n}$ and $\overline {\mathscr {M}^r}$ show a nearly linear increase as $\overline {\mathscr {M}^f}$ decreases, indicating that the impact of shear and rortices on density becomes more pronounced with intense turbulence. The reason behind the deviation of $\boldsymbol {\nabla }\rho$ from $\hat {\boldsymbol {d}_3}$ or $\hat {\boldsymbol {f}}$ at higher $\theta \mbox { {Re}}$ can be attributed to more intense mixing and overturning.

In figure 5(c), we find that $\overline {\mathscr {M}_1^r}>\overline {\mathscr {M}_1^n}>\overline {\mathscr {M}_1^f}$, indicating that the component of the density gradient along the direction of maximum dissipation $\boldsymbol {\nabla } \rho _{1}$ aligns preferentially with rotation $\hat {\boldsymbol {r}}$, and secondarily with non-rotational shear $\hat {\boldsymbol {n}}$. As $\overline {\mathscr {M}_1^f} \le 20\,\%$, it is likely that $d_1$ lies close to the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane. If $\hat {\boldsymbol {r}}\parallel \hat {\boldsymbol {d}_1}$, both a rotation and stretching occur along the $\hat {\boldsymbol {d}_1}$ direction, strengthening the velocity induced by the non-rotational shear, as depicted in figure 2(a). This suggests that the stretching of the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane, especially along the $\hat {\boldsymbol {r}}$ direction, is closely related to the dissipation.

In figure 5(d), we see that the preferential alignment is reversed for $\boldsymbol {\nabla } \rho _3$: $\overline {\mathscr {M}_{3}^f}>\overline {\mathscr {M}_{3}^n}>\overline {\mathscr {M}_{3}^r}$, i.e. $\boldsymbol {\nabla}\rho_3$ is primarily aligned with $\hat {\boldsymbol {f}}$. In the most turbulent flows T2 and T3, $\overline {\mathscr {M}_{3}^f} \approx \overline {\mathscr {M}_{3}^n} \approx 2 \overline {\mathscr {M}_{3}^r}$, i.e. $\boldsymbol {\nabla } \rho _3$ is preferentially in the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {f}}$ plane and perpendicular to $\hat {\boldsymbol {r}}$. We recall from §§ 3.1–3.2 that high-shear regions have approximately $\boldsymbol {\nabla } \rho \parallel \hat {\boldsymbol {f}} , \hat {\boldsymbol {d}_3} , \hat {\boldsymbol {g}}$, thus $\boldsymbol {\nabla } \rho _{3}$ is the component of the gradient most likely to be disturbed during mixing. Having $\boldsymbol {\nabla }\rho _3 \perp \hat {\boldsymbol {r}}$ signals stirring and overturning, a requirement for diapycnal mixing by small-scale diffusion.

As $\theta$ increases at nearly constant $\mbox { {Re}}$, the $\hat {\boldsymbol {d}_3}$–$\boldsymbol {\nabla }\rho$ alignment decreases, as well as the $\hat {\boldsymbol {f}}$–$\boldsymbol {\nabla }\rho$ and $\hat {\boldsymbol {d}_3}$–$\hat {\boldsymbol {f}}$ alignments. This is observed in I3, I7, I8 and T3 with $\theta =2^\circ$–$6^\circ$ and $\mbox { {Re}}\approx 1000$. Conversely, increasing $\mbox { {Re}}$ at constant $\theta =5^\circ$ has a more limited impact on alignment, as seen in cases H2, H4, I5, I8 and T3. This suggests that $\theta$ exerts a more significant influence than $\mbox { {Re}}$ on $\boldsymbol {\nabla } \rho$ alignment, echoing the findings of Lefauve & Linden (2022a) that $\theta$ increases overturning more strongly than $Re$.

Figure 6(a) summarises the typical averaged SDGRs across all 15 flows and encapsulates the typical geometry of stratified turbulence in our experimental data. The left-hand branches (i.e. the three arrows under the blue shading indicated by $\overline {\mathscr {M}_i}$) show the three ‘dissipation’ SDGRs, while the right-hand branches show the three ‘structural’ SDGRs (indicated by $\overline {\mathscr {M}^\phi }$), and the middle branches show the nine ‘mixed’ SDGRs indicated by $\overline {\mathscr {M}_i^\phi }$. Note the approximate symmetry between the left- and right-hand branches, indicating a robust relation between diapycnal transport, structural and dissipation coordinates.

Finally, figure 6(b) shows the variation of seven key SDGRs along $z$ (in the duct frame of reference) in the turbulent flow T3. To examine the effect of stratification, we superimpose the averaged non-dimensional buoyancy frequency $\langle N^2\rangle =-Ri_b\partial _z \rho$ (thick grey line). We find that $\langle {\mathscr {M}}_{1}\rangle$ (thick red line) is maximum in the middle layer where $\langle N^2\rangle$ is minimum, indicating that $\boldsymbol {\nabla } \rho$ is more distorted and thus more aligned with $\hat {\boldsymbol {d}_1}$ in the turbulent mixing layer. Vice versa, $\langle \mathscr {M}_{1}\rangle$ is minimum at the upper and lower density interfaces of this mixing layer where $\langle N^2\rangle \gtrsim 0.2$. This trend is exactly reversed for $\langle {\mathscr {M}}_{3}\rangle$, which reaches a minimum value in the well-mixed region and a maximum value at the interfaces. Furthermore, $\langle \mathscr {M}_{2}\rangle$ is relatively uniform in $z$ and independent of the stratification. We also observe the same pattern of reversed behaviour between $\langle \mathscr {M}_1^r\rangle$ and $\langle \mathscr {M}_1^f\rangle$, as well as between $\langle \mathscr {M}_3^r\rangle$ and $\langle \mathscr {M}_3^f\rangle$. Recalling from § 3.2 that $\boldsymbol {\nabla }\rho \perp \hat {\boldsymbol {r}}$ at high shear, the finding here that stronger stratification favours the alignment of $\hat {\boldsymbol {r}}$ and $\hat {\boldsymbol {d}_1}$ (larger $\langle \mathscr {M}_1^r\rangle$) is consistent with the fact that more diapycnal transport occurs at the interface. This agrees with Jiang et al. (2022) who studied the interaction of hairpin-like structures and density gradients through $\boldsymbol {R}\times \boldsymbol {\nabla }\rho$, as well as with Riley, Couchman & de Bruyn Kops (2023), who found that most potential energy dissipation occurs near density interfaces.

3.4. Link with mixing coefficient

Figure 7 investigates the relationship between the SDGRs and the flux coefficient $\varGamma \equiv \bar {\mathcal {B}}/\bar {\epsilon '}$, the ratio of the globally averaged turbulent buoyancy flux $\mathcal {B} \equiv Ri_b\, w'\rho '$ and turbulent kinetic energy dissipation $\epsilon ' \equiv (2/Re)\, e'_{ij} e'_{ij}$, where the perturbations (prime variables) are taken with respect to the $(x,t)$ average (Lefauve & Linden 2022b, § 2.2). We restrict the analysis to the T3 case within the active turbulent lower layer $z\in [-0.5,0]$ (blue shaded region in figure 6b). Larger symbols indicate higher $|z|$, and darker colours indicate high $N^2$ (in this region, $N^2\propto |z|$).

First, the linear relations between $\varGamma$ and the SDGRs suggest a clear link between the efficiency of mixing, approximated by this flux coefficient (Caulfield 2020), and the geometry encapsulated in the SDGRs. Second, $\varGamma$ decreases in proportion to $\langle \mathscr {M}_1\rangle =\cos ^2\angle (\boldsymbol {\nabla }\rho,\hat {\boldsymbol {d}_1})$ and $\langle \mathscr {M}_1^f\rangle$ but increases in proportion to $\langle \mathscr {M}_{1}^r\rangle$ (figure 7a). Looking at the direction of minimal dissipation, the opposite is observed: $\varGamma$ is proportional to $\langle \mathscr {M}_3\rangle =\cos ^2\angle (\boldsymbol {\nabla }\rho,\hat {\boldsymbol {d}_3})$ and $\langle \mathscr {M}_3^f\rangle$ but inversely proportional to $\langle \mathscr {M}_{3}^r\rangle$. We interpret this by the expectation that low $\mathscr {M}_1$ or high $\mathscr {M}_3$ indicates the alignment of isopycnals along the direction of greatest dissipation (or stretch, see § 2.4), and hence an increase in the surface area of isopycnals and diffusive mixing. Third, figure 7 shows that closer to the density interface (larger darker symbols), $\hat {\boldsymbol {r}}$ aligns more with $\hat {\boldsymbol {d}_1}$ (i.e. $\langle \mathscr {M}_{1}^r\rangle$ becomes larger), while $\boldsymbol {\nabla } \rho$ aligns more with $\hat {\boldsymbol {f}}$ and $\hat {\boldsymbol {d}_3}$ (i.e. $\langle \mathscr {M}_{3}^f\rangle$ and $\langle \mathscr {M}_{3}\rangle$ become larger). This results in the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane aligning more with the $\hat {\boldsymbol {d}_1}$–$\hat {\boldsymbol {d}_2}$ plane, indicating stronger overturning and stretching leading to mixing.

4. Conclusions

In summary, we have explored the alignment of the density gradient $\boldsymbol {\nabla } \rho (\boldsymbol {x},t)$ in two local, instantaneous orthogonal bases determined from flow kinematics based on the velocity gradient tensor $\boldsymbol{\nabla} \boldsymbol{{u}}(\boldsymbol{{x}},t)$. Our goal was to shed light on the geometry of shear-driven, stably stratified turbulent mixing in a SID experiment. We used the datasets of increasingly turbulent flows investigated in Jiang et al. (2022), with a focus on two main flows representative of the Holmboe wave and turbulent regimes.

The first kinematic basis $(\hat {\boldsymbol {n}}, \hat {\boldsymbol {r}}, \hat {\boldsymbol {f}})$, constructed using the non-rotational shear $\hat {\boldsymbol {n}}$ and rigid-body rotation $\hat {\boldsymbol {r}}$, showed that in regions of weak shear, $\boldsymbol {\nabla }\rho$ preferentially aligns in the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane. However, increasing shear strength tends to reorient $\boldsymbol {\nabla } \rho$ normal to the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane, thus increasing the potential for destruction of stratification by rotational and shearing structures. The second kinematic basis $(\hat {\boldsymbol {d}_1}, \hat {\boldsymbol {d}_2}, \hat {\boldsymbol {d}_3})$, constructed using the principal directions of the pseudo-dissipation tensor (2.4), showed that $\boldsymbol {\nabla }\rho$ aligns preferentially with the direction of minimum dissipation $\hat {\boldsymbol {d}_3}$ and thus normal to the maximum dissipation $\hat {\boldsymbol {d}_1}$, a trend that is exacerbated in regions of high shear.

The SDGRs (see (2.6)–(2.9a–c)) quantified the alignment of $\boldsymbol {\nabla } \rho$ in both kinematic bases to reveal that approximately $\overline {\mathscr {M}_{1}}\approx 23\,\%$ of the total SDGR is along the maximum dissipation $\hat {\boldsymbol {d}_1}$, which is itself dominated by flow structures in the $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane. By contrast, a ratio $\overline {\mathscr {M}}_{3}\approx 50\,\%$ is along the minimum dissipation $\hat {\boldsymbol {d}_3}$, which aligns preferentially perpendicular to $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ and makes the isopycnals more susceptible to being distorted through overturning by vortices and hence more diapycnal mixing.

Focusing on the most turbulent flow, the variation of SDGRs across the shear layer explained why more efficient mixing, quantified by $\varGamma$, occurs at the edges of the mixing layer. At the edges (stronger stratification), $\hat {\boldsymbol {r}}$ aligns more with the maximum dissipation $\hat {\boldsymbol {d}_1}$, and $\boldsymbol {\nabla } \rho$ aligns more with the minimum dissipation $\hat{\boldsymbol {d}_3}$ and, therefore, is prone to overturning by rotational structures. The $\hat {\boldsymbol {n}}$–$\hat {\boldsymbol {r}}$ plane becomes more aligned with the $\hat {\boldsymbol {d}_1}$–$\hat {\boldsymbol {d}_2}$ plane, enhancing stirring and overturning and increasing flux coefficient. This mechanism explains the robust linear correlations between the mixing efficiency $\varGamma$ and the SDGRs shown in figure 7. Our findings also rationalise that $\varGamma$ decreases with increasing $\mathscr {M}_{1}$ and increases with increasing $\mathscr {M}_{3}$ because isopycnals aligned along the direction of maximum dissipation are more exposed to stretching and diffusion. These insights into the geometry of the density gradient and dissipation constitute a step towards a better physical understanding of turbulent mixing in stratified shear flows. Future research should study the applicability of these findings to stratified flows that differ from the present stably stratified, shear-driven inclined duct flows.