Stratified inclined duct: direct numerical simulations

Abstract

The stratified inclined duct (SID) experiment consists of a zero-net-volume exchange flow in a long tilted rectangular duct, which allows the study of realistic stratified shear flows with sustained internal forcing. We present the first three-dimensional direct numerical simulations (DNS) of SID to explore the transitions between increasingly turbulent flow regimes first described by Meyer & Linden (J. Fluid Mech., vol. 753, 2014, pp. 242–253). We develop a numerical set-up that faithfully reproduces the experiments and sustains the flow for arbitrarily long times at minimal computational cost. We recover the four qualitative flow regimes found experimentally in the same regions of parameter space: laminar flow, waves, intermittent turbulence and fully developed turbulence. We find good qualitative and quantitative agreement between DNS and experiments and highlight the added value of DNS to complement experimental diagnostics and increase our understanding of the transition to turbulence, both temporally (laminar/turbulent cycles) and parametrically (as the tilt angle of the duct and the Reynolds number are increased). These results demonstrate that numerical studies of SID – and deeper integration between simulations and experiments – have the potential to lead to a better understanding of stratified turbulence.

1. Introduction

Large-scale fluid motions in the ocean are almost always stably stratified in density due to differences in temperature and/or salinity at different depths. The transport of momentum and buoyancy by turbulence plays an important role in setting the large-scale structure and circulation of the ocean, with implications for the global climate. Consequently, the influence of stable stratification on turbulence and the resulting mixing has attracted much attention (Linden 1979; Riley & Lelong 2000; Gregg et al. 2018; Caulfield 2020, 2021; Dauxois et al. 2021).

Sustained stratified shear-driven flows are a particularly interesting and relevant class of flows to study this problem, since turbulence is produced internally within the flow by drawing energy from the background shear, and because turbulence persists over sufficiently long periods of time to allow for a statistically steady dissipative equilibrium. The stratified inclined duct (SID) experiment was developed to study these flows in a controlled laboratory environment (Meyer & Linden 2014). It establishes a two-layer exchange flow through a long, rectangular and slightly inclined duct connecting two large reservoirs containing fluids of different densities. The SID experiment revealed that the flow regime within the duct could be tuned by adjusting the tilt angle $\theta$ of the duct with respect to the horizontal, and/or the Reynolds number ${Re}$ based on the initial density difference and the height of the duct. The flow regimes are (ordered by increasing $\theta \,{Re}$): laminar two-layer flow, interfacial waves, intermittent turbulence with increased interfacial mixing and eventually full turbulence with significant mixing. Much insight has already been gained through experimental studies of these regimes and of their transitions (Meyer & Linden 2014; Lefauve & Linden 2020a, 2022a), of their energetics and mixing properties (Lefauve, Partridge & Linden 2019; Lefauve & Linden 2022b) and of their respective coherent structures (Lefauve et al. 2018; Jiang et al. 2022).

Despite vast technological improvements yielding unprecedented time-resolved, volumetric velocity and density data (Partridge, Lefauve & Dalziel 2019), experimental limitations remain. The SID experimental data do not yet cover the full length of the duct, do not yet achieve the spatial resolution required to fully quantify energy dissipation and mixing, and are not yet as instantaneous and accurate as we would ideally like (due to their reconstruction of volume by successive scanning of planes). In this paper, we present the first direct numerical simulations (DNS) of SID to help overcome these limitations and integrate experiments and simulations more deeply.

Previous DNS of stratified shear flows have considered more idealised problems, typically without any forcing to sustain the flow (e.g. Salehipour, Peltier & Caulfield 2018; Watanabe et al. 2019). The boundary conditions are usually idealised too, being typically periodic for velocity and density in the streamwise and spanwise directions. By contrast, experiments have revealed that the specific ‘natural’ forcing mechanisms in SID flows (a streamwise hydrostatic pressure gradient and the tilt angle $\theta$) and the lack of periodicity in the streamwise direction (i.e. the presence of reservoirs) are essential features that need to be modelled accurately in order to understand this canonical flow. For example, these features are thought to be closely linked to the notion of ‘hydraulic control’ of the exchange flow at high enough values of $\theta \,{Re}$, and to the ensuing transition to turbulence (Meyer & Linden 2014; Lefauve et al. 2019; Lefauve & Linden 2020a). The no-slip boundary conditions at the duct walls have also been shown to be important to the structures of instabilities (Lefauve et al. 2018; Ducimetière et al. 2021).

In § 2 we explain how we overcome the challenges of developing faithful DNS of SID. In particular, we discuss how we modelled the reservoirs with minimal computational cost, and how we handled technically challenging boundary conditions. In § 3 we validate this new DNS methodology by comparing different reservoir geometries and forcing with fully resolved computations that capture the reservoirs explicitly. In § 4 we describe the flow regimes and compare the DNS with experimental data, first from regime diagrams (i.e. the map of the observed qualitative flow regimes in the two-dimensional parameter space $\theta$–${Re}$) and then from shadowgraph visualisations of the density interfaces. Then in § 5 we describe further quantitative diagnostics from our DNS, generally inaccessible to experiments, and highlight their added value. These include the gradient Richardson number, the turbulent kinetic energy (TKE) and pressure fields along the entire length of the duct and the turbulent energy fluxes. Finally, in § 6 we conclude by summarising our results and outlook.

2. Methodology

2.1. Governing equations

Our simulation geometry in non-dimensional units is shown in figure 1(a,b). It replicates the experimental geometry (see e.g. figure 1 of Lefauve et al. (2019) in dimensional units), which consists of a duct of square cross-section with internal height $H$, width $W$ and length $L$ connecting two large reservoirs with fluids at densities $\rho _0\pm \Delta \rho /2$ (white and blue shaded areas in figure 1a). To match previous experimental studies of SID, we non-dimensionalise all lengths by the duct half-height $H/2$, making the duct non-dimensional length, height and width $2A\times 2B\times 2$, respectively, where $A\equiv L/H$ and $B=W/H$ are the streamwise and spanwise aspect ratios, respectively. We also non-dimensionalise (i) the velocities by the fixed buoyancy velocity scale $\Delta U/2 \equiv \sqrt {g^\prime H}$ (where $g^\prime =g\Delta \rho /\rho _0$ is the reduced gravity and $\rho _0$ is the reference density); (ii) the time by the advective time unit (ATU) $H/\Delta U$; (iii) the density variations around the reference $\rho _0$ by $\Delta \rho /2$; and (iv) the pressure by $\rho _0(\Delta U/2)^2$. Note that the $x$ axis (the streamwise direction) is aligned along the duct, whereas gravity points downwards at an angle $\theta$ from the $-z$ axis (the vertical direction in the frame of the duct); hence in these duct coordinates $\boldsymbol {g}=g {\hat {\boldsymbol {g}}}=g(\sin \theta, 0, -\cos \theta )$.

Figure 1. Schematics of SID geometry in non-dimensional units. (a) Overview of the rectangular simulation domain of dimensions $L_x,L_y,L_z$ within which immersed boundaries create a square duct of dimensions $2A\times 2\times 2$. (b) Detail of the duct geometry and coordinate system. (c) Shape of the different reservoirs considered in this paper (benchmark, AR, BR and SR), with a total domain length $L_x=2(A+L^r_x)$. All numerical parameters are summarised in table 1.

Table 1. Summary of the DNS. From left to right: reservoir geometry (case) as shown in figure 1(c); Reynolds number; tilt angle; duct streamwise aspect ratio; duct spanwise aspect ratio; reservoir size; grid size of the entire computational domain; and forcing parameters: streamwise length of forcing ($l_f$) and time scales controlling momentum ($\eta _u$) and density forcing ($\eta _\rho$) in (2.5)–(2.6). Bold font and superscripts denote the most used DNS.

The resulting non-dimensional governing equations for our DNS are the Navier–Stokes equations under the Boussinesq approximation:

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

(2.2)$$\begin{gather}\frac{{\rm D} \boldsymbol{u}}{{\rm D} t} ={-} \boldsymbol{\nabla}p + \frac{1}{Re} \nabla^{2}\boldsymbol{u}+ {Ri} \, \rho \, {\hat{\boldsymbol{g}}} - \boldsymbol{F}_u, \end{gather}$$

(2.3)$$\begin{gather}\frac{{\rm D} \rho}{{\rm D} t} = \frac{1}{{Re \,Pr}} \nabla^{2}\rho - F_\rho, \end{gather}$$

where the material derivative is ${\rm D}/{\rm D}t \equiv \partial _t+\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla }$, the velocity is $\boldsymbol {u}=(u,v,w)$ in the non-dimensional coordinate system $\boldsymbol {x}=(x,y,z)$ aligned with the duct, the non-dimensional pressure is $p$ and the non-dimensional density variation around the mean is $\rho$ (bounded between $-1$ and $1$). The forcing terms $\boldsymbol {F}_u$ and $F_\rho$ used to maintain the quasi-steady exchange flows are described in § 2.2.

The non-dimensional Reynolds, Richardson and Prandtl numbers are related to the dimensional experimental parameters as follows:

