Stratified inclined duct: two-layer hydraulics and instabilities

Abstract The stratified inclined duct (SID) sustains an exchange flow in a long, gently sloping duct as a model for continuously forced density-stratified flows such as those found in estuaries. Experiments have shown that the emergence of interfacial waves and their transition to turbulence as the tilt angle is increased appears to be linked to a threshold in the exchange flow rate given by inviscid two-layer hydraulics. We uncover these hydraulic mechanisms by (i) using recent direct numerical simulations (DNS) providing full flow data in the key flow regimes (Zhu et al., J. Fluid Mech., vol. 969, 2023, A20), (ii) averaging these DNS into two layers, and (iii) using an inviscid two-layer shallow-water and instability theory to diagnose interfacial wave behaviour and provide physical insight. The laminar flow is subcritical and stable throughout the duct and hydraulically controlled at the ends of the duct. As the tilt is increased, the flow becomes supercritical everywhere and unstable to long waves. An internal jump featuring stationary waves first appears near the centre of the duct, then leads to larger-amplitude travelling waves, and to stronger jumps, wave breaking and intermittent turbulence at the largest tilt angle. Long waves described by the (nonlinear) shallow-water equation are interpreted locally as linear waves on a two-layer parallel base flow described by the Taylor–Goldstein equation. This link helps us to interpret long-wave instability and contrast it with short-wave (e.g. Kelvin–Helmholtz) instability. Our results suggest a transition to turbulence in SID through long-wave instability relying on vertical confinement by the top and bottom walls.


Introduction
Buoyancy-driven exchange flows between water masses of different density are common in estuaries and straits, e.g. the straits of the Great Belt, Gibraltar, Bab el Mandab, and the Bosphorous ( Gregg & Özsoy 2002).These essentially hydrostatic and two-layer flows are known to exhibit hydraulic jumps (Farmer & Armi 1988;Thorpe et al. 2018b), which are important discontinuities in the internal flow properties (layer thickness and speed) and often lead to instability (Lawrence & Armi 2022).A local jump can influence the large-scale properties of the flow such as the exchange flow rate (Lawrence 1993).These non-local hydraulic flow features are often studied in an idealised 'shallow-water' setting consisting of 2 A. Atoufi, L. Zhu, and others fluid organised in two counter-flowing, frictionless layers of specified thickness and speed with constant densities (Armi 1986;Lawrence 1990;Dalziel 1991).
In this paper, we employ two-layer hydraulics as a diagnostic tool to derive insights from direct numerical simulation (DNS) data of the exchange flows in the stratified inclined duct (SID).SID is a canonical stratified shear flow through a long, square cross-section, tilted duct for which there is now ample data, both experimental (e.g.Partridge et al. 2019) and numerical (Zhu et al. 2022).The investigation of turbulence in two-layer shear flows through tubes dates back to the classic experiments of Reynolds (1883) and Thorpe (1968), who both used a closed tube.The opening of the tube into large reservoirs in SID is more recent and allows for interfacial waves and turbulence to be sustained for much longer time periods, and for various flow regimes to be distinguished.The successive transitions to increasingly turbulent regimes in SID, as the Reynolds number and tilt angle are increased, have been recognised since Macagno & Rouse (1961) and Kiel (1991), and have been much studied more recently (Meyer & Linden 2014;Lefauve et al. 2019; Lefauve & Linden 2020;Duran Matute et al. 2023).These transitions are underpinned by many fundamental features of stratified flows, including interfacial waves and turbulent intermittency.
Our first aim with this new analysis is to uncover internal hydraulic effects in order to explain some of these leading-order dynamics that DNS and experimental data have revealed but not yet explained.In particular, we will provide the first direct evidence for the existence of internal hydraulic jumps, where the layer thickness expands in the direction of the flow of each layer.We also show how the development of jumps and large-amplitude interfacial waves coincides with a plateau in the exchange flow rate through the duct, after an initial increase with increasing Reynolds number and tilt in the laminar regime.This upper bound is a remarkably robust feature of SID, which has deep ramifications for the flow energetics and transition to turbulence (Lefauve et al. 2019).Though the emergence of waves and turbulence has long been linked to the notion of 'hydraulic control' in the experimental SID literature, this link has not yet been studied in detail, primarily because experiments lack the full velocity, density and pressure data along the length of the duct.The recent availability of DNS data (Zhu et al. 2022) overcoming these limitations finally makes a rigorous hydraulics analysis of SID possible.
Our second aim is to link two-layer internal hydraulic effects to the growth of instabilities and to shorter (non-hydrostatic) waves, links which are rarely found explicitly in the literature.The streamwise variation of the base flow in the SID geometry, and generally in all exchange flows, distinguishes them from idealised parallel stratified shear flows.This variation, however, is essential to the formation of internal hydraulic jumps, to the ideas of hydraulic control and maximal exchange, and hence to the nature of SID turbulence.It is well known that under some flow conditions, the loss of hyperbolicity of the shallow-water equations renders long waves unstable.However much less is known about the consequences of this instability, and the relative importance of long-wave versus short-wave instability.We will clarify this by explaining how a certain range of unstable nonlinear shallow water waves associated with an internal jump can be interpreted as linear instabilities on a locally parallel base flow, and thus that the insights derived from stable hydraulic theory remain valid even under (moderately) unstable conditions.
To tackle these aims, we introduce our DNS datasets in §2, and the two-layer averaging in §3, using the averaged datasets to show evidence of jumps and maximal exchange.In §4, we adapt the two-layer shallow water theory to SID, summarise important results from the literature, and connect them to linear stability theory.In §5, we use this hydraulics and stability framework to analyse our DNS.Then, in §6, we explore the transition between long (hydrostatic) waves and short (non-hydrostatic) waves, the influence of molecular diffusion (Prandtl number) and of smoother flow profiles.Finally, we draw our conclusions in §7.

Direct numerical simulations
2.1.Methodology Our DNS solves the following non-dimensional Navier-Stokes equations under the Boussinesq approximation for the density-stratified flow in the SID setup sketched in figure 1: where Re =  *  * / is the input Reynolds number ( * = √   * is the characteristic buoyancy velocity,  * is the dimensional half-height of the duct,  is the kinematic viscosity,  is the reduced gravity); Ri =   * /(2 * ) 2 is the input bulk Richardson number, giving a fixed Ri = 1/4; and Pr ≡ / is the Prandtl number ( is the thermal diffusivity).The unit gravity vector in the coordinate system (, , ) aligned with the duct is ĝ ≡ [sin , 0, − cos ].The DNS is performed with the open-source solver Xcompact3D (Bartholomew et al. 2020).
The duct has a square cross-section of non-dimensional height and width 2 and of length 60 (giving a long aspect ratio of  = 30).No-slip boundary conditions for  and no-flux boundary conditions for  are applied on the duct walls in the spanwise ( = ±1) and vertical ( = ±1) directions with an immersed boundary method.The flow is driven by the density difference between the dense ( = 1) and light ( = −1) fluid in the left-and right-hand reservoirs, respectively, producing counter-flowing layers in the streamwise -direction.The experimental reservoirs are modelled by ad hoc forcing terms   and   , which damp flow in the reservoirs and restore the density field to  = ±1, allowing us to maintain a quasisteady exchange flow for arbitrarily long times at a minimal computational cost.Details about the numerical setup and validation against experiments and benchmark cases (with large reservoirs and   =   = 0) can be found in Zhu et al. (2022).
The DNS is started at  = 0 from lock-exchange initial conditions, after which two counterflowing gravity currents develop from  = 0, advance at absolute speed ≈ 1, and reach either end of the duct after ≈ 30 advective time unit (ATU).Shortly after, the statistically-steady exchange flow of interest in SID becomes established.Conservatively, we only retain  80 in the following analysis to discard any initial transients, and we run the simulation until  = 260.
When the setup is tilted by an angle  > 0 • , the streamwise component of gravity contributes Ri  sin  to the -component of the momentum and adds extra kinetic energy

