3.2.1. Analytical solution for the characteristics

The characteristics $\lambda$ of the three-layer model are obtained by solving for the eigenvalues in (3.5) that requires ${\det}(\boldsymbol{A}-\lambda \boldsymbol{C})=0$ for a non-trivial solution leading to the characteristic polynomial

(3.9) \begin{equation} \sum _{n=0}^{4} a_n \lambda ^n = 0, \end{equation}

with the coefficients

(3.10) \begin{align} a_4&=(h_0+\widehat {h}),\nonumber\\ a_3&=-2 \left [\widehat {h} \, u_0 + (u_1 + u_2) h_0+ \widehat {h}\, \widehat {u} \right ],\nonumber \\ a_2 &= \widehat {h} \, {u_0^2 } +4 \widehat {h}\, \widehat {u}\, u_0 + g_1 h_1 (h_0+h_2) \big(F_1^2-1\big) + g_2 h_2 (h_0+h_1) \big(F_2^2-1\big)\nonumber \\ & \quad +4\,h_0 \,u_1 \,u_2, \\ a_1&=-2 \widehat {h}\, \widehat {u} \,{u_0^2 } -2 g_1 h_1 h_2 \big(F_1^2-1\big)\,u_0-2 g_2 h_1 h_2 \big(F_2^2-1\big)\,u_0 -2\,h_0 \,{u_1 } \,u_2 (u_1+u_2)\nonumber \\ & \quad +2 \,h_0 \widehat {h}\, \widehat {u} \, \widehat {g}, \nonumber \\ a_0&={{\det}\; \unicode{x1D63C}}. \nonumber \end{align}

In these coefficients we introduced the height $\widehat {h}=h_1+h_2$ , velocity $\widehat {u}=(h_1 u_2 +h_2 u_1)/\widehat {h}$ and reduced gravity $\widehat {g}=(g_1 h_1 u_2 + g_2 h_2 u_1)/(\widehat {h} \,\widehat {u})$ . We also used the Froude numbers of the bottom, middle and top layers, respectively,

(3.11) \begin{eqnarray} F_2=\frac {u_2}{\sqrt { g_2 h_2}}, \quad F_0=\frac {u_0}{\sqrt {\frac {g_1 g_2}{g_1+g_2} h_0}}, \quad F_1=\frac {u_1}{\sqrt {g_1 h_1}}. \end{eqnarray}

The characteristics have the form

(3.12) \begin{equation} \begin{aligned} & \lambda _{1}=\,\bar {\!\lambda }-\delta _1\lambda - \delta _2 \lambda ,& \lambda _{2}= \,\bar {\!\lambda }-\delta _1\lambda + \delta _2 \lambda , \nonumber\\[-2pt] & \nonumber \lambda _{3}=\,\bar {\!\lambda }+\delta _1\lambda - \delta _3 \lambda ,& \lambda _{4}=\,\bar {\!\lambda }+\delta _1\lambda + \delta _3 \lambda ,\end{aligned} \end{equation}

where $\,\bar {\!\lambda }=-{a_3}/{4a_4}$ is the convective velocity and

(3.13) \begin{align} \begin{array}{c} \delta _1\lambda = \frac {1}{2}\sqrt {-\frac {2}{3} m_1 + \frac {1}{3 a_4} \left ( m_3 + \frac {\varDelta _0}{m_3} \right )}, \\ \\[-5pt] \delta _2 \lambda = \frac {1}{2} \sqrt {-4(\delta _1 \lambda )^2-2m_1+\frac {m_2}{\delta _1 \lambda }}, \quad \delta _3 \lambda = \frac {1}{2} \sqrt {-4\left (\delta _1 \lambda \right )^2-2m_1-\frac {m_2}{\delta _1 \lambda }} \end{array} \end{align}

are the phase speeds, where

(3.14) \begin{eqnarray} m_1&\!\!=\dfrac {8 a_4 a_2 - 3 a_3^2}{8 a_4^2},\quad m_2=\dfrac {a_3^2-4 a_4 a_3 a_2 + 8 a_4^2 a_1}{8 a_4^3},\quad m_3=\left ( \dfrac {\varDelta _1^2+\sqrt {\varDelta _1^2-4\varDelta _0^3}}{2} \right )^{\tfrac {1}{3}}\!, \quad \nonumber\\[-1.5pt] &&\\[-2.5pt]\nonumber \varDelta _0 &\!\!= a_2^2 - 3 a_3 a_1 + 12 a_4 a_0, \quad \varDelta _1 = 2 a_2^3 - 9 a_3 a_2 a_1 + 27 a_3^2 a_0 + 27 a_4 a_1^2 -72 a_4 a_2 a_0. \end{eqnarray}

Since $(\delta _1 \lambda ^R,\delta _2 \lambda ^R,\delta _3 \lambda ^R) \geqslant 0$ , we realise that $\lambda _1^R \leqslant \lambda _2^R$ , $\lambda _3^R \leqslant \lambda _4^R$ and $\lambda _1^R \leqslant \lambda _4^R$ , where the superscript $R$ denotes the real component. Consequently, in a reference frame moving at $\,\bar {\!\lambda }$ , an observer would see two oppositely propagating interfacial waves, which are associated with the eigenvalues $\lambda _{1,2}$ and $\lambda _{3,4}$ .

To justify this association, we note that the terms $\delta \lambda = -\delta _1\lambda \mp \delta _2 \lambda$ and $\delta \lambda = \delta _1\lambda \mp \delta _3 \lambda$ represent the effective phase speeds of interfacial waves in a three-layer flow. The structure of these expressions suggests that $\lambda _{1,2}$ and $\lambda _{3,4}$ correspond to long waves propagating along the upper and lower interfaces, respectively. Moreover, since the coefficients in (3.10) are real, if any characteristics become imaginary, they must appear in complex conjugate pairs. Complex conjugate eigenvalues correspond to physically linked interfacial waves, reinforcing the idea that each pair $(\lambda _1, \lambda _2)$ and $(\lambda _3, \lambda _4)$ can be naturally associated with a single interface. Thus, a key observation from (3.12) is that the eigenvalues $\lambda _{1,2}$ and $\lambda _{3,4}$ can be identified with long waves at the upper and lower interfaces, respectively. We note this labelling does not imply that the interfaces are completely decoupled, as their characteristics inherently depend on the velocities and heights of all three layers.

3.2.2. A special case: pure exchange flow with stagnant middle layer

Of particular interest is identifying the conditions under which interfacial waves become unstable. To make progress on this question, we examine a special asymmetric case where $u_2 = -u_1$ , $h_1 = h_2$ and $g_1 = g_2$ . The zero (barotropic) net flow condition implies that $u_0=-(u_1 h_1 + u_2 h_2)/h_0=0$ and the coefficient $a_3$ vanishes, leading to zero convective velocity ( $\,\bar {\!\lambda }=0$ ). This corresponds to the simplest three-layer exchange flow, where the upper and lower layers are symmetric in thickness and stratification, while the middle layer remains stationary. This case is also highly relevant to the SID flow, particularly near the centre of the duct. Under these conditions, the expressions for characteristics in (3.12) simplify significantly (since $F_1 = F_2)$ , yielding

(3.15) \begin{equation} \begin{aligned} \lambda _{1,2} &= \mp \frac {\sqrt {g_1 h_1 }}{\sqrt {h_0 +2\,h_1 }}\,\sqrt {h_0\big(1+{F_1^2 }\big)+h_1\big(1-{F_1^2 }\big) - \sigma}, \\ \lambda _{3,4} &= \mp \frac {\sqrt {g_1 h_1 }}{\sqrt {h_0 +2\,h_1 }}\sqrt {h_0\big(1+{F_1^2 }\big)+h_1\big(1-{F_1^2 }\big) +\sigma },\end{aligned} \end{equation}

where $\sigma =\sqrt {4 h_0 h_1 {F_1^2 }(1 - {F_1^2 } ) + 4 {h_0^2 } {F_1^2 } +{h_1^2 } (1-{F_1^2 } )^2}$ . Clearly, in this case, if we further restrict the system to $F_1^2=1$ then $\lambda _{1,2}=0$ and $\lambda _{3,4}=\mp {2 \sqrt {g_1 h_0 h_1}}/{\sqrt {h_0 +2\,h_1}}$ . Since waves in the upper and lower layers become decoupled in this limit, this result supports our approach of associating each eigenvalue pair $\lambda _{1,2}$ and $\lambda _{3,4}$ with a specific interface. However, we note that interfacial waves are coupled where $\,\bar {\!\lambda } \ne 0$ in general.

