Regime transitions in stratified shear flows: the link between horizontal and inclined ducts

Abstract

We present the analytical solution for the two-dimensional velocity and density fields within an approximation for laminar stratified inclined duct (SID) flows where diffusion dominates over inertia in the along-channel momentum equation but is negligible in the density transport equation. We refer to this approximation as the hydrostatic/gravitational/viscous in momentum and advective in density (HGV-A) approximation due to the leading balances in the governing equations. The analytical solution is valid for laminar flows in a two-layer configuration in the limit of long ducts. The non-dimensional volume flux within the HGV-A approximation is given by $Fr^* ={{Re}}_g/(AK)$, which is a control parameter with ${{Re}}_g$ the gravitational Reynolds number, $A$ the aspect ratio of the duct and $K$ a geometrical parameter that depends on the tilt of the duct and is obtained from the analytical solution. This analytical solution was validated against results from laboratory experiments, and allows us to gain new insight into the dynamics and properties of SID flows. Most importantly, constant values of $Fr^*$ describe, in both horizontal and inclined ducts, the transitions between increasingly turbulent flow regimes: from laminar flow, to interfacial waves, to intermittent turbulence and sustained turbulence.

1. Introduction

The understanding and improved modelling of stratified turbulent flows is of the utmost importance in several environmental situations. A particular flow configuration that attracts a lot of attention is the stratified shear flow. The interest in this flow is twofold. On the one hand, it is reminiscent of several flows occurring in the environment such as the estuarine gravitational circulation (Geyer & MacCready 2014), the lock exchange flow (Härtel, Meiburg & Necker 2000; Ottolenghi et al. 2016) and the exchange flow occurring in ocean straits (Gregg, Özsoy & Latif 1999). On the other hand, it presents rich dynamics encompassing the emergence of instabilities (Kaminski, Caulfield & Taylor 2014; Ducimetière et al. 2021), Holmboe waves (Holmboe 1962; Salehipour, Caulfield & Peltier 2016; Lefauve et al. 2018) and stratified turbulence (Salehipour, Peltier & Caulfield 2018; Lefauve & Linden 2020a; Smith, Caulfield & Taylor 2021; Lefauve & Linden 2022a,b).

The pre-eminent experimental set-up is the stratified inclined duct (SID) due to a high degree of control to explore different flow regimes and phenomena, and to the possibility of performing detailed measurements (Macagno & Rouse 1961; Meyer & Linden 2014; Lefauve, Partridge & Linden 2019a; Lefauve & Linden 2020a, 2022a,b). This set-up consists of two large tanks with fluid of different densities that are linked by an inclined, long duct (see figure 1). In recent years there has been vast progress in the understanding of the flow in SID experiments due to improved measurement capabilities that allow for simultaneous detailed measurements of the three-dimensional (3-D) density and velocity fields (Partridge, Lefauve & Dalziel 2019). A central research line has been the transitions between flow regimes: from laminar to the emergence of interfacial waves, to intermittently turbulent and to fully turbulent (Macagno & Rouse 1961; Meyer & Linden 2014; Lefauve et al. 2019a; Lefauve & Linden 2020a). Although these different regimes have been observed for 60 years, explaining them over a wide range of parameter values and determining the functional dependence of the regime transition on the governing parameters remains an unsolved problem. In fact, one of the unanswered questions is: ‘How to explain flow regime transitions in horizontal ducts or ducts inclined at a slightly negative angle?’ (Lefauve et al. 2019a).

Figure 1. Schematic representation of the side view of a SID experimental set-up. The duct of length $L$ and height $H$ is inclined at an angle $\theta$ with respect to the horizontal. The duct connects two large tanks: one with water with density $\rho =\bar {\rho }+{\rm \Delta} \rho /2$ and the other with density $\rho =\bar {\rho }-{\rm \Delta} \rho /2$. The internal angle of the duct is $\alpha =\arctan (H/L)$. The along-duct coordinate is $x$ and the coordinate perpendicular to the bottom and top of the duct is $z$. The origin $O$ of the coordinate system is located at the centre of the duct.

Lefauve et al. (2019a) distinguished between two situations: lazy and forced flows. To explain this distinction, it is necessary to define the internal angle of the duct $\alpha = \arctan (H/L)$, where $H$ is the height of the duct and $L$ its length. Lazy flows are defined as those occurring when the inclination angle $\theta$ of the duct satisfies $\alpha \gg \theta > -\alpha$, and forced flows as those occurring when $\theta >\alpha$. Between lazy and forced flows, a smooth transition occurs. The term forced refers to the increased importance of the gravitational forcing due to the duct tilt. Meyer & Linden (2014) and Lefauve et al. (2019a) have proposed two different criteria for the regime transitions in forced flows showing good agreement with experimental data. However, the criterion proposed by Meyer & Linden (2014) is not valid for $\theta \leq 0$, and the one proposed by Lefauve et al. (2019a) is not valid for lazy flows.

The current paper proposes a new criterion for the regime transitions in SID experiments that spans both lazy and forced flows (i.e. encompassing slightly negative inclinations, horizontal ducts and positive inclinations). This criterion is based on a formal perturbation analysis for long ducts and the analytical solution of the resulting simplified set of equations. This work builds upon the recent description by Kaptein et al. (2020) of what they called the high-advection/low-diffusion approximation for laminar flows in horizontal ducts. In this approximation, viscous diffusion dominates over inertia in the along-channel momentum equation while diffusion is negligible in the density transport equation. In fact, the flow gets organized into two layers with a sharp interface in between. For consistency with the work by Lefauve & Linden (2020a), we will refer here to this approximation as the HGV-A approximation because of the hydrostatic/gravitational/viscous (denoted by HGV) balance in the momentum equation and the dominance of advection (denoted by A) in the density transport equation (a detailed derivation is presented in § 3.1). This approximation is, in general, only possible if the Schmidt number (the ratio between the kinematic viscosity of the fluid and the diffusivity of the scalar responsible for the density differences) is much larger than unity, for example, for salt-stratified flows. Kaptein (2021) proposed new curves based on this approximation describing the regime transitions with a generalized Reynolds number. These curves show an overall better agreement with experimental results for salt-stratified flows than previously proposed criteria. However, they are based on an assumption of the slope of the interface as a function of the ratio $\theta /\alpha$, which has not been confirmed.

To provide a physical explanation for the transitions and propose the curves determining them, we derive, in the current paper, a two-dimensional (2-D) analytical solution for the velocity and density fields in the HGV-A approximation in both horizontal and slightly inclined ducts. This analytical solution provides several new insights into SID flows. In particular, it provides the non-dimensional parameter governing the regime transition for SID flows for horizontal and inclined ducts. The new understanding of SID flows should allow better targeted experimental campaigns to answer remaining questions due to the clearer view of the parameter space. Furthermore, the tilt angle in relatively small-scale, well-controlled SID experiments can be used to achieve more turbulent flows than in horizontal ducts with the same values of the other governing parameters. This means that understanding the link between horizontal and inclined ducts can be fruitful to connect the results of SID turbulence (which is primarily found in slightly inclined ducts) to larger-scale environmental flows (which are primarily horizontal).

2. Description of the system and background

