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 vice versa, 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 scalar 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 Reference Deacon1971, 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 Reference Deacon1971, chap. 7)). More recently, it has been recognised that nutrient transport from the Atlantic partially supported primary production in Mediterranean ecosystems (Estrada Reference Estrada1996). 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.
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 Reference Linden1999). 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 (Reference Reynolds1883, § 12) and Taylor (Reference Taylor1931). Multi-layered stratified shear flows have complex hydrodynamic stability and turbulent mixing properties (Caulfield Reference Caulfield1994; Peltier & Caulfield Reference Peltier and Caulfield2003). 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, Caulfield & Peltier Reference Salehipour, Caulfield and Peltier2016; Lefauve et al. Reference Lefauve, Partridge, Zhou, Caulfield, Dalziel and Linden2018); spatio-temporal turbulent intermittency (de Bruyn Kops Reference de Bruyn Kops2015; Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016; Taylor et al. Reference Taylor, Deusebio, Caulfield and Kerswell2016); and layering and mixing (Salehipour & Peltier Reference Salehipour and Peltier2015; Lucas, Caulfield & Kerswell Reference Lucas, Caulfield and Kerswell2017; Zhou et al. Reference Zhou, Taylor, Caulfield and Linden2017; Salehipour, Peltier & Caulfield Reference Salehipour, Peltier and Caulfield2018).
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 current state of knowledge on the behaviour of our three key variables 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 summarise and conclude in § 6.
2 The experiment
2.1 Set-up 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 $\unicode[STIX]{x1D70C}_{0}\pm \unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$, connected by a long rectangular duct that can be tilted at an angle $\unicode[STIX]{x1D703}$ 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 input parameters, as we discuss later in the paper.
Our notation is shown in figure 1(b) and largely follows that of Lefauve et al. (Reference Lefauve, Partridge, Zhou, Caulfield, Dalziel and Linden2018), Lefauve, Partridge & Linden (Reference Lefauve, Partridge and Linden2019). The duct has length $L$, height $H$ and width $W$. The streamwise $x$ axis is aligned along the duct and the spanwise $y$ axis is aligned across the duct, making the $z$ axis tilted at an angle $\unicode[STIX]{x1D703}$ from the vertical (resulting in a non-zero streamwise projection of gravity $g\,\sin \,\unicode[STIX]{x1D703}$). The angle $\unicode[STIX]{x1D703}$ 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 $\boldsymbol{u}(x,y,z,t)=(u,v,w)$ along $x,y,z$, and the density field is $\unicode[STIX]{x1D70C}(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
(2.1)$$\begin{eqnarray}Q\equiv \langle |u|\rangle _{x,y,z,t},\end{eqnarray}$$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|\,\text{d}x\,\text{d}y\,\text{d}z\,\text{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 over the duration of an experiment
(2.2)$$\begin{eqnarray}Q_{m}\equiv \frac{2}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}\langle (\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0})u\rangle _{x,y,z,t},\end{eqnarray}$$which is equivalent to a buoyancy flux up to a multiplicative constant $g$. By definition $0<Q_{m}\leqslant Q$. The first inequality holds since, in our notation, negatively buoyant fluid ($\unicode[STIX]{x1D70C}_{0}<\unicode[STIX]{x1D70C}\leqslant \unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$) flows on average to the right ($u>0$) and vice versa. The second inequality would be an equality in the absence of molecular diffusion or stirring inside the duct (i.e. if all fluid moving right had density $\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$ and vice versa). In any real flow, laminar (and potentially turbulent) diffusion at the interface are responsible for an interfacial layer of intermediate density $|\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0}|<\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$ of finite thickness $\unicode[STIX]{x1D6FF}>0$ (figure 1b).
2.2 Non-dimensionalisation
A total of seven parameters are believed to play important roles in the SID: four geometrical parameters: $L$, $H$, $W$, $\unicode[STIX]{x1D703}$; and three dynamical parameters: the reduced gravity $g^{\prime }\equiv g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}$ (under the Boussinesq approximation $0<\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}\ll 1$), the kinematic viscosity of water ($\unicode[STIX]{x1D708}=1.05\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$) and the molecular diffusivity of the stratifying agent (active scalar) $\unicode[STIX]{x1D705}$. In this paper, we will primarily consider salt stratification ($\unicode[STIX]{x1D705}_{S}=1.50\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$), but will also discuss temperature stratification ($\unicode[STIX]{x1D705}_{T}=1.50\times 10^{-7}~\text{m}^{2}~\text{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: $\unicode[STIX]{x1D703}$, and the aspect ratios of the duct in the longitudinal and spanwise direction, respectively,
We choose to non-dimensionalise lengths by the length scale $H/2$, defining the non-dimensional position vector as $\tilde{\boldsymbol{x}}\equiv \boldsymbol{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 on mixing in exchange flows: $\tilde{\unicode[STIX]{x1D6FF}}\equiv \unicode[STIX]{x1D6FF}/H$, such that $\tilde{\unicode[STIX]{x1D6FF}}\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$
Consequently, we non-dimensionalise the velocity vector as $\tilde{\boldsymbol{u}}\equiv \boldsymbol{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 (or Schmidt number), as the ratio of the momentum to active scalar diffusivity
where $\unicode[STIX]{x1D705}$ takes the value $\unicode[STIX]{x1D705}_{S}$ or $\unicode[STIX]{x1D705}_{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{\unicode[STIX]{x1D70C}}\equiv (\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0})/(\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2)$, such that $\tilde{\unicode[STIX]{x1D70C}}\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{\unicode[STIX]{x1D6FF}}$ in the five-dimensional space of non-dimensional input parameters $(A,B,\unicode[STIX]{x1D703},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 (2.4) using the velocity scale $\sqrt{g^{\prime }H}$, we show here 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 $\unicode[STIX]{x0394}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
In order to show that our Reynolds number (2.4) and our non-dimensionalisation of 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 $\unicode[STIX]{x0394}U/2\sim \sqrt{g^{\prime }H}$ and the non-dimensional function $f_{\unicode[STIX]{x0394}U}(A,B,\unicode[STIX]{x1D703},Re,Pr)\sim 1$. Although some aspects of this discussion can be found in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Caulfield, Dalziel and Linden2018, Reference Lefauve, Partridge and Linden2019), the importance of this dimensional analysis for this paper justifies the more detailed discussion that we offer below.
The velocity scale $\unicode[STIX]{x0394}U$ in quasi-steady-state results from a dynamical balance in the steady, horizontal momentum equation under the Boussinesq approximation (in dimensional units)
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 $\unicode[STIX]{x2202}_{x}p=g\cos \unicode[STIX]{x1D703}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/(4L)z$, driving a flow in opposite directions on either side of the neutral level $z=0$: $-(1/\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x2202}_{x}p=g^{\prime }\cos \unicode[STIX]{x1D703}/(4L)z$;
- (G)
a gravitational body force reinforcing the flow by 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 \unicode[STIX]{x1D703}>0$ (the focus of this paper), and vice versa when the tilt angle is negative.
Rewriting (2.7) in non-dimensional form and ignoring multiplicative constants, we obtain
where $\ell$ is the smallest length scale of velocity gradients ($\ell =\unicode[STIX]{x1D6FF}$ in laminar flows, and $\ell \ll \unicode[STIX]{x1D6FF}$ 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 $\unicode[STIX]{x0394}U$ (ignoring constants and assuming $\cos \unicode[STIX]{x1D703}\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 \unicode[STIX]{x1D703}\gg g^{\prime }L\sin \unicode[STIX]{x1D703}$, i.e. when the tilt angle of the duct $\unicode[STIX]{x1D703}$ is much smaller than its ‘geometrical’ angle
(2.9)$$\begin{eqnarray}0<\unicode[STIX]{x1D703}\ll \unicode[STIX]{x1D6FC},\end{eqnarray}$$where we define the geometrical angle as
(2.10)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}\equiv \tan ^{-1}(A^{-1}).\end{eqnarray}$$Second, we can neglect the viscous (V) term if $g^{\prime }H\gg \unicode[STIX]{x1D708}\unicode[STIX]{x0394}U\ell ^{-2}L$, i.e. if the Reynolds number is larger than $Re\gg HL/\ell ^{2}$. This corresponds to
(2.11)$$\begin{eqnarray}Re\gg A\end{eqnarray}$$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 $\unicode[STIX]{x0394}U\sim \sqrt{g^{\prime }H}$, i.e. $f_{\unicode[STIX]{x0394}U}\sim 1$, which corresponds to our choice in § 2.1.
- (IG)
The inertial–gravitational balance. Using analogous arguments, if $\unicode[STIX]{x1D703}\gg \unicode[STIX]{x1D6FC}$ and $Re\gg HL/\ell ^{2}$, we expect the scaling $\unicode[STIX]{x0394}U\sim \sqrt{g^{\prime }L\sin \unicode[STIX]{x1D703}}$, i.e. $f_{\unicode[STIX]{x0394}U}(A,\unicode[STIX]{x1D703})\sim \sqrt{A\sin \unicode[STIX]{x1D703}}\gg 1$.
- (HV)
The hydrostatic–viscous balance. If $\unicode[STIX]{x1D703}\ll \unicode[STIX]{x1D6FC}$ and $Re\ll A$, we expect $f_{\unicode[STIX]{x0394}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 $\unicode[STIX]{x1D703}\gg \unicode[STIX]{x1D6FC}$ and $Re\ll A$, we expect $f_{\unicode[STIX]{x0394}U}(B,\unicode[STIX]{x1D703},Re)\sim \sin \unicode[STIX]{x1D703}Re\ll A$.
Figure 2 summarises the above analysis and the following conclusions:
(i) The parameters $A$, $\unicode[STIX]{x1D703}$ and $Re$ play particularly important roles in SID flows, since the variation of $\unicode[STIX]{x1D703}$ and $Re$ above or below thresholds set by $A$ can alter the scaling of $\unicode[STIX]{x0394}U$ (i.e. $f_{\unicode[STIX]{x0394}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<\unicode[STIX]{x1D703}\ll \unicode[STIX]{x1D6FC}$, $f_{\unicode[STIX]{x0394}U}$ increases from $\ll 1$ when $Re\ll A$ to ${\sim}1$ when $Re\gg A$. At high enough $Re$, $f_{\unicode[STIX]{x0394}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 $\unicode[STIX]{x1D703}\gg \unicode[STIX]{x1D6FC}$ and Reynolds number $Re\gg A$, $f_{\unicode[STIX]{x0394}U}$ should increase well above $1$, and likely retains a dependence on $A,B,\unicode[STIX]{x1D703},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 ever satisfy $\unicode[STIX]{x1D703}\ll \unicode[STIX]{x1D6FC}$ or $\gg \unicode[STIX]{x1D6FC}$, we consider that in general $f_{\unicode[STIX]{x0394}U}=f_{\unicode[STIX]{x0394}U}(A,B,\unicode[STIX]{x1D703},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). The existence and value of the upper edge of this region, i.e. the $Re$ value at which viscous and diffusive effects are negligible (the ‘practical $Re=\infty$ limit’) are a priori unknown.
Although the above ‘two-way’ balances do not allow us to confidently guess the scaling of $f_{\unicode[STIX]{x0394}U}$ in the blue region, theoretical arguments and empirical evidence of hydraulic control support $f_{\unicode[STIX]{x0394}U}(A,B,\unicode[STIX]{x1D703},Re,Pr)\sim 1$ for IHV, IGHV and IGV flows.
Hydraulic control of two-layer exchange flows dates back to Stommel & Farmer (Reference Stommel and Farmer1953), Wood (Reference Wood1968, Reference Wood1970) and was formalised mathematically by Armi (Reference Armi1986), Lawrence (Reference Lawrence1990) and Dalziel (Reference Dalziel1991). 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 ($\unicode[STIX]{x2202}_{x}=0$), reads
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 end of the duct), and sets the maximal non-dimensional volume flux to $\tilde{Q}=1/2$.
In ‘plug-like’ hydraulic flows ($Re\rightarrow \infty$), the velocity in each layer $\unicode[STIX]{x0394}U/2$ is equal to its layer average $Q$, giving an upper bound $f_{\unicode[STIX]{x0394}U}=\tilde{Q}=1/2$. By contrast, in real-life finite-$Re$ flows, the peak $\unicode[STIX]{x0394}U/2$ is larger than the average $Q$ (typically by a factor ${\approx}2$), such that the upper bound is $f_{\unicode[STIX]{x0394}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 $\boldsymbol{u}$ by $\sqrt{g^{\prime }H}\approx \unicode[STIX]{x0394}U/2$ in order to have $|\tilde{\boldsymbol{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 Background
We sketch the current state of knowledge on the behaviour of flow regimes, mass flux and interfacial layer thickness with input parameters in § 3.1. We highlight the limitations of previous studies and the current open questions to motivate our study in § 3.2. A more thorough review of the literature supporting these conclusions is given in appendix A, and a synthesis is given in table 2.
3.1 Current state of knowledge
The flow regimes have been observed and classified in a relatively consistent way in the literature. Throughout this paper, we adopt the nomenclature of Meyer & Linden (Reference Meyer and Linden2014): $\mathsf{L}$ (laminar flow with flat interface), $\mathsf{H}$ (interfacial Holmboe waves), $\mathsf{I}$ (intermittently turbulent), $\mathsf{T}$ (fully turbulent). The consensus is that the flow becomes increasingly disorganised and turbulent with increasing $A$, $\unicode[STIX]{x1D703}$ and $Re$. At a fixed $\unicode[STIX]{x1D703}\geqslant 0^{\circ }$, all flow regimes ($\mathsf{L},\mathsf{H},\mathsf{I},\mathsf{T}$) can be visited by increasing $Re$, and vice versa at fixed $Re$ and increasing $\unicode[STIX]{x1D703}$ (Macagno & Rouse Reference Macagno and Rouse1961; Wilkinson Reference Wilkinson1986; Kiel Reference Kiel1991; Meyer & Linden Reference Meyer and Linden2014; Lefauve et al. Reference Lefauve, Partridge and Linden2019) (hereafter MR61, W86, K91, ML14 and LPL19, respectively). Both K91 and ML14 observed regime transitions scaling with $A\tan \unicode[STIX]{x1D703}=\tan \unicode[STIX]{x1D703}/\tan \unicode[STIX]{x1D6FC}$ (or $A\unicode[STIX]{x1D703}$ for small angles), i.e. $A$ controls the $\unicode[STIX]{x1D703}$ scaling. However, the scaling in $Re$ is subject to debate, and may change on either side of $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D6FC}$ (LPL19). These conclusions are illustrated schematically in figure 3(a) (the interrogation marks denote open questions).
The mass flux has a complex non-monotonic behaviour in $A,\unicode[STIX]{x1D703},Re$ sketched in figure 3(b). While the dependence on $Re$ is clear at $Re<500A$ (MR61, W86, ML14, LPL19) due to the influence of viscous boundary layers, it is still debated at $Re>500A$: Mercer & Thompson (Reference Mercer and Thompson1975) (hereafter MT75) and ML14 argued in favour of this dependence on $Re$ even above $500A$ whereas Leach & Thompson (Reference Leach and Thompson1975) (hereafter LT75) and K91 argued against it. The mass flux reaches a maximum $Q_{m}\approx 0.4$–$0.5$ at $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D6FC}/2$ and ‘high enough’ $Re$ (MT75, K91, ML14, LPL19) and decays for smaller/larger $\unicode[STIX]{x1D703}$ and $Re$ (W86, LPL19) in a poorly studied fashion.
The interfacial layer thickness has only been studied experimentally in K91, who observed monotonic increase of $\unicode[STIX]{x1D6FF}$ with both $A$ and $\unicode[STIX]{x1D703}$, good collapse with $A\tan \unicode[STIX]{x1D703}$ (reaching its maximum $\unicode[STIX]{x1D6FF}=1$ at $\unicode[STIX]{x1D703}\gtrsim 2\unicode[STIX]{x1D6FC}$) and independence of $Re$ (figure 3c). The behaviour of $\unicode[STIX]{x1D6FF}$ at low $Re<500A$ remains unknown.
3.2 Limitations of previous studies
Many aspects of the scaling of regimes, $Q_{m}$ and $\unicode[STIX]{x1D6FF}$ with $A,B,\unicode[STIX]{x1D703},Re,Pr$ remain open questions. For example, the effects of $Re$ on $\unicode[STIX]{x1D6FF}$, and the effects of $B$ and $Pr$ on all three variables have not been studied at all. Moreover, despite our efforts to unify their findings in § 3.1 and appendix A, these past studies of the SID experiment inherently provide a fragmented view of the problem due to the following limitations (made clear by table 2):
(i) they used slightly different set-ups and geometries (e.g. presence versus absence of free surfaces in the reservoirs, rectangular ducts versus 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}$, $\unicode[STIX]{x1D6FF}$), 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 set-ups in different regions (such as $Re<500A$ and ${>}500A$).
The experimental results and models in the next two sections 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 $\unicode[STIX]{x1D6FF}$ with $A,B,\unicode[STIX]{x1D703},Re,Pr$, we obtained a comprehensive set of experimental data using an identical set-up, measuring all three dependent variables with the same methodology (described in appendix B), 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 set-up 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~\text{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{mm}$, but keeps $A,B,Pr$ identical (this is done by scaling down $H,W,L$ by the same factor $100/45$ such that the mSID and LSID ducts remain geometrically similar). Note that the mSID and LSID configurations should yield identical data at identical $Re$ since $H$ should only play a role through the non-dimensional parameters $A,B,Re$. However, we will see in §§ 4.2–4.4 that this hypothesis is challenged by our data.
- tSID
(t for tall) which differs from the HSID primarily in its tall spanwise aspect ratio $B=1/4$ (and, secondarily, in a marginally smaller height $H=90~\text{mm}$);
- mSIDT
(m for mini and T 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 $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ achieved, reflected in the lower $Re$.
Table 1 also lists, for each data set, the range of variation of $\unicode[STIX]{x1D703}$ and $Re$, and the number of data points, i.e. distinct $(\unicode[STIX]{x1D703},Re)$ couples for which we have data on regime, $Q_{m}$ and $\unicode[STIX]{x1D6FF}$.
Note that the regime and $Q_{m}$ data of the top three data sets have already been published in some form by Meyer & Linden (Reference Meyer and Linden2014) (ML14, denoted by $^{\ast }$) and Lefauve et al. (Reference Lefauve, Partridge and Linden2019) (LPL19, denoted by $^{\dagger }$), as discussed in appendices A.1–A.2. However, ML14 plotted their LSID and HSID data together (see their figures 7–8) and did not investigate their potential differences, while LPL19 only commented in passing on a fit of the $Q_{m}$ data in the $(\unicode[STIX]{x1D703},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 (Reference Lefauve2018) (especially in chapters 3 and 5, and the detailed parameters of all experiments are tabulated in his appendix A). Most of the raw and processed data used in this paper are available on the repository doi:10.17863/CAM.48821 (more details in appendix B).
Our focus on long ducts, evidenced by our choice of $A=15$ and $30$, reflects our focus on flows relevant to geophysical and environmental applications, which are typically largely horizontal ($\unicode[STIX]{x1D703}\approx 0^{\circ }$) and stably stratified in the vertical (as opposed to the different case of vertical exchange flow with $\unicode[STIX]{x1D703}=90^{\circ }$). The SID experiment conveniently exhibits all possible flow regimes, including high levels of turbulence and mixing, between $\unicode[STIX]{x1D703}=0^{\circ }$ and a few $\unicode[STIX]{x1D6FC}$ at most (§ 3.1). In long ducts (large $A$), $\unicode[STIX]{x1D6FC}\equiv \tan ^{-1}(A^{-1})$ is therefore small enough to allow us to study all the key dynamics of sustained stratified flows while keeping $\unicode[STIX]{x1D703}$ small enough for these flows to remain largely horizontal, and thus geophysically relevant.
As a result of this focus on long ducts, in the remainder of the paper we make the approximation that
This approximation is accurate to better than $2\,\%$ for the angles considered in our data sets ($\unicode[STIX]{x1D703}\leqslant 10^{\circ }$). Unless explicitly specified, $\unicode[STIX]{x1D703}$ will now be expressed in radians (typically in scaling laws).
4.2 Flow regimes
The $\mathsf{L},\mathsf{H},\mathsf{I},\mathsf{T}$ flow regimes were determined following the ML14 nomenclature as in § B.1 (except for a new regime which we discuss in the next paragraph). Figure 4 shows the resulting regime maps in the $(\unicode[STIX]{x1D703},Re)$ plane corresponding to the five data sets.
First, we note the introduction of a ‘new’ $\mathsf{W}$ regime in the tSID and mSIDT data (panels d,e). This $\mathsf{W}$ (wave) regime is similar to the $\mathsf{H}$ (Holmboe) regime, but describes interfacial waves which were not recognised as Holmboe waves in shadowgraphs. These waves were of two types. First, in the tSID geometry at positive angles $\unicode[STIX]{x1D703}>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 mSIDT 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 (Reference Lefauve2018, §§ 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 $\mathsf{H}$ and $\mathsf{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 $(\unicode[STIX]{x1D703},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 mSIDT data (panels c,e) is due to $Pr$.
However, one of the most surprising differences is that between LSID and mSID data (panels a,c), due to the dimensional height of the duct $H$ (already somewhat visible in LPL19, figure 2). It is reasonable to expect that this $H$-effect is responsible for the main differences between the LSID/HSID/tSID data and the mSID/mSIDT data. In other words, it appears that the dimensional $H$ is the main reason why the LSID/HSID/tSID transitions curves lie well above those for mSID/mSIDT, 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, such a length scale is missing in our dimensional analysis (§ 2.2) because we are unable to think of an additional length scale (such as the thickness of the duct walls or the level of the free surfaces in the reservoirs) that could play a significant dynamical role in the SID experiment.
We conclude that the transitions between flow regimes can be described by hyper-surfaces depending on all five parameters $A,B,\unicode[STIX]{x1D703},Re,Pr$ because their projections onto the $(\unicode[STIX]{x1D703},Re)$ plane for different $A,B,Pr$ do not overlap. This dependence of flow regimes on all five parameters is interesting because it was not immediately obvious from our dimensional analysis which concerned the scaling of the velocity $f_{\unicode[STIX]{x0394}U}$ alone (§ 2.3 and figure 2). Furthermore, the existence of another non-dimensional parameter involving $H$ and a ‘missing’ length scale is a major result that could not be predicted by physical intuition, and which this paper unfortunately does not elucidate.
Let us now investigate in more detail the scaling of regime transitions with respect to $\unicode[STIX]{x1D703}$ and $Re$, for which we have much higher density of data than for $A,B,Pr$. In figure 5, we replot the $\unicode[STIX]{x1D703}>0$ data of figure 4 using a log–log scale (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 $\unicode[STIX]{x1D703}Re=$ const. scaling, and dotted lines with a $\unicode[STIX]{x1D703}Re^{2}=$ const. scaling. We also show using grey shading special values of interest: $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FC}$ and $Re=50A$. The former was highlighted as particularly relevant in our scaling analysis (§ 2.3) and literature review (§ 3.1), notably as the boundary between lazy and forced flows (LPL19, § A.1). 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<50A$ (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 ) 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 § 4.3.
Figure 5 shows that regime transitions scale with $\unicode[STIX]{x1D703}Re^{2}=$ const. (dotted lines) in LSID, tSID and mSIDT (panels a,d,e), and with $\unicode[STIX]{x1D703}Re=\text{const.}$ (dashed lines) in HSID (panel b). In mSID (panel c), these two different scalings coexist: $\unicode[STIX]{x1D703}Re^{2}$ for $\unicode[STIX]{x1D703}\lesssim \unicode[STIX]{x1D6FC}$ (lazy flows) and $\unicode[STIX]{x1D703}Re$ for $\unicode[STIX]{x1D703}\gtrsim \unicode[STIX]{x1D6FC}$ (forced flows), as previously observed by LPL19, who physically substantiated the $\unicode[STIX]{x1D703}Re$ scaling in forced flows, but not the $\unicode[STIX]{x1D703}Re^{2}$ scaling in lazy flows. Furthermore, these five data sets show that this dichotomy in scalings between lazy and forced flows in mSID does not extend to all other geometries: lazy flows in the HSID exhibit a $\unicode[STIX]{x1D703}Re$ scaling and forced flows in the mSIDT exhibit a $\unicode[STIX]{x1D703}Re^{2}$ scaling. These observations further highlight the complexity of the scaling of regime transitions with $A,B,\unicode[STIX]{x1D703},Re,Pr$.
4.3 Mass flux
Mass fluxes were determined using the same salt balance methodology as ML14 described in § B.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 $\unicode[STIX]{x1D703}$ (from $\unicode[STIX]{x1D703}=-1^{\circ }$ in panel a to $\unicode[STIX]{x1D703}=3.5^{\circ }$ in panel j). The shape and colour of each symbol denote 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}$, as explained in § B.2). 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 prior to this work).
At low angles $\unicode[STIX]{x1D703}\lesssim 1^{\circ }<\unicode[STIX]{x1D6FC}$ (where $\unicode[STIX]{x1D6FC}\approx 2^{\circ }$ in mSID and $4^{\circ }$ in tSID) we observe low values $Q_{m}\approx 0.2$–0.3 in the $\mathsf{L}$ and $\mathsf{H}$ regimes. At intermediate angles $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D6FC}-2\unicode[STIX]{x1D6FC}$ 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 $\mathsf{I}$ and $\mathsf{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$). At higher angles $\unicode[STIX]{x1D703}\gtrsim \unicode[STIX]{x1D6FC}\approx 2^{\circ }$, $Q_{m}$ drops with $Re$ while remaining fairly constant with $\unicode[STIX]{x1D703}$.
As in the regime data, the mSID and tSID $Q_{m}$ data do not collapse with $Re$: all the tSID data (open symbols) are shifted to larger $Re$ compared to the mSID data (full symbols) suggesting again that a Reynolds number based on a ‘missing’ 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:
This is the general equation of a conic section, where $\unicode[STIX]{x1D71E}$ is commonly referred to as the matrix of the quadratic equation. It is well suited to describe the non-monotonic behaviour observed above, despite the fact that the non-monotonicity in $\unicode[STIX]{x1D703}$ (i.e. the decay of $Q_{m}$ at large $\unicode[STIX]{x1D703}$ widely observed in the literature) cannot be clearly confirmed by our data.
These contours describe hyperbolas ($\det \unicode[STIX]{x1D71E}<0$) for LSID, HSID and mSID (panels a–c), and concentric ellipses ($\det \unicode[STIX]{x1D71E}>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–c), or at the centre of the ellipses (panel d), and, encouragingly, no $Q_{m}=0.6$ contour exists here.
We again note that these four data sets do not collapse in the $(\unicode[STIX]{x1D703},Re)$ plane. For example, the angle at which this maximum $Q_{m}$ is achieved is a modest $\unicode[STIX]{x1D703}=0.3\unicode[STIX]{x1D6FC}$ in mSID (panel c) but appears much larger in tSID. The eigenvectors of $\unicode[STIX]{x1D71E}$ for each data set reveal that the major axis of these conic sections has equation $\unicode[STIX]{x1D703}Re^{\unicode[STIX]{x1D6FE}}$ where $\unicode[STIX]{x1D6FE}=2.6,0.3,1.5,1.2$ respectively for panels (a–d) (a larger exponent $\unicode[STIX]{x1D6FE}$ represents a larger dependence on $Re$, hence a more horizontal axis).
The exponent $\unicode[STIX]{x1D6FE}$ characterising the slope of the horizontal major axis is roughly of the same order as the exponent characterising the lines of regime transition (which is 1 for the $\unicode[STIX]{x1D703}Re$ scaling, and 2 for the $\unicode[STIX]{x1D703}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 $\unicode[STIX]{x1D6FE}=1.5$ is precisely the average of the two different regime transition exponents. This general lack of correlation suggests that the relationship between regimes and $Q_{m}$ is not straightforward and is 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 $|\unicode[STIX]{x2202}_{zz}\unicode[STIX]{x1D70C}|$ which we define as the edges of the interfacial density layer, as explained in § B.3.
In figure 7 we plot $\unicode[STIX]{x1D6FF}$ for our four duct geometries (rows) and three particular angles (representing only a subset of our data) $\unicode[STIX]{x1D703}=1^{\circ },2^{\circ },3^{\circ }$ (columns). In figure 8 we plot a quadratic fit (black contours) to all the available $\unicode[STIX]{x1D6FF}$ data (represented by the symbols) in the $(\log \unicode[STIX]{x1D703},\log Re)$ plane following (4.2). We also added in grey shading the $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FC}$ and $Re=50A$ values of interest for comparison between panels. In both figures, the shape and colour of the symbol denote the flow regimes as in figures 4–6.
In figures 7 and 8, $\unicode[STIX]{x1D6FF}$ monotonically increases with both $\unicode[STIX]{x1D703}$ and $Re$, starting from values as low as $\unicode[STIX]{x1D6FF}\approx 0.05$ in the $\mathsf{L}$, $\mathsf{H}$ and $\mathsf{W}$ regimes (see figure 13a for an illustration with $\unicode[STIX]{x1D6FF}=0.069$), and ending with values as high as $\unicode[STIX]{x1D6FF}\approx 0.8$ in the $\mathsf{T}$ regime (see figure 13c for an illustration with $\unicode[STIX]{x1D6FF}=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 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, $\unicode[STIX]{x1D6FF}_{99}\approx 10A^{1/2}(Re\,Pr)^{-1/2}\approx 0.03-0.1$ in the range $Re\in [300,6000]$ where the $\mathsf{L},\mathsf{H},\mathsf{W}$ regimes are found.
Figure 7 also shows a greater scatter of data points in the $\mathsf{I}$ and $\mathsf{T}$ regimes than in the $\mathsf{L}$ and $\mathsf{H}$ regime. This scatter cannot be attributed to measurement artefacts caused by turbulent fluctuations in the streamwise or spanwise position of the mixing layer (§ B.3), but rather demonstrates the inherent physical variability and limited reproducibility of $\mathsf{I}$ and $\mathsf{T}$ flows.
Both figures show the role of the dimensional parameter $H$ in ‘shifting’ the LSID/HSID/tSID data to higher $Re$ than the mSID data and hindering their overlap, hinting at a ‘missing’ length scale, as already discussed in the regime and $Q_{m}$ data. Note that $A$ and $B$ play additional, more subtle roles as shown by the differences between the LSID and HSID data and between the HSID and tSID data, respectively.
Finally, figure 8 shows good agreement between iso-$\unicode[STIX]{x1D6FF}$ contours and ‘iso-regime’ curves, or regime transitions curves (not shown for clarity, but easily visualised by the different symbols). This suggests that $\unicode[STIX]{x1D6FF}$ is more closely correlated to regime 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}$ with a frictional two-layer hydraulic model. In § 5.3, we tackle the scaling of $\unicode[STIX]{x1D6FF}$ in the $\mathsf{I}$ and $\mathsf{T}$ regimes by a variety of turbulence mixing models.
5.1 Volume-averaged energetics
The simultaneous volumetric measurements of the density and three-component velocity fields of LPL19 confirmed their theoretical prediction that, in forced flows, ($\unicode[STIX]{x1D703}\gtrsim \unicode[STIX]{x1D6FC}$) the time- and volume-averaged norm of the strain rate tensor (non-dimensional dissipation) followed the scaling $\langle \unicode[STIX]{x1D668}^{2}\rangle _{x,y,z,t}\sim \unicode[STIX]{x1D703}Re$ (see § A.1 for a review). They further decomposed the dissipation into:
(i) a ‘two-dimensional’ component $\unicode[STIX]{x1D668}_{2d}^{2}$ (based on the $x-$averaged velocity $\boldsymbol{u}_{2d}\equiv \langle \boldsymbol{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 $\unicode[STIX]{x1D668}_{2d}^{2}\sim \langle (\unicode[STIX]{x2202}_{z}u_{2d})^{2}\rangle _{x,y,z,t}=O(1)$. They indeed observed that $\langle \unicode[STIX]{x1D668}_{2d}^{2}\rangle _{x,y,z,t}$ plateaus at ${\approx}4$ in the $\mathsf{I}$ and $\mathsf{T}$ regimes due to the hydraulic limit;
(ii) a complementary ‘three-dimensional’ part $\unicode[STIX]{x1D668}_{3d}^{2}=\unicode[STIX]{x1D668}^{2}-\unicode[STIX]{x1D668}_{2d}^{2}$ which, as a consequence of the plateau of $\unicode[STIX]{x1D668}_{2d}^{2}$, takes over in the $\mathsf{I}$ and $\mathsf{T}$ regime and explains the $\unicode[STIX]{x1D703}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 ${\sim}10A^{1/2}Re^{-1/2}\ll 1$, and therefore volume-averaged two-dimensional dissipation to be higher $\unicode[STIX]{x1D668}_{2d}^{2}\sim \langle (\unicode[STIX]{x2202}_{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
by conjecturing that they correspond to threshold values of
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 $\unicode[STIX]{x1D703}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 ($\unicode[STIX]{x1D703}<\unicode[STIX]{x1D6FC}$, see LPL19 figure 8a) does not allow us to predict the rate of energy dissipation ($\unicode[STIX]{x1D668}^{2}$) from the rate of energy input (${\sim}\unicode[STIX]{x1D703}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 9a) 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 $\unicode[STIX]{x1D702}(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 $\unicode[STIX]{x1D702}(x)=0$ for $x\in [-A,A]$ and a volume flux $Q=u_{1}=-u_{2}=1/2$ as already mentioned in § 2.3 (see § C.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 (Reference Schijf and Schönfled1953), Assaf & Hecht (Reference Assaf and Hecht1974) and Anati, Assaf & Thompson (Reference Anati, Assaf and Thompson1977), it was formalised by Zhu & Lawrence (Reference Zhu and Lawrence2000), Gu (Reference Gu2001) and Gu & Lawrence (Reference Gu and Lawrence2005), 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 ($\unicode[STIX]{x1D703}>0$) and friction at the top and sidewalls. The full derivation of this model can be found in (L18, § 5.2) and we offer a summary in appendix C. Some of its conclusions were briefly discussed in LPL19 (their § 4.3.1), e.g. the distinction between lazy/forced flows. Below we provide a 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 $\text{d}x\times 2B\times 2$ centred around $x$ is subject to horizontal, resistive stresses at the bottom wall $\unicode[STIX]{x1D70F}_{1}^{Z}(x)$, top wall $\unicode[STIX]{x1D70F}_{2}^{Z}(x)$ (in blue), sidewalls $\unicode[STIX]{x1D70F}_{1,2}^{Y}(x)$ (respectively in the bottom and top layers, in green) and interface $\unicode[STIX]{x1D70F}^{I}$ (in red). The inclusion of these stresses in the evolution of Bernoulli potential along the duct (see § C.1) yields a nonlinear differential equation for the slope of the interface along the duct of the form
(see (C 24) for the full expression). Here, $f_{Z},f_{Y},f_{I}$ are constant friction factors parameterising, respectively, the top and bottom wall stresses, the sidewall stress and the interfacial stress (they can be determined a posteriori from any finite-$Re$ flow profile, see § C.2 and (C 21)). For any set of parameters $\unicode[STIX]{x1D703},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 $\unicode[STIX]{x1D702}(x)$ (§ C.3). The volume flux $Q$ generally increases with the forcing $\unicode[STIX]{x1D703}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 $\unicode[STIX]{x1D702}^{\prime }(x=0)$, whose simplified expression shows the balance between the forcing $\unicode[STIX]{x1D703}Re$ and the ‘composite friction parameter’ $F$
and $r_{Y}\equiv B^{-1}f_{Y}/f_{Z}$ and $r_{I}\equiv f_{I}/f_{Z}$ are respectively the sidewall friction ratio and the interfacial friction ratio.
We further note that our model has three properties: (i) the interface must slope down everywhere ($\unicode[STIX]{x1D702}^{\prime }(x)<0$) since the lower layer accelerates convectively from left to right ($u_{1}u_{1}^{\prime }(x)>0$) and vice versa ($u_{2}u_{2}^{\prime }(x)<0$); (ii) the interface must remain in the duct $|\unicode[STIX]{x1D702}(x=\pm A)|<1$; (iii) $\unicode[STIX]{x1D702}^{\prime }$ always reaches a maximum ($|\unicode[STIX]{x1D702}^{\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
and, therefore, using (5.4), we obtain the following bounds:
The first inequality in (5.5) 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 \unicode[STIX]{x1D6FC}$. The second inequality comes from (i) and (iii) and means that the mid-duct interfacial slope must be negative for $\unicode[STIX]{x1D702}(x)$ to be negative everywhere.
When suitably rescaled by $2Q\in [0,1]$, the combined friction parameter $F$ must therefore follow a $\unicode[STIX]{x1D703}Re$ scaling, strictly bounded from below. The upper bound in (5.6) is loose ($b>0$) in lazy flows, and tight ($b\rightarrow 0$) in forced flows ($A\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D703}/\unicode[STIX]{x1D6FC}\gg 1$ and $Q\rightarrow 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’ $\unicode[STIX]{x1D703}Re$ (i.e. at $\unicode[STIX]{x1D703}\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 $\unicode[STIX]{x1D703}Re$. The friction $F$ is independent of $\unicode[STIX]{x1D703}Re$ and the flow is typically laminar ($\mathsf{L}$ regime). The interface has a noticeable slope all along the duct $\unicode[STIX]{x1D702}^{\prime }(0)\ll 0$, associated with a small volume flux $Q\ll 1/2$ (see (C 25)). Such lazy flows are underspecified, and the scaling of $Q$ and $F$ with $\unicode[STIX]{x1D703}Re$ is therefore impossible to predict a priori.
(ii) At moderate $\unicode[STIX]{x1D703}Re$ ($\unicode[STIX]{x1D703}>0$): $2QF$ increases above its ‘default’ $\unicode[STIX]{x1D703}=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 ($\mathsf{L}\rightarrow \mathsf{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 $\mathsf{L}\rightarrow \mathsf{H}$ (or $\mathsf{L}\rightarrow \mathsf{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),\unicode[STIX]{x1D70C}(z)$ responsible for triggering the Holmboe instability, and the different scalings of $Q$ and $F$ that could allow the product $2QF$ to follow a $\unicode[STIX]{x1D703}Re$ scaling.
(iii) At high $\unicode[STIX]{x1D703}Re$: the hydraulic limit is reached, the upper bound is tight ($b\approx 0$), the interface is mostly flat ($\unicode[STIX]{x1D702}(x)\approx 0$ everywhere) and the inequality (5.6) becomes $2QF\approx F\approx \unicode[STIX]{x1D703}Re$. In such forced flows, the friction parameter $F$ alone must precisely balance the forcing. Arbitrarily large $\unicode[STIX]{x1D703}Re$ requires arbitrarily large $F$, which we conjecture is largely achieved by turbulent interfacial friction (increase in $r_{I}$ responsible for the $\mathsf{H}\rightarrow \mathsf{I}$ and eventually the $\mathsf{I}\rightarrow \mathsf{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 we explain in § C.2, can be written $r_{I}=1+K_{I}$, where $K_{I}$ is a turbulent momentum diffusivity (non-dimensionalised by the molecular value $\unicode[STIX]{x1D708}$) parameterising interfacial Reynolds stresses (see (C 19)). 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)
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
Comparing (5.7)–(5.8) to (5.1)–(5.2) 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 $\unicode[STIX]{x1D703}Re^{1/2}$ or $\unicode[STIX]{x1D703}^{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 $\mathsf{I}\rightarrow \mathsf{T}$ transition, undermining its usefulness to predict the $\mathsf{H}\rightarrow \mathsf{I}$ and $\mathsf{I}\rightarrow \mathsf{T}$ transitions.
From implication (ii), we understand why the volume flux $Q$, and hence the mass flux $Q_{m}$, both increase with $\unicode[STIX]{x1D703}$ and $Re$ in the $\mathsf{L}$ and $\mathsf{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 \unicode[STIX]{x1D703}Re$), which does not allow us to obtain the value of the exponent $\unicode[STIX]{x1D6FE}$ in the scaling $Q\sim \unicode[STIX]{x1D703}Re^{\unicode[STIX]{x1D6FE}}$. From implication (iii), we conjecture two potential reasons for the decrease of the mass flux $Q_{m}$ in the $\mathsf{T}$ regime (labelled ‘a’ and ‘b’ in figure 10b,c). In scenario ‘a’, $Q_{m}$ decreases due to increasing mixing despite the volume flux $Q$ staying relatively constant ($2QF\sim F\sim \unicode[STIX]{x1D703}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 $\unicode[STIX]{x1D703}Re$). Accurate $Q$ and $Q_{m}$ data obtained by volumetric measurements of velocity and density in L18 (figure 5.12b) support scenario ‘b’ up to $\unicode[STIX]{x1D703}Re=132$, but additional data are required to draw general conclusions.
The above frictional hydraulic model assumes a two-layer flow without any form of mixing, and thus ignores the behaviour of the interfacial thickness $\unicode[STIX]{x1D6FF}$, which is the subject of the next section.
5.3 Mixing models
The importance and difficulty of modelling interfacial mixing in exchange flows have long been recognised (Helfrich Reference Helfrich1995; Winters & Seim Reference Winters and Seim2000). However, despite the existence of hydraulic models for multi-layered or continuously stratified flows (Engqvist Reference Engqvist1996; Lane-Serff, Smeed & Postlethwaite Reference Lane-Serff, Smeed and Postlethwaite2000; Hogg & Killworth Reference Hogg and Killworth2004), 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}$, $\unicode[STIX]{x1D6FF}$, and their relation to fundamental stratified turbulence properties such as diapycnal diffusivity and mixing efficiency.
5.3.1 Turbulent diffusion models
Cormack, Leal & Imberger (Reference Cormack, Leal and Imberger1974a) 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 ($\unicode[STIX]{x1D6FF}=1$) in which turbulent mixing dominates to such an extent that the exchange flow is only weakly stratified in the vertical (i.e. $\langle |\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}|\rangle _{z}<1$ because $|\unicode[STIX]{x1D70C}(z=\pm 1)|<1$) and becomes stratified in the horizontal (i.e. $|\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D70C}(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. Cormack et al. (Reference Cormack, Leal and Imberger1974a) solved this problem analytically and found
Hogg et al. (Reference Hogg, Ivey and Winters2001) built on the above results and developed a model with linear velocity and density profiles within an interfacial layer of thickness $\unicode[STIX]{x1D6FF}<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. (Reference Cormack, Leal and Imberger1974a) $AK_{T}>1/15$ high-mixing limit (the ‘HM’ balance where turbulent diffusion dominates over inertia, $\unicode[STIX]{x1D6FF}=1$ and (5.9) holds) and the $AK_{T}\rightarrow 0$ hydraulic limit (the ‘IH’ balance where inertia dominates over mixing, $\unicode[STIX]{x1D6FF}=0$ and $Q=Q_{m}=1/2$ holds). Hogg et al. (Reference Hogg, Ivey and Winters2001) argued that $\unicode[STIX]{x1D6FF}$ would increase diffusively during the ‘duct transit’ advective time scale $A$, and integrated the linear velocity and density profiles across the interfacial layer to find
In order to use these models to explain the scaling of $Q_{m}$ and $\unicode[STIX]{x1D6FF}$ with $A,B,\unicode[STIX]{x1D703},Re,Pr$, we need to (i) extend them to the more complex ‘IHGM’ balance of SID flows in the $\mathsf{I}$ and $\mathsf{T}$ regimes in which gravitational forcing is present ($\unicode[STIX]{x1D703}>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$–$\unicode[STIX]{x1D716}$ model, yet these require either a length scale or the turbulent kinetic energy, which are both unknown (only the rate of dissipation is known, see (5.1)–(5.2)).
Second, borrowing from the frictional hydraulics results of § 5.2, it seems natural to conjecture that the ‘Reynolds stresses’ interfacial diffusivity $K_{I}$ in the $\mathsf{I}$ and $\mathsf{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 (5.7)–(5.8) with (5.10) would suggest
5.3.2 Previous mixing efficiency measurements and models
In this section we discuss two studies of the interfacial layer thickness $\unicode[STIX]{x1D6FF}$ 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. (Reference Prastowo, Griffiths, Hughes and Hogg2008) studied exchange flows through short ($A\approx 2{-}3$), wide ($B\gg 1$), horizontal ($\unicode[STIX]{x1D703}=0$) contractions. Their measurements suggest an approximately constant interfacial thickness $\unicode[STIX]{x1D6FF}\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, Imberger & Corcos Reference Sherman, Imberger and Corcos1978, 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 Reference Thorpe and Liu2009; Smyth & Moum Reference Smyth and Moum2013). Assuming a linear profile for $u(z)$ and $\unicode[STIX]{x1D70C}(z)$ across the mixing layer yields $Ri_{g}\approx \unicode[STIX]{x1D6FF}\approx 0.25$.
Prastowo et al. (Reference Prastowo, Griffiths, Hughes and Hogg2008) 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}$ is 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 $\unicode[STIX]{x1D6FF}$ 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 \unicode[STIX]{x1D6FF}/2\approx 0.125$.
Hughes & Linden (Reference Hughes and Linden2016) studied horizontal lock exchange gravity currents, which behave similarly to our exchange flows for part of their life cycle. They measured $\unicode[STIX]{x1D6FF}\approx 0.33$ in the range $Re\in [10^{4},10^{5}]$. Using similar measurements to Prastowo et al. (Reference Prastowo, Griffiths, Hughes and Hogg2008), 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\unicode[STIX]{x1D6FF}^{2}/3)(1-2\unicode[STIX]{x1D6FF}/3)(1-\unicode[STIX]{x1D6FF}/2)^{-2}\approx \unicode[STIX]{x1D6FF}^{2}\approx 0.08$.
However, we have seen that exchange flows in inclined ducts have $\unicode[STIX]{x1D6FF}$ monotonically increasing not only with $A$ and $Re$, but also with $\unicode[STIX]{x1D703}$. In addition, much higher values of $\unicode[STIX]{x1D6FF}\gg 0.3$ (up to 0.8, and even 1 in K91) can be achieved even at moderate values of $\unicode[STIX]{x1D703}$ of a few $\unicode[STIX]{x1D6FC}$ and $Re<10^{4}$. Therefore, the above models supporting values of $\unicode[STIX]{x1D6FF}=0.2{-}0.3$ and ${\mathcal{M}}=0.08{-}0.12$ in the $\mathsf{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 ($\unicode[STIX]{x1D703}>0$).
5.3.3 New mixing efficiency model
To address this, we propose a different 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: 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. Reference Winters, Lombard, Riley and D’Asaro1995).
As explained in LPL19 (see their §§ 4.1–4.3 and figure 8b), forced flows have, to a good approximation, the following quasi-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 $\unicode[STIX]{x1D6F7}_{P_{A}}^{adv}\approx Q_{m}\unicode[STIX]{x1D703}/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 \unicode[STIX]{x1D668}^{2}\rangle _{x,y,z,t}$. When turbulent mixing is neglected, all these fluxes have equal magnitude, and $D\approx (1/8)Q_{m}\unicode[STIX]{x1D703}$. 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 $\unicode[STIX]{x1D6F7}^{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 $|\unicode[STIX]{x1D6F7}_{P_{B}}^{adv}|=|\unicode[STIX]{x1D6F7}^{d}|$. The mixing efficiency quantifies the percentage of total time- and volume-averaged power throughput $\unicode[STIX]{x1D6F7}_{P_{A}}^{adv}$ that is spent to irreversibly mix the density field inside the duct
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)=\unicode[STIX]{x1D70C}(z)$ at either end of the duct as a minimal model to estimate the magnitude of $\unicode[STIX]{x1D6F7}_{P_{B}}^{adv}$ as a function of the interfacial layer thickness $\unicode[STIX]{x1D6FF}$, and eventually link it to input parameters $A,\unicode[STIX]{x1D703},Re$. We consider that fluid comes from the external reservoirs into the duct unmixed below (respectively above) the interfacial mixing layer at the left (respectively right) end of the duct, and leaves the duct mixed with a linear profile, going from 0 at the bottom (respectively top) edge of the mixed layer to $-1$ (respectively 1) at the top (respectively bottom) edge of the mixed layer. (In more central sections of the duct, mixing smooths 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:
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 (5.12)–(5.13) and $\unicode[STIX]{x1D6F7}_{P_{A}}^{adv}\approx Q_{m}\unicode[STIX]{x1D703}/8$, we now deduce
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 $\mathsf{I}$ and $\mathsf{T}$ regimes, we conjecture that most of the dependence on $Re$ observed in the $\unicode[STIX]{x1D6FF}$ data of figure 7 is due to the underlying monotonic increase of ${\mathcal{M}}(Re)$, which remains a priori unknown, but consistent with the observations of Prastowo et al. (Reference Prastowo, Griffiths, Hughes and Hogg2008) and Hughes & Linden (Reference Hughes and Linden2016). The observation of K91 and figure 7(a,b) that $\unicode[STIX]{x1D6FF}$ primarily scales with the group $A\unicode[STIX]{x1D703}$ at $Re\gg 500A$ (as sketched in figure 3c) would suggest that ${\mathcal{M}}$ asymptotes to a constant value at high $Re$, which is also consistent with the observations of Prastowo et al. (Reference Prastowo, Griffiths, Hughes and Hogg2008) and Hughes & Linden (Reference Hughes and Linden2016) at $ARe>10^{5}$ and $Re>10^{5}$ respectively. Assuming their high-$Re$ asymptotic value of ${\mathcal{M}}\approx 0.1$, we obtain
This gives, for example, $\unicode[STIX]{x1D6FF}\approx 0.4$ when $\unicode[STIX]{x1D703}/\unicode[STIX]{x1D6FC}\approx 1/2$. This value agrees with the K91 data (see § A.3 and figure 12f,g) and our LSID data (figure 7a, at $Re>10^{4}$ and $ARe>10^{5}$). However, this value does not agree well with our HSID, tSID and mSID $\unicode[STIX]{x1D6FF}$ data (figure 7b–d), in which $\unicode[STIX]{x1D6FF}$ 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 $\unicode[STIX]{x1D6FF}$ data and (5.15) 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 $\unicode[STIX]{x1D6FF}$ as the five non-dimensional input parameters were varied: the duct longitudinal aspect ratio $A$, spanwise aspect ratio $B$, tilt angle $\unicode[STIX]{x1D703}$, Reynolds number $Re$ and Prandtl number $Pr$.
Dimensional analysis (figure 2) and the experimental literature (figure 3, appendix A and table 2) showed that the rich dynamics of these sustained stratified shear flows were accessible for a wide range of $Re$ and for $\unicode[STIX]{x1D703}$ of at most a few duct aspect ratios $\unicode[STIX]{x1D6FC}=\tan ^{-1}(A^{-1})$. Our focus on ‘long’ ducts ($A\gg 1$) allowed us to explore these dynamics while keeping $\unicode[STIX]{x1D703}=O(\unicode[STIX]{x1D6FC})$ small enough to remain relevant to largely horizontal, stably stratified geophysical flows and turbulence, which are our ultimate motivation.
To overcome the limitations of previous studies of the problem, we presented extensive experimental results for all three variables of interests (regimes, $Q_{m}$ and $\unicode[STIX]{x1D6FF}$) in the $(\unicode[STIX]{x1D703},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 4–8) confirmed the conclusions of past studies: that increasingly disordered and turbulent regimes are found as $A,\unicode[STIX]{x1D703},Re$ are increased, that $Q_{m}$ is non-monotonic in $\unicode[STIX]{x1D703}$ and $Re$, and that $\unicode[STIX]{x1D6FF}$ is monotonic in $A,\unicode[STIX]{x1D703},Re$. Second, our data revealed the existence and importance of at least one additional non-dimensional input parameter involving the dimensional height of the duct $H$ and a ‘missing’ length scale, because our regime, $Q_{m}$ and $\unicode[STIX]{x1D6FF}$ data at the same $A,B,\unicode[STIX]{x1D703},Re,Pr$ but different $H$ do not collapse. This missing length scale has been an enduring puzzle that remains unsolved. Third, our data highlighted the complex dependence of all variables on all five parameters $A,B,\unicode[STIX]{x1D703},Re,Pr$. Regime transition, iso-$Q_{m}$, and iso-$\unicode[STIX]{x1D6FF}$ curves are not only shifted in the $(\unicode[STIX]{x1D703},Re)$ plane at different $A,B$, or $Pr$, but they also generally exhibit different power law scalings in $\unicode[STIX]{x1D703}$ and $Re$ at different $A,B,Pr$.
Given the breadth of our observations summarised above, and the relative richness of our data in the $(\unicode[STIX]{x1D703},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 $\unicode[STIX]{x1D703}Re^{\unicode[STIX]{x1D6FE}}=$ const. governing the regime transitions curves and the major axis of hyperbolas best fitting $Q_{m}$ in the $(\log \unicode[STIX]{x1D703},\log Re)$ plane. Even within this specific focus, we discovered that $\unicode[STIX]{x1D6FE}$ 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) $\unicode[STIX]{x1D6FE}$ is generally different for the regime data ($\unicode[STIX]{x1D6FE}=1$ or 2) and the $Q_{m}$ data ($0.3<\unicode[STIX]{x1D6FE}<2.6$) implying that regimes and $Q_{m}$ are not well correlated (whereas regimes and $\unicode[STIX]{x1D6FF}$ are); and (ii) in one regime data set, $\unicode[STIX]{x1D6FE}$ even takes two different values (of 1 and 2) in different regions of the $(\unicode[STIX]{x1D703},Re)$ plane.
6.3 Modelling results and outlook
To provide a modelling framework to understand the above observations (i)–(ii), we first split the $(\unicode[STIX]{x1D703},Re)$ plane into four quadrants delimited by $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FC}$ (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 boundary layers are fully developed across the duct). We then discussed three families of candidate models.
In § 5.1 we considered the volume-averaged energetics framework of Lefauve et al. (Reference Lefauve, Partridge and Linden2019) (LPL19). LPL19 physically explained the $\unicode[STIX]{x1D703}Re=$ const. scaling of regime transitions of forced ($\unicode[STIX]{x1D703}\gtrsim \unicode[STIX]{x1D6FC}$), low-$Re$ ($Re\lesssim 50A$), salt-stratified ($Pr=700$) flows as being caused by threshold values of the three-dimensional kinetic energy dissipation (5.1). We carried out the natural extension of their argument to high-$Re$ ($Re\gg 50A$) flows, by accounting for two-dimensional, laminar boundary layer dissipation. However, the resulting predicted scaling in $\unicode[STIX]{x1D703}Re-A^{-1/2}Re^{1/2}$ (5.2) 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 (Reference Gu and Lawrence2005) 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 (5.6). This model explained the qualitative behaviour of $Q_{m}$ with $\unicode[STIX]{x1D703}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 $\unicode[STIX]{x1D703}Re$ (5.7) was identical to that predicted by the energetics model and correct (at least for $Pr=700$), the high-$Re$ scaling in $A^{1/2}\unicode[STIX]{x1D703}Re^{1/2}$ (5.8) did not agree with our regime data.
Neither the energetics nor the frictional hydraulics model could predict the observed scalings in $\unicode[STIX]{x1D703}Re$ or $\unicode[STIX]{x1D703}Re^{2}$ observed in lazy flows ($\unicode[STIX]{x1D703}\lesssim \unicode[STIX]{x1D6FC}$) 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 \unicode[STIX]{x1D703},\log Re)$ plane break down for lazy flows at slightly negative angles ($-\unicode[STIX]{x1D6FC}\lesssim \unicode[STIX]{x1D703}\lesssim 0$), which we largely ignored in this paper.
In § 5.3, we focused on the scaling of $\unicode[STIX]{x1D6FF}$ underpinned by turbulent mixing. We first considered a model with constant turbulent diffusivity imposed throughout the domain (Cormack et al. Reference Cormack, Leal and Imberger1974a; Hogg et al. Reference Hogg, Ivey and Winters2001). We attempted to link this diffusivity to input parameters following insights gained from frictional hydraulic theory, but the resulting scalings (5.11) did not agree with our data. We then explained why previous measurements and models of $\unicode[STIX]{x1D6FF}$ in related stratified shear flows (Prastowo et al. Reference Prastowo, Griffiths, Hughes and Hogg2008; Hughes & Linden Reference Hughes and Linden2016) 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, (5.13)). The resulting expression for $\unicode[STIX]{x1D6FF}$ (5.14)–(5.15) is qualitatively consistent with our observations, but it involves ${\mathcal{M}}$ (not measured in our experiments) whose scaling on $Re$ is critical. Our model and data indirectly support previous observations of Prastowo et al. (Reference Prastowo, Griffiths, Hughes and Hogg2008) and Hughes & Linden (Reference Hughes and Linden2016) 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}}$ are needed 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,\unicode[STIX]{x1D703},Re$, our experimental observations have raised an even larger number of questions which remain open. Among these are the elusive existence of a sixth non-dimensional input parameter, the influence of the spanwise aspect ratio $B$ and Prandtl number $Pr$, and the scaling of the mixing efficiency ${\mathcal{M}}$.
Acknowledgements
A.L. was 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 C. Meyer and S. Vincent for collecting the regime and mass flux data in the LSID and HSID geometries in 2012–2014, and we thank L. 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 D. Kiel from Coanda Research & Development Corporation (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 and expertise of S. Dalziel, J. Partridge and the technicians of the G. K. Batchelor Laboratory. The data associated with this paper can be downloaded from the repository doi:10.17863/CAM.48821.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Literature review
In this appendix we support and complement the conclusions of § 3 by reviewing the experimental literature on the questions of flow regimes (§ A.1), mass flux (§ A.2) and interfacial layer thickness (§ A.3). We limit the discussion to the results that are most relevant to this paper, and give further details about the geometry and parameters used in each study in table 2.
A.1 Flow regimes
Macagno & Rouse (Reference Macagno and Rouse1961) (MR61) constitutes, to our knowledge, the first experimental study in a set-up 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_{\ast },R_{\ast })$, where $F_{\ast }\approx 2\sqrt{2}Q$ and $R_{\ast }\approx 4QRe$ (using our notation) are ‘effective’ Froude and Reynolds numbers. Arguing that flows above a certain $R_{\ast }$ and $F_{\ast }$ would be unstable, they proposed and verified experimentally that the ‘transition curves’ separating the flow regimes and the iso-curves of interfacial shear and mixing coefficients scaled with $R_{\ast }\,F_{\ast }\approx Re\,Q^{2}$ (i.e. these curves are $Re\,Q^{2}=~\text{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 set-up (which they attribute to Helmholtz) differs from the SID in that they were 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,\unicode[STIX]{x1D703}$ to reach target values of $R_{\ast },F_{\ast }$, without apparently realising that the flow was hydraulically controlled and that $\unicode[STIX]{x1D703}$ and $Re$ were the relevant independent input parameters.
Wilkinson (Reference Wilkinson1986) (W86) used shadowgraph and observed regime transitions similar to MR61 in a horizontal, circular pipe: ‘shear-induced instabilities $[\ldots ]$ 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 (Reference Kiel1991) (K91) used shadowgraph and laser sheet visualisations at larger $Re$, and classified the regimes differently: laminar; turbulent with $\unicode[STIX]{x1D6FF}<1$; and turbulent with $\unicode[STIX]{x1D6FF}=1$. Using a semi-empirical model based on the ratio of ‘IG’ to ‘IH’ kinetic energy scales (see § 2.3), he proposed regime transitions scaling with a ‘geometric Richardson number’ $Ri_{G}\sim A\tan \unicode[STIX]{x1D703}$ (more details in § A.2) independently of $Re$, i.e. the opposite of W86.
Meyer & Linden (Reference Meyer and Linden2014) (ML14) used shadowgraph visualisations, and (unaware of MR61) described essentially the same four regimes of MR61:
- $\mathsf{L}$
laminar flow with a thin, flat density interface;
- $\mathsf{H}$
mostly laminar flow with quasi-periodic waves on the density interface, identified as Holmboe waves;
- $\mathsf{I}$
spatio-temporally intermittent turbulence with small-scale structures and noticeable mixing between the two layers;
- $\mathsf{T}$
statistically steady turbulent flow with a thick interfacial density layer.
Interestingly, the only difference between the MR61 and ML14 nomenclatures lies in the letter characterising the wavy regime (W in MR61 and $\mathsf{H}$ in ML14), simply because MR61 observed Holmboe waves (see their figure 5) before they were explained by Holmboe (Reference Holmboe1962). ML14 mapped these regimes in the $(\unicode[STIX]{x1D703},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 $\unicode[STIX]{x1D703}>0$ (i.e. the square of the ‘IG’ velocity scaling $g^{\prime }L\sin \unicode[STIX]{x1D703}$) should be dissipated turbulently. By non-dimensionalising this excess energy by $(\unicode[STIX]{x1D708}/H)^{2}$, ML14 proposed and verified that regime transitions scale with a Grashof number (see their equation (4.4))
This scaling has two limitations: the ‘IG’ energy does not explain the transitions at $\unicode[STIX]{x1D703}=0$, and its non-dimensionalisation by $(\unicode[STIX]{x1D708}/H)^{2}$ lacks a physical basis.
Lefauve et al. (Reference Lefauve, Partridge and Linden2019) (LPL19) repeated the shadowgraph observations of ML14 in a smaller duct ($H=45~\text{mm}$ versus $H=100~\text{mm}$) with otherwise equal parameters $(A,B,Pr)=(30,1,700)$ and mapped the regimes in the $(\unicode[STIX]{x1D703},Re)$ plane. LPL19 observed two distinct scalings: a $\unicode[STIX]{x1D703}Re^{2}$ scaling for $\unicode[STIX]{x1D703}\lesssim \unicode[STIX]{x1D6FC}$ (in agreement with ML14), and a $\unicode[STIX]{x1D703}Re$ scaling for $\unicode[STIX]{x1D703}\gtrsim \unicode[STIX]{x1D6FC}$ (not observed in ML14). They developed from first principles 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 (for visualisations of flow fields in all four regimes, see their figures 2–3). They showed that for $\unicode[STIX]{x1D703}\gtrsim \unicode[STIX]{x1D6FC}$ (for so-called ‘forced flows’), the time- and volume-average rate of dissipation of kinetic energy could be predicted a priori as
where $\unicode[STIX]{x1D668}$ is the non-dimensional strain rate tensor. Because the magnitude of streamwise velocities and wall shear stresses are bounded to $O(1)$ by hydraulic control, the requirement of high strain rates at high $\unicode[STIX]{x1D703}Re$ caused transitions to increasingly three-dimensional (turbulent) flow regimes with smaller-scale gradients. The $\unicode[STIX]{x1D703}Re$ scaling of energy dissipation matched the observed regime transitions in ‘forced flow’ ($\unicode[STIX]{x1D703}\gtrsim \unicode[STIX]{x1D6FC}$), but the $\unicode[STIX]{x1D703}Re^{2}$ transition scaling in ‘lazy flows’ ($\unicode[STIX]{x1D703}\lesssim \unicode[STIX]{x1D6FC}$) remains unexplained.
A.2 Mass flux
Leach & Thompson (Reference Leach and Thompson1975) (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 CO$_{2}$/air and salt/fresh water). Surprisingly, they observed no dependence on $A,Re,Pr$.
Mercer & Thompson (Reference Mercer and Thompson1975) (MT75) reported dramatic non-monotonicity of $Q_{m}(A,\unicode[STIX]{x1D703})$: $Q_{m}\approx 0.2{-}0.3$ at $\unicode[STIX]{x1D703}=0^{\circ }$ (in agreement with LT75), increasing to $Q_{m}\approx 0.4$ at $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D6FC}/2$ and decreasing to $Q_{m}\approx 0.01-0.1$ at $\unicode[STIX]{x1D703}=90^{\circ }$ (we reproduce some of their data in figure 12b). In a small set of experiments at $\unicode[STIX]{x1D703}=30^{\circ }$ in a larger pipe ($Re=2\times 10^{4}$ versus $2\times 10^{3}$, and $A=6$), they reported 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=\unicode[STIX]{x03C0}/8\approx 0.39$ (non-dimensionalised using the pipe diameter), making the analogy with the hydraulic control arguments in Wood (Reference Wood1970) 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 \unicode[STIX]{x1D703}+A\sin \unicode[STIX]{x1D703}}$. In agreement with our dimensional analysis in § 2.3, this expression predicts a transition from an ‘IH’ balance at $0<\unicode[STIX]{x1D703}\ll \unicode[STIX]{x1D6FC}$ with $Q\approx 2/3$ to an ‘IG’ balance at $\unicode[STIX]{x1D703}\gg \unicode[STIX]{x1D6FC}$ with $Q\approx \sqrt{A\sin \unicode[STIX]{x1D703}}$. 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 \unicode[STIX]{x1D703}$’, leading to his definition of ‘geometric Richardson number’ whose inverse we interpret as being the square ratio of the ‘potential’ (‘IG’) to the ‘maximal’ (‘IH’) volume flux $Q$:
K91’s unpublished data are reproduced in figure 12. K91 obtained good collapse of his and MT75’s $Q_{m}$ data with $Ri_{G}\sim A\tan \unicode[STIX]{x1D703}$, with a peak at $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D6FC}/2$, and a decay at larger $\unicode[STIX]{x1D703}$ for a range of $A$ and $\unicode[STIX]{x1D703}$ (figure 12b–e corresponding to his figures 2.6 and 5.2). Further, in agreement with W86’s arguments, K91 reported independence of his results with $Re$ above $Re>400A$ (figure 12a) and intentionally focused on these high $Re$ throughout.
ML14 observed monotonic increase of $Q_{m}(\unicode[STIX]{x1D703})$ with $Q_{m}\approx 0.2{-}0.3$ at $\unicode[STIX]{x1D703}=0^{\circ }$ and $Q_{m}\approx 0.5$ at $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FC}/2$. They did not comment on the hint of non-monotonic behaviour suggested by their data for $\unicode[STIX]{x1D703}\gtrsim 2\unicode[STIX]{x1D6FC}$.
LPL19 observed (in passing) non-monotonic behaviour of $Q_{m}(\unicode[STIX]{x1D703},Re)$. Their data are well fitted by a hyperbolic paraboloid in the $\log \unicode[STIX]{x1D703}-\log Re$ plane, where $Q_{m}=$ const. curves are hyperbolas, with $Q_{m}\approx 0.5$ along the major axis $\unicode[STIX]{x1D703}Re^{3/2}=100$ ($\unicode[STIX]{x1D703}$ in radians), and $Q_{m}$ decays on either side of it.
A.3 Interfacial layer thickness
K91 performed conductivity probe measurements, observed monotonic increase of $\unicode[STIX]{x1D6FF}$ with both $A$ and $\unicode[STIX]{x1D703}$ (figure 12f), good collapse with $Ri_{G}\sim A\tan \unicode[STIX]{x1D703}$ (figure 12g) and independence of $Re$ (presumably because of his focus on $Re>500A$). The interfacial mixing layer is turbulent and thick ($\unicode[STIX]{x1D6FF}\approx 0.4-0.7$) at $\unicode[STIX]{x1D703}\approx \unicode[STIX]{x1D6FC}/2$ and fills the whole duct height ($\unicode[STIX]{x1D6FF}=1$) at $\unicode[STIX]{x1D703}\gtrsim 2\unicode[STIX]{x1D6FC}$. At even larger tilt angles, the mean vertical density gradient $|\unicode[STIX]{x1D70C}(z=1)-\unicode[STIX]{x1D70C}(z=-1)|/2$ drops below 1 (this ‘extreme’ turbulent scenario falls outside the scope of this paper).
Appendix B. Experimental methodology
B.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 of ML14 (see their § 3.1 and figure 3). For a schematic of the shadowgraph set-up, see 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 fluorescence, since more detailed measurements of the velocity and density fields (incompatible with simultaneous shadowgraph) have been performed in this geometry (Lefauve et al. Reference Lefauve, Partridge, Zhou, Caulfield, Dalziel and Linden2018, Reference Lefauve, Partridge and Linden2019; Partridge, Lefauve & Dalziel Reference Partridge, Lefauve and Dalziel2019).
All raw video data, including those obtained by other experimenters (acknowledged at the end of the paper), were reprocessed by the first author 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).
Most of the shadowgraph data (still images and movies) are available on the repository doi:10.17863/CAM.48821, and some of the velocity and density data are available on the repository doi:10.17863/CAM.41410 (linked to LPL19).
B.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 $\unicode[STIX]{x1D70C}_{1}^{i}=\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$ and finally at a well-mixed density $\unicode[STIX]{x1D70C}_{1}^{f}$ and ‘2’, initially at $\unicode[STIX]{x1D70C}_{0}-\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/2$ and finally at $\unicode[STIX]{x1D70C}_{2}^{f}$, giving the following two estimations:
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 misadjustment 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 (B 1) 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$ increases the relative duration of initial transients (typically fixed) and therefore the uncertainty about $T$.
Note that measurements of $Q_{m}$ in temperature-stratified experiments (mSIDT data set) could not be performed due to the practical impossibility to control the heat loss occurring through the walls of the reservoirs and the free surface.
For more details on these measurements, see L18, § 2.2.
All mass flux data (including $Q_{m,1}$ and $Q_{m,2}$, and for mSID and tSID the upper and lower bound estimations using $T=t^{d}-t^{a}$ or $t^{c}-t^{b}$) are available on the repository doi:10.17863/CAM.48821.
B.3 Interfacial layer thickness
The interfacial density layer thickness $\unicode[STIX]{x1D6FF}$ 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 grey scale 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 _{y=-B}^{B}\unicode[STIX]{x2202}_{zz}\unicode[STIX]{x1D70C}\,\text{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 methodology, illustrated in figure 13:
(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 grey scale 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 $\unicode[STIX]{x1D6FF}$.
All images were processed by the first author to ensure consistency, and yielded a total of 351 values of $\unicode[STIX]{x1D6FF}$ for all four duct geometries (table 1).
This methodology has at least two potential sources of error which we estimated to be relatively modest by performing an ad hoc set of additional measurements as we explain next.
First, the determination of $\unicode[STIX]{x1D6FF}$ from averages of shadowgraph images may give artificially large results in flows with significant streamwise variations in the vertical position of the interfacial layer. To quantify this effect, we compared measurements of $\unicode[STIX]{x1D6FF}$ made using a single snapshot to those made using an average of three snapshots in five $\mathsf{H}$ flows and five $\mathsf{T}$ flows (including those in figure 13a,c), and performing each measurement at 10 random $x$ locations (rather than one) to increase statistical robustness (total: $(5+5)\times 2\times 10=200$ measurements). We found that $\mathsf{H}$ flows were most prone to this effect (as expected from waves distorting the interface), with $\unicode[STIX]{x1D6FF}$ being estimated in three-snapshot averages an average of $14\,\%$ above its single-snapshot value, compared to a more modest average of $7.5\,\%$ in $\mathsf{T}$ flows.
Second, the determination of $\unicode[STIX]{x1D6FF}$ at a single $x$ location may not give representative results in flows with significant streamwise variations in the thickness of their interfacial layer. To quantify this effect, we investigated the streamwise variability of $\unicode[STIX]{x1D6FF}$ using the same set of 200 additional measurements by focusing on the standard deviation (spread) around the mean of each set of 10 different $x$ locations. We found that $\mathsf{H}$ flows had an average streamwise spread of $12\,\%$ of their mean, compared to a more modest $8.5\,\%$ in $\mathsf{T}$ flows.
Note that measurements of $\unicode[STIX]{x1D6FF}$ in temperature-stratified experiments (mSIDT 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)$ (sufficient to determine the flow regime but not $\unicode[STIX]{x1D6FF}$).
All interfacial thickness data (including a large number of still images, the code to determine $\unicode[STIX]{x1D6FF}$ and the quantification of errors) are available on the repository doi:10.17863/CAM.48821
Appendix C. Frictional two-layer hydraulic model
In this appendix 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 (Reference Gu2001), Gu & Lawrence (Reference Gu and Lawrence2005) but includes non-zero tilt angles and a wider range of frictional stresses suited to the SID. We cover the model formulation in § C.1, the parameterisation of frictional effects in § C.2, and the solution to the full problem in § C.3.
C.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 ($\unicode[STIX]{x2202}_{y,z}u_{1,2}=0$), while implicitly acknowledging and parameterising the effects of viscous stresses resulting from $\unicode[STIX]{x2202}_{y,z}u_{1,2}\neq 0$. This model is, however, internally 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 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)$
where $\langle \cdot \rangle _{z_{1,2}}$ denotes averaging over the lower/upper layer ($z\in [-1,\unicode[STIX]{x1D702}]$ and $z\in [\unicode[STIX]{x1D702},1]$ respectively, see figure 9a), and the velocity distribution coefficient $\unicode[STIX]{x1D706}_{1,2}$ (also called kinetic energy correction coefficient or Coriolis coefficient) is defined as
respectively in the lower and upper layer (see e.g. Chow Reference Chow1959, §§ 2.7–2.8 and Chanson Reference Chanson2004, § 3.2.2). The greater the non-uniformity of the velocity profile $u$, the larger $\unicode[STIX]{x1D706}$ is. For the SID flows considered in this paper, volumetric velocity measurements showed that $\unicode[STIX]{x1D706}$ varies over a relatively small range $1<\unicode[STIX]{x1D706}\lesssim 2$ (see L18, § 5.5.2). To simplify the following discussion, and since the effects of $\unicode[STIX]{x1D706}$ are not central here (they quantitative rather than qualitative), we make the approximation that $\unicode[STIX]{x1D706}_{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
Second, the conservation of volume and zero-net flux conditions are expressed all along the duct as
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 Reference Armi1986; Lawrence Reference Lawrence1990):
where the second equality uses (C 4) and the third equality is the control condition.
In horizontal, frictionless ducts, $E(x)=0$, hence $\unicode[STIX]{x1D702}=0$ and $u_{1}=-u_{2}=Q=1/2$ all along the duct.
When the combined effects of a small positive tilt angle $\unicode[STIX]{x1D703}>0$ and frictional stresses are added, the slope of the internal energy becomes
(this is the two-layer equivalent of single-layer ideas found in (Henderson (Reference Henderson1966), §§ 4.4–4.5)). By analogy with the topographic slope $\unicode[STIX]{x1D703}$, the ‘frictional’ slope $S(x)$ is computed by a balance of all the stresses acting on an infinitesimal slice of thickness $\text{d}x$ (figure 9b):
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), sidewall 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 side of the interface $|\unicode[STIX]{x1D70F}_{1}^{I}|=|\unicode[STIX]{x1D70F}_{2}^{I}|\equiv \unicode[STIX]{x1D70F}^{I}$. Following figure 9(b) and after elementary algebra, the balance in (C 7) can be rewritten as
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.
C.2 Parameterisation of shear stresses
We now tackle the relation between the stresses $\unicode[STIX]{x1D70F}_{i}^{j}$ and the underlying ‘real’ flow profiles $u(x,y,z)$. We start by considering the bottom wall stress of the lower layer $\unicode[STIX]{x1D70F}_{1}^{Z}$ in order to introduce the key concepts and definitions, before extending them to the other stresses. Using non-dimensional variables for $\unicode[STIX]{x1D70F}_{1}^{Z}$ and $u(x,y,z)$, we first write the dimensional equation for this stress as a simple function of the local shear
where the $\unicode[STIX]{x0394}U/2$ and $H/2$ factors come from non-dimensionalising $\unicode[STIX]{x1D70F}_{1}^{Z},u,z$, and simplify to
In order to correctly parameterise $\unicode[STIX]{x1D70F}_{1}^{Z}(x)$ and all other relevant stresses using well-defined, constant friction coefficients, we follow the following five steps.
(i) First, we define the cross-sectional ‘shape’ ($y$-$z$ dependence) of the local velocity profile in the lower layer as
(C 11)$$\begin{eqnarray}\hat{u} _{1}(x,y,z)=\frac{u(x,y,z)}{u_{1}(x)},\end{eqnarray}$$such that $\langle \hat{u} _{1}(x,y,z)\rangle _{y,z_{1}}=1$. This decomposition allows us to rewrite (C 10) as(C 12)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{1}^{Z}(x)=\frac{1}{Re}u_{1}(x)\left\langle \left|\frac{\unicode[STIX]{x2202}\hat{u} _{1}(x,y,z)}{\unicode[STIX]{x2202}z}\right|_{z=-1}\right\rangle _{y},\end{eqnarray}$$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.(ii) Second, we define a ‘layer-rescaled’ coordinate $\hat{z}_{1}$ as
(C 13)$$\begin{eqnarray}\hat{z}_{1}:=\frac{z}{1+\unicode[STIX]{x1D702}}=z\frac{Q}{u_{1}(x)},\end{eqnarray}$$in which layer 1 always has thickness one ($\hat{z}_{1}\in [-1,0]$), giving us(C 14)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{1}^{Z}(x)=\frac{1}{Re}\frac{u_{1}^{2}(x)}{Q}\left\langle \left|\frac{\unicode[STIX]{x2202}\hat{u} _{1}(x,y,\hat{z}_{1})}{\unicode[STIX]{x2202}\hat{z}}\right|_{\hat{z}_{1}=-1}\right\rangle _{y}.\end{eqnarray}$$(iii) Third, we define a constant, bottom friction parameter $f_{Z_{1}}$ to parameterise the stress
(C 15)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{1}^{Z}(x)=\frac{f_{Z_{1}}}{Re}\frac{u_{1}^{2}(x)}{Q}\quad \text{with }f_{Z_{1}}\equiv \left\langle \left|\frac{\unicode[STIX]{x2202}\hat{u} _{1}(x,y,\hat{z}_{1})}{\unicode[STIX]{x2202}\hat{z}}\right|_{\hat{z}_{1}=-1}\right\rangle _{x,y}.\end{eqnarray}$$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 $\unicode[STIX]{x2202}\hat{u} _{1}(x,y,z)/\unicode[STIX]{x2202}\hat{z}|_{\hat{z}_{1}=-1}$ must now be averaged not only over $y$ but over $x$ and $y$, as shown in (C 15). We also note that the $u_{1}^{2}(x)/Q$ factor in (C 15) 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 $\unicode[STIX]{x1D70F}^{I}$ too. However, the equivalent formulation to (C 14) for the sidewall stress in layer 1, $\unicode[STIX]{x1D70F}_{1}^{Y}$, is(C 16)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{1}^{Y}(x)=\frac{1}{Re}u_{1}(x)\left\langle \left|\frac{\unicode[STIX]{x2202}\hat{u} _{1}(x,y,z)}{\unicode[STIX]{x2202}y}\right|_{y=\pm 1}\right\rangle _{z_{1}},\end{eqnarray}$$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_{1}^{2}(x)/Q$ for $z$ derivatives).(iv) Fourth, we generalise the above definitions of $\hat{u} _{1}$ and $\hat{z}_{1}$ to both layers by defining a global $\hat{u}$ as
(C 17)and a global $\hat{z}$ as(C 18)(v) Fifth, we consider the role of turbulence at the interface, caused by Reynolds stresses which we parameterise, by analogy with (C 10), as follows
(C 19)$$\begin{eqnarray}\langle -\hat{u} ^{\prime }w^{\prime }\rangle _{x,y,z_{I},t}=\frac{1}{Re}K_{I}\left\langle \frac{\unicode[STIX]{x2202}\langle \hat{u} \rangle _{x,y,t}}{\unicode[STIX]{x2202}\hat{z}}\right\rangle _{z_{I}},\end{eqnarray}$$where $\boldsymbol{u}^{\prime }\equiv \boldsymbol{u}-\langle \boldsymbol{u}\rangle _{t}$ is the perturbation around the temporal mean and $K_{I}$ the turbulent momentum diffusivity non-dimensionalised by the molecular value $\unicode[STIX]{x1D708}$. Under these conditions, the total (molecular $+$ turbulent) interfacial stress $\unicode[STIX]{x1D70F}^{I}$ can be expressed precisely as(C 20)$$\begin{eqnarray}\unicode[STIX]{x1D70F}^{I}(x)=\frac{1+K_{I}}{Re}\frac{(u_{1}(x)-u_{2}(x))^{2}}{Q}\left\langle \left|\frac{\unicode[STIX]{x2202}\hat{u} (x,y,\hat{z})}{\unicode[STIX]{x2202}\hat{z}}\right|\right\rangle _{y,\hat{z}_{I}},\end{eqnarray}$$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
C.3 Key equations and solution method
We can now rewrite the frictional slope $S(x)$ in (C 8) using (C 21) and (C 4) as
By combining this expression for $S(x)$ with the expression for the composite Froude number $G^{2}(x)$ in (C 5) we finally obtain an expression for the differential equation governing the evolution of the interfacial slope $\unicode[STIX]{x1D702}^{\prime }(x)$ in (C 6)
where the spanwise friction ratio $r_{Y}$ and interfacial friction ratio $r_{I}$ are defined under (5.4). This equation was simplified to (5.3) for the discussion in § 5.2.
The idea behind the solution to this kind of problem can essentially be found in Gu & Lawrence (Reference Gu and Lawrence2005). However, contrary to their model (which had no tilt angle and no top and sidewall friction, i.e. $\unicode[STIX]{x1D703}=f_{Z_{1}}=f_{Y}=0$), ours does not allow us to find an analytical solution to (C 24). We therefore resort to an iterative numerical approach which we briefly outline below.
By symmetry of the problem (guaranteed under the Boussinesq approximation), $\unicode[STIX]{x1D702}$ is an odd function of $x$. We impose the boundary condition $\unicode[STIX]{x1D702}(0)=0$ and need only solve (C 24) in half of the domain (say $x\in [0,A]$).
However, since the volume flux $Q$ in (C 24) 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 $^{\ast }$)
where $\unicode[STIX]{x1D702}^{\ast }$ is the result of the forward integration of (C 24)
The coupled problem for $\unicode[STIX]{x1D702}(x)$ and $Q$ for any given set of forcing and friction parameters $(\unicode[STIX]{x1D703},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 (C 24) from $x=0$ to $A$ to get $\unicode[STIX]{x1D702}^{\ast }$ as in (C 26).
(iii) Get the $Q$ corresponding to this $\unicode[STIX]{x1D702}^{\ast }$ by the criticality condition (C 25).
(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 $(\unicode[STIX]{x1D703},Re,f_{Z},r_{Y},r_{I})$ from an experiment in the $\mathsf{L}$ regime, and quantitative agreement with $\unicode[STIX]{x1D702}(x)$ and $Q$ measurements was found L18, § 5.5.3.