Figure 5.

The onset of long-wave instability in the pure exchange flow case visible where (a) $\lambda _1^I\gt 0$ and (b) $\lambda _4^I\gt 0$ . Here, $F$ denotes the Froude number of the upper layer and $r$ represents the ratio of middle-to-upper/lower layer heights.

Equation (3.15) reveals that the onset of complex characteristics, indicative of long-wave instability, can be anticipated by identifying the condition under which $\sigma$ vanishes. Setting $\sigma = 0$ yields the critical Froude number thresholds

(3.16) \begin{equation} F_{\varDelta \mp }^2 = 1 + \frac {2r(r - 1) \mp 2r^{3/2} \sqrt {r + 2}}{4r - 1}, \end{equation}

where $r = h_0/h_1$ denotes the ratio of the middle layer depth to the upper (or lower) layer depth. Figure 5 shows $\lambda _1^I$ and $\lambda _4^I$ as functions of $F_1^2$ for the pure exchange flow. The critical Froude numbers predicted by (3.16) are shown as red dashed lines. We observe that, for $F_1^2 \geqslant F_{\varDelta +}^2$ , the characteristics become complex, and the largest values of $\lambda _4^I$ occur within the range $F_{\varDelta +}^2 \leqslant F_1^2 \leqslant F_{\varDelta -}^2$ . Moreover, the magnitudes of $\lambda _1^I$ and $\lambda _4^I$ change abruptly across the $F_{\varDelta -}^2$ curve. We identify $F_{\varDelta -}^2$ as the bifurcation Froude number and $F_{\varDelta +}^2$ as the marginal-stability Froude number for the pure exchange flow case.

The velocity difference between the upper and middle layers is given by $u_1^2 = F_1^2 g_1 h_1$ , since $u_0 = 0$ in the pure exchange flow. Hence, for $F_1^2 \geqslant F_{\varDelta +}^2$ , the velocity shear between the upper and middle layers, and between the middle and lower layers, becomes sufficiently strong to trigger long-wave instability on each interface separately. As $r$ increases, the middle layer thickens relative to the upper layer, causing the upper layer to thin and accelerate. This leads to larger velocity differences across the interfaces, particularly between the faster upper layer and the stagnant middle layer, leading to instability once $F_1^2 = F_{\varDelta +}^2$ . This behaviour is analogous to that observed in two-layer flows, where sufficiently strong interfacial shear also leads to long-wave instability.

The condition $\sigma = 0$ corresponds to $\lambda _1 = \lambda _3$ and $\lambda _2 = \lambda _4$ , which implies that the interfacial waves on the lower and upper interfaces have equal phase speeds and can resonate. Instability may therefore arise from this resonance. In the more general case with non-zero convective velocity, however, the identification of instability is less straightforward, a point revisited below and addressed in detail in § 4.

3.2.3. Criticality of interfaces

In a two-layer configuration, the sign of the real component of the characteristics sets the direction of information propagation and the hydraulic regime of the flow, i.e. subcritical, critical or supercritical regime (Long Reference Long1956; Dalziel Reference Dalziel1991; Atoufi et al. Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023). The flow is subcritical when $\,\bar {\!\lambda }\lt \delta \lambda ^R$ and information propagates in both leftward (towards decreasing $x$ ) and rightward (towards increasing $x$ ) directions, supercritical when $\,\bar {\!\lambda }\gt \delta \lambda ^R$ and information propagates only either rightward or leftward and critical when $\,\bar {\!\lambda }=\delta \lambda ^R$ .

In a three-layer configuration, complications arise as the upper and lower layers can individually exhibit different hydraulic regimes. As we found above, $\lambda _{1,2}$ and $\lambda _{3,4}$ are respectively associated with the upper and lower interface. Therefore, the necessary condition for the upper (respectively lower) layer to be supercritical is $\lambda _1^R \lambda _2^R \gt 0$ (respectively $\lambda _3^R \lambda _4^R \gt 0$ ). Information on each interface thus propagates in a particular direction, depending on the signs of the characteristics (Sannino et al. Reference Sannino, Carillo and Artale2007, Reference Sannino, Pratt and Carillo2009).

The characteristic pairs $\lambda _1^R \lambda _2^R$ and $\lambda _3^R \lambda _4^R$ determine criticality with respect to the individual layers, not the full three-layer system. For example, the flow may be supercritical in both layers if $\lambda _1^R \lambda _2^R \gt 0$ and $\lambda _3^R \lambda _4^R \gt 0$ , yet remain subcritical overall, e.g. if $\lambda _1^R, \lambda _2^R \lt 0$ and $\lambda _3^R, \lambda _4^R \gt 0$ . The three-layer flow is fully supercritical only when all layers are hydraulically supercritical and the corresponding characteristics have the same sign, so that information propagates in a single direction. It is fully subcritical when all layers are subcritical, or when the layers are supercritical but the characteristic signs differ, indicating opposing directions of information propagation. A critical condition arises when at least one characteristic pair satisfies $\lambda _i^R \lambda _j^R = 0$ , at which point information propagation is blocked and internal hydraulic control is established.

3.2.4. Critical middle layer thickness

The emergence of imaginary components, as conjugate pairs, in the characteristics implies the onset of instability of long waves as established in Appendix A. Consequently, we can predict critical values of the layer velocities and heights that cause interfaces to be unstable to disturbances of the given wavenumber $k$ by examining the coefficients $b_i$ of the quartic stability equation (A10) given below,

(3.17) \begin{align} b_4 &= \sinh \left (k\,{\left (h_0 +\widehat {h} {\kern2pt}\right )}\right ), \nonumber\\ b_3&=\sinh \left (k\,(h_0-\widetilde {h})\right ){\left (u_0 -u_1 \right )}- b_4 (u_1+u_2+2 u_0) + \sinh \left (k (h_0+\widetilde {h})\right ){\left (u_0 -u_2 \right )}, \nonumber\\ b_2 &= u_0^2 (\alpha _1 + \alpha _2 + 6 \alpha _4) + u_1^2 (\alpha _1 + \alpha _3) + u_2^2 (\alpha _2 + \alpha _3) + 4 u_0 (u_1 \alpha _1 + u_2 \alpha _2) + 4 u_1 u_2 \alpha _3 \nonumber \\ & \quad + \gamma _1 + \gamma _2, \\ b_1 &= \big(1 - 2 \sigma _1^2\big)(\gamma _{7} + \gamma _9) - \sigma _1 \sigma _2 \gamma _{8}, \nonumber\\ b_0 &= u_0^2 \gamma _3 + 2 u_0^4 \sigma _1 \beta _1 \sigma _2 + 2 u_1^2 u_2^2 \beta _4 \sigma _1 \sigma _2 + u_2^2 \gamma _4 + u_1^2 \gamma _5 + \gamma _6,\nonumber \end{align}