The SID set-up, mentioned earlier and sketched in figure 1, consists of two tanks with fluid at densities $\bar {\rho }\pm {\rm \Delta} \rho /2$ (due to differences in, for example, salt concentration or temperature), joined by a duct. The duct has length $L$ and height $H$, and it is inclined at an angle $\theta$ with respect to the horizontal. The fluid is considered to have uniform and constant viscosity $\nu$. It is convenient to define the buoyancy velocity scale $U_g \equiv \sqrt {g' H}$, where $g'\equiv g {\rm \Delta} \rho /\bar {\rho }$ with $g$ the gravitational acceleration. Besides the inclination angle of the duct $\theta$, the system can be described by three non-dimensional parameters: the aspect ratio of the duct $A \equiv \cot \alpha = L/H$, the gravitational Reynolds number

(2.1)\begin{equation} {{Re}}_g \equiv\frac{H U_g}{2\nu}= \frac{H \sqrt{g' H}}{2\nu} \end{equation}

and the Schmidt number $Sc \equiv \nu /\kappa$, with $\kappa$ the diffusivity of salt (or heat, in which case the Schmidt number is referred to as the Prandtl number). We consider long ducts ($A \gg 1$) for which $A^{-1}\approx \alpha$. For ducts with finite width $W$, we must introduce an additional parameter $B\equiv W/H$.

Meyer & Linden (2014) proposed an empirical criterion for the transition between different flow regimes by defining the Grashof number as ${{Gr}} \equiv 2A {{Re}}_g^2\sin \theta$, which quantifies the ratio of the buoyancy force to the viscous force. They proposed the critical value ${{Gr}}= 4\times 10^7$ for the transition between the intermittently turbulent and fully turbulent regimes showing good agreement with experimental results. Lefauve et al. (2019a) proposed that the transitions between different regimes for a SID set-up with a given $A$ value occur at constant $\theta {{Re}}_g$ values. Lefauve & Linden (2020a) checked the proposed transitions at constant $\theta {{Re}}_g$ values against several experimental data sets including those of Meyer & Linden (2014). They remarked particularly good agreement with experiments for forced flows ($\theta >\alpha$) when, in addition, $A^{-1} {{Re}}_g \lesssim 50$. However, the comparison was inconclusive for other values of $\theta$ and $A^{-1} {{Re}}_g$.

Both previously mentioned criteria (those proposed by Meyer & Linden 2014; Lefauve et al. 2019a) have a crucial shortcoming: they are not valid for $\theta \leq 0$ because the proposed governing parameters (${{Gr}}$ and $\theta {{Re}}_g$) are then equal or smaller than zero. This would mean that, for $\theta \leq 0$, the flow does not transition away from laminar, which is in disagreement with experimental results (see e.g. figure 4 by Lefauve & Linden 2020a). Moreover, it is well known that the Grashof number defined as ${{Gr}} = {{Re}}_g^2$ is the governing parameter in the case of a horizontal duct (see e.g. Härtel et al. 2000; Hogg, Ivey & Winters 2001). More precisely, the governing parameter is ${{Gr}} A^{-2} = {{Re}}_g^2 A^{-2}$ according to Hogg et al. (2001). Hence, the definition of ${{Gr}}$ proposed by Meyer & Linden (2014) is inconsistent with the relevant definition for horizontal ducts. For these reasons, there is still a need to find a physical explanation and a generalized governing parameter determining the transitions in positively inclined, horizontal and negatively inclined ducts.

For horizontal ducts, it is known that diffusion dominates and $U/U_g\propto {{Re}}_g A^{-1}$ (with $U$ denoting the typical magnitude of the velocity) for ${{Re}}_g A^{-1} \ll (180/ {{Sc}})^{1/2}$ (Kaptein et al. 2020). This approximation is known as the viscous advective–diffusive (VAD) solution (Cormack, Leal & Imberger 1974; Hogg et al. 2001), the hydrostatic–viscous balance (Lefauve & Linden 2020a) or the diffusion-dominated regime (Kaptein et al. 2020). In the opposite case of ${{Re}}_g A^{-1}\gg (180/ {{Sc}})^{1/2}$, Kaptein et al. (2020) found two distinct approximations that arise depending on ${{Sc}}$: one for $Sc\approx 1$ and another for $Sc\gg 1$.

For ${{Sc}}\approx 1$, the flow tends to the hydraulic limit, in which $U/U_g\propto 1$ (Hogg et al. 2001). This is the theoretical limit for a steady, inviscid, irrotational, hydrostatic flow in which the two layers are of equal thickness all along the duct (Gu & Lawrence 2005; Lefauve & Linden 2020a). For ${{Sc}}\gg 1$, Kaptein et al. (2020) found that $U/U_g\propto {{Re}}_g A^{-1}$, similarly to the VAD solution but with a different proportionality constant. This behaviour would correspond to what we refer here to as the HGV-A approximation. However, the observed scaling cannot continue indefinitely as the value of ${{Re}}_g$ increases, since the hydraulic limit also exists for $Sc\gg 1$, and hence, there is the upper bound $U/U_g=1$. Most probably, Kaptein et al. (2020) did not observe the transition towards the hydraulic limit because they only modelled laminar, steady flows. In fact, the limit of the parameter space that they explored was established by the emergence of waves and instabilities at the interface. This suggests that the HGV-A approximation should hold for flows with ${{Sc}}\gg 1$ between the VAD solution (which is not hydraulically controlled) and a hydraulically controlled flow. It is this hypothesis that we will endeavour to verify in this paper.

3. Analytical description of the HGV-A approximation

In this section we focus on the HGV-A approximation and derive its consequences for regime transitions in inclined ducts. Although scaling analysis for a SID flow has been previously done by Lefauve & Linden (2020a), the HGV-A approximation was not considered. Furthermore, we take a slightly different approach in the non-dimensionalization that allows us to simplify the equations in a mathematically formal way by using an asymptotic analysis of the momentum and density transport equations describing the flow in a long duct ($A\gg 1$) (see Van Dyke 1975, for background theory). Such an analysis has been proven to be a powerful tool to analyse problems where the geometrical shape of the domain introduces a small parameter (see, e.g. Rienstra & Chandra 2001; Duran-Matute et al. 2012), which in this case is $A^{-1}$. In fact, Cormack et al. (1974) used a similar approach to derive the VAD approximation in a closed container.

First, in § 3.1 we present the non-dimensional governing equations that are later used, in § 3.2, to derive the HGV-A approximation using asymptotic analysis for a long duct. In § 3.3 we present the analytical solution for the 2-D velocity and density fields. Then, in § 3.4 we discuss its implications for the regime transitions.

3.1. Governing equations

We consider the 2-D, steady flow in a $(x,z)$ cross-section of an infinitely wide duct, where $x$ is the along-channel coordinate and $z$ is the coordinate going from the bottom to the top of the duct. The fluid velocity is $\boldsymbol {v}=(u,0,w)$. The flow is described by the continuity and the steady Navier–Stokes equations with the Boussinesq approximation for an incompressible fluid,

(3.1)$$\begin{gather} \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z} = 0, \end{gather}$$

(3.2)$$\begin{gather}u\frac{\partial u}{\partial x}+w\frac{\partial u}{\partial z} ={-}\frac{1}{\bar{\rho}}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial z^{2}}\right)+g \frac{\rho'}{\bar{\rho}}\sin \theta, \end{gather}$$

(3.3)$$\begin{gather}u\frac{\partial w}{\partial x}+w\frac{\partial w}{\partial z} ={-}\frac{1}{\bar{\rho}}\frac{\partial p}{\partial z}+\nu \left(\frac{\partial^{2} w}{\partial x^{2}}+\frac{\partial^{2} w}{\partial z^{2}} \right) - g \frac{\rho'}{\bar{\rho}}\cos \theta, \end{gather}$$

where $p$ is the pressure and $\rho '=\rho -\bar {\rho }$ is the variable part of the total density, with $\bar {\rho } \gg \rho '$. In the reference frame of the inclined duct, the gravity vector is $\boldsymbol {g}=(g \sin \theta, 0, -g \cos \theta )$. A linear equation of state relates the density $\rho$ to, for example, salt concentration or fluid temperature. In this way, the density $\rho$ is also governed by a steady transport equation

(3.4)\begin{equation} u\frac{\partial \rho}{\partial x}+ w\frac{\partial \rho}{\partial z}= \kappa \left( \frac{\partial^{2} \rho}{\partial x^{2}} + \frac{\partial^{2} \rho}{\partial z^{2}} \right). \end{equation}

We define the non-dimensional variables denoted by a tilde such that

(3.5a–f)\begin{equation} u=U \tilde{u},\quad w=\frac{U}{A} \tilde{w},\quad x=\frac{L}{2}\tilde{x},\quad z=\frac{H}{2}\tilde{z},\quad \rho'={\rm \Delta}\rho \tilde{\rho}',\quad p=\frac{2\bar{\rho} U \nu L}{H^2} \tilde{p}. \end{equation}

Using the non-dimensional variables defined in (3.5a–f), (3.1)–(3.4) can be written as

(3.6)$$\begin{gather} \frac{\partial \tilde{u}}{\partial \tilde{x}}+\frac{\partial \tilde{w}}{\partial \tilde{z}}=0, \end{gather}$$

(3.7)$$\begin{gather}A^{{-}1} {{Re}}_g\, {{Fr}} \left(\tilde{u}\frac{\partial \tilde{u}}{\partial \tilde{x}}+\tilde{w} \frac{\partial \tilde{u}}{\partial \tilde{z}}\right) ={-}\frac{\partial \tilde{p}}{\partial \tilde{x}}+A^{{-}2} \frac{\partial^{2} \tilde{u}}{\partial \tilde{x}^2}+ \frac{\partial^{2} \tilde{u}}{\partial \tilde{z}^2}+ \frac{{{Re}}_g}{2 {{Fr}}} \tilde{\rho}'\sin \theta, \end{gather}$$

(3.8)$$\begin{gather}A^{{-}3} {{Re}}_g {{Fr}} \left(\tilde{u}\frac{\partial \tilde{w}}{\partial \tilde{x}}+\tilde{w} \frac{\partial \tilde{w}}{\partial \tilde{z}}\right) ={-}\frac{\partial \tilde{p}}{\partial \tilde{z}}+A^{{-}4} \frac{\partial^{2} \tilde{w}}{\partial \tilde{x}^{2}}+A^{{-}2} \frac{\partial^{2} \tilde{w}}{\partial \tilde{z}^{2}} - A^{{-}1}\frac{{{Re}}_g}{2{{Fr}}} \tilde{\rho}' \cos \theta, \end{gather}$$

(3.9)$$\begin{gather}A^{{-}1}{{Re}}_g\,{{Sc}}\,{{Fr}} \left(\tilde{u}\frac{\partial \tilde{\rho}'}{\partial \tilde{x}}+ \tilde{w}\frac{\partial \tilde{\rho}'}{\partial \tilde{z}}\right)= A^{{-}2}\frac{\partial^{2} \tilde{\rho}'}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{\rho}'}{\partial \tilde{z}^{2}}, \end{gather}$$

where the Froude number is defined as

(3.10)\begin{equation} {{Fr}}\equiv \frac{U}{U_g} = \frac{1}{HU_g}\int_{{-}H/2}^{H/2}|u|\,{\rm d} z. \end{equation}

We do not set, a priori, $U_g$ as the typical velocity scale (as done by, for example, Cormack et al. 1974) because we do not know if $U$ scales with $U_g$, and it is imperative to use a representative velocity scale to fully profit from the asymptotic analysis. In fact, determining the relationship between $U$ and $U_g$ (i.e. determining $Fr$) as a function of the control parameters of the problem is a crucial step in deriving the analytical solution in the HGV-A approximation. Hence, the Froude number represents the non-dimensional volume flow rate through the duct, and it is related to the function $f_{{\rm \Delta} U}$ discussed by Lefauve & Linden (2020a) for which they used the peak-to-peak velocity as the typical velocity scale.

3.2. Asymptotic analysis for long ducts

We now make use of the fact that $A^{-1}\ll 1$ to describe the main characteristics and balances in the HGV-A approximation using an asymptotic analysis. We introduce the symbol $O$ to denote the mathematical order of a function and $O_s$ to denote the sharp order (see Appendix A for formal definitions). Note that the symbol $O$ does not denote the physical order of magnitude because no account is kept of constants of proportionality (Van Dyke 1975). Our approach differs in three ways from that of Cormack et al. (1974). First, as already mentioned, we do not consider a priori $U_g$ as a representative velocity scale. Second, we simplify the problem by considering only the lowest order terms in the expansion. Third, we complicate the problem by not assuming that the interface is parallel to the bottom and top of the duct.

The starting point of the asymptotic analysis is to write the velocity, the pressure and the variable part of the density as asymptotic expansions on $A^{-1}$, such that

(3.11a–d)\begin{equation} \tilde{u}=\sum_{n=0}^\infty A^{{-}n}\tilde{u}_n, \quad \tilde{w}=\sum_{n=0}^\infty A^{{-}n}\tilde{w}_n, \quad \tilde{p}=\sum_{n=0}^\infty A^{{-}n}\tilde{p}_n, \quad \tilde{\rho}'=\sum_{n=0}^\infty A^{{-}n}\tilde{\rho}'_n. \end{equation}

As a second step, it is necessary to express the response parameters, in this case, $Fr$ as a function of $A^{-1}$. We write

(3.12)\begin{equation} {{Fr}} = C A^{{-}l}, \end{equation}

where $C$ is a finite constant such that ${{Fr}} =O_s(A^{-l})$ with $l \in \mathbb {N}$. We now substitute (3.11a–d) and (3.12) into (3.7)–(3.9), and extract the equations for the lowest (zeroth-order) terms.

We start the analysis with the $x$ and $z$ components of the Navier–Stokes equations, (3.7) and (3.8), respectively. Considering $A^{-1}\ll 1$ immediately yields, as expected for long ducts, that horizontal viscous diffusion (terms with $\partial ^2/\partial x^2$) can be neglected from these equations when compared with vertical viscous diffusion (terms with $\partial ^2/\partial z^2$) because of the additional $A^{-2}$ factor. The lowest order term of vertical viscous diffusion in (3.8) is $O(A^{-2})$ and can also be immediately neglected with respect to the lowest order terms of the pressure gradient and gravity forces that are $O(A^{0})$.

Different balances can be obtained by considering different values of $l$ in (3.12). However, here we are interested only in the laminar flow for which we assume that the nonlinear terms in (3.7) and (3.8) can be neglected, and we will verify this a posteriori by comparing with laboratory experiments. Neglecting the nonlinear terms formally means that we consider here only the case where ${{Re}}_g {{Fr}} A^{-1}=O(A^{-1})$. Equation (3.8) yields then, at lowest order, the hydrostatic balance

(3.13)\begin{equation} \frac{\partial \tilde{p}_0}{\partial \tilde{z}} ={-}\frac{K}{2} \tilde{\rho}'_0\cos \theta, \end{equation}

with

(3.14)\begin{equation} K \equiv \dfrac{{{Re}}_g}{A\,{{Fr}}}. \end{equation}

However, the hydrostatic balance (3.13) can only hold if $K\cos \theta$ is a finite constant (formally, if $K\cos \theta =O_s(A^0)$), so that the terms on both sides are of the same order. Furthermore, $\theta \ll 1$ so $\cos \theta \approx 1$, and thus, $K$ must be a finite constant ($0< K<\infty$) in the limit $A^{-1}\to 0$ (formally, $K=O_s(A^0$)).

We now turn our attention to (3.7), where the remaining terms yield

(3.15)\begin{equation} \frac{\partial^{2} \tilde{u}_0}{\partial \tilde{z}^{2}}= \frac{\partial \tilde{p}_0}{\partial \tilde{x}}-\frac{K}{2} \tilde{\rho}'_0 A \sin \theta. \end{equation}

Here, we can clearly see that the hydrostatic pressure and gravity are balanced by viscous momentum diffusion. This is the balance that gives rise to the HGV acronym. Here, the condition that $\theta \ll 1$ is also important so that $A \sin \theta$ remains finite when $A\gg 1$. Note that taking $\sin \theta =0$ yields the governing equation for a horizontal duct or what Lefauve & Linden (2020a) call the hydrostatic/viscous balance. Finally, the $z$ component of the velocity, $\tilde {w}$, is given by the continuity equation

(3.16)\begin{equation} \frac{\partial \tilde{u}_0}{\partial \tilde{x}}+ \frac{\partial \tilde{w}_0}{\partial \tilde{z}}=0. \end{equation}

The simplified equations (3.13), (3.15) and (3.16) were already introduced by Macagno & Rouse (1961) as the governing equations for SID flows. However, the perturbation analysis yields the additional important result that $K=O_s(A^0)$. This means that ${Fr=O_s(A^{-1})}$, i.e. that ${{Fr}}=U/U_g = {{Re}}_g /( A\,K)$ (in agreement with the numerical results of Kaptein et al. 2020), with $K^{-1}$ a proportionality constant in the limit of $A^{-1}\to 0$. We can see from (3.15) that the balance and, hence, the value of $K$ still depends on the value of $A \sin \theta$ that remains a finite parameter in the simplified set of equations. In short, the perturbation analysis so far tells us that, to find the dimensionless volumetric flow rate $Fr= {{Re}}_g /(K\, A)$, we must find the value of $K$ as a function of $A \sin \theta$ since the values of $Re_g$ and $A$ are known.

We now continue our analysis with the density transport equation (3.9), which yields at lowest order

(3.17)\begin{equation} \frac{\partial^{2} \tilde{\rho}'_0}{\partial \tilde{z}^{2}} = 0, \end{equation}

meaning that the vertical diffusion of density is equal to zero.

In short, the HGV-A approximation is characterized by a hydrostatic balance in the vertical, while in the horizontal, the flow is driven by a pressure gradient force, and by gravity in the along-duct direction if $\theta \neq 0$. Finally, we have the hallmark of the HGV-A approximation: diffusion dominates over inertia in the along-channel momentum equation but is equal to zero in the density transport equation. This combination requires that $Sc$ is large (as observed by Kaptein et al. 2020), which is the case, for example, for water flows where the density difference is caused by a difference in salt concentration (${{Sc}}\approx 700$) but not by a difference in temperature ($Sc\approx 7$).

3.3. Analytical solution in the HGV-A approximation

We consider now the governing equations at lowest order (3.13) and (3.15)–(3.17), referred to as the HGV-A approximation. For simplicity in the notation, we drop the tildes. The goal now is to find an analytical solution for this set of equations. We split the procedure in two steps. First, we determine the expressions for the along-channel velocity $u_0$ and the density $\rho '_0$. Second, we determine the value of $Fr$ that, as mentioned earlier, is equivalent to finding the value of $K$ as a function of $A \sin \theta$.

3.3.1. Determining the along-channel velocity and the density

To find an expression for the along-channel velocity $u_0$ and the density $\rho '_0$ in the HGV-A approximation, we begin with the transport equation (3.17). For the solution to this equation, we propose

(3.18)\begin{equation} \rho_0'(x,z) = \frac{1}{2} - \mathcal{H}(z-\eta(x)), \end{equation}

where $\mathcal {H}(z)$ is the Heaviside function defined as $\mathcal {H}(z>0)=1$, $\mathcal {H}(z<0)=0$ and $\mathcal {H}(z=0)=1/2$. The physical meaning of this solution is that the density gets organized in two layers with a sharp interface located at $z=\eta (x)$ with $\eta (0)=0$ due to the definition of the coordinate system.

The proposed solution satisfies $\partial ^2 \rho '_0/\partial z^2=0$ in each of the layers. However, formally, the solution is not defined at the interface, and it is convenient to approach the problem by treating each layer separately, with the interface as a boundary between them. We use the index $\zeta =\pm 1$ with $\zeta =+1$ referring to the top layer and $\zeta =-1$ to the bottom layer. Using this notation, $\rho '_{0,\zeta }=-\zeta /2$. Macagno & Rouse (1961) also solved the same governing equations for a two-layer system, but they considered the simplified problem where the interface is parallel to the top and bottom of the duct, i.e. $\eta (x)=0$.

We now turn our attention to (3.13) to derive the pressure distribution and the horizontal pressure gradient that drives the flow. Integrating this equation with respect to $z$ yields an expression for the pressure,

(3.19)\begin{equation} p_{0,\zeta}(x,z)=K\frac{\zeta}{4}z\cos \theta + \gamma_\zeta(x), \end{equation}

with $\gamma _\zeta (x)$ as integration constants that are functions of $x$ to be determined using the boundary conditions. At $z=1$, the pressure is an unknown function of $x$, but it is convenient to define, without loss of generality,

(3.20)\begin{equation} p_{0,+1} (x,z=1)={-}K[f(x)- \eta(x) \cos\theta]/4, \end{equation}

with $f(x)$ an unknown function giving the $x$ dependence of the barotropic pressure. Using (3.19) and (3.20), one obtains an expression for $\gamma _{+1}(x)$ as a function of $\eta (x)$ and $f(x)$, and applying continuity of pressure at $z=\eta (x)$ yields an expression for $\gamma _{-1}(x)$. The pressure as a function of $x$ and $z$ is then given by

(3.21)\begin{equation} p_{0,\zeta}(x,z)=\zeta\frac{K}{4}[\cos\theta(z-\eta(x)-\zeta)+\zeta f(x)] \end{equation}

so that the along-channel pressure gradient is

(3.22)\begin{equation} \frac{\partial p_{0,\zeta}}{\partial x}=\zeta\frac{K}{4}\left(\cos\theta \frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}-\zeta \frac{{\rm d} f(x)}{{\rm d}\kern0.06em x}\right), \end{equation}

where we can see that it is composed of a baroclinic part due to the sloping interface and a barotropic part given by the gradient of $f(x)$.

Substituting (3.22) into (3.15) yields

(3.23)\begin{equation} \frac{\partial^2 u_{0,\zeta}}{\partial z^2}={-}\zeta\frac{K}{4} F_\zeta(x), \end{equation}

with

(3.24)\begin{equation} F_\zeta(x)=\cos\theta \frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}-\zeta\frac{{\rm d} f(x)}{{\rm d}\kern0.06em x}-A\sin \theta. \end{equation}

Integrating twice with respect to $z$ and applying the boundary conditions $u_0(z=\pm 1) =u_0(z=\eta (x))=0$ gives

(3.25)\begin{equation} u_{0,\zeta}(x,z) ={-}\zeta\frac{K}{8}F_\zeta(x)(z-\zeta)(z-\eta(x)). \end{equation}

Since the barotropic pressure gradient is responsible for an equal volume flow rate in both directions, we use the condition of zero mean flow through the duct,

(3.26)\begin{equation} \int_{{-}1}^{1} u_0(x,z)\,{\rm d} z=0, \end{equation}

to determine ${\rm d} f(x)/{{\rm d}\kern0.06em x}$. We do this by integrating $u_{0,-1}(x,z)$ from $z=-1$ to $z=\eta (x)$ and $u_{0,+1}(x,z)$ from $z=\eta (x)$ to $z=1$, yielding

(3.27)\begin{equation} \frac{{\rm d} f(x)}{{\rm d}\kern0.06em x}={-}\frac{\eta(x)[3+\eta^2(x)]}{1+3\eta^2(x)} \left(\cos\theta \frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}-A\sin \theta \right), \end{equation}

