Parameterization of mixing in stratified open channel flow

Abstract

The dynamics and parameterization of mixing in temporally evolving turbulent open-channel flow is investigated through direct numerical simulations as the flow transitions from an initially neutral state to stable stratification. We observe three distinctly different mixing regimes separated by transitional values of turbulent Froude number $Fr$: a weakly stratified regime for $Fr >1$; an intermediate regime for $0.3< Fr<1$; and a saturated regime for $Fr<0.3$. The mixing coefficient $\varGamma =B/\epsilon _K$, (where $B$ is the buoyancy flux and $\epsilon _K$ is the dissipation rate of kinetic energy), is well predicted by the parameterization schemes of Maffioli et al. (J. Fluid Mech., vol. 794, 2016) and Garanaik & Venayagamoorthy (J. Fluid Mech., vol. 867, 2019, pp. 323–333) across all three regimes through instantaneous measurements of $Fr$ and the ratio $L_E/L_O$, where $L_E$ and $L_O$ are the Ellison and Ozmidov length scales, respectively. The flux Richardson number $R_f = B/(B+\epsilon _K)$ shows linear dependence on the gradient Richardson number $Ri_g$ up to a transitional value of $Ri_g =0.25$, past which it saturates again to a constant value independent of $Fr$ or $Ri_g$. By examining the flow as a balance of inertial, shear and buoyancy forces, we derive physically based scaling relationships to demonstrate that $Ri_g \sim Fr^{-2}$ and $Ri_g \sim Fr^{-1}$ in the weakly and moderately stratified regimes and that $Ri_g$ becomes independent of $Fr$ in the saturated regime. Our results suggest that the $L_E/L_O \sim Fr^{-1}$ scaling of Garanaik & Venayagamoorthy (J. Fluid Mech., vol. 867, 2019, pp. 323–333) in the intermediate regime manifests due to the influence of mean shear. The differences in the relationships between $Fr$ and $L_E/L_O$ for non-sheared flows within this regime are discussed.

1. Introduction

Due to the ubiquity of geophysical turbulent flows under the effect of stable density stratification, investigation into the mechanisms and accurate quantification of diapycnal mixing has been and remains a fundamental area of research within the stratified turbulence community. Reviews of past work include Ivey, Winters & Koseff (2008), Peltier & Caulfield (2003) and Caulfield (2020). In particular, mixing within wall-bounded stratified flows creates a unique and interesting set of mixing dynamics as the inherent spatial inhomogeneity allows for the simultaneous coexistence of various energetic mixing regimes within a single flow (Taylor, Sarkar & Armenio 2005; García-Villalba & del Álamo 2011; Williamson et al. 2015). In this study we investigate the mixing dynamics within temporally evolving open channel flow through its transition from an initially neutral to a stably stratified state. The emphasis of our study falls on two key ideas: the accurate parameterization of mixing efficiency in the context of our spatiotemporally inhomogeneous flow; and a robust investigation into the dynamics and relationship between key non-dimensional parameters within the varying flow regimes.

Central to the quantification and estimation of mixing in stratified flows are the diapycnal diffusivity $K_P$ and the mixing efficiency coefficient $\varGamma$, which are linked through the relation

(1.1)\begin{equation} K_P = \varGamma \frac{\epsilon_K}{N^{2}}, \end{equation}

where $\varGamma = R_f/(1-R_f)$, $R_f$ is the mixing efficiency or flux Richardson number, $\epsilon _K$ is the dissipation rate of turbulent kinetic energy and $N$ is the buoyancy frequency. Historically, Osborn (1980) argued that under equilibrium conditions $\varGamma$ can be assumed to have a constant value of $0.2$, however, it has since been demonstrated that $\varGamma$ can significantly vary with respect to the energetic state of the flow (Shih et al. 2005; Ivey et al. 2008; Venayagamoorthy & Koseff 2016). As such, numerous parameterization schemes have been proposed in the literature to estimate $\varGamma$ based on relevant non-dimensional parameters. However, as summarized in Gregg et al. (2018), a significant challenge within the study of stratified turbulence is the plethora of varied parameterization schemes for $\varGamma$ and the subsequent ambiguity in the relationship between the different parameters. Further, as outlined by Caulfield (2021), a limitation of numerous parameterization frameworks is that they are derived under the assumptions of homogeneity and stationarity and further tested within idealized triply periodic domains where such homogeneity and stationary is enforced and the statistics are correlated after appropriate spatial and temporal averaging. Although such flows are extremely useful for evaluating flow properties in a precise parameter range, real flows, however, can exhibit significant variability in time and space, often resulting in disparity between instantaneous correlations of flow properties to their respective non-dimensional parameters relative to a statistically stationary case. In this context, our spatiotemporally inhomogeneous channel flow in which no local parameters are externally enforced or known a priori, presents a robust testing ground for such schemes as well as a distinct opportunity to investigate the relationships and similarities between the varied parameterization frameworks.

Historically, due to the relative ease of measuring mean gradients, parameterization of the mixing efficiency in wall-bounded and shear flows has focused on the gradient Richardson number $Ri_g = N^{2}/S^{2}$, where $S$ is the mean shear. Throughout numerous studies it has been repeatedly shown that for stationary shear flows and within the upper limit of $Ri_g \lesssim 0.25$, the mixing efficiency displays a monotonic, essentially linear dependence on $Ri_g$ (Armenio & Sarkar 2002; Taylor et al. 2005; García-Villalba & del Álamo 2011; Chung & Matheou 2012; Deusebio, Caulfield & Taylor 2015; Karimpour & Venayagamoorthy 2015; Zhou, Taylor & Caulfield 2017a). Although the deviation of a monotonic relationship between $\varGamma$ and $Ri_g$ at $Ri_g = 0.25$ is conceptually consistent with the idea of critical gradient Richardson number $Ri_{g,c} = 0.25$ proposed in the seminal work of Miles (1961) as the idealized threshold for the formation of local shear instabilities, however, it remains unclear if stability is in fact the mechanism for the departure from monotonic behaviour or the value of $Ri_g = 0.25$ is simply ‘fortuitous’ (Salehipour, Peltier & Caulfield 2018). The critical value itself as applied to real three-dimensional flows is an area of debate in itself (Galperin, Sukoriansky & Anderson 2007). As summarized in Caulfield (2021), there, however, exists no singular behaviour for $\varGamma$ in the very stable limit of high $Ri_g$, with multiple evolution paths available in so-called ‘right flank’ of the flux-gradient curve of Phillips (1972).

In a recent study Maffioli, Brethouwer & Lindborg (2016) (henceforth referred to as MBL16) presented scaling arguments to propose an alternative parameterization framework through the turbulent Froude number $Fr = \epsilon _K/ N E_K$, where $E_K$ is the turbulent kinetic energy. In their paper and under the assumption of sufficiently high Reynolds number, they classify stratified turbulence into the ‘strongly stratified’ – $Fr \ll O(1)$ – and ‘weakly stratified’ – $Fr \gg O(1)$ – regimes with $\varGamma \sim \text {const.}$ and $\varGamma \sim Fr^{-2}$ scaling relationships, respectively. Garanaik & Venayagamoorthy (2019) (henceforth referred to as GV19) expand upon this idea to propose a separate ‘moderately stratified’ – $Fr = O(1)$ – regime where both buoyancy and inertial forces are significant, and derive scaling arguments to propose a novel $\varGamma \sim Fr^{-1}$ relationship within this regime. Expanding on similar concepts from past works such as Ivey & Imberger (1991) and Smyth, Moum & Caldwell (2001), they further propose that $Fr$ and subsequently $\varGamma$ can be inferred from the ratio of $L_E/L_O$ across all three regimes; where $L_E$ is the overturning Ellison length scale and $L_O$ is the Ozmidov length scale, with the underlying concept being that $L_E$ has been proven to correlate directly to the easily measurable overturning Thorpe length scale $L_T$ (Smyth & Moum 2000; Mater, Schaad & Venayagamoorthy 2013), thus inferring the mixing efficiency through field measurements of $L_T/L_O$ becomes a conceivably easier task rather than directly through $Fr$. As $Fr$ is a parameter composed of fundamental turbulent flow properties that inherently exist in stratified turbulence irrespective of physical boundaries or mean shear, MBL16 and GV19 hence both argue that an $Fr$-based framework may present a degree of universality across a broad range of stratified flows. However, the testing of its applicability to wall-bounded and shear flows remains relatively limited, in particular for stronger stratification levels of $Fr<1$.

The concept of a single unifying parameterization scheme becomes somewhat complicated when considering that in the presence of mean shear, $Ri_g$ and $Fr$ may not be independent parameters, in fact it has been suggested that a multiparameter framework may be necessary to accurately describe the mixing dynamics when both shear and buoyancy forces are present within the flow (Mater & Venayagamoorthy 2014). The two frameworks are reconciled for weakly stratified flow ($Fr \gg O(1)$ or $Ri_g \ll O(1)$) as it can readily be shown that $Ri_g \sim Fr^{-2}$ (Zhou et al. 2017a). However, in the moderately and strongly stratified regimes the relationship between $Ri_g$ and $Fr$ remains unclear. The underlying basis for the scaling in the strongly stratified regime of MBL16 draws directly from established strongly stratified turbulence theory of a regime defined by $Fr \ll O(1)$ and $Re_B \gg O(1)$ (Billant & Chomaz 2001; Riley & de BruynKops 2003; Lindborg 2006; Brethouwer et al. 2007), where $Re_B = \epsilon _K/ N^{2} \nu$ is the buoyancy Reynolds number and $\nu$ is the kinematic viscosity. As summarized in the review of Caulfield (2021), it is, however, still an open question whether sheared stratified turbulence can access this regime in the sense described. For instance, Thorpe & Liu (2009) hypothesize that sheared stratified turbulence inherently self-regulates within a loop between states of marginal stability and instability. Recent studies have shown support for such self-regulatory behaviour that appears to drive sheared Holmboe instability or ‘scouring’ driven turbulence towards a state described by the classic Miles and Osborne estimates of a critical gradient Richardson number $Ri_{g,c} \approx 0.25$ and mixing efficiency $\varGamma _{c} \approx 0.2$ (Salehipour et al. 2018; Smyth, Nash & Moum 2019). The matter becomes even further complicated if we consider that in weakly stratified flows, $Ri_g$ and $Re_B$ can also become deeply correlated such that $Ri_g \sim Re_B^{-1}$ (Riley & de BruynKops 2003; Hebert & de Bruyn Kops 2006; Chung & Matheou 2012), where $Re_B$ itself has been a parameter widely used to parameterize mixing (Shih et al. 2005); thus creating a complex multiparameter space in which $\varGamma,Ri_g,Fr,Re_B$ may all conceivably be interdependent in varying ways across varying energetic regimes.

In particular, the dynamics and relationships between these key parameters in the intermediate $Fr = O(1)$ regime remains relatively uninvestigated in the literature. Motivated by oceanic and atmospheric flows, recent high-resolution body-forced numerical studies have predominantly focused on the ‘strongly stratified’ regime (Brethouwer et al. 2007; Maffioli et al. 2016; Portwood et al. 2016; Maffioli 2017; Taylor et al. 2019). Meanwhile, studies with temporally evolving simulations tend to traverse this regime temporarily between states of weak and strong stratification, or vice versa (Shih et al. 2005; Maffioli & Davidson 2016; Garanaik & Venayagamoorthy 2018) rather than in a forced quasi-stationary state. However, a key study by Portwood, de Bruyn Kops & Caulfield (2019) demonstrated that in stationary sheared stratified turbulence with a broad range of $Re_B$, and where $Fr$ and $Ri_g$ were allowed to evolve as free parameters, the flow ‘tuned’ to fixed values of $Ri_g \approx 0.16$, $Fr \approx 0.5$ and $\varGamma \approx 0.2$, independent of $Re_B$. In the context of the self-regulatory behaviour described previously, such results suggest that the $Fr = O(1)$ regime may have significant relevance across a variety of sheared stratified flows and warrants deeper investigation.

In light of the discussion presented above, we summarize the main concepts and subsequent open questions we aim to address within this study.

1. (i) In the context of our highly spatiotemporally inhomogeneous flow, can $\varGamma$ be accurately parameterized through instantaneous measurements of $Fr,Ri_g,L_E/L_O$ or $Re_B$?