and to long waves by examining $a_i$ in (3.10) where $\widetilde {h}={h_1-h_2}$ and $\alpha _i, \sigma _i$ and $\gamma _i$ are defined in Appendix A.1. For example, the necessary condition guaranteeing only real roots is that $b_i$ should be all positive or all negative (based on the Routh–Hurwitz stability criterion (Shinners Reference Shinners1998). As wavenumbers $k \geqslant 0$ , then $b_4 \geqslant 0$ and this condition leads to one of the $b_3, b_2, b_1, b_0 \lt 0$ for all $k$ for instability to happen.

In fact, the number of sign changes in $b_i$ yields the number of complex roots and potential instability. As $\lambda$ (i.e. the long-wave limit of the phase speed $c$ of disturbance waves) appears in the form of complex conjugates, there are four complex roots and $a_{3,2,1,0}$ change sign (note $a_4$ is constant inside the duct). In the limit of long waves, using the $b_3=a_3 \lt 0$ condition and incorporating the zero net flow condition $Q=u_0 h_0 + u_1 h_1 + u_2 h_2 = 0$ to eliminate $u_0$ in $a_3$ yields a restriction for the height of the middle layer that $ -{2{(h_0 +\widehat {h})}{ (h_0 (u_1 + u_2) -h_1 u_1 -h_2 u_2 )}}/{h_0 } \lt 0$ for instability to occur.

This leads to

(3.18) \begin{equation} h_0^{cr}=\frac {u_1 h_1 + u_2 h_2}{u_1 + u_2} \end{equation}

as the critical height for the middle layer. Therefore, when $u_2 \lt {\lvert {u_1} \rvert }$ (i.e. an interface with a positive slope) and $h_0 \gt h_0^{cr}$ , the three-layer exchange flow is become unstable to long waves and may be supercritical. In the limit of short waves, where $k \gg 1$ , the sign of $b_3$ depends on the coefficients of the largest $k$ (positive) when expanding $b_3$ with respect to $k$ . Specifically, when $u_2 h_0 + u_1 h_1 + u_2 h_2 \gt 0$ , we have $b_3 \lt 0$ , satisfying the necessary condition. This leads to the condition $h_0 \gt {(|u_1| h_1 - u_2 h_2)}/{u_2}$ , which provides a smaller lower bound for $h_0$ and a smaller critical middle layer thickness than in (3.18). Therefore, it is expected that short waves will become unstable first. Nonetheless, it is the long waves (hydraulic effects) that set the interface slopes, and thus, the upper and lower layer flow rates $u_1 h_1$ and $u_2 h_2$ , which then determine the value of $h_0$ that results in $b_3 \lt 0$ and the onset of the short-wave instabilities.

3.2.5. Application to the data

The three-layer model, yielding layer characteristics, is applied to the layer-averaged DNS data in figure 6. The space–time evolution of $\lambda _1$ and $\lambda _3$ in the TW case are shown with $x{-}t$ plots. The dashed lines are reproduced from the layer-averaged quantities in figure 3, indicating the drastic changes in the layer statistics. Noticeably, the left and right lines are, respectively, in alignment with $\lambda _1^R$ and $\lambda _3^R\approx 0$ , which indicates the emergence of critical points in the upper and lower layer interfaces. The flow becomes critical on these lines, which is also evidenced by the appearance of hydraulic jumps as illustrated by the sudden changes of the corresponding layer heights in figure 3. In addition, the three-layer model predicts one control point near the inlet/outlet of the duct due to geometry changes, which agrees with the two-layer model in Atoufi et al. (Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023).

Figure 6.

Selected eigenvalues of the three-layer model for the TW case ( $\textit{Re} = 650, \theta = 6^\circ$ ): (a,d) real and imaginary components of the characteristic $\lambda _1$ for the upper interface, (b,e) $\lambda _3$ for the lower interface and (c,f) product of the real and imaginary components of $\lambda _1$ and $\lambda _3$ . The dashed and dash-dotted lines correspond to those in figure 3.

In figure 6(c), outside of the critical lines, the real components of the characteristics have the same sign and are positive/negative on the left/right side of the duct, respectively. Therefore, the information of the three-layer system can only propagate from the critical point near the inlet/outlet toward the duct centre. In between the critical lines, $\lambda _1^R\lambda _3^R\lt 0$ and the flow is subcritical. Waves in this region can propagate in either direction. This region is where the thickening of the middle layer occurs (see figure 3).

The imaginary components of characteristics correspond to the growth rate of the long interfacial waves (as discussed in Appendix A). As indicated by figure 6(d,e), the long waves at the upper interface become unstable ( $\lambda _1^I\gt 0$ ) in the left and middle regions of the duct, whereas the long waves at the lower interface are unstable in the middle and right regions ( $\lambda _3^I \gt 0$ ) where the superscript I refers to the imaginary component of $\lambda$ . While the stability of individual interfaces can be deduced from the associated characteristics, the stability of the three-layer flow system is less straightforward and requires consideration of resonance between interfaces, which will be discussed in § 4.

Figure 7.

Instantaneous characteristics of the interfacial waves on the upper and lower interfaces in DNS, diagnosed using the three-layer hydraulic analysis. The real components of the characteristics, representing the phase speed of the upper and lower interfacial waves ( $\lambda _1^R$ and $\lambda _3^R$ ), are shown in (a,c). The imaginary components, representing the growth rate of the upper and lower interfacial waves ( $\lambda _1^I$ and $\lambda _3^I$ ), are shown in (b,d).

The real and imaginary components of the instantaneous characteristics are shown in figure 7 at $t = 120$ for the L, SW, TW and I cases, to complement the space–time distribution of the characteristics, which was previously given only for the TW case. The instantaneous $\lambda ^R$ and $\lambda ^I$ of all non-laminar cases follow a similar trend. The imaginary component in figure 7(b,d) shows that in the L case, $\lambda ^I \approx 0$ for both the upper and lower interfaces, whereas in the SW, TW and I cases, $\lambda _1^I, \lambda _3^I \gt 0$ in most places inside the duct.

The application of the three-layer theory to the DNS data suggests that three-layer characteristics are not only useful in quantifying information propagation but also in informing us about regions where turbulence and mixing occur. We detail this connection between characteristics and turbulence next.

3.3.1. Critical exchange flow regime

To define hydraulic control, we use the steady-state solution of (3.1). For $\boldsymbol{\mathcal{F}}=0$ , the steady solution of (3.1) can be found by multiplying the inverse matrix $\unicode{x1D63C}^{-1}={\rm adj} (\unicode{x1D63C})/{\det}\; \unicode{x1D63C}$ on both sides leading to

(3.19) \begin{equation} \frac {\partial \boldsymbol{q}}{\partial x}=\frac {1}{{\det}\; \unicode{x1D63C}} {\rm adj} (\unicode{x1D63C}) \unicode{x1D63D}\frac {\partial \boldsymbol{s}}{\partial x}+\frac {1}{{\det}\; \unicode{x1D63C}} {\rm adj} (\unicode{x1D63C}) {\tau } \boldsymbol{q}, \end{equation}

where ${\rm adj} (\unicode{x1D63C})$ represent the transpose of the co-factor matrix of $\unicode{x1D63C}$ . A non-trivial steady solution $\boldsymbol{q}(x)$ of (3.1) requires

(3.20) \begin{equation} {\det}\; \unicode{x1D63C}= g_1 g_2 h_1 h_0 h_2 ({G}-1) \ne 0, \end{equation}

where $G$ is the composite Froude number defined based on the individual layers

(3.21) \begin{equation} {G}=1+{\big (\epsilon _1 {F_1^2 } + \epsilon _2 {F_2^2 } -1\big )}\,{F_0^2 } +{\big ({F_1^2 } -1\big )}\,{\big ({F_2^2 } -1\big )}, \end{equation}

where $\epsilon _1=g_1/(g_1+g_2)$ , $\epsilon _2=g_2/(g_1+g_2)$ are the ratio of cross-interface to total density differences.

The three-layer flow is hydraulically controlled when $h_0=0$ or ${G}=1$ (i.e. ${\det}\; \unicode{x1D63C}=0$ ). The condition $h_0=0$ is trivial as it removes the third layer. Generally, the hydraulic control condition incorporating (3.21) defines a critical state in the three-layer exchange flow, where

(3.22) \begin{equation} F_0^2=\frac {\left(F_1^2-1\right)\left(F_2^2-1\right)}{\left(1- \epsilon _1 F_1^2 - \epsilon _2 F_2^2\right)}, \end{equation}

yielding $G=1$ . This equation imposes a fundamental constraint on the middle layer, determining its velocity and thickness based on the Froude numbers of the upper and lower layers, thereby fully characterising its dynamics under hydraulic control. Determining whether a three-layer exchange flow beyond the critical state is subcritical or supercritical requires linking (3.21) to the characteristics of long interfacial waves, as discussed in the next section. Nonetheless, (3.22) allows the exploration of the possible combinations of supercritical and subcritical flows in each layer even though the overall three-layer system remains in a critical state. Figure 8(a) shows the middle layer Froude number $F_0^2$ at the critical state, computed from (3.21) for $\epsilon _1 = \epsilon _2 = 0.5$ , as a function of $F_1^2$ and $F_2^2$ . The structure is similar to figure 3 in Sannino, Pratt & Carillo (Reference Sannino, Pratt and Carillo2009), based on a three-layer model for the Strait of Gibraltar. Figure 8(b) shows the corresponding contour plot in the $F_1^2$ – $F_2^2$ plane. Following Sannino et al. (Reference Sannino, Pratt and Carillo2009), Sánchez-Garrido et al. (Reference Sánchez-Garrido, Sannino, Liberti, García Lafuente and Pratt2011), four distinct regimes are identified and labelled using the notation $(F_1^2, F_2^2, F_0^2)$ , where each symbol ( $\gt$ or $\lt$ ) denotes whether the corresponding layer is hydraulically supercritical or subcritical relative to the others. These labels illustrate how individual layers may exhibit distinct hydraulic states, even though the overall flow remains critical. For instance, $(\gt ,\lt ,\gt )$ represents a critical regime where the upper and middle layers are supercritical, while the lower layer is subcritical. Sannino et al. (Reference Sannino, Pratt and Carillo2009) referred to the region marked $(\lt ,\lt ,\lt )$ as a provisionally subcritical regime, in which all layers are subcritical and wave motions in the upper and lower layers are coupled. They referred to the region $(\gt ,\gt ,\lt )$ as a provisionally supercritical regime, where the upper and lower layers are supercritical while the middle layer remains subcritical.

Figure 8.

Critical state of the exchange flow. (a) Surface plot of the middle layer Froude number $F_0^2$ as a function of the upper layer ( $F_1^2$ ) and lower layer ( $F_2^2$ ) Froude numbers, satisfying the critical condition given in (3.22). (b) Contour plot of the same critical surface, showing level curves of $F_0^2$ in the $F_1^2$ – $F_2^2$ plane. Regions are labelled using triple inequality notation (e.g. $\gt ,\lt ,\lt$ ), indicating whether each layer is supercritical or subcritical relative to the others. These combinations illustrate how individual layers may depart from criticality, while the overall three-layer flow remains in a critical state.

In the special case where the middle layer is stagnant ( $u_0 = 0$ , so that $F_0 = 0$ ) and either the upper layer or lower layer Froude number reaches unity, i.e. ${\lvert {F_1} \rvert } = 1$ or ${\lvert {F_2} \rvert } = 1$ , a hydraulically controlled region emerges – bounded by ${\lvert {F_1} \rvert } = 1$ or ${\lvert {F_2}\rvert } = 1$ – in contrast to the localised control points typically found in two-layer flows (Atoufi et al. Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023).

To extend this classification beyond the critical state, the characteristics must be related to a composite Froude number $G$ , as introduced in the following section, that reflects the direction of information propagation.

3.3.2. Three-layer hydraulic regimes identifier

The link between characteristics and $G$ can be highlighted by noting $\Pi _{i=1}^4 \lambda _i = {a_0}/{a_4}$ combined with (3.20), yielding the identity

(3.23) \begin{eqnarray} {G} = 1 + \dfrac {(h_1 + h_0 + h_2)}{g_1 g_2 h_1 h_0 h_2} \Pi _{i=1}^4 \lambda _i. \end{eqnarray}

Although the condition $G=1$ in a three-layer flow can be used to identify control points, neither (3.21) nor (3.23) can provide insights into whether the three-layer flow is subcritical or supercritical. The reason is that when characteristic speeds of long interfacial waves at both interfaces have imaginary components, the multiplication of all complex characteristic speeds mixes up information propagation associated with each interface and typically gives $G \geqslant 1$ . As only the real components of the characteristics determine the hydraulic regimes (imaginary components quantify the long-wave growth rate, as shown in Appendix A), based on the definition of $G$ , we define a modified composite Froude number

(3.24) \begin{eqnarray} \widetilde {G}=1+ \dfrac {(h_1+h_0+h_2)}{g_1 g_2 h_1 h_0 h_2} \textrm{sign}{\big(\lambda ^R_1 \lambda ^R_4\big )} \Pi _{i=1}^4 \lambda ^R_i \end{eqnarray}

that is capable of identifying hydraulic regimes in three-layer flows with complex characteristics. Note that $\widetilde {G}$ can be negative, and $\widetilde {G}=1$ at control points, but generally $G \geqslant \widetilde {G}$ . When the real components of the four characteristics have the same sign, then $\widetilde {G}\gt 1$ , indicating a supercritical regime. Conversely, when $\widetilde {G}\lt 1$ , the real components have opposite signs, indicating a subcritical regime. The condition $\widetilde {G}=1$ indicates critical flow since the real component of the characteristic associated with one of the interfaces vanishes.

Figure 9.

(a–c) Normalised TKE $k^\prime _m/\langle k^\prime _m\rangle _{\mathcal{V},t}$ averaged over the duct cross-section ( $\langle \cdot \rangle _{\mathcal{V},t}$ denotes the spatial and temporal average). (d–f) Modified composite Froude number for the three-layer model, $\widetilde {G}$ , as defined in (3.24). (g–i) Modified composite Froude number for the two-layer model, $\widetilde {G}^{2L}$ , as defined in (3.28). The first row represents the SW case, the second row represents the TW case and the third row represents the I case. The dashed and dashed-dotted green curves are identical to those in figure 3.

Thus, $\widetilde {G}$ provides a diagnostic based on the direction of information propagation by long interfacial waves, rather than a strict indicator of a hydraulic regime under instability. The modified composite Froude number $\widetilde {G}$ is specifically defined to involve only the real parts of the characteristic speeds, avoiding the complications introduced by long-wave instability. As discussed in Appendix D, the real part $\lambda _i^R$ continues to predict the direction of information propagation by long interfacial waves to leading order, provided the spatiotemporal variation in $\lambda$ is weak. Therefore, while $\widetilde {G}$ should not be interpreted as a strict indicator of a hydraulic regime in the classical sense under instability, it remains a useful diagnostic for identifying regions where long interfacial wave propagation becomes predominantly unidirectional (i.e. supercritical regime), bidirectional (i.e. subcritical regime) or blocked (i.e. critical regime), even when the three-layer flow contains unstable long-wave modes.

3.3.3. Link to turbulence generation

Having established $\widetilde {G}$ as an indicator of the hydraulic regime, we now aim to show the link between transition in hydraulic regimes and the turbulence generation and compare predictions from the three-layer model to those of the two-layer model of Atoufi et al. (Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023). Following Zhu et al. (Reference Zhu, Atoufi, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023), we diagnose the local turbulent kinetic energy (TKE)

(3.25) \begin{eqnarray} k^\prime _{m}\equiv \frac {1}{2}\boldsymbol{u}^\prime _m \cdot \boldsymbol{u}^\prime _m, \end{eqnarray}

computed from the perturbations around a low-pass filtered velocity field

(3.26) \begin{align} \boldsymbol{\bar {u}}_m\!\left (x,y,z,t\right )&\equiv \frac {1}{\Delta L}\int _{-\Delta L/2}^{\Delta L/2}{\boldsymbol{u}\!\left (x-s,y,z,t\right ){\rm d}s},\end{align}

(3.27) \begin{align} \boldsymbol{u}^\prime _m(x,y,z,t)&\equiv \boldsymbol{u}-\boldsymbol{\bar {u}_m}.\end{align}

We choose the same constant filter width $\Delta L=10$ to maximise the time- and duct-volume-averaged kinetic energy. Figure 9(a,b,c) shows the normalised TKE averaged over the duct cross-section and figure 9(d,e,f) shows $\widetilde {G}$ . For comparison, we also show the modified composite Froude number for the two-layer hydraulic model,

(3.28) \begin{equation} \widetilde {G}^{2L}=1+\frac {h_1^{2L}+h_2^{2L}}{\big({\rho }_2^{2L}-{\rho }_1^{2L}\big) {Ri} \cos {\theta } \ h_1^{2L} h_2^{2L}} \ \lambda _1^{R^{2L}} \ \lambda _2^{R^{2L}} \end{equation}

in figure 9(g,h,i), where the superscripts ‘2L’ denote variables based on the two-layer averaging of Atoufi et al. (Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023) where the expressions for $\lambda _1$ and $\lambda _2$ are also given in their (4.10).

As expected, $\widetilde {G}$ in figure 9(e) is consistent with figure 6(a,b) for the TW in the sense that it indicates control points as well as the hydraulic regime. In all three cases, the flow experiences a transition from supercritical to subcritical in the central region of the duct where $\widetilde {G}$ changes sign. The same auxiliary dashed guidelines (i.e. space–time variation of middle layer height) presented in figures 3 and 6 agree well with $\widetilde {G} \approx 1$ . We observe a similar agreement between $\widetilde {G}$ and the characteristics in the SW and I cases, although their $x{-}t$ plots are not shown for brevity. Therefore, $\widetilde {G}$ identifies hydraulic regimes and control (i.e. flow regularisation when $G=\tilde {G}=1$ discussed in § 3.4) even in the presence of mixing and beyond the limit where (3.1) is no longer hyperbolic. Although not shown here, the contours of $G$ have qualitatively similar patterns to $\widetilde {G}$ except that $G \geqslant 1$ in the entire $x{-}t$ plots for the SW, TW and I cases due to $\lambda ^I \ne 0$ . Reduced values of $G$ are also observed at the same locations and times where $\widetilde {G} \lt 0$ in figure 9(d,e,f).

In the centre of the duct between the auxiliary curved lines (corresponding to the upper and lower interface locations) in figure 9, the averaged TKE tends to be elevated in the TW and I cases (figure 9 b,c), suggesting a correlation between enhanced turbulence and hydraulic regime transitions identified by $\widetilde {G} = 1$ (figure 9 e,f). In these cases, large values of $\langle k_m^{\prime } \rangle _{y,z}$ coincide with regions where $\widetilde {G}$ transitions from $\widetilde {G} \gt 1$ (supercritical) to $\widetilde {G} \lt 1$ (subcritical), particularly near the duct exits and within central interior regions. In the SW case (figure 9 a) this trend is less apparent where the TKE appears visually weak across most of the space–time domain. This reflects a lower level of turbulence rather than a lack of correlation. Nevertheless, localised turbulence is still present, especially near the duct ends and at the central region, which align with changes in $\widetilde {G}$ (figure 9 d). These observations indicate that, although the spatiotemporal correspondence between TKE and hydraulic regime transitions is strongest in the TW case and more distributed in the I case, it remains qualitatively evident across all three cases. In the I case, small-scale shear instabilities – such as Kelvin–Helmholtz roll-up and breakdown – embedded within the non-hydrostatic pressure gradient force $\boldsymbol{\mathcal{F}}$ also contribute to turbulence generation, in addition to hydraulic effects that are primarily governed by the hydrostatic pressure distribution. As shown in Atoufi et al. (Reference Atoufi, Zhu, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023), the influence of non-hydrostatic pressure is more pronounced in the I case than in the SW and TW cases.

Near both ends, we typically observe short waves forming and breaking leading to turbulence (Zhu et al. Reference Zhu, Atoufi, Lefauve, Taylor, Kerswell, Dalziel, Lawrence and Linden2023), suggesting that the flow is unstable to both long- and short-wave instabilities there. As we approach the centre of the duct, the contribution of long waves to the destabilisation of the flow becomes increasingly significant.

For the two-layer model, $\widetilde {G}^{2L} = 1$ is more weakly correlated with $\langle k^\prime _m\rangle _{y,z}$ . A key result of this paper is that the three-layer $\widetilde {G}$ outperforms the two-layer $\widetilde {G}^{2L}$ in predicting regions (not only points) where the flow is most unstable and turbulent. Therefore, regime transitions identified by changes in sign of $\widetilde {G}$ (i.e. changes in the direction of information propagation) provide excellent predictions of where and when turbulence is generated and strong mixing is expected in stratified exchange flows. The mechanism underlying this agreement is explained in § 4.

3.4.1. Inviscid control mechanism

The inviscid regularity condition in (3.29) implies that

(3.30) \begin{align} \frac {1}{3} {\rm adj} (\unicode{x1D63C}) \unicode{x1D63D}\frac {\partial \boldsymbol{s}}{\partial x}= \left (\begin{array}{c} g_1 u_1 h_2 \left ( {u_0^2 } - g_2 h_0 \left (1-F_2^2 \right )\right ) \left (\tan {\theta }-\dfrac {\text{d} b}{\text{d} x} \right ) -g_2 h_2 {u_0^2 } u_1 \\[10pt] g_1 g_2 h_1 h_2 u_0 \left (1-F_1^2\right ) +\left (g_1 h_1\right )^2 u_0 \left (1-F_1^2\right ) \left (\tan {\theta } -\dfrac {\text{d} b}{\text{d} x} \right )\\[12pt] g_2 u_2 \left (h_1 {u_0^2 } - g_1 h_0 h_1 \left (1-F_1^2\right )\right ) -g_1 h_1 {u_0^2 } u_2 \left (\tan {\theta }-\dfrac {\text{d} b}{\text{d} x}\right ) \\[12pt] g_2 h_1 h_2 {u_0^2 } -g_1 h_1 \left (h_2 {u_0^2 } - g_2 h_0 h_2 \left (1-F_2^2 \right )\right ) \left (\tan {\theta } -\dfrac {\text{d} b}{\text{d} x} \right ) \\[12pt] g_1 g_2 h_0 h_1 h_2 \left (1-F_2^2 \right ) \left (\dfrac {\text{d} b}{\text{d} x}-\tan {\theta } \right ) -g_1 g_2 h_0 h_1 h_2 \left (1-F_1^2\right ) \\[12pt] g_1 h_1 h_2 {u_0^2 } \left (\tan {\theta }-\dfrac {\text{d} b}{\text{d} x} \right )-g_2 h_2 \left (h_1 {u_0^2 } - g_1 h_0 h_1 \left (1-F_1^2\right )\right ) \end{array}\right )=\boldsymbol{0} \end{align}

at the control points. We recall that two conditions lead to ${\det}\; \unicode{x1D63C}=0$ and hydraulic control: either (i) $h_0=0$ or (ii) $h_0 \neq 0$ but $G=1$ , which further requires ${\lvert {F_1 } \rvert }=1 \ \text{or} \ {F_2 }=1$ and $F_0=0$ (i.e. $u_0=0$ ).

For condition (i), (3.30) enforces $u_0=0$ . Therefore, both the criticality and inviscid control mechanisms imply that the middle layer velocity vanishes near the jumps. This is consistent with the sudden reduction of $\lvert {u_0} \rvert$ near the duct ends and close to the control point in figure 4(d). Note that near the duct ends where ${\lvert {{\text{d} b}/{\text{d} x}} \rvert } \gg \tan {\theta }$ the flow is also hydraulically controlled as $h_0 \approx 0$ (condition (i)) and $\lvert {u_0} \rvert$ drops sharply in figure 4(a) to satisfy (3.30). Inside the duct ${\text{d} b}/{\text{d} x}=0$ and $\lvert {u_0} \rvert$ drops more smoothly.

As stated above, one possibility for $G=1$ is $F_0=0$ and either ${\lvert {F_1} \rvert }=1$ or ${\lvert {F_2} \rvert }=1$ at control points. For condition (ii), (3.30) enforces that ${\lvert {F_1} \rvert }=1$ (critical upper layer) yields ${\lvert {F_2} \rvert }=1$ (critical lower layer), assuming the middle layer has non-zero velocity. However, this implication does not hold in general. For example, when $u_0 = 0$ , $F_1 = 0$ and ${\text{d} b}/{\text{d} x} = \tan \theta$ , all terms in (3.30) vanish, the condition is trivially satisfied and does not constrain $F_2$ , meaning that both interfaces need not be simultaneously critical. This reflects the fact that when the middle layer is stationary near the control point, the system may exhibit decoupled dynamics at the upper and lower interfaces. Criticality may occur in one interface without enforcing it in the other as also discussed for the special case of pure exchange flow in § 3.2.2. This type of decoupling, particularly relevant for $ N$ -layer ( $N\gt 2$ ) flows with stagnant interior layers, was analysed by Engqvist (Reference Engqvist1996), who showed that under such conditions, the upper and lower groups of non-stagnant layers can be governed by independent control criteria.

In summary, in addition to flow regularisation, condition (i) implies a thin middle layer and a two-layer exchange flow at both ends. Under condition (ii), when the middle layer is dynamically active, the two oppositely propagating flows at critical velocity ( ${\lvert {F_1} \rvert } = {\lvert {F_2} \rvert } = 1$ ) cause vigorous mixing ( $h_0 \ne 0$ ) in the middle of the duct (due to symmetry). However, when the middle layer is stationary near the control, the upper and lower interfaces may become decoupled, and criticality may occur independently at one or both interfaces without requiring symmetry.

3.4.2. Viscous control mechanism

The effects of viscosity are of particular importance as they allow for enhanced dissipation in hydraulic jumps, and are connected to turbulent mixing when the hydraulic regime transitions from supercritical to subcritical.

In figure 10 we computed the viscous control term ${\tau } \boldsymbol{q}$ (non-zero elements) to verify that it indeed dominates the right-hand side of (3.1) in the central region of the duct where $\widetilde {G} = 1$ (at a representative time $t=110$ ). Panels (a,b) present the viscous terms and the sum of all the terms on the right of (3.1), respectively, and show that the viscous term dominates, particularly in the central region of the duct. For the wave cases, the viscous control mechanism, represented by $\boldsymbol{\tau q}$ , becomes $\lesssim 0.01$ in the centre of the duct to regulate the flow (§ 3.4) at the internal control points. The summation of these terms also reaches a local minimum near the control points. Note that the viscous term and sum also become approximately zero at the duct ends, where geometric control is established. The regularity condition is therefore satisfied for all control mechanisms at the duct ends.

The large contribution of viscous terms ${\tau } \boldsymbol{q}$ suggests that information propagation along the interfaces is predominantly modulated by viscous effects (see § 3.2). By comparing figures 7 and 10, we observe that along the upper interface, viscous damping occurs on the left side of the duct, where $\zeta _1 = (\tau q)_1 \gt 0$ (with the subscript 1 indicating the term from the velocity difference equation between the middle and upper layers, i.e. the first row of (3.1)) and $\lambda _1^R \gt 0$ , indicating viscous dissipation. On the right side of the duct, where $\zeta _1 = (\tau q)_1 \lt 0$ and $\lambda _1^R \lt 0$ , the viscous forcing reinforces the direction of wave propagation. This results in a local amplification of interfacial waves due to enhanced mixing with the lower layer. Similarly, along the lower interface, viscous damping is observed on the right side, where $\zeta _2 = (\tau q)_2 \lt 0$ and $\lambda _3^R \gt 0$ , again indicating dissipation. Here, the subscript 2 refers to the velocity difference equation between the lower and middle layers (the second row of (3.1)). On the left side, where $(\tau q)_2 \gt 0$ and $\lambda _3^R \lt 0$ , viscous forcing contributes constructively to wave propagation, leading to amplification through mixing with the upper layer.

Figure 10.

Control mechanisms that dictate how information propagates along the characteristics in the three-layer model at $t = 110$ for all cases: (a) viscous control and (b) the sum of all terms on the right-hand side of (3.1). Solid lines ( $\zeta _1$ ) and dashed lines ( $\zeta _2$ ) represent the first and second momentum equations (first and second rows, respectively) in (3.1). The lines are shown for the SW, TW and I cases in light purple, dark purple and yellow, respectively.

Thus far our focus has been primarily on the physical significance of characteristics in hydraulic control and regime transition embedded in $\lambda ^R$ . Although we showed a strong correlation between turbulence generation and $\lambda ^R$ , to reveal their true connections, we need to consider instabilities through the imaginary part $\lambda ^I$ , as we do next.

Figure 11.

Space–time plots of the long-wave growth rate simultaneity, $\lambda _2^I \lambda _4^I$ , for (a) SW, (b) TW and (c) I cases.

4.3.1. Diagnosing long-wave resonance across interfaces from DNS

In matrix perturbation theory (Stewart & Sun Reference Stewart and Sun1990), the term spectral gap is often used to describe the difference between the spectra of two matrix pairs. Given a matrix pair $(\unicode{x1D63C}, \unicode{x1D63E}{\kern2pt})$ , where we consider the generalised eigenvalue problem $\unicode{x1D63C} \boldsymbol{v} = \lambda \unicode{x1D63E} \boldsymbol{v}$ , a perturbation of the matrix pair results in a perturbed generalised eigenvalue problem $(\unicode{x1D63C} + \Delta \unicode{x1D63C}, \unicode{x1D63E} + \Delta \unicode{x1D63E}{\kern2pt})$ . The spectral gap in this context quantifies the variation in eigenvalues caused by the perturbation of both matrices in the pair. This measure quantifies the effects of small changes to both matrices on the generalised eigenvalues and, consequently, the stability and dynamics of the system. In the context of perturbed upper and lower layers (represented by matrix pairs $\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}$ and $\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}$ , respectively), the spectral gap measures how the spectra of these perturbed layers deviate from each other and from the isolated middle layer represented by $\unicode{x1D655}$ . By examining the spectral gap, we gain insight into how resonant interactions between perturbed layers influence the characteristics of the full three-layer system. Further explanations and mathematical descriptions for the spectral gap in the context of perturbed layers are provided in Appendix B.