and $F_\zeta (x)$ can be rewritten as

(3.28)\begin{equation} F_\zeta(x)= \frac{[1+\zeta\eta(x)]^3}{1+3\eta^2(x)} \left(\cos \theta \frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}-A\sin \theta\right). \end{equation}

Replacing this expression for $F_\zeta$ into (3.23) clearly shows the physical origin of the two drivers of the flow: (i) a baroclinic pressure gradient due to a sloping interface between the two layers with homogeneous density, and (ii) gravity due to the tilt of the duct. For a given value of $A$, the relative importance of these two forcing terms varies with the angle $\theta$ and the slope of the interface. For example, for horizontal ducts, the pressure gradient will be the only forcing, while it would be expected that gravity takes over with increasing $\theta$ values (particularly, if ${\rm d} \eta (x)/{{\rm d}\kern0.06em x}$ decreases simultaneously). Hence, considering these two forcing terms without neglecting a priori any of them is crucial to link our knowledge of horizontal and inclined ducts. The difficulty here is that although the value of $\theta$ is known since it is a control parameter, the value of ${\rm d} \eta (x)/{{\rm d}\kern0.06em x}$ is not.

Now, the only missing part of the solution is to determine the shape of the interface $\eta (x)$. For this, we consider the flow rate through each of the layers yielding