2. (ii) Are these frameworks interconnected and if so how are they and the relationships between their relative parameters reconciled across the different mixing regimes?

3. (iii) What are the limitations on their applicability to open channel flow?

To that end, the remainder of the paper is structured as follows. In § 2 we present our flow configuration and numerical method. In § 3 we present a brief qualitative overview of our flows evolution and demonstrate the local parameter range of $Fr$ and $Re_B$ available within our flow. In § 4 we demonstrate the initial time-dependence of our flow on the development of the buoyancy field and demonstrate that after this period, the mixing properties within the flow become insensitive to global temporal effects and can be described by local processes. In § 5 we examine the applicability of $Fr, Ri_g, L_E/L_O$ and $Re_B$-based parameterization frameworks for both the mixing efficiency as well as the energetic state of the flow itself and subsequently derive relationships between all four non-dimensional parameters across the varying flow regimes within the channel. Finally, in § 6 we discuss our main findings within the study and their direct implications to the parameterization of mixing within stratified turbulence.

2. Numerical methodology

2.1. Problem formulation

We formulate our flow configuration in the same framework as Williamson et al. (2015) and their simulations of stationary radiatively heated open channel flow. A similar framework has been used for the destratifying open channel simulations of Kirkpatrick et al. (2019) and Kirkpatrick et al. (2020). A schematic of the flow is presented in figure 1. The flow is periodic in the streamwise ($x$) and spanwise ($y$) directions with no-slip and free-slip adiabatic bottom and top boundaries, respectively, and driven by a constant pressure gradient in $x$. The flow is heated through a depth varying volumetric heat source $q_I(z)$ on the principle of Beer–Lambert's law, defined as

(2.1)\begin{equation} q_I(z) = I_s\alpha {\rm e}^{(z-\delta)\alpha}, \end{equation}

where $I_s$ is the radiant surface heat flux, $\alpha$ is the absorption coefficient and $\delta$ is the channel height. Hence we can define domain-averaged mean heat source and normalized heat flux terms

(2.2a,b)\begin{equation} \langle q \rangle = \frac{1}{\delta}\int_0^{\delta}q_I(z)\,{\rm d}z, \quad q_N = \frac{1}{\delta^{2}}\int_0^{\delta}(\langle q \rangle-q_I(z))(z-\delta)\,{\rm d}z. \end{equation}

The dimensional temperature $\varTheta$ at time $t$ is decomposed into the statistically steady temperature fluctuation deviating from a domain-averaged mean, defined as

(2.3)\begin{equation} \varTheta(\boldsymbol{x},t) = \varTheta(\boldsymbol{x},t)' + \langle \varTheta(t) \rangle \end{equation}

and, under the assumption that no heat is lost through the boundaries, it follows that

(2.4)\begin{equation} \frac{\partial \langle \varTheta(t) \rangle}{\partial t} = \frac{\langle q \rangle}{\rho_0 C_p}, \end{equation}

where $\rho _0$ is the reference fluid density and $c_p$ is the specific heat. By defining the initial equilibrium friction velocity $u_{\tau,0}$ and channel depth $\delta$ as the characteristic velocity and length scales, respectively, we can then define a characteristic temperature scale

(2.5)\begin{equation} \varTheta_N = \frac{q_N\delta}{\rho_0 c_p u_{\tau,0}}. \end{equation}

Hence we can define a non-dimensional temperature and heat source

(2.6a,b)\begin{equation} \theta(\boldsymbol{x},t) = \frac{\varTheta(\boldsymbol{x},t) - \langle \varTheta(t) \rangle}{\varTheta_N}, \quad q(z) = \frac{q_I(z) - \langle q \rangle}{q_N}. \end{equation}

Our flow is then fully defined by four non-dimensional parameters: the initial friction Reynolds number $Re_{\tau,0}$; the Prandtl number $Pr$; the turbidity profile $\alpha \delta$; and an initial bulk stability parameter $\lambda _0$, defined as

(2.7a–d)\begin{equation} Re_{\tau,0} = \frac{u_{\tau,0}\delta}{\nu}, \quad Pr = \frac{\nu}{\kappa}, \quad \alpha \delta, \quad \lambda_0 = \frac{\delta}{\mathcal{L}}, \end{equation}

where $\nu$ and $\kappa$ are the kinematic viscosity and thermal diffusivity, respectively. Also note that in our framework, $\alpha \delta \Rightarrow \infty$ is analogous to a heat flux at the top surface similar to the simulations of Taylor et al. (2005). The stability of our flow is defined in the Monin–Obhukov framework as the ratio of the domain confinement scale $\delta$ to bulk Obhukov length $\mathcal {L}$ defined as

(2.8)\begin{equation} \mathcal{L} = \frac{u_{\tau,0}^{3}}{g\beta I_s/\rho_0 C_p} \left( \frac{1}{2} - \frac{1}{\alpha\delta}\right)^{{-}1}, \end{equation}

where $g$ is the gravitational acceleration and $\beta$ is the coefficient of thermal expansion. Note that $\mathcal {L}$ is analogous to the standard definition of the Obhukov length $L = u_{\tau }^{3}/\kappa b_*$ where $b_*$ is the surface buoyancy flux (Flores & Riley 2011). Thus we define our non-dimensional buoyancy scale as

(2.9)\begin{equation} b = \lambda_0 \theta. \end{equation}

We non-dimensionalize time around the initial advection time scale $T_{\tau,0} = \delta / u_{\tau,0}$ such that

(2.10)\begin{equation} t = T/T_{\tau,0}, \end{equation}

where $T$ is the dimensional time. Thus, under the Oberbeck–Boussinesq assumption, the non-dimensional governing equations for our flow are

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

(2.12)\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-} \boldsymbol{\nabla} p + \frac{1}{Re_{\tau,0}}\nabla^{2}\boldsymbol{u} + b\boldsymbol{e}_z +\boldsymbol{e}_x, \end{gather}

(2.13)\begin{gather}\frac{\partial b}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} b = \frac{1}{Pr Re_{\tau,0}}\nabla^{2}b + \lambda_0 q, \end{gather}

with our boundary and initial conditions explicitly defined as

(2.14)\begin{gather} z = 0; \quad u=v=w=0, \quad \frac{\partial b}{\partial z} = 0, \end{gather}

(2.15)\begin{gather}z = \delta; \quad \frac{\partial u}{\partial z}=\frac{\partial v}{\partial z}=w=0, \quad \frac{\partial b}{\partial z} = 0, \end{gather}

(2.16)\begin{gather}t = 0; \quad b = 0. \end{gather}

We have selected this open channel flow configuration for three reasons. Firstly, the use of an adiabatic bottom boundary has been shown to ensure that while the bulk flow becomes stratified, the near-wall turbulence structure remains relatively unchanged by the effects of buoyancy (Taylor et al. 2005; Williamson et al. 2015), allowing for simulations to be run at significantly higher levels of buoyancy strength than isoflux or fixed buoyancy boundary conditions where the relaminarization and collapse of turbulence inherently occurs at the wall (Flores & Riley 2011; García-Villalba & del Álamo 2011; Deusebio et al. 2015; Zhou et al. 2017a). Secondly, relative to a surface flux boundary condition, use of the volumetric heat source shifts the pycnolcine deeper into the channel away from the upper boundary, creating an appreciable region in the central bulk flow of significantly stronger stratification (Williamson et al. 2015). This subsequently allows us to access regimes of lower $Fr$ farther away from the top boundary, where confinement effects may influence mixing dynamics (Flores, Riley & Horner-Devine 2017). Thirdly, our direct numerical simulation (DNS) configuration of a stratified open channel flow heated through radiative surface heating is a canonical representation of stratified river flow and in particular the regulated river flows in inland Australia, where the accurate estimation and prediction of diapycnal mixing remains an important task (Williamson et al. 2015; Kirkpatrick et al. 2019).

Figure 1. Schematic diagram of the flow configuration; the domain is periodic in $x$ and $y$.

2.2. Notation and presentation of statistics

Within our temporally varying and vertically inhomogeneous flow, we define that any flow variable at a spatial location $\boldsymbol {x}$ and time $t$ can be decomposed into a horizontally averaged mean and fluctuating components denoted with an overbar and prime, respectively, such that

(2.17)\begin{equation} f(\boldsymbol{x},t) = \bar{f}(z,t) + f'(\boldsymbol{x},t) \end{equation}

and where the mean at a vertical location of $z$ is calculated through a volumetric integral across the horizontal plane at time $t$, as follows:

(2.18)\begin{equation} \bar{f}(z,t) = \frac{1}{L_x L_y} \int_0^{L_x} \int_0^{L_y} f(\boldsymbol{x},t) \,{\rm d}\kern0.7pt x\,{\rm d}y. \end{equation}

Similarly, we define that unless otherwise explicitly stated, it is implicit that flow statistics composed out of the velocity and buoyancy fields (e.g. $\epsilon _K,E_K,S,N$ etc) are presented as horizontal averages of that quantity at location $z$ and time $t$ as defined in (2.18).

2.3. DNS

Equations (2.11), (2.12), (2.13) were solved using a fractional-step finite-volume solver as outlined in Armfield et al. (2003) and Williamson et al. (2015). The advective spatial derivatives are discretized using fourth-order central differencing, whilst other spatial derivatives are calculated using second-order central differencing. Cell-face velocities are calculated using Rhie–Chow momentum interpolation and the time advancement is performed using a second-order Adams–Bashforth scheme.

We present our simulation parameters in table 1. Our simulations cover two friction Reynolds numbers, $Re_{\tau,0}=400,900$. As we are primarily interested in the low $Fr$ (high $Ri_g)$ regimes, our simulations cover a range of stability parameters, from moderately to extremely stable $(\lambda _0 = 0.5 - 5)$. We limit all simulations to $Pr = 1$ for numerical efficiency and set $\alpha \delta = 8$ for all simulations with a single control case of $\alpha \delta = 16$.

Table 1. List of DNS performed, and relevant parameters.

Following past studies of stratified wall-bounded turbulence (García-Villalba & del Álamo 2011; Deusebio et al. 2015; Williamson et al. 2015), we discretize our domain as follows. For all simulations the streamwise and spanwise grid size in initial viscous wall units is kept constant at $\Delta x_0^{+} = 5$ and $\Delta y_0^{+} = 2.5$. The vertical grid size for the $Re_{\tau,0}=400$ simulations is logarithmically stretched from $\Delta z_0^{+} = 0.4$ at the wall to $\Delta z_0^{+} = 4$ at $z = 0.25$ where it stays constant to the half-channel height $z = 0.5$. The vertical grid spacing in the top half of the channel is then set as symmetrical about the midpoint axis to ensure accurate resolution of viscous near-surface mechanics (Calmet & Magnaudet 2003). A similar procedure for the vertical grid size of the $Re_{\tau,0}=900$ simulations was employed with a further refinement of $\Delta z_0^{+} = 2.5$ in the bulk of the flow. To maintain accurate resolution of the viscosity affected near-wall and near-surface regions, we ensure that we have more than 10 grid points within a $\Delta z_0^{+} = 10$ distance from either boundary.

We keep the domain size constant at $L_x \times L_y \times L_z = 2 {\rm \pi}\delta \times {\rm \pi}\delta \times \delta$ across all simulations. We acknowledge that the size of the domain may affect the intermittent regime where laminar and turbulent patches coexist, as it has been shown that a smaller domain often leads to earlier laminarization for the same set of bulk parameters and can exhibit strong temporal intermittency (Flores & Riley 2011; García-Villalba & del Álamo 2011; Brethouwer, Duguet & Schlatter 2012; Deusebio et al. 2015). In this study, however, we are primarily interested in instantaneously correlating local turbulent mixing properties to appropriate local non-dimensional parameters within fully turbulent regions of the flow. Furthermore, our adiabatic bottom boundary condition ensures that the near-wall region remains fully turbulent, hence we do not expect the domain size to significantly influence the results presented in this study.

