2.1. Field observations

Over the past decades, field observations have provided much data and insight on a variety of geophysical shear-driven stratified turbulent flows.

River plumes are outflows of buoyant water into the coastal ocean primarily forced by freshwater runoff (McPherson, Stevens & O'Callaghan Reference McPherson, Stevens and O'Callaghan2019) and/or wind (Yoshida et al. Reference Yoshida, Ohtani, Nishida and Linden1998). The strength and spatial heterogeneity of turbulent mixing between these two water masses impact the physical, chemical and biological properties of the developing coastal current (MacDonald, Carlson & Goodman Reference MacDonald, Carlson and Goodman2013).

Exchange flows between reservoirs of fluids at different densities are also highly relevant and occur on a variety of scales. At small scales, Lawrence et al. (Reference Lawrence, Pieters, Zarembaa, Tedford, Gu, Greco and Hamblin2004) investigated the exchange flow through a shallow ship canal connecting a small harbour to a lake undergoing seasonal, wind-driven cool upwelling, and the effects of this exchange on lake-shore pollution. At larger scales, the strongly stratified exchange flows in estuaries are primarily forced by periodic tides (Geyer & Smith Reference Geyer and Smith1987; Peters & Bokhorst Reference Peters and Bokhorst2000; MacDonald & Horner-Devine Reference MacDonald and Horner-Devine2008; Tedford et al. Reference Tedford, Carpenter, Pawlowicz, Pieters and Lawrence2009; Geyer et al. Reference Geyer, Lavery, Scully and Trowbridge2010). At even larger scales, the relatively steady baroclinic exchange flows through straits are weakly modulated by tides and influenced by the Earth's rotation, such as the much-studied strait of Gibraltar (Armi & Farmer Reference Armi and Farmer1988; Farmer & Armi Reference Farmer and Armi1988; Wesson & Gregg Reference Wesson and Gregg1994; Macias et al. Reference Macias, Garcia, Navas, Vazquez-Lopez-Escobar and Mejias2006).

In the deep Atlantic ocean, sill overflows of cold, Antarctic bottom water (AABW) through fractures such as the Romanche Trench are responsible for significant transport and mixing across ocean basins (Ferron et al. Reference Ferron, Mercier, Speer, Gargett and Polzin1998; van Haren et al. Reference van Haren, Gostiaux, Morozov and Tarakanov2014). In the upper-equatorial Pacific Ocean, deep periodic turbulent mixing events are caused by the interaction of a sustained vertical shear (between the wind-driven surface current and the opposing deep equatorial under-current) with a stable stratification modulated by diurnal solar heating (Smyth, Moum & Nash Reference Smyth, Moum and Nash2011; Smyth et al. Reference Smyth, Moum, Li and Thorpe2013, Reference Smyth, Pham, Moum and Sarkar2017).

Although field observations yield the most ‘realistic’ data that one can hope for, they come at the cost of a limited control over the flow parameters, of great complexity in geometry and external forcing (wind, sun, tides, buoyancy, rotation), and of limited measurement abilities, all of which add up to make their general understanding challenging.

2.2. Numerical simulations

A complementary approach is to isolate physical mechanisms by controlling the flow parameters, geometry and forcing conditions in direct numerical simulations (DNSs) of the three-dimensional governing equations.

One of the key idealised models is the ‘stratified shear layer’, or unforced parallel shear flow with hyperbolic tangent profiles for the velocity $\boldsymbol {u} = (u(z),0,0)$ and density $\rho (z)$, free slip in the vertical direction and periodicity in the streamwise directions (see e.g. Smyth & Moum Reference Smyth and Moum2000). Such mixing layers are prone to a range of linear instabilities, even in the presence of a single density interface; in particular the Kelvin–Helmholtz instability, whose initially two-dimensional billows undergo a ‘zoo’ of secondary three-dimensional instabilities mediating the transition to turbulence at $Re=O(10^3)$ (Caulfield & Peltier Reference Caulfield and Peltier2000; Mashayek & Peltier Reference Mashayek and Peltier2012a,Reference Mashayek and Peltierb, Reference Mashayek and Peltier2013). Mixing layers have complicated turbulent and mixing properties dependent on parameters such as the Reynolds, bulk Richardson and Prandtl numbers (Salehipour & Peltier Reference Salehipour and Peltier2015; Salehipour, Peltier & Mashayek Reference Salehipour, Peltier and Mashayek2015; Salehipour, Caulfield & Peltier Reference Salehipour, Caulfield and Peltier2016; Watanabe, Riley & Nagata Reference Watanabe, Riley and Nagata2017; Salehipour, Peltier & Caulfield Reference Salehipour, Peltier and Caulfield2018). The lack of forcing in these studies means that the turbulence is (at best) quasi-steady during a relatively short time before the initial kinetic energy is dissipated.

More recent studies focused on continuously forced turbulence, using boundary forcing in the stratified plane Couette flow (Zhou et al. Reference Zhou, Taylor, Caulfield and Linden2017b; Zhou, Taylor & Caulfield Reference Zhou, Taylor and Caulfield2017a), and using a relaxation of the mean profiles to initial conditions in the stratified shear layer flow (Smith, Caulfield & Taylor Reference Smith, Caulfield and Taylor2021).

Although these studies provided exceptionally detailed quantitative insight, it remains challenging to perform continuously forced simulations of flows satisfying all criteria in § 1 with parameters relevant to field observations (in particular the Reynolds and Prandtl numbers, see Appendix A, table 2). More fundamentally, simulations are approximations of idealised equations, which typically assume (among others): no rotation; incompressibility; the Boussinesq approximation; a linear equation of state; a single active scalar; spatially homogenous momentum; and scalar diffusivity with idealised values of $Pr$; as well as simplistic geometry, initial conditions, boundary conditions, and forcing.

2.3. Laboratory experiments

Laboratory experiments offer a valuable intermediate approach, by allowing more control over flow parameters, geometry, forcing and measurements than in the field, while retaining some of the inherent complexity of ‘real’ flows discarded in simulations.

A few laboratory flows satisfying all seven criteria in § 1 have been studied (see Appendix A), typically using a combination of qualitative flow visualisations and quantitative single-plane velocity/density data at relatively low resolution in space and/or time. Strang & Fernando (Reference Strang and Fernando2001) studied the entrainment at a weakly turbulent interface in a closed-loop recirculating flume driven by disk pumps (known as the Kovasznay flume after Odell & Kovasznay (Reference Odell and Kovasznay1971)). Odier et al. (Reference Odier, Chen, Rivera and Ecke2009), Odier, Chen & Ecke (Reference Odier, Chen and Ecke2014) and Odier & Ecke (Reference Odier and Ecke2017) studied the similar problem of entrainment and mixing of a turbulent wall jet developing into a sloping gravity current over a dense quiescent layer. Meyer & Linden (Reference Meyer and Linden2014) and Lefauve & Linden (Reference Lefauve and Linden2020) studied the transitions between flow regimes in exchange flows taking place in an inclined duct.

The added value of laboratory experiments such as those cited above in the three-pronged (observational, numerical, experimental) approach has so far been somewhat limited by the challenge of obtaining high-resolution, three-dimensional measurements of turbulent flow fields.

However, such measurements are now becoming available. The novel scanning stereo particle image velocimetry–planar laser-induced fluorescence (PIV–PLIF) system introduced in Partridge, Lefauve & Dalziel (Reference Partridge, Lefauve and Dalziel2019) achieves simultaneous measurements of the density and three-component velocity fields in a three-dimensional volume. Using this novel system in the stratified inclined duct (SID) geometry, Lefauve et al. (Reference Lefauve, Partridge, Zhou, Caulfield, Dalziel and Linden2018) studied the three-dimensional stability properties of interfacial Holmboe waves, and Lefauve, Partridge & Linden (Reference Lefauve, Partridge and Linden2019a) studied the time- and volume-averaged energy budgets of 16 data sets spanning a range of flows on either side of the turbulent transition.

However, although it is in principle ‘easier’ to generate high-$Re$, high-$Pr$ flows in the laboratory than in numerical simulations, it is in practice extremely challenging to obtain sufficiently resolved density and velocity measurements of such flows.

2.4. Approach

To achieve the objectives set out in § 1, we will further analyse the 16 experimental data sets of Lefauve et al. (Reference Lefauve, Partridge and Linden2019a,b). These cutting-edge density and velocity data are ideally suited for our purpose since they are non-intrusive, three-dimensional and three-component, simultaneous, high-resolution (in space and time) and accurate.

The key methodological differences between this paper and Lefauve et al. (Reference Lefauve, Partridge and Linden2019a) are: (i) our focus, in this paper, on turbulent fluctuations and statistics inside the shear layer (as opposed to volume-averages including wall effects); (ii) our analysis of these data in a framework consistent with the observational and numerical literature on stratified shear layers (in particular the non-dimensional notation), allowing for more direct comparison and added value to the general community. We introduce this methodology in the next section.