Database
We use the DNS database recently acquired by Zhu et al. (2022), which shows good agreement with experiments when all five non-dimensional parameters Re, Ri, Pr, ,  are matched.This database provides the complete set of flow variables all along the duct, which is a requirement to properly study hydraulic processes in SID.
We consider flows at a fixed Reynolds number Re = 650 and fixed Prandtl number Pr = 7, corresponding to temperature-stratified water.Four cases are examined, corresponding to tilt angles  = 2 • , 5 • , 6 • , 8 • , denoted by B2, B5, B6, B8, respectively in Zhu et al. (2022).Each of these four cases covers specific flow regimes: B2 is laminar (L), B5 has stationary waves (SW), B6 has travelling waves (TW), and B8 has intermittent turbulence (I).In the following, we refer to these datasets as L, SW, TW, and I respectively.We only briefly touch upon more turbulent cases (e.g.B10 having Re = 1000 and  = 10 • , referred to here as T), as they deviate significantly from the assumptions of the model in this paper, that the flow is primarily composed of two layers with a hydrostatic pressure field (discussed in §A).For more details about the flow regimes, statistics and spatio-temporal dynamics, see Zhu et al. (2022).

Layer averaging procedure
We seek to reduce the dimensionality of our DNS datasets to a set of two layers in order to interpret their dynamics using simple two-layer hydraulics, as sketched in figure 2. To do this, we will define the interface that separates the layers as the height  = (, ) where  = 0.The streamwise velocities, heights and densities of each layer are   , ℎ  , and   (where  = 1, 2 correspond to the upper and lower layer, respectively), are then obtained by averaging the DNS data over the -direction and in  over the height of each layer.Specifically, the flow properties of the layers are obtained by where the top-layer average is ∫  −1 •  and and the spanwise average is Recall that  = 1 and  = −1 are the non-dimensional height of the top and bottom walls, respectively.Figure 3 shows a single snapshot at time  = 110 of  (colours) and  (contours) at the  = 0 midplane Focus on Fluids articles must not exceed this page length for the two datasets L and TW, highlighting the  = 0 density interface with a thick green contour.
3.2.Layer-averaged DNS data and evidence of jumps Figure 4 illustrates the results obtained after applying this layer averaging to the DNS database.We show  −  diagrams of the lower-layer height ℎ 2 (top row) and lower-layer velocity  2 > 0 (bottom row) for the laminar (L), stationary wave (SW), travelling wave (TW), and intermittently turbulent (I) cases.The layer averaging (in , ) was performed in the range || 32, which includes the duct || 30 and extends slightly into the reservoirs where the layers flow in and out.
The L flow regime exhibits sudden changes in depth and speed only at the ends of the duct  = ±30, as is expected from the flow exiting into the deep reservoirs.Figure 4(a) and an instantaneous snapshot of the flow in figure 3(a) confirms that the interface () is steady and gently sloping down ( < 0) throughout the duct, and is symmetric in the duct centre, i.e. (, ) = () = −(−).
In contrast, for the SW and TW flow regimes, the upper layer thickness ℎ 2 suddenly increases along the direction of the flow (purple for  < 0 to orange for  > 0 in figure 4), which is accompanied by a sudden drop in  2 (dark to light red).This discontinuity indicates the presence of what is commonly called an 'internal hydraulic jump' (Baines 2016; Thorpe et al. 2018a), as can be seen in figure 3(b), where both layers experience a sudden expansion and deceleration in their respective flow direction.In the TW flow (figure 3(b)) the interface is sloping up ( > 0) in the vicinity of the jump, and is negative ( < 0) throughout most of the left-hand side of the duct, and vice versa, which is the opposite of what is found in the L flow (figure 3(a)).
Figure 5 shows the temporal evolution of the interface in all four flows, with (, ) plotted at intervals separated by one ATU and vertically stacked.In the SW case, the jump remains in the narrow interval  = ±1 whereas in the TW case, it oscillates over a much larger interval and sends off waves in either direction, hence the distinction between the 'stationary' and 'travelling' wave regime.In the I case, moving jumps are observed in the quiet phase (150 𝑡 200), being initiated near both ends of the duct at  ≈ 150 and progressively moving toward the centre of the duct.These jumps merge at around  = 200 and then stay

Flow rate and evidence of maximal exchange
To a good approximation, the flow in SID has zero net (barotropic) volume flow rate but a non-zero exchange volume flow rate (or volume flux) (, ) and mass flow rate (or mass flux)   (, ) defined as The approximation in (3.5) comes from the fact that the layer averages of the product  is not exactly equal to the product     of the layer averages.Also, recall that the non-dimensional density is defined such that the mean density is 0 and the minimum and maximum are −1 and 1 respectively.Hydraulic jumps in two-layer flows are often connected to the notion of maximal exchange, i.e. of an upper bound in the exchange volume flux  and mass flux   .This means that flows lacking such jumps have a lower , and that no realisable flow may have a higher .While hydraulic jumps have not yet been investigated in detail in the experimental SID literature, the mass flux   has, due to the simplicity with which it can be measured.
In figure 6 we show the dependence of the mass flux   with RiRe sin  ≈ (1/4)Re using its temporal mean (symbols) and extreme values (error bars) in 15 different DNS from Zhu et al. (2022), containing the four main datasets L, SW, TW and I, as well as 11 others including one turbulent dataset (T).We find that   increases approximately linearly with Re in the L regime (where little mixing ensures that  ≈   ), in agreement with the laminar analytical solution of Duran Matute et al. (2023).For higher values of Re it reaches an upper bound   ≈ 0.5 in the W regime, before dropping below 0.5 in the I and T regimes.These observations are consistent with the corresponding experimental data of Lefauve & Linden (2020) (their figures 5-6), if we exclude their data at small angles  < 2 • (which behave slightly differently and we did not simulate).This apparent upper bound is an evidence of maximal exchange, since further increase in  does not lead to an increase in the velocity difference between the upper and lower layer.It is a significant departure from the laminar solution based on the balance of gravitational forcing and viscous drag, suggesting   ∝ Re (Lefauve & Linden 2020).As it will be shown later (section 5.3) this   ≈ 0.5 is the threshold of the long-wave instabilities and onset of hydraulic transitions and a further increase in   beyond this limit is taxed by turbulent dissipation.
This large body of experimental and numerical evidence in favour of maximal exchange in SID and   ≈ 0.5 at the transition between the L and W regime, combined with our observations in §3.2 of the existence of an internal jump, is strongly suggestive that hydraulic effects dominate the flow in the W regime onwards.
In the next section, we develop the two-layer hydraulics framework to study these jumps and maximal exchange and their relation to the transition from laminar flow to waves and to turbulence in SID.