To quantify such resonant interactions in DNS data, we measure the relative spectral gap between (characteristics of) perturbed upper and lower layers with respect to the (characteristics of) isolated middle layer as

(4.8) \begin{eqnarray} \mathcal{R}(\unicode{x1D652},\unicode{x1D654}{\kern2pt}) = \max _{i}\, \max _{j} \dfrac {{\big\lvert {\, {\big\lvert {\Delta {\lambda ^{\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}}_i}}\big\rvert }\sqrt {1+{\big\lvert {\lambda _j({\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}}{\kern2pt})}\big\rvert }^2}-{\big\lvert {\Delta {\lambda ^{\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}}_j}}\big\rvert }\sqrt {1+{\big\lvert {\lambda _i({\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}}{\kern2pt})}\big\rvert }^2} \ } \big\rvert }}{{\big\lvert {\Delta {\lambda ^{\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}}_i}} \big\rvert }\sqrt {1+{\big\lvert {\lambda _j({\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}}{\kern2pt})} \big\rvert }^2}+{\big\lvert {\Delta {\lambda ^{\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}}_j}} \big\rvert }\sqrt {1+{\big\lvert {\lambda _i({\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}}{\kern2pt})} \big\rvert }^2}}, \end{eqnarray}

where $\Delta {\lambda ^{\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}}_i}= [\lambda ^R_i(\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}{\kern2pt}) -\lambda (\unicode{x1D655}{\kern2pt}) ]+ I \lambda ^I_i(\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}{\kern2pt})$ and $\Delta {\lambda ^{\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}}_i}= [\lambda ^R_i(\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}{\kern2pt}) -\lambda (\unicode{x1D655}{\kern2pt}) ]+ I \lambda ^I_i(\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}{\kern2pt})$ . The superscripts $R$ and $I$ refer to real and imaginary components of the corresponding characteristics and $i,j=1,\ldots, 4$ are indices counting individual characteristics. Note that $\mathcal{R}(\unicode{x1D652},\unicode{x1D654}{\kern2pt})=\mathcal{R}(\unicode{x1D654},\unicode{x1D652}{\kern2pt})$ . The reader is referred to Appendix B.2 to see a further description and mathematical derivation of (4.8).