(2.4a–c)\begin{equation} {Re}\equiv \frac{\dfrac{\Delta U}{2}\dfrac{H}{2}}{\nu} \equiv \dfrac{\sqrt{g'H}H}{2\nu}, \quad {Ri}\equiv \frac{\dfrac{g}{\rho_0}\dfrac{\Delta\rho}{2}\dfrac{H}{2}}{\left(\dfrac{\Delta U}{2}\right)^2}\equiv \frac{1}{4}, \quad {Pr}\equiv \frac{\nu}{\kappa} \equiv 7, \end{equation}

where $\nu$ is the kinematic viscosity and $\kappa$ is the mass diffusivity. Previous studies of SID showed that the streamwise velocity scales with $\Delta U/2$, motivating this definition of Reynolds number. The Richardson number is always equal to 1/4 due to the definition of $\Delta U$. The Prandtl number in all simulations was set to ${Pr}=7$, approximately representative of temperature stratification in water at room temperature. For a given duct and reservoir geometry, there are two remaining free non-dimensional parameters: the tilt angle $\theta$ and the Reynolds number ${Re}$ (based on the driving density difference $\Delta \rho$).

2.2. Artificial restoring of the exchange flow

The exchange flow in the duct is driven by the hydrostatic longitudinal pressure gradient and by the longitudinal gravitational acceleration $g\sin \theta$. In the context of a two-layer flow, the along-duct component of gravity accelerates the heavier layer rightwards (downhill) and the lighter layer leftwards (uphill). The role of the hydrostatic pressure gradient turns out to be more intricate and is examined in § 5.3. In the experiments, the flow inside the duct is sustained over long time periods (typically several hundred ATUs) until the discharged fluids accumulated in the large reservoirs have reached the level of the duct. Simulating such large reservoirs would be prohibitively expensive. In the simulations, we use smaller reservoirs and add ad hoc forcing terms $\boldsymbol {F}_u, F_\rho$ in the momentum and buoyancy equations (2.2) and (2.3), respectively:

(2.5)$$\begin{gather} \boldsymbol{F}_u \equiv F_u \boldsymbol{u} \equiv \left[ \dfrac{1-\tanh\left(\dfrac{2}{\varDelta}\left(x+\dfrac{L_x-l_{f}}{2}\right)\right)}{\eta_u}+ \dfrac{1+\tanh\left(\dfrac{2}{\varDelta}\left(x-\dfrac{L_x-l_{f}}{2}\right)\right)}{\eta_u} \right]\boldsymbol{u} , \end{gather}$$

(2.6)$$\begin{gather}F_\rho\equiv \dfrac{1-\tanh\left(\dfrac{2}{\varDelta}\left(x+\dfrac{L_x-l_{f}}{2}\right)\right)} {\eta_\rho}(\rho-1)+\dfrac{1+\tanh\left(\dfrac{2}{\varDelta}\left(x-\dfrac{L_x-l_{f}}{2}\right)\right)}{\eta_\rho}(\rho+1) , \end{gather}$$

where $l_{f}$ is the streamwise length of influence of the forcing and $\varDelta =2 l_f/L_f$ (with a fixed $L_f\equiv 8$) defines the steepness of the transition from the forced to the unforced regions. The density forcing term restores the density of the fluid in the reservoir to the prescribed value (i.e. $\pm 1$), and the momentum forcing term acts to dampen motion in the reservoir. The time scales $\eta _u$ and $\eta _\rho$ control the momentum and density forcing terms, respectively. Compromise values of these time scales must be found, as large values are too slow to sufficiently damp reservoir motion and restore density, while small values are too fast and overreact, threatening numerical stability. These parameters were optimised with the size of the reservoirs in order to minimise their influence on the large-scale flow in the duct compared with the benchmark cases with large reservoirs and without forcing (see § 2.6). Tests revealed little variation in the range $l_{f}\in [0.3L_x^r ,0.7L_x^r]$; therefore we set ${l_f=0.5L_x^r}$, confining the forcing region to half the reservoir (greyed out in figure 1a). The time scales $\eta _u$ and $\eta _\rho$ should then be smaller than the times for a discharging flow (with non-dimensional speed 1) to pass through the forcing region, i.e. $\approx l_f$. Practically, we set $2.5 \lesssim \eta _u \lesssim 5$ and $0.1 \lesssim \eta _\rho \lesssim 0.5$, depending on $l_f$.

Physically, $\boldsymbol {F}_u$ decelerates the fluid entering the reservoir until it comes to rest, and $F_\rho$ ensures that the density of this fluid matches that of the reservoir before it re-enters the duct. This forcing thus effectively mimics the action of infinitely large reservoirs with a finite-sized, computationally feasible domain.

2.3. Solver

The DNS were performed with the open-source solver Xcompact3D (Bartholomew et al. 2020), which uses fourth-order and sixth-order compact finite-difference schemes for the first and second spatial derivatives, respectively, and a third-order Adams–Bashforth scheme (Peyret 2002; Zhu & Xi 2020) for the time integration with a time step ${\delta _t=0.001}$. The pressure field is obtained from a conventional Poisson equation (based on applying a divergence operator on (2.2) and employing continuity (2.1)), which is then solved numerically using the fast Fourier transform with modified wavenumbers. For more details about the core of the code (Incompact3D), see Laizet & Lamballais (2009) and Laizet & Li (2011), and for the application of Xcompact3D to stratified turbulent flows, see Frantz et al. (2021). We modified Xcompact3D to include the forcing terms $\boldsymbol {F}_u,F_\rho$ discussed above.

2.4. Domain and boundary conditions

The computational domain had dimensions $L_x$, $L_y=2$ and $L_z$ along $x$, $y$ and $z$, respectively (see figure 1a,b). On the boundaries of this domain, we applied a no-slip condition for $\boldsymbol {u}$ and a no-flux condition for $\rho$ as in Laizet & Lamballais (2009). To represent the duct and reservoir geometry within this computational domain, we applied the immersed boundary method (IBM) in Xcompact3D to the yellow-shaded region in figure 1(a).

The IBM treatment of $\boldsymbol {u}$ (no slip) uses a direct forcing method described in Mohd-Yusof (1997) and specifically for Xcompact3D in Laizet & Lamballais (2009) and Gautier, Laizet & Lamballais (2014), which imposes $\boldsymbol {u}=\boldsymbol {0}$ in the solid regions. The pressure $p$ in the solid region is treated by reducing the Poisson equation to a Laplace equation (Laizet & Lamballais 2009). The IBM allows relatively simple implementation of complex geometries in scalable codes such as Xcompact3D that are built upon Cartesian coordinates and rectangular computational domains.

The IBM treatment of $\rho$ (no flux) required a slightly different approach to minimise the modifications of Xcompact3D and maintain the consistency between $\boldsymbol {u}$ and $\rho$. A $\tanh$ function was used to reconstruct points inside the solid region:

(2.7)\begin{equation} \rho_i = \frac{1}{2}(\rho_{l}+\rho_{r})+\frac{1}{2}(\rho_{r}+\rho_{l})\tanh \left[\frac{L_{w}}{\xi_{r}-\xi_{l}}\left(\xi_i-\frac{1}{2}(\xi_{r}+\xi_{l})\right)\right], \end{equation}

as shown in figure 2. Here $\xi$ is the wall-normal coordinate (horizontal or vertical), the subscript $i$ is the index of grid points and $r$ and $l$ are the right and left grid points (in the fluid region), respectively, adjacent to the solid walls. The $\tanh$ function ensures a zero-flux boundary condition at the wall while maintaining a smooth change of the density through the solid region. A similar approach using a polynomial reconstruction has been used to treat the Dirichlet and Neumann boundary conditions in Gautier et al. (2014) and Frantz et al. (2021). The length scale $L_w=10$ was chosen to ensure a smooth change of density inside the solid region while maintaining exponentially small density flux at SID walls.

Figure 2. Schematic of the IBM to implement the no-flux boundary condition on density $\rho$ at fluid–solid boundaries in the $x$ or $z$ direction. The curved red line represents the fictitious density profile across the solid region.

2.5. Initial conditions

All simulations were initialised at $t=0$ with a density $\rho = \tanh (x/L_I)$ (where $L_I=0.1$) at the centre of the duct, simulating ‘lock exchange’ conditions with a sharp but continuous change from densest fluid on the left-hand side to lightest fluid on the right-hand side. A zero-mean uniformly distributed random noise with non-dimensional amplitude $\varsigma =0.5$ was applied to the initial velocity $\boldsymbol {u}_n$ to break symmetry and initiate instabilities inside the duct. Note that smaller perturbation amplitudes (e.g. $\varsigma =0.005$) can be applied, but we verified that doing so did not influence the main features of the flow (see supplementary material S1 available at https://doi.org/10.1017/jfm.2023.502).

Shortly after $t=0$, a gravity current formed at the centre of the duct ($x=0$) and propagated in both directions towards the ends of the duct. After a typical duct transit time of order $t\approx A$ (transiting at non-dimensional velocity $\approx$1 over a non-dimensional length $A$), the exchange flow was established.

2.6. Parameters, duct and reservoir geometries

In order to investigate the various flow regimes we varied the Reynolds number ${Re}$ in the range 400–1250 and the duct tilt angle $\theta$ in the range $1^\circ \unicode{x2013}10^\circ$. As mentioned above, the Prandtl number $Pr=7$ and Richardson number ${Ri}=1/4$ were fixed. The suite of DNS is summarised in table 1.

Most DNS were run with a long duct of streamwise and spanwise aspect ratios $A=30$ and $B=1$, respectively, for direct comparison with the experiments in Lefauve & Linden (2020a) (the ‘mini SID Temperature’ dataset abbreviated ‘mSIDT’). However, two DNS (cases ‘BRw’ in table 1) were run in a longer and wider duct at $A=44$, $B=2$ to compare with a new experimental set-up.

To validate the performance of our forcing to sustain a realistic exchange flow, we ran a benchmark DNS without forcing ($F_u=F_\rho =0$) but with large reservoirs ($L_x^r\times L_z=30\times 8$). This benchmark had a combined reservoir volume of four times that of the duct ($60\times 8/(60 \times 2)=4$), which is sufficient for our validation but still much smaller than that of the experiments (volume ratio $\approx$30). We show in § 3 that the different reservoirs do not seem to influence the flow statistics within SID. This is expected from the knowledge that the flow in SID is hydraulically controlled (Meyer & Linden 2014), i.e. that information from the reservoirs cannot travel into the duct because of ‘control’ regions at the inlet and outlet, where the convective flow speed is faster than interfacial waves (Lawrence 1990). This conveniently ensures that different reservoir geometries and conditions do not influence the flow within the duct, as long as unmixed and quiescent fluids are available at either end of the duct.

All the other DNS had non-zero forcing and smaller, more computationally affordable reservoirs. To test the impact of reservoir size, we used the three following reservoirs sketched in figure 1(c): the A-reservoir (‘AR’) of dimensions $L_x^r\times L_z=10\times 8$, which is a third of the length of benchmark but equally tall; the B-reservoir (‘BR’) of $L_x^r\times L_z=5\times 4$, which is half the length and half the height of the A-reservoir; and finally, the smallest S-reservoir (‘SR’) of $L_x^r\times L_z=10\times 2$. Note that (unlike the experiments) all reservoirs have the same spanwise width as the duct $L_y$. The set of forcing parameters $(l_f, \nu _u,\nu _\rho )$ for $\boldsymbol {F}_u,F_\rho$ that we found to have minimal impacts on the duct for each case are also listed in table 1.

The bold font for the nine cases at ${Re}=650$ and $1000$ highlight the DNS that we analyse in more detail in this paper, with the superscripts giving their shorthand names (B2, B5, B6, B8, B10, S6, S8, W3 and W5). The full flow data of the five cases B2-B10 can be downloaded from the repository (Zhu et al. 2023). The other DNS were used for validation and for plotting the regime diagrams in the $(\theta,{Re})$ plane.

Finally, we adopt a uniform grid size $N_x \times N_y\times N_z$ for the entire domain $L_x\times L_y\times L_z$, which has the advantage of helping maintain numerical stability near the immersed boundaries. The grid size was small enough to capture the Kolmogorov turbulent length scale, and $2\unicode{x2013}3$ times the Batchelor length scale in our most turbulent dataset B10 (discussed in more detail in § 5.4), ensuring adequate resolution of the kinetic and scalar energy spectra.

3. Validation

In figure 3 we assess the ability of the forcing introduced in (2.5)–(2.6) to sustain the exchange flow by comparing, for the B-reservoir, a standard DNS with forcing (‘forced’) and a DNS without forcing (‘unforced’), i.e. $F_u=F_\rho =0$.

Figure 3. Demonstration of the effects of the forcing terms $F_u,F_\rho$ in finite-sized reservoirs (here in a B-reservoir, ‘BR’) at $({Re},\theta )=(650,4^\circ )$. (a) Time series of the volume and mass flux. Instantaneous mid-duct slices of $\rho (x,y=0,z,t=100)$ for (b) forced and (c) unforced DNS, showing only the rightmost quarter of the duct and the right-hand reservoir.

Figure 3(a) shows the time series of the volume flux $Q$ and the mass flux $Q_m$ defined as

(3.1)$$\begin{gather} Q(t) \equiv \langle \vert u \vert\rangle_\mathcal{V}, \end{gather}$$

(3.2)$$\begin{gather}Q_m(t) \equiv \langle \rho u \rangle_\mathcal{V}, \end{gather}$$

where $\langle \cdot \rangle _\mathcal {V} \equiv (1/8A)\int _{-1}^1\int _{-1}^1 \int _{-A}^{A}\, {{\rm d}x} \, {{\rm d}y} \, {\rm d}z$ denotes an average over the entire volume of the duct. Note that since $\vert \rho \vert \leqslant 1$, by definition $Q_m\leqslant Q$. A more diffuse interface and turbulent mixing can cause $Q_m$ to be significantly lower than $Q$. The values of $Q(t)$ (solid lines) and $Q_m(t)$ (dashed lines) in the forced (red) and unforced (green) DNS are identical in the initial stage of accelerating gravity current ($0 < t \lesssim 60$). They remain equal until the exchange flow approaches a steady state at $Q\approx 0.5$ and $Q_m\approx 0.45$ ($60\lesssim t \lesssim 100$). However, from $t\approx 100$, the unforced time series drops sharply, signalling that the flow slows down (see $Q(t)$) and becomes overall more mixed inside the duct (as $Q_m(t)$ decays faster than $Q(t)$). By contrast, the forced time series remains steady until the end ($t\approx 160$) of the simulation.

Figures 3(b) and 3(c) show $x$–$z$ slices of the density field on the rightmost quarter of the computational domain at $t=100$ for the forced DNS and unforced DNS, respectively. In the unforced DNS, the dense, right-flowing bottom layer (in blue) has filled over half of the B-reservoir. The large kinetic energy of this layer has led to mixing inside the reservoir. This dense fluid contaminates the exchange flow as it is entrained back into the duct by the left-flowing buoyant layer (in red). In the forced DNS, this does not happen; the outflowing layer is slowed down and its density is gradually converted to that of the inflowing fluid. This allows an infinitely-long quasi-steady exchange flow to be maintained inside the duct.

In figure 4 we compare the statistics of the established exchange flow in the benchmark case (very large reservoirs, unforced), and in progressively smaller, but forced, reservoirs: BR and SR. We compare two different flows: a laminar regime at $({Re},\theta )=(400,5)$ (red, blue and green curves) and a wave regime at $({Re},\theta )=(650,6)$ (purple, pink and cyan curves). Figure 4(a) shows $\langle u \rangle (z)$ (where $\langle \cdot \rangle \equiv \langle \cdot \rangle _{x,y,t}$ is the average over the entire duct length, width and time series), figure 4(b) shows $\langle \rho \rangle (z)$ and figure 4(c) shows the time series of the total kinetic energy $\langle \bar {k} \rangle _\mathcal {V}(t)$ (where $\bar {k}\equiv |\boldsymbol {u}|^2/2$).

Figure 4. Comparison of the effects of reservoir sizes on the (a) streamwise velocity and (b) density profiles, and (c) kinetic energy time series in both the laminar and wave regimes.

Comparing the benchmark, BR and SR cases, we find excellent agreement between the vertical profiles and the time series of kinetic energy. Minor temporal discrepancies in the wave regime after $t\gtrsim 150$ may be caused by variations in the initial random noise, but have negligible influence on the flow statistics and dynamics. Overall, our forcing method faithfully models the effects of reservoirs as far as the flow inside the duct is concerned, even with very small reservoirs. We provide additional evidence in supplementary material S1 that the key flow dynamics in all flow regimes in SID are largely independent of the reservoirs. We compare spatio-temporal diagrams of the TKE for $({Re},\theta )=(650,6^\circ )$ for the benchmark, AR, BR, as well as BR with a reservoir wider than the duct and show that details of wave motion and occasional turbulence over 200 ATUs do not vary more than they would under different initial noise conditions. The spatio-temporal diagrams of TKE for $({Re},\theta )=(650,8^\circ )$ for the benchmark and BR also agree well in terms of turbulent intermittency and wave propagation.

Our typical computation at ${Re}=650$ required $45\times 10^6$ points in the AR, but only $17\times 10^6$ in the BR (a reduction of 60 %). This explains why, in the following, we use the BR for more detailed analyses requiring longer time series of order $\approx$200 ATU. We use the even more affordable SR more sparingly in this paper, since our main goal is to compare the SR results with the BR to investigate the ability of the SR to reproduce the key flow physics.

4. Comparison between DNS and experiments

4.1. Regimes: observations

Figure 5 shows snapshots of the DNS density field illustrating the different flow regimes. All cases use the B-reservoir, with a duct aspect ratio of $A=30$, as highlighted in table 1 and named B2, B5, B6, B8 and B10. The first four cases B2–B8 were at ${Re}=650$, the last one B10 was at ${Re}=1000$, and the numbers 2, 5, 6, 8, 10 indicate the respective value of $\theta$ in degrees. Slices through the density field at the middle $y=0$ plane (figure 5a–e), and through a cross-sectional $y$–$z$ plane at $x=0$ in the duct are shown. The full temporal evolution of these five cases can be seen in our supplementary movies.

Figure 5. Snapshots of the density field in the mid-plane $y=0$ (a–e) and in the duct cross-section ( f–j) in five representative flows: (a, f) laminar (B2), (b,g) stationary wave (B5), (c,h) travelling wave (B6), (d,i) intermittent turbulence (B8, active phase) and (e,j) fully developed turbulence (B10). These cases are highlighted in bold font in table 1.

We recover, in our DNS, the same four key flow regimes (laminar, wave, intermittently turbulent and fully turbulent) identified in experimental studies of SID, in particular Meyer & Linden (2014), Lefauve (2018), Lefauve et al. (2019) and Lefauve & Linden (2020a), which is the first key result of this paper. Moreover, the DNS allow us, for the first time, to observe three-dimensional instantaneous snapshots along the entire domain, including the duct and the in- and out-flow in the reservoirs, which were not accessible to experiments. We describe each regime in turn.

4.1.1. Laminar regime

First, in B2 (figure 5a, f) we observe a simple laminar (L) flow, which is largely parallel and steady without any observable waves or turbulent fluctuations. Molecular diffusion creates a relatively thin interface of intermediate density (in white). This density interface slopes at an angle since the two counter-flowing layers (in blue or red) get thinner as they accelerate along the duct. This convective acceleration $u\partial _x u$ in each layer is caused by the pressure gradient $-\partial _x p$, gravity ${Ri} \sin \theta \rho$ and opposed by the viscous term ${Re}^{-1}{\nabla }^2 u$.

4.1.2. Wave regime(s)

Second, in B5 and B6 (figure 5b,c,g,h) we observe a wave (W) flow, including large-scale waves with streamwise wavelength of order $O(1)\unicode{x2013}O(A)$ (i.e. $O(1)\unicode{x2013}O(30)$). These waves are triggered by disturbances within the duct; information (waves) does not travel from the reservoirs into the duct. In this regime, small-scale, weakly turbulent structures of limited spatial extent may occasionally be generated by the breakdown of large-scale waves but they always dissipate rapidly. We note that most of the previous SID experiments were done with salt stratification ($Pr\approx 700$), in which case the much thinner density interface supports Holmboe waves; hence these studies called this wave regime the ‘Holmboe’ regime. However, with temperature stratification ($Pr\approx 7$, as simulated here) Lefauve (2018, § 3.4.2) and Lefauve & Linden (2020a) highlighted that Holmboe waves were never found on the thicker interface. At $Pr\approx 7$ they found the same wave regime as that observed here, with interfacial gravity waves on the edges of a thicker density interface.

We find that W flows can feature waves that are either stationary (B5; figure 5b,g) or travelling (B6; figure 5c,h) in the streamwise direction $x$. Increasing $\theta$ tends to first decrease the slope of the interface (relative to the duct) and accelerate the flow until the interface is parallel to the axis of the duct (reducing $u\partial _x u$ and $-\partial _x p$), at which point the gravitational term $Ri \sin \theta \rho$ can no longer be balanced by laminar diffusion alone. This appears to coincide with the creation of a third, partially mixed layer (in light red, white and light blue) that is neutrally buoyant and thus reduces the gravitational forcing $Ri \sin \theta \rho$. This third layer, often located near the centre of the duct (i.e. around $|x|\approx 0$) rather than near the ends ($|x|\approx A$) supports both stationary and travelling interfacial waves. We note that, in contrast to L flow, in these W flows the in-flowing layers (before reaching the ‘wavy’ area in the centre of the duct) are thinner than the out-flowing layers (after going through the ‘wavy’ area), which is reminiscent of an internal hydraulic jump.

Travelling waves (figure 5c) tend to travel along the two density interfaces (between the top and middle layers, and between the middle and bottom layers) in a specific fashion. Left-going waves are most often found on the right quarter of the duct ($x \gtrsim A/2$), travelling towards the centre. Vice versa, right-going waves are most often found on the left quarter of the duct ($x \lesssim -A/2$), travelling towards the centre. Once they reach the central region ($|x| \lesssim A/2$), both types of waves usually end up decaying. This observation suggests that the flow may be ‘supercritical’ outside of the central region, i.e. that information transported by interfacial waves can only propagate in one direction (towards the centre but not towards the ends of the duct).

4.1.3. Intermittently turbulent regime

Third, in B8 (shown in figure 5d,i) we observe an intermittently turbulent (I) flow, which becomes more chaotic and in which patterns of individual waves become indistinguishable. Small-scale turbulent structures (of typical non-dimensional scale $\ll 1$) are generated, often by a breakdown of waves akin to the ‘bursting’ events of turbulent boundary layers (Robinson 1991; Jiménez & Simens 2001; Zhu & Xi 2020). This interfacial turbulence, which persists for much longer times than in the W regime, enhances interfacial mixing and creates a third partially mixed layer (shown in white) over an increasingly long streamwise extent (as compared with the W regime). The interfacial turbulence sometimes extends along the full length of the duct. The combination of the decreasing magnitude of the gravitational forcing $Ri \sin \theta |\rho |$ by the increasingly mixed layer and the increasing smaller-scale viscous dissipation are presumably the key ingredients that keep the flow steady as $\theta$ is increased from $2^\circ,5^\circ,6^\circ$ to $8^\circ$ in B2, B5, B6, B8 (at constant $Re=650$).

The defining characteristic of the I regime is that the turbulence identified by small-scale structures is temporally intermittent; turbulence occasionally decays and the flow ‘relaminarises’ before transitioning to turbulence again. These cycles are described in § 5.2. The time scales associated with the transition to turbulence and its decay, and the advection of perturbations along the length of the duct, occasionally make this turbulence also spatially intermittent in $x$.

4.1.4. Fully turbulent regime

Fourth, in B10 (shown in figure 5e,j) we observe a fully turbulent (T) regime in which turbulence is sustained in time and is more vigorous than in the I regime. Although the intensity of the turbulence can fluctuate in time, the flow in this regime never fully relaminarises. The central partially mixed layer typically covers the entire length of the duct and at least a third of its height.

4.2. Regime diagrams

In figure 6 we map the flow regimes described above in 12 DNS with the AR (figure 6a) and in 15 DNS with the BR (figure 6b) for a range of $\theta$ and $Re$.

Figure 6. Regime transitions in $\theta$–$Re$ parameter space. The large symbols are our DNS data in the (a) A-reservoir and (b) B-reservoir, and the small markers are the temperature-stratified experimental data of Lefauve & Linden (2020a) (see their figure 4e) with matched non-dimensional parameters.

These 27 DNS data points, shown as large symbols, are compared with the 148 experimental data points of Lefauve & Linden (2020a) taken from their figure 4(e) and displayed here as smaller, fainter symbols using the same colour coding for the different flow regimes. The experimental data points were obtained using the same aspect ratios ($A=30, B=1$) and with temperature stratification ($Pr\approx 7$). The regimes were identified by shadowgraph visualisation, often over a small streamwise extent of the duct (their movies can be downloaded from Lefauve & Linden (2020b)). Figure 6 therefore represents the first direct comparison of DNS results with experimental results in SID, with all non-dimensional control parameters matched.

First, we find a general agreement between DNS and experiments in the location of flow regimes in the $(\theta,{Re})$ plane, as evidenced by the fact that most large symbols (DNS) are of the same type as the smaller symbols (experiments). This is the second key result of this paper, because it confirms that our DNS, with small computationally efficient reservoirs, can reproduce the key physics of SID, encapsulated in the flow regimes.

The minor exception to this agreement is found near the L/W transition, where some of our DNS found the W regime whereas the experiments found the L regime. This may be a genuine difference, but we suspect that this may be due to the fact that the weak stationary waves found near the L transition may have been missed in the experiments. This is because the experimental shadowgraphs were visualised over a limited extent of the duct and because low-amplitude waves in low-${Re}$ temperature-stratified flows produce very small changes in refractive index and thus weak shadowgraph signals.

Second, the AR and BR yield consistent results (figure 6a,b), confirming that the smallest ‘true’ reservoir (excluding the SR for now) is indeed sufficient to reproduce the experiments. These results offer strong further support to the preliminary validation of our suite of DNS in § 3.

4.3. Shadowgraphs

We now turn to a side-by-side comparison of shadowgraph visualisations of the flow in DNS and experiments within a particular flow regime.

Experimental shadowgraph movies are obtained by the projection onto a semi-transparent screen of initially parallel light rays that have travelled through the duct along the spanwise $y$ direction. Any variations in the curvature (normal to the rays) of the density field $\rho$ (and hence refractive index field $n$) cause the rays to focus or defocus, varying the intensity that reaches the screen (Weyl 1954). In the limit of weak variations, the intensity of the image formed is (see e.g. Lefauve 2018, § 2.1)

(4.1)\begin{equation} I(x,z,t) = \beta I_0(x,z) \int_{{-}B}^{B} (\partial_{xx}+\partial_{zz})\rho(x,y,z,t) \,{{\rm d}y}. \end{equation}

Here $\beta$ depends on $(\rho _0/n_0)\partial n/\partial \rho$ and the experimental geometry, while $I_0$ is the (approximately) uniform background intensity of the illumination. This field is thus particularly suited to detect density interfaces, and is a simple and efficient proxy to compare the structure of interfacial density waves and small-scale turbulence in DNS and experiments.

Figure 7 compares false colour instantaneous snapshots of $I(x,z)$ in the I regime over a central portion of the duct $|x|<9$. The DNS shadowgraphs reconstructed from the calculated density fields (assuming $\beta I_0 =1$) are shown in figure 7(a–c) (W3 and W5; see table 1), and the matching experimental shadowgraph images of $I/I_0$ are shown in figure 7(d–f), all at $Re=650$ and ${Pr}=7$. We show a single snapshot at $\theta =3^\circ$, at the boundary between the W and I regimes (figure 7a) and two snapshots at $\theta =5^\circ$, well into the I regime, where the flow is in a quiet laminar phase (figure 7b) and in an active turbulent phase (figure 7c). The full temporal evolution of these four shadowgraphs can be found in our supplementary movies.

Figure 7. Snapshots of shadowgraph comparing the normalised intensity $I(x,z)$ in DNS (a–c) to the matching experiments (d–f) in two cases W3 (a,d) and W5 (b,e,c, f). Magnitudes (colour bar limits) are naturally different due to the unknown experimental $\beta$ factor in (4.1). The times at which these snapshots were taken are shown by the vertical lines in the spatio-temporal diagrams of figure 8.

Figure 8. Spatio-temporal diagrams of shadowgraph in W3 (a,c) and W5 (b,d) comparing DNS (a,b) with experiments (c,d). The vertical black solid lines indicate the time of the snapshots in figure 7.

Note that these shadowgraphs were obtained in a new experimental apparatus having a wide duct $B=2$, with a regular straight rectangular section of length $A=40$. In these experiments, the reservoirs were tilted in line with the duct (the whole apparatus pivots about one end), which exactly matches the DNS. The apparatus had trumpet-shaped expansions at either end (over an additional 10 % of its length) for a smoother connection to the reservoirs. While we did not model the trumpet ends in our DNS, we used $B=2$ and the total length $A=44$ to reproduce the geometry as faithfully as possible. Trumpet ends were first used in Meyer & Linden (2014) who reported no visible impact in shadowgraphs when compared to straight ends. With the parameters $A$ and $B$ increased with respect to cases B2–B10, the I regime is found at smaller $\theta$ values than would be expected from figure 6 (see Lefauve & Linden 2020a).

We find good agreement in the structure of interfacial waves, somewhat reminiscent of Kelvin–Helmholtz billows, in DNS and experiments (compare figure 7a,b with figure 7d,e). These waves have higher amplitude than the stationary waves previously found in B5 (see figure 5b) because the flow is more energetic and prone to the growth of stratified shear instabilities at $B=2$ than at $B=1$, due to a weaker influence of the no-slip sidewalls (see Ducimetière et al. 2021, § IIIc). These waves tend to break into weak and short-lived turbulence at $\theta =3^\circ$ (placing it borderline in the I regime), and into stronger and longer-lived turbulence at $\theta =5^\circ$ (placing it well into the I regime). We also find good agreement in the overall appearance of small-scale turbulence in the ‘active’ phase (compare figures 7c and 7f). Active turbulence in the experiment extends slightly closer to the top and bottom boundaries than in the DNS. This may be a result of various factors including the non-zero thermal conductivity of the experimental duct walls, spurious reflections of light and excessive cropping of near-wall regions caused by the difficulty in locating the exact location of the walls in the shadowgraph images.

Figure 8 illustrates these temporal dynamics with the corresponding $z$–$t$ spatio-temporal diagrams in DNS (figure 8a,b) and in experiments (figure 8c,d). We find again good agreement, both in the vertical growth and decay of the waves and in the alternation and approximate period of the quiet and active phases.

These shadowgraphs show that our DNS faithfully reproduce not only the qualitative flow regimes and their distribution in $\theta$–$Re$ space, but also details of their spatial structures and temporal dynamics, which is a third key result of this paper.

5. Added value of DNS

In this section we examine quantitative DNS diagnostics which, because they are difficult or impossible to obtain in experiments, add value to the experimental study of SID.

5.1. Vertical profiles and gradient Richardson number

Figure 9 shows, for the five flow regimes previously shown in figure 5, the $x,y,t$-averaged velocity $\langle u \rangle (z)$, density $\langle \rho \rangle (z)$ and gradient Richardson number ${Ri}_g$ based on the gradients of these mean flows:

(5.1)\begin{equation} {Ri}_g(z)\equiv{-}{Ri}\,\frac{\partial_z \langle \rho \rangle}{(\partial_z \langle u \rangle)^2}. \end{equation}

Figure 9. Vertical profiles of the mean (a) streamwise velocity $\langle u \rangle$, (b) density $\langle \rho \rangle$ and (c) gradient Richardson number $Ri_g$ in the five flows of figure 5. We also include the I and T experimental profiles of Lefauve et al. (2019) at $Pr\approx 700$. The vertical dashed lines in (c) denote $Ri_g=$ 0.1 and 0.25.

Such simultaneous velocity and density diagnostics are available in salt-stratified experiments (at ${Pr}\approx 700$), and we superimpose on figure 9 the mean profiles in the I and T regimes from Lefauve et al. (2019) (their figure 4f,l). However, these diagnostics cannot be accurately obtained in temperature-stratified experiments (at $Pr\approx 7$) to match our DNS for two main reasons. First, the velocity field measurements rely on particle image velocimetry in a refractive index matched fluid, which is impossible without the introduction of another stratifying agent (necessarily having a much smaller diffusivity than temperature). Second, the density field measurements rely on laser-induced fluorescence with a dye having a much smaller diffusivity than temperature, and therefore ‘tagging’ it poorly (temperature-sensitive fluorescent dyes exist but such measurements are more difficult and less accurate).

In figure 9(a), the velocity in the L regime (B2) adopts an approximately sinusoidal profile with a low amplitude ($\max |u| \approx 0.3$), whereas in the SW, TW, I and T regimes (B5, B6, B8 and B10, respectively), the mean velocity varies nearly linearly with height for $|z|\lesssim 0.5$, and $\max |u| \approx 1$. The vertical locations of the peaks in velocity, initially around $z \approx \pm 0.5$ in the L regime, shift slightly towards the top and bottom walls $z \approx \pm 0.7$ in the I and T regimes. These observations agree qualitatively with the experimental profiles of Lefauve et al. (2019) in the four regimes (see their figures 3f,l and 4f,l). However, exact agreement should not be expected as their $\theta$, ${Re}$ and ${Pr}$ values differ from ours.

In figure 9(b), the density profile resembles an error function in the L regime (B2) whereas it has a partially mixed layer in the W and I regimes (B5, B6, B8 and I, Exp) identified by a central region of reduced gradient (a layer) flanked by two regions of enhanced gradient (two interfaces). These three profiles are almost identical, with the nuance that the intermediate layer becomes slightly thicker from B5 to B6 to B8, as expected. Finally, in the T regime (B10 and T, Exp.), the middle layer becomes noticeably thicker, and the two interfaces flanking it become less sharp, leading to a profile approaching a uniform stratification.

In figure 9(c), the L flow has $Ri_g\approx 0.5$ throughout the central quarter of the height of the duct, flanked by steeply increasing values. The W and I flows have ${Ri}_g\approx 0.1$ near $z=0$. The T flow has a broader minimum with ${Ri}_g\approx 0.1\unicode{x2013}0.15$. The ${Ri}_g$ profiles in the I and T flows are qualitatively consistent between DNS and experiments, showing that, despite the difference in parameters, some key dynamical features of turbulence are not sensitive to the fluid properties.

Note that ${Ri}_g<0.25$ at $z=0$ in the W, I and T flows, but not in the L flow. Therefore, the mean profiles in the non-laminar flows reassuringly satisfy the Miles–Howard criterion necessary for the development of instabilities in a steady, inviscid, Boussinesq, parallel stably stratified shear flow. Moreover, in the T regime ${Ri}_g(z) \approx {Ri}_e \approx 0.1\unicode{x2013}0.15$ (i.e. robustly below the Miles–Howard criterion of 0.25). This agrees with the experimental conclusions of Lefauve & Linden (2022a) (see their figure 5), drawn from a wider dataset of 16 flows with increasing levels of turbulence. The reasons for this particular equilibrium, originally suggested by Turner (1973) and much observed since in numerical, experimental and observational data, are still debated. Other authors (Thorpe 2010; Smyth & Moum 2013; Salehipour et al. 2018) have called it ‘self-organised criticality’ or ‘marginal stability’, and found values ranging from ${Ri}_e\approx 0.07$ to $0.25$ in various stratified shear flows that differ (sometimes significantly) from SID (Lefauve & Linden 2022a).

5.2. Kinetic energy

We now study the spatio-temporal dynamics of kinetic energy along the entire length of the duct in our DNS. Similar experimental diagnostics are not yet available since the resolution of video cameras and geometry of the laser sheet limit us to shorter windows spanning only a limited part of the duct length.

We start by decomposing the velocity into mean and turbulent (fluctuating) components. Lefauve & Linden (2022b) defined the mean as the $x,t$ average. However, as our DNS data are available along the entire length of the duct, the flow (especially $u$) becomes noticeably inhomogeneous in $x$ (see figure 5a–c). A simple $x$ average would therefore make the ‘turbulent’ component artificially large by incorporating a significant non-parallel – but laminar – component. To resolve this, we define the mean and fluctuations using a moving average:

(5.2)$$\begin{gather} \bar{\boldsymbol{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{gather}$$

(5.3)$$\begin{gather}\boldsymbol{u}^\prime_m(x,y,z,t)\equiv\boldsymbol{u}-\bar{\boldsymbol{u}}_{\boldsymbol{m}}, \end{gather}$$

and the respective ‘moving’ mean kinetic energy (MKE) and TKE are

(5.4)$$\begin{gather} \bar{k}_{m}\equiv\tfrac{1}{2}\bar{\boldsymbol{u}}_m \boldsymbol{\cdot} \bar{\boldsymbol{u}}_m, \end{gather}$$

(5.5)$$\begin{gather}k^\prime_{m}\equiv\tfrac{1}{2}\boldsymbol{u}^\prime_m \boldsymbol{\cdot} \boldsymbol{u}^\prime_m. \end{gather}$$

The length of the averaging stencil $\Delta L=10$ was chosen to maximise the time- and duct-volume-averaged MKE $\langle \bar {k}_{m}\rangle _{\mathcal {V},t}$, as shown in figure 10(d).

Figure 10. Snapshots of DNS B6 at $t=160$ showing the (a) density, (b) moving-averaged TKE $k^\prime _m$ defined in (5.5) and (c) ‘naive’ TKE based on $\boldsymbol {u}'=\boldsymbol {u}-\langle \boldsymbol {u}\rangle _{x,t}$. (d) Bulk MKE $\langle \bar {k}_m \rangle _{\mathcal {V},t}$ as a function of the stencil length $\Delta L$. The dashed line corresponds to our choice in the remainder of the paper $\Delta L=10$.

In figure 10(a–c) we demonstrate the use of this moving average with a snapshot in the travelling wave regime (B6, as in figure 5c). The underlying turbulence ‘hotspots’ visualised by the density field (figure 10a) are faithfully captured by our moving-averaged TKE $k^\prime _{m}$ (figure 10b), whereas they are greatly overestimated by the ‘naive’ TKE based on the $x,t$-averaged velocity (figure 10c), which is equivalent to setting $\Delta L=2A=60$ and $x=0$ in (5.2).

Figure 11 shows time series of the duct-volume-averaged MKE $\langle \bar {k}_{m} \rangle _\mathcal {V}$ (figure 11a) and TKE $\langle k^\prime _{m} \rangle _\mathcal {V}$ (figure 11b) for the five regimes shown in figure 5. The laminar (B2) and the stationary wave (B5) cases quickly reach a steady state with constant or near-constant MKE (figure 11a) and zero or near-zero TKE (figure 11b). We note that the flow in B5 is faster than in B2 as their MKE plateau at $\approx$0.2 and $\approx$0.06, respectively. The travelling wave (B6) case shows more fluctuations in the MKE (but also $\approx$0.2), and a larger TKE fluctuating between $\approx$0.001 and 0.005. The intermittent and turbulent cases (B8 and B10) show much larger fluctuations in MKE and TKE, and a significantly larger TKE than in the travelling wave case. In these two cases, the MKE and TKE fluctuations are of comparable magnitude to their temporal mean. The TKE fluctuations are particularly striking, showing that the flow, even when averaged over the entire duct volume, is alternating between phases of intense and weak turbulence. These fluctuations appear to be quasi-periodic with a period of $\approx$100 ATU, corresponding to approximately one-and-a-half full-duct transit times at advective speed 1. Although long known from experiments, the mechanisms responsible for these fluctuations remain poorly understood and beyond the scope of this paper. The intermittent regime (B8) differs from the turbulent regime (B10) in that its TKE occasionally drops to zero for extended periods of time (here $150\lesssim t \lesssim 210$); i.e. the flow relaminarises in the duct. This never happens in the turbulent regime, although it does feature cycles of weaker and stronger turbulence.

Figure 11. Time series of volume-averaged (a) MKE and (b) TKE for the five cases B2–B10 of figure 5. A snapshot of density and TKE in B6 at $t=160$ is shown in figure 10(a,b).

Our MKE and TKE time series are approximately similar to their salt-stratified experimental counterparts in Lefauve & Linden (2022b) (see their figures 1(a,b) and 3( f,i,o,r), noting that they correspond to $Pr\approx 700$). This agreement between the DNS and experiments extends from the values of MKE $\approx$0.2 in the W/I/T regimes to the mean TKE $\approx 0.01\unicode{x2013}0.02$ in the I/T regimes, corresponding to a typical turbulent/mean velocity ratio of $\sqrt {k^\prime _m/\bar {k}_m} \approx 20\,\%\unicode{x2013}30$ % (although the ratio appears to be slightly higher in the DNS B10 than in the experiments T2 and T3). The large temporal fluctuations in B10 (which is at the limit of our computational resources) are typical of a flow near the I/T regime transition rather than well into the T regime, as was already clear from its location on the regime diagram (figure 6b). The experimental time series of TKE (see Lefauve & Linden (2022b), figure 3(o,r) in their datasets T1 and T3) suggest that these temporal fluctuations would significantly decrease in a more highly turbulent flow at higher ${Re}$ and $\theta$ (i.e. that the TKE would become increasingly steadily sustained at a higher level). Our B8 has a time series similar to their dataset T1 (both being near the I/T transition), whereas their more turbulent dataset T3 is further away from the I/T transition.

Our MKE and TKE time series are generally out of phase in the I and T regimes (figure 11), although the negative correlation is less clear in the T regime. In other words, in B8, the MKE tends to decrease as the TKE increases (i.e. turbulence slows down the mean flow), and vice versa (the MKE tends to increase when the TKE decreases or is zero). This behaviour, wherein the mean flow and the turbulence appear to regulate one another, supports the ideas of ‘self-organised criticality’ and ‘marginal stability’ discussed previously. In B10, this negative correlation holds until $t\approx 200$, at which point the TKE increases rapidly while the MKE continues increasing, leading both TKE and MKE to peak approximately at the same time. In the T regime the mean flow thus appears closer to a turbulent threshold, such that perturbations grow more readily without ‘waiting’ for the mean flow to fully accelerate. Equivalently, in the T regime (which has the highest $\theta \,{Re}$), the mean flow is able to keep accelerating despite the growing turbulence, presumably due to a higher forcing (because it is proportional to $\theta$) and a lower TKE dissipation than in the I regime (because it is inversely proportional to ${Re}$) .

Finally, figure 12 shows $x$–$t$ diagrams of TKE (averaged along $y$ and $z$) for B2 to B10 (left to right) after the initial transients have decayed ($t>80$). In the L regime (B2, figure 12a), the TKE is negligible except very near the ends of the duct $|x|\approx 30$, where tiny fluctuations are found where the exchange flow discharges into the reservoirs. In the stationary W regime (B5, figure 12b), the TKE at both ends is higher, extends a little further into the duct and is also occasionally visible near the centre of the duct $|x|\approx 0$ where it appears to remain stationary. Similarly, the ‘end waves’ do not propagate into the duct and may be swept into the reservoirs, implying that their phase speed is smaller than the convective speed of the flow (i.e. that the flow may be critical or supercritical in these regions).

Figure 12. Spatio-temporal diagram of TKE $\langle k^\prime _m \rangle _{y,z}(x,t)$ for $t\in [80, 260]$ (after the initial transients) in DNS (a) B2, (b) B5, (c) B6, (d) B8 and (e) B10. Note the colourbar is in log scale.

In the travelling W regime (B6; figure 12c), larger TKE develops and it sometimes appears to propagate along the duct. These waves appear to be generated within the duct rather than travelling from the ends. In the I regime (B8, figure 12d), a laminar phase develops between $120 \lesssim t \lesssim 220$, lasting over a full duct transit time (taking ${\approx }2A=60$ ATU at the maximum flow speed $\approx$1). The boundary between laminar and turbulent phases appears to propagate from one end of the duct to the other. In the T regime, a quiescent patch develops near the centre of the duct just after $t=200$. The quiescent phase ends when energetic turbulent regions move in from both ends of the duct.

5.3. Pressure

We now analyse the pressure field throughout the interior of the duct, which is inaccessible to experiments.

Figure 13 shows representative snapshots of the spanwise-averaged pressure $\langle p \rangle _y$ and five equally spaced isopycnals for B2, S6, B6, S8 and B8. Note that our definition of the non-dimensional density $\rho$ as the perturbation around the reference $\rho _0$ (see § 2.1) implicitly subtracts the hydrostatic pressure due to $\rho _0$ in the reservoirs. The pressure distribution is qualitatively similar in the W and I regimes (S6, B6, S8, B8; figure 13b–e), but different in the L regime (B2; figure 13a).

Figure 13. Spanwise-averaged pressure field snapshots (colours) superimposed with five isopycnals (lines) $\rho =0,\pm 0.4, \pm 0.8$ in DNS (a) B2, $t=100$, (b) S6, $t=108$, (c) B6, $t=108$, (d) S8, $t=235$ and (e) B8, $t=250$.

In the L flow inclined at $\theta =2^\circ$ (figure 13a), the pressure conforms to what we expect from an exchange flow in a horizontal duct where $\theta$ does not play a major role (see the sketch in Lefauve (2018), figure 1.4). Essentially, each layer experiences a favourable streamwise pressure gradient: $-\partial _x p > 0$ in the right-flowing lower layer where $u>0$, causing a convective acceleration along the duct, $u \partial _x u >0$; and vice versa, $-\partial _x p < 0$ in the left-flowing upper layer where $u<0$, causing the expected $u \partial _x u <0$ ($\partial _x u >0$ in both layers). This is achieved by a high-pressure zone in the bottom left and top right reservoirs (in red), and a low-pressure zone in the top left and bottom right reservoirs (in blue), as would be naturally obtained by the hydrostatic equilibrium of two solutions having different densities and required to match hydrostatic pressures at mid-height.

However, in flows inclined at $\theta =6^\circ$ and $8^\circ$ (figure 13b–e), the pressure has a large-scale global minimum in the centre of the duct (in blue). This low-pressure zone causes both layers to experience a favourable pressure gradient over approximately the first half of their course, causing the fluid to accelerate as it flows towards the centre of the duct, but an adverse pressure gradient over the second half of their course, causing the fluid to decelerate from the centre of the duct as it flows away from the centre of the duct. This adverse pressure gradient confirms the predictions of Lefauve & Linden (2022a) (see their § 4.3) who, without having direct access to the pressure field, found that in most datasets the Reynolds-averaged budget implied the existence of an adverse pressure gradient.

Furthermore, these features of the pressure distribution are found in both the BR and SR geometries (compare figures 13(b,c) and 13(d,e), despite the difference in the location of the low-pressure region in the I cases) as the outflowing layers must decelerate in both geometries. In the BR (and by extension in AR and benchmark), this occurs since the streams encounter fluid at rest in a large reservoir; in the SR, this occurs since the streams encounter the artificial forcing region where the flow is brought to rest. This suggests that the SR cases appear to qualitatively mimic the benchmark case despite the absence of reservoirs, which may help reduce computational costs further in future studies.

Finally, we address the impact of the adverse pressure gradient on the density field and its interface(s). The isopycnals in figure 13(a–d) (the yellow lines denoting the lower, densest layer, and the dark blue lines denoting the upper, lightest layer) show that the central low-pressure zone coincides with an increasing depth of each layer in the direction of flow. This is expected from mass conservation along a straight duct: an accelerating layer (because of a favourable pressure gradient or gravitational forcing) must become thinner along its course, and vice versa, a decelerating layer (caused by an adverse pressure gradient) must become thicker. We also see from the isopycnals that this thickening coincides with the emergence of displaced isopycnals (in B6, S6) and turbulence (in B8, S8). Note that this interface thickening implies the occurrence of internal hydraulic ‘jump’ and the set-up of hydraulic control (Meyer & Linden 2014; Lefauve & Linden 2020a) in the middle of the duct. This effect of internal hydraulics is beyond the scope of this paper and will be revisited in more detail in our future work.

5.4. Turbulent energy fluxes

We conclude this section with an analysis of the turbulent energy fluxes in datasets B5, B6, B8 and B10 and use it to demonstrate how DNS data allow us to overcome the current experimental limitations identified in Lefauve & Linden (2022b) (their Appendix B) and improve our physical understanding of stratified turbulent mixing.

To do so, we adopt the ‘shear-layer’ non-dimensional framework of Lefauve & Linden (2022a) (§ 3.3) by rescaling all velocities such that the $(x,t)$-averaged $\langle u \rangle _t(y=0,z)$ has extrema $\pm 1$, and rescaling all spatial variables such that the $z$ locations of these extrema are located at $\pm 1$ (whereas previously the top and bottom walls were located at $\pm 1$); we call this central region of non-dimensional height $=2$ the shear layer. This rescales the effective Reynolds number and bulk Richardson number of the flow, which we now denote $Re^s$ and $Ri_b^s$, respectively, and allows for more meaningful comparison of datasets with one another as well as with the literature. Following Lefauve & Linden (2022a) we further remove from our analysis all data outside the shear layer, i.e. exclude the top and bottom near-wall boundary layers in $z$, as well the boundary layers in $y$ (where the peak $|u|$ is less than 0.7), in order to focus on the ‘core’ region with turbulent activity.

We then define the non-dimensional time- and volume-averaged MKE $\bar {K} = \langle \bar{k}_m \rangle _\mathcal {S}$ and TKE $K' = \langle k'_m \rangle _\mathcal {S}$, as well as the mean scalar variance $\bar {K}_\rho \equiv Ri^b_s \langle \bar {\rho }_m^2 /2 \rangle _\mathcal {S}$ and turbulent scalar variance $K_\rho ' \equiv Ri^b_s \langle \rho '_m{}^2 /2 \rangle _\mathcal {S}$. Note that the subscript $m$ denotes the moving average introduced in § 5.2, and that the multiplying factor $Ri_b^s$ allows us to interpret $K_\rho,K_\rho '$ as proxies for potential energy under linear stratification. The simple bracket averaging $\langle \cdot \rangle _\mathcal {S}$ denotes a combined three-dimensional volume averaging over the shear-layer region and along the central two-thirds of the duct (by excluding one averaging window $\Delta L$ on either end) and time averaging over $t\in [80, 280]$ (focusing on the established dynamics as in figure 12).

Considering the evolution equations of these four energy reservoirs under a set of ‘safe’ approximations in SID, Lefauve & Linden (2022b) derived the following approximate balances between energy fluxes in a statistically steady state:

(5.6a)$$\begin{align}\mathcal{P} &\approx \mathcal{F} - \bar{\epsilon} \,\quad \text{(production of } K' = \text{forcing} - \text{laminar dissipation)}, \end{align}$$

(5.6b)$$\begin{align}\mathcal{E} &\approx \mathcal{P} - \mathcal{B} \quad \text{(dissipation of } K' = \text{production} - \text{buoyancy flux)}, \end{align}$$

(5.6c)$$\begin{align}\mathcal{P}_\rho &\approx \varPhi^{\bar{K}_\rho} \quad\quad \text{(production of } K'_\rho = \text{boundary net flux of unmixed fluid)}, \end{align}$$

(5.6d)$$\begin{align}\chi &\approx \mathcal{P}_\rho \ \ \quad\quad \text{(dissipation of } K'_\rho \text{ i.e. mixing} = \text{production)}, \end{align}$$

where the eight fluxes are

(5.7)\begin{gather} \left.\begin{gathered} \mathcal{P} \equiv{-} \langle u'v'\partial_y \bar{u}+ u'w' \partial_z \bar{u} \rangle_\mathcal{S}, \quad \mathcal{F} \equiv Ri_b^s \sin \theta \langle \bar{u}\bar{\rho}\rangle_\mathcal{S} , \\ \bar{\epsilon} \equiv \frac{2}{Re^s} \langle||\bar{\textsf{\textit{\textbf{s}}}}||^2\rangle_\mathcal{S}, \quad \mathcal{E} \equiv \frac{2}{Re^s} \langle ||{\boldsymbol{\mathsf{s}}}'||^2 \rangle_{\mathcal{S}}, \end{gathered}\right\} \end{gather}

(5.8)\begin{gather} \left.\begin{gathered} \varPhi^{\bar{K}_\rho} \equiv{-}Ri_b^s \left\langle \frac{\left.\left( u\dfrac{\rho^2}{2}\right)\right|^{x^+}_{x^-}}{x^+{-} x^-} \right\rangle_\mathcal{S}, \quad\mathcal{B} \equiv Ri_b^s \langle w' \rho' \rangle_\mathcal{S}, \\ \mathcal{P}_\rho \equiv{-} Ri_b^s \langle w' \rho' \partial_z \bar{\rho} \rangle_\mathcal{S}, \quad \chi \equiv \frac{Ri_b^s}{Re^s \, Pr}\langle |\boldsymbol{\nabla} \rho'|^2\rangle_\mathcal{S}, \end{gathered}\right\} \end{gather}

where $\bar {{\boldsymbol{\mathsf{s}}}}, {\boldsymbol{\mathsf{s}}}'$ are the strain rate tensors of the mean and turbulent velocity fields, respectively. From (5.7) and henceforth we omit the subscripts $m$, using the moving average in all bar and prime quantities. Equations (5.6b) and (5.6d) are the classical balances of Osborn (1980) and Osborn & Cox (1972), respectively, while (5.6a) and (5.6c) are specific to SID: the mean kinetic energy is sustained through $\mathcal {F}$ by the gravitational acceleration and the mean scalar variance is sustained through $\varPhi ^{\bar {K}_\rho }$ by the inflow of unmixed fluids from the reservoirs into the volume of interest.

Figure 14(a–d) demonstrates how closely the four balances (5.6a–d) hold in DNS, where the thick diagonal solid line denotes equality between the left-hand side (vertical axis) and right-hand side (horizontal axis) of each equation. The open symbols in figure 14(a,b,d) denote the value of the left- and right-hand side exactly as in (5.6a,b,d), and are generally in balance, except in a few cases. In these few cases, the balance is improved by the filled symbols obtained after adding the following boundary fluxes:

(5.9a)\begin{align} \varPhi^{\bar{K}} &\equiv{-} \left\langle \frac{\left.\left(u \left(\dfrac{{\boldsymbol{u}}^2}{2}+p\right)\right)\right|^{x^+}_{x^-}}{x^+{-} x^-} \right\rangle_{\mathcal{S}}\ \text{to the right-hand side of } {(5.6a)},\end{align}

(5.9b)\begin{align} \varPhi^{K'} &\equiv{-} Ri_b^s \left\langle \frac{\left.\Biggl(u' \Biggl(\dfrac{{\boldsymbol{u}'}^2}{2}+p'\Biggr)\Biggr)\right|^{x^+}_{x^-}}{x^+{-} x^-} + \dfrac{\left.\Biggl(v' \Biggl(\dfrac{{\boldsymbol{u}'}^2}{2}+p'\Biggr)\Biggr)\right|^{y^+}_{y^-}}{y^+{-} y^-} \right.\nonumber\\ &\quad \left. + \dfrac{\left.\Biggl(w' \Biggl(\dfrac{{\boldsymbol{u}'}^2}{2}+p'\Biggr)\Biggr)\right|^{z^+}_{z^-}}{z^+{-} z^-} \right\rangle_\mathcal{S} \ \text{to the right-hand side of } {(5.6b)}, \end{align}

(5.9c)\begin{align} \varPhi^{K'_\rho} &\equiv{-} Ri_b^s \left\langle \frac{\left.\left(u'\dfrac{\rho' {}^2}{2}\right)\right|^{x^+}_{x^-}}{x^+{-} x^-} + \frac{\left.\left(v' \dfrac{\rho'{}^2}{2}\right)\right|^{y^+}_{y^-}}{y^+{-} y^-} \right.\nonumber\\ &\quad \left. + \frac{\left.\left(w' \dfrac{\rho'{}^2}{2}\right)\right|^{z^+}_{z^-}}{z^+{-} z^-} \right\rangle_\mathcal{S} \text{to the right-hand side of } {(5.6d)},\end{align}

where $x^\pm, y^\pm, z^\pm$ are the edges of our shear-layer averaging domain, and we neglected in (5.9a) the spanwise and vertical mean transport, in (5.9a)–(5.9b) the work of viscous forces and in (5.9c) the transport by molecular diffusion, in order to focus on the dominant contributions. The fact that these fluxes improve the balance only slightly demonstrates that they can, to a reasonable approximation, be neglected in SID energetics. This result was hypothesised in Lefauve et al. (2019) and Lefauve & Linden (2022b) but their experiments could not clearly confirm it due to excessive noise in the computation of (5.9) and the lack of the pressure. After correction, the only remaining discrepancy is found in $\chi$ in the most turbulent flow B10 (figure 14d), where it is $33\,\%$ below the expected value to balance $\mathcal {P}_\rho$.

Figure 14. Correlation between some of the time- and volume-averaged energy fluxes (5.7)–(5.8) in B5, B6, B8 and B10. (a–d) Verification of the energy balances in (5.6) (left-hand side versus right-hand side), filled symbols denoting the minor correction of boundary fluxes (5.9). (e, f) Verification of the empirical relations (5.10). (g) Determination of the empirical flux coefficient (5.11) suggesting $\varGamma =0.1$.

Part of this discrepancy ($13\,\%$ out of $33\,\%$) is explained by our neglect of molecular diffusion $Ri_b^s/(Re^s Pr)|\boldsymbol {\nabla }\bar {\rho }|^2$ in (5.6d), which at these values of $Re^s$ and $Pr$ is ‘only’ a factor of five smaller than $\chi$. The remaining $20\,\%$ of this discrepancy is most likely be due to under-resolution, as we explain next. Our spatial grid in B10 approaches the time- and volume-averaged Batchelor length scale computed in the shear layer $[\Delta x, \Delta y, \Delta z] = [3.3, 2.3, 2.3] \, \ell _B$, where $\ell _B \equiv \mathcal {E} ^{-1/4}(Re^s)^{-3/4}Pr^{-1/2}= 0.01$ in non-dimensional shear-layer units. We verified that $\chi$ was nearly identical on the coarser grid (by a factor of $\approx$2 in $x$ and $y$; see table 1). However, the time series of the volume-averaged $\mathcal {E}$ undergoes cycles (much like the TKE in figure 11b) with peaks being three times larger than the time average used for $\ell _B$. Using this peak value, our grid only resolves $\approx [4, 3, 3] \, \ell _B$, which may be too coarse to capture the entire spectrum of $\boldsymbol {\nabla }\rho '$ during the most turbulent events. This should be kept in mind for future DNS which aim to fully resolve mixing.

Having largely verified (5.6), we now examine in figure 14(e, f) two robust empirical relations from the SID experimental literature in very turbulent flows ($\mathcal {E} \gg \bar {\epsilon }$):

(5.10a)$$\begin{gather} \mathcal{P_\rho} \approx \mathcal{B} \quad \text{because } \partial_z \langle \rho\rangle_{x,t} \approx{-}1 \text{ in the the shear layer}, \end{gather}$$

(5.10b)$$\begin{gather}\mathcal{E} + \bar{\epsilon} \approx \mathcal{E} \approx 0.035 \theta \quad \text{because of hydraulic control and } Ri^b_s \approx 0.10\unicode{x2013}0.15, \end{gather}$$

where $\theta$ is in radians (Lefauve & Linden 2022b). Our DNS data generally confirm (5.10a) (figure 14e) and (5.10b) (figure 14f) but with two reservations. First, in B10, $\mathcal {P}_\rho$ is $25\,\%$ below the expected value $\mathcal {B}$ as a result of the mean vertical density gradient being slightly weaker in our DNS (at $Pr=7$) than in the experiments (at $Pr=700$). Second, although the total dissipation $\mathcal {E}+\bar {\epsilon }$ corrected with the appropriate boundary fluxes (filled symbols) follows (5.10b) (in particular since our $Ri_b^s$ indeed converges to 0.12–0.14 in all datasets B5–B10), the agreement is less good for the turbulent dissipation alone $\mathcal {E}$ (smaller open symbols). In other words, our ability to fully capture $\mathcal {E}$ and all boundary fluxes in DNS allows us to quantify the relative importance of the (subdominant) terms $\varPhi ^{\bar{K}},\varPhi ^{K'}$ and $\bar {\epsilon }$ in SID energetics. We hypothesise that stronger turbulence at higher values of $\theta \,Re^s$ would see these subdominant terms plateau, and that $\mathcal {E}$ would follow $0.035 \theta$ increasingly closely.

We now move to the ultimate goal of this energetics analysis, and more broadly of research into turbulent mixing, which is to eventually connect all turbulent fluxes to known non-dimensional parameters of the flow in a closed system of equations. In the asymptotic ‘strong SID turbulence’ scenario under which $\bar {\epsilon }$ becomes subdominant, we have seven key turbulent fluxes in (5.7), (5.8) and six independent equations: the four equations (5.6) expressing conservation of energy, and the two robust empirical equations (5.10), one of which crucially involves the single input parameter $\theta$. To close the system, we require a seventh independent equation, which we choose to be the classical flux parameter in the ocean mixing literature:

(5.11)\begin{equation} \varGamma \equiv \frac{ \mathcal{B} }{\mathcal{E}}. \end{equation}

Combining these seven equations in matrix form, and inverting this linear system, we deduce all fluxes in closed form as

(5.12) \begin{align} &\begin{bmatrix} 1 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & -1 \\ 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & \varGamma & -1 & 0 & 0 & 0 \\ \end{bmatrix} \begin{bmatrix} \mathcal{F} \\ \mathcal{P} \\ \mathcal{E} \\ \mathcal{B} \\ \mathcal{P}_\rho \\ \varPhi^{\bar{K}_\rho} \\ \chi \end{bmatrix} = \begin{bmatrix} 0 \\ 0\\ 0 \\ 0\\ 0\\ 0.0035\theta \\ 0 \end{bmatrix} \nonumber\\ &\quad \Longrightarrow\ \begin{bmatrix} \mathcal{F} \\ \mathcal{P} \\ \mathcal{E} \\ \mathcal{B} \\ \mathcal{P}_\rho \\ \varPhi^{\bar{K}_\rho} \\ \chi \end{bmatrix} = 0.035\theta \begin{bmatrix} 1+\varGamma \\ 1+\varGamma \\ 1 \\ \varGamma \\ \varGamma \\ \varGamma \\ \varGamma \end{bmatrix} . \end{align}

This expression highlights the importance of knowing the value of $\varGamma$, and its potential dependence on any of the non-dimensional flow parameters such as $\theta, Re$ or $Pr$, as the keystone to turbulent mixing in SID. Figure 14(g) shows $\mathcal {B}$ versus $\mathcal {E}$ (both of which are fully resolved in our DNS) which strongly supports the constant value of $\varGamma \approx 0.1$ in the two most turbulent datasets (0.097 and 0.096 in B8 and B10, respectively). The experimental data Lefauve & Linden (2022b) also found $\varGamma \approx 0.1$ but their insufficient spatial resolution to fully capture $\mathcal {E}$ led them to conjecture a slightly lower value in the range 0.05–0.07. Our DNS allow us to confirm that $\varGamma \approx 0.1$ is a robust estimate, at least for turbulence at $Pr=7$ in this narrow region of the $(\theta,Re)$ space. This value is smaller than the classical $\varGamma \approx 0.2$ found in most of the literature, potentially because of the shear-driven nature of SID turbulence, sustained by gravitational forcing, inducing a different type of equilibrium from that in other canonical flows.

6. Conclusions

In this paper, we performed and interpreted DNS of stratified shear flows in a SID, that is, a long rectangular duct connecting two reservoirs. The flow is continuously forced by gravity by a modest positive tilt angle $\theta =0^\circ{-}10^\circ$, and has Reynolds number ${Re}=400{-}1250$, bulk Richardson number ${Ri}=0.25$ and Prandtl number ${Pr}=7$. Our results are summarised as follows.

6.1. An efficient numerical paradigm for SID

In § 2 we presented a new numerical set-up (figure 1) designed to closely mimic the experimental set-up of the SID. We introduced a new forcing term in the reservoirs (figure 1) that allows the exchange flow to be sustained indefinitely with smaller reservoirs than the experiments (Lefauve et al. 2019), thus focusing our computational resources on the flow of interest within the duct. We also implemented an IBM (figure 2) to enforce the boundary conditions on the duct walls that match the experiments, i.e. no-slip for velocity and no-flux for density.

In § 3 we validated this numerical configuration. First, we showed that our artificial forcing in the reservoirs was necessary to sustain the exchange flow by ‘refreshing’ any finite-sized reservoirs beyond the short time scale over which they would otherwise fill with mixed fluid (figure 3). Second, we showed that small reservoirs combined with the appropriate forcing were sufficient to reproduce the flow of a benchmark case having very large reservoirs and no forcing (figure 4). Additional evidence that the key flow kinetics (e.g. wave propagation and turbulent intermittency) and statistics (e.g. mean flow quantities) in the SID are largely independent of reservoir size, spanwise length and perturbation magnitude is provided in supplementary material S1. This independence has been explained by the establishment of hydraulic control (Meyer & Linden 2014; Lefauve 2018).

In § 4 we described the properties of increasingly disorganised and turbulent flow regimes found by increasing ${Re}$ and $\theta$ (figure 5). These regimes are similar to those found in experiments where the stratification is achieved by temperature (approximately matching our ${Pr}=7$), which further validates the relevance and accuracy of our DNS to faithfully reproduce experimentally realisable flows. These flow regimes are generally found in the same region of $\theta$–${Re}$ parameter space as the experiments (figure 6), with very little difference between the larger reservoir (AR) and the smaller reservoir (BR). This agreement between DNS and experiments carries over to more detailed flow characteristics visualised by instantaneous shadowgraph snapshots (figure 7) and spatio-temporal diagrams (figure 8).

In § 5 we studied quantitative DNS diagnostics that complement experimental diagnostics. We first investigated the vertical velocity and density profiles. The gradient Richardson number (figure 9) displayed the same turbulent ‘equilibrium’ with nearly uniform ${Ri}_g \approx 0.10{-}0.15$ across the shear layer as in the experiments, despite the difference in ${Pr}$. We then moved to the mean and turbulent kinetic energies (figure 10), introducing a moving-average (in $x$) definition suited to our DNS data, focusing on the intermittency of the turbulence (figure 11). We investigated the spatio-temporal behaviour of the TKE (figure 12), exploiting the fact that our DNS data are available along the full length of the duct. We contrasted stationary and travelling waves, described where waves originate from and how turbulence or relaminarisation sometimes occurs synchronously along the duct, and sometimes in ‘waves’ propagating at the advective speed. Next, we investigated the pressure field (figure 13) and discovered a large-scale low-pressure zone inside the duct in all non-laminar flows, which was previously conjectured with experimental data, but only proven now with DNS data. We showed that our ad hoc forcing, even in the smallest geometry SR, was sufficient to reproduce this behaviour observed with real reservoirs. This low-pressure zone creates a favourable pressure gradient in both layers over roughly the first half of their transit, allowing them to accelerate as they flow in, but, crucially, an adverse pressure gradient over roughly the second half of their transit, causing them to decelerate before they flow out. This result suggests a potentially new mechanism for hydraulic control in exchange flows tilted at a favourable angle $\theta$, which requires further study. Finally, we have largely confirmed the simplified model in (5.6) for the steady-state kinetic and scalar energy fluxes in SID turbulence, as well as the empirical relations in (5.10) (figure 14) hypothesised from experimental data. This allowed us to express in (5.12) all seven fluxes fully characterising the time- and volume-averaged turbulent energetics and mixing in SID as functions of $\theta$ and the flux coefficient $\varGamma$. Our data suggest $\varGamma \approx 0.1$ in the most turbulent flows, a value lower than the classical value of $0.2$ used in most of the ocean mixing literature. The physics behind this difference warrants further study.

6.2. Outlook

This paper introduced a computationally efficient way to simulate realistic shear-driven stratified turbulence over long time periods, which shows excellent agreement with the experiments, with all non-dimensional parameters being matched (at $Pr=7$). We consider this comprehensive agreement between highly nonlinear numerical and experimental fluid dynamics to be the major result of this paper, and a milestone in the study of SID and stratified turbulence. Furthermore, the numerics add considerable value to the experiments by providing accurate and highly resolved data over the entire domain, and by allowing arbitrarily long integration times with the addition of forcing terms in the reservoirs.

There is significant scope to build on this study, by overcoming technical challenges. For example, improved experimental technology is needed to obtain more accurate, higher-resolution data in more highly turbulent flows at ${Re}=O(10^3{-}10^4)$. Increased computational power is also needed to match such ${Re}$ and ${Pr}$ to tackle the differences between temperature $(Pr=O(1))$ and salt stratification $({Pr}=O(10^2))$. Finally, studying the slow, quasi-periodic dynamics of intermittent turbulence requires large physical reservoirs and data acquisition as well as long integration times and costly simulations. Nevertheless, as experimental technology improves and computational power increases, we anticipate that these techniques will be able to cover a much larger range in parameter space and answer questions previously inaccessible to theory, observations, experiments or simulations alone.