3.1. The SID experiment

We consider the stratified inclined duct (SID) experiment sketched in figure 1(a). We study the steady-state exchange flow sustained inside a long ($L=1350$ mm) duct of square cross-section ($H=45$ mm), inclined at a small angle $\theta$, connecting two large reservoirs initially filled with aqueous salt solutions ($Pr=700$) of different densities $\rho _0 \pm \Delta \rho /2$. This exchange flow naturally achieves continuously forced, shear-driven, strongly stratified turbulence at the interface, i.e. away from the solid duct boundaries (a good approximation to free shear).

Figure 1. Set-up and notation. (a) The SID experiment (see § 3.1). (b) Measured duct subvolume, dimensional coordinate system $(x^d,y^d,z^d)$ and flow variables $\rho ^d, u^d$, with three key schematic flow profiles (in grey, red and blue) at $y^d=0$ and $z^d_{max}$. We then compare these key flow profiles (c) in dimensional units; (d) after the hydraulic non-dimensionalisation in (3.1a–d); (e) after the shear-layer rescaling in (3.3a–d), yielding profiles comparable to canonical $\tanh$ shear layers (we discard the dashed line profiles outside the main shear layer $|y^s|>L_y, |z^s|>1$).

The SID experiment has been studied in detail in prior publications, and we refer the reader to these for further details about the set-up: Meyer & Linden (Reference Meyer and Linden2014) (hereafter ML14, see their § 2); Lefauve et al. (Reference Lefauve, Partridge, Zhou, Caulfield, Dalziel and Linden2018) (hereafter LPZCDL18, see their § 3); Lefauve et al. (Reference Lefauve, Partridge and Linden2019a) (hereafter LPL19, see their § 1–2); and Lefauve & Linden (Reference Lefauve and Linden2020) (hereafter LL20, see their § 2).

3.2. Hydraulic non-dimensionalisation