Equation (4.8) yields that under the condition where $\lambda ^R(\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}{\kern2pt}) =\lambda ^R(\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}{\kern2pt}) =\lambda (\unicode{x1D655}{\kern2pt})$ and $\lambda ^I(\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}{\kern2pt})=\lambda ^I(\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}{\kern2pt}) \gt 0$ , then $\mathcal{R}(\unicode{x1D652},\unicode{x1D654}{\kern2pt})=0$ , which implies that long waves can resonantly interact since they are phase locked and have an equal growth rate. Therefore, $\mathcal{R} \rightarrow 0$ (note that $0 \leqslant \mathcal{R} \leqslant 1$ ) is a condition for near resonance, whereas when $\mathcal{R} \rightarrow 1$ , long waves are not phase locked and do not resonate. Note that for the special case of pure exchange flow discussed in § 3.2.2, the $F_1^2=F_{\varDelta +}^2$ condition that defines the neutral stability curve for long waves also leads to $\lambda (\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}{\kern2pt})=\lambda (\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}{\kern2pt})$ , which corresponds to the resonance condition $\mathcal{R} = 0$ (since   $\lambda_1=\lambda_3$ and $\lambda_2=\lambda_4$ ). Equivalently, $\mathcal{R} = 0$ marks the neutral stability of the long waves.

