Buoyancy-driven exchange flows in inclined ducts

TL;DR

This study systematically investigates buoyancy-driven exchange flows in inclined ducts, analyzing how flow regimes, mass flux, and interfacial thickness depend on five key parameters through extensive experiments.

Contribution

It provides the first comprehensive experimental dataset exploring the combined effects of five parameters on stratified exchange flows in inclined ducts.

Findings

Identified scaling laws for flow regime, flux, and interfacial thickness.

Discovered complex dependencies of flow characteristics on all parameters.

Revealed the influence of an additional sixth parameter on flow behavior.

Abstract

Buoyancy-driven exchange flows arise in the natural and built environment wherever bodies of fluids at different densities are connected by a narrow constriction. In this paper we study these flows in the laboratory using the canonical stratified inclined duct experiment, which sustains an exchange flow in an inclined duct of rectangular cross-section over long time periods (Meyer & Linden, J. Fluid Mech., vol. 753, 2014). We study the behaviour of these sustained stratified shear flows by focusing on three dependent variables of particular interest: the qualitative flow regime (laminar, wavy, intermittently turbulent, or fully turbulent), the mass flux (net transport of buoyancy between reservoirs), and the interfacial thickness (thickness of the layer of intermediate density between the two counter-flowing layers). Dimensional analysis reveals five non-dimensional independent input parameters: the duct aspect ratios in the longitudinal direction $A$ and spanwise direction $B$, the tilt angle $\theta$, the Reynolds number $Re$ (based on the initial buoyancy difference driving the flow), and the Prandtl number $Pr$ (we consider both salt and temperature stratifications). After reviewing the literature and open questions on the scaling of regimes, mass flux, and interfacial thickness with $A,B,\theta,Re,Pr$, we present the first extensive, unified set of experimental data where we varied systematically all five input parameters and measured all three output variables with the same methodology. Our results in the $(\theta,Re)$ plane for five sets of $(A,B,Pr)$ reveal a variety of scaling laws, and a non-trivial dependence of all three variables on all five parameters, in addition to a sixth elusive parameter. We further develop three classes of candidate models to explain the observed scaling laws: ...

Tables

Table 1. Table 1 : The five data sets used in this paper, using four duct geometries (abbreviated LSID, HSID, mSID, tSID) with different dimensional heights H 𝐻 H , lengths L = A ​ H 𝐿 𝐴 𝐻 L=AH and widths W = B ​ H 𝑊 𝐵 𝐻 W=BH , and two types of stratification (salt and temperature). We emphasise in bold the resulting differences in the ‘fixed’ non-dimensional parameters A , B , P ​ r 𝐴 𝐵 𝑃 𝑟 A,B,Pr with respect to the ‘control’ geometry (top row). We also emphasise the difference in H 𝐻 H between LSID and mSID, to test whether or not H 𝐻 H plays a role other than through the non-dimensional parameters A , B , R ​ e 𝐴 𝐵 𝑅 𝑒 A,B,Re . We also list the range of θ , R ​ e 𝜃 𝑅 𝑒 \theta,Re explored, and the number of regime, Q m subscript 𝑄 𝑚 Q_{m} and δ 𝛿 \delta data points obtained in the ( θ , R ​ e ) 𝜃 𝑅 𝑒 (\theta,Re) plane. Some of these data have been published or discussed in some form in ML14 (denoted by ∗ ) and LPL19 (denoted by † ) and are reused here with their permission for further analysis. Measurements of Q m subscript 𝑄 𝑚 Q_{m} and δ 𝛿 \delta were not practical with heat stratification (hence the - symbol, see text for more details). Total: 886 individual experiments and 1545 data points.

	Duct scale		Fixed params.			Varied params.		Number of data points
Name	$H$	Cross-section	$A$	$B$	$Pr$	$\theta {(}^{\circ })$	$Re(\times {10}^3)$	regime	$Q_m$	$\delta$
LSID	100		30	1	700	$[-1,4]$	$[2,20]$	173∗	20∗	115
HSID	100		15	1	700	$[0,4]$	$[1,20]$	74∗	34∗	58
mSID	45		30	1	700	$[-1,6]$	$[0.3,6]$	360†	162†	91
tSID	90		15	$1/4$	700	$[-1,3]$	$[3,15]$	131	92	87
mSID ​ ​	45		30	1	7	$[0,10]$	$[0.3,1.5]$	148	-	-

Equations

$$Q\equiv ⟨|u|{⟩}_{x,y,z,t},$$

$$Q_m\equiv \frac{2}{\Delta \rho }⟨(\rho -{\rho }_0)u{⟩}_{x,y,z,t},$$

$$A\equiv \frac{L}{H}\text{and}B\equiv \frac{W}{H},$$

$$Re\equiv \frac{\sqrt{g^{\prime }H}H}{2\nu }=\frac{\sqrt{g{H}^3}}{2\nu }\sqrt{\frac{\Delta \rho }{{\rho }_0}}.$$

$$Pr\equiv \frac{\nu }{\kappa }.$$

$$\frac{\Delta U}{2}=\sqrt{g^{\prime }H}{f}_{\Delta U}(A,B,\theta ,Re,Pr).$$

$$\underset{\mathrm{inertial}(\mathrm{I})}{\underbrace{\text{u}\cdot \nabla u}}=\underset{\mathrm{hydrostatic}(\mathrm{H})}{\underbrace{-(1/{\rho }_0){\text{p}}_{x}p}}+\underset{\mathrm{gravitational}(\mathrm{G})}{\underbrace{g\sin \theta (\rho -{\rho }_0)/{\rho }_0}}+\underset{\mathrm{viscous}(\mathrm{V})}{\underbrace{\nu {\nabla }^{2}u}},$$

$$\underset{\mathrm{I}}{\underbrace{(\Delta U{)}^2\tilde{\text{u}}\cdot \tilde{\nabla }\tilde{u}}}\sim \underset{\mathrm{H}}{\underbrace{(g^{\prime }H\cos \theta )\tilde{z}}}+\underset{\mathrm{G}}{\underbrace{(g^{\prime }L\sin \theta )\tilde{\rho }}}+\underset{\mathrm{V}}{\underbrace{(\nu \Delta U{\ell }^{-2}L)\widetilde{{\nabla }^2}\tilde{u}}},$$

$$0<\theta \ll \alpha$$

$$\alpha \equiv {\tan }^{-1}(A^{-1})$$

$$Re\gg A$$

$$G^2=4\frac{⟨u^2{⟩}_{x,y,z,t}}{\sqrt{g^{\prime }H}}=1⟹⟨|\tilde{u}|{⟩}_{x,y,z,t}=\tilde{Q}=\frac{1}{2}.$$

$$Gr\equiv \frac{g^{\prime }L\sin \theta }{(\nu /H{)}^2}\approx 4A\theta R{e}^2,$$

$$⟨{\mathsf{s}}_{ij}{\mathsf{s}}_{ij}{⟩}_{x,y,z,t}\approx \frac{1}{16}\theta Re,$$

$$Ri_{G}^{-1}\equiv\bigg{(}\frac{\sqrt{(4/9)\cos\theta+A\sin\theta}}{(1/2)\sqrt{\cos\theta}}\bigg{)}^{2}=\frac{16}{9}+4A\tan\theta=\frac{16}{9}+4\frac{\tan\theta}{\tan\alpha}.$$

$$\cos \theta \approx 1\text{and}\sin \theta \approx \theta ,$$

$$Q_m(\theta ,Re)$$

$$⟨{\mathsf{s}}_{3d}^2{⟩}_{x,y,z,t}\sim \theta Re\text{for}Re<50A,$$

$$⟨{\mathsf{s}}_{3d}^2{⟩}_{x,y,z,t}\sim \theta Re-A^{-1/2}R{e}^{1/2}\text{for}Re\gg 50A,$$

$${\eta }^{\prime }(x)={\eta }^{\prime }(\eta ,Q,\theta ,Re,f_Z,f_Y,f_I)$$

$${\eta }^{\prime }(0)=\frac{\theta Re-2QF}{Re(1-4Q^2)}\text{where}F\equiv f_Z(1+2r_Y+8r_I),$$

$$-A^{-1}<{\eta }^{\prime }(0)<0,$$

$$\theta Re<2QF<(1+b)\theta Re\text{where}b(A,\theta ,Q)\equiv \frac{1-4Q^2}{A\theta }$$

$$K_I\sim \theta Re\text{for}Re<50A.$$

$$K_I\sim A^{1/2}\theta R{e}^{1/2}\text{for}Re\gg 50A.$$

$$\mathcal{M}=\frac{{\Phi }^d}{D+{\Phi }^d}=\frac{{\Phi }_{P_B}^{\text{adv}}}{{\Phi }_{P_A}^{\text{adv}}}.$$

$$\Phi^{\textrm{adv}}_{P_{B}}=\frac{1}{4A}\langle z\rho u\rangle_{z}|_{L-R}=\frac{2}{4A\delta}\int^{\delta/2}_{-\delta/2}z\big{(}z+\frac{\delta}{2}\big{)}^{2}\,\mbox{d}z=\frac{\delta^{3}}{24A},$$

$$\delta \approx (3A\theta Q_m\mathcal{M}{)}^{1/3}.$$

$$\delta\rightarrow 0.5\Big{(}\frac{\theta}{\alpha}\Big{)}^{1/3}.$$

$${\tilde{Q}}_{m,1}=\frac{-({\rho }_1^f-{\rho }_1^i)V_1}{\Delta \rho (H^2/2)\sqrt{g^{\prime }H}T}\text{and}{\tilde{Q}}_{m,2}=\frac{({\rho }_2^f-{\rho }_2^i)V_2}{\Delta \rho (H^2/2)\sqrt{g^{\prime }H}T},$$

Full text

1 Introduction

Buoyancy-driven exchange flows naturally arise where relatively large bodies of fluid have different densities on either side of a relatively narrow constriction. In a gravitational field, this difference in buoyancy, usually in the horizontal direction, results in a horizontal hydrostatic pressure gradient along the constriction, of opposite sign above and below a ‘neutral level’, a height at which the pressures on either side of the constriction are equal. This pressure gradient drives a counter-flow through the constriction, in which fluid from the negatively-buoyant reservoir flows below the neutral level towards the positively-buoyant reservoir, and conversely, with equal magnitude. Such buoyancy-driven exchange flows result in little to no net volume transport, but crucially, in a net buoyancy transport between the reservoirs which tends to homogenise buoyancy differences in the system (i.e. towards equilibrium). In addition, irreversible mixing often occurs across the interface between the two counter-flowing layers of fluid, creating an intermediate layer of partially mixed fluid, and partially reducing the buoyancy transport. The net transport and mixing of the active scalar field (e.g. heat, salt, or other solutes) and of other potential passive scalars fields having different concentrations in either reservoirs (e.g. pollutants or nutrients) have a wide range of consequences of interest. For this reason, the study of buoyancy-driven exchange flows has a rich history. (The primary role of buoyancy being implicit throughout the paper, we will simply refer to these flows as ‘exchange flows’.)

Aristotle offered the first recorded explanation of the movement of salty water within the Mediterranean Sea (Deacon, 1971, pp. 8-9). Ever since, exchange flows through the straits of Gibraltar and the Bosphorus have driven much speculation and research, due to their crucial roles in the water and salt balances of the Mediterranean Sea, countering its evaporation by net volume transport and allowing its very existence (as first demonstrated experimentally by Marsigli in the 1680s (Deacon, 1971, Chap. 7)). More recently, it has been recognised that nutrient transport from the Atlantic partially supported primary production in Mediterranean ecosystems (Estrada, 1996). The quantification, modelling, and discussion of the past and current impact of exchange flows in straits, estuaries, or between lakes continues to generate a vast literature (as an example, the dozens of papers with titles containing ‘exchange’ and ‘Gibraltar’ have attracted over 1,000 citations in total.)

Exchange flows of gases also have a great variety of perhaps even more tangible and ancient applications to society in the ‘natural ventilation’ of buildings (Linden, 1999). It would be surprising indeed if some ice-age prehistoric Homo Sapiens did not ponder the inflow of cold outside air and the outflow of heat or fire combustion products when choosing a cave suitable for living. More recently, engineering problems of air flow through open doorways or ventilation ducts, or the escape of gases from ruptured industrial pipes, have stimulated further research.

More fundamentally, exchange flows are stably-stratified shear flows, a canonical class of flows widely used in the mathematical study of stratified turbulence, dating back at least to (Reynolds, 1883, § 12) and Taylor (1931). Multi-layered stratified shear flows have complex hydrodynamic stability and turbulent mixing properties (Caulfield, 1994; Peltier & Caulfield, 2003). The straightforward and steady forcing of exchange flows make them ideal laboratory stratified shear flows because of the ability to sustain, over long time periods, high levels of turbulent intensity and mixing representative of large-scale natural flows.

The aim of this paper is to carry out a thorough review and exploratory study of buoyancy-driven exchange flows in inclined ducts. To do this, we will focus on the behaviour of three key variables:

(i)    

the qualitative flow regime (e.g. laminar, wavy, intermittently or fully turbulent);

(ii)   

the mean buoyancy transport;

(iii)  

the mean thickness of any potential interfacial mixing layer.

The above three variables are particularly relevant in applications to predict exchange rates of active or passive scalars (e.g. salt, heat, pollutants, nutrients) between two different fluid bodies (e.g. rooms in a building, seas or lakes on either sides of a strait).

However, our primary motivation in this paper is to contribute to a larger research effort into the fundamental properties of turbulence in sustained stratified shear flows of geophysical relevance. The above three variables have thus been chosen for their particular ability to be readily captured by simple laboratory techniques while encapsulating several key flow features that are currently the subject of active research, such as: interfacial ‘Holmboe’ waves (Salehipour et al., 2016; Lefauve et al., 2018); spatio-temporal turbulent intermittency (de Bruyn Kops, 2015; Portwood et al., 2016; Taylor et al., 2016); and layering and mixing (Salehipour & Peltier, 2015; Zhou et al., 2017; Lucas et al., 2017; Salehipour et al., 2018).

To achieve this aim, the remainder of the paper is organised as follows. In § 2 we introduce a canonical experiment ideally suited to study the rich dynamics of exchange flows, and analyse the a priori importance of its non-dimensional input parameters. In § 3 we review the relevant experimental literature on the scaling laws of our variables of interest and highlight the limitations of previous research in order to motivate our study. In § 4 we present our experimental results and scaling laws. In § 5 we explain some of these results with a variety of models, and we conclude in § 6.

2 The experiment

2.1 Setup and notation

