Experimental properties of continuously forced, shear-driven, stratified turbulence. Part 1. Mean flows, self-organisation, turbulent fractions

Abstract

We study the experimental properties of exchange flows in a stratified inclined duct, which are simultaneously turbulent, strongly stratified by a mean vertical density gradient, driven by a mean vertical shear, and continuously forced by gravity. We focus on the ‘core’ shear layer away from the duct walls, where these flows are excellent experimentally realisable approximations of canonical hyperbolic-tangent stratified shear layers, whose forcing allows mean and turbulent properties to reach quasi-steady states. We analyse state-of-the-art data sets of the time-resolved density and velocity in three-dimensional subvolumes of the duct in 16 experiments covering a range of flow regimes (Holmboe waves, intermittent turbulence, full turbulence). In this Part 1 we first reveal the permissible regions in the multidimensional parameter space (Reynolds number, bulk Richardson number, velocity-to-density layer thickness ratio), and their link to experimentally controllable parameters. Reynolds-averaged balances then reveal the subtle momentum forcing and dissipation mechanisms in each layer, the broadening or sharpening of the density interface, and the importance of the streamwise non-periodicity of these flows. Mean flows suggest a tendency towards self-similarity of the velocity and density profiles with increasing turbulence, and gradient Richardson number statistics support prior ‘internal mixing’ theories of ‘equilibrium Richardson number’, ‘marginal stability’ and ‘self-organised criticality’. Turbulent volume fractions based on enstrophy and overturn thresholds quantify the nature of turbulence between different regimes in different regions of parameter space, while highlighting the challenges of obtaining representative statistics in spatiotemporally intermittent flows. These insights may stimulate and assist the development of numerical simulations with a higher degree of experimental realism.

1. Introduction

In this two-part study we present experimental results relevant to a wide class of geophysical flows that are simultaneously as follows.

1. (i) Turbulent, i.e. inherently three-dimensional and unsteady, possessing a range of dynamically active scales, and in which momentum and scalar diffusion occurs primarily through macroscropic fluctuations (e.g. Reynolds stresses for momentum). This is quantified by a large Reynolds number (to be defined in § 3, typically $Re \gg 10^3$) reflecting the overwhelming importance of inertial forces over viscous forces.

2. (ii) Strongly stratified, i.e. flows in which the stable density stratification (typically in the form of a relatively sharp density interface) plays a significant role, for example through interfacial waves and the energetic cost of mixing the active scalar. This is quantified by a non-negligible bulk Richardson number (to be defined in § 3, typically $Ri_b = O(0.1\text {--}1)$), reflecting the non-negligible ratio of potential to kinetic energy in the system.

3. (iii) Shear-driven, i.e. flows in which turbulent kinetic energy is primarily extracted from a large-scale, largely parallel, mean shear flow, away from solid boundaries. This configuration is implicit in the definition of the bulk Richardson number mentioned above, and excludes stratified turbulence forced by moving boundaries, internal waves and other forms of spectral forcing (common in simulations using periodic geometry).

4. (iv) Continuously forced, i.e. flow in which a continuous, steady flux of energy and unmixed fluid balance the turbulent dissipation and irreversible mixing, respectively. This allows a statistically steady state of vigorous turbulence to be sustained for long periods of time (e.g. $\gg 10^2$ advective time units), as in many flows of geophysical interest (excluding horizontal gravity currents which are inherently transient).

As is often implicit in most of the geophysically oriented literature on continuously forced, shear-driven and strongly stratified turbulence, we further reduce the scope of this paper to flows that are as follows.

1. (v) Boussinesq, i.e. in which density differences are small enough (typically $\ll 5\,\%$ of the mean density) that they only play a relevant role through the acceleration of the reduced gravity.

2. (vi) High Prandtl number, i.e. in which the ratio of momentum to scalar diffusion $Pr \equiv \nu /\kappa$ (also called the Schmidt number) is typically large, as is the case of temperature and salt stratification in water (where $Pr=7$ and $700$, respectively). As a result, the region of mean shear in which the mean-to-turbulent kinetic energy transfer occurs – commonly referred to as the shear layer – is typically thicker than the density interface and embeds it (i.e. the ratio of shear layer to density interface thickness is $R >1$).