This equivalence highlights a connection between the resonance mechanism and the onset of long-wave instability in pure exchange flows. For more general three-layer exchange flows, where the neutral stability curve cannot be expressed analytically, we expect the resonance measure $\mathcal{R}$ to continue to serve as a useful diagnostic for identifying resonant interactions between interfacial waves. To confirm this, we apply the $\mathcal{R}$ measure to DNS data in the following section.

Figure 12.

Contours of the long-wave resonance indicator $\mathcal{R}$ plotted in space–time $(x, t)$ for the SW (a), TW (b) and I (c) cases. The dark regions indicate the minimum spectral gap, as defined in (4.8), between the perturbed upper and lower layers, derived from the building-block matrix decomposition of the three-layer system into perturbed upper and lower layer subsystems in (4.1).

4.3.2. Application to the data

The contours of the long-wave resonance indicator $\mathcal{R}$ are shown in figure 12 for the SW, TW and I cases. All cases qualitatively correspond to the $\widetilde {G}$ contours in figure 9(d–f), particularly in the sense that the regime transitions from a supercritical hydraulic regime ( $\widetilde {G} \gt 1$ ) to a subcritical regime ( $\widetilde {G} \lt 1$ ) when $\mathcal{R} \approx 0$ . This similarity is noteworthy and not self-evident as the characteristics of the perturbed upper layer and those of the upper interface in the original three-layer flow are different by definition, i.e. $\lambda (\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}{\kern2pt}) \neq \lambda _{1,2}(\,\widehat {\!{\kern-1pt}\unicode{x1D653}}{\kern2pt})$ , and, similarly, $\lambda (\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}{\kern2pt}) \neq \lambda _{3,4}(\,\widehat {\!{\kern-1pt}\unicode{x1D653}}{\kern2pt})$ in general. Therefore, one would expect that $\lambda (\,\widetilde {\!{\kern-1pt}\unicode{x1D654}}{\kern2pt}) + \lambda (\,\widetilde {\!{\kern-1pt}\unicode{x1D652}}{\kern2pt}) \neq \lambda (\,\widehat {\!{\kern-1pt}\unicode{x1D653}}{\kern2pt})$ . The striking resemblance between the spatiotemporal distributions of $\mathcal{R}$ and $\widetilde {G}$ suggests that the building-block decomposition introduced in § 4.2 provides valuable insight not only into the pathway for growth of interfacial waves but also their interactions and three-layer hydraulics. Moreover, it suggests that the hydraulic regime transition, driven by changes in the direction of information propagation along the interfaces, is also accompanied by long-wave resonance.