The stratified inclined duct experiment (hereafter abbreviated ‘SID’) is sketched in figure 1(a). This conceptually simple experiment consists of two reservoirs initially filled with aqueous solutions of different densities $\rho_{0}\pm\Delta\rho/2$, connected by a long rectangular duct that can be tilted at an angle $\theta$ from the horizontal. At the start of the experiment, the duct is opened, initiating a brief transient gravity current. Shortly after, at $t=0$, an exchange flow starts and is sustained through the duct for long periods of time, until the accumulation of fluid of a different density from the other reservoir reaches the ends of the duct and the experiment is stopped at $t=T$ (typically after several minutes and many duct transit times). This exchange flow has at least four qualitatively different flow regimes, based on the experimental parameters, as we discuss in more detail in § 3.1.

Our notation is shown in figure 1(b) and largely follows that of Lefauve et al. (2018, 2019). The duct has length $L$, height $H$, and width $W$. The streamwise $x$ axis is aligned along the duct, the spanwise $y$ axis is aligned across the duct, making the $z$ axis tilted at an angle $\theta$ from the vertical (resulting in a non-zero streamwise projection of gravity $g\sin\theta$). The angle $\theta$ is defined to be positive when the bottom end of the duct sits in the reservoir of lower density, as shown here. The velocity vector field is $\textbf{u}(x,y,z,t)=(u,v,w)$ along $x,y,z$, and the density field is $\rho(x,y,z,t)$. All spatial coordinates are centred in the middle of the duct, such that $(x,y,z,t)\in[-L/2,L/2]\times[-W/2,W/2]\times[-H/2,H/2]\times[0,T]$.

Next, we define two integral scalar quantities of particular interest in exchange flows:

(i)   $Q$ 

the volume flux as the volumetric flow rate averaged over the duration of an experiment

[TABLE]

where $\langle|u|\rangle_{x,y,z,t}\equiv 1/(LWHT)\int_{0}^{T}\int_{-H/2}^{H/2}\int_{-W/2}^{W/2}\int_{-L/2}^{L/2}|u|\,\mbox{d}x\,\mbox{d}y\,\mbox{d}z\,\mbox{d}t$. The volume flux $Q>0$ measures the magnitude of the exchange flow between the two reservoirs. It is different from the net (or ‘barotropic’) volume flux $\langle u\rangle_{x,y,z,t}\approx 0$, since, to a good approximation, the volume of fluid in each reservoirs is conserved during an experiment (assuming the levels of the free surface in each reservoir are carefully set before the start of the experiment).

(ii)   $Q_{m}$

the mass flux as the net flow rate of mass averaged of the duration of an experiment

[TABLE]

which is equivalent to a buoyancy flux up to a multiplicative constant $g$. By definition $0<Q_{m}\leq Q$. The first inequality holds since, in our notation, negatively-buoyant fluid ($\rho_{0}<\rho\leq\rho_{0}+\Delta\rho/2$) flows on average to the right ($u>0$) and conversely. The second inequality would be an equality in the absence of molecular diffusion inside the duct (i.e. if all fluid moving right had density $\rho_{0}+\Delta\rho/2$ and conversely). In any real flow, laminar (and potentially turbulent) diffusion at the interface are responsible for an interfacial layer of intermediate density $|\rho-\rho_{0}|<\Delta\rho/2$ of finite thickness $\delta>0$ (figure 1(b)).

2.2 Non-dimensionalisation

A total of seven parameters are believed to play important roles in the SID: four geometrical parameters: $L$, $H$, $W$, $\theta$, and three dynamical parameters: the reduced gravity $g^{\prime}\equiv g\Delta\rho/\rho_{0}$ (under the Boussinesq approximation $0<\Delta\rho/\rho_{0}\ll 1$), the kinematic viscosity of water ($\nu=1.05\times 10^{-6}$ m2 s*-1*) and the molecular diffusivity of the stratifying agent (active scalar) $\kappa$. In this paper, we will primarily consider salt stratification ($\kappa_{S}=1.50\times 10^{-9}$ m2 s*-1*), but will also discuss temperature stratification ($\kappa_{T}=1.50\times 10^{-7}$ m2 s*-1*). From these seven parameters having two dimensions (of length and time), we construct five independent non-dimensional parameters below.

The first three non-dimensional parameters are geometrical: $\theta$, and the aspect ratios of the duct in the longitudinal and spanwise direction, respectively:

[TABLE]

We choose to non-dimensionalise lengths by the length scale $H/2$, defining the non-dimensional position vector as $\tilde{\textbf{x}}\equiv\textbf{x}/(H/2)$ such that $(\tilde{x},\tilde{y},\tilde{z})\in[-A,A]\times[-B,B]\times[-1,1]$. As an exception, we choose to non-dimensionalise the typical thickness of the interfacial density layer by $H$, for consistency with other definitions in the literature: $\tilde{\delta}\equiv\delta/H$, such that $\tilde{\delta}\in[0,1]$.

The last two non-dimensional parameters are dynamical. We define an ‘input’ Reynolds number based on the velocity scale $\sqrt{g^{\prime}H}$ and length scale $H/2$:

[TABLE]

Consequently, we non-dimensionalise the velocity vector as $\tilde{\textbf{u}}\equiv\textbf{u}/\sqrt{g^{\prime}H}$, and time by the advective time unit $\tilde{t}\equiv 2\sqrt{g^{\prime}/H}t$ (hereafter abbreviated ATU). We define our last parameter, the Prandtl number, as the ratio of the momentum to active scalar diffusivity:

[TABLE]

where $\kappa$ takes the value $\kappa_{S}$ or $\kappa_{T}$ depending on the type of stratification (salt or temperature, giving respectively $Pr=700$ and $Pr=7$). Finally, we define the non-dimensional Boussinesq density field as $\tilde{\rho}\equiv(\rho-\rho_{0})/(\Delta\rho/2)$, such that $\tilde{\rho}\in[-1,1]$.

We now reformulate the aim of this paper (introduced in § 1) more specifically as: exploring the behaviour of flow regimes, mass flux $\tilde{Q}_{m}$, and interfacial layer thickness $\tilde{\delta}$ in the five-dimensional space of non-dimensional input parameters $(A,B,\theta,Re,Pr)$.

In the next section we address the dimensional scaling of the velocity in the experiment. By discussing the a priori influence of the input parameters identified above on the velocity scale in this problem, we will provide a basis for subsequent scaling arguments in the paper.

2.3 Scaling of the velocity

Having constructed our Reynolds number (4) using the velocity scale $\sqrt{g^{\prime}H}$, here we show that it is the relevant velocity scale to use in such exchange flows. As sketched in figure 1(b), we define the typical peak-to-peak velocity as $\Delta U$. This velocity scale is not set by the experimenter as an input parameter, rather it is chosen by the flow as an output parameter. From dimensional analysis, we write

[TABLE]

In order to show that our Reynolds number (4) and our non-dimensionalisation the velocity by $\sqrt{g^{\prime}H}$ are relevant (and such that $\tilde{u}\in[-1,1]$), we will show below that we indeed expect $\Delta U/2\sim\sqrt{g^{\prime}H}$ and $f_{\Delta U}(A,B,\theta,Re,Pr)\sim 1$. Although some aspects of this discussion can be found in Lefauve et al. (2018, 2019), the importance of this dimensional analysis for this paper justifies the more detailed discussion that we offer below.

The velocity scale $\Delta U$ in quasi-steady state results from a dynamical balance in the steady, horizontal momentum equation under the Boussinesq approximation (in dimensional units)

[TABLE]

In addition to the standard inertial (I) and viscous (V) terms, this equation highlights the two distinct ‘forcing’ mechanisms in SID flows:

(H)   

a hydrostatic longitudinal pressure gradient, the minimal ingredient for exchange flow, resulting from each end of the duct sitting in reservoirs at different densities. This hydrostatic pressure in the duct increases linearly with depth $\p_{x}p=g\cos\theta\Delta\rho/(4L)z$, driving a flow in opposite directions on either side of the neutral level $z=0$: $-(1/\rho_{0})\p_{x}p=g^{\prime}\cos\theta/(4L)z$;

(G)   

a gravitational body force reinforcing the flow by the the acceleration of the positively-buoyant layer upward (to the left in figure 1) and of the negatively-buoyant layer downward (to the right) when the tilt angle is positive $g\sin\theta>0$ (the focus of this paper), and conversely when the tilt angle is negative.

Rewriting (7) in non-dimensional form and ignoring multiplicative constants, we obtain

[TABLE]

where $\ell$ is the smallest length scale of density gradients ($\ell=\delta$ in laminar flows, and $\ell\ll\delta$ in turbulent flows).

To simplify this complex ‘four-way’ balance, it is instructive to consider the four possible ‘two-way’ dominant balances to deduce four possible scalings for $\Delta U$ (ignoring constants and assuming $\cos\theta\approx 1$ since the focus of this paper is on small angles).

(IH)  

The inertial-hydrostatic balance. First, we can neglect the gravitational (G) term with respect to the hydrostatic (H) term if $g^{\prime}H\cos\theta\gg g^{\prime}L\sin\theta$, i.e. when the tilt angle of the duct $\theta$ is much smaller than its ‘geometrical’ angle:

[TABLE]

where we define the geometrical angle as

[TABLE]

Second, we can neglect the viscous (V) term if $g^{\prime}H\gg\nu\Delta U\ell^{-2}L$, i.e. if the Reynolds number is larger than $Re\gg HL/\ell^{2}$. This corresponds to

[TABLE]

in laminar flow (ignoring the case $B\ll 1$ for simplicity), and to a larger lower bound in turbulent flows. Under these conditions, balancing I and H gives the scaling $\Delta U\sim\sqrt{g^{\prime}H}$, i.e. $f_{\Delta U}\sim 1$, which corresponds to our choice in § 2.1.

(IG)  

The inertial-gravitational balance. Using analogous arguments, if $\theta\gg\alpha$ and $Re\gg HL/\ell^{2}$, we expect the scaling $\Delta U\sim\sqrt{g^{\prime}L\sin\theta}$, i.e. $f_{\Delta U}(A,\theta)\sim\sqrt{A\sin\theta}\gg 1$.

(HV)  

The hydrostatic-viscous balance. If $\theta\ll\alpha$ and $Re\ll A$, we expect $f_{\Delta U}(A,B,Re)\sim A^{-1}Re\ll 1$ (some dependence on $B$ being unavoidable in such a viscous flow).

(GV)  

The gravitational-viscous balance. If $\theta\gg\alpha$ and $Re\ll A$, we expect $f_{\Delta U}(B,\theta,Re)\sim\sin\theta Re\ll A$.

Figure 2 summarises the above analysis and the following conclusions.

(i)    

The parameters $A$, $\theta$ and $Re$ play particularly important roles in SID flows, since the variation of $\theta$ and $Re$ above or below thresholds set by $A$ can alter the scaling of $\Delta U$ (i.e. $f_{\Delta U}$). The parameter $B$ appears less important in this respect (except in narrow ducts where $B\ll 1$ and the $Re$ threshold becomes $AB^{-2}$).

(ii)   

At low tilt angles $0<\theta\ll\alpha$, $f_{\Delta U}$ increases from $\ll 1$ when $Re\ll A$ to $\sim 1$ when $Re\gg A$. At high enough $Re$, $f_{\Delta U}$ likely retains a dependence on $A,B,Re$ due to turbulence (the constant ‘IH’ scaling being a singular limit for $Re\rightarrow\infty$).

(iii)  

At high tilt angles $\theta\gg\alpha$ and Reynolds number $Re\gg A$, $f_{\Delta U}$ should increase well above $1$, and likely retains a dependence on $A,B,\theta,Re$ (the ‘IG’ scaling being a singular limit for $Re\rightarrow\infty$).

(iv)  

The blue rectangle in figure 2 represents the region of interest in most exchange flows of practical interest and in this paper. In this region, three or four physical mechanisms must be considered simultaneously (IHV, IGHV or IGV). Since few flows satisfy $\theta\ll\alpha$ or $\gg\alpha$, we consider that in general $f_{\Delta U}=f_{\Delta U}(A,B,\theta,Re,Pr)$ (the $Pr$ dependence reflects the fact that the active scalar can no longer be neglected at high $Re$ due to its effect on turbulence and mixing).

(v)   

Although the above ‘two-way’ balances do not allow us to confidently guess the scaling of $f_{\Delta U}$ in the blue region, theoretical arguments and empirical evidence support $f_{\Delta U}(A,B,\theta,Re,Pr)\sim 1$ for IHV, IGHV and IGV flows. This in agreement with the IH scaling and with (ii), but in disagreement with (iii).

The last conclusion (v) is supported by the theoretical concept of hydraulic control of two-layer exchange flows, which dates back to Stommel & Farmer (1953); Wood (1968, 1970) and was formalised mathematically by Armi (1986); Lawrence (1990); Dalziel (1991). In steady, inviscid, irrotational, hydrostatic (i.e. ‘IH’) exchange flows, the ‘composite Froude number’ $G$ is unity, which using our notation and assuming streamwise invariance of the flow ($\p_{x}=0$), reads:

[TABLE]

Such exchange flows are called maximal: the phase speed of long interfacial gravity waves $\sqrt{g^{\prime}H}$ ‘controls’ the flow at sharp changes in geometry (on either ends of the duct), and sets the maximal non-dimensional volume flux to $\tilde{Q}=\nicefrac{{1}}{{2}}$.

In ‘plug-like’ hydraulic flows ($Re\rightarrow\infty$), the velocity in each layer $\Delta U/2$ is equal to its layer-average $Q$, giving an upper bound $f_{\Delta U}=\tilde{Q}=\nicefrac{{1}}{{2}}$. By contrast, in real-life finite-$Re$ flows, the peak $\Delta U/2$ is larger than the average $Q$ (typically by a factor $\approx 2$), such that the upper bound is $f_{\Delta U}\approx 2\tilde{Q}\approx 1$. This upper bound remains approximately valid throughout the blue region of figure 2. We thus answer the question motivating this section: our choice of non-dimensionalising u by $\sqrt{g^{\prime}H}\approx\Delta U/2$ in order to have $|\tilde{\textbf{u}}|\lesssim 1$ is indeed relevant to SID flows.

Henceforth, we drop the tildes and, unless explicitly stated otherwise, use non-dimensional variables throughout.

3 Literature review

In this section we review the experimental literature on the questions of flow regimes (§ 3.1), mass flux (§ 3.2), and interfacial layer thickness (§ 3.3). For conciseness, we limit the discussion to the most relevant results of each study, and give further details about the geometry and parameters $A,B,\theta,Re,Pr$ used by each in appendix A (table A provides a synthesis of this literature review). We then highlight the main limitations of previous research to further motivate the paper in § 3.4.

3.1 Flow regimes

Macagno & Rouse (1961) (MR61) constitutes, to our knowledge, the first experimental study in a setup similar to the SID. MR61 used dye visualisations to describe four qualitatively different regimes:

L    

‘uniform laminar motion with straight streamlines’;

W  

‘laminar motion with regular waves’;

I    

‘incipient turbulence, with waves which break and start to show irregularity and randomness’;

T    

‘pronounced turbulence and active mixing across the interface’

