Long-wave instabilities of sloping stratified exchange flows

Abstract We investigate the linear instability of two-layer stratified shear flows in a sloping two-dimensional channel, subject to non-zero longitudinal gravitational forces. We reveal three previously unknown instabilities, distinct from the well-known Kelvin–Helmholtz instability and Holmboe wave instability, in that they have longer wavelengths (of the order of 10 to $10^3$ shear-layer depths) and often slower growth rates. Importantly, they can grow in background flows with gradient Richardson number $\gg 1$, which offers a new mechanism to sustain turbulence and mixing in strongly stratified flows. These instabilities are shown to be generic and relatively insensitive to Reynolds number, Prandtl number, base flow profile and boundary conditions. The nonlinear evolution of these instabilities is investigated through a forced direct numerical simulation, in which the background momentum and density are sustained. The growth of long unstable waves in background flows initially stable to short wave causes a decrease in the local gradient Richardson number. This leads to local nonlinear processes that result in small-scale overturns resembling Kelvin–Helmholtz billows. Our results establish a new energy exchange pathway, where the mean kinetic energy of a strongly stratified flow is extracted by primary unstable long waves and secondary short waves, and subsequently dissipated into internal energy.


Introduction
The study of stratified flows has attracted considerable attention over the past few decades due to their importance in many environmental and industrial processes.In the oceans, stratification occurs due to differences in salinity and/or temperature, leading to mostly stably stratified flows.Turbulence in these flows plays a significant role in the transport of momentum and mass and is crucial in shaping the global climate (Linden 1979;Riley & Lelong 2000;Gregg et al. 2018;Caulfield 2020).An interesting open question concerns the maintenance of turbulence and its associated irreversible turbulent mixing under strong stable stratification, which tends to suppress turbulence.
When stratification is relatively weak, stably-stratified flows can be linearly unstable.It is well known that linear shear instabilities, such as Kelvin-Helmholtz instability (KHI) (Hazel 1972;Smyth et al. 1988) and Holmboe wave instability (HWI) (Holmboe 1962), can cause transition of a laminar stratified flow to turbulence, inducing strong mass and momentum transport (Caulfield 2021).Over the past 50 years, numerous studies have been carried out to understand these instabilities and their relation to mixing (Thorpe 1968;Smyth et al. 1988;Carpenter et al. 2010;Salehipour et al. 2015;Zhou et al. 2017).In most of these studies, the density isopycnals are perpendicular to the direction of gravity, which does not explicitly drive the flows.
However, in many natural systems, density isopycnals are not exactly perpendicular to gravity, in which case nonzero streamwise gravity forces come into play and may partially drive the flow.One notable example is the internal tide interacting with the sloping bottom topography of the oceanic continental shelf (Garrett & Kunze 2007).At a critical slope, the internal tide provides an additional energy production pathway that leads to turbulent mixing of temperature, salinity, and other tracers (Gayen & Sarkar 2010).Similarly, many engineering flows occur along an inclined boundary.Examples can be found in building ventilation systems (Linden 1999), where indoor/outdoor air is often exchanged through inclined ventilation ducts, producing mixing and dispersion of heat and indoor pollutants.In gas-cooled nuclear reactors, carbon dioxide and air are exchanged through inclined coolant ducts, which can result in the depressurization and damage of the reactors in case of failure (Leach & Thompson 1975;Mercer & Thompson 1975).
Studies on the influence of longitudinal gravitational forcing on the onset of turbulence in stratified exchange flows remain limited.One notable recent body of work is the Stratified Inclined Duct (SID) experiment (Meyer & Linden 2014;Lefauve et al. 2019;Lefauve & Linden 2020).These studies investigated the transition and turbulent mixing of the exchange flow in an inclined duct that connected two reservoirs with fluids at different densities or temperatures.To understand the mechanism of transition in SID, Lefauve et al. (2018) conducted a linear stability analysis using a base state extracted from the SID experiment.Subsequently, Ducimetière et al. (2021) systematically investigated the three-dimensional unstable modes in inclined ducts, focusing on the effects of side wall confinement.These studies focused primarily on HWI (and secondarily on KHI), which have wavelengths comparable to the thickness of the shear layers.Interestingly, Ducimetière et al. (2021) observed a secondary instability at significantly longer wavelengths than KHI and HWI and attributed it to the effect of the inclination angle.Recently, Atoufi et al. (2023) studied the mechanism of transition by applying shallow water equations as a diagnostic tool to analyse a new numerical database of SID (Zhu et al. 2023).They suggested that the instability of long shallow water waves (long-wave KHI in the presence of top and bottom solid boundaries) may cause turbulence in the SID.Although the longitudinal gravitational forcing was included in the numerical simulation data, it was not included explicitly in the shallow water model.
In this paper, we explore explicitly the impact of longitudinal gravitational forces on the instability of long waves and on potential new pathways toward turbulence, restricting ourselves to a two-dimensional geometry.In § 3, we examine the linear instabilities in inclined channels and conduct a thorough exploration of the parameter space.We identify three new distinct families of long-wave instabilities distinct from the well-known HWI and KHI, and map in parameter space these long-wave instabilities that dominate the flow.In § 4, we then investigate the evolution of these new instabilities by conducting two-dimensional forced direct numerical simulations (DNS), and discuss their impact on turbulence and energy transfers.Finally, we conclude in § 5.