Two-layer equations: characteristics and instabilities
In this section, we aim to define and relate the characteristic velocity of the two-layer flows in the context of SID to the hydraulic regime of the flow.We also aim to provide a physical implication of when this characteristic velocity is not purely real.We start by simplifying the inviscid Navier-Stokes equations with shallow water (long wave) theory in §4.1-4.3 before comparing it to the Taylor-Goldstein (linear wave) theory in §4.4, and interpreting our findings in §4.5.

Shallow water equations: nonlinear long waves
The validity of this hydrostatic assumption is verified in our DNS data in appendix A. The upper layer then obeys and the lower layer obeys Here,   (, ) is the pressure at the upper wall and () is the elevation of the bottom wall.Conveniently, the bottom wall in SID is fixed at () = −1.We subtract the momentum equations in (4.1) to (4.2) to remove   and reduce the number of unknowns.The variation of density of the upper and lower layers in  is also neglected compared to variations in height of the layers.The shallow water equations become , with two auxiliary equations to satisfy the no-net (barotropic) flow condition and geometric constraint inside the duct (/ = 0) By taking the derivative of (4.5) with respect to , these four equations can be written compactly as where the state vector  and coefficient matrices A, C are The shallow water equation (SWE) (4.6) is quasilinear, i.e. linear in the derivatives of  but with the coefficient matrix A dependent on .The quasi-constant local density difference between layers is defined as Δ(, ) ≡  2 −  1 ∈ (0, 2) (specified by our layer-averaged DNS data).This equation does not assume that interfacial waves have infinitesimal amplitudes, but it does assume, through the hydrostatic assumption, that waves are long with respect to the layer height (i.e. that their non-dimensional wavenumber  1).In the following, we neglect the forcing  = [Ri sin  Δ, 0, 0, 0]  , (4.8) recalling that sin  1.We will study (4.6) in the homogeneous limit (  = 0), with a focus on the eigenvalues of A, in which  plays no role.The role of forcing in shallow water theory is an interesting question left for future work.However, the forcing proportional to sin  does influence the DNS, and is thus implicitly taken into account in the state vector  obtained after layer averaging the DNS data.

Characteristic curves and propagation of information
Consider a left eigenvector  and eigenvalue  associated with the matrix pair (A,C) such that   A =   C, where  denotes Hermitian transpose.Multiplying (4.6) by   yields The eigenvalues  are called characteristic velocities, since they define characteristics curves  in the (, ) plane along which the partial differential equation (4.6) reduces to an ordinary differential equation (Whitham 2011).For this homogeneous equation, the combinations of flow variables   C / = 0 is conserved.These characteristic velocities, henceforth simply referred to as 'characteristics', can be complex  =   +   ∈ C. Their real part   ≡ () represents the phase speed of shallow water waves, describing the trajectories .
Their imaginary part   ≡ () represents any potential growth (  > 0) or decay (  < 0) of these waves.The direction of information propagation is set by the sign of   : when   > 0, information propagates rightward (towards increasing ), and vice versa, whereas   = 0, indicates stationary waves.The characteristics  are given by the two distinct solutions of det(  − ) = 0: where (4.13) Here, (4.11) and (4.13) use the volume flux  > 0 in (3.4) (which is an approximation relying on (3.3)) and the interface position instead of the layer heights and velocities, as well as the fact that Ri = 1/4 and ℎ 1 + ℎ 2 = 2 in SID.We see that the characteristics consist of a convective velocity λ and a phase speed , which can be imaginary, depending on the value of the 'stability Froude number'  2 Δ (Long 1956; Lawrence 1990; Dalziel 1991).Note that (4.7), (4.10) and (4.12) are slightly modified versions of those given in previous studies (Long 1956;Armi 1986;Lawrence 1993) adapted to SID flows.
If  2 Δ < 1 the two characteristics  1,2 = λ ±  are real and information propagates in both directions relative to the convective velocity λ.The absolute direction of propagation is given by the sign of  1,2 , which we shall return to in §4.3.
If  2 Δ > 1 the characteristics become complex conjugates  1,2 =   ±   = λ ± ||, indicating that the system is no longer hyperbolic and that waves grow temporally unstable.The condition  2 Δ > 1 is sometimes known as Long's instability criterion (Long 1956), although it is quoted in Lamb (1932) and likely dates back to Helmholtz.The real part is the convective velocity, i.e. information only propagates in one direction, and the growth rate is 7(a-c) shows how  1,2 vary with the volume flux  and interface position  following (4.11),assuming no mixing (Δ = 2) and a horizontal duct (cos  = 1).We compare the case where the interface is locally symmetric ( = 0, panel b) to the cases where it is asymmetric, i.e. above or below the mid-level ( = ±0.5 in panels a and c, respectively).Panel b shows that with a symmetric interface  1,2 = ± √︁ 1 − 4 2 are real (red curves) for  0.5 and become purely imaginary (blue curves) for  > 0.5.However, when the interface is either below or above the midplane ( = ±0.5, panels a,c), this transition to instability occurs at a lower critical volume flux  >   ≈ 0.4.In summary, instability is caused, for a given volume flux , by an increasingly asymmetric interface ||, and vice versa, for a given interface position, by an increasing volume flux .This offers a possible explanation for the transition from L to W flow observed in figure 6.

Composite Froude number and hydraulic control
To further stress the importance of the characteristics  1,2 , we return to the original SWE (4.6), and note that a non-trivial steady solution () requires det A ≠ 0, i.e.
where  2 is the squared composite Froude number defined with the Froude numbers of the upper and lower layers,  1 , and  2 , respectively, Points at which  2 = 1 are called control points (Armi 1986;Lawrence 1990;Dalziel 1991).At control points, A is non-invertible and a regularity condition must exist to recover a steady solution.
The link between characteristics and the composite Froude number can be highlighted with the identity From this expression, we deduce the following, illustrated in figure 8.If the waves are stable ( 2 Δ < 1), the characteristics  1,2 are real.The flow is called subcritical and information propagates in both directions (along positive and negative ) i.e.  1  2 < 0 ⇔  2 < 1.In other words, the absolute phase speed || is larger than the absolute convective velocity λ in (4.10).Vice versa, the flow is called supercritical when information propagates in only one direction i.e.  1  2 > 0 ⇔  2 > 1, and the absolute phase speed || is smaller than the absolute convective velocity λ.Note that for supercritical flow, the direction of propagation associated with  1 and  2 is the same as that given by the convective velocity λ, i.e. the waves are swept downstream.
In the next sections, we clarify the interpretation of unstable shallow water waves using linear stability analysis around a locally parallel base flow.In §4.4 we take the long-wave limit of solutions of the Taylor-Goldstein equations, and in §4.5 we linearise the shallow water equations assuming the waves are sufficiently (but not excessively) long.