The initial simulation field is of fully developed statistically stationary neutral open channel flow at a given $Re_{\tau,0}$. For all simulations we then initialize the isothermal buoyancy field $b = 0$ at $t = 0$ and simultaneously switch on the volumetric heat source and the effects of buoyancy. All simulations have been run for $t_{{final}} = 10$ time units. For all simulations, transient data has been recorded at intervals of $\Delta t = 0.02$ time units to ensure accurate representation of temporal effects within the flow.

3. Flow overview

3.1. Qualitative description and local parameter range

A defining feature of our DNS configuration is that as the channel transitions from a neutral to stratified state, the local flow at any given vertical location $z$ evolves through a broad range of non-dimensional parameters such that, at any point in time, the flow contains regions of both high and low $Fr$ simultaneously. This creates a data set that traverses a variety of energetic regimes within a single simulation and where the flow, both mean and fluctuating, evolves in a relatively natural way without external imposition. We briefly present the local parameter range available within our data set through figure 2 showing the temporal evolution of $Fr$ and $Re_B$ as a function of $z$ for case R900L2 which shows behaviour typical of our flow, where we redefine

(3.1a,b)\begin{equation} Fr = \frac{\epsilon_K}{N E_K}, \quad Re_B = \frac{\epsilon_K}{N^{2} \nu}, \end{equation}