(3.29)\begin{equation} -\int_{\eta(x)}^{1} u_{0,+1}(x,z)\,{\rm d} z=1 =\int_{{-}1}^{\eta(x)}u_{0,-1}(x,z)\,{\rm d} z, \end{equation}

where the value of one is due to the way the velocity was made dimensionless. In this way, we obtain an autonomous differential equation for $\eta (x)$,

(3.30)\begin{equation} \cos\theta \frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}={-}\frac{48}{K} \frac{1+3\eta^2(x)}{[1-\eta^2(x)]^3}+A\sin \theta. \end{equation}

Equation (3.30) can be further used to rewrite $F_\zeta (x)$ in (3.28) yielding

(3.31)\begin{equation} F_\zeta(x)={-}\frac{48}{K}\frac{1}{[1-\zeta\eta(x)]^3}, \end{equation}

such that (3.23) is now written as

(3.32)\begin{equation} \frac{\partial ^2 u_{0,\zeta}}{\partial z^2}={-}\frac{12\zeta}{[1-\zeta\eta(x)]^3}. \end{equation}

Finally, the along-duct velocity given by (3.25) can be written as

(3.33)\begin{equation} u_{0,\zeta}(x,z) =\frac{6 \zeta}{[1-\zeta\eta(x)]^3}(z-\zeta)(z-\eta(x)), \end{equation}

where we can notice that the velocity profile at $x=0$,

(3.34)\begin{equation} u_{0,\zeta}(0,z) =6 \zeta (z-\zeta)z, \end{equation}

always has the same shape consisting of two parabolas. To reach the complete solution, it is still necessary to determine the position of the interface, which is equivalent to finding the value of $K$ as a function of $A \sin \theta$, which can be done by solving (3.30) for $\eta (x)$.

3.3.2. Determining $K$

Since the analytical solution for horizontal ducts is tractable, we focus first on this case for which (3.30) simplifies to

(3.35)\begin{equation} \frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}={-}\frac{48}{K}\frac{1+3\eta(x)^2}{[1-\eta(x)^2]^3}, \end{equation}

yielding

(3.36)\begin{equation} x={-}\frac{K}{19\,440}[320\sqrt{3} \arctan (\sqrt{3} \eta(x))-555 \eta(x)+150 \eta(x)^3 - 27 \eta(x)^5]. \end{equation}

The value of $K$ is obtained by imposing boundary conditions. Finding the slope of the interface from an autonomous equation – similar to (3.30) – was previously done for horizontal ducts by Gu & Lawrence (2005) and for inclined ducts by Lefauve & Linden (2020a) who derived the equation using internal hydraulics. Gu & Lawrence (2005) found in fact a similar expression to (3.36), but the constant to be determined from the boundary conditions was the composite Froude number given a certain magnitude of the frictional effects. When using internal hydraulics, the boundary condition is proposed using maximum exchange flow theory (Armi & Farmer 1986). This theory states that the maximum flow rate possible is such that the flow becomes critical at the ends of the channel (i.e. that the composite Froude number becomes unity at $x=\pm 1$; see e.g. Dalziel 1991; Zaremba, Lawrence & Pieters 2003; Gu & Lawrence 2005 for details). In such a case, the flow is said to be hydraulically controlled. However, enforcing the composite Froude number to be equal to unity at the edges is inconsistent with the results of the perturbation analysis leading to the HGV-A approximation. On the one side, for internal hydraulics to be valid, inertia in the $x$ component of the momentum equation should not be negligible. On the other hand, setting the value of the composite Froude number at $x=\pm 1$ would mean that the position of the interface at these locations depends on the Froude number. Imposing this boundary condition to find the value of $K$ would, in turn, make $K$ depend on the Froude number. However, this is inconsistent with the results from the asymptotic analysis that $K=O_s(A^0)$ while ${{Fr}}=O_s(A^{-1})$, which is a condition needed for the hydrostatic balance (3.13) to hold.

We propose then a different boundary condition inspired by the results of Kaptein et al. (2020) for the simulation used to exemplify the HGV-A approximation (${{Re}}_g=500$, ${A=60}$, ${{Sc}}=300$, $\theta =0^\circ$) and shown in figure 2. In this figure we can see that the interface curves sharply when reaching the end of the duct, with the light fluid turning upwards at $x=-1$ and the dense fluid turning downwards at $x=1$. The currents at the ends of the duct must have a dimensionless thickness $\epsilon$ in the $x$ direction with $\epsilon < A^{-1}\ll 1$, meaning that they are thin with respect to the length of the duct. We propose then to impose $\eta [x=\pm ( 1+\epsilon )]=\mp 1$. For long ducts, we can assume a small error of order $A^{-1}$, and enforce, instead $\eta (x=\pm 1)=\mp 1$, meaning also that ${\rm d}\eta /{{\rm d}\kern0.06em x}|_{x=\pm 1}=\mp \infty$. Note that, even though the flow inside the duct is not controlled by the composite Froude number being equal to one at the edges, the boundary conditions at the edges remain crucial in determining the slope of the interface, and hence, the volume flow rate. Taking $\eta (\pm 1)=\mp 1$ yields $K\approx 131$ and ${\rm d}\eta /{{\rm d}\kern0.06em x} |_{x=0}=-48/K \approx -0.366$ for a horizontal duct. Note that the value of ${\rm d}\eta /{{\rm d}\kern0.06em x} |_{x=0}$ is close to the slope of $-1/3$ obtained empirically and numerically by Kaptein et al. (2020) and shows good agreement with the density field shown in figure 2.