Taylor-Goldstein equations: linear short and long waves
We relax the previous restriction that waves must be long ( 1) by studying waves of possibly larger  but of infinitesimal amplitude developing on a parallel base flow described by a velocity profile U () and a density profile R ().The perturbation streamfunction ψ() exp  ( − ) describing the evolution of these two-dimensional linear waves is given by the inviscid Taylor-Goldstein equation (TGE) This can be analysed following standard methods (e.g.Drazin 2002; Smyth & Carpenter 2019) with details given in appendix B. Note that we assumed small tilt angles 0 < sin  cos , although the more general TGE in appendix B shows that sin  has a destabilising effect (ignored here).In the ordinary differential equation above, ' ' denotes differentiation with respect to , and  ∈ C is the phase speed of the plane waves, akin to the characteristics  of the SWE.However, we use a different notation for the characteristics of shallow water nonlinear waves and the phase speed of TG linear waves to emphasise that while the former implicitly assumes  1, the latter implicitly assumes   −1 .In other words, the TG linear waves we investigate should be shorter than the duct length  for the local analysis on a parallel base flow to be sensible.In other words,   −1 ensures that the waves do not 'feel' the streamwise variations of the base flow along the duct.
To make analytical progress and obtain solutions for , we assume a two-layer flow bounded by solid walls, with fixed layer heights ℎ 1 , ℎ 2 , velocities  1 ,  2 , and an interface at  = 0, consistent with the two-layer model adopted throughout this paper Note that ℎ 1 + ℎ 2 = 2 and, by a simple vertical shift, the above model is equivalent to having a domain restricted  = ±1 and an interface at an arbitrary  =  (giving ℎ 1,2 = 1 ± ).
By enforcing the matching conditions for the stream function and pressure at the interface (see appendix B), we derive the dispersion relation for the complex phase speeds , (4.20) These waves are the well-known Kelvin-Helmholtz (KH) waves supported by a single vortex sheet (see e.g.Smyth & Carpenter (2019) § 4.6.1).They are modulated by stratification; increasing Ri cos Δ always stabilises them.Importantly, the dispersion relation (4.20) describes waves in a domain bounded by the top and bottom walls, which as we will see, strongly affect waves that have a wavelength comparable to or longer than the domain height.
For 'short' waves, i.e. having a wavelength of the order of the duct height  =  (1) or shorter  > 1, the dispersion relation (4.20) cannot be simplified further, and these waves are dispersive.These will be referred to simply as KH waves, as in the short wave limit they are identical to those found in vertically unbounded domains.
For long waves (which, for clarity, we do not call KH waves), i.e. having a wavelength much longer than the duct height, but still shorter than the duct length  −1  1, the dispersion relation (4.20) simplifies as sinh  ℎ −→  ℎ and cosh  ℎ −→ 1 and (4.20) reduces to (4.10), i.e.
The dispersion relation (4.20) shows that there exists a smooth transition between short KH waves and long shallow water waves, as those waves differ by their wavelength but not the underlying physical mechanism.We return to this smooth transition and plot the dispersion relation in §6.
Although the limit (4.21) can be inferred from the PhD thesis of Gu (2001) ( § 3.3) and is briefly alluded to in Boonkasame & Milewski (2012) ( § 3), it does not appear to be widely disseminated in the hydraulics literature.

Link between complex characteristics and instability
This natural link between  1,2 and  1,2 in the long-wave limit can be further understood by considering another limit.It is possible to directly linearise the SWE (4.6) to study the evolution of infinitesimal perturbations  q(, ) (0 <  1) on a parallel, steady base flow  0 = ( 1 ,  2 , ℎ 1 , ℎ 2 ) akin to (4.19).We perform a first-order Taylor expansion of the coefficient matrix from (4.7) A() = A( 0 +  q) and obtain At order  we have the linear, local SWE where A 0 ≡ A( 0 ) is the local constant coefficient matrix.Importantly, the right-hand side coming from the Taylor expansion, and acting as a forcing term, vanishes when we assume the base flow to vary slowly.Substituting the plane wave ansatz q = q exp  ( − ) gives the phase speed  as the eigenvalue of the matrix pair (A 0 , C).The two distinct solutions  1,2 of det(A 0 − C) = 0 then become identical to the local characteristics  1,2 (at a fixed , ) derived in (4.10).
In other words, the nonlinear shallow water wave ( 1) characteristics  1,2 can be interpreted as the linear shallow water waves phase speed  1,2 propagating on a locally parallel base flow (  −1 ), which are themselves the long wave, non-dispersive limit case of the Taylor-Goldstein two-layer phase speeds  1,2 .This result ultimately stems from the fact that, simply put, the linearisation and the long wave limit commute.In particular, the potential positive imaginary component of characteristics   > 0 can be interpreted as the exponential growth rate of such unstable waves satisfying  −1  1.This interpretation is summarised in figure 9.The long aspect ratio  of SID (in this paper  = 30), ensures that the range of waves  −1  1 exists, and therefore that this interpretation is useful, unlike in shorter geometries having  10.We also note in passing that this interpretation can be generalised to any number of layers greater than two.

Two-layer hydraulics applied to DNS
In this section, we use the modelling results of §4 to understand the observations of hydraulic jumps and maximal exchange of §3.In §5.1 we study the stability Froude number flagging long-wave instability, in §5.2 we focus on the characteristics and composite Froude number to diagnose internal jumps and hydraulic control, and in §5.3 we explain the observed flow rate with notions of maximal exchange, viscous friction and mixing.

Stability Froude number and instability
In figure 10 we plot the spatio-temporal diagram of the stability Froude number  2 Δ (, ) given by (4.12) in our four datasets L, SW, TW, and I (all at  = 650) to diagnose any long wave instability.
The laminar flow (L) has  2 Δ < 1 everywhere (in blue in figure 10a), hence long waves are stable at  = 2 • .In all other cases,  2 Δ > 1 (in red in figure 10) in most of the duct, hence long waves are unstable at  = 650 for  >   where   ∈ [2 • , 5 • ].In the wave flows (SW and TW), the waves are most unstable (maximum  2 Δ , deep red) near the centre of the duct.In the I flow, long waves are very unstable ( 2 Δ 1) throughout most of the duct.These  2 Δ (, ) diagrams correspond to the absence of a jump in the duct in the L flow (figure 4a) and the existence of a jump in the SW, TW, and I flows (figure 4b-d), although the reasons for the jump remain to be explained.These diagrams also show that as the tilt angle  is modestly increased between 5 • and 8 • , further changes take place as the flow becomes increasingly unstable to long waves.Next, we delve deeper into the hydraulics analysis to understand jumps by focusing on the characteristics in each flow.