where $\epsilon _K = \nu (\partial u_i'/\partial x_j)^{2}$, $E_K = 1/2(\overline {u_i'u_i'})$ and $N = (\partial \bar {b}(z)/\partial z)^{1/2}$. From figure 2(a) we observe that the flow obtains an $Fr$ range that encapsulates all three energetic mixing regimes proposed by GV19. Of particular interest is that within the central portion of the channel, $Fr$ stabilizes within approximately one eddy turnover unit ($t \approx 1$) and obtains an appreciable depth range for the so-called ‘moderately stratified regime’ where $Fr = O(1)$ in an energetically ‘quasi-stationary’ state such that turbulent properties adapt rapidly to the evolving buoyancy gradient $N$ in order to obtain energetic equilibrium. As outlined in § 1, in previous numerical studies this regime is only often temporarily traversed as a transient state and subsequently there is a lack of data points in the available literature within this regime. As such our flow and extensive data set of $Fr = O(1)$ allows us to examine the mixing properties and scaling relationships within this regime in much greater detail than has been previously reported.

Figure 2. Evolution in time of key local parameters as a function of $z$ for case R900L2. (a) Turbulent Froude number $Fr$, vertical dotted lines left to right represent $Fr = 0.02$ (as the upper limit for the strongly stratified regime outlined in Lindborg 2006) and $Fr = 1$; (b) buoyancy Reynolds number $Re_B$, vertical dotted line represents $Re_B = 1$.

Although the flow achieves $Fr<0.02$ close to the free surface, our simulations are not able to access the ‘strongly stratified regime’ in the same sense as outlined in Billant & Chomaz (2001) or Brethouwer et al. (2007). This is made evident in figure 2(b), which demonstrates that for our flow configuration as we approach the free surface, the stratification grows stronger and a reduction in $Fr$ inevitably leads to a reduction in $Re_B$, leading to the relaminarization of the flow. Subsequently at our parameter range, regions within our flow of very low $Fr$ are more indicative of a diffusive regime as described by Brethouwer et al. (2007), where $Re_B <1$ and the flow is dominated by viscous effects and strong vertical shearing at large scales.

This is made clear if we consider the instantaneous visualizations of the buoyancy $b$ and enstrophy $|\boldsymbol {\omega }^{2}|$ fields in the vertical ($x,z$) plane at $t = 5$ in figure 3 for case R900L2. We can clearly observe that the flow is separated into three distinct regimes relative to the vertical coordinate $z$: a lower near-wall regime, where due to the adiabatic boundary condition the boundary layer turbulence structure remains relatively unchanged by the effects of stratification; a central region within which turbulent structures exhibit the characteristic ‘lean’ of sheared stratified turbulence and where we observe the formation of distinctly classic Kelvin–Helmholtz instability (KHI) overturning structures within the shear layer, causing highly vigorous and energetic mixing; an upper quasi-laminar or diffusive regime where turbulence is essentially suppressed by the effects of buoyancy and separated by a sharp buoyancy interface that experiences sporadic overturning by turbulent structures. Thus, the turbulence structure visually observed in figure 3 corresponds closely to the three regimes defined by $Fr$ and $Re_B$ in figure 2.

Figure 3. Instantaneous flow visualizations in the vertical ($x,z$) plane at $t = 5$ for case R900L2 of (a) the buoyancy field $b$ and (b) enstrophy field $|\boldsymbol {\omega }^{2}|$. Flow is moving left to right.

3.2. A note on the mixing efficiency

Throughout this study we defer to a definition of the instantaneous mixing efficiency through $R_f$ and $\varGamma$ as defined in Ivey & Imberger (1991),

(3.2a,b)\begin{equation} R_f = \frac{B}{B+\epsilon_K}, \quad \varGamma = \frac{R_f}{1-R_f} = \frac{B}{\epsilon_K}, \end{equation}

where $B = \overline {-b'w'}$ is the buoyancy flux. The vertical buoyancy flux, however, not only incorporates the small-scale irreversible mixing rate (the quantity of interest) but also large-scale reversible stirring processes (Caulfield & Peltier 2000; Peltier & Caulfield 2003). As such, for linearly stratified flows a more robust definition of the irreversible mixing rate is usually defined through the destruction rate of buoyancy variance $\chi$ (Maffioli et al. 2016; Venayagamoorthy & Koseff 2016; de Bruyn Kops & Riley 2019; Howland, Taylor & Caulfield 2020), where

(3.3)\begin{equation} \chi = \frac{\kappa}{N^{2}}\left(\overline{\frac{\partial b'}{\partial x_j}}\right)^{2}. \end{equation}

For $\chi$ to accurately represent the irreversible conversion of kinetic to available potential energy, a fundamental condition is that the local buoyancy period $N(z)$ must be invariant in both space and time (Caulfield 2020), a condition that is inherently unsatisfied within our temporally evolving inhomogeneous channel flow.

An alternative framework is through the adiabatic resorting of the buoyancy in the $z_*$ coordinate space of Winters et al. (1995) to obtain an irreversible mixing rate $\mathcal {M}$ and subsequently an irreversible mixing efficiency $\eta = \mathcal {M}/(\mathcal {M}+\epsilon _K)$ (Caulfield & Peltier 2000). In past studies, Zhou et al. (2017b) and Smith, Caulfield & Taylor (2021) have shown that the framework may be also applied to inhomogeneous shear flows to evaluate irreversible mixing across a midplane shear interface through appropriate spatial and temporal integration over the shear layer. However, in stratified open-channel flow where we are interested in investigating the correlation between mixing efficiency and local flow parameters across a broad range of vertical locations, rather than a single central shear layer, the $z_*$ framework presents obvious limitations. We note, the aim of this study is not to quantify an exact measure of the irreversible mixing efficiency, but rather to investigate its behavioural trends across different mixing regimes and the subsequent implications on the relationships between varying non-dimensional parameters at a given vertical location. Furthermore, as shown in Venayagamoorthy & Koseff (2016), we expect that in the weakly and moderately stratified regimes the differences in the definitions of mixing efficiency to be relatively small and the qualitative behaviour to remain similar. As such, our definition of mixing efficiency through (3.3) still accurately captures the dynamics of interest in our study.

4. Initial time dependence

We first briefly consider the initial time dependence exhibited by our flow properties related to the buoyancy field due to our idealized initial condition of $b=0$ by plotting $\varGamma$ as a function of time across a range of vertical locations and simulations in figure 4. From the results we observe that $\varGamma$ initially grows proportional to $t^{2}$ and shows clear time dependence up to approximately one eddy turnover time unit ($t \approx 1$). We consider the transport equation for the horizontally averaged mean buoyancy $\bar {b}$, which under the assumptions of homogeneity in $x$ and $y$, becomes

(4.1)\begin{equation} \frac{\partial \bar{b}}{\partial t} = \frac{\partial B}{\partial z} + \kappa \frac{\partial^{2} \bar{b}}{\partial z^{2}} + \lambda_0 q. \end{equation}

We can further take the vertical derivative $\partial / \partial z$ to obtain the evolution equation for $N^{2}$ as follows:

(4.2)\begin{equation} \frac{\partial N^{2}}{\partial t} = \frac{\partial^{2} B}{\partial z^{2}} + \kappa \frac{\partial^{3} \bar{b}}{\partial z^{3}} + \frac{\partial}{\partial z} \lambda_0 q. \end{equation}

For our simulations with initial condition $b=0$ we make the assumption that the turbulent and diffusive terms (first and second terms on the right-hand side) are negligible at the start of the simulation and the equation reduces to

(4.3)\begin{equation} \frac{\partial N^{2}}{\partial t} = \frac{\partial}{\partial z} \lambda_0 q. \end{equation}

Integrating forward in time from the initial reference time of $t = 0$ we arrive at an estimate of $N^{2}(t)$ at a given horizontal plane

(4.4)\begin{equation} \int_0^{t} \frac{\partial N^{2}}{\partial t} {\rm d}t = \int_0^{t} \frac{\partial}{\partial z}\lambda_0 q(z)\,{\rm d}t \quad \Rightarrow \quad N^{2}(t)= \frac{\partial}{\partial z}\lambda_0 qt. \end{equation}

Hence it is clear that initially $N^{2} \sim t$. For the remainder of this section, in the interest of simplification, we drop the notation $(t)$, however, the dependence on time of flow properties remains implicit. In the same framework as GV19, consider the turbulent vertical displacement of a fluid parcel $L_{disp} = w't$. The fluctuating buoyancy $b'$ can therefore be estimated as

(4.5)\begin{equation} b' \sim L_{disp}N^{2} = w'N^{2}t = w' \frac{\partial}{\partial z}\lambda_0 qt^{2} \sim t^{2}. \end{equation}

If we take one further assumption that buoyancy initially acts as a passive scalar, then it follows that flow properties related to the velocity field are not a strong function of $t$ and remain unchanged by the introduction of the buoyancy field. We can thus construct an initial time dependant expression for any flow property that incorporates the buoyancy field. Consider for example the buoyancy flux such that

(4.6)\begin{equation} B = \overline{-b'w'} \sim t^{2}. \end{equation}

By similar logic we can obtain an expression for $\varGamma$ such that

(4.7)\begin{equation} \varGamma = \frac{B}{\epsilon_K} \sim t^{2}, \end{equation}

as clearly demonstrated in figure 4 for all cases where $\varGamma \sim t^{2}$, until approximately one tenth of the characteristic eddy turnover time unit ($t \approx 0.1$), corresponding to the estimate for the time taken for energy injected at large scales to travel down the energy cascade to the dissipative range and hence affect the flow field (see Pope 2000). Past this time scale, buoyancy begins to affect the flow, nullifying our assumption of buoyancy acting as a passive scalar. Meanwhile we expect the turbulent and diffusive terms in (4.2) to become appreciable and influence the growth of $N^{2}$, causing it to diverge from a linear $N^{2} \sim t$ growth. This creates a set of complex dynamics, causing the buoyancy field to exhibit nonlinear time-dependence as the flow adjusts to the sudden imposition of buoyancy. This time dependence lasts of the order of one eddy turnover time unit ($t \approx 1$) across all simulations and vertical locations, suggesting the adjustment to the buoyancy field is a global rather than local process. For $t \gtrsim 1$ temporal variability becomes negligible and $\varGamma$ begins to evolve at a quasi-steady rate. We note similar dependence on the initial eddy turnover time scale has been observed in previous studies with distinctly different flow configurations, yet with a similar $b=0$ initial condition (Métais & Herring 1989; Venayagamoorthy & Stretch 2006; Maffioli & Davidson 2016), suggesting some universality on the eddy turnover time scale for the nonlinear adjustment of the flow. Subsequently we can expect that for $t \gtrsim 1$, $B$ and hence $\varGamma$ become independent of time and evolve relative to local processes. Thus it follows that for $t \gtrsim 1$ we can expect locally based parameterization schemes for $\varGamma$ to become applicable to our flow.

Figure 4. Mixing coefficient $\varGamma$ plotted against non-dimensional time $t$ (a) at varying vertical locations for case R900L2 (b) at a vertical location of $z = 0.5$ for all simulations. Vertical dashed lines in both figures represent $t= 0.1$ and $t = 1$; the diagonal dashed line represents a line proportional to $t^{2}$.

5. Parameterization of mixing efficiency, applicability and comparison

We now turn to the main theme of this study. In this section we explore the parameterization of mixing efficiency in our temporally evolving inhomogeneous flow through the $Fr, Ri_g, L_E/L_O$ and $Re_B$ methods. We then use the $Fr$-based framework of MBL16 and GV19 as a base case scenario to further investigate the relationships and similarities between the varying frameworks and relevant non-dimensional parameters across the varying mixing regimes.

5.1. $Fr - \varGamma$ framework and moderately stratified regime scaling

To begin we first explicitly test the novel $\varGamma - Fr^{-1}$ scaling relationship derived in GV19 for the $Fr = O(1)$ regime. The authors argue that the pertinent time scales for $B$ and $\epsilon _K$ are the buoyancy and inertial time scales $T_N = 1/N$ and $T_L = E_K/ \epsilon _K$, respectively. Under the assumption of a stationary linearly stratified flow they obtain

(5.1a,b)\begin{equation} B \sim \chi \sim \frac{w'^{2}}{T_N}, \quad \epsilon_K \sim \frac{w'^{2}}{T_L}. \end{equation}

Thus they obtain the new scaling for the irreversible mixing coefficient $\varGamma = B/\epsilon _K \sim \chi /\epsilon _K \sim T_L/T_N = Fr^{-1}$. For the purpose of generality we make one further assumption, that at $Fr = O(1)$ the separation of vertical and horizontal velocity scales is still relatively small and hence we can approximate $w'^{2} \sim E_K$. We can thus rewrite (5.1) in the classic form of velocity and length scales $\epsilon \sim u^{3}/l$ as

(5.2a,b)\begin{equation} B \sim \chi \sim \frac{E_K}{T_N} \sim \frac{E_K^{3/2}}{L_{N}}, \quad \epsilon_K \sim \frac{E_K}{T_L} \sim \frac{E_K^{3/2}}{L_{I}}, \end{equation}

where

(5.3a,b)\begin{equation} L_{N} = \frac{E_K^{1/2}}{N}, \quad L_{I} = \frac{E_K^{3/2}}{\epsilon_K}. \end{equation}

Here $L_{N}$ is an energetic buoyancy length scale constructed as a result of dimensional analysis and can be taken to represent the conceptual size of large energy-containing eddies when the effects of buoyancy are significant. Now $L_N$ has also been shown to be an accurate indicator of the size of overturns when $Fr < 1$ (Mater et al. 2013). Here $L_{I}$ is the well known inertial energy-containing turbulent length scale (Pope 2000). Such that for the moderately stratified regime we obtain $\varGamma = B/\epsilon _K \sim \chi /\epsilon _K \sim L_{I}/L_{N} = Fr^{-1}$, the same scaling as GV19. If considered in the context of the strongly stratified turbulence theory of Billant & Chomaz (2001) and Lindborg (2006) it can be readily seen that $L_{N}$ and $L_{I}$ have physical analogues in the vertical and horizontal integral scales $l_v \sim u_h/N$ and $l_h \sim u_h^{3}/\epsilon _K$, respectively, where $u_h$ is the horizontal velocity scale.

As discussed in § 3, the definition of $\chi$ as an irreversible mixing rate in the derivation above is not equivalent to that of an instantaneous local $\chi (z)$ as measured within our flow, hence, and without loss of generality, we explore these scaling proposals through the buoyancy flux $B$ instead. Figure 5(a) shows $B$ normalized by $E_K^{3/2} /L_{N}$ as function of time at a vertical location of $z = 0.5$, corresponding to a location in the flow where $Fr = O(1)$ and $Re_B \gg O(1)$ for all simulations. For the proposed scaling to hold, we expect the former ratio to approach a constant value of $O(1)$. We further note that it is sufficient to show the $B$ scaling alone to validate the $\varGamma \sim Fr^{-1}$ assumption; as it can readily be shown that $\epsilon _K / (E_K^{3/2}/L_{I}) = (\epsilon _K /E_K^{3/2}) /(\epsilon _K /E_K^{3/2}) = 1$.

Figure 5. Buoyancy flux $B$ normalized by $E_K^{3/2} /L_{N}$, (a) as a function of time $t$ for all simulations at a vertical location of $z = 0.5$, (b) as a function of the turbulent Froude number $Fr$ presented as a two-dimensional probability density function (p.d.f.) for $t>1$ and $0.2 \leqslant z \leqslant 0.8$. Vertical dashed lines indicate $Fr = 0.3$ and $Fr = 1$. Dashed horizontal line in both figures indicates an empirical constant of 0.08.

In agreement with our analysis in § 4 we observe that the scaling does not hold for $t < 1$ during the transient adjustment period, past which the results clearly show an agreement with the scaling in (5.2) such that $B / (E_K^{3/2} /L_{N})$ approaches a constant value of approximately $0.08$ across all simulations and appears invariant in time. When considered in the context of our temporally inhomogeneous flow where not only the turbulent properties but also the background stratification $N^{2}$ is evolving in time, the result implies that the above scaling for $B$ in (5.2) is both valid and relatively robust. This is further made evident in figure 5(b) which shows the two-dimensional p.d.f. of $Fr$ and the ratio $B / (E_K^{3/2} /L_{N})$ for all simulations. Due to the initial time dependence observed in § 4 we exclude the data for $t<1$ in the construction of the p.d.f. Furthermore, as the confinement effects of the top and bottom boundaries are outside the scope of this study we also exclude the data for $z<0.2$ and $z>0.8$. The vertical limits of $z=0.2$ and $z=0.8$ have been chosen as the approximate values within which fully developed open channel flow has been shown to obtain local energetic equilibrium (Williamson et al. 2015; Kirkpatrick et al. 2019). It is indeed clear from the results that $B / (E_K^{3/2} /L_{N}) \approx 0.08$ and remains constant only for a regime within the bounds of $0.3 \lesssim Fr \lesssim 1$, presenting further strong evidence for the argument of GV19 for the existence of a separate intermediate energetic regime that exists in a quasi-stationary state and not as a transient transition between the weakly and strongly stratified regimes. We note that the value of $Fr = 0.3$ is in direct agreement with the transitional value observed in MBL16 towards the asymptotic $\varGamma$ regime within their study, suggesting some level of universality to this transitional value.

We now turn our attention to the $Fr$-based parameterization of the mixing efficiency. Figure 6(a) shows the two-dimensional p.d.f. of $Fr$ and $\varGamma$ for all simulations constructed out of the data within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$.

Figure 6. (a) Two-dimensional p.d.f. of the turbulent Froude number $Fr$ and the mixing coefficient $\varGamma$ constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. The axes on the insert within the figure are presented on a linear scale. (b) The $Fr$ bin-averaged mixing coefficient $\langle \varGamma \rangle$ plotted against bins of corresponding turbulent Froude number $\langle Fr \rangle$ for all data points within $t>1$ and $z>0.2$. Solid lines indicate the proposed scaling of MBL16 and GV19 as well as empirically observed $\varGamma = 0.3$. Vertical dashed lines indicate $Fr = 0.3$ and $Fr = 1$.

From the results we observe that $\varGamma$ collapses well across all three energetic regimes and respective scaling relationships proposed by MBL16 and GV19. Although we only access a very small section of the weakly stratified regime, we observe a distinct $Fr^{-2}$ slope developing in the high $Fr$ range similar to that observed in the Couette flow studies of Zhou et al. (2017a,b). Crucially, in agreement with our analysis above, we observe that within a range of $0.3 \lesssim Fr \lesssim 1$ where the majority of our data is concentrated, $\varGamma$ collapses across all simulations and depths with a distinct $Fr^{-1}$ scaling, providing further strong evidence for the existence of a distinctly separate intermediate regime as proposed by GV19. Unlike MBL16, however, who find that $Fr = 0.3$ corresponds to a peak in mixing efficiency, we observe no non-monotonic behaviour in $\varGamma$, rather the transition at $Fr = 0.3$ corresponds to a saturation of the mixing efficiency, with $\varGamma$ displaying independence of $Fr$ and trending towards an empirically observed asymptotic value of $\varGamma \approx 0.3$. The asymptotic value being significantly higher than the upper limit of $0.2$ proposed by Osborn (1980), is more reminiscent of the higher values of cumulative mixing efficiency seen in studies of KHI-induced mixing (Salehipour & Peltier 2015; Mashayek, Caulfield & Peltier 2017; Salehipour et al. 2018), seemingly in agreement with our visual observations of highly energetic KHI-driven mixing events within our flow.

With further decreasing $Fr$ we observe a significant amount of scatter in the data. This can be explained as the flow transitions towards the diffusive regime and both $B$ and $\epsilon _K$ grow very small and a minor change in either of the two quantities causes large variation in $\varGamma$. Additionally, as has been reported in the study García-Villalba & del Álamo (2011), channel flow within this regime becomes increasingly affected by the propagation of internal gravity waves as well as counter-gradient fluxes resulting from convective instabilities in the flow (Taylor et al. 2005). This results in instantaneous measurements of $B$ becoming strongly susceptible to contamination from these reversible processes, subsequently causing our measurements of $\varGamma$ within this regime to show significant scatter about its asymptotic mean value. Furthermore, as discussed in Venayagamoorthy & Koseff (2016) the presence of the counter-gradient fluxes may cause the observed saturated mean value of $\varGamma = 0.3$ to underestimate the true irreversible mixing efficiency. However, as can be observed from the insert in figure 6(a) which presents the same data on a linear–linear axes, the observed frequency of counter-gradient fluxes is relatively negligible in comparison with the full data set. As such the significant takeaway from the results presented within this regime is not the precise value for the saturated mixing efficiency but rather the observation that $\varGamma$ appears to grow independent of $Fr$ within this regime.

To investigate this further we sort and average our instantaneous measurements of horizontal plane measurements of $\varGamma (z,t)$ across bins of constant $Fr$ to obtain the binned data set of $\langle \varGamma \rangle$. Avoiding near-wall mechanics we consider all data for $z>0.2$ and $t>1$. We note that for this case we include the data near the free-surface boundary to highlight the departure point from the saturated regime towards the diffusive regime. Figure 6(b) shows $\langle \varGamma \rangle$ plotted against bins of corresponding $\langle Fr \rangle$ for all simulations under the conditions described above. A colour bar depicting $\langle Re_B \rangle$ is also shown for reference. When presented in this manner, it becomes clear that for $Fr \lesssim 0.3$ the mixing efficiency indeed saturates towards a constant value for all simulations, independent of $Fr$, a detail that is somewhat less clear in the scattered data presented in figure 6(a). Furthermore, it is clear that the transition from a saturated mixing efficiency to the $\varGamma \sim Fr^{-1}$ regime at $Fr \approx 0.3$ occurs irrespective of $Re_{\tau }, \lambda, \alpha \delta$ or vertical location $z$ within the channel, again suggesting a degree of universality for this transitional value as argued by MBL16. Conversely no singular value of $Fr$ appears for the transition away from the constant $\varGamma$ regime within the ‘left flank’ of the curve. Rather, as seen from the results, this corresponds to a transition to the diffusive regime at low $Re_B$ as will be demonstrated in further detail in § 5.5.

We note that within the ‘left flank’ of figure 6(b), following the saturated regime, case R900L5 displays a further increase in $\varGamma$ with a clear peak before dropping away. The exact mechanism of this is unknown, and may be attributed to a number of influencing factors that fall outside the scope of this study such as large vertical displacement of fluid parcels through internal gravity waves, surface confinement effects and strong spatiotemporal intermittency of turbulence as the flow approaches the diffusive regime. Furthermore, as can be observed from the colour bar, this behaviour occurs for flows where $Re_B \ll O(1)$, which can be considered essentially laminar such that the vertical transfer of buoyancy through the turbulent flux $B$ is significantly smaller than that of molecular diffusion. Subsequently, our definition of a ‘turbulent’ mixing efficiency through $\varGamma$ begins to lose relevance. We further note that similar behaviour has been observed at low $Re_B$ in the studies of Zhou et al. (2017b) and Smith et al. (2021), suggesting that as $\epsilon _K$ goes to zero the diapycnal flux is not fully suppressed within this regime.

5.2. $Ri_g$ framework

In this section we briefly examine the behaviour of $\varGamma$ relative to the more established framework based on $Ri_g$, where

(5.4)\begin{equation} Ri_g = \frac{N^{2}}{S^{2}}, \end{equation}

where $S = \partial \bar {u}(z)/\partial z$. Figure 7 shows the two-dimensional p.d.f. of $Ri_g$ and $R_f$ constructed analogously to that of figure 6(a). For reference we also plot the line corresponding to the linear relationship $R_f = Ri_g$ (dotted line) as well as the empirical fit corresponding to $R_f = 0.25(1-{\rm e}^{-7Ri_g})$ of Venayagamoorthy & Koseff (2016) (solid lined). From the results it is clear that, in agreement with numerous studies of sheared stratified turbulence (see § 1), the mixing efficiency displays linear monotonic dependence on $Ri_g$ up to a critical value of $Ri_g \approx 0.25$, where it can be seen that the majority of our data is concentrated. The results add to the mounting evidence that in weak/moderate levels of stratification, the diapycnal and momentum diffusivities are relatively equal, such that the turbulent Prandtl number $Pr_T \sim Ri_g/R_f$ is approximately unity. A further key observation from the results is that within our flow we clearly observe sustained turbulence and vigorous KHI-induced mixing at $Ri_g$ values significantly higher than the critical ‘Miles–Howard’ limit of $Ri_{g,c} = 0.25$, suggesting exercising caution in assuming a strict upper limit of KHI stability based on $Ri_{g,c}$ for more complex and initially turbulent flows, as suggested by Galperin et al. (2007).

Figure 7. Two-dimensional p.d.f. of the gradient Richardson number $Ri_g$ and the flux Richardson number $R_f$ constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. Solid line indicates the proposed empirical fit model of Venayagamoorthy & Koseff (2016). Dotted line represents a linear relationship of $R_f = Ri_g$. Vertical dashed line indicates $Ri_g = 0.25$. The axes on the insert within the figure are presented on a linear scale.

Within the ‘right flank’ of the curve we again observe no non-monotonic behaviour in $R_f$ as observed in the high $Pr$ Couette flow simulations of Zhou et al. (2017b), rather the saturation of the mixing efficiency seems very well described by the empirical fit of Venayagamoorthy & Koseff (2016) showing clear asymptotic behaviour in the limit of high $Ri_g$. We note in our study that we maintain $Pr = 1$ for all simulations, rather than higher values of $Pr = 6 - 7$ which are typical of stably stratified river or oceanic flows that provide the motivation behind this study. Past studies of shear instability driven mixing (Smyth et al. 2001; Salehipour & Peltier 2015; Salehipour, Peltier & Mashayek 2015) have shown that at higher $Pr$ the dynamics of a KHI mixing event can vary significantly to that at lower $Pr$, particularly the secondary instabilities that form at smaller scales. Furthermore, we note our definition of the mixing efficiency through $R_f$ has fundamentally different behaviour at high $Ri_g$ compared with its irreversible counterpart used in the aforementioned studies, making a direct comparison somewhat difficult. It remains unclear how higher $Pr$ values for our flow configuration would affect the results presented in this study, and presents direction for future work.

5.3. Influence of mean shear and $Fr - Ri_g$ scaling

A distinct key observation from the above analysis is that the transition to the saturated mixing efficiency regime is described well by both $Fr \approx 0.3$ and $Ri_g \approx 0.25$. However, for $0.3 \lesssim Fr \lesssim 1$ and $Fr \gtrsim 1$ we have observed two distinctly different mixing regimes with separate scaling for $\varGamma$, in contrast to a simple linear dependence of the mixing efficiency on the gradient Richardson number for $Ri_g \lesssim 0.25$. Such results imply a more complex relationship between $Fr$ and $Ri_g$ within our flow than the $Ri_g \sim Fr^{-2}$ relationship for weakly stratified flow (Zhou et al. 2017a). To investigate this further, we consider that the dynamics of our channel flow can be described in the conceptual framework of Mater & Venayagamoorthy (2014) as a competition of the inertial, buoyancy and shear forces within the flow, or analogously as a competition of the three respective time scales $T_L, T_N, T_S$, where $T_S= 1/S$ is the shear time scale. Thus we consider our channel flow within their $S^{*} - Fr^{-1}$ regime map, where $S^{*}$ is the non-dimensional shear rate defined as

(5.5)\begin{equation} \textit{Note that }S^{*} = \frac{S E_K}{\epsilon_K}.\end{equation}

$S^{*}$ has been frequently used in the literature to describe sheared turbulence, both with and without the presence of buoyancy (Rogallo 1981; Shih et al. 2000; Chung & Matheou 2012). Furthermore, It can be readily seen that $S^{*} = T_L/T_S$, and hence $S^{*}$ represents the competition between the shear and inertial time scales in the same manner that $Fr = T_N/T_L$ analogously represents the competition of buoyancy and inertial time scales. We note that by convention, $S^{*}$ is often defined by $S^{*} = Sq/\epsilon _K$ where $q = 2E_K$ is twice the turbulent kinetic energy. For this study, however, we defer to the definition in (5.5) to maintain consistency in the comparison of the respective three time scales.

In this sense, Mater & Venayagamoorthy (2014) propose that stratified shear flow can be divided into three distinct regimes. An inertia dominated regime where flow tends to revert to isotropy and is defined by both $S^{*}< S_c^{*}$ and $Fr^{-1} < Fr_c^{-1}$, where $S_c^{*}$ and $Fr_c$ are some critical values at which shear and buoyancy forces become significant relative to inertia and begin to distort the isotropic flow. Meanwhile the buoyancy and shear dominated regimes are defined by $S^{*}> S_c^{*}$ and $Fr^{-1} > Fr_c^{-1}$ and separated by the diagonal line of $Ri_g = 0.25$ as it can be readily shown that $Ri_g =(T_S/T_N)^{2} = (Fr^{-1}/S^{*})^{2}$. We adopt the critical values defined in Mater & Venayagamoorthy (2014) of $S_c^{*} = 3.3$ and $Fr_c^{-1} = 1.65$. Furthermore, we note that the limit of $Ri_g = 0.25$ delineating the shear and buoyancy dominated regime is not a strict limit as in the sense of the ‘Miles–Howard’ criterion but rather an empirical choice due to the evidence that in forced stratified sheared turbulence $Ri_g$ tends to an upper bound of $0.25$ in a stationary state (Rohr et al. 1988; Chung & Matheou 2012).

Figure 8 shows the two-dimensional p.d.f. of our flow within the $S^{*} - Fr^{-1}$ regime map described above, constructed in the same manner as figure 6. From the results we can again observe three distinct regimes of behaviour separated by $Fr = 1$ and $Fr = 0.3$.

Figure 8. Two-dimensional p.d.f. of the inverse of the turbulent Froude number $1/Fr$ and non-dimensional shear rate $S^{*}$ in the $ST_L - NT_L$ regime map of Mater & Venayagamoorthy (2014) constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. Solid lines indicate the separation of the proposed inertia, buoyancy and shear dominated regimes as outlined in § 5.2. Vertical dashed lines indicate $1/Fr = 1$ and $1/Fr = 3.33$ corresponding to $Fr = 1$ and $Fr = 0.3$, respectively.

In the weakly stratified regime for $1/Fr < 1$, buoyancy acts essentially as a passive scalar and the flow travels along a horizontal path of $S^{*} \approx S_c^{*}$ where the shear and inertial forces are in balance, much like in a neutral channel flow.

Within the moderately stratified regime we observe a transitional regime where buoyancy begins to change the dynamics of the flow. As $Fr$ decreases and buoyancy begins to restrict the vertical velocity fluctuation $w'$ and subsequently the turbulent shear stress $\overline {-u'w'}$, an imbalance appears in the flow between the total local shear stress and the driving pressure gradient. This causes the flow to accelerate in an effort to obtain energetic equilibrium, translating into an increase in a local measure of $S$, such that $S^{*}$ increases accordingly. Thus, it is clear that for our flow this intermediate regime where both $S^{*} = O(1)$ and $Fr = O(1)$ represents a critical state of the flow where inertia, buoyancy and shear all play a significant and interconnected role in the dynamics of the flow.

Conversely in the saturated regime for $Fr < 0.3$ or $Ri_g > 0.25$ we observe that with further increasing stratification, the shear increases accordingly such that $S^{*}$ grows large and $Fr$ grows small. Within this regime we can thus expect inertia to become insignificant relative to the effects of shear and stratification such that $N \sim S$ and $Fr$ becomes decoupled from $Ri_g$ as the effects of inertia become negligible relative to buoyancy and shear. Accordingly, although the data in this regime shows some scatter we can observe that within this regime the flow tends to settle down and evolve along diagonal lines of constant $Ri_g$, with the majority of the data in the ‘right flank’ of the figure again being concentrated around the ubiquitous value of $Ri_g = 0.25$. The scatter in the data can be directly explained by the relatively slow process that is the acceleration of the mean flow from its initial state. We note that at $t_{{final}} = 10$ all the simulations are still significantly distant from their final equilibrium states. As such the mean shear $S$ is still increasing to obtain shear stress equilibrium such that $Ri_g$ is still evolving towards its stationary state. Meanwhile, the turbulent properties have adapted to the growing background stratification and $Fr$ has essentially stabilized as observed in figure 2, causing the scattered trajectories of the flow in $S^{*} - Fr^{-1}$ space within this regime.

From the observations above it is clear that the relationship between the inertial, shear and buoyancy forces within the flow vary significantly along the three energetic regimes of the flow. We can hence derive scaling arguments to directly relate $Fr$ to $Ri_g$ within our flow across the three regimes.

Within the weakly stratified, $Fr > 1$ regime we can relate $Fr$ and $Ri_g$ with a simple mixing length argument for $S$ such that

(5.6)\begin{equation} S \sim \frac{U_*}{L_*},\end{equation}

where $U_*$ and $L_*$ are some characteristic velocity and length scales pertinent to energetic mixing within the flow. Observing a balance between shear and inertial forces we can estimate $U_* \sim E_K^{1/2}$ and $L_* \sim L_{I}$ such that $S \sim \epsilon _K/E_K$ or conversely $T_S \sim T_L$. Thus we can obtain

(5.7)\begin{equation} Ri_g = \frac{N^{2}}{S^{2}} \sim \frac{N^{2} E_K^{2}}{\epsilon_K^{2}} = Fr^{{-}2};\end{equation}

the same scaling as derived in MBL16 or Zhou et al. (2017a).

Within the moderately stratified, intermediate $0.3< Fr<1$ regime we can adopt a similar mixing length argument for $S$ as in (5.6). Having hypothesized that in this critical regime inertia, buoyancy and shear are all significant, the dimensional group that could influence the dynamics becomes $(N,S,E_K,\epsilon _K)$. Analogously to our scaling of $B$ within this regime in (5.2), we make the assumption that the pertinent length scale is the energetic vertical buoyancy length scale such that $L_* \sim L_{N}$. Hence, dimensional analysis dictates that a velocity scale constructed out of the remaining dimensional quantities becomes $U_* \sim (\epsilon _K/N)^{1/2}$, where $(\epsilon _K/N)^{1/2}$ is the velocity scale analogue of the Ozmidov length scale (Shih et al. 2005). Hence we obtain

(5.8)\begin{equation} S \sim \frac{U_*}{L_*} \sim \frac{(\epsilon_K/ N)^{1/2}}{E_K^{1/2} /N} = \left(\frac{\epsilon_K N}{E_K}\right)^{1/2}. \end{equation}

Conversely this can be seen as a comparison of time scales such that

(5.9)\begin{equation} T_S \sim (T_L T_N)^{1/2} \end{equation}

and we can obtain a scaling relation between $Fr$ and $Ri_g$ within this regime

(5.10)\begin{equation} Ri_g = \frac{N^{2}}{S^{2}} \sim \frac{N^{2} E_K}{N \epsilon_K} = \frac{N E_K}{\epsilon_K} = Fr^{{-}1}. \end{equation}

We note the scaling proposed in (5.8)–(5.10) are new scaling relationships for the $Fr=O(1)$ regime and thus directly reconcile the observed $\varGamma \sim Fr^{-1}$ and $R_f \sim Ri_g$ scaling for mixing efficiency within this regime.

Finally, in the $Fr < 0.3$ saturated regime we can directly imply that as $T_N$ and $T_S$ both become much smaller than $T_L$, the effects of inertia become negligible such that $N \sim S$ and $Ri_g$ becomes decoupled from $Fr$.

We first test the new scaling in (5.9) explicitly by plotting the ratio $S / (\epsilon _K N /E_K)^{1/2}$ as a function of time at $z = 0.5$ in figure 9(a) and as a function of $Fr$, presented in the form of a two-dimensional p.d.f., in figure 9(b). From the results it is clear that $S / (\epsilon _K N /E_K)^{1/2}$ approaches a constant value of approximately $0.3$ that appears invariant with respect to both time and space within the defined regime of $0.3 \lesssim Fr \lesssim 1$. Again in the context of our evolving and inhomogeneous flow, such results suggest that the scaling is reasonably robust and may pertain to a wider range of stratified shear flows.

Figure 9. Mean shear rate $S$ non-dimensionalized by $(\epsilon _K N /E_K)^{1/2}$ (a) plotted against time $t$ for all simulations, horizontal dashed lined indicates empirically observed constant of 0.3, (b) presented in the form of a two-dimensional p.d.f. with the turbulent Froude number $Fr$. The p.d.f. is constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. Vertical dashed lines indicate $Fr = 0.3$ and $Fr = 1$.

To confirm the different scaling across all regimes, we plot the two-dimensional p.d.f. of $Fr$ and $Ri_g$ in figure 10(a) for our DNS results as well as from the stationary homogeneous sheared turbulence studies of Shih et al. (2000), Chung & Matheou (2012) and Portwood et al. (2019). From our instantaneous results it is clear that the flow is again well divided into the three distinct regimes divided by $Fr = 1$ and $Fr = 0.3$, with clear $Ri_g \sim Fr^{-2}$ and $Ri_g \sim Fr^{-1}$ behaviour in the weakly and moderately stratified regimes, respectively. The results of the three aforementioned studies with distinctly different forcing mechanics also show clear support for the derived scaling showing excellent agreement across both regimes, suggesting a degree of universality in their application. Meanwhile the transition to the saturated regime at $Fr \approx 0.3$ and $Ri_g \approx 0.25$ corresponds to a decoupling between $Fr$ and $Ri_g$ as the variables become uncorrelated, although this is somewhat obscured by the scatter in the data.

Figure 10. (a) Two-dimensional p.d.f. of the turbulent Froude number $Fr$ and the gradient Richardson number $Ri_g$ constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. (b) The $Fr$ bin-averaged gradient Richardson number $\langle Ri_g \rangle$ plotted against bins of corresponding $\langle Fr \rangle$ for all data points within $t>1$ and $z>0.2$. Large blue circles show the data of the stationary runs in table 2 of Shih et al. (2000). Large diamonds (cyan in (a), black in (b)) show data of Chung & Matheou (2012). Large green ‘X’ shows ‘tuned’ values of Portwood et al. (2019). Solid lines indicate the proposed scaling in (5.7) and (5.10). Vertical dashed lines indicate $Fr = 0.3$ and $Fr = 1$. Horizontal dashed lines indicate $Ri_g = 0.25$.

To demonstrate this result more clearly we plot the $Fr$ bin-averaged data set of $\langle Ri_g \rangle$ plotted against corresponding bins of $\langle Fr \rangle$ in figure 10(b) in a similar manner to that of the process described for figure 6(b). A colour bar showing $\langle Re_B \rangle$ is again included for reference. From the results it is again clear that the scaling derived for the weakly and moderately stratified regimes in (5.7) and (5.10) distinctly collapses on one line and holds irrespective of external flow parameters of vertical location in the channel, with a clear transition at $Fr \approx 0.3$ or $Ri_g \approx 0.25$ at which point each simulation continues its own unique trajectory in $Fr-Ri_g$ space. It is also worth noting the long individual ‘tails’ in the ‘left flank’ of the curve, where $Ri_g \gg O(1)$ correspond to low $Re_B$ flow that is essentially laminar and where the scaling derived above on the assumption of turbulent flow becomes invalid. Furthermore, the data of the aforementioned three studies is again plotted for visual clarity to highlight the validity of the scaling.

In addition to the stationary runs plotted in figures 10(a) and 10(b) in which $Ri_g$ is free to evolve to its stationary state, Shih et al. (2000) also ran non-stationary simulations (see table 1 in Shih et al. (2000)) in which $Ri_g$ and subsequently the mean gradients $S,N$ are kept fixed for a given simulation. These simulations do not follow the proposed scaling but instead evolve along constant lines of $Ri_g$ across all three regimes. Conversely, in our study and the other data sets shown in figures 10(a) and 10(b), neither $Fr$ or $Ri_g$ are known a priori and subsequently the both the mean and turbulent flow properties adapt to the proposed scaling.

The findings discussed in the analysis above present implications for the parameterization of stratified shear flow. Firstly, it is clear that for $Fr >0.3$ or $Ri_g<0.25$, in the weakly and moderately stratified regimes, $Fr$ and $Ri_g$ become interchangeable in any parameterization scheme by using the scaling relations in (5.7) and (5.10), respectively. Such interchangeability allows for more flexible parameterization as relationships derived through turbulent properties and $Fr$ may be inferred in real flows from relatively simple measurements of the mean buoyancy and velocity field. Secondly, although very appealing, the simple division of stratified shear flow into two regimes along a line of $Ri_g = 0.25$ does not accurately capture the subtleties and the differences in dynamics between the distinctly different $0.3< Fr<1$ and $Fr>1$ regimes. Lastly, care should be taken in any assumptions on the state of the flow by inferring relationships between $Ri_g$ and $Fr$ in the $Fr < 0.3$ regime; in particular in temporally evolving flow where the mean fields may exhibit appreciably long adjustment periods as the flow transitions to energetic equilibrium. Furthermore, in the context of the strongly stratified turbulence theory of Brethouwer et al. (2007), the apparent decoupling of $Fr$ and $Ri_g$ within this regime suggests that even if an upper limit fundamentally exists on the stationary value of $Ri_g$, it does not translate to a lower bound on $Fr$, suggesting the possibility of accessing the strongly stratified $Fr \ll O(1)$ regime within stratified channel flow.

5.4. Overturning length scale framework

In this section we investigate the applicability of the $L_E/L_O - \varGamma$ framework and hence the more easily measurable $L_T/L_O - \varGamma$ framework of GV19 to our temporally evolving channel flow. Where we define

(5.11a,b)\begin{equation} L_E = \frac{b'_{{rms}}}{N^{2}}, \quad L_O = \left(\frac{\epsilon_K}{N^{3}}\right)^{1/2}, \end{equation}

where $L_E$ is the well known Ellison overturning length scale and represents the large energy containing overturns within the flow. Meanwhile the Ozmidov length $L_O$ represents the maximum conceptual size of an isotropic eddy that is not confined by stable stratification (Smyth & Moum 2000). In this section we use the overturning length scales $L_E$ and $L_T$ interchangeably when referencing different studies, due to their linear relationship in fully turbulent flow as shown in Mater et al. (2013).

Figure 11 shows the two-dimensional p.d.f. of $L_E/L_O$ and $\varGamma$ for the same data set as in figure 6. It is clear that $\varGamma$ is accurately described through instantaneous measurements of $L_E/L_O$ across all three regimes and scaling relationships described in GV19. Although accessing only a very small section of the regime, we observe a $\varGamma \sim (L_E/L_O)^{4/3}$ relationship in the weakly stratified regime in agreement with the observational oceanic study of Ijichi & Hibiya (2018). However, as $L_E/L_O$ becomes increasingly small we note that the results appear to deviate slightly from the proposed scaling for the weakly stratified regime. We shall return to this shortly in the analysis to come. Again, we observe a distinct collapse of the data in the moderately stratified regime, showing a $\varGamma \sim (L_E/L_O)^{1}$ scaling that holds for almost an entire decade of $(L_E/L_O)$. A key observation is that the transition towards the saturated regime occurs at precisely $L_E/L_O \approx 1$. Taking the approximation of $L_E \approx L_T$ and considering the visual observation of vigorous KHI-induced mixing in figure 3, this is conceptually consistent with the work of Mashayek et al. (2017) who found that mixing efficiency in a KHI mixing event peaks when $L_T \sim L_O$. That is, when energy being injected into the downscale energy cascade through the overturning of the KHI is at a scale corresponding to the upper end of the inertial subrange such that it is not constricted by the mean stratification. In this sense $L_E/L_O \approx 1$ may offer a somewhat more conceptually appealing transitional value to the saturated regime in our flow, rather than the empirically observed $Fr \approx 0.3$ or $Ri_g \approx 0.25$.

Figure 11. Two-dimensional p.d.f. of the length scale ratio $L_E/L_O$ and the mixing coefficient $\varGamma$ constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. Solid lines indicate the proposed scaling of GV19 as well as empirically observed $\varGamma = 0.3$. Vertical dashed line indicates $L_E/L_O = 1$.

An important observation from the above results is the appreciably large region of the channel with a distinct $\varGamma \sim (L_E/L_O)^{1}$ scaling within a quasi-stationary state as proposed by GV19. This is in direct contrast to the recent study of Howland et al. (2020) who instead found that $\varGamma \sim (L_E/L_O)^{2}$ within this regime. As the $L_E/L_O - \varGamma$ framework is essentially an indirect $Fr$-based framework, we consider the scaling arguments of GV19 for the relationship between $Fr$ and $L_E/L_O$.

In the weakly stratified $Fr \gg O(1)$ regime, it is expected that the overturning length scale is well approximated by the inertial energy containing scale $L_T \sim L_E \sim L_{I}$, hence GV19 shows that

(5.12)\begin{equation} L_T/L_O \sim L_E/L_O \sim L_{I}/L_O= \frac{E_K^{3/2} / \epsilon_K}{\epsilon_K^{1/2}/N^{3/2}} = \frac{E_K^{3/2} N^{3/2}}{ \epsilon_K^{3/2}} = Fr^{{-}3/2}. \end{equation}

Conversely in the limit of strong stratification where $Fr \ll O(1)$ and the effects of buoyancy strongly influence flow dynamics, the overturning scale will be expected to scale as the vertical buoyancy scale such that $L_T \sim L_E \sim L_{N}$, and GV19 obtains

(5.13)\begin{equation} L_T/L_O \sim L_E/L_O \sim L_{N}/L_O = \frac{E_K^{1/2}/N}{\epsilon_K^{1/2}/N^{3/2}} = \frac{E_K^{1/2} N^{1/2}}{ \epsilon_K^{1/2}} = Fr^{{-}1/2}. \end{equation}

However, in the intermediate regime the relationship between $Fr$ and $L_E/L_O$ becomes somewhat ambiguous. For instance, in the decaying homogeneous (unsheared) stably stratified DNS study of Mater et al. (2013), no such intermediate regime was observed for any appreciable range of $Fr$, rather the flow was divided into the two regimes described by (5.12) and (5.13) with a single critical crossover point of $L_T/L_O \sim Fr \sim 1$.

We note that the scaling in GV19 makes no assumptions about the presence of mean shear. However, let us consider that in stratified shear flow it has been shown that the overturning length scale is well approximated by the turbulent shear length such that $L_T \sim L_E\sim L_S$ (Venayagamoorthy & Stretch 2010; Mater & Venayagamoorthy 2014), where

(5.14)\begin{equation} L_S = \frac{E_K^{1/2}}{S}. \end{equation}

Here $L_S$ is the shear analogue of $L_N$ and can be directly related back to the mixing length $L_M = \overline {-u'w'}/S$. We note that the $L_E\sim L_I$, $L_E\sim L_N$ and $L_E\sim L_S$ scaling relationships can also be derived in the framework of dominant time scales by considering the vertical displacement of a fluid parcel from its background stratification similar to the analysis presented in § 4, such that

(5.15)\begin{equation} b' \sim w'N^{2} T_*, \end{equation}

where $T_*$ is some time scale pertinent to the evolution of $b'$ due to overturning of the buoyancy field. For instance, if we consider that for the buoyancy dominated regime the pertinent time scale is $T_N$, and by taking the approximation that $w' \sim E_K^{1/2}$, we obtain

(5.16)\begin{equation} L_E = \frac{b'}{N^{2}} \sim \frac{E_K^{1/2}N^{2} T_N}{N^{2}} \sim E_K^{1/2}T_N = \frac{ E_K^{1/2}}{N} = L_N. \end{equation}

It can hence be readily shown that a substitution of $T_L$ or $T_S$ into (5.15) analogously yields $L_E \sim L_I$ and $L_E \sim L_S$, respectively. In this sense and under the assumption that $L_E \sim L_T$, we present physical scaling arguments that provide support to past studies that have shown correlation between measurements of $L_T$ and the three respective energetic length scales $L_I,L_N,L_S$ across varying flow regimes (Mater et al. 2013; Mater & Venayagamoorthy 2014; Ijichi & Hibiya 2018).

To investigate this further, we plot the two-dimensional p.d.f.s of $Fr$ and the ratios $L_E/L_S$, $L_E/L_I$, $L_E/L_N$ and $L_E/L_O$ in figure 12.

Figure 12. Two-dimensional p.d.f.s of the turbulent Froude number $Fr$ and the length scale ratios: (a) the ratio of the Ellison length $L_E$ to turbulent shear length scale $L_S$; (b) the ratio of $L_E$ to the inertial turbulent length scale $L_I$; (c) the ratio of $L_E$ to vertical buoyancy length scale $L_N$; (d) the ratio of $L_E$ to Ozmidov length scale $L_O$. All p.d.f.s constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. Horizontal dashed line for all figures indicates a ratio of unity. Vertical dashed lines indicate $Fr = 0.3$ and $Fr = 1$.

From figure 12(a) it is clear that within the moderately stratified and saturated regimes for $Fr < 1$, the $L_S$ or alternatively the $T_S$ scaling becomes valid and the ratio $L_E/L_S$ becomes a constant of approximately unity. In particular, the excellent correlation within the moderately stratified regime is conceptually consistent with our visual observations in figure 3 of vigorous KHI-driven mixing arising from shear instabilities. This is further supported by the observation that in the $S^{*} - Fr^{-1}$ regime map in figure 8 the flow predominantly lies within the ‘shear dominated’ regime and is reflected in the $T_S \sim (T_L T_N)^{1/2}$ scaling derived in § 5.3. The scatter in the far ‘left flank’ of the saturated regime can again be attributed to the relatively slow process that is acceleration of the mean flow and development of the mean shear profile in our time-varying simulations. For $Fr > 1$ where we observed the slight disagreement between the $\varGamma \sim (L_E/L_O)^{4/3}$ scaling and our results, we observe a small transitional regime where $L_E/L_S$ grows with $Fr$ before again appearing to plateau to a constant of $O(1)$; although this behaviour is not fully developed at our parameter range. To explain this we consider our assumption that for $Fr > 1$ we expect the inertial scaling $L_E \sim L_I$ such that

(5.17)\begin{equation} L_E/L_S \sim L_I/L_S = \frac{E_K^{3/2}/\epsilon_K }{E_K^{1/2}/S} = \frac{S E_K}{\epsilon_K} = S^{*} \end{equation}

and multiplying through by $N/N$ we obtain

(5.18)\begin{equation} S^{*} \frac{N}{N} = \frac{S E_K}{\epsilon_K} \frac{N}{N} = Ri_g^{{-}1/2} Fr^{{-}1}. \end{equation}

Now recalling $Ri_g \sim Fr^{-2}$ in the limit of high $Fr$ as derived in (5.7) we obtain

(5.19)\begin{equation} L_E/L_S \sim S^{*} = Ri_g^{{-}1/2} Fr^{{-}1} \sim Fr^{1}Fr^{{-}1} = \text{const}. \end{equation}

This is in agreement with our observations in figure 8 that $S^{*}$ remains a constant of $O(1)$ for $Fr>1$.

From figure 12(b), however, we observe that within the $Fr > 1$ regime we do not observe the expected inertial scaling of $L_E \sim L_I$ for an appreciable range of $Fr$. Rather, we see continued growth of the ratio $L_E/L_I$ with increasing $Fr$. We can explain this under two considerations. Firstly, for the data presented our highest measurements of $Fr$ are still of $O(1)$. For comparison, in the study of Mater et al. (2013) (their figure 6), the $L_T \sim L_I$ scaling only appears for flow where $Fr = O(10)$. Secondly, as seen in figure 2, by nature of our DNS configuration the data for high $Fr$ inevitably occurs close to the bottom wall. At the lower vertical extent of $z = 0.2$ for which our p.d.f. results are presented, $L_E$, $L_I$, $L_N$ and $L_O$ are all larger than the geometric confinement length scale defined by the distance from the wall $z$ for all simulations. As discussed in the study Taylor et al. (2005), this creates an additional confinement scale that changes and further complicates the relationship between the varying length scales. In this sense we find the absence of a clear $L_E \sim L_I$ scaling unsurprising and we hypothesize that a $L_E \sim L_I$ regime would manifest for our DNS configuration similar to that of Mater et al. (2013) if we were able to access a flow regime where $Fr \gg O(1)$, and simultaneously all the relevant length scales were smaller than the physical confinement scale $z$. As the confinement effects of physical boundaries on the parameterization of mixing fall outside the scope of this study – we leave this for future work.

Within the moderately stratified regime, we can take the $L_E \sim L_S$ scaling derived above to again obtain

(5.20)\begin{equation} L_E/L_I \sim L_S/L_I = 1/S^{*} = Ri_g^{1/2} Fr^{1}, \end{equation}

similarly to the derivation in (5.18). Taking the new scaling of $Ri_g = Fr^{-1}$ derived in (5.10) we obtain

(5.21)\begin{equation} L_E/L_I \sim Ri_g^{1/2} Fr^{1} \sim Fr^{{-}1/2}Fr^{1} = Fr^{1/2}. \end{equation}

Our results show clear support for this scaling, with a clear region of $L_E/L_I \sim Fr^{1/2}$ developing in the intermediate regime for an entire decade of $Fr$. This falls in direct contrast to the study of Mater et al. (2013) where in the absence of mean shear this scaling does not manifest, but rather a $L_T/L_I \sim Fr^{1}$ relationship is observed for $Fr < 1$. We similarly observe an $Fr^{1}$ scaling relationship for our results, but only in the far ‘left flank’ of the figure within the saturated regime, suggesting a buoyancy dominated regime such that $L_E \sim L_N$ and as derived in Mater et al. (2013) we can obtain

(5.22)\begin{equation} L_E/L_I \sim L_N/L_I = \frac{E_K^{1/2}/N }{E_K^{3/2}/\epsilon_K} = \frac{\epsilon_K}{N E_K} = Fr^{1}. \end{equation}

From figure 12(c) we observe direct support for this with a clear trend that in the saturated regime $L_E/L_N$ becomes a constant of order unity. Furthermore, we note that it can clearly be shown that

(5.23)\begin{equation} L_S/L_N = \frac{E_K^{1/2}/S}{E_K^{1/2}/N} = \frac{N}{S} = Ri_g^{1/2}. \end{equation}

Hence, within the saturated regime for $Fr < 0.3$ it is clear that both $L_E/L_S$ and $L_E/L_N$ become constant in agreement with the scaling in (5.13) and our observation that as the flow evolves towards stationarity, $Ri_g$ will trend towards a constant critical value becoming independent of $Fr$. Within the moderately stratified regime, using $L_E \sim L_S$ and $Ri_g \sim Fr^{-1}$ as derived in (5.10), we can obtain

(5.24)\begin{equation} L_E/L_N \sim L_S/L_N = Ri_g^{1/2} \sim Fr^{{-}1/2}. \end{equation}

This intermediate scaling that arises due to the presence of mean shear is again strongly supported by our results, with a clear collapse of the data across a decade of $Fr$.

Finally for the weakly stratified regime, as shown in Mater et al. (2013), we can substitute the expected inertial scaling $L_E\sim L_I$ to obtain

(5.25)\begin{equation} L_E/L_N \sim L_I/L_N = \frac{E_K^{3/2}/\epsilon_K}{E_K^{1/2}/N} = \frac{E_K N}{\epsilon_K} = Fr^{{-}1}. \end{equation}

At our simulation parameter range, we cannot, however, confirm this scaling for our results within this regime for the reasons discussed above.

In light of the results and derivation above we now consider the ratio of $L_E/L_O$ plotted against $Fr$ in figure 12(d). Indeed, it is clear that in the saturated regime we observe the respective $L_E/L_O \sim Fr^{-1/2}$ as proposed by GV19 in (5.13) and shown in the derivation above. As we barely access the weakly stratified regime we again are unable to definitively test the $L_E/L_O \sim Fr^{-3/2}$ scaling at our parameter range. For the moderately stratified regime we now consider that we have explicitly shown $L_E \sim L_S$ and also recall our scaling of $S \sim (\epsilon _K N /E_K)^{1/2}$ in (5.8) to obtain

(5.26)\begin{equation} L_E/L_O \sim L_S/L_O = \frac{E_K^{1/2}/S}{\epsilon_K^{1/2}/ N^{3/2}} \sim \frac{E_K /\epsilon_K^{1/2} N^{1/2}}{\epsilon_K^{1/2}/ N^{3/2}} = \frac{E_K N}{\epsilon_K} = Fr^{{-}1}, \end{equation}

which is the scaling proposed by GV19 but explicitly derived for our shear driven flow and shown to hold only for the intermediate range of $0.3 \lesssim Fr \lesssim 1$ within our flow. Hence we obtain their scaling for $\varGamma$ as follows:

(5.27)\begin{equation} L_E/L_O \sim Fr^{{-}1} \Rightarrow \varGamma \sim Fr^{{-}1} \sim (L_E/L_O)^{1}; \end{equation}

as demonstrated in our results in figure 11.

Furthermore, due to the relationships between $Ri_g$ and $Fr$ derived in (5.7) and (5.10) it can be readily shown that $L_E/L_O$ is directly related to $Ri_g$ across the three regimes. Within the weakly stratified regime we expect

(5.28)\begin{equation} L_E/L_O \sim Fr^{{-}3/2} \sim Ri_g^{3/4}.\end{equation}

Within the moderately stratified regime we expect

(5.29)\begin{equation} L_E/L_O \sim Fr^{{-}1} \sim Ri_g^{1}.\end{equation}

Within the saturated regime, where $Ri_g$ and $Fr$ become decoupled, we expect that $L_E/L_O$ will become independent of $Ri_g$, similar to the results of Rohr et al. (1988) for high $Ri_g$. Figure 13(a) shows the two-dimensional p.d.f. of $Ri_g$ and $L_E/L_O$, with reasonable collapse for the scaling of the three regimes presented above. Due to the scatter in the ‘right flank’ of the curve we illustrate this point more clearly by plotting the $Ri_g$ bin-averaged set of $\langle L_E/L_O \rangle$ against corresponding bins of $\langle Ri_g \rangle$ in figure 13(b). A colour bar depicting $\langle Re_B \rangle$ is again included for reference. From the results it becomes clear that for $Ri_g \lesssim 0.25$ or $L_E/L_O \lesssim 1$ the data collapses well along lines of the proposed scaling derived in (5.28) and (5.29). Meanwhile, in the ‘right flank’ of the curve we observe that $L_E/L_O$ grows independent of $Ri_g$ in agreement with the analysis above, showing separate horizontal trajectories for each simulation.

Figure 13. (a) Two-dimensional p.d.f. of the gradient Richardson number $Ri_g$ and the length scale ratio $L_E/L_O$ constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. (b) The $Ri_g$ bin-averaged $\langle L_E/L_O \rangle$ plotted against bins of corresponding $\langle Ri_g \rangle$ for all data points within $t>1$ and $z>0.2$. Solid lines indicate the proposed scaling in (5.28) and (5.29). Vertical dashed line indicates $Ri_g = 0.25$. Horizontal dashed line indicates $L_E/L_O = 1$.

In this sense, for our quasi-steady shear driven flow it becomes clear that the $Fr, Ri_g$ and $L_E/L_O$ frameworks for the parameterization of the mixing efficiency can all be directly reconciled across the weakly and moderately stratified regimes, with a clear transition to the saturated regime at $Fr \approx 0.3$, $Ri_g \approx 0.25$ or $L_E/L_O \approx 1$, with it being implicit that these three values are interchangeable for our flow. In light of the above analysis we find that the results of Mater et al. (2013) or Howland et al. (2020) do not contradict ours, as their studies inherently have zero mean shear and hence the $L_E \sim L_S$ scaling within the intermediate regime becomes invalid. As such, for non-sheared flows we expect that the buoyancy scaling of $L_E \sim L_N$ to hold across both the intermediate and saturated regimes for $Fr \lesssim 1$ as in Mater et al. (2013), and accordingly may explain the $\varGamma \sim Fr^{-1} \sim (L_E/L_O)^{2}$ scaling observed in Howland et al. (2020) for the moderately stratified regime. Hence in flows where the mean shear is not significant, we expect no appreciable range of $L_E/L_O$ to develop a $\varGamma \sim (L_E/L_O)^{1}$ scaling as in the aforementioned studies. Such findings suggest that the inference of $Fr$ and hence the state of turbulence and mixing through direct field measurements of the overturning length scale within the $Fr = O(1)$ regime may prove a problematic task as it would require additional information as to the state of the flow.

5.5. $Re_B$ framework and transition to diffusive regime

Since the important work of Shih et al. (2005), recent studies have shown that $Re_B$ may not be an optimal parameter in the parameterization of mixing efficiency as it does not truly describe the strength of stratification within the flow in the same sense as $Fr$ or $Ri_g$ (Maffioli et al. 2016; Scotti & White 2016; Garanaik & Venayagamoorthy 2019; Portwood et al. 2019). Rather, it is argued that since $Re_B = (L_O/L_K)^{4/3}$, where $L_K = (\nu ^{3} /\epsilon _K)^{1/4}$ is the well known Kolmogorov scale, its use should be restricted to a measure of the size of the inertial subrange or how ‘energetic’ the flow is.

We explore this by plotting the two-dimensional p.d.f. of $Re_B$ and $\varGamma$ in figure 14(a). From the results we again observe that our flow is divided into three mixing regimes. Again we observe a clear saturated regime where $\varGamma$ trends towards constancy and a distinct regime where $\varGamma \sim Re_B^{-1/2}$, in agreement with the ‘intermediate’ and ‘energetic’ regimes proposed by Shih et al. (2005). Furthermore, we observe a region of the flow corresponding to the weakly stratified regime with a clear $\varGamma \sim Re_B^{-1}$ scaling, in agreement with the Couette flow results and scaling derived in Zhou et al. (2017a). The results suggest the validity of an $Re_B$-based parameterization approach of mixing efficiency for our flow at the parameter range available, albeit with slightly more scatter in the high $Re_B$ ‘right flank’ of the plot with two distinctly separate ‘tails’ emerging in the results. However, as described in MBL16, $Re_B$ can be directly linked to the horizontal Reynolds number $Re_h$ and Froude number $Fr_h$ through the expression $Re_B = Re_h Fr_h^{2}$. By considering $E_K^{1/2}$ as the velocity scale instead of $u_h$, a similar expression can be constructed for turbulent Reynolds number $Re_T$ and $Fr$ such that

(5.30)\begin{equation} Re_B = Re_T Fr^{2}, \end{equation}

where $Re_T = E_K^{2}/(\nu \epsilon _K)$ is the turbulent Reynolds number (Mater & Venayagamoorthy 2014). Hence the parameterization of $\varGamma$ through $Fr$ can be directly related back to $Re_B$. For the weakly stratified regime it is clear that

(5.31)\begin{equation} \varGamma \sim Fr^{{-}2} \Rightarrow \varGamma \sim \left( \frac{Re_B}{Re_T}\right)^{{-}1}. \end{equation}

For the moderately stratified regime we obtain

(5.32)\begin{equation} \varGamma \sim Fr^{{-}1} \Rightarrow \varGamma \sim \left( \frac{Re_B}{Re_T}\right)^{{-}1/2}. \end{equation}

For the saturated regime, provided the flow remains turbulent, $\varGamma$ will become independent of both parameters. As discussed by MBL16, in this sense and under the assumption that $\varGamma$ is fundamentally linked to $Fr$ rather than $Re_B$, an $Re_B$-based parameterization inherently contains an $Re_T$ dependence within itself. Furthermore, taking the estimation that in weakly and moderately stratified flow $E_K^{1/2} \sim u_{\tau }$ and $\epsilon _K \sim u_{\tau }^{3} /\delta$, it can be clearly shown that

(5.33)\begin{equation} Re_T = \frac{E_K^{2}}{\nu \epsilon_K} \sim \frac{u_{\tau}^{4} \delta}{\nu u_{\tau}^{3}} = \frac{u_{\tau \delta}}{\nu} = Re_{\tau}. \end{equation}

Hence, provided the scaling in (5.31)–(5.32) is valid, we expect the flow should show $Re_{\tau }$ sensitivity in the parameterization of mixing efficiency through $Re_B$ within these regimes.

Figure 14. (a) Two-dimensional p.d.f. of the buoyancy Reynolds number $Re_B$ and the mixing coefficient $\varGamma$ constructed out of the instantaneous data of all simulations within the range of $t>1$ and $0.2 \leqslant z \leqslant 0.8$. (b) The $Re_B$ bin-averaged mixing coefficient $\langle \varGamma \rangle$ plotted against bins of corresponding buoyancy Reynolds number $\langle Re_B \rangle$ for all data points within $t>1$ and $z>0.2$. Solid lines indicate scaling lines of $\varGamma \sim Re_B^{-1/2}$ and $\varGamma \sim Re_B^{-1}$ as well as empirically observed $\varGamma = 0.3$. Vertical dashed line indicates $Re_B = 1$.

To investigate this we consider an $Re_B$ bin-averaged data set of $\langle \varGamma \rangle$ plotted against corresponding bins of $\langle Re_B \rangle$ in figure 14(b) in the same manner as for the $Fr$ sorted data in figure 6(b). From the results it is clear that there is no singular transitional value of $Re_B$ from the saturated regime of constant mixing efficiency to the $\varGamma \sim Re_B^{-1/2}$ regime. Rather, two separate evolution paths develop for the $Re_{\tau,0} = 400$ and $Re_{\tau,0}=900$ cases, respectively, showing a clear Reynolds number dependence on the transition. Similarly, no clear singular $Re_B$ value emerges for the transition to the $\varGamma \sim Re_B^{-1}$ regime. In contrast to the $Fr$ averaged results in 6(b), it becomes clear that an $Re_B$-based parameterization of mixing efficiency is inherently dependant on $Fr$, while the reverse is untrue.

In the ‘left flank’ of figure 14(b) we can observe that the transition away from a constant $\varGamma$ regime to the diffusive regime is well approximated by $Re_B \approx 1$ in agreement with the theory of Brethouwer et al. (2007). The exception is the data for case R400L05 (green diamonds), where due to the relatively low level of stability, the flow remains turbulent up to the free surface. In this case the free surface itself rather than buoyancy acts to confine and modify the turbulence properties (Calmet & Magnaudet 2003; Flores et al. 2017), causing deviation from the constant $\varGamma$ regime. The exact mechanics of this is relatively unknown and is an area of study within itself, and are subsequently outside the scope of this study.

6. Concluding remarks

In this study we have investigated temporally evolving stratified open channel flow through DNSs as the flow transitions from a neutral to a stably stratified state, with the emphasis of the study being on the parameterization of mixing across varying energetic regimes within stratified channel flow and the subsequent analysis of the relationship between the relevant non-dimensional mixing diagnostics.

We find that after an initial transient adjustment period of approximately one eddy turnover time unit ($t \approx 1$), the turbulent flow within the channel is distinctly divided into weakly stratified, moderately stratified and saturated mixing regimes separated by transitional values of $Fr =1$ and $Fr =0.3$ across all simulations. Within the three regimes we find that instantaneous measurements of the mixing coefficient $\varGamma$ are predicted well through both the $Fr$ and $L_E/L_O$ parameterization frameworks as outlined in MBL16 and GV19. To our knowledge, ours is the first DNS study to extensively test both the $Fr$ and $L_E/L_O$ parameterization frameworks for stratified wall-bounded flow across a wide range of $Fr$. Considering the strong inherent vertical inhomogeneity within our sheared flow due to the depth varying flux profiles as well as the spatiotemporally evolving mean gradients $S$ and $N$, the remarkable collapse of the results from purely instantaneous measurements of $Fr$ within our flow presents a very strong argument in favour of the case put forward by MBL16 and GV19 for the applicability of an $Fr$-based approach to parameterization of mixing across a variety of stratified flows.

A defining characteristic of our flow is that the majority of the channel evolves into an energetically quasi-stationary state of $Fr = O(1)$. Within this regime we are able to explicitly verify the novel ‘moderately stratified’ scaling of GV19 by showing that $B \sim E_K^{3/2}/L_{N}$, invariant in time and only within the range of $0.3< Fr<1$. By considering our flow within the inertia–shear–buoyancy regime map of Mater & Venayagamoorthy (2014), we find this regime describes a critical state where the inertial, shear and buoyancy forces are all significant in describing the energetic state of the flow. We subsequently provide physically based scaling arguments to show that $T_S \sim (T_N T_L)^{1/2}$ and $Ri_g \sim Fr^{-1}$ within this regime, hence reconciling the concept of a separate intermediate $\varGamma \sim Fr^{-1}$ scaling with the established evidence that in sheared flow and for $Ri_g<0.25$, the turbulent Prandtl number is approximately unity resulting in a linear relationship between the mixing efficiency and $Ri_g$. In contrast we find that for $Fr > 1$ in the weakly stratified regime, buoyancy becomes negligible and a simple balance between shear and inertial forces leads to $Ri_g \sim Fr^{-2}$ in agreement with the arguments presented by MBL16. By considering the mixing length scaling in shear flow of $L_E \sim L_S$ within the intermediate regime, we demonstrate that the remarkable collapse of the data with a distinct $\varGamma \sim (L_E/L_O)^{1}$ scaling for $0.3< Fr<1$ as proposed by GV19 comes as a direct result of the $T_S$ scaling presented above. However, our analysis suggests that in flows devoid of mean shear, the $L_E/L_O - Fr$ and hence $L_E/L_O - \varGamma$ scaling within this regime may differ, implying limited applicability to a wider range of stratified flows under a single framework.

As such, the results suggest that for weakly and moderately stratified quasi-stationary shear flow, the three parameterization schemes become equivalent. As the unification of parameterizing mixing across various fields and applications in the study of stratified turbulence remains a pressing challenge (Caulfield 2020), the results of our study present strong evidence of universal mixing behaviour that appears invariant under differing frameworks in these ubiquitous shear driven flows. Furthermore, as our DNS configuration is an idealization of stratified river flows (Williamson et al. 2015; Kirkpatrick et al. 2019), the present results suggest that under sufficient levels of stratification an appreciable region will inevitably develop where $Fr = O(1)$ and the separate scaling relationships derived in this study for the moderately stratified regime become physically relevant.

For flow that remains turbulent and within the regimes described by the equivalent transitional values $Fr < 0.3$, $Ri_g> 0.25$ or $L_E/L_O >1$, we find that the mixing efficiency saturates to a constant asymptotic value, seemingly in agreement with the strongly stratified scaling of MBL16, however, we note our lowest values of $Fr$ obtained in this study are still significantly higher than the theoretical upper limit of the strongly stratified regime. We further find that in agreement with the theory of Brethouwer et al. (2007), the transition away from a saturated mixing efficiency into the diffusive regime occurs at $Re_B \approx 1$ for our flow, with the caveat that the transition occurs sufficiently far from the free-surface boundary.

Furthermore, the DNS results further suggest that in the saturated regime $Fr$ and $Ri_g$ become decoupled. This suggests that the strongly stratified regime in the context of the stratified turbulence theory of Billant & Chomaz (2001) or Brethouwer et al. (2007) may be accessible within open channel flow even under the assumption that some upper limit pertains to $Ri_g$ to maintain local equilibrium. Whether this result pertains not only to open channel flow but also to wider range of shear flows remains an important open question.