3. (vii) Nearly horizontal, i.e. in which the normal to the mean shear flow and density interface is inclined with respect to the direction of gravity at most by a small angle (e.g. $\theta < 10^\circ$), such that the main dynamics are horizontal (thus excluding plumes and exchange flows on steep slopes).

We will derive insights from recently available, three-dimensional velocity and density experimental data on exchange flows that satisfy the above conditions (ii)–(vii) and belong to four different flow regimes (from laminar, to wavy, to intermittently turbulent and fully turbulent). While only the latter two regimes satisfy (i), the former two regimes are near the turbulent transition and thus provide valuable information.

The remainder of the paper is organised as follows. We motivate this study and explain our approach in § 2, and introduce our methodology (experiment, notation and data sets) in § 3. We will then make progress on the following sets of questions, to each of which we devote a section.

1. § 4 What are the key non-dimensional parameters ($Re,Ri_b,R$), the mean profiles, the forcing and dissipative mechanisms characterising these flows in various flow regimes? How do these compare with similar flows in other observational, experimental and numerical studies?

2. § 5 What is the distribution of the gradient Richardson number – a key non-dimensional measure of the flow stability – in various regimes? Does vigorous turbulence tends towards ‘self-organisation’, i.e. a kind of self-sustaining weakly stratified equilibrium observed in other studies?

3. § 6 How to measure quantitatively and characterise the character of intermittent or sustained turbulence using the concept of turbulent fraction with simultaneous velocity and density data? How do various data sets, spanning a range of non-dimensional parameters, compare and why?

Finally, we conclude in § 7 and distil the key insights gained for the modelling of continuously forced, shear-driven, stratified turbulence. In the companion Part 2 paper (Lefauve & Linden 2022a), we tackle the energetics, anisotropy and parameterisation challenges.

2. Context

To provide context and motivation for our study, we discuss relevant field observations, numerical simulations and laboratory experiments in §§ 2.1–2.3 (for a summary table of the most recent and data-rich studies, see Appendix A, table 2). We then show where our study fits in and explain our approach in § 2.4.

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 2019) and/or wind (Yoshida et al. 1998). 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 2013).

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. (2004) 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 1987; Peters & Bokhorst 2000; MacDonald & Horner-Devine 2008; Tedford et al. 2009; Geyer et al. 2010). 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 1988; Farmer & Armi 1988; Wesson & Gregg 1994; Macias et al. 2006).

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. 1998; van Haren et al. 2014). 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 2011; Smyth et al. 2013, 2017).

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 2000). 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 2000; Mashayek & Peltier 2012a,b, 2013). Mixing layers have complicated turbulent and mixing properties dependent on parameters such as the Reynolds, bulk Richardson and Prandtl numbers (Salehipour & Peltier 2015; Salehipour, Peltier & Mashayek 2015; Salehipour, Caulfield & Peltier 2016; Watanabe, Riley & Nagata 2017; Salehipour, Peltier & Caulfield 2018). 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. 2017b; Zhou, Taylor & Caulfield 2017a), and using a relaxation of the mean profiles to initial conditions in the stratified shear layer flow (Smith, Caulfield & Taylor 2021).

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 (2001) 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 (1971)). Odier et al. (2009), Odier, Chen & Ecke (2014) and Odier & Ecke (2017) 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 (2014) and Lefauve & Linden (2020) 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 (2019) 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. (2018) studied the three-dimensional stability properties of interfacial Holmboe waves, and Lefauve, Partridge & Linden (2019a) 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. (2019a,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. (2019a) 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. Methodology

We introduce our experimental set-up in § 3.1, the hydraulic non-dimensionalisation of variables in § 3.2 and the non-dimensional rescaling of experimental data suited to comparison with canonical stratified shear layers in §§ 3.3–3.4. Finally, we introduce our data sets in § 3.5.

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 (2014) (hereafter ML14, see their § 2); Lefauve et al. (2018) (hereafter LPZCDL18, see their § 3); Lefauve et al. (2019a) (hereafter LPL19, see their § 1–2); and Lefauve & Linden (2020) (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 2018, § 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 2019b). 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. (2019) (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 2022b) 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. Flow parameters and Reynolds averages

In this section we further characterise our data sets with three key pieces of information: the output flow parameters in § 4.1; the mean flow profiles in § 4.2; and the Reynolds-averaged balances sustaining these mean flows in § 4.3.

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 2015; Salehipour et al. 2016; Smith et al. 2021). 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 (2018) § 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. Gradient Richardson number and self-organisation

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