Characteristics, composite Froude number and control
In figure 11 we plot the characteristic velocities  1,2 () given by (4.10) (real parts in panel a and imaginary parts in panel b) and the composite Froude number  2 () given by (4.16) (panel c) at time  = 110, to diagnose the propagation of information and criticality of the flow, respectively.In figure 12    Note that we only show data inside the duct up to  = ±30.
opposite signs throughout most of the duct (  1   2 < 0 ⇔  2 < 1), allowing characteristics to cross.Information propagates in both directions and the flow is subcritical inside the duct.The speed of propagation is of order 0.2 − 0.4, i.e. significantly lower than the advective velocity 1.However, at the ends of the duct, the characteristics vanish locally (  1   2 = 0 ⇔  2 = 1 at || ≈ 30) signalling that long waves become stationary.Figure 11(a) shows that immediately outside the duct, information propagates only in one direction ( 2 > 1 at || 30), in fact, leftward on the left-hand side (  1 ,   2 < 0) and rightward on the right-hand side (  1 ,   2 > 0), i.e. always away from the duct.The ends of the duct, therefore, act as control points, in the sense that no information from the reservoirs can propagate into the duct.The existence of two such hydraulic control points, with their respective characteristics pointed outwards, means that the interior of the duct is 'fully controlled' in the hydraulic sense, isolating the flow from hydrostatic disturbances within the reservoirs to either side.This is the first direct evidence that SID flows in the L regime are hydraulically controlled.
In contrast, in the unstable SW, TW, and I cases, the roots are complex conjugates throughout most of the duct (figure 11(a,b)), a consequence of instability ( 2 Δ > 1), causing supercriticality ( 2 > 1, figure 11(c)).These unstable waves always move at the local convective velocity of the flow   = λ(, ), which we recall from (4.11) is non-zero if the interface is not at mid-depth (, ) ≠ 0.
In the SW and TW cases, the unstable waves are initially carried rightward ( λ > 0) throughout most of the left-hand side of the duct, and leftward ( λ < 0) throughout most of the right-hand side of the duct, as seen in figure 12(b,c).Thus, all unstable waves are carried towards the centre ( = 0) where their characteristics converge, creating the hydraulic jump observed in figure 5.The largest values of   = || are found in the region where the   = λ convective components converge (see figure 11(b) and green-yellow shades in figure 12(b,c)).Their growth rate is fast (  ≈ 0.2 − 0.5) and slightly higher in TW than in SW.Such a jump, bounded by supercritical regions on either side, is distinguished from the standard hydraulic jumps which make the flow transition from a supercritical to a subcritical state.It can be viewed as the limit of the length of the subcritical region tending to zero.The jumps in SW and TW can be called 'undular jumps' because of their moderate 'upstream' Froude numbers  2  ≈  2 /2 of each layer (between 1 and 2) and the small energy they dissipate, compared to direct, breaking hydraulic jumps.Some jumps in I are stronger, as evidenced by their locally higher  2 values and their visibly higher dissipation.
This pattern of unstable characteristics converging toward the centre to form a jump can be summarised by sign   = −sign .This can be understood first by sign   = sign λ = −sign  from (4.11), i.e. the waves are carried at the local convective velocity, and second, by sign  = sign , i.e. the interface does not slope down as in L but is instead lowered on the left-hand side of the duct, and lifted on the right-hand side (central jump).
The main difference between SW (stationary wave regime) and TW (travelling wave regime) lies in the behaviour of   near the jump around  = 0 (figure 11(a)).While   () goes smoothly through zero in SW, it oscillates more in TW, suggesting that the location of the jump is prone to oscillations.This is confirmed by comparing the  −  trajectory of the locus of the maximum  2 Δ in figure 10(b,c) or the maximum   in figure 12(b,c) as a proxy to the location of the jumps.
Finally, we turn to the I flow, which exhibits more vigorous interfacial turbulence and wave instability (especially between  = 150 − 250) and greater variability in  and  than SW and TW.The characteristics in figure 12(d) converge quickly to form a large number of local jumps around  ≈ 100, which then organise into three main clusters: a central stationary cluster flanked by a left and a right cluster which themselves converge to form discrete jumps around  ≈ 160 while being carried to the centre, eventually converging into a single jump  200.The convergence of characteristics tends to coincide with the maximum instability (  ≈ 0.5).During the more stable (transitional) period at  = 110 the multiple sign reversals of   () (figure 11(a)) hinder a straightforward interpretation of wave propagation along .
This pattern by which unstable waves are carried in SW, TW and I flows differs so greatly from the classical picture of hydraulic control in L flow that it prompts the question: since information travels toward the duct centre, does it travel from the reservoirs into the duct, and is the flow still hydraulically controlled? Figure 11 shows that close to the duct exits (|| ≈ 28), the composite Froude number  2 (panel c) of all the cases becomes 1.Meanwhile, immediately outside the duct 30 < || < 32 the waves are stable (  = 0, panel b), and both the curves of the real characteristics (panel a) and of the composite Froude number  2 (panel c) closely match those in the L flow.In other words, the SW and TW flows also have a A. Atoufi, L. Zhu, and others control point ( 2 = 1), flanked by narrow regions of subcriticality ( 2 < 1).We conclude that the flow within the duct remains isolated from the reservoirs , and hence that it is also hydraulically controlled.This represents the first direct evidence that SID flows in the W and I regimes, in addition to having a (supercritical-to-supercritical) jump in the centre of the duct, are also hydraulically controlled at the ends of the duct.

Maximal exchange and critical flow rate
In the previous section, we showed that all four flows cases (L, SW, TW and I) were hydraulically controlled in the sense that control points at the ends of the duct isolated the flow within the duct from processes in the reservoirs and prevented the flow inside the duct from reaching velocities exceeding a maximal volume flux .In this section, we seek to explain the differences in the value of this critical volume flux (and by extension mass flux) observed between the L, W, and I regimes in figure 6.
In the stable L flow, we find a time-averaged  ≈ 0.31 well below the absolute upper bound of 0.5 for instability given by (4.13).This is understood by the frictional hydraulic theory of Gu & Lawrence (2005), subsequently adapted to SID in Lefauve & Linden (2020) (their § 5.2).In short, the relatively low values of the Reynolds number in these low-tilt L flows mean that viscous friction at the duct walls and at the interface must be parameterised in the shallow water equations.This parameterisation allows a correct prediction of the sloping interface () (a consequence of viscous friction, i.e. loss of momentum along the flow of each layer), which in turn allows prediction of  by imposing the criticality condition at the ends of the duct  2 ( = ±) = 1.Simply speaking, the lower  and the longer the duct aspect ratio , the more friction occurs along the duct, the more offset the interface || becomes at the ends of the duct, and the lower the volume flux  becomes to satisfy  2 (||, ) = 1 (since  2 is an increasing function of both || and ).
In the unstable SW, TW and I flows, despite the existence of viscous friction, we find a remarkably consistent time-averaged  ≈ 0.51 − 0.53, slightly above the critical   = 0.5 upper bound for frictionless two-layer hydraulics and a flat interface  = 0 (figure 7).These values require a different explanation.Although the Reynolds number of SW, TW, and I are identical to L, their larger tilt angle  pushes these flows beyond the instability threshold  2 Δ = 1, corresponding to the transition between 'lazy' and 'forced' flows (Lefauve et al. 2019).However,  does not continue to increase with .Rather, beyond the instability threshold, (4.13) suggests that, in the centre of the duct (where  ≈ 0),  = 0.5 √︁ (Δ)/2 cos  Δ , i.e. a linear increase with the stability Froude number.We deduce that, since  never greatly exceeds 0.5 (the value reached the instability threshold), the subsequent increase in  Δ > 1 (indirectly caused by the forcing ∝  in the DNS) must be compensated by a decrease in Δ/2 ∝ 1/ 2 Δ , i.e. by increased mixing.The data show that the average Δ/2 indeed decreases from 0.79 in SW, to 0.76 in TW, to 0.68 in I.This mixing in turn explains why   (roughly ≈ (Δ/2)) stays robustly below 0.5 in SID, and indeed decreases from the W to the I regime (figure 6).