Figure 2. Density field from the numerical simulation by Kaptein et al. (2020) for ${{Re}}_g=500$, $A=60$, ${{Sc}}=300$ and $\theta =0$. The white dashed lines represent the limits of the duct ($x=\pm 1$). The black dashed line represents the interface given by (3.36) with $K=131$.

To find the value of $K$ for the inclined ducts, we follow a similar approach as for the horizontal ducts, but we do the calculations numerically. The value of $K$ is obtained by solving (3.30) while enforcing $\eta (\pm 1)=\mp 1$. We further simplify the problem by taking $\cos \theta \approx 1$ since $\theta \ll 1$. It was already discussed in § 3.2 that $K$ is a constant for $A^{-1}\to 0$ but that it depends on $A \sin \theta$, and this can be seen again in (3.30). The solution to (3.30) for three different values of $A \sin \theta$ is shown in figure 3. Several properties of the shape of the interface expected from previous work (see e.g. Gu & Lawrence 2005; Lefauve et al. 2019a; Kaptein et al. 2020) and observed in figure 2 are reproduced by the analytical solution. For example, the solution naturally yields a slope of the interface that is constant to a good approximation over a large portion of the duct around $x=0$, and that it bends up or down as it approaches the edges of the duct. Furthermore, the slope of the interface at $x=0$ is given by

(3.37)\begin{equation} S\equiv \left.\frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}\right|_{x=0}\approx\cos \theta \left.\frac{{\rm d} \eta(x)}{{\rm d}\kern0.06em x}\right|_{x=0}={-}\frac{48}{K}+A \sin\theta \leq 0, \end{equation}

and it goes towards zero for increasing values of $\theta$. The values of $K$ and $S$ as a function of $A\sin \theta$ are shown in figure 4. Note that, due to the way the variables where made dimensionless, the ‘real’ slope (e.g. as observed in the experiments) is given by $S'=S/A$.

Figure 3. Shape of the interface for three different values of $A\sin \theta$ as obtained from solving (3.30).

Figure 4. (a) Value of the parameter $K$ defined in (3.14) and (b) value $S$, the slope of the interface at $x=0$, as a function of $A \sin \theta$. The value of $K$ and $S$ are such that the solution to the autonomous equation for $\eta (x)$ (3.30) satisfies $\eta (\pm 1)=\mp 1$ assuming $\cos \theta \approx 1$. The dotted line in (b) represents the empirical approximation used by Kaptein (2021).

Previous work has suggested a change in behaviour around $A\sin \theta =1$ where the flow transitions from lazy to forced (Lefauve et al. 2019a). In particular, the slope of the interface is considered relatively flat throughout the duct ($S\approx 0$) for $A \sin \theta >1$. Although here $S$ does not tend to zero for $A \sin \theta >1$, the values of $S'$ are quite small for typical values of $A$ used in the experiments. For example, in the case of $A=30$, the height of the interface varies about 3 mm over 1 m, which could be imperceptible by the eye. Furthermore, we can see that, for $A \sin \theta >1$, the value of $K$ varies little when compared with the variation for $A \sin \theta <1$.

The fact that we obtain the value of $S$ as a function of $A \sin \theta$ is a critical difference with respect to the parametrization used by Kaptein (2021) to determine the regime transition curves. In that work, the variation of $S$ with $A \sin \theta$ was assumed based on previous results by Lefauve et al. (2019a) and Kaptein et al. (2020). In particular, it was assumed that $S= (A \sin \theta -1)/3$ for $A\sin \theta <1$ and that $S=0$ for $A\sin \theta \geq 1$. Although certain features and trends of the regime transition curves might be reproduced using these assumptions, differences are also expected when the value of $S$, as shown in figure 4(b), is considered.

3.4. Implications of the HGV-A approximation for the regime transition

The derivation of the analytical solution in the HGV-A approximation yielded that $K$ is a constant for a given value of $A\sin \theta$ by solving (3.30) while enforcing $\eta (\pm 1)=\mp 1$. For the upcoming discussion and an easier comparison with the work by Lefauve et al. (2019a), it is convenient to use the fact that, for a long duct with a small tilt angle ($\theta,A^{-1}\ll 1$), $A\sin \theta \approx \theta /\alpha$. In such a case, the solution in the HGV-A approximation yields that $K = {{Re}}_g/ ({{Fr}} A)$ must be a constant for a given value of $\theta /\alpha$. Both ${{Re}}_g$ and $A$ are control parameters of the problem, while ${{Fr}}$ is a response parameter equivalent to the non-dimensional volume flow rate. Hence, if the value of ${{Re}}_g A^{-1}$ is increased while keeping $\theta /\alpha$ fixed, the value of ${{Fr}}$ should also increase keeping the value of $K$ constant.

It is, here, convenient to define the Froude number in the HGV-A approximation as

(3.38)\begin{equation} {{Fr}}^* \equiv \frac{{{Re}}_g }{A K}\approx \frac{{{Re}}_g \alpha}{K}. \end{equation}

Note that ${{Fr}}^*$ can be seen as either a response parameter or a control parameter. It is a response parameter because it represents the non-dimensional volumetric flow rate that results from the choice of control parameters: ${{Re}}_g$, $\theta$ and $A$. In addition, it is a control parameter because ${{Re}}_g$ and $A$ are directly imposed for a given experiment, and $K$ is a geometrical parameter set by imposing the value of $\theta /\alpha$. The value of $K$ is determined from the analytical solution in the HGV-A approximation, and it is, to a very good approximation for $\theta \ll 1$, the value shown in figure 4(a). In this way, the value of ${{Fr}}^*$ is a known quantity for a given experiment. Within the HGV-A approximation, ${{Fr}} = {{Fr}}^*$, but this is not the case if the approximation does not hold. Furthermore, since $K\approx 131$ is a constant for $\theta = 0$, saying that ${{Fr}}^*$ is the control parameter for horizontal ducts is equivalent to saying that ${{Re}}_g\,A^{-1}$ is the control parameter as shown by Hogg et al. (2001).

To study the limit of validity of the HGV-A approximation, we consider a hypothetical SID set-up with given values for $A$ and $\theta$ satisfying $A\gg 1$ and $\theta \ll 1$. As just mentioned, the value of $K$ is known then. In this hypothetical set-up, we first take a sufficiently small value of $Re_g$ so that ${{Fr}}^*$ is sufficiently small to neglect inertia from the $x$ component of the momentum equation (3.7). Furthermore, $Sc\gg 1$ so that the flow gets organized (too good approximation) in a two-layer configuration and the HGV-A approximation holds. In such a case, $Fr={{Fr}}^*$. We now do a series of experiments increasing ${{Re}}_g$, which is equivalent to increasing ${{Fr}}^*$. We know from internal hydraulic theory that there is a maximum flow rate possible through the duct, i.e. a maximum possible value for $Fr$. In the frictionless case, this maximum value is $Fr=1$ that is known as the hydraulic limit (Hogg et al. 2001; Lefauve et al. 2019a), while friction reduces the maximum possible value of $Fr$ (Gu & Lawrence 2005). Hence, there must be a critical value ${{Fr}}^*_c$ at which the HGV-A approximation becomes invalid. For $Fr^*$ values larger than this critical value, the actual Froude number is such that $Fr< Fr^*$. Approaching this limit of validity of the HGV-A approximation, it is expected that a transition (i.e. a qualitative change in the flow) occurs.

Now, the next step is to determine the range of validity of the HGV-A approximation by determining ${{Fr}}_c^*$ since the asymptotic analysis does not give information about the physical order of magnitude at which the different terms become relevant (Van Dyke 1975). Hence, to determine the range of validity of the HGV-A approximation, it is necessary to use numerical simulations or laboratory experiments. In the following section we study this range of validity and the consequence for regime transitions in four different experimental set-ups.

4. Experimental verification

It was already discussed that several of the properties of the analytical solution in the HGV-A approximation agree with the numerical results of Kaptein et al. (2020). Now, we use results from laboratory experiments to verify the theoretical results derived in the previous section. We have two main aims: (i) to show that the solution in the HGV-A approximation described in § 3.3 does exist and is observed experimentally in a region of the parameter space, and (ii) to show that curves of ${{Fr}}^*=\text {const.}$ describe the transition between different regimes.

4.1. Description of the data