This long–long-wave resonance will modify the base state (mean flow) by the long waves themselves (figure 2). Such changes in the base state may lead to the onset of short-wave instabilities and result in the TKE generation at later times, as also found in Zhu et al. (Reference Zhu, Atoufi, Lefauve, Kerswell and Linden2024). For all cases, the contours of $\mathcal{R}$ are consistent with the TKE contours presented in figure 9 in the sense that strong TKE events are linked with interfacial wave resonances (i.e. the strong TKE occurs generally in places with the smallest $\mathcal{R}$ ). A more detailed comparison between instantaneous profiles of TKE and $\mathcal{R}$ is presented in Appendix C.

To properly establish the link between long-wave resonance and TKE generation, we also identify the resonance between long-wave and short-wave packets, eventually leading to the generation of small scales, which will be discussed next.

4.4.1. Linking long-wave and short-wave instability: hydraulically controlled instability

In the previous section we examined the link between internal hydraulics and the instability of long waves. We now show that the stability of the entire wave spectrum – including short and intermediate waves – is likewise governed by the internal hydraulic state, characterised by the layer Froude numbers. This mechanism, which we refer to as hydraulically controlled instability, highlights that although the system is a stratified shear flow both the shear and stratification are prescribed by hydraulic conditions, which in turn control stability across all wavenumbers.

To illustrate this, we again consider the pure exchange flow case discussed in § 3.2.2, which simplifies the analysis. The imaginary components of the phase speeds of all waves in the spectrum, derived from the linear theory developed in Appendix A, are shown in figure 13. The top row corresponds to the bifurcation Froude number $F_{\varDelta -}^2$ and the bottom row to the marginal-stability Froude number $F_{\varDelta +}^2$ . As the Froude number increases from $F_{\varDelta -}^2$ to $F_{\varDelta +}^2$ , the wave instability structure changes qualitatively across all wavenumbers – not just in the long-wave limit.

Figure 13(a,b) shows that when $F_1^2 = F_{\varDelta -}^2$ , unstable modes span a broad wavenumber range – from long to short waves – up to $k \lesssim 10$ , particularly in three-layer flows with a thick middle layer ( $r \geqslant 1$ ), where $r$ denotes the ratio of middle to upper layer thickness. Moreover, the phase speeds $c_1$ and $c_4$ exhibit similar imaginary components across wavenumbers, indicating that interfacial waves at both interfaces grow at nearly the same rate. In contrast, when $F_1^2 = F_{\varDelta +}^2$ , corresponding to marginal long-wave stability, instability is confined to short waves, which persist across a wide range of $r$ , as shown in figure 13(c,d). In this case, the $c_1^I$ and $c_4^I$ differ significantly, both in the values of $r$ that support instability and in the wavenumber ranges associated with peak $c_1^I$ and $c_4^I$ , reflecting a decoupling of instability mechanisms across the layers.

These results suggest that the layer Froude numbers dictate both the onset and extent of instability across the wave spectrum. The observed similarity in the $c_1^I$ and $c_4^I$ in the $F_1^2=F_{\varDelta -}^2$ case between long- and short-wave behaviour raises the question of whether these waves can interact resonantly – a possibility we now explore.

Figure 13.

Instability map from contours of the imaginary component of the phase speed, $c^I$ , for linear disturbance waves in a three-layer pure exchange flow with a stagnant middle layer introduced in § 3.2.2. Both the bifurcation Froude number case $F_1^2 = F_{\varDelta -}^2$ and the marginal-stability Froude number case $F_1^2 = F_{\varDelta +}^2$ defined in (3.16) are shown, isolating the dependence of wave amplification on the layer depth ratio $r$ and wavenumber $k$ . For $F_{\varDelta -}^2$ , the interfacial modes $c_1$ and $c_4$ exhibit similar values of $c^I$ , indicating comparable growth rates at both interfaces. In contrast, for $F_{\varDelta +}^2$ , the growth characteristics differ between interfaces, with distinct ranges of $r$ and $k$ contributing to instability.

Figure 14.

Imaginary component of the phase speed of the lower interfacial wave, $c_4^I$ , for pure exchange flow with a stagnant middle layer introduced in § 3.2.2 at selected wavenumbers $k = 0.001$ (a), $1$ (b) and $10$ (c), plotted as a function of layer thickness ratio $r$ and Froude number $F_1^2$ .

4.4.2. Resonance between long and short waves

To assess the potential for resonance between long and short waves, figure 14 shows $c_4^I$ for the pure exchange flow at selected wavenumbers ( $k = 0.001$ , $1$ and $10$ ), across a range of layer thickness ratios $r$ and Froude numbers. At $k = 0.001$ shown in figure 14(a), $c_4^I$ closely matches its value in the long-wave limit ( $k \to 0$ ) shown previously in figure 13(b), confirming that the solution represents a long-wave mode. As can be seen from figure 14(a–c), for $k \lt 1$ , the variation in phase speed as the wavenumber changes is relatively weak, while for $k \gt 1$ , the phase speed continues to change significantly. This suggests that dispersion is strongest around $k = 1$ , where small differences in wavenumber lead to substantial differences in propagation speed. This behaviour implies that wave packets centred near $k = 1$ are highly dispersive – that is, their constituent wavenumbers propagate at different phase speeds and, consequently, different group velocities. The implications of this for resonance and energy exchange will be explored in the following section. Additionally, increasing $k$ broadens the region of instability as shown in figure 14(a–c), extending to lower Froude numbers ( $F_1^2\lt F_{\varDelta +}^2$ ) and a wider range of $r$ , indicating that short-wave instability becomes more prominent where $F_1^2\lt F_{\varDelta +}^2$ .

For systematic identification of resonance between long and short waves, let us recall that a packet of linear waves with phase velocity $c$ propagates at the group velocity,

(4.9) \begin{equation} c_g= \frac {\partial (ck)}{\partial k}=c + k \frac {\partial c}{\partial k}=c+k d_g, \end{equation}

where $d_g = \partial c/\partial k$ measures the wave dispersivity. The kinetic energy is carried by the group velocity $c_g$ as the wave packet propagates, and $c$ is the phase velocity, whose mathematical expression is defined in Appendix A (from (A10)). Of particular interest is a packet of short waves with a group velocity close to the phase speed of the long waves (i.e. $c_g = \lambda$ ), since long and short waves can then resonate (Ma Reference Ma1981).

In general, $c \in \mathbb{C}$ , so $c_g$ also contains real and imaginary parts. The real part $c_g^R$ represents the propagation speed of the wave envelope, while the imaginary part $c_g^I$ corresponds to its temporal growth. Resonance between wave packets and long interfacial waves is most likely when both $c_g^R \approx \lambda ^R$ and $c_g^I \approx \lambda ^I$ .

We utilise the general form of the spectral gap, as defined in (B13), to assess the relative group velocity of short wave packets compared with the phase speed of long waves at both the lower and upper interfaces. For a given wavenumber $k_0$ , we interpret the group velocity $c_g|_{(k=k_0)} \in \mathbb{C}$ as the generalised eigenvalue of a perturbed matrix pair involving the three-layer matrix pair $\,\widehat {\!{\kern-1pt}\unicode{x1D653}} = (\unicode{x1D63C}, \unicode{x1D63E}{\kern2pt})$ , expressed as

(4.10) \begin{equation} {c_g|_{(k=k_0)} = \lambda (\,\widehat {\!{\kern-1pt}\unicode{x1D653}} + {\unicode{x1D648}{\kern2pt}''}),} \end{equation}