Applicability of unstable hydraulics
In this section, we study the applicability of the previous results to problems not usually considered in two-layer hydraulics.In §6.1 we study the transition between long waves (the propagation of which is identical to the local characteristics of the SWE) and shorter Kelvin-Helmholtz waves (only predicted by the TGE).In §6.2 we study the waves diagnosed from DNS run at different values of the Prandtl number to investigate their indirect dependence on scalar diffusion and the thickness of the density interface.In §6.3 we study how the growth  of long waves is impacted by smooth, diffuse (i.e.not strictly two-layer) density and velocity profiles, which are expected in all real-world flows (having a finite Re and Pr).
6.1.Long vs short waves Figure 13 shows the dispersion relation from the TGE (4.20) with the phase speed ( 1 ) =   1 (blue to red contours) and growth rate ( 1 ) =   1 (colour map with deep blue being stable) as functions of the wavenumber  (vertical axis) and the volume flux  (horizontal axis).We compare a symmetric interface ( = 0, panel a), an asymmetric interface ( = −0.5, panel b) and a case with a symmetric interface but without solid top and bottom boundaries (panel c), whose dispersion relation (B 13) is derived analytically in §B.3.
We recall that for  1 the phase speed   and growth rate   of TGE become identical to   and   of SWE, respectively.In this case the  = 10 −2 data of figure 13(a,b) become indistinguishable from those plotted in figure 7(b,c), respectively.
From figure 13(a-b) we recover the results from previous sections.First, long waves are non-dispersive (the phase speed contours do not depend on ).Second, they become unstable (lighter shades of blue and green) above a critical volume flux  >   , equal to 0.5 for a symmetric interface (panel a) and lower than 0.5 for an asymmetric interface (panel b).Third, for a symmetric interface, all unstable waves (long and short) are stationary (absence of contours), but all stable waves are travelling (presence of contours).For an asymmetric interface, even unstable waves are travelling in the reference frame of the duct, as we explained in §5.2.Fourth, the transition between short and long waves is smooth, i.e. there is a continuity between the long shallow water waves controlling the hydraulics of two-layer flows and the short KH waves.
Panels (a,b) of figure 13 also give new results.First, short waves ( 1) become unstable at smaller values of  compared to the long wave threshold   .This transition to short waves becomes noticeable from  10 −0.5 ≈ 0.3 and is clear for  > 1.The shortest waves shown here ( = 10 2 ) are predicted to become unstable above a very small volume flux  0.1.However, we note that this threshold would be closer to the long-wave   if we included viscosity in the TGE, as viscosity would significantly damp the growth of short waves.These panels show that for a given value of , the growth rate increases monotonically with  (i.e. the shortest waves are the most unstable).As is often the case, this 'ultraviolet catastrophe' would be regularised by viscosity, with the possible existence of a maximum growth rate at an intermediate  for intermediate values of Re.
In figure 13(c), the absence of solid walls does not affect short waves (the colours and contours at  3 are identical to panel a), because they do not 'feel' the presence of the walls.However, the absence of solid walls precludes the existence of long waves  0.3, because this setup (despite being bounded at  = ±1) approximates layers of infinite depth, compared to which all waves are 'short'.In other words, this analysis explicitly shows that the presence of solid walls in SID are crucial to explain the leading order dynamics in the DNS by allowing long-wave instability.Adding walls (panel a) creates the long waves on which hydraulic effects rely, an long waves transition smoothly into short waves as  increases.
Figure 14 shows the linear growth rate   obtained by substituting into the TGE dispersion relation (4.20) (function of ) the two-layer properties (as functions of ) extracted from the DNS.Diagnostics are shown for the L flow (panel a) and the TW flow (panel b) at time  = 110.The TW flow is, unlike the L flow, unstable to long waves  1, with the maximal growth rate found near the centre of the duct, as previously seen in figure 11(b) as we know from (4.21) that   ( 1) −→   .However, we also find that the L flow appears mildly unstable to short waves (especially very short waves  10), and the TW flow appears even more strongly unstable to them.However, we know that in the DNS the L flow is visibly stable and does not have observable interfacial waves, while the TW flow is primarily unstable to long waves, whereas short waves play a more minor role.This suggests that, at least for the present values of Re = 650 and Pr = 7, the growth of short waves is sufficiently damped by viscosity, mass diffusion and/or other effects not taken into account in this inviscid two-layer model.

Low vs high Prandtl numbers
Although Pr does not appear explicitly in the SWE, the DNS dynamics from which the twolayer properties are extracted certainly depend on Pr.In applications, three typical values are of particular interest: Pr ≈ 1 (representative of temperature stratification in air), Pr ≈ 7 (representative of temperature stratification in water, as studied in this paper) and Pr ≈ 700 (representative of salt stratification in water).
To study the impacts of diffusion, we carried out two additional DNS with parameters identical to the TW flow (with Pr = 7) at Pr = 1 and Pr = 28 (the latter requiring a much higher spatial resolution, hindering the study of higher Pr).In figure 15 we compare the characteristics curves  () and growth rates   at these three different values of Pr using the same visualisation as figure 12 (where the Pr = 7 data was already shown as TW).We find In figure 16 we compare the characteristics   () (panel a) and   () (panel b) as well as the composite Froude number  2 () at  = 110 (the Pr = 7 data was already shown in figure 11).All curves (solid red, dashed blue, and dotted green) have essentially the same qualitative features described earlier in figure 11.However, as noted in figure 15, the growth rate appears to decrease slightly with Pr.
The synoptic features of the flow governed by long waves, therefore, appear relatively unaffected by Pr.This can be rationalised by the fact that Pr will primarily influence the thickness of the density interface separating the two layers, rather than its location  (the locus of  = 0) or the speed of the flow (), which are the two key model variables in the SWE.
We expect short waves to be more strongly influenced by a decreasing thickness of the density interface with increasing Pr, and vice versa.However, the short waves observed in our DNS at Pr = 28 and in experiments at  = 700 are Holmboe waves, not the KH waves supported by our two-layer model.Unlike the KH instability caused by a single vortex sheet, the Holmboe instability is caused by the resonance of vorticity waves (e.g. on the edges of a diffuse velocity interface) with a non-collocated gravity wave (on a sharper density interface) (Carpenter et al. 2011).Lefauve et al. (2018) performed a linear stability analysis on the experimentally measured mean flow (at Re = 440, Pr = 700), including viscosity and scalar diffusion.They found that intermediate 1  2 Holmboe waves were most unstable (see their figure 6a).However, tackling Holmboe waves -and their presumed coexistence with the long waves governing hydraulic processes at high Pr SID -would require a three-layer model (for velocity) mixed with a two-layer model (for density).