MR61 mapped the above regimes together with measurements of the interfacial stress and mixing coefficients in the plane $(F_{*},R_{*})$, where $F_{*}\approx 2\sqrt{2}Q$ and $R_{*}\approx 4QRe$ (using our notation) are ‘effective’ Froude and Reynolds numbers. Arguing that flows above a certain $R_{*}$ and $F_{*}$ would be unstable, they proposed and experimentally verified that the ‘transition curves’ separating the flow regimes and the iso-curves of interfacial shear and mixing coefficients scaled with $R_{*}\,F_{*}\approx Re\,Q^{2}$ (i.e. these curves are $Re\,Q^{2}=$ const.). As we have seen in § 2.3, $Q$ is in reality a dependent variable, not an input parameter. This confusion in MR61 comes from the fact that their setup (which they attribute to Helmholtz) differs from the SID in that they were (to some extent) able to prescribe the volume flux $Q$ by controlling the inflow of salt water by a piston communicating with one of the reservoirs (their system was closed, i.e. it had no free surface). They varied $Q$ together with $\theta$ in non-trivial ways in order to reach target values $R_{*},F_{*}$. They did not appear to realise that the flow was hydraulically controlled and that $\theta$ and $Re$ were the relevant independent input parameters.

Wilkinson (1986) (W86) used shadowgraph and observed regime transitions similar to MR61 in a horizontal, circular pipe: ‘shear-induced instabilities […] initially in the form of cusp-like waves, but as the shear was further increased, Kelvin-Helmholtz billows were seen to grow and collapse creating a turbulent shear layer’. He suggested a scaling in $Re$ alone, independent of $A$: laminar flow under $Re<2450$, ‘interfacial waves radiating in both directions’ for $Re\in[2600,2700]$, and turbulence for $Re>2700$, but his experiments were limited in number ($\approx 18$).

Kiel (1991) (K91) used shadowgraph and laser sheet visualisations at larger $Re$, and classified the regimes differently: laminar; turbulent with $\delta<1$; and turbulent with $\delta=1$. Using a semi-empirical model that we will describe in § 3.2, he proposed regime transitions scaling with $A\tan\theta$, independently of $Re$, i.e. the opposite of W86.

Meyer & Linden (2014) (ML14) used shadowgraph visualisations, and (unaware of MR61) described the following four regimes (see their figure 3):

L     

laminar flow with a thin, flat density interface;

H     

mostly laminar flow with quasi-periodic waves on the density interface, identified as Holmboe waves;

I      

spatio-temporally intermittent turbulence with small-scale structures and noticeable mixing between the two layers;

T     

statistically-steady turbulent flow with a thick interfacial density layer.

This is the nomenclature adopted by LPL19 and that we adopt in this paper. The only difference between the MR61 and ML14 nomenclatures lies in the letter characterising the wavy regime (W in MR61 and H in ML14), simply because MR61 observed Holmboe waves before they were explained by Holmboe (1962). ML14 mapped these regimes in the $(\theta,Re)$ plane for two different $A=15,\,30$. They argued that, because the flow was hydraulically controlled, the ‘excess kinetic energy’ gained by the flow at $\theta>0$ should be dissipated turbulently (we identify this energy as the square of the ‘IG’ velocity scaling $g^{\prime}L\sin\theta$ of § 2.3). By non-dimensionalising this excess ‘IG’ energy by $(\nu/H)^{2}$, ML14 proposed and verified experimentally that regime transitions scale with the Grashof number

[TABLE]

assuming $\sin\theta\approx\theta$. This scaling has two limitations: the ‘IG’ energy does not explain the transitions at $\theta=0$, and its non-dimensionalisation by $(\nu/H)^{2}$ lacks a physical basis.

Lefauve et al. (2019) (LPL19) repeated the shadowgraph observations of ML14 in a smaller duct ($H=45$ mm vs $H=100$ mm) with otherwise equal parameters $(A,B,Pr)=(30,1,700)$ and mapped the regimes in the $(\theta,Re)$ plane. Surprisingly, their observations do not support the $Gr=4\times 10^{7}$ curve proposed by ML14 for the I to T transition. This suggests that $H$ plays a major role other than through the non-dimensional parameters $A,B,Re$, and therefore the existence of at least an additional non-dimensional parameter that we have hitherto overlooked. LPL19 observed two distinct scalings: a $\theta Re^{2}$ scaling for $\theta\lesssim\alpha$, in agreement with ML14), and a $\theta Re$ scaling for $\theta\gtrsim\alpha$, which they explained using a physical model. They developed from first principles detailed energy budgets which they applied to 16 experiments in which the full density field and three-component velocity field were simultaneously measured in a three-dimensional volume of the duct at high spatio-temporal resolution (for visualisations of flow fields in all four regimes, see their figures 2-3). They showed, theoretically and experimentally, that for $\theta\gtrsim\alpha$ (for so called ‘forced flows’), the time- and volume-average rate of dissipation of kinetic energy could be predicted a priori as

[TABLE]

where $\mathsf{s}$ is the non-dimensional strain rate tensor and $\sin\theta\approx\theta$ was assumed since their focus was on small angles. They showed that because the magnitude of streamwise velocities and wall shear stresses were bounded by hydraulic control, the requirement of high strain rates at high $\theta Re$ caused transitions to increasingly three-dimensional (turbulent) flow regimes. The $\theta Re$ scaling of energy dissipation matched the observed regime transitions in ‘forced flow’ ($\theta\gtrsim\alpha$), but the $\theta Re^{2}$ transition scaling in ‘lazy flows’ ($\theta\lesssim\alpha$ and in particular $\theta=0$) remains unexplained.

In summary, there is good consensus in the literature on the nomenclature of flow regimes observed in the SID and, in particular, on the fact that the flow becomes increasingly disorganised and turbulent with $A$, $\theta$ and $Re$. At a fixed $\theta\geq 0^{\circ}$, all flow regimes ($\textsf{L},\textsf{H},\textsf{I},\textsf{T}$) can be visited by increasing $Re$, and conversely at fixed $Re$ and increasing $\theta$ (W86, K91, ML14, LPL19). Both K91 and ML14 observed regime transitions scaling with $A\tan\theta=\tan\theta/\tan\alpha$ (or $A\theta$ for small angles), i.e. $A$ controls the $\theta$ scaling. However, the scaling in $Re$ is subject to debate, and may change on either side of $\theta\approx\alpha$ (LPL19). These conclusions are illustrated schematically in figure 3(a) (the interrogation marks denote open questions)

3.2 Mass flux

Leach & Thompson (1975) (LT75) measured $Q_{m}=0.23$ in horizontal circular pipes for high Reynolds number $Re=O(10^{4}-10^{5})$, and $Pr=1$ and $700$ (respectively CO2/air and salt/fresh water). Surprisingly, they observed no dependence on $A,Re,Pr$.

Mercer & Thompson (1975) (MT75) reported dramatic non-monotonicity of $Q_{m}(A,\theta)$ (we reproduce some of their data in appendix B, figure B.1(b)). They measured $Q_{m}\approx 0.2-0.3$ at $\theta=0^{\circ}$ (in agreement with LT75), increasing to $Q_{m}\approx 0.4$ at $\theta\approx\alpha/2$, and decreasing to $Q_{m}\approx 0.01-0.1$ at $\theta=90^{\circ}$. In a small additional set of experiments at $\theta=30^{\circ}$ in a larger pipe ($Re=2\times 10^{4}$ vs $2\times 10^{3}$, and $A=6$), they observed dependence on $Re$ even in the ‘very high’ range $Re\in[300A,3000A]$ (though it might be due to subtle differences in apparatus).

W86 developed a Bernoulli model in a horizontal circular pipe which predicted an upper bound of $Q=\pi/8\approx 0.39$ (non-dimensionalised using the pipe diameter), making the analogy with the hydraulic control arguments in Wood (1970) who predicted $Q=0.5$ in rectangular ducts. Including viscous boundary layers at the circular walls, he predicted and verified experimentally a monotonic increase of $Q$ with $A^{-1}Re$ (as the thickness of boundary layers decreases): $Q_{m}=0.13$ at $Re\approx 20A$ to $Q_{m}=0.35$ at $Re\approx 500A$ (larger than LT75), in agreement with the dimensional analysis of § 2.3 (conclusion (ii)).

K91 developed an inviscid Bernoulli model in an inclined duct for two counter-flowing layers of equal thickness and predicted $Q\approx\sqrt{(4/9)\cos\theta+A\sin\theta}$. In agreement with our dimensional analysis in § 2.3, this expression predicts a transition from an ‘inertial-hydrostatic’ (IH) balance at $0<\theta\ll\alpha$ with $Q\approx 2/3$ to an ‘inertial-gravitational’ (IG) balance at $\theta\gg\alpha$ with $Q\approx\sqrt{A\sin\theta}$. K91 showed, however, that this ‘IG’ scaling could only be observed experimentally when communication and mixing between the two counter-flowing layers was artificially suppressed by a rigid ‘splitter plate’ along the duct. He argued that the non-realisation of the IG scaling was due to a turbulent transition occurring when the IG scaling for $Q$ ‘that potentially exists’ exceeds a threshold dependent on the ‘stabilising effect of $g^{\prime}\cos\theta$’ (i.e. the maximal $Q=\nicefrac{{1}}{{2}}\sqrt{\cos\theta}$). He summarised this argument in a ‘geometric Richardson number’ $Ri_{G}$, whose inverse we interpret as being the square ratio of the ‘potential’ to the maximal $Q$

[TABLE]

K91 obtained a good empirical collapse of MT75’s $Q_{m}$ data (figure B.1(c)), as well as of his own $Q_{m}$ data (figure B.1(e-f)) with $Ri_{G}$. The mass flux peaks at $Ri_{G}\approx 0.2$, i.e. $\theta\approx\alpha/2$, and decays at larger tilt angles as the flow regime is increasingly turbulent. In agreement with W86’s arguments, K91 reported independence of his results with $Re$ above $Re>400A$ (figure B.1(a)), and intentionally focused on these high $Re$ throughout.

ML14 observed monotonic increase of $Q_{m}(\theta)$ with $Q_{m}\approx 0.2-0.3$ at $\theta=0^{\circ}$ and $Q_{m}\approx 0.5$ at $\theta=\alpha/2$. They did not comment on the hint of non-monotonic behaviour ($Q_{m}\lesssim 0.5$) suggested by their data at $\theta\approx 2\alpha$.

LPL19 observed (in passing) clear non-monotonic behaviour of $Q_{m}(\theta,Re)$. Their data are well fitted by a hyperbolic paraboloid in the $\log\theta-\log Re$ plane, where $Q_{m}=$ const. curves are hyperbolas, with $Q_{m}\approx 0.5$ along the major axis $\theta Re^{3/2}=$ 100 ($\theta$ in radians), and $Q_{m}$ decays on either side of it.

In summary, we can infer that the mass flux has a complicated non-monotonic behaviour in $A,\theta,Re$ as sketched in figure 3(b). While the dependence on $Re$ is clear at $Re<500A$ (MR61, W86, ML14 and LPL19) due to the influence of viscous boundary layers, it is still debated at $Re>500A$ (MT75 and ML14 argue in favour, K91 argues against). The mass flux $Q_{m}$ reaches a maximum $Q_{m}\approx 0.4-0.5$ at $\theta\approx\alpha/2$ and ‘high enough’ $Re$ (MT75, K91, ML14, LPL19), and decays on either side for smaller/larger $\theta$ and $Re$ (W86, LPL19) in a poorly-studied fashion.

3.3 Interfacial layer thickness

K91 is, to our knowledge, the only thorough experimental study of $\delta(A,\theta)$. K91 performed conductivity probe measurements and observed monotonic increase with both $A$ and $\theta$ (figure B.1(f)) and independence on $Re$. He found good collapse with $Ri_{G}$ (figure B.1(g)). The interfacial mixing layer is turbulent and thick ($\delta\approx 0.4-0.7$) when the maximal $Q_{m}$ is reached at $Ri_{G}\approx 0.2$ ($\theta\approx\alpha/2$), and it fills the whole duct height ($\delta=1$) at $Ri_{G}<0.1$ ($\theta\gtrsim 2\alpha$), see figure 3(c). At even larger tilt angles, the mean vertical density gradient $|\rho(z=1)-\rho(z=-1)|/2$ drops below 1 (this ‘extreme’ turbulent scenario falls outside the scope of this paper). The behaviour of $\delta$ at low $Re<500A$ has not yet been studied.

3.4 Limitations of previous research

The above experimental literature (summarised in figure 3) provides much quantitative insight to complement the qualitative predictions of the scaling analysis in § 2.3 (summarised in figure 2).

However, many aspects of the scaling of regimes, $Q_{m}$ and $\delta$ with $A,B,\theta,Re,Pr$ remain open questions (symbolised by the interrogation marks in figure 3). As extreme examples, the effects of $Re$ on $\delta$, and the effects of $B$ and $Pr$ on all three variables have not been studied at all.

Moreover, despite our above efforts to unify their findings, these past studies of the SID experiment inherently provide a somewhat fragmented view of the problem due to the following limitations (made clear by table A):

(i)    

they used slightly different setups and geometries (e.g. presence vs absence of free surfaces in the reservoirs, rectangular ducts vs circular pipes), and slightly different measuring methodologies (e.g. for $Q_{m}$);

(ii)   

only one study (K91) addressed the interdependence of the three variables of interest (regime, $Q_{m}$, $\delta$), while the remaining studies measured either only regimes (MR61), only $Q_{m}$ (LT75, MT75), or both (ML14, LPL19);

(iii)  

they focused on the variation of a single parameter (MR61), two parameters (W86, K91, LPL19), or at most three parameters (MT75, ML14) in which case the third parameter took only two different values;

(iv)  

they studied limited regions of the parameter space, and it is difficult to confidently interpolate results obtained by different setups in different regions (such as $Re<500A$ and $>500A$).

The experimental results in the next section attempt to overcome the above limitations by providing a more unified view of the problem.

4 Experimental results

In order to make progress on the scaling of flow regimes, $Q_{m}$ and $\delta$ with $A,B,\theta,Re,Pr$, we obtained a comprehensive set of experimental data using an identical setup, measuring all three dependent variables with the same methodology (described in appendix C), and varying all five independent parameters in a systematic fashion. We introduce the different duct geometries and data sets used in § 4.1, and present our results on flow regimes in § 4.2, on mass flux in § 4.3, and on interfacial layer thickness in § 4.4.

4.1 Data sets

All experimental data presented in the following were obtained in the stratified inclined duct (SID) setup sketched in figure 1. We used four different duct geometries and two types of stratification (salt and temperature) to obtain the following five distinct data sets, listed in table 1:

LSID     

(L for Large) with height $H=100$ mm, and $A=30$, $B=1$;

HSID     

(H for Half) which only differs from the LSID (the ‘control’ geometry) in that it is half the length: $A=15$ (highlighted in bold in table 1);