5.5. Discussion

Our $\overline {\overline {Ri_g}}(z)$ data in ${\rm H}$/${\rm I}$ flows bear similarities to the deep-sill ocean overflow data of van Haren et al. (2014), especially to their figure 2(b). They reported long trains of Kelvin–Helmholtz overturning billows in a sustained stratified shear flow with intermittent levels of dissipation (see Appendix A for their $Re^s,Ri_b^s$ values).

Our $\overline {\overline {Ri_g}}(z)$ data in ${\rm T}$ flows are also consistent with the growing body of evidence on self-organisation in turbulent stratified shear flows subject to ‘internal mixing’, as opposed to ‘external mixing’ imposed by boundary forcing external to the shear layer (Turner 1973). The evidence suggests that a self-similar equilibrium adjustment of $\bar {u},\bar {\rho }$ occurs such that the gradient Richardson number based on the mean flows is approximately uniform across the shear layer.

This ‘equilibrium Richardson number’ hypothesis dates back at least to Turner (1973) (see his § 10.2), who quoted equilibrium values in the literature in the range $\overline {\overline {Ri_g}}(z) \approx Ri_e = 0.06\text {--}0.3$. This hypothesis is also supported by the Monin–Obhukov similarity theory, which assumes a constant buoyancy flux and derives self-similar $\bar {u},\bar {\rho }$ far enough away from any solid boundary (see Turner (1973), § 5.1), a regime verified numerically in stratified plane Couette flows (Deusebio, Caulfield & Taylor 2015; Zhou et al. 2017a).

A related ‘marginal instability’ hypothesis was also formulated in Thorpe & Liu (2009) that turbulence maintains itself on the edge of instability flagged by the linear Miles–Howard criterion of $\overline {\overline {Ri_g}}=0.25$. This was supported by the Pacific equatorial undercurrent data and calculations of Smyth & Moum (2013) (see the p.d.f. in their figure 2). A further (related) ‘self-organised criticality’ hypothesis was put forward by Salehipour et al. (2018) that strongly stratified Holmboe wave turbulence is continuously attracted to a critical value of $\overline {\overline {Ri_g}}=0.25$ (see the p.d.f. in their figure 5), making a connection to the scale-invariant energy ‘avalanches’ in the original sandpile toy model of Bak, Tang & Wiesenfeld (1988).

Comparing our $\overline {\overline {Ri_g}}$ data in flows I6–T3 with $(Re^s,Ri_b^s,R,Pr)\approx (10^3, 0.15, 2, 700)$ to the canonical stratified shear layer DNSs of Salehipour et al. (2018) (their figure 13), we find that our ‘peak’ value $Ri_e \approx 0.10\text {--}0.15$ is lower than their $Ri_e \approx 0.20$ found in ‘critical Holmboe wave turbulence’ with $(Re^s,Ri_b^s,R,Pr)=(6000, 0.16, 10, 8)$, but comparable to their $Ri_e \approx 0.10$ found in ‘subcritical Kelvin–Helmholtz turbulence’ (with much lower $Ri_b^s=0.04, R=1$). However, we note that their flow is a ‘rundown’ from an initial condition, not forced as in our experiments. Comparing with data in gravity currents forced by a positive $\theta =10^\circ$ slope, our value is compatible with $Ri_e \approx 0.1$ in the experiments of Krug et al. (2015) (their figure 8b with $(Re^s,Ri_b^s,Pr) \approx (4000, 0.30, 700)$), but higher than $Ri_e \approx 0.07$ in the DNSs of van Reeuwijk, Holzner & Caulfield (2019) (see their figure 3a with $(Re^s,Ri_b^s,Pr) \approx (4000, 0.10, 1)$).

Our data thus appear to support the recent findings of Salehipour & Peltier (2019) that ‘Kelvin–Helmholtz’ and ‘Holmboe’ turbulence share key similarities that can be ‘learned’ from large DNS data sets by deep learning methods. This provides an implicit demonstration that, once fully developed, stratified turbulence appears relatively independent of the mechanistic route taken to get to it.

6. Turbulent fractions