Sharp vs smooth two-layer flow profiles
Finally, we study the influence of smooth density and velocity profiles U () and R () on the growth rate of long waves.To do so, we solve the eigenvalue problem from the Taylor-Goldstein equation (4.18) before the two-layer base flow ansatz (4.19).Numerical solutions for the growth rate   are shown in figure 17 as functions of .In panel a, we show the results for hyperbolic-tangent U ()/ = R () = tanh / where the interface thickness is progressively decreased from 1/80 (almost exactly two layers, solid lines) to 1/8 (dashed lines) to 1/4 (thicker interface, dotted lines).We compare these 'smooth tanh' growth rates (in black) to the equivalent 'sharp two-layer' growth rates (in red) obtained from the analytical dispersion relation (4.20) by layer-averaging the tanh profiles (in which case   =   ).In panel b, we keep the same density profiles but use U ()/ = − sin , which is a good approximation of the mean velocity at these relatively low values of Re =  (10 2 − 10 3 ).
Both panels a and b show that the   = 0.5 threshold for long-wave instability is virtually unchanged by smooth profiles, with only a slight increase of a few percent for the thickest interface  = 1/4.This result supports the relevance of two-layer hydraulics even in 'realworld' flows which depart significantly from the two-layer model.
However, the 'smooth tanh' growth rates are always lower than the corresponding 'sharp two-layer' growth rates.The thicker the interface, the lower the 'smooth tanh' growth rate.Comparing the vertical scale in panels a and b, we conclude that the sinusoidal velocity profile is more stable (by approximately a factor of 10) than the tanh profile.The combination of a sinusoidal velocity and a diffuse velocity interface (dotted blue line in panel b) yields the slowest growth as  increases.As such profiles are a better approximation of the mean flow of the DNS at low Pr (e.g.Pr = 1) than sharp two-layer profiles, these results warn us not to interpret the large growth rates   =  (0.1) found in this paper too literally.
In other words, although the qualitative predictions of two-layer hydraulics are robust (in particular the long-wave instability threshold   ) when the underlying data is not exactly two-layer, the quantitative growth rates predictions are over-estimated when the interface is diffuse, as in low-Pr flows.