mSID     

(m for mini) which only differs from the LSID in its height $H=$45\text{,}\mathrm{mm}$$ (in bold), but should yield identical data to LSID for identical $\theta,Re,Pr$ since $H$ should only play a role through the non-dimensional parameters $A$, $B$ (identical in LSID and mSID) and $Re$. Our emphasis on ‘should’ denotes the fact that we have already discussed in the literature review § 3.1 that this was already observed not to be the case: the I to T regime transition curve of ML14 in LSID did not overlap with that of LPL19 in mSID;

tSID      

(t for tall) which differs from the HSID primarily in its tall spanwise aspect ratio $B=\nicefrac{{1}}{{4}}$ (and, secondarily, in a marginally smaller height $H=90$ mm);

mSID ​ ​  

(m for mini and ​ ​for temperature) which differs from the mSID in that the stratification was achieved by different reservoir water temperatures (hence $Pr=7$), as opposed to different salinities in the above data sets (where $Pr=700$). This limited the density difference $\Delta\rho$ achieved, reflected in the lower $Re$.

Table 1 also lists, for each data set, the range of variation of $\theta$ and $Re$, and the number of data points, i.e. distinct $(\theta,Re)$ couples for which we have data on regime, $Q_{m}$, and $\delta$.

Note that the regime and $Q_{m}$ data of the top three data sets have already been published in some form by Meyer & Linden (2014) (ML14, denoted by ∗) and Lefauve et al. (2019) (LPL19, denoted by †), as discussed in the literature review § 3.1-3.2. However, ML14 plotted their LSID and HSID data together (regimes in the $(\theta,Gr)$ plane, see their figure 7 and $Q_{m}$ as a function of $\theta$ alone, see their figure 8) and did not investigate their potential differences, while LPL19 only commented in passing on a fit of the $Q_{m}$ data in the $(\theta,Re)$ plane (see their figure 9). The individual reproduction and thorough discussion of these data alongside more recent data using a unified non-dimensional approach will be key to this paper. All five data sets have been used in the PhD thesis of Lefauve (2018) (especially in Chapters 3 and 5), and the detailed parameters of all experiments are tabulated in his appendix A for completeness.

Our focus on long ducts, evidenced by our choice of $A=15$ and $30$, reflects our primary motivation, i.e. flows relevant to geophysical and environmental applications, which are typically largely horizontal ($\theta\approx 0^{\circ}$) and stably stratified in the vertical (as opposed to the different case of vertical exchange flow with $\theta=90^{\circ}$). As we have seen in § 3, the SID experiment conveniently exhibits all possible flow regimes, including high levels of turbulence and mixing, between $\theta=0^{\circ}$ and a few $\alpha$ at most. In long ducts (large $A$), $\alpha\equiv\tan^{-1}(A^{-1})$ is therefore small enough to allow us to study all the key dynamics of sustained stratified flows while keeping $\theta$ small enough for these flows to remain largely horizontal and relevant to our motivation.

As a result of this focus on long ducts, in the remainder of the paper we make the approximation that

[TABLE]

as in our previous discussion of the results of ML14 and LPL19 in § 3. This approximation is accurate to better than $2~{}\%$ for the angles considered in our data sets ($\theta\leq 10^{\circ}$). Unless explicitly specified by the ∘ symbol, $\theta$ will now be expressed in radians (typically in scaling laws).

4.2 Flow regimes

The $\textsf{L},\textsf{H},\textsf{I},\textsf{T}$ flow regimes were determined following the ML14 nomenclature as in appendix C.1 (except for a new regime which we discuss in the next paragraph). Figure 4 shows the resulting regime maps in the $(\theta,Re)$ plane corresponding to the five data sets.

First, we note the introduction of a ‘new’ W regime in the tSID and mSID data (panels (d,e)). This W (wave) regime is similar to the H (Holmboe) regime, but describes interfacial waves which could not be recognised as Holmboe waves in shadowgraphs. These waves were of two types. First, in the tSID geometry at positive angles $\theta>0$, the waves did not exhibit the distinctive ‘cusped’ shape of Holmboe waves and the waves appeared to be generated at the ends of the duct and to decay as they travel inside the duct. Second, in the mSID ​ ​larger-amplitude, tilde-shaped internal waves were observed across most of the height of the duct, contrary to Holmboe waves which are typically confined to a much thinner interfacial region. Further discussion of these waves falls outside the scope of this paper, but can be found in (Lefauve, 2018, §§ 3.2.3-3.2.4) (hereafter abbreviated L18). This new observation highlights the richness of the flow dynamics in the SID experiment. However, for the purpose of this paper, the H and W regimes are sufficiently similar in their characteristics (mostly laminar flow with interfacial waves) that we group them under the same regime for the purpose of discussing regime transitions.

The main observation of figure 4 is that the transitions between regimes can be described as simple curves in the $(\theta,Re)$ plane that do not overlap (or ‘collapse’) between the five data sets. The slope and location of the transitions varies greatly between panels: the difference between the LSID and HSID data (panels (a,b)) is due to $A$, the difference between the HSID and tSID data (panels (b,d)) is due to $B$, and the difference between the mSID and mSID ​ ​data (panels (c,e)) is due to $Pr$.

However one of the most surprising differences is that between LSID and mSID data (panel (a,c)), due to the dimensional height of the duct $H$, already observed in LPL19. It is reasonable to expect that this $H$-effect is responsible for the main differences between the LSID/HSID/tSID data and the mSID/mSID ​ ​data. In other words $H$ is the main reason why the LSID/HSID/tSID transitions curves lie well above those for mSID/mSID ​ ​, i.e. the same transitions occur at higher $Re$ for larger $H$. The factor of $\approx 2$ quantifying this observation suggests that a Reynolds number built using a length scale identical in all data sets (rather than $H/2$) would better collapse the data; however our dimensional analysis does not provide such a length scale.

We conclude that the transitions between flow regimes can be described by hyper-surfaces depending on all five parameters $A,B,\theta,Re,Pr$ because their projections onto the $(\theta,Re)$ plane for different $A,B,Pr$ do not overlap. This dependence of flow regimes on all five parameters was not immediately obvious from our dimensional analysis which concerned the scaling of the velocity $f_{\Delta U}$ alone (§ 2.3 and figure 2). Furthermore, the existence of at least another important non-dimensional parameter related to $H$ is a major result that could not be predicted a priori.

Let us now investigate in more detail the scaling of regime transitions with respect to $\theta$ and $Re$, for which we have much higher density of data than for $A,B,Pr$. We replot the $\theta>0$ data of figure 4 using a log-log scale in figure 5 (each panel corresponding to the respective panel of figure 4). To guide the eye to the two main types of regime transition scalings observed in these data, we also plot two families of lines: dashed lines with a $\theta Re=$ const. scaling, and dotted lines with a $\theta Re^{2}=$ const. scaling. We also show using grey shading special values of interest: $Re=50A$, and $\theta=\alpha$. The latter was highlighted as particularly relevant in the literature review § 3 (notably as the boundary between lazy and forced flows in LPL19). Although W86 and K91 quoted $Re=500A$ as a threshold beyond which the effects of viscosity should be negligible on the turbulence in the SID, we believe that $Re=50A$ is a physically justifiable threshold beyond which the influence of the top and bottom walls of the duct becomes negligible. In the absence of turbulent diffusion, laminar flow in the duct is significantly affected by the top and bottom walls if the interfacial and wall 99 % boundary layers overlap in the centre of the duct ($x=0$), which occurs for $Re<50$ (L18, § 5.2.3). If, on the other hand, $Re\gg 50A$ ($Re=500A$ being a potential threshold), the top and bottom wall laminar boundary layers (as well as the side wall laminar boundary layers, assuming that $B\centernot\ll 1$) do not penetrate deep into the ‘core’ of the flow (however at these $Re$, we expect interfacial turbulence to dominate the core of the flow). Note that black contours representing a fit of the $Q_{m}$ data are superimposed in panels (a-d); these will be discussed in the next section.

We observed regime transitions scaling with $\theta Re^{2}=$ const. (dotted lines) in LSID, tSID and mSID ​ ​(panels (a,d,e)), and with $\theta Re=$ const. (dashed lines) in HSID (panel (b)). In mSID (panels (c)), we observe these two different scalings: $\theta Re^{2}$ for $\theta\lesssim\alpha$ (lazy flows) and $\theta Re$ for $\theta\gtrsim\alpha$ (forced flows), as previously observed by LPL19. LPL19 physically substantiated the $\theta Re$ scaling in forced flows, but not the $\theta Re^{2}$ scaling in lazy flows. Furthermore, these five data sets show that this dichotomy in scalings between lazy and forced flows does not extend to all other geometries; indeed lazy flows in the HSID exhibit a $\theta Re$ scaling and forced flows in the mSID ​ ​exhibit a $\theta Re^{2}$ scaling. These observations further highlight the complexity of the scaling of regime transitions with $A,B,\theta,Re,Pr$.

4.3 Mass flux

Mass fluxes were determined using the same salt balance methodology as ML14 described in appendix C.2.

In figure 6 we plot the $Q_{m}$ data for mSID (full symbols) and tSID (open symbols) as a function of $Re$ for all the available $\theta$ (from $\theta=-1^{\circ}$ in panel (a) to $\theta=3.5^{\circ}$ in panel (j)). The colour of each symbol denotes the regime as in figures 4-5 and the error bars denote the uncertainty about the precise duration $T$ of the ‘steady’ flow of interest in an experiment (used to average the volume flux and obtain $Q_{m}$, see appendix C.2 for more details). We do not plot the LSID and HSID data in this figure because they are sparser and do not have error bars (these data were collected by ML14).

At low angles $\theta\lesssim 1^{\circ}<\alpha$ (where $\alpha\approx 2^{\circ}$ in mSID and $4^{\circ}$ in tSID) we observe low values $Q_{m}\approx 0.2-0.3$ in the L and H regimes. At higher angles $\theta\approx\alpha-2\alpha$ we observe convergence to the hydraulic limit $Q_{m}\rightarrow 0.5$ (denoted by the dashed line), as discussed in § 2.3, which coincides with the I and T regimes. We also note that this hydraulic limit is not a strict upper bound in the sense that we observe values up to $Q_{m}=0.6$ in some experiments (some error bars even going to $0.7$).

Moreover, the non-monotonic behaviour of $Q_{m}$ with respect to both $\theta$ and $Re$ is particularly clear in the mSID data: at ‘large’ angles $\theta\gtrsim\alpha\approx 2^{\circ}$, $Q_{m}$ drops with both $\theta$ and $Re$. All of these observations are in qualitative agreement with the literature (§ 3.2). However, similarly to the regime data of § 4.2, we also note that the mSID and tSID data do not collapse with $Re$: all tSID data (open symbols) are shifted to larger $Re$ compared to mSID data (full symbols). This again suggests that $H$ plays an important role not presently captured by our five non-dimensional parameters, and that a Reynolds number based on a length scale independent of $H$ would better collapse the data.

To gain more insight into the scaling of $Q_{m}$ and its relation to the flow regimes, we superimpose on the regime data of figure 5(a-d) black contours representing the least-squares fit of our four $Q_{m}$ data sets using the following quadratic form:

[TABLE]

This is the general equation of a conic section, well suited to the non-monotonic behaviour observed above, and $\bm{\Gamma}$ is commonly referred to as the matrix of the quadratic equation.

These contours describe hyperbolas ($\det\bm{\Gamma}<0$) for LSID, HSID and mSID (panels (a,b,c)), and concentric ellipses ($\det\bm{\Gamma}>0$) for tSID (panel (d)). The hydraulic limit $Q_{m}\approx 0.5$ is reached either at the saddle point of the hyperbolas (panels (a,b,c)), or at the centre of the ellipses (panels (d)), and, encouragingly, no $Q_{m}=0.6$ contour exists here.

We again note that these four data sets do not collapse in the $(\theta,Re)$ plane. For example, the angle at which this maximum $Q_{m}$ is achieved is a modest $\theta=0.3\alpha$ in mSID (panel (c)) but appears much larger in tSID. The eigenvectors of $\bm{\Gamma}$ for each data set reveal that the major axis of these conic sections has equation $\theta Re^{\gamma}$ where $\gamma=2.6,0.3,1.5,1.2$ respectively for panels (a,b,c,d) (a larger exponent $\gamma$ represents a larger dependence on $Re$, hence a more horizontal axis).

The exponent $\gamma$ characterising the slope of the major axis is roughly of the same order as the exponent characterising the lines of regime transition (which is 1 for the $\theta Re$ scaling, and 2 for the $\theta Re^{2}$ scaling), suggesting that both phenomena (regime transition and non-monotonic behaviour of $Q_{m}$) are linked. However, this agreement is not quantitative except in mSID (panel (c)) where $\gamma=1.5$ is precisely the average of the two different regime transition exponents. This general lack of correlation suggest that the relationship between regimes and $Q_{m}$ in the SID is not straightforward and dependent on the geometry.

4.4 Interfacial layer thickness

Interfacial layer thickness was determined using the non-intrusive shadowgraph imaging technique (in salt experiments only). Shadowgraph is particularly suited to detect peaks in the vertical curvature of the density field $\p_{zz}\rho$ which we define as the edges of the interfacial density layer, as explained in appendix C.3.

In figure 7 we plot $\delta$ for our four duct geometries (rows) and three particular angles (representing a subset of our data) $\theta=1^{\circ},2^{\circ},3^{\circ}$ (columns). In figure 8 we plot a quadratic fit (black contours) to all the available data (represented by the symbols) in the $(\log\theta,\log Re)$ plane following (17). We also added in grey shading the $\theta=\alpha$ and $Re=50A$ values of interest for comparison between panels. In both figures, the colour of the symbol denotes the flow regimes as in figures 4-6.

We observe in both figures that $\delta$ monotonically increases with both parameters $\theta$ and $Re$, starting from values as low as $\delta\approx 0.05$ in the L, H, and W regimes (see figure C.1(a) for an illustration with $\delta=0.069$), and ending with values as high as $\delta\approx 0.8$ in the T regime (see figure C.1(a) for an illustration with $\delta=0.47$). The upper bound corresponds to the turbulent mixing layer filling 80 % of the duct height, with unmixed fluid only filling the remaining top and bottom 10 %. We substantiate this value of the lower bound by the thickness of the 99 % laminar boundary layer resulting from the balance between streamwise advection and vertical diffusion of an initially step-like density field. This calculation gives, at any point in the duct, $\delta_{99}\approx 10A^{1/2}(Re\,Pr)^{-1/2}\approx 0.03-0.1$ in the range $Re\in[300,6000]$ where the $\textsf{L},\textsf{H},\textsf{W}$ regimes are found.