Problem formulation and governing equations
In this section, we present the equations required for linear stability analysis (LSA) of a stratified exchange flow between two fluid layers having density  0 ± Δ/2 (where  0 is the reference density and 0 < Δ ≪  0 is the density difference) in a two-dimensional stratified inclined channel (SIC, see fig.1(a)).Following the SID experimental literature, lengths are nondimensionalized by the half-channel height  * , velocity by the buoyancy-velocity scale  * ≡ √  ′  * (where  ′ = Δ/ 0 is the reduced gravity), time by the advective time unit  * / * , pressure by  0  * 2 , and density variations around  0 by Δ/2, respectively.The nondimensional continuity, Navier-Stokes and scalar governing equations under the Boussinesq approximation are where u = (, , ) is the non-dimensional velocity in the three-dimensional coordinate system x = (, , ), where -,-,-axis are the longitude, spanwise and wall-normal direction of the channel respectively.In this coordinate system gravity g is pointing downward at a angle  to the − axis, i.e., g =  ĝ =  [sin , 0, − cos ], and  and  are the non-dimensional pressure and density, respectively.The dimensionless parameters are the Reynolds number Re ≡  *  * / ( is the kinematic viscosity), the Prandtl number Pr ≡ / ( is the scalar diffusivity), and Richardson number Ri ≡  ′  * /(2 * ) 2 = 1/4 (fixed here because of the buoyancy velocity scale).

Formulation of linear stability analysis
We now apply a linear stability analysis (LSA) (Drazin & Reid 2004;Smyth & Carpenter 2019) to the SIC, noting that in agreement with Squire's theorem, Lefauve et al. (2018);Ducimetière et al. (2021) have shown that the fastest-growing mode is two dimensional (2D).We impose infinitesimal 2D perturbations to a 1D base state.The velocity, density, and pressure fields are thus decomposed as (2.4) where capital letters and superscript prime represent the mean and perturbation components of quantities, respectively.A normal mode perturbation of the form (, , ) = φ() exp ( + ), (2.7) is adopted.The base flows are obtained by solving for the numerical solution of the laminar exchange flow following Thorpe (1968), which will be introduced in §2.3.Substituting (2.4)-(2.6)into (2.1)-(2.3) and linearising yields the same system as Lefauve et al. (2018), i.e.
where 0 and I are the zero and identity matrices, respectively and where Δ =  2 −  2 (the operator  = / and  2 =  2 / 2 ).At the top and bottom boundaries ( = ±1), no-slip and no-flux boundary conditions are applied for velocity and density, respectively.We also demonstrate the negligible effect of choosing a free-slip boundary condition for velocity in appendix A. To obtain the unstable modes, we solve the linear system (2.8) numerically using a second-order finite-difference discretization method described in Smyth & Carpenter (2019).The spatial resolution is chosen based on the sharpness of the interface and is (150,150,250,400) grid points for Pr = (1, 7, 28, 70), respectively.A sensitivity analysis for resolution ensured convergence of the results.