Like all exchange flows, we expect the flow in the SID to be forced by a mean streamwise pressure gradient of opposite directions in each layer, even when the duct is horizontal. This streamwise pressure gradient results from the expectation that the pressure is constant along the plane of neutral density $\rho =\rho _0$ and that the ends of the duct sit in reservoirs of different densities (see Lefauve Reference Lefauve2018, § 1.2.2). The resulting two-layer hydraulic flow has (dimensional) peak-to-peak velocity jump set (approximately) by the buoyancy velocity scale $\Delta U \equiv 2\sqrt {g'H}$, where $g'\equiv g\Delta \rho /\rho _0$ is the reduced gravity.

As commonly done is the hydraulic community, LPZCDL18, LPL19 and LL20 used halves of the total density difference ($\Delta \rho /2$), duct height ($H/2$) and velocity jump ($\Delta U/2$) to non-dimensionalise all variables:

(3.1a–d)\begin{equation} \rho^h \equiv \frac{\rho^d-\rho_0}{\Delta\rho/2}, \quad \boldsymbol{u}^h = \frac{\boldsymbol{u}^d}{\Delta U/2}, \quad \boldsymbol{x}^h \equiv \frac{\boldsymbol{x}^d}{H/2}, \quad t^h \equiv \frac{t^d}{H/\Delta U}, \end{equation}

where $\rho$ and $\boldsymbol {u}=(u,v,w)$ are the density and velocity fields, respectively, $\boldsymbol {x}=(x,y,z)$ is the position vector in the coordinate system defined in figure 1(b) and $t$ is the time. The superscripts ${}d$ and ${}h$ denote, respectively, dimensional and non-dimensional hydraulic variables.

This hydraulic non-dimensionalisation leads to the natural definitions of ‘input’ Reynolds and bulk Richardson numbers (i.e. depending only on parameters set by the experimentalist) that we refer to in this paper as ‘hydraulic’ Reynolds number and ‘hydraulic’ bulk Richardson numbers:

(3.2a,b)\begin{equation} Re^h \equiv \frac{\dfrac{ \Delta U}{2}\dfrac{H}{2}}{\nu} = \frac{\sqrt{g'H}H}{2\nu} = 1.42\times 10^4 \sqrt{\frac{\Delta \rho}{\rho_0}}, \qquad Ri^h_b \equiv \frac{\dfrac{g}{\rho_0}\dfrac{\Delta \rho}{2}\dfrac{H}{2}}{\left(\dfrac{\Delta U}{2}\right)^2} = \frac{1}{4}. \end{equation}

The flow in the SID is not only forced by a streamwise pressure gradient, but also by the projection of gravity $\boldsymbol {g}$ along $x$ due to the tilt angle $\theta >0$ of the duct (in this paper between $1^\circ$ and $6^\circ$, as sketched in figure 1a). These two forcing mechanisms yield a variety of flow regimes: from laminar flow with flat interface (${\rm L}$ regime); to mostly laminar Holmboe waves propagating at the interface (${\rm H}$ regime); to intermittent turbulent (${\rm I}$ regime); to fully developed turbulence (${\rm T}$ regime). These flow regimes and their transitions have been mapped in the ($\theta,Re^h$) plane and discussed extensively in ML14, LPL19 and LL20.

One of the key conclusions of these past studies of the SID experiment is that, while the dimensional peak-to-peak velocity scale of $u^d$ is primarily set as $\Delta U \equiv 2\sqrt {g'H}$ by the longitudinal pressure gradient (hydraulic scaling), the actual (measured) non-dimensional peak-to-peak magnitude of $u^h$ (in an $x$- and $t$-averaged sense) is a complicated $O(1)$ function of $Re^h$ and $\theta$.

To illustrate this point, we define three key profiles in the duct subvolume of figure 1(b): the vertical profiles of density $\rho ^d$ (in grey), the streamwise velocity $u^d$ (in red) in the vertical plane of maximum velocity (the midplane $y^d=0$), as well as the spanwise profile of $u^d$ (in blue) in the horizontal plane of maximum velocity ($z^d=z^d_{max}$). These three profiles are drawn schematically in dimensional variables in figure 1(c) and after the hydraulic non-dimensionalisation (3.1a–d) in figure 1(d).

Figure 1(d) shows that while the duct height, duct width and the magnitude of the total density jump are always $2$, the peak-to-peak velocity $\delta u$ and the height between the velocity peaks $h$ are both a priori unknown.

3.3. Shear-layer rescaling

In order to analyse our data in a non-dimensional framework quantitatively consistent with most of the literature on stratified shear layers, we define the following shear-layer rescaling, using halves of the ‘output’ (measured) velocity jump $\delta u$ and shear-layer depth $h$:

(3.3a–d)\begin{equation} \rho^s \equiv \rho^h, \quad \boldsymbol{u}^s \equiv \frac{\boldsymbol{u}^h}{\delta u/2}, \quad \boldsymbol{x}^s \equiv \frac{\boldsymbol{x}^h}{h/2}, \quad t^s \equiv \frac{t^h}{h/\delta u}, \end{equation}

where the superscripts $h$ and $s$ denote, respectively, the hydraulic non-dimensional variables defined in (3.1a–d) and the new shear-layer variables.

Figure 1(e) shows that the rescaled total velocity jump and shear-layer depth are now always 2. Since the symmetry of the flow with respect to $y^s, z^s=0$ is sometimes approximate, we further shift the $y^s,z^s$ axes to centre them such that the bounds of the shear layer are exactly $|y^s| \le L_y, |z^s| \le 1$. The total shear-layer width $2L_y$ is the smallest width in which both profiles $u^s(y^s,z^s=\pm 1)$ are at least $70\,\%$ of their extremum value (typically $2L_y \approx 3$ in our data).

We also define the velocity-to-density thickness ratio $R\equiv 1/h_\rho$, where $2h_\rho$ is the typical non-dimensional density layer thickness defined as spacing between the points at which $\rho ^s = \pm \tanh (1) = \pm 0.76$ (giving typically $R>1$ when $Pr\gg 1$).

The dashed lines in figure 1(e) denote flow outside the shear layer, where velocities decay to zero to satisfy the no-slip boundary condition at the four duct walls. In the remainder of this paper we ignore wall effects by deliberately discarding data outside the shear layer. This choice has shortcomings, e.g. it rules out the study of turbulent entrainment across the shear layer. However, it is advantageous for the analyses presented in this paper, and consistent with our focus on the forcing and dissipative mechanisms, the self-organisation and the turbulent fractions inside stratified shear layers.

This rescaling leads to the definitions of the following ‘shear’ Reynolds number and ‘shear’ bulk Richardson number:

(3.4a,b)\begin{equation} Re^s \equiv Re^h \, \frac{\delta u \, h}{4}, \quad Ri_b^s \equiv Ri_b^h \, \dfrac{(h/2)}{(\delta u/2)^2} \equiv \dfrac{h}{2(\delta u)^2}. \end{equation}

Note that our $Ri_b^s$ is sometimes called $Ri_0$ or $J$ in the literature.

In the remainder of this paper, unless specified otherwise, we use the shear-layer variables defined in (3.3a–d) and drop the superscript $s$ (except in $Re^s$ and $Ri_b^s$, for clarity).

The corresponding governing equations for momentum and density in shear-layer variables under the Boussinesq approximation are then

(3.5a)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

(3.5b)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla} p + Ri_b^s \, \rho \, (-\boldsymbol{\hat{z}} + \sin \theta \, \boldsymbol{\hat{x}}) + \frac{1}{Re^s} \boldsymbol{\nabla}^{2}\boldsymbol{u} , \end{gather}$$

(3.5c)$$\begin{gather}\frac{\partial \rho}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \rho = \frac{1}{Re^s \, Pr} \boldsymbol{\nabla}^{2} \rho, \end{gather}$$

where we assumed that $\cos \theta \approx 1$ in nearly horizontal flows (accurate to better than $1\,\%$ in this paper). We discuss boundary conditions next.

3.4. Comparison with canonical shear layers

The above rescaling makes our data (figure 1e), non-dimensional parameters (3.4) and governing equations (3.5) comparable to those found in studies of canonical stratified shear layers defined by the initial ($t=0$) profiles:

(3.6a,b)\begin{equation} u(z)={-}\tanh(z), \quad \rho(z)={-}\tanh(R z). \end{equation}

Note the minus signs, typically absent in the literature, but retained here for historical reasons and of minor significance (note that some studies prefer to use the buoyancy field, here simply equal to $-\rho$). A relatively small number of studies opt for a non-dimensionalisation based on the total (as opposed to half) velocity jump $\delta u$ and shear-layer depth $h$, making their Reynolds number four times as large as ours, and their bulk Richardson number half as large as ours. The values of $Re^s$ and $Ri_b^s$ in Appendix A have been estimated and/or converted from various studies to match our definitions consistent with the governing equations (3.5) and the canonical $\tanh$ model (3.6). Note that while most numerical studies report the initial flow profile, observational and experimental studies more often report instantaneous or mean flow profiles.

We find at least five interesting differences between our rescaled SID flows and most canonical $\tanh$ shear layers.

1. (i) Our rescaled profiles in figure 1(e) are understood as ‘mean flows’ averaged in the horizontal direction and over a long-time equilibrium, as opposed to carefully designed initial conditions.

2. (ii) Our velocity at the top and bottom boundaries of the shear layer reaches approximately $\pm 1$ (in the midplane $y=0$) and $\pm 0.8\text {--}0.9$ (when averaged in $y$ across the layer), as opposed to the more modest $\tanh (1)\approx \pm 0.76$.

3. (iii) Our vertical shear at the top and bottom boundaries of the shear layer is zero $\partial _z u (z \approx \pm 1)=0$, because of the influence of the nearby top and bottom walls, as opposed to the typical free-slip boundary conditions at $z\rightarrow \pm \infty$.

4. (iv) Our spanwise velocity gradient is non-zero at the spanwise edges of the shear layer $\partial _y u (y = \pm L_y)\neq 0$, because of the influence of the nearby sidewalls, as opposed to the typical periodic boundary conditions in $y$.

5. (v) Our long-time equilibrium is achieved by a gravitational body force along $x$ ($Ri_b^s \sin \theta \rho$) and by non-periodic boundary conditions along $x$ responsible for both a mean horizontal pressure gradient and a mean horizontal buoyancy flux (continuously replacing partially mixed fluid in the duct by unmixed fluid from the reservoirs), all of which are typically absent in canonical shear-layer simulations.

3.5. Data sets

We use 16 sets of time-resolved, volumetric data of the density and three-dimensional, three-component velocity $(u,v,w,\rho )(x,y,z,t)$ freely available online (Lefauve, Partridge & Linden Reference Lefauve, Partridge and Linden2019b). These were obtained by successive $x$–$z$ planar measurements of stereo particle image velocity (sPIV) and laser induced fluorescence (LIF) performed simultaneously in a rapid, continuous, back-and-forth scanning motion across $y$ to reconstruct successive three-dimensional volumes. In all experiments, the duct streamwise aspect ratio was $30$, the duct spanwise aspect ratio was $1$ (square) and the Prandtl number was $Pr \approx 700$ (NaNO$_3$/NaCl salt solutions with matched refractive indices). For more information on the set-up, scanning technique and postprocessing (including imposing $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} =0$ in all volumes), we refer the reader to Partridge et al. (Reference Partridge, Lefauve and Dalziel2019) (their §§ 3 and 4), LPZCDL18 (their §§ 3.3 and 3.4) and LPL19 (their § 3.1 and 3.2).

To suit the objectives of the present paper, these data sets were modified in the following two ways. First, as explained in § 3.3, we only retain data in the shear layer (by discarding the near-wall data dashed in figure 1e). The final size of each volume $(2L_x,2L_y,2L_z)$ is given in Appendix B, table 3 (in shear-layer units, where by definition $2L_z=2$), together with the total remaining number of grid points in each direction $(n_x,n_y,n_z)$, and the resulting resolution $(\Delta x, \Delta y, \Delta z) \equiv (2L_x/n_x, 2L_y/n_y, 2/n_z)$). Second, small errors in the initial levels of free surfaces in each reservoir (figure 1a) caused small barotropic (net) flow oscillations between the two reservoirs, which decayed exponentially with time. Data sets showing these early-time damped oscillations were cropped in time to keep only the later time, statistically steady part of the flow. The resulting length of each data set $L_t$ (in shear-layer advective time units) is given in Appendix B together with the total number of volumes $n_t$ and the temporal resolution $\Delta t\equiv L_t/n_t$ (or time taken to scan a volume from wall to wall, i.e. $y^h=\pm 1$).

The resulting final 16 data sets and accompanying movies have been made available online (Lefauve & Linden Reference Lefauve and Linden2022b) and their key properties are shown in table 1. One flow belongs to the laminar (${\rm L}$) regime (named L1), four flows to the Holmboe wave (${\rm H}$) regime (named H1–H4), eight flows to the intermittently turbulent (${\rm I}$) regime (named I1–I8) and three flows to the fully turbulent (${\rm T}$) regime (named T1–T3). These data sets are ordered by increasing values of the product of input parameters $\theta Re^h$, as in LPL19 (see their table 2) who showed that $\theta Re^h$ controlled the time- and volume-averaged kinetic energy dissipation and thus the transitions between flow regimes (note $\sin \theta \approx \theta$ in our nearly horizontal flows). The output parameters $\delta u, h, R, Re^s$ and $Ri^b_s$ were determined as explained in § 3.3, where the key profiles drawn in figure 1 were interpreted as $x$- and $t$-averages over the data set (i.e. over $x\in [0,2L_x]$, and $t\in [0,L_t]$). We discuss the values of these output parameters next.

Table 1. List of the 16 volumetric data sets used, with input parameters $\theta$ and $Re^h$ (note $Ri_b^h = 1/4$), sorted by increasing $\theta Re^h$, and thus by flow regime ${\rm L}$, ${\rm H}$, ${\rm I}$, ${\rm T}$ (as in LPL19 table 2). The output parameters follow the shear-layer rescaling in § 3.3 and figure 1(e). The four parameters, $\theta, Re^s, Ri^s_b, R$, are necessary and sufficient to describe the model in (3.5)–(3.6).

4.1. Output parameters

In figure 2 we plot maps of all 16 data sets of table 1 in the space of input parameters $(Re^h,\theta )$ (figure 2a) and in the space of our three independent output parameters $(Re^s,Ri^s_b,R)$ (figure 2b,c). We also show the power law regressions of the output parameters with respect to the input parameters (figure 2d, f). The search for such power laws is motivated by LPL19 and especially LL20, who demonstrated that many features of these flows (e.g. the flow regime, flow rate, and interfacial density layer thickness) could be understood and predicted theoretically as products of power laws of the input parameters (see LL20, § 5).

Figure 2. Maps of all 16 data sets of table 1 in the space of (a) input parameters $(Re^h,\theta )$ and (b,c) output parameters $(Re^s,Ri^s_b,R)$. The dashed curves in (a) are the regime transitions in $\theta Re^h =$ const. previously given by LPL19. The dashed rectangles in (b,c) highlight the fact that regimes also occupy distinct regions in the output space. (d–f) Best power law fit (least squares linear regression in log–log space) of the output parameters $(Re^s,Ri^s_b,R)$ by the input parameters $\theta,Re^h$ (fit versus actual value, the dashed line denoting equality). In power law scalings, $\theta$ is always expressed in radians.

First, we see that the bulk Richardson number $Ri_b^s$ and the velocity-to-density thickness ratio $R$ typically decrease as the Reynolds number $Re^s$ increases, and appear to reach asymptotic values in the turbulent regime (figure 2b,c). In other words, the relatively wide and uniformly sampled region of the input space (figure 2a) is mapped by the mean flow dynamics into a relatively narrow and specific region of the output space (figure 2b,c). The flow dynamics also have an inherent degree of randomness making them not generally repeatable, because we see that near-identical input parameters can be mapped into fairly different output parameters (e.g. compare the couples H2/H4, I4/I6 and I5/I8 in figure 2a and figure 2b,c). The above two observations mean that the experimentalist has only a limited (and not fully understood yet) ability to control the output parameters from the input parameters. This is likely due to the inherent chaotic nature of such high-$Re$ flows, combined with the limited control over initial conditions (procedure and care with which the duct was opened), and limited duration of the time-averaging intervals $L_t$ (especially in the ${\rm I}$ and ${\rm T}$ regimes).

Second, we see that different qualitative flow regimes ${\rm L},{\rm H}, {\rm I}, {\rm T}$ (in blue, green, yellow and red, respectively) occupy distinct regions in the $(Re^s,Ri_b^s,R)$ output space (sketched in figure 2b,c by the dashed rectangles), although they are generally better separated in $Re^s$ than in $Ri_b^s$ or $R$. This result, which implies that the different flow regimes reflect different physics, could not simply be predicted a priori from the previously known result that regimes occupy distinct regions in the $(Re^h,\theta )$ input space (sketched in figure 2a by the dashed curves of LPL19).

In the output space (figure 2b,c), the transition from stable laminar flow to regular Holmboe waves (${\rm L} \rightarrow {\rm H}$) is correlated with $Re^s \gtrsim 100-200$, $Ri^s_b \lesssim 0.6-0.8$ and $R\lesssim 12$; values that are consistent with the triggering of Holmboe instability. The transition to intermittent turbulence (${\rm H}\rightarrow {\rm I}$) is correlated with $Re^s\gtrsim 500$, $Ri^s_b\lesssim 0.2$ and $R \lesssim 7$, while the transition to sustained turbulence (${\rm I} \rightarrow {\rm T}$) is correlated with $Re^s \gtrsim 1000$, and the asymptotic values $Ri^s_b \approx 0.15$ and $R\approx 2$.

Third, we observe that the maps in the output space are not entirely consistent with the use of $\theta Re^h$ as a proxy for flow regimes and as a means to quantitatively order flows within regimes (based on their closeness to another regime), as was done in LPL19 and in our nomenclature of the data sets. For example, we see in figure 2(b,c) that I2/I3 are closer to ${\rm T}$ flows than I6/I8 are, and that T2/T3 are closer to ${\rm I}$ flow than T1 is, whereas our nomenclature suggests otherwise in both cases. We also see in figure 2b that, although the five flows I4–I8 have near-identical $\theta Re^h$, they stretch all the way from the ${\rm H}$ transition to the ${\rm T}$ transition.

Fourth, we note that vigorous turbulence can be sustained even at relatively low $Re^s\sim 1000$ due to the continuously forced nature of SID flows. Indeed, our largest value $Re^s \approx 1500$ in the ${\rm T}$ regime is a factor three to four lower than the values of $4000$–$6000$ investigated in the latest numerical simulations of stratified shear layers (Salehipour & Peltier Reference Salehipour and Peltier2015; Salehipour et al. Reference Salehipour, Caulfield and Peltier2016; Smith et al. Reference Smith, Caulfield and Taylor2021). Although much higher $Re^s \approx Re^h = O(10^4-10^5)$ can readily be achieved in the SID experiment (see LL20), they are not shown here because they remain out of reach of detailed quantitative measurements due to limitations in the spatiotemporal resolution of the scanning sPIV/LIF technique (discussed in LPL19, Appendix A).

Fifth, we study the power law regression (best fit) of the output parameters with respect to the input parameters. The scaling $Re^s \propto \theta ^{0.73} (Re^h)^{1.4}$ (figure 2d) is an excellent fit, since most symbols lie close to the dashed line (the coefficient of determination is $r^2=0.98$). This shows that $\theta$ plays a key role in setting the non-dimensional scales $\delta u, h$, and thus $Re^s$ (remembering from (3.4) that $Re^s\equiv \delta u \, h \, Re^h/4$). However, $Ri_b^s \propto \theta ^{-0.91} (Re^h)^{-0.42}$ (figure 2e) is a poorer fit ($r^2=0.73$); although $Ri_b^s$ tends to decrease with both $\theta$ and $Re^h$, the data have more variability than can be explained by a simple power law. Finally, $R \propto \theta ^{-1.14} (Re^h)^{-0.96}$ (figure 2 f) is a good fit ($r^2=0.88$), showing that the non-dimensional density layer thickness $2h_\rho =2R^{-1}$ tends to increase slightly more strongly with $\theta$ than with $Re^h$. This is consistent with the findings of LL20 (see their figures 7 and 8), who applied a similar (though higher-order) fitting to density layer thickness data obtained by shadowgraph image analysis in various duct geometries, across hundreds of experiments covering a wider range of $\theta \text{ and } Re^h$ than in the present paper.

4.2. Mean flows

We now turn to mean flows. Here, and in the remainder of this paper, we define the averages for any flow variable $\phi$ as follows:

(4.1a)$$\begin{gather} \bar{\phi}(y,z) \equiv \langle \phi \rangle_{x,t} \equiv \frac{1}{2L_t L_x} \int_0^{L_t} \int_0^{2L_x}\phi \ \, \text{d}\kern0.06em x\, \, \text{d} t, \end{gather}$$

(4.1b)$$\begin{gather}\langle \phi \rangle \equiv \langle \phi \rangle_{x,y,z,t} \equiv \frac{1}{4L_y} \int_{{-}L_y}^{L_y} \int_{{-}1}^{{-}1}\bar{\phi} \ \, \text{d} y\, \, \text{d} z, \end{gather}$$

where $\langle \phi \rangle _{i}$ denotes averaging with respect to any coordinate $i$; $\bar {\phi }$ denotes specifically $x$- and $t$-averaging (what we usually call the ‘mean’); and $\langle \phi \rangle$ denotes time- and volume-averaging. All averaging is performed using accurate trapezoidal numerical integration.

Figure 3 shows the mean streamwise velocity $\bar {u}$ and density $\bar {\rho }$ from all 16 data sets. Each panel (a–p) corresponds to a data set; the top subpanels show vertical profiles (both the midplane velocity maximum $\bar {u}(y=0)$ and $y$-averages $\langle \bar {u}\rangle _{y}$, $\langle \bar {\rho }\rangle _{y}$, across the whole shear layer $|y|\le L_y$), while the bottom subpanels show the spanwise profiles at the top and bottom edges of the shear layer $\bar {u} (z=\pm 1)$.

Figure 3. Mean flows in all 16 data sets of table 1: vertical profiles (top subpanels) and spanwise profiles (bottom subpanels), see legend for details. Legend and axes limits are identical in all respective subpanels.

First, we see that the horizontal profiles $\bar {u}(z=\pm 1)$ (black dashed and dash–dotted curves) show excellent spanwise symmetry (about the $y=0$ plane), as expected from the symmetry of the duct. In the shear-layer region plotted here ($|y|\le L_y$), where we recall that by definition velocities are at least 70 $\%$ of their extrema, we see a fairly extended flat region where $\partial _y \bar {u} \approx 0$. This region typically occupies at least $|y|\le 1$, and is slightly wider in some data sets, with no obvious dependence on flow parameters (not even on $Re^s$, surprisingly). This suggests that despite the existence of sidewalls, SID flows exhibit shear layers whose mean flows exhibit very little spanwise variation over an extent at least as large as the vertical extent ($|z|\le 1$). Closer to the spanwise edges of the shear layer, the mean flows have $\partial _y \bar {u} \neq 0$ and the resulting effects of this spanwise shear on the turbulence can in principle be investigated (which is not possible in simulations with periodic boundary condition in $y$).

Second, we see in some data sets that the vertical profiles of $\bar {\rho }$ (red solid) and $\bar {u}$ (black solid) are ‘offset’ with respect to one another, i.e. the $\bar {\rho }=0$ and $\bar {u}=0$ levels are not collocated and $\bar {\rho }\bar {u}<0$ (this is particularly visible in figure 3c,e,g,h,i,j,k).

In the Holmboe wave regime, where the density interface is sharp ($R>7$) and $\tanh$-like, this offset gives rise to asymmetric (i.e. one-sided) Holmboe waves (in H2 and H4). (For further empirical observations of this offset, see Lefauve (Reference Lefauve2018) § 3.2.2; and for visualisations and explanation of these waves in H4, see LPZDCL18.) By contrast, the absence of offset gives rise to symmetric (i.e. two-sided) Holmboe waves (in H1 and H3). (For a visualisation of these waves in H1, see LPL19 figure 3g-j.) This offset is inconsistent with the effects of gravitational forcing alone (see the term $Ri^s_b \sin \theta \rho$ in (3.5b)), and thus suggests the existence of a horizontal pressure gradient with a more complicated $z$ profile than hitherto assumed.

In the ‘weakly’ intermittent regime (I2–I6), the density interface is broader ($R\approx 2\text {--}5$) and this offset appears correlated with unequal entrainment and mixing (i.e. asymmetry) on either side of the $\bar {\rho }=0$ level. Further observation of the vertical profiles in figure 3(g–l) reveals that the density is indeed better mixed above its $0$ level, and that the density interface lies below the velocity interface. This is consistent with the fact that the measured duct volume lies nearer the end sitting in the $\rho =1$ reservoir (i.e. on the ‘left’, as sketched in figure 1(a); see LPL19, table 2 for the precise locations). Assuming that mixing occurs uniformly across the length of the duct, the bottom layer with initial density $\rho =1$ (coming from the ‘left’) has therefore travelled less, and thus presumably experienced less mixing, than the top layer with initial density $\rho =-1$ (coming from the ‘right’). This slight but important non-periodicity along $x$ is an important aspect of SID flows, which appears necessary to obtain continuously forced exchange flows in the laboratory.

In the more strongly turbulent regime (T2 and especially T3), the vertical density and velocity profiles become similar ($\bar {u}(z) \approx \bar {\rho } (z)$), and closer to $\tanh$/linear. The vertical symmetry of ${\rm T}$ flows and their lower thickness ratio $R \lesssim 2$ result from a more intense and sustained mixing than in ${\rm I}$ flows.

4.3. Reynolds-averaged balances

We now explain the quasi-steady maintenance of these mean flows $\bar {u},\bar {\rho }$ by analysing the steady Reynolds-averaged $x$-momentum and density equations,

(4.2a)$$\begin{gather} \underbrace{\ - \, \overline{\partial_x p} \, }_{\substack{\text{mean pressure} \\ \text{gradient} \, \equiv \, \varPi }} \underbrace{+ \ Ri_b^s \, \sin \theta \, \bar{\rho}}_{ \text{body force}} \ \underbrace{+ \ (Re^s)^{{-}1} (\partial_{yy} \bar{u} + \partial_{zz} \bar{u})}_{\text{molecular diffusion}} \ \underbrace{-\ \partial_y(\overline{u'v'}) - \partial_z (\overline{u'w'})}_{\text{turbulent diffusion}} \approx 0, \end{gather}$$

(4.2b)$$\begin{gather}\underbrace{\ - \, \overline{\partial_x(u \rho)} \, }_{\substack{\text{mean advective} \\ \text{buoyancy flux} \, \equiv \, \varLambda}} \, \underbrace{+ \ (Re^s \, Pr)^{{-}1} \partial_{zz} \bar{\rho}}_{\text{molecular diffusion}} \ \ \underbrace{-\ \partial_y(\overline{v'\rho'}) - \partial_z (\overline{w'\rho'})}_{\text{turbulent diffusion}} \approx 0, \end{gather}$$

where flow fluctuations are defined as $\phi ' \equiv \phi - \bar {\phi }$. We used incompressibility $\partial _x u +\partial _y v + \partial _z w = 0$ (imposed at all times) and the (good) approximations that $\bar {v},\bar {w},\partial _{yy} \bar {\rho }\approx 0$ and that mean flows are steady (i.e. $\overline {\partial _t u} \approx \overline {\partial _t \rho } \approx 0$).

The slight non-periodicity in $x$ gives rise to two previously mentioned important forcing terms: the mean streamwise pressure gradient denoted $\varPi (y,z) \equiv -\overline {\partial _x p}$, and the mean streamwise advective buoyancy flux denoted $\varLambda (y,z) \equiv -\overline {\partial _x (u \rho )}$ (continuously replacing partially mixed fluid in the duct by unmixed fluid from the reservoirs).

Figure 4 shows the vertical structure of each term in (4.2a) (figure 4a–e) and (4.2b) (figure 4 f–j) for five representative data sets spanning the ${\rm H}, {\rm I}$ and ${\rm T}$ regimes. Derivatives were computed using second-order-accurate finite differences, and we only plot the $y$-average of all terms, ignoring their (weak) spanwise structure. Note that we cannot measure directly the mean pressure gradient $\varPi$ in panels (a–e); instead we plot its indirect estimation $\varPi ^{estim}$ assuming a perfect balance of the three remaining terms in (4.2a). Similarly, although we measured the mean advective buoyancy flux as $\varLambda (z) = (2L_x)^{-1} [ \langle u \rho \rangle _{y,t} (x=0) - \langle u \rho \rangle _{y,t} (x=2L_x)]$, we also plot for comparison its indirect estimation $\varLambda ^{estim}$ assuming a perfect balance of the two remaining terms in (4.2b).

Figure 4. Vertical profiles of the $y$-averaged terms in the Reynolds-averaged (a–e) $x$-momentum balance (see (4.2a)) and ( f–j) density balance (see (4.2b)) in five representative data sets H1, H4, I2, I6, T3. Legends and axes limits are identical in all panels. The horizontal black dashed lines denote the velocity (a–e) and density ( f–j) interfaces, i.e. the location where $\langle \bar {u}\rangle _y,\langle \bar {\rho }\rangle _y=0$, respectively (see figure 3b,e,g,k,p). The colouring of quadrants indicate regions of diffusion (diff.) and antidiffusion (antidiff.) of these profiles (see e,j for the legend). ‘Unexpected’ behaviour occurs where the line colour does not the match quadrant colour.

In this two-layer exchange flow, terms in the momentum balance (4.2a) that are positive above the $\bar {u}=0$ level (thin black dashed lines in figure 4a–e) and terms that are negative below this level are both diffusive in the sense that they tend to weaken the flow in each layer and thus decrease $\bar {u}$. These two ‘diffusive quadrants’ are shaded in light blue in figure 4(a–e) and the terms that are expected to be diffusive (molecular and turbulent diffusion) have a similar light blue line colour. Vice versa, terms that have opposite values on either side of the $\bar {u}=0$ level are antidiffusive in the sense that they tend to strengthen the flow in each layer and thus increase $\bar {u}$. These two ‘antidiffusive’ quadrants and the terms expected to be antidiffusive are coloured purple. We extend this diffusive/antidiffusive distinction to the density balance (4.2b) and figure 4( f–j) using pink and maroon, respectively. As a result, unexpected behaviour occurs in regions where line and quadrant colours do not match, which is the focus of the discussion below.

First, we see that molecular (laminar) diffusion of momentum (dotted blue) and density (dotted pink) is negligible in all flows (the lines are barely distinguishable from 0), at least in the shear-layer region ($|z|\le 1$). By contrast, turbulent diffusion (dash–dotted blue and pink) is important in this region, reaching locally absolute values of order $O(0.01)$, which would be responsible for $O(1)$ changes over $O(100)$ advective time units in the absence of counter-acting mechanisms (i.e. over $O(L_t)$, the total time captured in our data sets). Turbulent diffusion behaves diffusively as expected (i.e. these lines are in the quadrant matching their colour, blue and pink, respectively), except in the Holmboe regime where these terms are strikingly antidiffusive in the vicinity of their respective $\bar {u}=0$ and $\bar {\rho }=0$ interfaces, and diffusive farther away from them (figure 4a,b, f,g). This means that the fluctuations of Holmboe waves effectively sharpen, or ‘scour’ both the velocity and density interface. This sharpening occurs symmetrically on either side of the interfaces in H1 (figure 4a, f) and asymmetrically (only above the interfaces) in H4 (figure 4b,g). This is consistent with the previously mentioned fact that H1 sustains symmetric (both upward- and downward-pointing) Holmboe waves, while H4 sustains asymmetric (upward-pointing only) Holmboe waves.

Second, the gravitational body force (solid purple) is, as expected, antidiffusive almost everywhere (i.e. sustaining $\bar {u}$), except in the regions where velocity and density interfaces are offset (figure 4b–d) as discussed in the previous section. However, an unexpected result of figure 4(a–e) is that the estimated mean pressure gradient $\varPi ^{estim}$ (dashed purple) is diffusive (i.e. adverse) almost everywhere. In ‘offset’ regions where $\bar {u} \bar {\rho }<0$, this unexpected adverse pressure gradient may provide an explanation for the sustained offset of interfaces (fluid is forced by the pressure gradient to flow against the natural direction suggested by gravitational forcing). However, in ‘regular’ regions where $\bar {u}\bar {\rho }>0$, this unexpected adverse pressure gradient is contrary to our intuition derived from horizontal ($\theta =0^\circ$) exchange flows where $\varPi$ is necessarily antidiffusive (i.e. favourable), as it is the only forcing causing the flow. In exchange flows inclined at even small angles (e.g. $\theta =1^\circ$ in H1) and thus forced by gravity, our results suggest that the particular equilibrium enforced by hydraulic control within the duct results in an adverse pressure gradient (at least throughout most of the shear layer in $z$, and away from the in-flow at either ends of the duct in $x$ where a relatively favourable pressure gradient is expected). This adverse $\varPi$ is consistent with the fact that hydraulic control enforces a relatively low velocity scaling $\Delta U\propto \sqrt {g'H}$ (inertial–hydrostatic balance), instead of the much larger $\sqrt {g'H}\sin \theta Re^s \approx \sqrt {g'H}\theta Re^s$ expected in an infinitely long or periodic tilted duct (gravitational–viscous balance), as explained in LL20 § 2.3. What is not yet understood is the underlying structure of the pressure field which would allow this combination of in-flowing favourable $\varPi$ and out-flowing adverse $\varPi$, while matching with the far-field hydrostatic pressure distribution in the reservoirs.

Third, we see that the measured value of the mean streamwise advective buoyancy flux $\varLambda$ (maroon dashed, not closing the density balance) and its estimated value $\varLambda ^{estim}$ (maroon solid, closing the density balance) are only significantly different in H1 and I2 (figure 4 f,h). In H1, turbulent antidiffusion (scouring) near the density interface requires $\varLambda ^{estim}$ to be (unexpectedly) diffusive, but direct measurement of $\varLambda$ suggests otherwise, which is not presently understood. In I2, $\varLambda ^{estim}$ apparently underestimates $\varLambda$, possibly due to limitations in the spatiotemporal resolution of these measurements (see Appendix B and quantification of this effect in LPL19, figure 12). This suggests that in some data sets we could use the measured $\varLambda$ as a proxy for turbulent diffusion of density (rather than the other way around). However, doing so would require trust in $\varLambda$ and in the exact balance of (4.2a), and we have seen above that at least one of these could be questionable (see H1).

5.1. Definitions

The gradient Richardson number $Ri_g(\boldsymbol {x},t) \equiv N^2/S^2$ is the ratio of the squared buoyancy frequency $N^2(\boldsymbol {x},t)\equiv -Ri_b^s \partial _z \rho$ to the square of the vertical shear frequency $S^2(\boldsymbol {x},t) \equiv (\partial _z u)^2$ (in non-dimensional shear-layer units, recalling that $\partial _z \bar {u}, \partial _z \bar {\rho }<0$ throughout the shear layer). It gives a pointwise measure of the stability of stratified shear flows, since stratification (high $N^2$) tends to stabilise the flow, whereas shear (high $S^2$) tends to destabilise it.

However, in order to work with more tractable (lower-dimensional and smoother) statistics, we consider instead the buoyancy frequency, shear frequency and the gradient Richardson number based on the mean flow:

(5.1a–c)\begin{equation} \overline{\overline{N}}^2(y,z) ={-}Ri_b^s \, \partial_z\bar{\rho}, \quad \overline{\overline{S}}^2(y,z) \equiv (\partial_z \bar{u})^2, \quad \overline{\overline{Ri_g}}(y,z) \equiv \frac{\overline{\overline{N}}^2}{\overline{\overline{S}}^2}\equiv{-}Ri_b^s \frac{\partial_z \bar{\rho}}{(\partial_z \bar{u})^2} . \end{equation}

Note that we use this ‘double overbar’ notation to avoid confusion with the single overbar notation implying the different quantities $\langle N^2\rangle _{x,t}$, $\langle S^2\rangle _{x,t}$, $\langle Ri_g \rangle _{x,t}$, which are noisier and not discussed here.

5.2. Vertical profiles

In figure 5(a–p) we plot the vertical structure of this ‘mean’ $\overline {\overline {Ri_g}}$ in all 16 data sets (log–lin scale). We show averages in $y$ across the shear layer (thick black line) together with the total spread across all $y$ locations (grey shading).

Figure 5. Gradient Richardson number based on the mean flow. (a–p) Vertical profiles in all 16 data sets (thick lines denote $\langle \overline {\overline {Ri_g}} \rangle _y$, while grey shadings denote the spread over the entire range $y\in [-L_y,L_y]$). Note the log scale in $\overline {\overline {Ri_g}}$ spanning three decades $0.01$–$10$. (q–t) Correlations between $\overline {\overline {N}}^2$ and $\overline {\overline {S}}^2$ in H1, I2, I6 and T3. Symbol shade denotes the absolute vertical position $|z|$. The dashed and dotted lines correspond to a ratio of $0.1$ and $1$, respectively. (u) Probability distribution functions (p.d.f.s) stacked with successive offsets of $+1$.

First, we note that the spread in $y$ is generally modest (less than an order of magnitude), especially near the interface ($|z|\lesssim 0.5$). (For a visualisation of the $y$ dependence across the whole duct cross-section in data set T3, see Partridge et al. (Reference Partridge, Lefauve and Dalziel2019) figure 7c,i.)

Second, focusing on the $y$-averages, we observe that the ${\rm L}$ and ${\rm H}$ profiles (figure 5a–e) tend to have two minima of order $0.02$–$0.1$ on either side of the sharp density interface, and a distinct hump of order $0.2$–$2$ around the interface. Overall, $\overline {\overline {Ri_g}}$ values tend to monotonically decrease with increasing forcing (i.e. from L1 to T3), except near the edges of the shear layer $z \approx \pm 1$ where large values are always expected since $\partial _z \bar {u}=0$ by definition.

Third, a clear change in structure occurs in ${\rm I}$ profiles, where the single (dromedary) hump of ${\rm L}$ and ${\rm H}$ profiles breaks into a double (camel) hump on either side of the growing interfacial layer of mixed fluid. A final change in structure occurs in the stronger ${\rm I}$ and in all ${\rm T}$ profiles, where the double hump flattens and $\overline {\overline {Ri_g}}$ becomes nearly constant at $\approx 0.1-0.2$ across most of the shear layer.

5.3. Gradient correlations

In order to understand this last observation that $\overline {\overline {Ri_g}} \rightarrow 0.1\text {--}0.2$ in the turbulent shear layer, we investigate in figure 5(q–t) correlations between the numerator $\overline {\overline {N}}^2$ (vertical axis) and the denominator $\overline {\overline {S}}^2$ (horizontal axis). We plot, for four representative data sets (H1, I2, I6 and T3), the cloud of all $n_y n_z$ data points (those visible within those axis limits), and denote the absolute $|z|$ position in shades of grey (white representing data near the edges, of lesser interest).

We see in H1 a ‘comma’-shaped cloud, having a flat low-gradients part corresponding to an asymptote in $\overline {\overline {N}}^2$, and a steep straight high-gradients part corresponding to local values $\overline {\overline {Ri_g}} \approx 0.1\text {--}1$ (see the dashed and dotted guide lines). This structure is more or less conserved in I2 (weak ${\rm I}$ regime) although lower local values $\overline {\overline {Ri_g}} \ll 0.1$ (below the dashed line) are found at midheights (grey colour), in agreement with the profile in figure 5(g). However, very low and very high density gradients disappear in I6 (strong ${\rm I}$ regime) and T3, where the cloud becomes increasingly small and compact around $\overline {\overline {N}}^2\approx 0.04\text {--}0.4$ and $\overline {\overline {S}}^2\approx 0.2\text {--}2$, while remaining tangent to the $\overline {\overline {Ri_g}} = 0.1$ scaling (dashed line).

5.4. Histograms

To complement the above observations, we plot in figure 5(u) estimates of the p.d.f. of $\overline {\overline {Ri_g}}$ in all 16 data sets, stacked vertically for visualisation purposes. These p.d.f.s are essentially histograms based on $n_y n_z$ points, normalised such that $\int _0^1 \text {p.d.f.} \, \, \text {d} \overline {\overline {Ri_g}} = 1$ (note that we ignore the large $\overline {\overline {Ri_g}}>1$ values at the edges of the shear layer).

This figure shows that the relatively broad and/or multipeaked p.d.f.s of ${\rm L}$ and ${\rm H}$ flows progressively become narrower and single-peaked in late ${\rm I}$ and ${\rm T}$ flows. Intense turbulent flows are thus characterised by mean gradient Richardson numbers largely in the range $0.1-0.2$, with a sharp peak near $0.10\text {--}0.15$ in each case.

6.1. Perturbation enstrophy fraction

We start by defining the perturbation enstrophy as

(6.1)\begin{equation} \omega'^2(\boldsymbol{x},t) \equiv \|\boldsymbol{\nabla} \times \boldsymbol{u}'\|^2, \end{equation}

where we recall that $\boldsymbol {u}' \equiv \boldsymbol {u}(\boldsymbol {x},t) - \overline {\boldsymbol {u}}(y,z)$ (as defined and used in (4.1a) and (4.2)). This measure ignores the shear associated with the mean flows (figure 3) in order to capture perturbations away from it, representative of waves or turbulence. Note that the use of our shear-layer-rescaled velocity field (implicit throughout since § 3.3) ensures that all data sets can be meaningfully compared side-by-side.

First, we plot in figure 6(a) the time- and volume-averaged $\langle \omega '^2 \rangle$ (as defined in (4.1b)) for all 16 data sets (ordered following the nomenclature of table 1 based on $\theta Re^h$). We also plot the standard deviation in time of this volume average (shown as error bars) to highlight temporal variability. The average $\langle \omega '^2 \rangle$ increases from $\approx 0$ in L1 to $\approx 0.1\text {--}0.3$ in ${\rm H}$ flows, to $\approx 0.2\text {--}0.8$ in ${\rm I}$ flows, to $\approx 0.5\text {--}1.5$ in ${\rm T}$ flows, with some overlap between regimes. This increase is not entirely monotonic; the symmetric Holmboe wave flows H1 and H3 have slightly higher values and temporal variability than the asymmetric Holmboe wave flows H2 and H4, and those values are comparable to the weaker intermittent flows I1–I4, while the stronger intermittent flows I6–I8 are comparable to the weaker turbulent flow T1. Although absolute temporal variability roughly follows a similar pattern, the ratio of standard deviation to mean (sometimes called the coefficient of variation) is fairly constant at $\approx 15\,\%\text {--}30\,\%$ in most ${\rm H}$, ${\rm I}$ and ${\rm T}$ flows, except in I5, I7, I8, T1 where it reaches $\approx 35\,\%\text {--}45\,\%$. In other words, those four flows could be considered the ‘most intermittently turbulent’, although the remaining ${\rm I}$ flows do exhibit turbulent and more quiescent events, and the remaining ${\rm T}$ flows do exhibit temporal variability in the amplitude of their turbulence. Overall, these results confirm the expectation that higher values of $\omega '^2$ and of relative temporal variability, respectively, represent higher levels of turbulence and intermittency. Therefore, both provide a quantitative basis generally consistent (but not exactly coincident) with the earlier qualitative flow regime classification.

Figure 6. Turbulent fractions based on enstrophy and density overturns in all 16 data sets. (a) Time-and volume-average of the perturbation enstrophy $\langle \omega '^2 \rangle$ (error bars show $\pm$ one standard deviation in time of the volume average). (b) cumulative distribution function (c.d.f.) of the perturbation enstrophy (using all $n_x n_y n_z n_t$ points), highlighting the threshold values 0.5 (dotted line) and 2 (dashed line). (c) Enstrophy fraction corresponding to threshold values $\omega '^2>0.5$ (small empty symbols) and $\omega '^2>2$ (large full symbols), plotted against the averages of (a). (d) Probability distribution function (p.d.f. or normalised histogram) of the density (using all $n_x n_y n_z n_t$ points), separating ${\rm L}$/${\rm H}$, ${\rm I}$ and ${\rm T}$ data for greater clarity. Individual p.d.f.s are shown in thin lines, and the average p.d.f.s for each subpanel are shown in thick lines. (e) Cumulative distribution function of the negative density gradient $-\partial _z \rho$ limited to points where $|\rho |<0.9$, highlighting the ‘overturn threshold’ value $-0.1$ (dashed line). ( f) Overturn fraction $-\partial _z \rho <-0.1$ (with $|\rho |<0.9$), plotted against the enstrophy fraction of (c). Only non-negligible fractions $>0.1\,\%$ are shown. All gradients are computed by second-order finite differences. To remove outliers caused by these gradient computations (for the purpose of this figure only) all of the $\boldsymbol {u}'$ and $\rho$ data were smoothed with a spatial filter having an isotropic three-dimensional Gaussian kernel with a modest standard deviation of one grid point and a tight window of $5\times 3 \times 5$ grid points in $x,y,z$.

Second, to go beyond averaged values, we plot in figure 6(b) the cumulative distribution function (c.d.f.) of $\omega '^2$, obtained by integration of the p.d.f.s (normalised histograms). All c.d.f.s have a similar sigmoidal shape, with an inflection point at $\approx 0.5$, but their relative position along the $\omega '^2$ axis differs widely, consistent with the pattern of generally increasing $\langle \omega '^2 \rangle$ observed in figure 6(a). Moreover, the proximity of some data sets in terms of their averages observed in figure 6(a) extends to their whole distribution; in particular H3/I2/I3/I4, H1/I5 and I7/I8/T1 have nearly identical (indistinguishable) c.d.f.s.

Third, we plot in figure 6(c) two sets of enstrophy turbulent fractions, defined as the ratio of data points above a certain threshold of $\omega '^2>0.5$ (small empty symbols) and $\omega '^2>2$ (large full symbols) corresponding to $1-\text {c.d.f.}(0.5)$ and $1-\text {c.d.f.}(2)$, respectively (these thresholds are highlighted by dotted and dashed lines in figure 6b). This enstrophy criterion is loosely based on ideas developed in Holzner et al. (Reference Holzner, Liberzon, Nikitin, Lüthi, Kinzelbach and Tsinober2008) and Krug et al. (Reference Krug, Holzner, Lüthi, Wolf, Kinzelbach and Tsinober2015) for the characterisation of the turbulent/non-turbulent interfaces, and more generally on the fact that turbulence is associated with extreme vorticity fluctuations (long ‘tail’ of the enstrophy p.d.f.s). These two sets of fractions are plotted against the average values of figure 6(a), and reveal an excellent correlation between all three measures. Focusing on the $\omega '^2>2$ fraction (a threshold value greater than any time- and volume-averages), we find that only the more energetic six data sets I6–T3 have non-negligible fractions $>1\,\%$ representative of significant turbulent events (reaching values of $\approx 20\,\%$ for T2).

6.2. Density overturn fraction

Before investigating density overturns, we plot in figure 6(d) the p.d.f.s of the full density field $\rho$, segregating the ${\rm L}$/${\rm H}$ data (top subpanel), ${\rm I}$ data (middle subpanel), and ${\rm T}$ data (bottom subpanel). Individual p.d.f.s (thin lines) and subpanel averages (thick lines) show a similar trend: ${\rm L}$/${\rm H}$ flows have a roughly bimodal distribution $|\rho | \approx 0.9-1$; ${\rm I}$ flows develop an extra middle peak flanked by two flat and $\approx 0$ intermediate plateaus; and ${\rm T}$ flows strengthen and broaden the middle peak and increase the value of the intermediate plateaus. This increasingly broad distribution of the density field (from an initially bimodal $\rho =\pm 1$ distribution in the external reservoirs) owes to the increasing intensity of mixing. Furthermore, the asymmetry of the middle peak, almost systematically between $-0.5<\rho <0$ rather than around $0$, reveals a stronger/more efficient mixing above the density interface ($\rho <0$) than below it. This is consistent with our observations on the mean density profiles in figure 3, which we explained in § 4.2 by the non-periodicity of the flow along $x$ and the asymmetrical location of our measuring volume with respect to the duct length. However, note that the c.d.f.s (not shown here) corresponding to these p.d.f.s are not mathematically equivalent to the mean density profiles $\langle \bar {\rho }\rangle _y(z)$ of figure 3; instead the c.d.f.s of $\rho$ at any given time $t$ would yield the instantaneous background density field used to calculate the background potential energy of the flow (Winters et al. Reference Winters, Lombard, Riley and D'Asaro1995), which is beyond the scope of this paper.

Since turbulent mixing is caused by a combination of large-scale stirring and small-scale diffusion, we proceed by investigating in figure 6(e) the c.d.f.s of the vertical gradients of density $-\partial _z \rho$ (the negative sign is added for convenience). (Note the use of lin–log axes in figure 6e, as opposed to the log–lin axes in figure 6b, preventing a direct comparison of the shapes of the c.d.f.s between figure 6b and figure 6e.) Because of inherent noise in the density field, aggravated by the computation of gradients, we restricted these c.d.f.s to points where $|\rho |<0.9$, i.e. where the density field was at least partially mixed. We find that ${\rm H}$ flows (particularly H3) tend to have sharper stable gradients $-\partial _z \rho \gg 1$ (their c.d.f. plateaus to unity at higher values than most ${\rm I}$ and ${\rm T}$ flows). However, ${\rm I}$ and ${\rm T}$ flows (particularly I4, I6, I8, T1–T3) tend to have much more unstable gradients $-\partial _z \rho < 0$ (left of the dashed line), which signal density overturns. This is consistent with higher levels of turbulent mixing and the earlier qualitative flow regime classification.

Finally, we plot in figure 6( f) the overturn turbulent fraction defined as the ratio of data points having $-\partial _z \rho <-0.1$ and $|\rho |<0.9$ (corresponding to $\text {c.d.f.}(-0.1)$). These thresholds were chosen to avoid noisy and spurious gradient values caused either by clearly unmixed fluid ($|\rho |> 0.9$) or by very well-mixed fluid ($|\rho |<0.9$ but $\partial _z\rho \lesssim 0$). This overturn criterion is loosely based on ideas developed in Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) for the identification of dynamically distinct regions in stratified turbulence. The overturn fraction is plotted against the enstrophy fraction of figure 6(c), and the axis limits $>0.1\,\%$ hide the least turbulent flows of lesser interest. Overall, overturn fractions tend to be fairly low ($<6\,\%$), and lower than enstrophy fractions. Moreover, we find a very good correlation between both fractions (most points follow a linear scaling), with the exception of I4, whose overturn fraction is an order of magnitude above that expected (based on its enstrophy fraction, and on the neighbouring flows I2, I5 which we recall have very similar $\theta Re^h$ values).

6.3. Flow visualisations

To delve deeper into the above observations, figure 7 offers visualisations of these turbulent fractions in four ‘case studies’ using data sets I2, I4, I7 and T2, which are representative of the four main clusters in figure 6( f). We plot three types of information. First, for each data set (highlighted by the three yellow boxes and one red box around figure 7a–l) we plot a snapshot of the underlying full enstrophy $\omega ^2 \equiv \|\boldsymbol {\nabla } \times \boldsymbol {u}\|^2$ and the simultaneous density $\rho$ in the midplane $y=0$ (for T2 only, we also plot the midplanes $z=0$ and $x=-17.6$). Second, we identify in black contours the regions exceeding the perturbation enstrophy threshold $\omega '^2>2$ and the overturn threshold $-\partial _z\rho <-0.1$ (with $|\rho |<0.9$) as discussed in figure 6. The corresponding turbulent fractions in each plane (relative area as a percentage) are displayed in the top right corner of each panel. Third, we plot in figure 7(m–p) the time series of these turbulent fractions averaged over the whole volume (recall that the time-average of these two series was shown in figure 6f). The vertical dashed lines in figure 7(m–p) denote the time of the respective snapshots in figure 7(a–l), proving that our choice of snapshots represents typical (rather than extreme) values. We recall that the mean flows corresponding to these four data sets were shown previously in figure 3(g,i,l,o). We now describe each flow in turn to highlight their salient features.

Figure 7. Snapshots and time series of turbulent fractions based on perturbation enstrophy and density overturn for (a,b,m) I2, (c,d,n) I4, (e, f,o) I7 and (g–l,p) T2. For each data set, we show a single snapshot in time of the (a,c,e,h) midplane total enstrophy $\omega ^2(y=0)$ and (b,d, f,j) midplane total density $\rho (y=0)$. For T2 only, we also show the midplanes (g,i) $x=-17.6$ and (k,l) $z=0$ (the dashed lines in (g–l) denote the location of these plane cuts). Black contours show the regions exceeding the respective turbulent thresholds, and their respective fractions in each plane are given as a percentage. For each data set we also show in (m–p) the time series of the volume-averaged fractions (the dashed lines denote the time of the snapshots in (a–l)). Time does not start at $t=0$ because all data sets were cropped to remove any early-time net flow oscillations, as explained in § 3.5. Although the determination of turbulent regions and fractions is based on gradients computed on smoothed $\boldsymbol {u}',\rho$ fields (see caption of figure 7), the underlying $\omega ^2,\rho$ snapshots plotted here in colour are not smoothed. Colour bars are identical for all panels.

6.4. ‘Flavours’ of stratified turbulence

Now that we have described in more details the turbulence in ${\rm I}$–${\rm T}$ flows, we seek to clarify its relation to non-dimensional parameters. The question is: how could the distribution of turbulent fractions in roughly four clusters in figure 6( f) be predicted from the maps of input or output parameters in figure 2(a–c)?

To mention only a few apparent paradoxes prompting this question: the turbulent fractions of I4–I8 are scattered among all four clusters despite having a nearly equal product of input parameters $\theta Re^h \approx 80$–$85$ (where $\theta$ is in radians, as in all scaling laws); T1 is much ‘less turbulent’ than T2 and T3 despite having similar input $\theta Re^h \approx 120$–$130$ and the largest output $Re^s$ of all flows; and I7 is much less turbulent than T2 despite having similar output $Re^s,Ri_b^s,R$.

The power law regressions of figure 2(d–f) demonstrated the major role of $\theta$ in the (approximate) scaling of the three output parameters $Re^s,Ri_b^s,R$. It is also natural to expect that $\theta$, as key non-dimensional parameter in the governing equations (3.5), also plays a major role in the turbulent fractions, which is not captured by $Re^s,Ri_b^s,R$ alone.

To confirm this, we plot in figure 8 the power law regression (best fit) of the turbulent fractions with respect to $\theta$ and $Re^s$, using the nine data sets with turbulent fractions $>0.1\,\%$ plotted in figure 6( f). A multivariate regression including the additional two parameters $Ri_b^s$ and $R$ was also performed, but it provided little additional predictive power, probably because these two parameters are fairly constant across all nine data sets.

Figure 8. Power law fit of turbulent fractions with $\theta$ and $Re^s$ (least squares linear regression in log–log space) for the data with $>10^{-3}=0.1\,\%$ fraction of figure 6( f). (a) Enstrophy fraction and (b) overturn fraction (fit versus measured value, dashed line denoting equality). (c) Map of data sets in the $(\theta,Re^s)$ plane and contours of the enstrophy fraction and overturn fraction fits shown in (a,b). We plot the $0.1\,\%$ contour, and contours with increments of $3\,\%$ for the enstrophy fraction and of $1\,\%$ for the overturn fraction. Data sets with $<0.1\,\%$ fraction (not used for fitting) are plotted as small open symbols.

First, we see in figure 8(a) that the enstrophy fraction follows an approximate scaling $\propto \theta ^{2.7} (Re^s)^{2.8}$, and we see in figure 8(b) that the overturn fraction follows an approximate scaling $\propto \theta ^{3.2} (Re^s)^{1.8}$. This shows that both turbulent fractions increase steeply (superlinearly) with both $\theta,Re^s$, at least in the ‘low-fraction’ ($\lesssim 20\,\%$) regions investigated here. This also highlights the fact that the enstrophy fraction scales equally strongly with $\theta$ and $Re^s$, while the overturn fraction scales more strongly with $\theta$ than $Re^s$. This latter finding is consistent with the power law regression results of figure 2(d, f) from which we deduce $R\propto \theta ^{-1.6}(Re^s)^{-0.67}$, or equivalently $h_\rho \propto \theta ^{1.6}(Re^s)^{0.67}$, i.e. mixing scales more strongly with $\theta$ than with $Re^s$.

Second, we show in figure 8(c) the contours resulting from the two turbulent fraction scalings in the $(\theta,Re^s)$ log–log plane (enstrophy fraction fit in solid dark lines, and overturn fraction fit in dash–dotted light lines), together with the location of data sets in this plane. We see that the fits correctly predict the $0.1\,\%$ fraction ‘thresholds’ since all data sets used for the fitting (full symbols) are indeed located to the right of both $0.1\,\%$ contours (the remaining data sets are shown as small empty symbols, and are located to the left of at least one $0.1\,\%$ contour). These two sets of superposed contours also illustrate the existence of different clusters (or quadrants) in figure 6( f): I4, I6 have low $Re^s$ and high $\theta$, and thus achieve relatively high overturn fractions with respect to the enstrophy fraction, and vice versa for I7. This can be further quantified by the ratio of overturn-to-enstrophy fraction which follows a scaling $\propto \theta ^{0.5} (Re^s)^{-1.1}$, thus increasing with $\theta$ and decreasing with $Re^s$, and hinting at a general change in the ‘type’ or ‘flavour’ of stratified turbulence. Recalling from LPL19 that the product $\theta Re^h$ (or $\theta Re^s$, as we will show in Part 2) approximately controls the rate of turbulent kinetic energetic dissipation, and that high-$Re^s$ flows have a wider spectral inertial subrange, a tentative explanation for this change of ‘flavour’ is that high-$\theta$, low-$Re^s$ flows dissipate comparatively more energy in large-scale (overturning) eddies, whereas low-$\theta$, high-$Re^s$ flows dissipate comparatively more energy in small-scale eddies.

Finally, we note that the above empirical trends rely on relatively poor fits (figure 8(a,b) have $r^2=0.43$ and $0.56$, respectively). Therefore, they have been used to generate qualitative insight rather than detailed quantitative predictions. They would benefit from being verified by further experimental data sets across a broader range of parameters, and from receiving a theoretical basis. This would allow us to more confidently extrapolate experimental results to a parameter space more relevant to geophysical flows, especially in Reynolds number.