As already discussed in the regimes and $Q_{m}$ data, we also observe a clear role of the dimensional parameter $H$ in ‘shifting’ the LSID/HSID/tSID data to higher $Re$ than the mSID data and hindering their overlap. Note that $A$ and $B$ play additional, more subtle roles shown by the differences between the LSID and HSID data and between the HSID and tSID data respectively.

Finally, we observe in figure 8 that there is good agreement between iso-$\delta$ contours and ‘iso-regime’ curves, or regime transitions (not shown for clarity, but easily visualised by the colours of the symbols). This suggests that $\delta$ is more closely correlated to regimes than $Q_{m}$ is.

5 Models and discussion

In this section, we attempt to explain some of the above observations with three particular classes of models, whose prior success in the literature make them natural candidates to tackle this problem.

In § 5.1 we attempt to explain the scaling of regime transitions at high $Re\gg 50A$ by generalising the time- and volume-averaged energetics analysis of LPL19. In § 5.2, we investigate the scaling of regimes and $Q_{m}$ by developing a frictional two-layer hydraulic model. In § 5.3, we tackle the scaling of $\delta$ in the I and T regimes by a variety of turbulence mixing models.

5.1 Volume-averaged energetics

We recall from § 3 that the simultaneous, volumetric measurements of the density and three-component velocity fields of LPL19 confirmed their theoretical prediction that in forced flows ($\theta\gtrsim\alpha$) the time- and volume-averaged norm of the strain rate tensor (non-dimensional dissipation) followed the scaling $\langle\mathsf{s}^{2}\rangle_{x,y,z,t}\sim\theta Re$ (14). They further decomposed the dissipation into:

(i)    

a ‘two-dimensional’ component $\mathsf{s}^{2}_{2d}$ (based on the $x-$averaged velocity $\textbf{u}_{2d}\equiv\langle\textbf{u}\rangle_{x}$). LPL19 measured flows in the mSID geometry at $Re<2500$, i.e. $Re\not\gg 50A=1500$, in which case the viscous interfacial and top and bottom wall boundary layers are well or fully developed and $\mathsf{s}^{2}_{2d}\sim\langle(\partial_{z}u_{2d})^{2}\rangle_{x,y,z,t}=O(1)$. They indeed observed that $\langle\mathsf{s}^{2}_{2d}\rangle{x,y,z,t}$ plateaus at $\approx 4$ in the I and T regimes due to the hydraulic limit;

(ii)   

a complementary ‘three-dimensional’ part $\mathsf{s}^{2}_{3d}=\mathsf{s}^{2}-\mathsf{s}^{2}_{2d}$ which, as a consequence of the plateau of $\mathsf{s}^{2}_{2d}$, takes over in the I and T regime and explains the $\theta Re$ scaling of regime transitions for forced flows in mSID.

In flows at $Re\gg 50A$ (well above the horizontal grey shading in figures 5, 8) we expect the 99 % viscous boundary layers to be of typical thickness $\delta_{99}\sim 10A^{1/2}Re^{-1/2}\ll 1$, and therefore volume-averaged two-dimensional dissipation to be higher $\mathsf{s}^{2}_{2d}\sim\langle(\partial_{z}u_{2d})^{2}\rangle_{x,y,z,t}\sim 10^{-1}A^{-1/2}Re^{1/2}\gg 1$. Therefore, we extend the prior results of LPL19 that regime transitions correspond to threshold values of

[TABLE]

by conjecturing that they correspond to threshold values of

[TABLE]

which introduces $A$ and a different exponent to $Re$ into the scaling.

Unfortunately we have little regime data for forced flows at $Re\gg 50A$ (upper right quadrants of each panel in figure 5) except in LSID (panel (a)). Nevertheless, it does not appear that this conjectured scaling would be able to explain the observed $\theta Re^{2}$ scaling. Detailed flow measurements would be required in this geometry to confirm or disprove the above two assumptions that two-dimensional dissipation follows a different scaling, and that regime transitions are tightly linked to three-dimensional dissipation.

Furthermore, we recall that the under-determination of the energy budgets of lazy flows ($\theta<\alpha$, see LPL19 figure 8*(a)*) does not allow us to predict the rate of energy dissipation ($\mathsf{s}^{2}$) from the rate of energy input ($\sim\theta Re$) and therefore to substantiate the transition scalings in lazy flows (left two quadrants of each panel in figure 5).

5.2 Frictional two-layer hydraulics

We introduce the fundamentals of this model in § 5.2.1, before examining the physical insight it provides in § 5.2.2, and its implications for the scaling of regime transitions and mass flux in § 5.2.3.

5.2.1 Fundamentals

The two-layer hydraulic model for exchange flows (figure 9(a)) assumes two layers flowing with non-dimensional velocities $u_{1}(x)>0$ (lower layer) and $u_{2}(x)<0$ (upper layer), and separated by an interface of non-dimensional elevation $\eta(x)\in[-1,1]$ above the neutral level $z=0$.

In the idealised inviscid hydraulic model (i.e. in the absence of viscous friction) the conservation of volume and of Bernoulli potential, and the requirement of hydraulic control yield a horizontal and symmetric interface $\eta(x)=0$ for $x\in[-A,A]$ and a volume flux $Q=u_{1}=-u_{2}=\nicefrac{{1}}{{2}}$ as already mentioned in § 2.3 (see appendix D.1 for more details).

The frictional hydraulic model is of more relevance to SID flows at finite $Re$. This model parameterises the effects of viscous friction while retaining the hydraulic assumptions (hydrostatic, steady, two-layer flow with uniform velocities $u_{1,2}(x)$). Dating back to Schijf & Schönfled (1953); Assaf & Hecht (1974); Anati et al. (1977)), it was formalised by Zhu & Lawrence (2000); Gu (2001); Gu & Lawrence (2005), who considered the effects of friction at the interface and bottom wall only, with applications to wide, open, horizontal channels. Here we further develop this model to add the effects of gravitational forcing ($\theta>0$) and friction at the top and side walls. The full development of this model can be found in (L18, § 5.2) and we offer a summary in appendix D. Some of its conclusions were briefly discussed in LPL19 § 4.3.1 (e.g. the distinction between lazy/forced flows). Below we provide a concise self-contained presentation of the key results of this model regarding the particular problem of the scaling of regimes and $Q_{m}$.

As sketched in figure 9(b), we consider that each infinitesimal duct sub-volume $dx\times 2B\times 2$ centred around $x$ is subject to horizontal, resistive stresses at the bottom wall $\tau^{Z}_{1}(x)$, top wall $\tau^{Z}_{2}(x)$ (in blue), side walls $\tau^{Y}_{1,2}(x)$ (respectively in the bottom and top layers, in green), and interface $\tau^{I}$ (in red). The inclusion of these stresses in the evolution of Bernoulli potential along the duct (see § D.1) yields a nonlinear differential equation for the slope of the interface along the duct of the form

[TABLE]

(see (52) for the full expression). Here $f_{Z},f_{Y},f_{I}$ are constant friction factors parameterising respectively the top and bottom wall stress, the side wall stress, and the interfacial stress, (they can be determined a posteriori from any finite-$Re$ flow profile, see § D.2 and (D.2)). For any set of parameters $\theta,Re,f_{Z},f_{Y},f_{I}$, this dynamical equation can be combined with the hydraulic control condition and solved numerically using an iterative method to yield a unique solution for $Q$ and $\eta(x)$ (§ D.3). The volume flux $Q$ generally increases with the forcing $\theta Re$, and decreases with friction $f_{Z},f_{Y},f_{I},$ and $A$.

5.2.2 Physical insight

We now consider the mid-duct slope $\eta^{\prime}(x=0)$, whose simplified expression shows the balance between the forcing $\theta Re$ and the ‘composite friction parameter’ $F$:

[TABLE]

and $r_{Y}\equiv B^{-1}f_{Y}/f_{Z}$ and $r_{I}\equiv f_{I}/f_{Z}$ are respectively the side wall friction ratio and the interfacial friction ratio.

We further note that our model has three properties: (i) the interface must slope down everywhere ($\eta^{\prime}(x)<0$) since the lower layer accelerates convectively from left to right ($u_{1}u^{\prime}_{1}(x)>0$) and conversely ($u_{2}u^{\prime}_{2}(x)<0$); (ii) the interface must remain in the duct $|\eta(x=\pm A)|<1$; (iii) $\eta^{\prime}$ always reaches a maximum ($|\eta^{\prime}|$ reaches a minimum) at the inflection point $x=0$.

From these properties we deduce that the existence of a solution requires the mid-duct interfacial slope to satisfy

[TABLE]

and therefore, using (21), we obtain the following bounds:

[TABLE]

The first inequality in (22) comes from property (ii) and means that the mid-duct interfacial slope must not be too steep compared to the duct geometrical slope $A^{-1}\approx\alpha$. The second inequality comes from (i) and (iii) and means that the mid-duct interfacial slope must be negative for $\eta(x)$ to be negative everywhere.

When suitably rescaled by $2Q\in[0,1]$, the combined friction parameter $F$ must therefore follow a $\theta Re$ scaling, strictly bounded from below. The upper bound in (23) is loose ($b>0$) in lazy flows, and tight ($b\rightarrow 0$) in forced flows ($A\theta\approx\theta/\alpha\gg 1$ and $Q\rightarrow\nicefrac{{1}}{{2}}$).

5.2.3 Implications for regimes and $Q_{m}$

Combining the above physical insight with our experimental observations, we conjecture the following behaviour about regimes and $Q_{m}$, summarised in figure 10.

(i)    

At zero or ‘low’ $\theta Re$ (i.e. at $\theta\approx 0$, since $Re$ must be large for hydraulic theory to hold) due to the inevitable presence of wall and interfacial friction ($F>0$) and the looseness of the upper bound $b$, $2QF$ is typically well above the forcing $\theta Re$. The friction $F$ is independent of $\theta Re$ and the flow is typically laminar (L regime). The interface has a noticeable slope all along the duct $\eta^{\prime}(0)\ll 0$, associated with a small volume flux $Q\ll\nicefrac{{1}}{{2}}$ (see (53)). Such lazy flows are underspecified, and the scaling of $Q$ and $F$ with $\theta Re$ is therefore impossible to predict a priori.

(ii)   

At moderate $\theta Re$ ($\theta>0$): $2QF$ increases above its ‘default’ $\theta=0$ value. This is achieved, on one hand, through an increase in $Q$ (and therefore $Q_{m}$), making the flow approach the hydraulic limit (panel (b)), and on the other hand, through an increase in $F$, in particular through laminar interfacial shear ($r_{I}$), rendering the flow unstable to Holmboe waves above a certain threshold ($\textsf{L}\rightarrow\textsf{H}$ transition, panel (c)). The phenomenology of this transition agrees with that proposed by the energetics of LPL19 (see their § 6.2-6.3). The fact that the $\textsf{L}\rightarrow\textsf{H}$ (or $\textsf{L}\rightarrow\textsf{W}$) transition exhibits different scalings in our different data sets is not presently understood. It may come from the complex, individual roles of $Q$ and $F$ in the precise flow profiles $u(y,z),\rho(z)$ responsible for triggering the Holmboe instability, and the different scalings of $Q$ and $F$ that could allow $2QF$ to follow a $\theta Re$ scaling.

(iii)  

At high $\theta Re$: the hydraulic limit is reached, the upper bound is tight ($b\approx 0$), the interface is mostly flat ($\eta(x)\approx 0$ everywhere), and the inequality (23) becomes $2QF\approx F\approx\theta Re$. In such forced flows, the friction parameter $F$ alone must precisely balance the forcing. Arbitrarily large $\theta Re$ requires arbitrarily large $F$, which we conjecture is largely achieved by turbulent interfacial friction (increase in $r_{I}$ responsible for the $\textsf{H}\rightarrow\textsf{I}$ and eventually the $\textsf{I}\rightarrow\textsf{T}$ transition).

From implication (iii), it is natural to conjecture that these two transitions are also caused by threshold values of the interfacial friction ratio $r_{I}$, which, as explain in appendix D.2, can be written $r_{I}=1+K_{I}$, where $K_{I}$ is a turbulent momentum diffusivity (non-dimensionalised by the molecular value $\nu$) parameterising interfacial Reynolds stresses (see (49)). Assuming that all wall shear stresses are similar ($r_{Y}\approx 1$), and that interfacial Reynolds stresses eventually dominate over laminar shear ($K_{I}\gg 1$), we have $K_{I}\approx F/(8f_{Z})$. For $Re<50A$, fully-developed boundary layers yield $f_{Z}\sim 1$, implying regime transitions scaling with (ignoring pre-factors)

[TABLE]

For $Re\gg 50A$, thin top and bottom wall boundary layer arguments similar to those of § 5.1 yield $f_{Z}\sim A^{1/2}Re^{-1/2}$, implying regime transitions scaling with

[TABLE]

Comparing (24)-(25) to (18)-(19) we see that the $Re<50A$ scaling obtained with frictional hydraulics is identical to that obtained by the energetics. However, the $Re\gg 50A$ scaling is different, and unfortunately it does not allow us to explain the regime transitions data (a $\theta Re^{1/2}$ or $\theta^{2}Re$ scaling is never observed). In addition, direct estimations of friction coefficients using three-dimensional, three-component velocity measurements in all flow regimes (L18, § 5.5) suggest a posteriori that the assumption that $K_{I}\gg 1$ might only hold beyond the $\textsf{I}\rightarrow\textsf{T}$ transition, undermining its usefulness to predict the $\textsf{H}\rightarrow\textsf{I}$ and $\textsf{I}\rightarrow\textsf{T}$ transitions.

From implication (ii), we understand why the volume flux $Q$, and hence the mass flux $Q_{m}$, both increase with $\theta$ and $Re$ in the L and H regimes, as observed in § 4.3. However lazy flows are under-specified; only one equation governs both the volume flux and friction ($2QF\sim\theta Re$), which does not allow us to obtain the value of the exponent $\gamma$ in the scaling $Q\sim\theta Re^{\gamma}$. From implication (iii), we conjecture two potential reasons for the decrease of the mass flux $Q_{m}$ in the T regime (labelled ‘a’ and ‘b’ in panels b,c). In scenario ‘a’, $Q_{m}$ decreases due to increasing mixing despite the volume flux $Q$ staying relatively constant ($2QF\sim F\sim\theta Re$). In scenario ‘b’, $Q_{m}$ decreases partly due to mixing, and partly due to a decrease in $Q$ (compensated by $F$ increasing faster than $\theta Re$). Accurate $Q$ and $Q_{m}$ data obtained by volumetric measurements of velocity and density in L18 (figure 5.12*(b)*) support scenario ‘b’ up to $\theta Re=132$, but additional data are required to draw general conclusions.

Finally, because this frictional hydraulic model assumes a two-layer flow without any form of mixing, it does not allow us to discuss the behaviour of the interfacial thickness $\delta$. We thus discuss several models that explicitly address mixing in the next section.

5.3 Mixing models