Base flows
The base state for density in our exchange flow is taken as a hyperbolic tangent (figure 1 (2.10) The interfacial thickness is  = 1/(2 √ Pr) to approximate the effect of diffusion (Smyth & Peltier 1991).The typical model (Smyth et al. 1988, e.g.) considers a shear layer driven by an arbitrary, controllable background shear.A similar procedure is applied to our SIC by modifying the laminar solution developed by Thorpe (1968) and imposing a background body force F = − (where  is a variable to control the magnitude of the force).This decouples the base velocity from the inclination angle in SIC, allowing for the exploration of the  −  space, as if being influenced by arbitrary external tidal forces or pressure gradients.The mean velocity profile  () of the steady laminar exchange flow is obtained by integrating the 2D momentum equation (2.11) where −/ = 0 to satisfy the zero-flux condition of SIC.This yields the following laminar base state for the forced SIC where  2 is the polylogarithm function of order 2. The constants  1 and  2 are computed given the no-slip boundary condition at the walls  ( = ±1) = 0 and are (2.14) This solution  () is sinusoidal-like (figure 1(b)), much like those observed in experiments and simulations (Lefauve et al. 2018;Zhu et al. 2023).The magnitude of the base velocity depends on Re, , and , while the shape depends more on .In addition to the base state described by (2.12), we also conducted a LSA with a tanh-shape velocity profile in appendix A, to compare with the standard stratified free-shear layer model (Smyth et al. 1988).These results were qualitatively consistent with those in the remainder of the paper, in terms of the existence of the same long-and short-wave families in SIC.

Results: new families of linear instabilities in SIC
Here we present the results from the LSA of SIC.We explore the parameter space of  −  and map out three new families of long-wave instabilities in addition to the well-known short-wavel HWI and KHI.We also investigate the impacts of Re and Pr in order to further understand the importance of these newly discovered long waves in the laminar-turbulence transition.

Five families of instabilities
We first fix (Re, Pr, Ri) = (1000, 7, 0.25) and vary the inclination angle  from −10 • to 10 • .When  > 0, the SIC slopes downward, the streamwise gravity energises the mean flow and vice versa.We vary the forcing factor , on which two important physical quantities depend: the interfacial background Richardson number   , defined as the gradient Richardson number of the background flow at the density interface  = 0, i.e., and the mass flux (or flow rate of buoyancy), which is given by The Richardson number   is an important measure of the relative importance of stratification compared with shear, which is critical for stratified shear flow stability (Caulfield 2020).The mass flux   is closely associated with the hydraulic control of exchange flows; a threshold value of   ≈ 0.5 indicates the emergence of an internal hydraulic jump (Meyer & Linden 2014;Lefauve et al. 2019) which Atoufi et al. (2023) demonstrated to be equivalent to a relatively long KHI (requiring the existence of a top and bottom boundaries).
Figure 2(a,c) shows the distribution of the growth rate and wave frequency of the fastestgrowing modes in the parameter spaces (,   ) and (,   ), respectively.Examining the contour lines reveals five distinct families of unstable modes, shown schematically in figure 2(b,d).To better understand these modes, we show the dispersion relation of five based on their longer wavelengths ( (10 ∼ 10 4 )) compared to the 'short' HWI and KHI ( (10 −1 ∼ 10)).To the best of our knowledge, these unstable modes have not previously been investigated in the literature.We find that the features of these instabilities are generally insensitive to the shapes of base profile and boundary conditions, despite adopting a base profile (2.12) and no-slip boundary in this section.To support this, we show in appendix A that these instabilities are found using a tanh-shape base state and free-slip boundary condition, as used by Smyth & Winters (2003).This suggests that these instabilities can exist in a wide range of stratified exchange flows along a slope.In the following sections, we characterise the five families of unstable modes in more detail.

Holmboe wave instability (HWI)
The HWI (Holmboe 1962) occurs when the density interface is thinner than the shear layer and results from the resonance between vorticity waves at the edges of the shear layer and internal gravity waves at the density interface (Caulfield 1994;Carpenter et al. 2010).It gives rise to a pair of counter-propagating growing modes on either side of the density interface.
In SIC, the regime dominated by HWI exists from  = −10 • to 2 • and   = 0.3 to 4 (  = 0.1 to 0.3) in figure 2. The dispersion relation of HWI is shown in figure 3, where HWI has a pair of complex conjugate eigenvalues with non-zero phase speed  = −  /.Despite the well-known feature that HWI can exist in horizontal flows at   values significantly higher than 0.25 (Miles 1961;Howard 1961), we notice that HWI can also be induced over a wide range of .More interestingly, the HWI-dominated regime gradually shrinks from  < 0 to  ≈ 2, beyond which HWI ceases to exist.This indicates that increasing downward slopes have a negative effect on HWI, a phenomenon that has not been previously discussed in the literature and constitutes a new result.

Kelvin-Helmholtz instability (KHI)
The KHI arises due to the interaction of vorticity waves at two edges of finite shear layers, leading to a sequence of stationary vortex billows that roll up the denser fluids and cause significant mixing (Hazel 1972;Smyth et al. 1988).However, unlike these previous studies (with the exception of the recent Atoufi et al. (2023)) the KHI observed here in the SIC geometry is bounded by no-slip solid boundaries at  = ±1.
In SIC, KHI has a zero phase speed and a characteristic wavelength of , consistent with previous studies by Smyth & Carpenter (2019); Caulfield (2021); Smyth & Peltier (1991).KHI dominates the flow at small   ≲ 0.25, in agreement with the Miles-Howard criterion.Interestingly, like HWI, the longitudinal gravity force can affect the regimes of KHI.The upper bound of the KHI-dominant regime in figure 2(a) increases linearly from   = 0.15 to 0.25 as  increases from −10 to 10.This suggests an enhancement of KHI by a downward slope, which we believe to be an additional new result.

Long-wave instability (LWI)
Of the three new instabilities that arise with slopes, the novel long-wave instability (LWI) dominates the flow at large downward slopes ( > 4 • ) and a weak shear (strong stratification).In contrast to KHI and HWI, the LWI has a longer wavelength ( (10 − 10 2 )).Note that the LWI discussed in this paper is distinct from the long waves supported by shallow-water (hydraulic) theory (Lawrence 1990;Atoufi et al. 2023) which are essentially KH waves with a large  (satisfying the hydrostatic approximation) and which can exist at  = 0. LWI, on the other hand, specifically requires  ≠ 0. As depicted in figure 3, its phase speed is near-zero.This instability can be triggered at   ≫ 1, at which the shear-induced HWI and KHI vanish.Note that the presence of a mean shear can affect LWI by modifying its growth rate and phase speed.In terms of wave interaction, since vorticity waves vanish as   → 0, we hypothesise that LWI is a result of the interaction between two gravity waves at the density interface whose symmetry is broken by the non-zero slope.However, the   = 0 condition may be arbitrary when subjected to a non-zero slope, as it requires the gravity and pressure forces to be precisely cancelled by external body forces F in (2.11).In practice, such a precisely balanced condition is expected to be rarely observed.

Downslope very-long-wave instability (VLWI-DS)
The new VLWI-DS shares similarities with the LWI, in that it can exist at weak shear (strong stratification) and has a long wavelength.However, VLWI-DS dominates the flow under different conditions, namely when 2 • <  < 5 • and   > 0.25 (  < 0.5).It is also characterized by very long wavelengths of  (10 2 − 10 3 ) (wave numbers  =  (10 −3 ) ∼  (10 −2 )) and, interestingly, a pair of eigenmodes with complex conjugate phase speeds (figure 3).As with the HWI, we thus expect a pair of unstable VLWI-DS modes propagating with opposite phase speeds.The evolution of these unstable long waves and their connections to the onset of turbulence will be further discussed in §4.

Upslope very-long-wave instability (VLWI-US)
At a negative inclination angle ( < 0), i.e. for upward slopes, another type of very-long-wave instability (VLWI-US) appears with wavelengths ⩾ 10 2 (wave numbers  <  (10 −2 )) and a zero phase speed (figure 3).This instability is similar to LWI and VLWI-DS in that it requires a slope ( ≠ 0) and can exist in a strongly stratified environment.Contrary to the usually significantly smaller growth rate of the long waves compared with the corresponding short waves, the VLWI-US has in fact a comparable growth rate as HWI; this will be further discussed in §3.3.
Importantly, these long-wave instabilities have the potential to trigger and sustain turbulence in strongly-stable stratified flows, which are a priori regarded as stable.In §4 we will show that these new instabilities can indeed destabilise the flow at   ≫ 1, eventually resulting in nonlinear bursting and a transition to turbulence and mixing.It is also important to note that figure 2 only shows the fastest growing modes, whereas multiple families of instabilities can coexist in certain regions, as shown in figure 3.As a result, the regions of instability overlap, and the neutral boundary of each instability cannot be identified from figure 2. In §3.3, we will address this challenge by introducing an unsupervised clustering technique to isolate the neutral boundary of each family.Furthermore, in figure 2, we include a black line computed from  = 0, i.e. the natural convective 'Thorpe' base state with forcing F = 0.Under the parameters discussed so far (Re = 1000, Pr = 7), this line does not overlap with the regimes of long-wave instabilities in parameter space.Nonetheless, it is important to note that different Re and Pr or boundary conditions can modify the regimes of the long wavelength instability and interact with the base flow.An example is demonstrated in §3.4 for Pr = 28.
3.2.Eigenfunctions Further insights into these SIC instabilities can be gained by examining their eigenfunctions expressed in (2.7) for representative cases (see figure 2 and table 1).In figure 4, we present the vorticity (first row) and density (second row) eigenfunctions of the fastest growing modes for cases I, ..., V, marked in figure 2, each of which represents one of the five branches of instabilities: HWI, KHI, LWI, VLWI-DS, and VLWI-US, respectively.Note that the -axis in these cases has been re-scaled to compare modes having very different wavelengths.In figure 4, the wavelengths of HWI and KHI are ≈ 4, LWI is ≈ 70, VLWI-DS is ≈ 300, and VLWI-US is ≈ 420.
The density eigenfunctions of all modes are concentrated near the interface, indicating the critical role of stratification.Near the walls, the intensity of vorticity eigenfunctions is large due to the no-slip effects of the walls.(Note that, with a free-slip velocity boundary condition, the corresponding modes do not exhibit this intense vorticity at the wall, see appendix A).In the shear layer, one of the HWI modes plotted here (left-propagating) exhibits two pairs of counter-rotating roll cells centred at  ≈ 0.5.For KHI, the vorticity and density eigenfunctions are highly concentrated at the interface, leaving a weaker bulk region in the rest of channel.By contrast, the vorticity eigenfunctions of LWI, VLWI-DS, and VLWI-US fill the channel and are asymmetric with respect to  = 0.

Neutral boundaries of instabilities
As mentioned in §3.1, different families of instabilities can coexist at the same parameters, making it difficult to determine the neutral boundary of each family from the distribution of fastest growing modes in figure 2. To identify the different neutral boundaries we employ an unsupervised machine learning algorithm called DBSCAN (density-based spatial clustering of applications with noise) (Ester et al. 1996).The DBSCAN algorithm clusters the local maxima of the dispersion relation (figure 3(a)) of all the cases in figure 2 using ,   , and   as input variables.These variables are first logarithmically transformed and normalised before being fed into DBSCAN for clustering.Note that the DBSCAN groups the local optimal modes of LWI and VLWI-US together in a single cluster due to their similarity in ,   , and   .An additional step is taken to distinguish between the two branches by using the fact that LWI occurs when  > 0, while VLWI-US occurs when  < 0. The clustering analysis in figure 5(b-f) reveals the regimes of different families of instabilities, which could not have been identified by simply looking at the distribution of fastest amplifying modes (figure 5(a)).The KHI regime (panel (c)) exactly matches the distribution of the fastest amplifying modes (panel (a)), while other modes (LWI, VLWI-DS, VLWI-US) that overlap with KHI are omitted.This suggests that KHI always has the fastest growth rate.For HWI (panel (b)), increasing  clearly decreases the growth rate while shrinking its 'territory', causing it to disappear when  > 2 • .When  is fixed, the fastest growing HWI appears at   ≈ 1, while the growth rate decreases as   departs from 1.The territory of HWI overlaps with VLWI-US (panel (f)) which can exist when  < −0.5 • .The growth rate of these two modes is comparable so that figure 5(a) cannot display the neutral boundaries of these two modes properly.As for LWI (panel (d)), it generally persists at large positive  except for   ≲ 0.2.The critical  for the appearance of VLWI-DS (panel (e)) is ≈ 2.5 • .It overlaps with KHI and LWI at large  and small   , respectively, but is mostly omitted in the plot of the fastest growing mode due to its relatively small growth rate.
In general, these long-wave families of instabilities can persist across a wide range of   , ranging from   ≪ 0.25 (especially for VLWI-DS and VLWI-US) to   ≫ 1.Consequently, we anticipate their widespread presence in sloping stratified exchange flows.

Effect of Reynolds and Prandtl numbers
In this section, we study the impacts of Re and Pr on these different families of instabilities.

Reynolds number effects
Figure 6 shows the   −  parameter space of the fastest growing modes at a lower Re = 650 (panel (a)) and higher Re = 5000 (panel (b)) than the standard case discussed in §3.1.It is anticipated that in the inviscid limit Re → ∞ the critical  will approach 0 • .Therefore, it is a reasonable speculation that these gravity-induced long waves may be generic in high-Re natural water bodies subjected to shear, stratification and even the most shallow slope.

Prandtl number effects
Figure 7 displays the   −  parameter space of the fastest growing modes at Pr = 1, 28, and 70, respectively, corresponding to the increasingly sharper interface of the density base state, following (2.10).To ensure convergence, the grid resolution for the LSA was set to 150, 250, and 400, respectively.As Pr increases, the influence of the slope  on KHI becomes more significant, resulting in a wider upper boundary of KHI, which can be triggered at   > 0.25 for large downward slopes  ≳ 5 • .Meanwhile, HWI is also significantly affected by Pr.At Pr = 1, HWI does not appear due to the thick density interface determined by (2.10).However, as Pr increases, the region of HWI expands significantly towards larger .The long-wave families exist at all Pr.As Pr increases from Pr = 1 to 28, the territory of the long waves converges towards  = 0.However, the changes in the territory become less significant from Pr = 28 to 70, indicating a potential convergence of the wave regime at moderate Pr.However, due to the dominance of HWI at high Pr, the long-wave families are largely omitted by the fastest growing HWI at   ≲ 10 in figure 7(c).Interestingly, at Pr = 28, the profile of the Thorpe exchange flow (F = 0 in (2.12)) passes sequentially through the HWI, VLWI-DS, LWI, and KHI dominated regimes.This provides an example where VLWI-DS and LWI can dominate Thorpe's SIC flow.

Nonlinear evolution of unstable modes
To gain insight into the subsequent nonlinear evolution of these unstable modes we conduct forced two-dimensional direct numerical simulations (DNS).We describe our DNS in §4.1 and discuss the evolution and breakdown of the unstable flows in §4.2.The instantaneous flow kinetics of these unstable waves and the mechanisms leading to their breakdown are discussed in §4.3 and §4.4,respectively.

Forced DNS formulation
To simulate the growth of linear unstable perturbations on the desired base state, we add to the right-hand sides of (2.2) and (2.3) the two forcing terms respectively.In this way, the mean velocity and density of the DNS are forced towards the targeted base profile of  () and ().These terms can be regarded as enforcing a pressuredriven exchange flow under a sustained stratification.Similar approaches that apply body forces to the stratified flows were introduced in Taylor et al. (2016) and Smith et al. (2021).
We perform the simulations using the open-source solver Dedalus (Burns et al. 2020) employing a Fourier-Chebyshev pseudo-spectral scheme for spatial discretisation and a 3rdorder, 4-stage diagonally-implicit+explicit Runge-Kutta scheme (Ascher et al. 1997) for time stepping.We imposed periodic boundary conditions in the streamwise  direction, while we applied no-slip and no-flux boundary conditions for velocity and density, respectively, to the solid walls at  = ±1, as in the LSA.The streamwise length   of the channel was set equal to the wavelength of the fastest growing mode, while the channel height   was fixed at 2. We employed a uniform grid for the  direction and a Chebyshev grid for the  direction.The simulation resolution was determined by the geometrical and physical parameters of the problem.We initialised the simulations by superimposing on the base state the eigenfunctions of the LSA unstable modes with a perturbation magnitude .The parameters of the production runs are listed in table 1.

Temporal evolution
In this section, we focus on the temporal evolution of the fastest growing modes of each instability family, I, II, III, IV, and V, as marked in figure 2. Figure 8 shows the temporal behaviour of the unstable modes through the time series of the mass flux   () (3.2) and the spatially-averaged vertical velocity of perturbations ⟨ 2 ⟩(), where ⟨•⟩ denotes averaging over , .The magnitude  of the perturbation was chosen differently for each mode in order  to obtain a reasonably long linear growth period.The forcing magnitude  is determined so that the base velocity matches the selected cases in figure 2. The exchange flow is simulated by forcing the background flow in time using (4.1) and allowing the perturbations to grow.
In figure 8(a), the background state is controlled by the body forces (4.1) so that   initially remains constant and consistent with the targeted base state (as marked in figure 2) until the perturbations are significantly amplified and the flow enters the nonlinear stage.This initially constant   value indicates the effectiveness of the forcing method to maintain a sustained background state before the intense nonlinear dynamics set in.The evolution of the disturbance amplitudes is shown in figure 8(b), with all cases exhibiting a clear exponential growth period for  2 , with growth rates matching the corresponding linear unstable modes.For KHI, LWI, and VLWI-US, following the exponential growth period, an intense nonlinear bursting process is caused by the breakdown of the primary waves, leading to intense mixing and changes in   and  2 .In contrast, the sudden changes in   do not appear for HWI and VLWI-DS since their primary waves do not break down.HWI and VLWI-DS have a pair of conjugate modes, represented by oscillating  2 profiles, due to the synchronization of complex-conjugate modes, as discussed in Yang et al. (2022).Interestingly, after the nonlinear bursting at  = 1250, the nonlinear HWI still maintains the oscillating pattern (Lefauve et al. 2018).The time series of   and ln⟨ 2 ⟩ pinpoint the critical time when the nonlinear effects become prominent.Specifically, this occurs when   deviates from its constant level or when ln⟨ 2 ⟩ no longer shows exponential growth after reaching a certain amplitude.Note that the critical amplitude for nonlinear bursting remains independent of the initial amplitude of perturbations.However, it varies for each individual unstable mode, as illustrated in figure 8.The critical time may vary depending on the particular unstable mode, the growth rate, and the magnitude of the initial perturbation.
Figure 9 shows the  −  diagrams of ln⟨ 2 ⟩  , where ⟨•⟩  indicates −averages.In case HWI (panel (a)), we observe left-going waves from the spatial-temporal diagram, while its conjugate pair is omitted as only one mode of the pair is imposed as initial perturbation in the DNS.Nonlinear effects become significant at  ≈ 1250 (ln ⟨ 2 ⟩ ≈ −9) owning to a relatively small growth rate (0.0038), as indicated by the saturation of the exponential growth of  2 in figure 8(b).Interestingly, the spatial-temporal pattern of HWI does not change significantly after  ≈ 1250 in figure 9(a).This implies that nonlinear effects only halt the linear growth of HW structures, which maintain their forms as nonlinear HW, as observed in experiments and nature (Meyer & Linden 2014;Cudby & Lefauve 2021).By contrast, KHI generates strong secondary instabilities (Mashayek & Peltier 2012) after the onset of nonlinear effects at  ≈ 200 (ln ⟨ 2 ⟩ ≈ −5).Consequently, the KH billows break up, leading to highly chaotic flow stages with small-scale structures.
In figure 9(c-e), the  −  diagrams of ln⟨ 2 ⟩  reveal that for the three new families of long waves, small-scale structures emerge in the latter stages of the transitions, characterised by highly fluctuating contour lines.In the case of LWI (figure 9(c)), the onset of an intense chaotic flow period occurs at  = 710, during which small-scale structures are initially generated at  = 25 and  = 55 where ln⟨ 2 ⟩  of the linear wave peaks.These structures then propagate towards the quiet regions and ultimately trigger a disorganised flow field across the channel.Interestingly, we have observed from figure 8(a) that the nonlinear effects set in at  = 630 (ln ⟨ 2 ⟩ ≈ −14) as   significantly deviates from the original constant level.At this stage, the two peaks of the linear disturbance approach each other while contour curves twist.The ln⟨ 2 ⟩  distribution no longer maintains its shape as is in the linear growing period.As we will show later, this nonlinear dynamics is the breakdown of the long waves.
In case VLWI-DS, as shown in figure 9(d), the unstable wave moves leftward and grows exponentially until nonlinear dynamics set in at  = 2860 (ln ⟨ 2 ⟩ ≈ −19), which is identified by a jump in ⟨ 2 ⟩() in figure 8(b).Small-scale waves/structures are formed at a local peak of the long wave  2 , which propagate both leftward and rightward, creating strong mixing.Despite the intense bursting of the flows, the long wave does not break down like LWI, presumably due to its low growth rate.It continues to propagate at the same phase speed as the linear wave energy that was previously used to amplify the long wave and is then fed to the small-scale waves, which eventually break down and dissipate, allowing the long wave to persist for a long period of time and propagate over a long distance.
For VLWI-US (shown in figure 9(e)), local nonlinear bursting and small-scale structures are directly created at  = 500 (ln ⟨ 2 ⟩ ≈ −16) and  ≈ 150 and 350 on top of the long wave.Similar to VLWI-DS, a preliminary breakdown of the long-wave is not observed.Soon after, intense secondary instabilities fill the entire channel and the long waves are no longer distinguishable.
In summary, the evolution of these linear long waves eventually leads to the appearance of nonlinear dynamics and intense secondary short-wavelength structures.The formation of these small-scale structures is a result of perturbation amplification, which alters the base state and allows the growth of short-wave instabilities.We will explore this mechanism in more detail in §4.3-4.4.

Features of flow kinematics
In this section we present instantaneous flow fields corresponding to the key stages of evolution for the long wave instabilities.The kinematics of the short waves, i.e., HWI and KHI, have been well-documented in the literature (Smyth & Winters 2003;Salehipour et al. 2015;Mashayek & Peltier 2012, 2013;Lefauve et al. 2018), and will not be repeated here.
Figure 10 shows snapshots of the total density  =  +  ′ (colours) and vertical velocity  ′ (lines) at four stages of LWI.We find two stages of nonlinear breakdown corresponding to the breakdown of the long wave and the generation of KH-like overturns.At the linear stage (figure 10(a)), the impacts of the density perturbation on the mean are barely observable and the interface is thin and flat.Later, the growth of LWI raise and drop on the left-and right-hand sides of  = 40, creating a large-scale jump (figure 10(b)).Meanwhile,  ′ becomes localized at  = 40.Further amplification of the unstable mode breaks down the jump, generating a chaotic region at  = 700.At this stage, the instability is no longer linear, as demonstrated in §4.2.At  = 760, a series of short overturns resembling to KH billows are formed inside and at two sides of the chaotic region at  = 40.These waves propagate away from the chaotic region and may eventually lead to (two-dimensional) turbulence.
Figure 11 shows snapshots of the flow fields at three stages of VLWI-DS.From  = 1500 (panel a) to  = 2500 (panel b), VLWI-DS amplifies, while light-blue (e.g., upper layer: 150 <  < 200) and light-red (e.g., bottom layer: 250 <  < 300) regions become distinguishable, indicating the enhancement of mixing in these regions which acts to dissipate the energy  Finally, in the VLWI-US case (figure 12), the amplification of the stationary waves directly induces localized KH-like billows characterised by strong fluctuations of  ′ in panel (c) at the interface without first breaking down as in LWI.This behaviour may be due to the slower growth rate, which prevents the formation of a distinct nonlinear 'jump' observed in LWI.
As discussed in §4.2, the evolution of these long wave families eventually lead to intense bursting processes that form strong small-scale KH-like overturns.These overturns are responsible for dissipating the kinetic energy injected by a positive slope that cannot be completely balanced by the dissipation of long waves.Note that in all the long waves cases, a short-wave instability (KHI and HWI) does not initially exist according to the LSA in §3.1.In the next section, we study how these short-wave KH-like overturns are induced by nonlinear long waves.
4.4.Breakdown mechanism In §3.1, we showed that the KHI can only occur when   ⪅ 0.25.In the new long wave cases considered in this study,   ≫ 1, hence KHI cannot be triggered.Instead, KH-like overturns are formed by the nonlinear evolution of these long waves.
To understand the cause of the formation of KH-like overturns, we computed the gradient Richardson number   at the total density interface  = 0, which is defined as where we recall that  ≡ 1/4 (see § 2.1).Note that the field   (, ) is based on the total velocity  and density , and thus differs from the constant   defined in (3.1), which is based on the initial base flow profiles  and .
In figure 13, we show the  −  diagrams of   , with  2 contour lines superimposed.In general, as the unstable long waves grow, the density gradient can decrease due to diffusivity which increases the mixing layer.Meanwhile, the velocity gradient can increase as the perturbations are amplified.These results in a decreasing interfacial   , until it reaches below 0.25 (shades of blue), with which small-scale structures are associated.
For each individual long-wave family, the process is slightly different.LWI (figure 13(a)) has two stages of nonlinear breakdown.The first stage appears at  = 680 when a 'jump' is formed at  = 45.This jump changes the density interface and generates two low   regions separated by a high   region.In these low   regions,   continuously decreases due to the amplification of long waves and eventually reduces below 0.25, which potentially allows the growth of the secondary KHI in these regions.Finally, overturns are formed in these regions, leading to the second stage of nonlinear breakdown.Similarly, the amplification of VLWI-DS (figure 13(b)) also causes a low   region that travels along with the waves.As soon as   < 0.25, intense overturns are formed in this region and later contaminate the entire duct.For VLWI-US, the overturns are first formed at the edges of the low   region (170 <  < 330).The close relation between the low   region and the onset of nonlinear short waves strongly suggests that the overturns are a consequence of the decreasing of local   caused by the nonlinear evolution of initially long waves.
From an energy budget perspective, the formation of short-wave overturns in these longwave simulations allows for more efficient dissipation of kinetic energy fed by external forces.In figure 14(a), we illustrate the pathways of turbulent kinetic energy  ′ as given schematically by ) where, , , , and Φ  ′ represent the production, dissipation, buoyancy flux, and transport terms of  ′ .The reader may refer to (Caulfield 2021;Lefauve & Linden 2022) for the definition and a more comprehensive discussion of the kinetic budget of stratified shear flows.When   ≫ 0.25, initially only the long waves are allowed to grow in the flow with strong stratification, gaining energy from the mean flow through production and buoyancy terms and losing it through dissipation.As the long waves are amplified, local shear is created and amplified by the growing velocity perturbations, leading to the decrease of local   .As   ≲ 0.25, the necessary condition for the growth of short waves (mostly KHI here) is satisfied.The short waves then grow, extracting  ′ from the long waves and dissipating it to internal energy.This opens a new energy pathway that allows flows with strong stratification (large   ) to dissipate energy by creating small-scale (turbulent) structures.When   ≪ 0.25 (figure 14(b)), long waves can coexist with short waves (e.g.case IV in figure 3) and may contribute to the energy dissipation.But they are often significantly weaker than short waves since short waves tend to have a faster growth rate.Meanwhile, the short wave directly gains most of the kinetic energy from the mean flow and converts it to internal energy.We also note that turbulence created by these unstable waves can also induce irreversible mixing which, in return, contributes to the production of internal energy.