In this section we seek to characterise the distinction between flow regimes in more quantitative ways than done hitherto in ML14, LPL19 and LL20. We introduce the concept of turbulent fractions, i.e. the ratio of spatial regions that are ‘turbulent’ with respect to two criteria, derived from our simultaneous measurements of the density field and of the three-dimensional, three-component velocity field. We first consider a criterion based on perturbation enstrophy in § 6.1, and then a criterion based on the overturning of the density field in § 6.2. We then plot and discuss flow visualisations of these fractions in § 6.3, and the dependence on non-dimensional parameters in § 6.4.

Note that most of our velocity and (especially) our density data sets are under-resolved in the sense that they do not capture the smallest dynamically active length scales. This under-resolution undoubtedly affects the statistics below – especially the extreme perturbations – in ways that are difficult to quantify in this Part 1. In Part 2 we address this by undertaking a more detailed study of turbulent statistics, length scales, under-resolution and other measurements artefacts.

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. (2008) and Krug et al. (2015) 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. 1995), 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. (2016) 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.3.1. I2 flow

First, we recall that I2 corresponds to the ‘bottom left quadrant’ of figure 6( f) (like I1, I3, I5), representing intermittent flows with the lowest enstrophy and overturn fractions ${\sim}0.1\,\%$. The enstrophy field (figure 7a) roughly exhibits two sets of vertically stacked, quasi-periodic ‘tilde-shaped’ structures (primarily due to $\partial _z u$) coinciding with the top and bottom edges of a partially mixed layer in the density field (figure 7b). These structures are spatially only weakly ‘turbulent’ in the sense that the enstrophy only deviates significantly from its long-time mean in a few small regions (in this plane a typical $1.1\,\%$), either due to the local weakening or strengthening of the ‘core’ shear or to the shedding of top and bottom ‘filaments’, occasionally coinciding with limited density overturns (in this plane a typical $0.2\,\%$). This weak turbulence is also temporally intermittent (alternating with quiescent periods), as evidenced by the time series in figure 7(m).

The enstrophy structures (representative of weaker ${\rm I}$ flows) can be described as more disorganised and intermittent cousins of the longer-lived tilde-shaped vorticity structures previously described in H4 by Lefauve et al. (2018) (named ‘confined Holmboe waves’) and more generally found in H1–H4. The density field of early ${\rm I}$ flows, as compared with ${\rm H}$ flows, also typically features a thicker layer of mixed fluid and interfacial waves of larger amplitude more likely to overturn.

6.3.2. I4 flow

Second, we recall that I4 is the only flow in the ‘top left quadrant’ of figure 6( f) having low enstrophy fraction $\sim 0.2\,\%$ but medium overturn fraction $\sim 2\,\%$. The enstrophy field of I4 (figure 7c) is generally of lower amplitude than that of I2, and regions exceeding the turbulent threshold remain very limited. The density field of I4 (figure 7d), exhibits a thicker intermediate mixed layer than that of I2, but with weaker gradients ($|\partial _z \rho |\approx 0$) causing widespread (but weak) overturns (in this plane $6.0\,\%$). Furthermore, I4 largely lacks the large-amplitude interfacial waves found in I2 on either edges of the mixed layer, which would normally be associated with perturbation enstrophy (through baroclinic torque), and which are ultimately required to mix the density field by entrainment.

The time series in figure 7(n) gives a clue to explain the apparent paradox of how ‘so much’ mixing (here overturn fraction) could be achieved with relatively ‘so little’ stretching or rotation (here enstrophy fraction). Until $t\approx 470$ (the time at which the snapshots are shown), there is indeed very little correlation between both fractions; the overturn fraction undergoes large oscillations while the enstrophy fraction remains close to zero.

Combining all this evidence on I4, we conclude that the majority of the mixing in I4 likely occurred in vigorously turbulent regions (large enstrophy and overturn fractions) located outside of the measurement volume, and that mixed fluid was subsequently advected into the more quiescent measurement volume (note that the measurement volume spans only $13\,\%$ of the duct length along $x$). This is consistent with our prior (unpublished) shadowgraphs observations of ${\rm I}$ flows along the whole length of the duct, which occasionally showed strong spatial intermittency, i.e. the coexistence and alternation of quiescent and vigorously turbulent pockets along $x$.

6.3.3. I7 flow