The importance and difficulty of modelling interfacial mixing in exchange flows has long been recognised (Helfrich, 1995; Winters & Seim, 2000). However, despite the existence of hydraulic models for multi-layered or continuously-stratified flows (Engqvist, 1996; Hogg & Killworth, 2004), to date there exists no ‘three-layer’ hydraulic model allowing for the exchange of momentum or mass between the two counter-flowing layers suitable to our problem (which would violate most hydraulic assumptions). Below we review some experimental, numerical, and theoretical work most relevant to the scaling of $Q_{m}$, $\delta$, and their relation to fundamental stratified turbulence properties such as diapycnal diffusivity and mixing efficiency.

5.3.1 Turbulent diffusion models

Cormack et al. (1974a) tackled natural convection in a shallow ($A\gg 1$) cavity with differentially heated walls. This problem is analogous to SID flows in the limit of maximum ‘interfacial’ thickness ($\delta=1$) in which turbulent mixing dominates to such an extent that the exchange flow is only weakly stratified in the vertical (i.e. $\langle|\partial_{z}\rho|\rangle_{z}<1$ because $|\rho(z=\pm 1)|<1$) and becomes stratified in the horizontal (i.e. $|\partial_{x}\rho(z=\pm 1)|>0$ and mean isopycnals are no longer horizontal). In their model, the horizontal hydrostatic pressure gradient is balanced only by a uniform vertical turbulent diffusion with constant $K_{T}$. Using the terminology of § 2.3, this balance could be called the hydrostatic-mixing (or ‘HM’) balance where ‘mixing’ plays a similar role to ‘viscosity’ in the ‘HV’ balance of § 2.3. Cormack et al. (1974a) solved this problem analytically and found: {subeqnarray} Q &= 5384(A K_T)^-1,

Q_m = 4 A K_T + 311451520(A K_T)^-3,

where we assumed a turbulent Prandtl number of unity for simplicity (i.e. the density equation has the same turbulent diffusivity). The above equations are adapted from equations (19) and (20) of Hogg et al. (2001) (in their review of the results of Cormack et al. (1974a)) to match our slightly different definitions of $Q,Q_{m},A$ and especially our definition of $K_{T}$ as being non-dimensionalised by the inertial scaling $\sqrt{g^{\prime}H}H/2$ (giving $K_{T}=(4Gr_{T})^{-1/2}$ where $Gr_{T}$ is their ‘turbulent Grashof number’). We also contrast the the uniform diffusivity $K_{T}$ in this model and the interfacial diffusivity $K_{I}$ in the frictional hydraulics model of § 5.2.3, which have different roles and different non-dimensionalisation ($\sqrt{g^{\prime}H}H/2$ vs $\nu$, hence ‘$K_{T}=K_{I}/Re$’). The predictions (5.3.1) were verified numerically and experimentally in two papers of the same series (Cormack et al., 1974b; Imberger, 1974), but only hold in the ‘high-mixing’ limit of $AK_{T}>\nicefrac{{1}}{{15}}$ below which inertia becomes noticeable and the assumptions start to break down (at $AK_{T}<\nicefrac{{1}}{{25}}$, $Q$ and $Q_{m}$ even exceed the hydraulic limit…).

Hogg et al. (2001) built on the above results and developed a model with linear velocity and density profiles within an interfacial layer of thickness $\delta<1$ and a uniform turbulent momentum and density diffusivity $K_{T}$. This models the ‘IHM’ balance, i.e. the transition between the Cormack et al. (1974a) $AK_{T}>\nicefrac{{1}}{{15}}$ high-mixing limit (the ‘HM’ balance where turbulent diffusion dominates over inertia, $\delta=1$, and (5.3.1) holds) and the $AK_{T}\rightarrow 0$ hydraulic limit (the ‘IH’ balance where inertia dominates over mixing, $\delta=0$, and $Q=Q_{m}=\nicefrac{{1}}{{2}}$ holds). Hogg et al. (2001) argued that $\delta$ would increase diffusively during the ‘duct transit’ advective timescale $A$, and integrated the linear velocity and density profiles across the interfacial layer to find

{subeqnarray}

δ&≈ 5(A K_T)^1/2,

Q ≈ 12-54(A K_T)^1/2, \slabelhogg-2

Q_m ≈ 12-53(A K_T)^1/2, \slabelhogg-3

where the prefactors $5,\nicefrac{{5}}{{4}},\nicefrac{{5}}{{3}}$ come from the imposed matching with the high-mixing solution (5.3.1). Hogg et al. (2001) validated these predictions with large eddy simulations and found good quantitative agreement for $Q,Q_{m},\delta$ across the range $AK_{T}\in[\nicefrac{{1}}{{2000}},\nicefrac{{1}}{{15}}]$, below which convergence to the inviscid hydraulic limit was confirmed.

In order to use these models to explain the scaling of $Q_{m}$ and $\delta$ with $A,B,\theta,Re,Pr$, we need to (i) extend them to the more complex ‘IHGM’ balance of SID flows in the I and T regimes in which gravitational forcing is present ($\theta>0$); (ii) have a model for the scaling of $K_{T}$ on input parameters (the above models prescribed $K_{T}$ as an input parameter, but it is a priori unknown in the SID). To do so, we propose to use insight gained by the energetics and frictional hydraulics models.

First, following the results of LPL19 and § 5.1 on the average rate of turbulent dissipation, it is tempting to model $K_{T}$ using a turbulence closure scheme like the mixing length or $K-\epsilon$ model, yet these require either a lenghtscale or the turbulent kinetic energy, which are both unknown (only the rate of dissipation is known, see (18)-(19)).

Second, following the frictional hydraulics results of § 5.2, we may conjecture that the ‘Reynolds stresses’ interfacial diffusivity $K_{I}$ in the I and T regimes may play a similar role to the uniform turbulent diffusivity in the present model. Recalling that by definition $K_{T}=K_{I}/Re$, combining the scalings (24)-(25) with (5.3.1) would suggest: {subeqnarray} δ∼12 - Q_m &∼ (A θ)^1/2       for Re ¡ 50A, \slabelwild-1

δ∼12 - Q_m ∼ (A^3 θ^2 Re^-1)^1/4  for  Re  ≫50A. \slabelwild-2

Unfortunately these scalings are not consistent with the observations of figures 5-8: $\delta$ is clearly a function of $Re$ for $Re<50A$ (less so at high $Re$ where the $A\theta$ scaling has indeed been observed by K91), and $\delta$ is clearly not a decreasing function of $Re$ for $Re\gg 50A$.

5.3.2 Previous mixing efficiency measurements and models

In this section we discuss two studies of the interfacial layer thickness $\delta$ and its relation to the Richardson number and mixing efficiency as a basis for the development of a more suitable model for SID flows in the next section.

Prastowo et al. (2008) studied exchange flows through short ($A\approx 2-3$), wide ($B\gg 1$), horizontal ($\theta=0$) contractions. Their measurements suggest an approximately constant interfacial thickness $\delta\approx 0.23-0.25$ across the range $Re\in[10^{4},10^{5}]$, in rough agreement with previously quoted estimates for shear-driven mixing flows (e.g. Sherman et al. (1978), p. 275 and references therein). They support this observation with ‘equilibrium’ or ‘marginally stable’ Richardson number arguments that the gradient Richardson number should be maintained near the Miles-Howard linear stability threshold, a phenomenon commonly observed subsequently in the observational literature on shear-driven mixing (Thorpe & Liu, 2009; Smyth & Moum, 2013). Assuming a linear profile for $u(z)$ and $\rho(z)$ across the mixing layer yields $Ri_{g}\approx\delta\approx 0.25$.

Prastowo et al. (2008) also measured the time-averaged mixing efficiency in their exchange flow using density profile measurements in the reservoirs at the end of the experiments, defined as $\mathcal{M}\equiv(P_{f}-P_{r})/(P_{i}-P_{r})\in[0,1]$, where $P_{i}$ the initial potential energy in the system (before the exchange flow starts), $P_{f}$ is the final measured potential energy in the system, and $P_{r}$ is the ‘reference’ or ‘minimum’ potential energy obtained by adiabatic (‘no-mixing’) rearrangement of fluid parcels from the initial conditions (i.e. $P_{i}-P_{r}$ is the initially available potential energy). They found collapse of the $\mathcal{M}$ data with $ARe$ and $\mathcal{M}\rightarrow 0.11$ for $ARe\rightarrow 10^{5}$ (using our notation). Finally, they supported this observation and linked $\mathcal{M}$ to $\delta$ by estimating mixing efficiency as the ratio of potential energy gain to kinetic energy deficit caused by the presence of a linear mixing layer, which yielded $\mathcal{M}\approx Ri_{g}/2\approx\delta/2\approx 0.125$.

Hughes & Linden (2016) studied horizontal lock exchange gravity currents, which for some of their life cycle, behave similarly to exchange flows. They measured $\delta\approx 0.33$ in the range $Re\in[10^{4},10^{5}]$. Using similar measurements to Prastowo et al. (2008), they found $\mathcal{M}\rightarrow 0.08$ asymptoting from below as $Re\rightarrow 10^{5}$. They supported this asymptotic value using a simple mixing model based on idealised linear profiles in the mixing layer, which yielded $\mathcal{M}=(2\delta^{2}/3)(1-2\delta/3)(1-\delta/2)^{-2}\approx\delta^{2}\approx 0.08$.

However, we have seen that exchange flows in inclined ducts have $\delta$ monotonically increasing not only with $A$ and $Re$, but also with $\theta$. In addition, much higher values of $\delta\gg 0.3$ (up to 0.8, and even 1 in K91) can be achieved even at moderate values of $\theta$ of a few $\alpha$ and $Re<10^{4}$. Therefore, the above models supporting values of $\delta=0.2-0.3$ and $\mathcal{M}=0.08-0.12$ in the T regime disagree with our data, despite (i) the similarity of SID flows to the flows assumed above (shear-driven mixing flows with the same ‘IH’ velocity scaling $-1\lesssim u\lesssim 1$) and (ii) the fact that these models would apparently not be modified by the presence of gravitational forcing ($\theta>0$).

5.3.3 New mixing efficiency model

To address this, we propose a different toy model of mixing based on the energetics framework of LPL19. As sketched in figure 11(a), we consider that the duct is composed of three volume-averaged energy reservoirs (in bold): potential energy $P$, kinetic energy $K$, and internal energy $I$ (heat). We further decompose the potential energy reservoir into an available potential energy $P_{A}$, and a background potential energy $P_{B}$ (such that $P=P_{A}+P_{B}$), as is customary in the study of mixing (see e.g. Winters et al. (1995)).

As explained in LPL19 (see their § 4.1-4.3 and figure 8*(b)*), forced flows have, to a good approximation, the following steady-state energetics: the external fluid reservoirs provide an advective flux of potential energy into the duct, which we identify here as being an advective flux of available potential energy $\Phi^{\textrm{adv}}_{P_{A}}\approx Q_{m}\theta/8$, which is then converted to kinetic energy by the horizontal buoyancy flux $B_{x}$, and to heat by the viscous dissipation $D\approx(2/Re)\langle\mathsf{s}^{2}\rangle_{x,y,z,t}$. When turbulent mixing is neglected, all these fluxes have equal magnitude, and $D\approx(1/8)Q_{m}\theta$. When turbulent mixing is included, a net vertical buoyancy flux $B_{z}$ converts part of $K$ back to $P_{A}$, and a net irreversible diapycnal flux $\Phi^{d}$ converts part of $P_{A}$ to $P_{B}$, at a steady-state rate equal to the advective flux of $P_{B}$ out of the duct, back into the external reservoirs $|\Phi^{\textrm{adv}}_{P_{B}}|=|\Phi^{d}|$. The mixing efficiency quantifies the percentage of total time- and volume-averaged power throughput $\Phi^{\textrm{adv}}_{P_{A}}$ that is spent to irreversibly mix the density field inside the duct

[TABLE]

It is expected that $\mathcal{M}\ll 1$ in such flows, as represented by the respective thickness of the arrows in figure 11(a), representing the order of magnitude of the fluxes.

As sketched in figure 11(b), we propose piecewise-linear flow profiles $u(z)=\rho(z)$ at either end of the duct as a minimal model to estimate the magnitude of $\Phi^{\textrm{adv}}_{P_{B}}$ as a function of the interfacial layer thickness $\delta$, and eventually link it to input parameters $A,\theta,Re$. We consider that fluid comes from the external reservoirs into the duct unmixed below (resp. above) the interfacial mixing layer at the left (resp. right) end of the duct, and leaves the duct mixed with a linear profile, going from 0 at the bottom (resp. top) edge of the mixed layer to $-1$ (resp. 1) at the top (resp. bottom) edge of the mixed layer. (In more central sections of the duct, mixing smoothes out the discontinuities at the edges of the mixing layer present at the ends, and we expect the continuous linear profile drawn in the centre, but it is irrelevant to the following calculations). The outflow of mixed fluid creates the following net flux of background potential energy out of the duct:

[TABLE]

where $|_{L-R}$ denotes the difference between the values at the left and right boundary, and the prefactor $1/(4A)$ comes from the non-dimensionalisation of the energy budget equations (see LPL19, equation (4.14a)). From (26)-(27) and $\Phi^{\textrm{adv}}_{P_{A}}\approx Q_{m}\theta/8$, we now deduce:

[TABLE]

Encouragingly, this estimation has the potential to be consistent with our data in the SID. Assuming that $Q_{m}\approx 0.5$ throughout most of the I and T regimes, we conjecture that most of the dependence on $Re$ observed in the $\delta$ data of figure 7 is due to the underlying monotonic increase of $\mathcal{M}(Re)$, which is a priori unknown, but consistent with the observations of Prastowo et al. (2008) and Hughes & Linden (2016). The observation of K91 and figure 7(a-b) that $\delta$ primarily scales with the group $A\theta$ at $Re\gg 500A$ (as sketched in figure 3(c)) would suggest that $\mathcal{M}$ asymptotes to a constant value at high $Re$, which is also consistent with the observations of Prastowo et al. (2008) and Hughes & Linden (2016) at $ARe>10^{5}$ and $Re>10^{5}$ respectively. Assuming their high-$Re$ asymptotic value of $\mathcal{M}\approx 0.1$, we obtain

[TABLE]