We use the experimental data sets by Meyer & Linden (2014), Lefauve et al. (2019a) and Lefauve & Linden (2020a). These data sets have been discussed and made available online (https://doi.org/10.17863/CAM.48821, https://doi.org/10.17863/CAM.41410) by Lefauve, Partridge & Linden (2019b) and Lefauve & Linden (2020a,b). Four ducts with different dimensions were used. We will refer to them as mSID (m for mini), tSID (t for tall), LSID (L for large), and HSID (H for half), in accordance with Lefauve & Linden (2020a). Table 1 summarizes the characteristics of all the set-ups. For each set-up, the inclination angle of the duct $\theta$ and the gravitational Reynolds number ${{Re}}_g$ were varied independently. The fluid used was a salt (NaCl) solution ($Sc\approx 700$), and the value of ${{Re}}_g$ was varied by changing the salt concentration in each of the tanks. In total, we used 738 data points to study the regime transitions.

Table 1. Characteristics of the experiments used in this paper. Four duct geometries (abbreviated mSID, LSID, HSID, tSID) are used (Lefauve & Linden 2020a). We list the values of the dimensionless numbers describing each duct geometry ($A$ and $B$), the value of $Sc$ for salt in water and the ranges of $\theta$ and ${{Re}}_g$ explored.

Meyer & Linden (2014) distinguished four different regimes: laminar (${\textsf {L}}$), Holmboe waves (${\textsf {H}}$), intermittently turbulent (${\textsf {I}}$), and turbulent (${\textsf {T}}$). Lefauve & Linden (2020a) introduced a ‘waves’ (${\textsf {W}}$) regime where waves other than Holmboe waves were observed. The different regimes were mostly identified by shadowgraph observations over a subsection of the duct, following the qualitative description of each regime by Meyer & Linden (2014). A schematic of the shadowgraph set-up was presented by Lefauve (2018). The observed regime and the mass flow rate as a function of the governing parameters are provided by Lefauve & Linden (2020b). A complication for the comparison with the experimental data is that, as pointed out by Lefauve & Linden (2020a), there is a surprising difference between the results from the LSID and the mSID set-ups: the regions in the $(\theta,{{Re}}_g)$ plane where the different regimes occur do not coincide even though the values of all the dimensionless parameters are identical.

To show that the proposed solution in the HGV-A approximation is observed experimentally and to observe the limit of validity, we analyse in the following section the along-duct velocity and the density fields for three experiments in the mSID duct all with $\theta =2^\circ$ (provided by Lefauve et al. 2019b). The first experiment falls within the ${\textsf {L}}$ regime with ${{Re}}_g=398$, the second within the ${\textsf {H}}$ regime with ${{Re}}_g=1059$ and the third within the ${\textsf {I}}$ regime with ${{Re}}_g=1466$. A detailed description of the experiments and the methodology is given by Lefauve et al. (2019a). With this choice of experiments, we have the same approach as mentioned for a hypothetical SID set-up in § 3.4, but with a real set-up and real experimental results.

4.2. Experimental confirmation of the HGV-A approximation

For the three experiments for which we analyse the velocity and density field, we first determine the value of $K$ numerically as explained in § 3.3.2. This value is given to good approximation in figure 4(a). For $\theta = 2^\circ$ and $A=30$, $A \sin \theta \approx \theta /\alpha \approx 1.0$ so $K\approx 39$. We then compute the value of the Froude number in the HGV-A approximation $Fr^*\equiv {{Re}}_g (KA)^{-1}$, giving ${{Fr}}^*\approx 0.34$ for the ${\textsf {L}}$ experiment, ${{Fr}}^*\approx 0.91$ for the ${\textsf {H}}$ experiment and ${{Fr}}^*\approx 1.25$ for the ${\textsf {I}}$ experiment. Here, we note already that the value of ${{Fr}}^*$ for the ${\textsf {I}}$ experiments is larger than ${{Fr}}=1$ corresponding to the hydraulic limit. In the following, we investigate if the key flow characteristics and working assumptions leading to the solution in the HGV-A approximation derived in § 3.3 hold in these experiments. We consider four aspects: (i) the density field, focusing on the slope of the interface; (ii) the $x$ component of the momentum equation (4.1), focusing on the fact that the HGV balance holds; (iii) the vertical profiles of the along-channel velocity component; and (iv) the prediction of the mass flux through the duct.

Figure 5 shows the experimentally obtained time-averaged density fields and the location of the interface as predicted by the HGV-A approximation for the three individual experiments considered. The slope of the interface close to $x=0$ is given to good approximation by (3.37), which yields $S\approx -0.20$ for these experiments. First, we see that the density is indeed organized into two layers with a sharp interface in between, with the ${\textsf {I}}$ experiment presenting a thicker interface. For the ${\textsf {L}}$ experiment, the semi-analytical solution for the slope of the interface agrees well with the experiment. This means that (3.30) with the boundary condition $\eta (\pm 1)=\mp 1$ predicts well the shape of the interface for this experiment and that the value of $K$ obtained from the analytical solution in the HGV-A approximation applies. The agreement is less for the ${\textsf {H}}$ experiment, with the slope of the interface being slightly less steep than predicted. Finally, for the ${\textsf {I}}$ experiment, the interface is fully parallel to the top and bottom of the duct, suggesting that the solution in the HGV-A approximation does not hold for this experiment.

Figure 5. Time-averaged density field in the $(x,z)$ plane at $y=0$ for three experiments in the mSID set-up. (a) Experiment within the ${\textsf {L}}$ regime ($\theta =2^\circ$, ${{Re}}_g=398$). (b) Experiment within the ${\textsf {H}}$ regime ($\theta =2^\circ$, ${{Re}}_g=1059$). (c) Experiment within the ${\textsf {I}}$ regime ($\theta =2^\circ$, ${{Re}}_g=1455$). The black dashed line represents the line $z=Sx$ with $S=-0.20$ as given by (3.37).

As a next step, we investigate if the HGV balance (3.32) is indeed the leading order balance in the $x$ component of the momentum equation (3.7) close to the centre of the duct ($x\approx 0$), which we rewrite here as

(4.1)\begin{align} \underbrace{A^{{-}1} {{Re}}_g {{Fr}} \left(u\frac{\partial u}{\partial x}+w \frac{\partial u}{\partial z}\right)}_{I} = \underbrace{A^{{-}2} \frac{\partial^{2} u}{\partial x^2}}_{II}+ \underbrace{\frac{\partial^{2} u}{\partial z^2}}_{III}+ \underbrace{\frac{{{Re}}_g }{2A{{Fr}}}\left( A \sin \theta\rho'-\cos\theta \int_{0}^{z}\frac{\partial\rho'}{\partial x}{\rm d} z'\right)}_{IV}, \end{align}

with $z'$ a dummy variable. Here, we have assumed that the hydrostatic balance (3.13) holds and followed the derivation by Kaptein et al. (2020) to write the pressure gradient term in this form. To see that the HGV balance (3.32) holds, vertical viscous diffusion (term III) must be equal to minus the forcing (term IV). In addition, we want to see that term I (inertia) is negligible with respect to term III and term IV. Figure 6 shows with black lines the vertical profiles of terms I, III and (minus) IV of (4.1) obtained from the time- and $y$-averaged velocity and density fields at $x\approx -0.2$. We do not show term II because it is practically equal to zero for all experiments. In general, we see that term I is also much smaller than the other two terms meaning that it can indeed be neglected. As can be seen, term III and (minus) term IV overlap well for all experiments, with the experiment in the ${\textsf {L}}$ regime showing better agreement, meaning that the HGV balance holds particularly well for the ${\textsf {L}}$ experiment. In the case of the ${\textsf {H}}$ and ${\textsf {I}}$ experiments, terms III and IV match in general, suggesting that the HGV balance is still the leading balance of the time-averaged flow, but some deviations are observed.

Figure 6. The vertical profiles of the time- and $y$-averaged terms of the horizontal momentum balance (4.1) as a function of $z$ and $x\approx -0.2$ for the same three experiments as in figure 5. The black dotted, solid and dashed lines denote the terms I, III and (minus) IV, respectively, as obtained from the laboratory experiments. Term II is always approximately equal to zero and not explicitly shown. The dashed grey line represents both terms III and (minus) IV as obtained from the analytic solution in the HGV-A approximation that assumes an infinitely sharp interface, and is given by the right-hand side of (3.32).

To understand these deviations, we make a comparison with the theory. Figure 6 also shows with dashed grey lines the value for terms III and IV according to the solution in the HGV-A approximation derived in § 3.3 and given by the right-hand side of (3.32), which assumes an infinitely thin interface. The interface in the experiments is not infinitely thin due to diffusion in the vertical direction that is fully neglected in the HGV-A approximation. This effect is most visible in the ${\textsf {I}}$ experiment due to turbulent diffusion across the interface. For all three experiments, the forcing corresponds well with forcing expected in the HGV-A approximation given by the sum of gravity and the baroclinic pressure gradient. It is the diffusion in the vertical direction that shows deviations for the ${\textsf {H}}$ and ${\textsf {I}}$ experiments. We note, in particular, an increase of the magnitude of vertical momentum diffusion in the bottom boundary layer, and just above and below the shear layer. Note that the reasons why the momentum balance is not closed for the ${\textsf {H}}$ and ${\textsf {I}}$ experiments is that (4.1) would need additional terms for the time-averaged experimental results. Most importantly, it misses the momentum diffusion in the across-duct direction and a term containing the time-averaged effect of temporal perturbation (similar to Reynolds stresses).

We now compare the velocity profiles. Figure 7 shows the vertical profiles of the $x$ component of the velocity at $x\approx -0.2$ and averaged in the $y$ direction for the same three experiments. Choosing different values of $x$ does not affect the results beside the fact that the height of the point with $u=0$ is slightly shifted. For the experiment in the ${\textsf {L}}$ regime, the flow is steady so the instantaneous and time-averaged velocity profiles coincide. Even for the ${\textsf {H}}$ and ${\textsf {I}}$ experiments, the variation in the instantaneous profiles is small compared with the magnitude of the time-averaged velocity. For the three experiments, we see a different level of disagreement in the amplitude of the velocity. Computing the value of ${{Fr}}$ for the experiments using the time- and $y$-averaged $x$-velocity component yields ${{Fr}}=0.40$ (compared with ${{Fr}}^*\approx 0.34$) for the ${\textsf {L}}$ experiment, ${{Fr}}=0.81$ (compared with ${{Fr}}^*\approx 0.91$) for the experiment in the ${\textsf {H}}$ regime and ${{Fr}}=0.90$ (compared with ${{Fr}}^*\approx 1.25$) for the experiment in the ${\textsf {I}}$ regime. As mentioned earlier, maximum exchange theory says that the maximum value is ${{Fr}}=1$, so we could have already expected that the Froude number for the ${\textsf {I}}$ experiment was not going to reach ${{Fr}}^*=1.25$, but we have here confirmation of this. It suggests again that the HGV-A approximation does not hold for this experiment. On the other hand, the value of the Froude number for the ${\textsf {L}}$ and ${\textsf {H}}$ experiments are well predicted (within 13 %) by the HGV-A analytical solution. We believe that this is a very good result, particularly if we consider, on the one side, experimental errors, and on the other, some of the assumptions made during the analytical derivation. For example, we assumed an infinitely wide duct, while the experiments have a cross-sectional aspect ratio $B=1$.

Figure 7. Vertical velocity profiles of the $x$ component of the velocity as a function of $z$ for $x\approx -0.2$ and averaged in the $y$ direction for the same three experiments as in figures 5 and 6. The grey lines represent the instantaneous velocity, the black dashed line the time-averaged velocity and the red line represents the velocity in the HGV-A approximation given by (3.33).

In the previous paragraphs, we studied the characteristics of only three experiments with equal value of $\theta /\alpha$ and, hence, equal value of $K$. In general, we saw that the ${\textsf {L}}$ experiment is very well reproduced by the HGV-A approximation, the ${\textsf {H}}$ experiment is slightly less so, and the ${\textsf {I}}$ experiment shows clear differences. In § 3.4 we discussed that we expect that, for ${{Fr}}^*$ values larger than a critical value ${{Fr}}^*_c$, the HGV-A approximation does not hold, and the actual Froude number must be smaller than ${{Fr}}^*$. The results suggest that ${{Fr}}^*_c$ is found somewhere between the ${\textsf {H}}$ and ${\textsf {I}}$ experiment. However, it is difficult to determine clear trends from only these three experiments. Although values of $Fr$ are not available for many experiments, the mass flux

(4.2)\begin{equation} Q_m\equiv\int_{{-}1}^{1}\int_{{-}1}^{1}|\rho' u|\,{\rm d} y\,{\rm d} z \end{equation}

is available for several experiments. The value of $Q_m$ is commonly used to compare (and, thus, characterize) the flow, with $Q_m=1/2$ corresponding to the hydraulic limit (Hogg et al. 2001; Lefauve et al. 2019a), while $Q_m={{Fr}}/2={{Fr}}^*/2$ in the HGV-A approximation. The factor $1/2$ arises from the fact that $\rho '=\pm 1/2$ in the HGV-A approximation following the way we made the density dimensionless.

The value of $Q_m$ as a function ${{Fr}}^*/2$ for the experiments in the mSID set-up are shown in figure 8. The values of $Q_m$ can then be compared with the expected values in the HGV-A approximation ($Q_m={{Fr}}^*/2$) and to the hydraulic limit ($Q_m=1/2$). For ${{Fr}}^*/2\lesssim 1/2$, $Q_m\approx {{Fr}}^*/2$, showing that the HGV-A approximation gives a good estimate for the value of $Q_m$. Although the values for some specific experiments are overestimated, and some others are underestimated, clear trends are visible. For example, for the experiments with $\theta \geq 1.0^\circ$, $Q_m\approx {{Fr}}^*/2$ up to ${{Fr}}^*/2=1/2$. Then, a sharp transition to the hydraulic limit occurs, and $Q_m\approx 1/2$ for ${{Fr}}^*/2\geq 1/2$. For experiments with $\theta = 0.5^\circ$, $Q_m$ seems to saturate reaching a maximum value $Q_m\approx 0.4<1/2$. The results for $\theta = 0^\circ$ suggest that $Q_m$ will similarly not reach the value $1/2$ for a horizontal duct. The impossibility of reaching $Q_m=1/2$ is in agreement with frictional hydraulic theory that states that frictional effects reduce the maximum possible value of $Q_m$ (Gu & Lawrence 2005). It is only for $\theta \geq 1.0^\circ$ that the additional forcing by gravity in the along-duct direction reduces the importance of frictional effects and makes the direct transition from the HGV-A approximation to the hydraulic limit ($Q_m=1/2$) possible.

Figure 8. Mass flux per unit width $Q_m$ as a function of the parameter ${{Fr}}^*/2$ for the experiments in the mSID set-up. The colour denotes the value of the angle $\theta$. The solid black lines represent $Q_m={{Fr}}^*/2$ (the expected value in the HGV-A approximation) and $Q_m=1/2$ (the expected value in the hydraulic limit). The experimental values tend to follow the trend of the prediction based on the HGV-A approximation $Q_m\approx Fr^*/2$ until ${{Fr}}^*/2\approx 0.5$. For larger values of ${{Fr}}^*/2$, a maximum, constant value of $Q_m$ is reached for each value of $\theta$ as expected from (frictional) two-layer hydraulics. The overall maximum value is $Q_m\approx 0.5$ as predicted for the hydraulic limit.

Figure 8 is reminiscent of similar plots by Hogg et al. (2001) where they plotted $Q_m$ as a function of $({{Re}}_g A^{-1})^2$ for simulations with $Sc=1$ in a horizontal duct. In that case, they observed the change from the VAD solution to the hydraulic limit. In our case with $Sc\approx 700$, we observe a similar change but from the HGV-A approximation to the hydraulic limit. In short, the flow changes from the HGV-A approximation where the value of $K$ is constant for a given value of $\theta /\alpha$, to a hydraulically controlled flow where the value of ${{Fr}}$ is constant for a given value of $\theta /\alpha$.

In general, the flows in the ${\textsf {L}}$ and ${\textsf {H}}$ regimes are very well described by the HGV-A approximation. For the ${\textsf {L}}$ experiment, ${{Fr}}=0.40$, ${{Re}}_g=398$ and $A=30$, meaning that ${{Fr}}{{Re}}_g A^{-1} \approx 5.3$ or $Re_g\approx 13.3 A$. We have seen in figure 6 that inertia can be neglected for this experiment. Furthermore, if we consider the series of experiments with $\theta =1^\circ$ ($K\approx 62$) shown in figure 8, we see that $Q_m\approx 0.5$ for ${{Fr}}^*/2>0.5$. This suggests that the limit of validity of the HGV-A (in which inertia is neglected) is around ${{Fr}}^*\approx 1$. For this set of experiments, ${{Fr}}^*=1$ corresponds to ${{Re}}_g\approx 62 A$. These two examples show that the restrictive condition ${{Re}}_g\ll A$ derived using a scaling analysis by Lefauve & Linden (2020a) to neglect advection terms can be drastically relaxed.

4.3. Regime transitions

In this section we finally discuss how the values of ${{Fr}}^*$ are related to regime transitions. Figure 9 shows the location of the different regimes in the parameter space $(\theta,{{Re}}_g)$ for all experimental set-ups with superimposed curves of constant ${{Fr}}^*$. It can be clearly seen that these curves describe the transitions between the different regimes through all values of $\theta /\alpha$. For $\theta /\alpha \gtrsim 1$, curves given by ${{Fr}}^*=\text {const.}$ are equivalent to curves given by $\theta {{Re}}_g=\text {const.}$ for a given duct. This is in agreement with the transitions derived for forced flows by Lefauve et al. (2019a) using an energetics budget approach and that Lefauve & Linden (2020a) derived using a two-layer frictional hydraulics model. However, the agreement of the curves of constant ${{Fr}}^*$ value, based on the HGV-A approximation, extends to lazy flows ($\theta /\alpha <1$) including $\theta /\alpha \leq 0$ where the curves given by $\theta {{Re}}_g=\text {const.}$ do not make physical sense.

Figure 9. Location of the different regimes in the (${{Re}}_g$,$\theta$) plane for the four different set-ups: mSID, LSID, HSID and tSID. The different symbols represent the different regimes, i.e. laminar (${\textsf {L}}$), Holmboe waves (${\textsf {H}}$), other waves (${\textsf {W}}$), intermittently turbulent (${\textsf {I}}$) and turbulent (${\textsf {T}}$). The solid lines represent curves of constant ${{Fr}}^*$ with the value indicated along the line. The dashed line represents $\theta /\alpha =1$. The dotted lines are examples of the transition between regimes for forced flows given by $\theta {{Re}}_g=\text {const.}$ as proposed by Lefauve et al. (2019a).

In the mSID set-up, the transitions between successive regimes are given by ${{Fr}}^*\approx 0.5,$ 1 and 2. The detailed analysis of the three experiments in the previous section (shown as red open squares in figure 9a) and of the mass flux as a function of ${{Fr}}^*$ coincides with this picture. We observe the emergence of turbulence at ${{Fr}}^*\approx 1$, where we find the limit of validity of the HGV-A approximation. This transition further represents a transition to hydraulically controlled flows. It is interesting now to note that experiments in the ${\textsf {I}}$ regime (i.e. for $1 \lesssim {{Fr}}^* \lesssim 2$) and $\theta \geq 2$ have values of $Q_m$ larger than in the hydraulic limit where $Q_m\approx 1/2$. For these experiments, $Q_m$ seems to overshoot following the trend of the HGV-A approximation ($Q_m={{Fr}}^*/2$). This suggests that the flow in the ${\textsf {I}}$ regime has mixed characteristics of the HGV-A approximation and the fully turbulent, hydraulically controlled flow. Finally, the characteristics of the HGV-A approximation are negligible when ${{Fr}}^*\gtrsim 2$ and the flow is in the ${\textsf {T}}$ regime.

Surprisingly, although curves of constant ${{Fr}}^*$ match the transitions for all ducts, the particular values of ${{Fr}}^*$ associated with the transitions differ. For the LSID set-up, the values of ${{Fr}}^*$ marking these transitions are about a factor two larger, for the HSID about a factor three larger, and for the tSID about a factor five larger. In spite of the differences in these values, there is a remarkable similarity: there is always a factor two between the values of ${{Fr}}^*$ describing successive regime transitions.

The difference in the transition values for the tSID set-up with respect to all the other set-ups could have been expected since there is one non-dimensional parameter ($B$) that is different. In fact, the tSID set-up is much thinner than the other set-ups and, hence, 3-D effects are more important. In particular, there is additional friction from the side walls, which results in a smaller value of ${{Fr}}$ for a given forcing magnitude. This, in turn, means that the successive transitions are delayed. Further investigation of the 3-D effects would require extension of the 2-D solution presented in § 3.3 to three dimensions, which we believe is possible, but outside the scope of the current paper.

On the other hand, the difference between the mSID and LSID (both with $A=30$), and the HSID (with $A=15$) set-ups needs further explanation since the HGV-A approximation accounts for having different $A$ values. It could be that $A=15$ is not long enough for the limit of long ducts to hold. However, we do not have enough information to test this hypothesis. Still, the most surprising difference is between the LSID and mSID set-ups, as already mentioned by Lefauve & Linden (2020a), since both set-ups have the same values for $A$, $B$ and ${{Sc}}$. If the values of these non-dimensional parameters are the same, the fact that two experiments with the same values of ${{Re}}_g$ and $\theta$ are in different regimes would go against the principle of dynamic similarity. Lefauve & Linden (2020a) already mentioned that there might be a missing length scale (and its accompanying non-dimensional parameter), but they could not think of a relevant one. If one is to consider just the duct itself, there are no other non-dimensional parameters that are needed to describe the problem. For example, the non-dimensional parameters used so far are the only parameters needed to set-up a numerical simulation such as those by Kaptein et al. (2020). This would suggest that there are certain characteristics of the set-up outside the duct that can delay or stimulate flow transitions. Unfortunately, there is much less data available for the LSID and HSID set-ups when compared with the mSID set-up. For example, there are no 3-D velocity and density fields, and there are fewer experiments for which the value of $Q_m$ is available. However, we saw that, for the mSID set-up, the different data give a self-consistent picture. Ultimately, new independent experimental or numerical efforts will be needed to clarify the differences between set-ups.

We believe that the key to understand further the regime transitions, and possibly, the differences between the mSID and LSID set-ups lies at the edges of the duct. There are three main reasons for this: (i) the flow is largely determined by flow details at the edges of the duct (e.g. the analytical solutions obtained using either two-layer hydraulics or the HGV-A approximation are highly dependent on the boundary conditions there); (ii) the HGV-A approximation formally breaks down immediately outside the duct; (iii) instabilities first emerge at the edges of the duct where the nonlinear terms are larger since both $u$ and $\partial u /\partial x$ grow as $|x|\to 1$ (this can be deduced from conservation of mass and the shape of the interface). Hence, further knowledge of the flow and the experimental set-ups at the edges of the duct is needed to clarify these points. Note, for example, that Meyer & Linden (2014) used two different shapes of edges in the duct, and any of these shapes are quite different from the solid long vertical walls in the simulations by Kaptein et al. (2020). Furthermore, more information about the flow at the edges of the duct can help us understand how it transitions from the HGV-A approximation to being hydraulically controlled.

In spite of the differences between the set-ups, figure 9 suggests that the parameter ${{Fr}}^*$ does quantify the forcing for a given set-up. Recall that this parameter groups together the control parameters: ${{Re}}_g$, $A$ and $\theta$ with the latter incorporated in the value of $K$. It was already known that ${{Re}}_g A^{-1}$ is the control parameter for horizontal ducts (Hogg et al. 2001). In this way, ${{Fr}}^*$ is an equivalent number that is only rescaled by dividing by $K$ to account for the relative inclination of the duct with respect to its internal angle (i.e. depending on $\theta /\alpha$). It would seem that $K$ is an intrinsic parameter of SID set-ups. A different way to analyse the problem is in terms of a generalized gravitational Reynolds number for inclined ducts. If the ${{Re}}_g$ quantifies the forcing due to the density differences, a generalized gravitational Reynolds number

(4.3)\begin{equation} {{Re}}_{g,\theta}= \frac{K(\theta=0)}{K}{{Re}}_g\approx \frac{131}{K} {{Re}}_g \end{equation}

can be defined, in a similar way as done by Kaptein (2021). In this way, ${{Re}}_{g,\theta }$ incorporates both aspects of the forcing: the density difference between the tanks and the tilting of the duct.

We have shown that the transition from laminar to turbulent in a (long) SID experiment with ${{Sc}}\gg 1$ depends solely on ${{Fr}}^*$. For weak forcing (small ${{Fr}}^*$ values), the flow is laminar and can be described by HGV-A approximation. As the forcing (the value of ${{Fr}}^*$) is increased, the flow develops nonlinear and time-dependent characteristics (waves, thinner boundary layers, turbulence). Once the flow becomes turbulent, the HGV-A approximation becomes invalid. The nonlinear characteristics result in additional vertical viscous momentum diffusion that, for ${{Fr}}^*>1$, compensates the increase in forcing. In this way, the value of ${{Fr}}$ can remain constant in agreement with (frictional) two-layer hydraulic theory (Gu & Lawrence 2005).

5. Conclusions

In the current paper we present an analytical solution for the laminar flow in SID experiments within the HGV-A approximation, where the acronym stands for hydrostatic/gravitational/viscous balance in momentum and advective for density. This approximation is derived for long ducts ($A=L/H\gg 1$) and small-inclination angles $\theta$ (including $\theta =0$ and slightly negative angles). The Reynolds number has to be small enough for viscosity to dominate over inertia in the along-channel momentum equation, but diffusion is negligible in the density transport equation. This combination is only observed if ${{Sc}} \gg 1$. Under these conditions, the flow is organized into two homogeneous layers and the non-dimensional volume flow rate is given by ${{Fr}}^*={{Re}}_g (KA)^{-1}$, where ${{Re}}_g$ is the gravitational Reynolds number and $K$ is a geometrical parameter that depends exclusively on $A\sin \theta \approx \theta /\alpha$ with $\alpha$ the internal angle of the duct. The parameter $Fr^*$ can be considered as a control parameter because its value can be set a priori, and it quantifies the relative magnitude of the forcing (gravity and the horizontal pressure gradient) against viscous forces.

A comparison with results from laboratory experiments has demonstrated that the analytical solution describes well the flow in the laminar regime and the time-averaged flow in the Holmboe waves regime. Furthermore, this comparison shows that the HGV balance holds and inertia can be neglected for values of ${{Re}}_g$ and $\theta$ well beyond what was previously thought. The HGV-A approximation presents the relevant balances in these two regimes, and it clearly exposes the two forcing mechanisms: a baroclinic pressure gradient due to the inclined interface and gravity in the along-duct direction due to the tilt of the duct. For $\theta /\alpha \gtrsim 1$ (forced flows), the inclination of the interface might be neglected, but for $\theta /\alpha \lesssim 1$ (lazy flows), it drives a crucial component of the forcing despite being small (e.g. $S'\approx 48/(131 A)$ for $\theta =0$). To make the link between horizontal and inclined ducts, it is crucial to include both of these terms.

Finally, the curves given by ${{Fr}}^*=\text {const.}$ describe the regime transitions observed experimentally for horizontal and slightly (both positively and negatively) inclined ducts. As the value of ${{Fr}}^*$ is increased, the forcing magnitude is increased. First, the HGV-A approximation breaks down as a critical ${{Fr}}^*$ value is approached and reached. Successive transitions from laminar flow, to interfacial waves, to intermittent turbulence and sustained turbulence are needed to increase momentum diffusion to keep the internal Froude number from exceeding the maximum value imposed by maximum exchange theory. In this way, the curves ${{Fr}}^*=\text {const.}$ provide a solution to the long-standing problem of finding the curves describing the regime transitions spanning horizontal and inclined ducts.