Third, I7 represents the intermediate flows in figure 6( f) (like I8, T1) having medium enstrophy and overturn fractions ${\sim}1\,\%$. The enstrophy field of I7 (figure 7e) exhibits similar sets of vertically stacked tilde-shaped structures to that of I2, but these are more disorganised and likely to break off and locally exceed the $\omega '^2>2$ threshold (in this plane $4.8\,\%$). The dynamics of these structures has been described in some qualitative detail in Lefauve (2018) § 3.3.1, using planar ($y=0$) two-dimensional two-component PIV/LIF measurements at high temporal resolution (see his figure 3.13). Essentially, a turbulent ‘event’ is typically initiated by a small defect in the lower (sharper) density interface, which grows and causes the interface to roll up. The corresponding ‘single’ tilde-shaped vorticity ($\partial _z u$) structure (typical of ${\rm H}$ flows) is then stretched by the mean shear, until it eventually splits into two smaller vertically stacked structures. These structures are in turn stretched and split to create vorticity at finer scales, until the flow is clearly ‘turbulent’. The corresponding density interfaces undergo successive stretching, ejection of fluid blobs and creation of thin filaments, which promote significant mixing, without large overturns (in this plane only $1.3\,\%$).

The time series in figure 7(o) highlights the intermittent character of such events, and confirms that the enstrophy and density snapshots were chosen towards the end of a turbulent event (see the vertical dashed line), after most of the initial stretching and splitting. Furthermore, the good correlation between the enstrophy and overturn fractions is consistent with the local dynamics summarised above (as opposed to the time series of I4 for $t<450$, which required us to invoke advection of mixed fluid from outside the volume).

6.3.4. T2 flow

Fourth, T2 corresponds to the ‘top right quadrant’ in figure 6( f) (like I6, T3) having the highest enstrophy and overturn fraction ${\sim}1\,\%$–$10\,\%$. The enstrophy field of T2 (figure 7h) has generally much higher values than that of I7 (even locally exceeding the colour bar limits) and thus higher enstrophy fraction ($30\,\%$ in the $y=0$ plane). The enstrophy field of sustained turbulent flows such as T2 is also much more disorganised than that of I7; tilde-shaped structures are barely detectable among the smaller-scale transient structures. The density field (figure 7j) has higher amplitude and more frequent roll-ups resulting in higher overturn fraction ($7.3\,\%$ in the $y=0$ plane) and in stronger mixing. The time series in figure 7(p) show that high turbulent fractions ($\gg 1\,\%$) are sustained, despite some unsteadiness, and that both fractions are somewhat correlated.

Finally, ${\rm T}$ flow structures are highly three-dimensional, as evidenced by the cross-sectional $y$–$z$ cut (figure 7g,i) and horizontal $x$–$y$ cut (figure 7k,l). A consequence of this three-dimensionality is the high variability of turbulent fractions and the lack of correlation between them in individual planes. For example, the enstrophy fraction is $19\,\%$ in figure 7(g) but $30\,\%$ in figure 7(h), while the overturn fraction is $3.6\,\%$ in figure 7(i) but $10\,\%$ in figure 7(l). This highlights the importance of three-dimensional, simultaneous data for the study of stratified turbulence in the laboratory.

We refer the reader interested in further three-dimensional visualisations of ${\rm T}$ flows to Partridge et al. (2019). Their figures 8 and 9 show a snapshot of $u,v,w,\omega ^2,\rho$ in three perpendicular cuts for flow T3, and include the whole duct cross-section $(y^h,z^h) \in [-1,1]^2$ (i.e. not only the shear layer $(y,z) \in [-1,1]\times [-L_y,L_y]$ as in the present paper).

We also refer the reader to movies of all datasets made available in Lefauve & Linden (2022b).

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.

7. Conclusions

In this Part 1 we presented some key ‘basic’ properties of continuously forced, shear-driven, stratified turbulence generated by exchange flow in a square duct inclined at a small angle $\theta$ (SID experiment). We analysed 16 data sets of the simultaneous density field and three-component velocity field in a three-dimensional subvolume of the duct, spanning a range of non-dimensional parameters and flow regimes. In § 3 we adopted a convenient shear-layer non-dimensional framework to focus on the core of the flow (discarding near-wall data), and to consistently define the effective ‘shear-layer’ Reynolds number $Re^s$, bulk Richardson number $Ri_b^s$ and interface thickness ratio $R$. This allowed for easy comparison of flow profiles and statistics across all data sets, as well as with other results on stratified shear layers to support the three-pronged (observational, numerical, experimental) effort outlined in § 2. Below we summarise the progress made on the three sets of questions raised in the end of § 1.