This gives, for example, $\delta\approx 0.4$ when $\theta/\alpha\approx 1/2$. This value agrees with the K91 data (see § 3.3) and our LSID data (figure 7(a), at $Re>10^{4}$ and $ARe>10^{5}$). However, this value does not agree well with our HSID, tSID, and mSID $\delta$ data (figure 7(b-d), in which $\delta$ remains dependent on $Re$. This is presumably due to the lower values of $A$ and/or $Re$ in these data sets, which remain below the asymptotic values of $Re>10^{5}$ and $ARe>10^{5}$. In other words, we believe that our $\delta$ data and (29) are consistent and provide further (albeit indirect) evidence for the monotonic increase of $\mathcal{M}$ with $Re$.

6 Conclusions

6.1 Problem and approach

In this paper, we investigated buoyancy-driven exchange flows taking place in inclined rectangular ducts (figure 1). We focused on the behaviour of three key dependent variables: the qualitative flow regime (laminar, wavy, intermittently turbulent, or fully turbulent), the non-dimensional mass (or buoyancy) flux $Q_{m}$, and the non-dimensional thickness of the interfacial layer $\delta$ as the five non-dimensional input parameters were varied: the duct longitudinal aspect ratio $A$, spanwise aspect ratio $B$, tilt angle $\theta$, Reynolds number $Re$, and Prandtl number $Pr$.

Dimensional analysis (figure 2) and the experimental literature (figure 3 and table A) showed that the rich dynamics of these sustained stratified shear flows were accessible for a wide range of $Re$ and for $\theta$ of at most a few duct aspect ratios $\alpha=\tan^{-1}(A^{-1})$. Our focus on ‘long’ ducts ($A\gg 1$) allowed us to explore these dynamics while keeping $\theta$ small enough to remain relevant to (largely-horizontal) stably-stratified geophysical flows and turbulence, which are our ultimate motivation.

To overcome the limitations of past studies of the problem, we presented extensive experimental results for all three variables of interests (regimes, $Q_{m}$, and $\delta$) in the $(\theta,Re)$ plane for five different data sets, between which $(A,B,Pr)$ were varied systematically (table 1).

6.2 Experimental results

First, our data (figures 5-8) confirmed the conclusions of past studies: that increasingly disordered and turbulent regimes are found as $A,\theta,Re$ are increased, that $Q_{m}$ is non-monotonic in $\theta$ and $Re$, and that $\delta$ is monotonic in $A,\theta,Re$. Second, our data confirmed beyond doubt the previously-supposed existence and importance of at least one additional elusive non-dimensional input parameter related to the dimensional height of the duct $H$, because our regime, $Q_{m}$, and $\delta$ data at the same $A,B,\theta,Re,Pr$ but different $H$ do not collapse. Third, our data highlighted the non-trivial dependence of all variables on all five parameters $A,B,\theta,Re,Pr$. Regime transition, iso-$Q_{m}$, and iso-$\delta$ curves are not only shifted in the $(\theta,Re)$ plane at different $A,B$, or $Pr$, but they also generally exhibit different power law scalings in $\theta$ and $Re$ at different $A,B,Pr$.

Given the breadth of our observations summarised above, and the relative richness of our data in the $(\theta,Re)$ plane compared to the few values of $A,B,Pr$ studied, we focused specifically on the very last observation above, i.e. on the various scalings of the form $\theta Re^{\gamma}=$ const. governing the regime transitions curves and the major axis of hyperbolas best fitting $Q_{m}$ in the $(\log\theta,\log Re)$ plane. Even within this specific focus, we discovered that $\gamma$ not only varies between data sets (at different $A,B,Pr$) but that it also varies within a given data set (at fixed $A,B,Pr$): (i) $\gamma$ is generally different for the regime data ($\gamma=1$ or 2) and the $Q_{m}$ data ($0.3<\gamma<2.6$) implying that regime and $Q_{m}$ are not well correlated (whereas regimes and $\delta$ are); and (ii) in one regime data set, $\gamma$ even takes two different values (of 1 and 2) in different regions of the $(\theta,Re)$ plane.

6.3 Modelling results and outlook

To provide a modelling framework, particularly to understand the above observation (ii), we split the $(\theta,Re)$ plane into four quadrants delimited by $\theta=\alpha$ (the ‘lazy/forced flow boundary’, based on the respective dominance of hydrostatic/gravitational forcing) and $Re=50A$ (the ‘low/high $Re$ boundary’, based on whether or not wall boundary layers are fully developed across the duct or remain thin). To explain the above observations (i)-(ii), we discussed three families of candidate models.

In § 5.1 we considered the volume-averaged energetics framework of Lefauve et al. (2019) (LPL19). LPL19 physically explained the $\theta Re=$ const. scaling of regime transitions of forced ($\theta\gtrsim\alpha$), low-Re ($Re\lesssim 50A$), salt-stratified ($Pr=700$) flows as being caused by threshold values of the three-dimensional kinetic energy dissipation (equation (18)). We extended their argument to high-$Re$ ($Re\gg 50A$) flows, by accounting for the scaling of two-dimensional, laminar boundary layer dissipation. However, the resulting predicted scaling in $\theta Re-A^{-1/2}Re^{1/2}$ (equation (19)) did not agree with any of our regime data. Detailed measurements of dissipation in these high-$Re$ flows (not found in LPL19) would be valuable to understand why this is the case, but they are very challenging to perform due to the required spatio-temporal resolution.

In § 5.2 we developed a two-layer frictional hydraulics model of SID flows (figure 9) from Gu & Lawrence (2005) and showed that the existence of a solution imposed a lower and upper bound on the product of the volume flux by a parameter quantifying wall and interfacial friction (equation (23)). This model explained the qualitative behaviour of $Q_{m}$ with $\theta Re$, and the fact that regimes and $Q_{m}$ could have different scalings (figure 10). This model also provided a quantitative scaling for the interfacial friction parameter and, in turn, for regime transitions, based on our conjecture that regime transitions were directly linked to interfacial turbulent stresses. Although the resulting low-$Re$ scaling in $\theta Re$ (equation (24)) was identical to that predicted by the energetics model and correct (at least for $Pr=700$), the high-$Re$ scaling in $A^{1/2}\theta Re^{1/2}$ (equation (25)) did not agree with our regime data.

Neither the energetics nor the frictional hydraulics model could predict the observed scalings in $\theta Re$ or $\theta Re^{2}$ observed in lazy flows ($\theta\lesssim\alpha$) because these flows are under-specified in either model (they have more unknowns than equations). In addition, scalings laws deduced from plots in the $(\log\theta,\log Re)$ plane break down for lazy flows at slightly negative angles ($-\alpha\lesssim\theta\lesssim 0$), which we largely ignored in this paper.

In § 5.3, we focused on the scaling of $\delta$ by modelling turbulent mixing. We first considered a model with constant turbulent diffusivity imposed throughout the domain (Cormack et al., 1974a; Hogg et al., 2001). We attempted to link this diffusivity to input parameters following insights gained from frictional hydraulic theory, but the resulting scalings (equation (5.3.1)) did not agree with our data. We then explained why previous measurements and models of $\delta$ in related stratified shear flows (Prastowo et al., 2008; Hughes & Linden, 2016) were inconsistent with our results on exchange flows in inclined ducts. We thus developed a new model that explicitly represents the rate of mixing in the energy budget analysis of LPL19, and quantifies this mixing as a function of known input parameters and an unknown mixing efficiency $\mathcal{M}$ using a simplified flow profile (figure 11, equation  (27)). The resulting expression for $\delta$ (equations (28)-(29)) is qualitatively consistent with our observations, but it involves $\mathcal{M}$ (not measured in these experiments) whose scaling on $Re$ is critical. Our model and data indirectly support previous observations of Prastowo et al. (2008) and Hughes & Linden (2016) that $\mathcal{M}$ monotonically increases with $Re$ to reach asymptotic values of $\mathcal{M}\approx 0.1$ at very high $Re$, but direct measurements of $\mathcal{M}$ would be desirable to confirm this.

While these models have allowed us to make significant progress by providing useful physical insights and partial quantitative results regarding scaling laws in $A,\theta,Re$, our experimental observations have raised an even larger number of questions which remain open. Among these are the elusive existence of an important sixth non-dimensional input parameter linked to the duct height $H$, the influence of the spanwise aspect ratio $B$ and Prandtl number $Pr$, and the scaling of the mixing efficiency $\mathcal{M}$.

Acknowledgements

AL is supported by an Engineering and Physical Sciences Research Council (EPSRC) Doctoral Prize Fellowship. We also acknowledge funding from EPSRC under the Programme Grant EP/K034529/1 ‘Mathematical Underpinnings of Stratified Turbulence’ (MUST) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation Grant No 742480 ‘Stratified Turbulence And Mixing Processes’ (STAMP). We thank Colin Meyer and Simon Vincent for collecting the regime and mass flux data in the LSID and HSID geometries in 2012-2014, and we thank Dr Ling Qu visiting from the South China Sea Institute of Oceanology (SCSIO, Guangzhou, China) for collecting part of the regime and mass flux data in the mSID and tSID in 2015-2016. We also thank Dr Darwin Kiel from Coanda Research and Development (Burnaby, Canada) for his reading of the manuscript, his suggestions, and his permission to reproduce some of his results. Finally, we are grateful for the invaluable experimental support of Prof Stuart Dalziel, Dr Jamie Partridge and the technicians of the G. K. Batchelor Laboratory.

Appendix A Summary table of the literature review

Summary of the literature review. For each paper, we specify the scale of the apparatus (duct height or pipe diameter $H$), the parameters that were either fixed or whose variation was not studied, those that were varied and studied, the key conclusions about the scaling of transitions between flow regimes (based on empirical or physical arguments), mass flux $Q_{m}$, and interfacial layer thickness $\delta$.

Name $H$ (mm) Fixed params. Varied params. Regimes Mass flux $Q_{m}$ Interfacial thickness $\delta$

MR61 150

$A=6$

$B=\nicefrac{{1}}{{3}}$

$Re\approx 2\times 10^{3}$

$Pr=700$

$\theta\in$ ?

$\sim Re\,Q^{2}$

(empirical)

$-$ $-$

LT75

$150-1140$

(pipe)

$B=1$

$\theta=0^{\circ}$

$A\in[0.5,20]$

$Re\in[15,75]\times 10^{3}$

$Pr=1,\,700$

$-$

Constant in $A,Re,Pr$

at $\theta=0^{\circ}$

$Q_{m}\approx 0.23$

$-$

MT75

$50$, $150$

(pipe)

$B=1$

$Pr=700$

$A\in[0.5,18]$

$Re=2,\,20\times 10^{3}$

$\theta\in[0^{\circ},90^{\circ}]$

$-$

Non-monotonic in $A,\,\theta$

$Q_{m}\approx 0.3$ at $\theta=0^{\circ}$

$Q_{m}\approx 0.4$ at $\theta\approx\nicefrac{{\alpha}}{{2}}$

$Q_{m}\rightarrow 0$ as $\theta\gg\alpha$

Dependence on $Re$

$-$

W86

$14-25$

(pipe)

$B=1$

$\theta=0^{\circ}$

$Pr=700$

$A=1.6,\,3.5,\,9.6$

$Re\in[0.2,3]\times 10^{3}$

$\sim Re$

(empirical)

Monotonic in $A^{-1}Re$

$Q_{m}\approx 0.15$ at $Re\approx 20A$

$Q_{m}\approx 0.35$ at $Re\approx 500A$

$-$

K91 50, 100

$B=1,\,2,\,4$

$Re\in[2,15]\times 10^{3}$

$Pr=700$

$A=1,\,2,\,4,\,8$

$\theta\in[-45^{\circ},90^{\circ}]$

$\sim A\,\tan\theta$

(empirical)

Non-monotonic in $A,\,\theta$

Collapse with $A\tan\theta$

Independence on $Re$

Monotonic in $A,\,\theta$

Max. $\delta=1$ at $\theta\approx 2\alpha$

Collapse with $A\tan\theta$

ML14 100

$B=1$

$Pr=700$

$A=15,30$

$\theta\in[-1^{\circ},4^{\circ}]$

$Re\in[1,20]\times 10^{3}$

$\sim A\,\theta\,Re^{2}$

(empirical)

Monotonic in $\theta$

$Q_{m}\rightarrow 0.5$ as $\theta>\nicefrac{{\alpha}}{{2}}$

$-$

LPL19 45

$A=30$

$B=1$

$Pr=700$

$\theta\in[-1^{\circ},6^{\circ}]$

$Re\in[0.3,6]\times 10^{3}$

$\sim\theta\,Re^{2}$ as $\theta\lesssim\alpha$

(empirical)

$\sim\theta\,Re$ as $\theta\gtrsim\alpha$

(physical)

Non-monotonic in $\theta,\,Re$ $-$

Appendix B Unpublished data from Kiel (1991)

In figure B.1 we reproduce some of the unpublished mass flux and interfacial thickness data in Kiel (1991) (K91) discussed in § 3.2-§ 3.3. In panel (a) K91 observed relative independence of $Q_{m}$ on $Re$ at large enough values $Re>400A=1600$ (his figure 4.10). In panel (b) we show his reproduction (his figure 2.6) of the $Q_{m}(A,\theta)$ in Mercer & Thompson (1975) (MT75, figure 5), and in panels (d,f), his $Q_{m}(A,\theta)$ and $\delta(A,\theta)$ data (his figures 5.2 and 5.13 respectively). In panels (c,e,f) we show the collapse with his proposed $Ri_{G}$ criterion (15).

Appendix C Experimental methodology

C.1 Flow regimes

Regimes were largely determined by shadowgraph observations over a subsection of the length of the duct, following the qualitative description of each regimes provided by ML14 (see their § 3.1 and figure 3). A schematic of the shadowgraph setup can be found in L18, § 2.1.

In the mSID data set, 48 out of 360 regime identifications were not made by shadowgraph, but rather by direct visualisation of the density field by planar laser induced fluorescense (PLIF), since more detailed measurements of the velocity and density fields (incompatible with simultanenous shadowgraph) have been performed in this geometry (Lefauve et al., 2018; Partridge et al., 2019; Lefauve et al., 2019).

All raw video data, including those obtained by other experimenters (acknowledged at the end of the paper), were reprocessed in an effort to ensure that regimes were identified as consistently as possible across all five data sets of table 1 (especially in the cases where the distinction between regimes can be subtle).

C.2 Mass flux

Mass fluxes were determined, as in ML14, by measuring the average initial (‘i’) and final (‘f’) density in each reservoir: reservoir ‘1’, initially at density $\rho_{1}^{i}=\rho_{0}+\Delta\rho/2$ and finally at a well mixed density $\rho_{1}^{f}$ and ‘2’, initially at $\rho_{0}-\Delta\rho/2$ and finally at $\rho^{f}_{2}$, giving the following two estimations

[TABLE]

where $V_{1},V_{2}$ are the (typically approximately equal) volumes of fluid in the respective reservoirs, and the tilde on $\tilde{Q}_{m}$ stresses the fact that they are non-dimensional (despite all quantities on the right side of the $=$ sign being dimensional). Experiments in which both estimates differed by more than $(Q_{m,1}-Q_{m,2})/(Q_{m,1}+Q_{m,2})>10~{}\%$ were rejected (typically due to an initial misadjustement of the free surfaces resulting in a net volume flux $\langle u\rangle_{x,y,z,t}\neq 0$). All data shown in this paper thus have near-zero net volume flux, and we only use the average value $Q_{m}\equiv(Q_{m,1}+Q_{m,2})/2$, .

We recall that $T$ in (30) is the (dimensional) duration of an experiment. The determination of the relevant $T$ was made carefully but remains subject to intrinsic uncertainties which affect $Q_{m}$ as we explain next. The duct is opened at time $t^{a}$ initiating a gravity current lasting until the exchange flow is considered fully established by shadowgraph visualisations at time $t^{b}$. The exchange flow of interest continues until the levels of the discharged fluids approach the ends of the duct, at which point one end of the duct is closed at time $t^{c}$, shortly before the other end of the duct is closed at $t^{d}$. To avoid under- and over-estimations of $Q_{m}$ by the intervals $t^{d}-t^{a}$ and $t^{c}-t^{b}$ (respectively), we choose to use the average of the two $T=(t^{d}-t^{a}+t^{c}-t^{b})/2$, and to use error bars to indicate the magnitude of the resulting uncertainty (the difference between the over- and under-estimation). Note that error bars tend to be larger at high $Re$ (figure 7) because the overall duration $T$ of an experiment is inversely proportional to the magnitude of the dimensional exchange velocities (scaling with $\sqrt{g^{\prime}H}$, and hence with $Re$) due to the finite size of the reservoirs. A smaller duration $T$ increasing the relative duration of initial transients (typically fixed) and therefore the uncertainty about $T$.

Note that measurements of $Q_{m}$ in temperature-stratified experiments (mSID ​ ​data set) could not be performed due to the practical impossibility to control the heat loss occurring through the boundaries of the reservoir and the free surface.

For more details on these measurements, see L18, § 2.2.

C.3 Interfacial layer thickness

The thickness of the interfacial density layer $\delta$ was estimated from shadowgraph images. To a reasonable approximation, the refraction of near-parallel light beams by inhomogeneities in the density field results in a recorded greyscale light intensity $I(x,z)$ proportional to the second vertical derivative of the density field integrated in the spanwise direction $I(x,z)\propto\int_{-B}^{B}\partial_{zz}\rho\,\mbox{d}y$ (for a full derivation and discussion of the approximations, see L18, § 2.1). This makes shadowgraphy particularly well-suited to detect the average location of large-scale curvatures in the density field, which are precisely the edges of the interfacial density layer.

Due to the nature of shadowgraph images, and to its sensitivity to air bubbles or scratches on the walls of the reservoirs, the identification of minima and maxima of $I(z)$ could only be semi-automated according to the following strict methodology, illustrated in figure C.1:

(i)    

A random sample of typically three to five snapshots per movie were selected and averaged (although in rare cases only one still image was available);

(ii)   

A randomly-generated location in the streamwise direction was selected (dotted blue lines) and they greyscale intensity profile $I(z)$ at this particular $x$ location was superimposed onto the image (solid red curves);

(iii)   

The profile $I(z)$ was carefully interpreted, and the local extrema representing the top and bottom duct boundaries (yellow crosses) and edges of the interface to measure (circles) were carefully selected by a click.

(iv)   

The ratio of pixel distances between the selected edges of the density interface and the top and bottom walls was computed to yield $\delta$.

All images were processed by the first author to ensure consistency, and yielded a total of 351 values of $\delta$ for all four duct geometries (table 1).

Note that measurements of $\delta$ in temperature-stratified experiments (mSID ​ ​data set) could not be performed since the refractive index of water is a weaker function of temperature than salinity at comparable density differences, resulting in insufficient contrast and thus noisier $I(z)$ (though sufficient to determine the flow regime).

Appendix D Frictional two-layer hydraulic model

In this section we give details of the two-layer frictional hydraulic model introduced in § 5.2 and sketched in figure 9. This model is based on Gu (2001); Gu & Lawrence (2005) but includes non-zero tilt angles and a wider range of frictional stresses suited to the SID. We cover the model formualtion in § D.1, the parameterisation of frictional effects in § D.2, and the solution to the full problem in § D.3.

D.1 Model formulation

The frictional hydraulic model appears at first inconsistent because it is based on velocities that are uniform in the cross-sectional plane ($\partial_{y,z}u_{1,2}=0$), while implicitly acknowledging and parameterising the effects of viscous stresses resulting from $\partial_{y,z}u_{1,2}\neq 0$. This model is however consistent provided that the departure from hydrostaticity is small (vertical and spanwise accelerations are negligible) and that viscous stresses are localised in relatively narrow boundary layers at the walls and interface ($Re\gg 50A$), rather than rather than through the whole volume ($Re<50A$).

Following standard hydraulic practice, the effective ‘hydraulic’ velocities $u_{1,2}(x)$ that will be used to compute the total Bernoulli head (kinetic energy) of each layer need to be defined in a way that accounts for the non-uniformity of the underlying ‘real’ velocity profile in the SID $u(x,y,z)$

[TABLE]

where $\langle\cdot\rangle_{z_{1,2}}$ denotes averaging over the lower/upper layer ($z\in[-1,\eta]$ and $z\in[\eta,1]$ respectively, see figure 9(a)), and the velocity distribution coefficient $\lambda_{1,2}$ (also called kinetic energy correction coefficient or Coriolis coefficient) is defined as

[TABLE]

respectively in the lower and upper layer (see e.g. Chow, 1959, § 2.7-2.8 and Chanson, 2004, § 3.2.2). The greater the non-uniformity of the velocity profile $u$, the larger $\lambda$ is. For the SID flows considered in this paper, volumetric velocity measurements showed that $\lambda$ varies over a relatively relatively small range $1<\lambda\lesssim 2$ (see L18, § 5.5.2). To simplify the following discussion, and since the effects of $\lambda$ are not central here (they quantitative rather than qualitative), we make the approximation that $\lambda_{1,2}(x)\approx 1$, effectively assuming that $u_{1,2}(x)=\langle u(x,y,z)\rangle_{y,z_{1,2}}$ in the following.

First, the conservation of Bernoulli potential in two-layer hydraulic flows is commonly expressed using the so-called ‘internal energy’ of the system

[TABLE]

Second, the conservation of volume and zero-net flux conditions are expressed all along the duct as

[TABLE]

The third important ingredient of two-layer hydraulics is the condition of hydraulic control, which requires that the composite Froude number $G$ is unity at sharp changes in geometry, i.e. at the duct ends (Armi, 1986; Lawrence, 1990):

[TABLE]

where the second equality uses (34) and the third equality is the control condition.

In horizontal, frictionless ducts, $E(x)=0$, hence $\eta=0$ and $u_{1}=-u_{2}=Q=\nicefrac{{1}}{{2}}$ all along the duct.

When the combined effects of a small positive tilt angle $\theta>0$ and frictional stresses are added, the slope of the internal energy becomes

[TABLE]

(this is the two-layer equivalent of single-layer ideas found in (Henderson, 1966, § 4.4-4.5)). By analogy with the topographic slope $\theta$, the ‘frictional’ slope $S(x)$ is computed by a balance of all the stresses acting on an infinitesimal slice of thickness $dx$ (figure 9(b)):

[TABLE]

The subscript $i=1,2$ represents respectively the bottom and top layers, the superscript $j=Z,Y,I$ represents the origins of the stresses in the model: top and bottom wall stresses ($Z$, shown in blue in the figure), side wall stresses ($Y$, in green) and interfacial stresses ($I$, in red), $A_{i}^{j}$ represents the surface area over which the respective stresses act, and $V_{i}$ the volume of each layer. Note that the interfacial stresses have equal magnitudes on either sides of the interface $|\tau_{1}^{I}|=|\tau_{2}^{I}|\equiv\tau^{I}$. Following figure 9(b) and after elementary algebra, the balance in (37) can be rewritten as:

[TABLE]

where all the stresses in this equation and henceforth are norms and have positive values. For further details about the development of this model from first principles, see L18, § 5.2.

D.2 Parameterisation of shear stresses

We now tackle the relation between the stresses $\tau_{i}^{j}$ and the underlying ‘real’ flow profiles $u(x,y,z)$. We start by considering the bottom wall stress of the lower layer $\tau_{1}^{Z}$ in order to introduce the key concepts and definitions, before extending them to the other stresses. Using non-dimensional variables for $\tau_{1}^{Z}$ and $u(x,y,z)$, we first write the dimensional equation for this stress as a simple function of the local shear

[TABLE]

where the $\Delta U/2$ and $H/2$ factors come from non-dimensionalising $\tau_{1}^{Z},u,z$, and simplify to

[TABLE]

In order to correctly parameterise $\tau_{1}^{Z}(x)$ and all other relevant stresses using well-defined, constant friction coefficients, we follow the following five steps.

First, we define the cross-sectional ‘shape’ ($y-z$ dependence) of the local velocity profile in the lower layer as

[TABLE]

such that $\langle\hat{u}_{1}(x,y,z)\rangle_{y,z_{1}}=1$. This decomposition allows to rewrite (40) as

[TABLE]

which is an exact expression for the local shear stress that does not require any assumptions about the value of the velocity gradient or flow profile. 2. 2.

Second, we define a ‘layer-rescaled’ coordinate $\hat{z}_{1}$ as

[TABLE]

in which layer 1 always has thickness one ($\hat{z}_{1}\in[-1,0]$), giving us

[TABLE] 3. 3.

Third, we define a constant, bottom friction parameter $f_{Z_{1}}$ to parameterise the stress:

[TABLE]

We note that despite the rescaling of $u(x,y,z)$ by $u_{1}(x)$ and the stretching of $z$ to $\hat{z}_{1}$ such that the interface is located at $\hat{z}_{1}(x)=0$, $\hat{u}_{1}$ still has a weak residual $x$ dependence. Since for simplicity, we choose to model $f_{Z_{1}}$ as independent of $x$, the velocity gradient $\p\hat{u}_{1}(x,y,z)/\p\hat{z}|_{\hat{z}_{1}=-1}$ must now technically be averaged not only over $y$ but over $x$ and $y$, as shown in (45). We also note that the $u_{1}^{2}(x)/Q$ factor in (45) results from the product of $u_{1}(x)$ (by definition of $\hat{u}_{1}$) by $u_{1}(x)/Q$ (by definition of $\hat{z}_{1}$). Physically, this quadratic dependence corresponds to the vertical shear being enhanced not only by the magnitude of $u_{1}$, but also by the enhanced vertical gradient due to the thinner layers where $u_{1}$ is larger. This $u_{1}^{2}(x)/Q$ scaling will be found in the interfacial stress $\tau^{I}$ too. However, the equivalent formulation to (44) for the side wall stress in layer 1, $\tau_{1}^{Y}$, is

[TABLE]

where we assume identical shear at $y=\pm 1$. We emphasise that since the $y$ derivative does not experience any rescaling due to the layer thickness, it follows a $u_{1}(x)$ scaling (as opposed to $u^{2}_{1}(x)/Q$ for $z$ derivatives). 4. 4.

Fourth, we generalise the above definitions of $\hat{u}_{1}$ and $\hat{z}_{1}$ to both layers by defining a global $\hat{u}$ as

[TABLE]

and a global $\hat{z}$ as

[TABLE] 5. 5.

Fifth, we consider the role of turbulence at the interface, caused by Reynolds stresses which we parameterise, by analogy with (40), as follows

[TABLE]

where $\textbf{u}^{\prime}\equiv\textbf{u}-\langle\textbf{u}\rangle_{t}$ is the perturbation around the temporal mean and $K_{I}$ the turbulent momentum diffusivity non-dimensionalised by the molecular value $\nu$. Under these conditions, the total (molecular $+$ turbulent) interfacial stress $\tau^{I}$ can be expressed precisely as:

[TABLE]

where $\hat{z}_{I}$ denotes averaging over the interfacial mixed layer.

Based on the five above steps, we propose the following parameterisation of frictional effects in the hydraulic model {subeqnarray} τ_1,2^Z(x) & = fZRe  u1,22(x)Q, \slabeldef-fz

τ_1,2^Y(x) = fYRe  —u_1,2(x)—, \slabeldef-fy

τ_1^I (x) = τ_2^I(x) = fIRe  (u1(x)-u2(x))2Q.\slabeldef-fi

where the vertical, spanwise and interfacial friction parameters are, respectively, {subeqnarray} f_Z &≡ ⟨— \p^u(x,y,^z)\p^z—_^z=±1 ⟩_xy,

f_Y ≡ ⟨— \p^u(x,y,^z)\py—_y±1 ⟩_xz,

f_I ≡ (1+K_I) ⟨— \p^u(x,y,^z)\p^z— ⟩_x,y,z_I.

The $y$ and $\hat{z}$ derivatives at $y,\hat{z}=\pm 1$ should be very similar, and the average of the two is implied. The three parameters can be computed from three-dimensional, three-component velocity measurements, as was done in L18, § 5.5.

D.3 Key equations and solution method

We can now rewrite the frictional slope $S(x)$ in (38) using (D.2) and (34) as

[TABLE]

By combining this expression for $S(x)$ with the expression for the composite Froude number $G^{2}(x)$ in (35) we finally obtain an expression for the differential equation governing the evolution of the interfacial slope $\eta^{\prime}(x)$ in (36)

[TABLE]

where the spanwise friction ratio $r_{Y}$ and interfacial friction ratio $r_{I}$ are defined under (21). This equation was simplified to (20) for the discussion in § 5.2.

The idea behind the solution to this kind of problem can essentially be found in Gu & Lawrence (2005). However, contrary to their model (which had no tilt angle and no top and side wall friction $\theta=f_{Z_{1}}=f_{Y}=0$), our model does not allow us to find an analytical solution to (52). We must therefore resort to an iterative numerical approach which we briefly outline below.

By symmetry of the problem (guaranteed under the Boussinesq approximation), $\eta$ is an odd function of $x$. We impose the boundary condition $\eta(0)=0$ and need only solve (52) in half of the domain (say $x\in[0,A]$).

However, since the volume flux $Q$ in (52) is a priori unknown, we must solve a coupled problem imposing the additional condition of hydraulic control at each duct end (denoted by the superscript ∗)

[TABLE]

where $\eta^{*}$ is the result of the forward integration of (52)

[TABLE]

The coupled problem for $\eta(x)$ and $Q$ for any given set of forcing and friction parameters $(\theta,Re,f_{Z},r_{Y},r_{I})$ can then be solved by the following iterative algorithm (illustrated in L18, figure 5.4).

(i)    

Guess $Q$;

(ii)   

Integrate numerically (52) from $x=0$ to $A$ to get $\eta^{*}$ as in (54);

(iii)   

Get the $Q$ corresponding to this $\eta^{*}$ by the criticality condition (53);

(iv)   

Compare this $Q$ with the initial guess and update the guess;

(v)    

Repeat until convergence of $Q$.

This model and its solution were validated using parameters $(\theta,Re,f_{Z},r_{Y},r_{I})$ from an experiment in the L regime, and quantitative agreement with $\eta(x)$ and $Q$ measurements was found L18, § 5.5.3.