where $\unicode{x1D648}{\kern2pt}''$ is a perturbation matrix pair that accounts for the finite-wavenumber modulation of the wave packet. We recall that $\lambda (\,\widehat {\!{\kern-1pt}\unicode{x1D653}}{\kern2pt})$ yields the phase speed and growth rate of long interfacial waves in the three-layer flow, while $\lambda (\,\widehat {\!{\kern-1pt}\unicode{x1D653}} + {\unicode{x1D648}{\kern2pt}''})$ provides the corresponding quantities for the amplitude envelope of a wave packet centred around $k_0$ . Using this interpretation, we apply the spectral gap definition from (B13) to compute

(4.11) \begin{equation} \mathcal{S} \equiv S_{\,\widehat {\!{\kern-1pt}\unicode{x1D653}}}({\unicode{x1D648}{\kern2pt}''}), \end{equation}

i.e. the spectral gap between long waves and short wave packets, by substituting $\unicode{x1D648} = \,\widehat {\!{\kern-1pt}\unicode{x1D653}}$ and ${\unicode{x1D648}{\kern2pt}'} = {\unicode{x1D648}{\kern2pt}''}$ in (B13).

The spectral gap $\mathcal{S}$ thus provides a diagnostic criterion for identifying near-resonant interactions between short-wave packets and long interfacial waves. Unlike the dispersion relation, which describes the behaviour of individual modes, $\mathcal{S}$ quantifies the extent to which a short-wave packet can propagate at the same speed as a long wave. When this condition is met, efficient energy transfer from the long-wave to the short-wave packet becomes possible. In the next section we examine how this resonance mechanism contributes to turbulence generation, particularly through the amplification of short-wave disturbances in regions where $\mathcal{S}$ is small.

4.4.3. The role of resonance in turbulence generation

The product $c_2^I c_4^I$ , shown in figure 17(c) for the TW and I cases, highlights regions in $(x, k)$ space where both interfacial modes are simultaneously unstable. These regions of modal amplification are especially pronounced for $k \lesssim 1$ and spatially coincide with zones of elevated TKE in figure 9. Although the peak values of $c_2^I$ and $c_4^I$ do not occur at the same $x$ for all wavenumbers, their product effectively identifies regions of overlap where both modes contribute to instability. This spatial alignment suggests that energy amplification is not solely attributable to isolated modal growth, but may be enhanced by interactions between wave modes on both interfaces. In particular, simultaneous instability at the upper and lower interfaces permits resonant coupling between long interfacial waves and dispersive short-wave packets, enabling efficient energy transfer across scales.

To understand the conditions that allow this energy transfer, we consider the limitations of the dispersion relation alone. While the dispersion relation derived in Appendix A plays a central role in identifying both long- and short-wave instabilities, it does not capture the interactions between modes at different scales. For instance, figure 17 shows that unstable modes exist on both interfaces over a broad wavenumber range. However, the dispersion relation does not indicate whether these modes are phase locked, whether their group velocities coincide or whether energy can be transferred between them.

This motivates the introduction of the spectral gap $\mathcal{S}$ , which quantifies the proximity between the group velocity of a short-wave packet and the phase speed of a long wave. When $\mathcal{S} \rightarrow 0$ , a near-resonant condition is established: the short-wave packet is not only unstable but can also efficiently extract energy from long energetic waves. This leads to resonance-driven amplification of short-wave disturbances, even in parameter regimes where their modal instability alone may not explain the observed turbulence.

4.4.4. Application to the data

To apply the theoretical framework developed above, we must select a representative wavenumber $k_0$ around which to evaluate the spectral gap $\mathcal{S}$ . This choice is guided by both dispersivity analysis and prior evidence of potential resonance. In particular, from the discussion in § 4.4.2, we have already seen that short-wave packets centred at $k_0 = 1$ exhibit phase speeds comparable to those of long interfacial waves, suggesting conditions favourable for resonance.

This is further supported by the dispersivity analysis shown in figure 18, which reveals that the wave dispersivity $|d_g|$ peaks at $k_0 = 1$ , indicating the strongest dispersion and the wavenumber, therefore, is effective for coherent wave packet formation that carries energy through the flow. It is noteworthy that the dispersivity of small-amplitude long waves tends to zero as $k \rightarrow 0$ , since these waves are non-dispersive. In contrast, for short-wave packets, the dispersivity peaks at a wavenumber $k_0$ , where $|d_g|$ is maximal. Importantly, at $k_0 = 1$ , the phase speed $c(k_0)$ remains $\mathcal{O}(1)$ , which is comparable to the long-wave phase speed $\lambda$ . Furthermore, since $k_0 |d_g| \ll |c|$ , the group velocity $c_g$ is close to $c$ , and hence, also comparable to $\lambda$ . Therefore, wave packets centred at $k_0=1$ travel at approximately the same speed as long waves, facilitating resonant coupling.

Figure 15 shows space–time contours of $\mathcal{S}$ and demonstrates it qualitatively corresponds to $\mathcal{R}$ and TKE for all cases despite their fundamentally different mathematical definitions. For example, figure 15(a,b) suggests that in the SW and TW cases, $\mathcal{S}$ is smallest in the central regions of the duct and, therefore, the relative velocity between packets of short waves and long waves is small in this region. Consequently, in SID, in the central region where long interfacial waves are in near resonance, they also resonate with energy-carrying short-wave packets, amplifying the growth of the unstable short waves, resulting in breaking and generating small-scale turbulence and TKE. A more detailed comparison of the spatial distributions of TKE, $\mathcal{R}$ and $\mathcal{S}$ for the SW, TW and I cases at various times is presented in Appendix C.

Therefore, $\mathcal{S}$ and $\mathcal{R}$ serve as complementary diagnostics to unravel the route to turbulence and small-scale generation in a stratified exchange flow such as SID. The next section summarises these mechanisms.

4.6.1. Implementations for adaptive resolution in ocean models

Recent advances in high-resolution numerical ocean modelling, particularly the lat-lon-cap (LLC) simulations such as LLC4320 (Gallmeier et al. Reference Gallmeier, Su, Rocha, Menemenlis and Klein2023), have enabled global ocean simulations at approximately $1/48^{\circ }$ horizontal resolution ( $\sim\! 2$ km at the equator) using the Massachusetts Institute of Technology general circulation model (MITgcm). Despite this relatively high spatial resolution, turbulence is far from being resolved, and mixing remains parameterised through local gradient-based closures.

The diagnostic measures $\mathcal{R}$ and $\mathcal{S}$ developed in this study quantify the conditions under which long waves resonate with dispersive short-wave packets. Since these diagnostics rely only on layer-averaged quantities – velocities, densities and interface displacements – they are directly applicable to model outputs.

In coastal and marginal seas, particularly in estuarine and strait-like geometries where the Rossby number is large and rotational effects are weak, the three-layer approximation proposed in this study is especially appropriate. Our results imply that even without directly resolving turbulence, ocean models can leverage these diagnostics to identify zones of intensified TKE generation. This insight may be valuable for implementing adaptive mesh refinement in regional models, where local grid refinement could be concentrated in dynamically active zones suggested by long–short wave resonance conditions. Because all diagnostic expressions derived here are analytical, this integration poses no additional computational burden and presents a promising route toward physics-informed resolution enhancement in estuarine and strait flows.

4.6.2. Implications for mixing parametrisation

While in this study we applied the theoretical framework to DNS datasets of SID flows, the methodology can be extended to ocean circulation models, which typically solve shallow water or layer-averaged primitive equations. These models do not resolve turbulence explicitly and usually assume hydrostatic balance, thus neglecting non-hydrostatic pressure gradients that play a key role in wave-driven turbulence and mixing.

The model output can be decomposed into three effective layers by discretising the vertical structure into bulk velocity, density and thickness in each layer via appropriate integration. These layer-averaged quantities serve as input to the theory developed here, enabling the computation of resonance diagnostics $\mathcal{R}$ and $\mathcal{S}$ , the linear wave speed $c$ , the imaginary part $c_i$ as a proxy for instability and the dispersivity $|d_g|$ , defined as the magnitude of the group velocity derivative with respect to wavenumber. In regions where $\mathcal{S}$ is locally minimum or $|d_g|$ is large, we suggest that turbulent viscosity in mixing parametrisations could adopt a velocity scale based on the group velocity of the most dispersive wave packet, $c_g(k_0)$ , where $k_0$ corresponds to the wavenumber at which $|d_g|$ peaks.

This provides a leading-order correction to existing mixing schemes that rely solely on local shear or stratification. It leverages the resolved mean state of an ocean model to estimate non-hydrostatic, wave-induced pathways to turbulence, which are otherwise absent in hydrostatic model formulations. Consequently, it offers a route to enhance parametrisation schemes with minimal computational overhead while capturing key non-local mixing mechanisms.