Conclusions
In this paper, we examined the effects of longitudinal gravitational forces on the stability of two-layer stratified exchange flows by conducting linear stability analyses and nonlinear forced DNS in a sloping channel with solid top and bottom boundaries.In addition to the well-known Holmboe and Kelvin-Helmholtz instabilities, we revealed the existence of three new families of long-wave instabilities subject to non-zero gravitational forces ( ≠ 0): • Long-wave instability (LWI), with wavelengths of the order 10 − 100 channel depths (wave number  =  (10 −2 − (10 −1 )) and a near-zero wave speed; • Downslope very-long-wave instability (VLWI-DS), with wavelengths of the order 100− 1000 channel depths (wave number  =  (10 −3 ) ∼  (10 −2 )), a non-zero wave speed, and complex conjugate eigenmodes implying travelling waves; • Upslope very-long-wave instability (VLWI-US), with wavelengths ≧ 100 channel depths (wave number  <  (10 −2 )) and a near-zero wave speed.
The LWI and VLWI-DS exist at a positive (favourable) slope where the along-slope component of gravity reinforces the pressure gradient, while VLWI-US emerges under a negative (adverse) slope condition.Interestingly, their onset is largely independent of the base flow speed.As a result, they can be triggered even at very high background gradient Richardson numbers   ≫ 1, and induce chaos and sustain (two-dimensional) turbulence and mixing in strongly stratified fluids.In a weakly stratified flow (low   ), they can also coexist with short-wave instabilities (KHI and HWI), but they generally have a lower growth rate.The short-wave HWI and KHI also exhibit interesting features under non-zero slopes.Increasing  tends to suppress the HWI regime while enhancing the KHI.Moreover, the neutral boundary of KHI increases linearly from   = 0.15 at  = −10 • to 0.25 at  = 10 • .
The long-wave families appear under broad flow conditions.To explore their dependence on flow parameters, we varied the Reynolds number Re, Prandtl number Pr, the base flow and boundary conditions.While increasing the Re does not significantly affect KHI, it does enhance the other instabilities.The range of HWI expands to larger   , while the range of the long-wave instabilities approaches  = 0. Therefore, it can be anticipated that as Re → ∞ (as is often the case in natural flows), the critical slope required to trigger these long instabilities approaches zero || → 0. Increasing Pr, or equivalently, decreasing the thickness of the density interface of the base flow, cause the unstable range of HWI to expand towards larger  and smaller   .Meanwhile, the range of the long-wave instabilities slowly approaches  = 0, indicating that these long waves can exist in both water (with Pr ranging from 7 to 700) and air (Pr ≈ 1).It should be noted that these instabilities are not limited to the sine-like base state and the no-slip boundary conditions used in this study.Instead, they can be triggered by, e.g., a tanh-shaped velocity base and free-slip (but impenetrable) velocity boundary conditions.
Finally, we studied the nonlinear evolution of the different instabilities in the inclined channel and their connections to turbulence using a two-dimensional forced DNS that maintains the base states.For all of the long-wave instabilities, the evolution eventually led to a nonlinear bursting process with significant small-scale secondary KH-like overturns and mixing.Specifically, the LWI exhibit two nonlinear stages where an initial breakdown of the long waves is followed by a secondary bursting process, creating multiple intense KH-like overturns.For VLWI-DS and VLWI-US, the long waves do not break down.Instead, they directly alter the base states and induce localized small-scale overturns.
The evolution of these long instabilities results in a decrease in the density gradient and an increase in the shear, which in turn reduces the local gradient Richardson number   .Our analysis reveals that the appearance of KH-like overturns is highly correlated with a local low   , which approaches the critical threshold of 0.25 (below which we find the KHI), substantiating the emergence of localised KH-like overturns.From a turbulent kinetic energy budget perspective, a new energy pathway allows the transfer of kinetic energy from the mean flow to the long waves (linearly) and then to the short waves (nonlinearly), eventually leading to the dissipation of turbulent kinetic energy, under conditions where short waves are linearly stable.
The circumstances under which turbulence can persist in strongly stratified flows remains a fascinating debate within the community (Caulfield 2021).We demonstrated that weakly unstable (very) long waves may trigger turbulence and mixing after long periods of time, even under initially very strongly stratified conditions (  ≫ 1).These results have particular relevance for high−Re flows in rivers (Yoshida et al. 1998) and straits (Gregg & Özsoy 2002), or any natural flow having even very shallow slopes  ≈ 0. A quantitative investigation of the turbulent transition and mixing associated with these long waves would require threedimensional direct numerical simulations, an endeavour left for future work.

Figure 1 :
Figure 1: (a) Schematic of the two-dimensional shear flow in a stratified channel inclined at an angle , and (b) base velocity  () and density () profiles computed from (2.10) and (2.11).

Figure 2 :
Figure 2: Parameter space projections of the fastest growing mode: (a) the growth rate   (colours) and wave frequency   (lines) and (b) the schematics of the   −  parameter space; (c) the growth rate   and wave frequency   and (d) the schematics of   −  parameter space.Markers represent the five cases I, ..., V in Table 1 for which the fastest growing mode is calculated.Black solid lines are the natural convective Thorpe base state (i.e.F = 0), and the horizontal dotted lines in (a) and (b) correspond to   = 0.25.

Figure 3 :
Figure 3: Dispersion relations for typical cases: (a) positive growth rate   versus wave number  and (b) positive growth rate   versus wave frequency   .Markers correspond to cases I, II, III, IV, and V from figure 2.

Figure 6 :
Figure 6: Effect of Reynolds number: fastest growing mode projected onto Ri b −  space for (a) Re = 650, and (b) Re = 5000.

Figure 13 :
Figure 13: Spatial-temporal diagrams of the gradient Richardson number   at the density interface in forced DNS for (a) LWI, (b) VLWI-DS, and (c) VLWI-US.The color maps show the values of   , while the lines represent ⟨ 2 ⟩  .

Figure 14 :
Figure 14: Pathways of turbulent kinetic energy in sloping exchange flows under (a) strong stratification   ≫ 0.25, where only very long waves are unstable, (b) weaker stratification   < 0.25 where both short and long waves coexist.

Table 1 :
Numerical parameter values used for the DNS runs.