Conclusions
In this paper, we employed a two-layer averaging procedure to extract a reduced-order representation of four direct numerical simulations (DNS) datasets in the stratified inclined duct (SID) at Re = 650 and Pr = 7 (with two supplementary datasets at Pr = 1 and Pr = 28.This two-layer representation revealed in §3 that the flow is stable in the laminar regime (L, tilt angle  = 2 • ), but develops an internal hydraulic jump (discontinuity in the layer properties) in the centre of the duct in the stationary wave regime (SW,  = 5 • ).This jump moves around in the travelling wave (TW,  = 6 • ) regime, and causes further disorganised wave breaking in the intermittently turbulent (I,  = 8 • ) regime.

Modelling results
To understand these findings, in §4 we adapted to SID DNS the well-known inviscid Boussinesq shallow water equations (SWE) governing the nonlinear evolution of long waves ( 1) at a sharp density interface.The SWE predict that information propagates along a pair of trajectories,  1,2 , that arise from the solution of a generalised eigenvalue problem and depend on the local () and instantaneous () state of the two-layer representation.The solutions can be written in the form  1,2 = λ ± , where the convective velocity λ is always real but the phase speed  may be either real or imaginary.When  1,2 are real ( ∈ R), they represent the propagation of two (neutrally stable) kinematic waves where λ can be interpreted as a convective velocity and  as the phase speed of waves relative to λ.The respective signs of  1,2 determine the direction of information propagation and whether the flow is subcritical (composite Froude number  2 < 1; product  1  2 < 0) with information propagating in both directions, or supercritical ( 2 > 1;  1  2 > 0) with information propagating only in the direction given by the sign of λ.When  1,2 are complex  1,2 =   ±   = λ ± ||, the real part still represents a convective velocity that carries information while the positive imaginary part indicates that the flow is unstable.Although in this unstable case the SWE are no longer hyperbolic, the flow may be viewed as supercritical ( 2 > 1) in the sense that information is propagated only in the direction given by λ.
To interpret the unstable SWE, we compared the characteristics with the dispersion relation from the inviscid Taylor-Goldstein equation (TGE) governing the linear normalmode stability of a two-layer base flow.We showed that the global nonlinear characteristics (, ) defined on the non-parallel base flow and sloping interface of the SWE could be interpreted locally as the phase speed and growth rate of linear waves in the long-wave limit (i.e. 1), propagating on a base flow that is locally assumed parallel.Importantly, this interpretation is only valid for waves that are much shorter than the duct length 2 (i.e.  −1 ).It provides a local, linear stability interpretation for unstable two-layer wave characteristics satisfying  −1  1, which is a relevant range in long ducts ( −1 1).The dispersion relation for the dispersive TGE waves () tend to the non-dispersive SWE characteristics  as  1, but they also allow us to explore the smooth transition to shorter ( 1), non-hydrostatic Kelvin-Helmholtz (KH) waves.

Physical results
Applying this two-layer hydraulics and instability framework to the two-layer-averaged DNS datasets yielded the main physical results of this paper in §5.We provided the first direct evidence that SID flows are, in all regimes (L, SW, TW and I), hydraulically controlled at the ends of the duct and thus in a state of maximal exchange.At these control points, the flow is locally supercritical; thus information from the reservoirs cannot enter the duct and influence the flow within it.In the SW, TW and I regime, the flow in the duct is always unstable to long waves ( 2 Δ > 1) and thus supercritical ( 2 > 1), explaining the existence of an undular jump within the duct, as a consequence of characteristic trajectories converging to a single point.In the I regime, multiple local jumps gradually merge into clusters and eventually into a single, stronger jump, resulting in greater instability.
The emergence of unstable, supercritical flow in SID beyond a certain tilt angle is rationalised by the fact that gravitational forcing continuously provides a surplus of kinetic energy which must be dissipated.From a hydraulics perspective, the required dissipation in a supercritical flow (i.e.having a surplus of kinetic energy compared to potential energy) must be associated with an decrease in kinetic energy and an increase in potential energy, hence a thickening of both layers downstream of the jump.The physical insight of energy surplus and dissipation dates back to Meyer & Linden (2014).It was later formalised by Lefauve et al. (2019) and Lefauve & Linden (2020), who showed using frictional two-layer hydraulics with a tilt  that the mid-duct interfacial slope obeyed  (0) ∝  − , where  represents viscous friction along the duct.The transition from subcritical to supercritical flow that we identified corresponds to the transition from 'lazy' to 'forced' flows, which they identified based on the relative importance of the tilt  and the duct geometric angle  = tan −1  −1 ≈ 2 • in this paper.'Lazy' flows are characterised by  < , and an interface gently sloping down.'Forced' flows are characterised by  > , and a relatively flat interface all along the duct.In forced flows, the tendency of  to tilt up the density interface exceeds the tendency of frictional losses  to tilt it down, thus  (0) > 0, which causes central jumps.
Next, we rationalised the values of the volume flux in all regimes.The value  ≈ 0.3 in the L regime (lazy flow) is explained by the offset of the interface at the ends of the duct where control ( 2 = 1) takes place, recalling that the sloping interface is caused by viscous friction along the duct.The robust values  ≈ 0.5 in the unstable SW, TW, and I regimes (forced flows) are all surprisingly close to the instability threshold   = 0.5 for a symmetric interface.Increasing instability from SW to TW to I as the tilt angle  is increased (which theory predicts should increase  above 0.5) appears balanced by increasing mixing in the layers (which decreases ).This explains why the maximal exchange threshold  = 0.5 predicted by inviscid long wave theory remains a remarkably robust feature of SID flows, even under turbulence.
Using the TGE analysis provided further physical insight into the applicability of unstable hydraulics in §6.We showed that short inviscid KH waves are always predicted to be more linearly unstable than long waves, despite the fact that long waves cause the internal hydraulic jumps observed in SW, TW and I and appear to dominate the dynamics of these SID flows.We explained this paradox by the neglect of viscosity in TGE, which would damp the shortest waves.We also showed that DNS at lower Pr = 1 or higher Pr = 28 showed qualitatively (but not quantitatively) similar two-layer long wave hydraulics to Pr = 7.We also highlighted that experimental observations at Pr = 700 of the simultaneous existence of unstable long waves (causing an internal jump and supercritical flow) with short finite-amplitude Holmboe waves could not be explained by the two-layer model, because it does not support Holmboe waves.Finally, we showed that the predictions of long wave instability (especially the threshold   = 0.5) were robust even in diffuse two-layer flows having a thick interface.

Outlook
These key 'hydraulic' features of SID flows, explaining the emergence of waves, increasingly supercritical jumps, and ultimately turbulence, result from long wave instability which (tautologically) relies on the existence of top and bottom solid boundaries confining the flow in a long, tilted duct.This 'long-wave' pathway to turbulence in SID appears a priori to differ from the classical 'short-wave' KH pathway in an unbounded stratified shear layer (see e.g.Caulfield & Peltier (2000); Mashayek & Peltier (2012)), often used as a paradigm for ocean mixing.Further work is needed to clarify the relative importance of long and short waves, and within short waves, of Kelvin-Helmholtz (two-layer) waves and Holmboe (three-layer) waves, and how they contribute to the transition to turbulence under varying , Re, Pr.
The current formulation of the two-layer SWE does not account for the mixing layer that develops and appears to be important beyond the wave regime.The omission of a mixed layer may reduce the accuracy when investigating detailed spatial features of the flow, such as the localization of unstable wave regions in the centre of the duct in TW and SW.This model is able to predict the formation of shocks and provide clues to the long-length-scale dynamics and how they govern the synoptic features of the flow, but further work is needed to accurately represent the non-hydrostatic processes of turbulent mixing itself.

Figure 1 :
Figure 1: A schematic of the stratified inclined duct (SID) numerical setup.The duct is oriented at angle  to the horizontal, which is equivalent to tilting the gravity vector.Densities are non-dimensionalised by the density differences between the reservoirs and lengths by the half-duct height.The duct's volume is (, , ) ∈ [−30, 30] × [−1, 1] × [−1, 1].

Figure 2 :
Figure 2: A schematic of two-layer shallow water flow.Flow in a part of the reservoir is considered.

Figure 3 :
Figure 3: Snapshot at  = 110 of streamwise velocity  (blue to red shading) and density  (colour contours) in the midplane ( = 0) for the () laminar (L) flow and () travelling wave (TW) flow, the latter showing evidence of an jump (see the interface  = 0 emphasised by the thick green line).

Figure 5 :
Figure 5: Temporal evolution of the density interface (, ) in (a) L, (b) SW, (c) TW, (d) I.The interface curves are stacked in time at intervals of one ATU.Jumps are revealed in (b)-(d) by the discontinuity in ().

Figure 6 :
Figure 6: The mass flux   from the 15 DNS of Zhu et al. (2022) increasing with  and  until the upper bound of ≈ 0.5.The symbols, coded by the respective regime, show the time-averaged value while the error bars depict the extreme values

Figure 7 :
Figure 7: Real and imaginary parts of the eigenvalues of the two-layer SWE in (4.11) at (a) positive (b) zero and (c) negative interface elevations  as functions of the volumetric flow rate .The eigenvalues are always complex for    = 0.5, and become complex for slightly lower critical values   ≈ 0.4 in the asymmetric cases (a,c).
Figure 9: (a) Summary of the interpretation of characteristics in §4.2 using linear theory, either by first linearising the NSE, and then taking the long-wave limit (TGE, §4.4) or vice versa (SWE, §4.5).Both approaches yield the same result for waves much shorter than the duct length but much longer than the duct height (range sketched in b).
we plot a set of discrete trajectories (white curves) by solving   =   (, ) for  > 80, (5.1) and initialising  ( = 80) as 61 equidistant points along the duct  ∈ [−30, 30].Additionally, we plot the local growth rate   (, ) in background colours (dark blue to yellow) for comparison.Figures 11(a) and 12(a) show that in the stable L flow,  has two distinct real roots of

Figure 11 :
Figure 11: Characteristics (a) real part   , (b) imaginary part   (only the positive values are shown), and (c) composite Froude number  2 of the L, SW, TW, I flows at  = 110.Note that we also show data immediately outside the duct, up to  = ±32.

Figure 12 :
Figure 12: Characteristic curves  () (in white) obtained by (5.1) for (a) L, (b) SW, (c) TW, (d) I.The curves (in white) originate from 61 equally spaced positions between  = −30 and 30.The colour contours show the growth rate   (, ) (which is zero in a).Note that we only show data inside the duct up to  = ±30.

Figure 13 :
Figure 13: Dispersion relation of all (long and short) inviscid two-layer waves: growth rate (colours) and phase speed (contours) solutions of the TGE varying with wavenumber  and volume flux .(a) Symmetric interface  = 0. (b) Asymmetric interface  = −0.5.(c) Symmetric interface but without solid top and bottom boundaries at  = ±1, in which case the long waves of SWE disappear.

Figure 14 :
Figure 14: Dispersion relation (growth rate only) predicted by the TGE two-layer model applied to the DNS (a) L and (b) TW flows.The dashed line at  = 1 represents the boundary between long and short waves.Short waves are predicted to be most unstable by this inviscid model but appear relatively stable in reality.

Figure 17 :
Figure 17: Growth rate   of long waves on smooth velocity and density profiles, obtained by a numerical solution of the TGE (4.18) as the volume flux  is increased (black curves).(a) Hyperbolic-tangent profiles for velocity U () and density R () profiles of increasing interface thickness.(b) Sinusoidal profiles for velocity U () (and same R () as in (a)).The growth is always slower than it would be using the sharp two-layer analytical solution (4.20) (red curves).