In § 4 we described the non-trivial mapping of SID flows from the space of input parameters $\theta,Re^h,Ri_b^h$ (set by the experimentalist or the numericist) to the space of output parameters $Re^s,Ri_b^s,R$ (set by the internal flow dynamics). We also highlighted that our flows sustain turbulence at much lower $Re^s$ than in unforced stratified shear layer simulations due to the forcing tilt angle $\theta$ and the continuous advection of unmixed fluid from the external reservoirs into the duct. Next, we investigated vertical profiles of the mean flows and of the Reynolds-averaged equations sustaining them. In particular, we found that Holmboe wave fluctuations actively sharpen (or ‘scour’) the density interface on which they rely, whereas intermittently turbulent and fully turbulent flows actively broaden it. We also discovered the wide-reaching influence of the large-scale streamwise inhomogeneity of the flow (or non-periodicity in the $x$ direction). Some of its effects were readily understood, e.g. the sharpening of the density interface by advection of unmixed fluid, or the asymmetric entrainment and mixing on either side of the interface along the duct. However, some of its effects remain unexplained, e.g. the density and velocity midpoints being substantially offset in some flows, or the role of the mean hydrostatic pressure gradient which decelerates the flow at $\theta >0$. Numerical simulations resolving the pressure field in the whole geometry (duct and external reservoirs) would help elucidate this question, which might be generic to two-layer, hydraulically controlled exchange flows inclined at an angle $\theta >0$.

In § 5 we showed that the vertical profiles of the mean gradient Richardson number $\overline {\overline {Ri_g}}(z)$ smoothly evolved from a single-hump structure due to the strongly stratified interface of Holmboe flows, to a double-hump structure due to increased mixing in intermittent flows; and finally to a broad plateau $\overline {\overline {Ri_g}}(z)\approx 0.1-0.2$ across most of the shear layer in fully turbulent flows. As the turbulent intensity increases within the shear layer, we showed that the mean density gradient (buoyancy frequency) and velocity gradients (shear frequency) become tightly linked by a single, near-constant gradient Richardson number, and the probability distribution function of $\overline {\overline {Ri_g}}$ becomes narrowly peaked around $\approx 0.15$. Our data are consistent with prior theories of ‘equilibrium Richardson number’, ‘marginal stability’ or ‘self-organised criticality’, and thus provide further evidence that continuously forced, shear-driven, stratified turbulence in the SID tends to self-organise in such a distinctive ‘internal mixing’ equilibrium. However, the precise value of this equilibrium Richardson number differs across the observational, numerical and experimental studies cited, and thus remains an open question.

In § 6 we quantified the differences between flow regimes by analysing their enstrophy, density and density gradient statistics. We defined two distinct and complementary turbulent fractions as the relative flow volume exceeding a threshold in perturbation enstrophy, or experiencing density overturning. This divided intermittently and fully turbulent flows into roughly four clusters based on their location in this enstrophy–overturn turbulent fraction space, and we investigated these differences using spatial and temporal visualisations of representative flows. This revealed two particular challenges in extracting converged turbulent statistics from our experimental data, acquired over a finite, inhomogeneous subvolume of the duct and over a finite time period. First, intermittently turbulent flows show cycles with various periods, some exceeding a hundred advective time units (of the order of our recording time). Second, well-mixed turbulent fluid can be suddenly advected into our subvolume causing spikes in turbulent fraction, not due to internal dynamics, but rather to the advection of spatially intermittent turbulent patches (typically along $x$, but possibly along $y$ and $z$ too since we excluded near-wall data). Despite these challenges, approximate scaling relations between the turbulent fractions and the two key non-dimensional parameters $\theta,Re^s$ suggest that turbulence at high $\theta$ (though here $\theta <10^\circ$) and low $Re^s$ is subject to larger-scale overturning and mixing but less extreme small-scale enstrophy compared to turbulence at low $\theta$ and high $Re^s$.