1. Introduction
Stably stratified turbulence may be thought of as a model flow that is potentially useful for understanding aspects of geophysical and engineering processes in many regions of the oceans and atmosphere. In particular, stratified turbulence describes the relatively small-scale dynamics, in length and time, at which turbulence and irreversible mixing occurs. Interpreting and modelling such idealised flows is a primary avenue for the development and calibration of global circulation and basin-scale models, for which mixing plays a leading-order role in global energy budgets (e.g. Ferrari & Wunsch Reference Ferrari and Wunsch2009; Jayne Reference Jayne2009; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018). However, even in idealised stably stratified flow configurations, parameterised modelling of mixing has emerged as a challenge due to the potential for dependence on a wide range of non-independent parameters (Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018; Ivey, Bluteau & Jones Reference Ivey, Bluteau and Jones2018; Caulfield Reference Caulfield2021).
Parallel to the growing recognition of the importance of modelling such small-scale processes in stratified turbulence for broader geophysical applications, stratified turbulence has been increasingly observed to exhibit certain quantitatively similar small-scale dynamics to isotropic turbulence in some parameter regimes, specifically when the buoyancy Reynolds number, $ {{Re_b}}$, is large enough (Gargett et al. Reference Gargett, Hendricks, Sanford, Osborn and Williams1981; Lindborg Reference Lindborg2006; de Bruyn Kops Reference de Bruyn Kops2015). Of course, this does not mean that stratified turbulence is equivalent in all respects, not least because stratification inevitably introduces anisotropy. We consider $ {{Re_b}}$ in more detail in the next section, but in simple terms it quantifies the dynamic range of turbulent length scales neither directly affected by large-scale stratification nor by molecular viscosity. Indeed, it is widely acknowledged that the details of the large-scale, or outer, mean-scale, dynamics do not strongly affect the small-scale dynamics provided there is sufficient scale separation between such large and small scales. While the small-scale dynamics is certainly not independent of large scales (Corrsin Reference Corrsin1958; Durbin & Speziale Reference Durbin and Speziale1991), the assumption that sufficient scale separation induces statistically independent small scales has proven to be a valuable tool in the modelling of anisotropic turbulent flows because it allows for the application of theoretical models based on statistical symmetries, e.g. isotropy and homogeneity, to dynamic modelling.
To take one example, Kolmogorov–Obukhov scaling depends on the assumptions of local isotropy and homogeneity (Kolmogorov Reference Kolmogorov1941; Oboukhov Reference Oboukhov1941) yet it has been observed to be a good approximation for anisotropic flows provided the scale separation is sufficiently large between anisotropic turbulence scales and the dissipative viscous–diffusive scales (Champagne, Harris & Corrsin Reference Champagne, Harris and Corrsin1970; Gargett, Osborn & Nasmyth Reference Gargett, Osborn and Nasmyth1984; Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994; Shen & Warhaft Reference Shen and Warhaft2002). However, not all stably stratified flows exhibit sufficient scale separation for the foregoing assumption to be made (Jackson & Rehmann Reference Jackson and Rehmann2014; de Bruyn Kops Reference de Bruyn Kops2015).
This phenomenology is useful because it motivates the adaptation of existing models for isotropic turbulence to stably stratified turbulence. For instance, by analysis of Obukhov–Corrsin and Kolmogorov scalings, Beguier, Dekeyser & Launder (Reference Beguier, Dekeyser and Launder1978) observed that the passive scalar time scale $\tau _\phi$ and turbulence time scale $\tau$ couple at sufficiently high Reynolds number such that they can be related by universal constants, i.e.
where $E_\phi$ is the turbulent scalar variance, $E_k$ is the turbulent kinetic energy, $\epsilon _\phi$ and $\epsilon$ are their respective dissipation rates and $\beta,C$ are the Obukhov–Corrsin and Kolmogorov constants, respectively. The assumption that the time-scale ratio is a constant is commonly applied to model scalar dissipation in Reynolds-averaged Navier–Stokes models (e.g. Newman, Launder & Lumley Reference Newman, Launder and Lumley1981; Ristorcelli Reference Ristorcelli2006).
In modelling mixing in stratified flows, the applicability of such a universal relation is interesting due to its implications for modelling the irreversible diapycnal flux, $j_b$. In such models, equilibrium assumptions are typically prescribed such that the irreversible diapycnal flux is equal to a particular definition of the available potential energy dissipation rate (Osborn & Cox Reference Osborn and Cox1972; Peltier & Caulfield Reference Peltier and Caulfield2003; Caulfield Reference Caulfield2021) such that coupling of the turbulent and scalar dynamics, as suggested by Beguier et al. (Reference Beguier, Dekeyser and Launder1978), can be an insightful relation for modelling the scalar dynamics.
The importance of accurately modelling such irreversible scalar fluxes in the ocean has motivated significant research on various idealised flows using an array of modelling frameworks. Perhaps most notable, estimating the irreversible flux from the scalar gradient via a turbulent diffusivity, $\kappa _b$, is the de facto standard approach (e.g. Ivey & Imberger Reference Ivey and Imberger1991; Barry et al. Reference Barry, Ivey, Winters and Imberger2001; Salehipour & Peltier Reference Salehipour and Peltier2015; Maffioli & Davidson Reference Maffioli and Davidson2016), with
where $N$ is an appropriately defined large-scale buoyancy frequency. A common parametrisation of $\kappa _b$ has been suggested by Osborn (Reference Osborn1980), which asserts that $\kappa _b$ be related to the irreversible rate of dissipation of kinetic energy via a turbulent flux coefficient $\varGamma _b$ such that
By analysis of the energy equations for a simple stratified flow model at stationary conditions, Osborn (Reference Osborn1980) suggested that $\varGamma _b$ is a constant $\leq 0.2$. More recently, the prescription of the parameter dependence for $\varGamma _b$ has emerged as a necessary improvement to turbulent diffusivity models due to increasing evidence that a constant coefficient for $\varGamma _b$ appears not to be at all suitable. Parametrisations of $\varGamma _b$ in this framework have been determined to be difficult due to a potential dependence on multiple parameters, which themselves may well be correlated (Caulfield Reference Caulfield2021), which difficulties lead to often contradictory subclosures (Monismith, Koseff & White Reference Monismith, Koseff and White2018). While alternative models exist, such as mixing length models (Odier et al. Reference Odier, Chen, Rivera and Ecke2009; Ivey et al. Reference Ivey, Bluteau and Jones2018) and flux transport models, they have not as yet gained widespread usage in large-scale simulations.
Therefore, as the first principal objective of this research, we consider formally adapting the Beguier et al. (Reference Beguier, Dekeyser and Launder1978) relationship (1.1) to stably stratified turbulence in order to model irreversible diapycnal mixing dynamically. To evaluate this hypothesis, we consider simulations of stationary homogeneous stratified and sheared turbulence (S-HSST). S-HSST is an idealised model flow configuration that enables easy adjustment of the range of length scales (the dynamic range) available for a locally isotropic subrange to form without changing other dimensionless flow parameters such as the characteristic Froude or Richardson numbers, as defined precisely in the next section. Such a fundamental flow configuration is well suited to the evaluation of the effects of locally isotropic scaling theories (Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000; Chung & Matheou Reference Chung and Matheou2012). To test the hypothesis that the Beguier et al. (Reference Beguier, Dekeyser and Launder1978) relation can be applied to stably stratified turbulence, we perform numerical experiments of S-HSST in parameter space extremes previously inaccessible to computation. Therefore, a second principal objective of this work is to report on the characteristics of S-HSST at very large Reynolds numbers and time scales, quantitatively described in subsequent sections.
To achieve these objectives, the rest of the paper is organised as follows. In § 2, we discuss a length-scale-based framework for stratified and sheared turbulence and use it to derive a mixing model hypothesis based on the arguments of Beguier et al. (Reference Beguier, Dekeyser and Launder1978). In § 3, we then describe numerical simulation experiments of homogeneous stratified and sheared turbulence (HSST). In a Reynolds number parameter space, we then present the ensemble-averaged dynamics in § 4, evaluate the mixing model in § 5 and finally draw our conclusions in § 6.
2. Theoretical background
2.1. Parametric framework
HSST is assumed to satisfy the Navier–Stokes equations subject to the non-hydrostatic Boussinesq approximation. The dimensional equations for the fluctuations relative to the planar means are
where $\boldsymbol {u} = (u_x, u_y, u_z)$ is the velocity vector in the coordinate system $(x,y,z)$ and functional dependencies have been dropped for convenience, and $\rho$ is the fluctuating density field so that the total density is
with $\bar {\rho }(z) = z \, {{\rm d}} \bar {\rho }/ {{\rm d}}z$ being the ambient density that is constant in time and has a uniform gradient antiparallel to gravitational acceleration $g_z$. Similarly, $\bar {u}_z= z \, {{\rm d}} \bar {u}_x/ {{\rm d}}z=z S$ is the mean velocity and $S\equiv {{\rm d}} \bar {u}_x/ {{\rm d}}z$ is the mean shear. The pressure decomposition is analogous to that of density and velocity so that $p$ is the fluctuating mechanical pressure relative to the temporally constant planar mean. Internal energy and scalar transport affecting density have been combined into a single evolution equation for $\rho$ with $D$ the molecular diffusivity, $\hat {\boldsymbol {x}}$ and $\hat {\boldsymbol {z}}$ are the unit vectors in the $x$- and $z$-directions respectively and $\nu$ is the kinematic viscosity.
2.1.1. Parametrisation
Dimensional analysis suggests that (2.1) can be described in terms of at least four non-dimensional parameters. As we wish to interpret the dynamics in terms of the dynamic range available for various aspects of the flow, we consider parametrisation in terms of the ratios of length scales.
Here, we use the convention of defining an approximate outer scale of turbulence, the large-eddy length scale, and the (viscous) Kolmogorov length scale as
where $E_k \equiv \langle \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u} \rangle /2$, $\epsilon \equiv \nu \langle \boldsymbol {\nabla } \boldsymbol {u} : \boldsymbol {\nabla } \boldsymbol {u} \rangle$, the operator ‘$:$’ denotes the double inner product and the notation $\langle \boldsymbol {\cdot } \rangle$ indicates an ensemble average.
The (buoyancy) Ozmidov length scale
defines the lower limit of length scales significantly affected by buoyancy forces, where $N$ is the (background) buoyancy frequency, defined here as
With an explicit shear scale applied in HSST, the Corrsin length scale (Corrsin Reference Corrsin1958) is thought to characterise the lower limit of scales affected by shear, and is defined as
Scales smaller than $L_C$ have been suggested to define a locally isotropic regime of length scales in the absence of additional smaller-scale affects of the mean flow (Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994).
Given the scales defined in (2.3a,b), (2.4) and (2.6), we choose to describe turbulence by the following three parameters:
in a fluid with a fixed Prandtl number $Pr\equiv \nu /D$ as the required fourth parameter (see Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) for more details on parameterisations). We choose the shear Reynolds number, $Re_s$, since it describes the range of length scales associated with isotropy for stationary flows with Richardson number, $Ri \leq 1$ (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014). Alternative Reynolds numbers can also be useful to consider for comparison with other flow regimes, where the turbulent Reynolds number and buoyancy Reynolds numbers are defined, respectively as
In flows driven by mean shear, the energetic stationarity of the flow is thought to be characterised by $Ri$ (Jacobitz, Sarkar & Van Atta Reference Jacobitz, Sarkar and Van Atta1997), measuring the scale separation between the driving shear scales and stabilising buoyancy scales. It is thought that turbulent flow configurations persist where $Ri<1$, where we distinguish between sustained turbulent flows and results derived from or interpreted in terms of linear hydrodynamic instability theory (see Zhou et al. (Reference Zhou, Taylor, Caulfield and Linden2017) for more discussion). Indeed, stratified turbulent simulations with prescribed shear scales have revealed that the characteristic gradient Richardson number associated with stationarity has been observed to take a Reynolds number insensitive value of
as observed in stratified homogeneous shear flow (Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000; Holt, Koseff & Ferziger Reference Holt, Koseff and Ferziger1992; Portwood, de Bruyn Kops & Caulfield Reference Portwood, de Bruyn Kops and Caulfield2019).
Furthermore, in shear-dominated flows, energy production due to shear naturally induces a strong coupling between shear and outer length scales such that the shear parameter $ {{S_\ast }}$, here defined in terms of $E_k$ (as opposed to in terms of $q\equiv \langle \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u} \rangle$) as
tends to a constant between 5 and 6 at sufficiently high Reynolds numbers (Jacobitz et al. Reference Jacobitz, Sarkar and Van Atta1997; Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000)
However, the foregoing turbulent parameters are likely insufficient to characterise active scalar dynamics, a specific point we wish to investigate in this paper. The smallest scales of the scalar are thought to be characterised by the Batchelor length scale $L_B$,
which coincides with $L_K$, at unity Prandtl number as considered in this work.
Obukhov–Corrsin similarity, which is either applicable at high Froude number or at scales smaller than those substantially affected by buoyancy, i.e. for scales smaller than $L_O$, implies an outer-scale scalar length scale
where $E_p$ is the available potential energy and $\chi$ is its irreversible dissipation rate, defined in the linearly stratified limit to be
respectively, demonstrating that $\chi$ in this context may also be interpreted as the scaled destruction rate of buoyancy variance (Caulfield Reference Caulfield2021). Therefore, the range of anisotropic scales of the scalar are characterised by
Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019) observed this parameter to approach unity in the high $Re_b$ limit for flows where $Pr=1$, and Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) suggest $N_\ast$ more descriptively parameterises turbulent mixing.
2.2. A mixing model constructed from analysis of subrange scaling
Local isotropy is a state wherein statistical symmetries of multi-point statistics, such as homogeneity, isotropy and stationarity, are present in a spatio-temporally localised region (Monin & Yaglom Reference Monin and Yaglom1975). In the presence of anisotropic integral scales, the conditions of local isotropy are defined, in the spatial sense, by a quasi-equilibrium range of scales that are sufficiently smaller than integral scales (Kolmogorov Reference Kolmogorov1941).
This is a prerequisite of classical inertial subrange scaling arguments (Kolmogorov Reference Kolmogorov1941; Oboukhov Reference Oboukhov1941; Corrsin Reference Corrsin1951). These state that scales within the quasi-equilibrium regime exist wherein turbulence is independent of the effects of viscosity. In shear flows, conditions for the application of local isotropy are thought to be valid at scales smaller than $L_C$ (Uberoi Reference Uberoi1957; Corrsin Reference Corrsin1958). Kolmogorov scaling (Kolmogorov Reference Kolmogorov1941, Reference Kolmogorov1962) in the presence of anisotropic outer scales has been revealed to coincide at scales consistent with local isotropy (Champagne et al. Reference Champagne, Harris and Corrsin1970; Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994). The equivalent local isotropy justification of the Kolmogorov subrange may be applied to Obukhov–Corrsin passive scalar scaling (Oboukhov Reference Oboukhov1949; Corrsin Reference Corrsin1951) such that it is fit to describe an active scalar below scales affected by stratification (Monin & Yaglom Reference Monin and Yaglom1975, p. 391). Therefore, the spectral density of potential energy is expected to be determined by the energy dissipation rates and a wavenumber, i.e.
where $\boldsymbol {k}$ is the wavenumber vector, $\beta$ is the Obukhov–Corrsin constant, the functions denoted by $f^{OC}_v$, $f^{OC}_L$ are the viscous and outer-scale correction functions that are necessary without applying assumptions of local isotropy. Furthermore, according to Oboukhov (Reference Oboukhov1941), the kinetic energy spectra
where $C$ is the Kolmogorov constant, and $f^{K}_L$ and $f^{K}_v$ are the correction functions. Following the analysis of Beguier et al. (Reference Beguier, Dekeyser and Launder1978), by integration of (2.15) and (2.16) when $L_C \gg |\boldsymbol {k}|^{-1} \gg L_K, L_B$, it is possible to obtain
Therefore, the explicit relation for modelling the irreversible flux from (2.17) is
with the limiting behaviour, at high Reynolds number, that $\varPi \rightarrow \beta /C$. A critical simplifying assumption in the development of (2.18) is that $ {{Fr}} \geq O(1)$, although the imposition of more complex scaling relations, which are thought to develop when $ {{Fr}} \ll O(1)$, may alternatively be imposed. Indeed, more generally,
We stress that the assumed model spectra (2.15) and (2.16) are approximations which are difficult to observe precisely in carefully controlled experiments or simulations. Whereas these two-point scaling models are observed in flows with anisotropic large scales, they are not claimed to be generally applicable for flows in arbitrary parameter regimes. However, the sensitivity of $\varPi$ with respect to measured deviations from the model spectra is unclear a priori, a point we discuss in § 4.
The relation (2.18) has profound implications for the calibration of a number of turbulent stratified mixing models. Generally, in gradient diffusion models,
while gradient diffusion subject to the model proposed by Osborn (Reference Osborn1980) leads to
Combining the turbulent viscosity model of Crawford (Reference Crawford1982) with the turbulent diffusivity model of Osborn (Reference Osborn1980), the turbulent Prandtl number may then be expressed as
Alternatively, the relevance of the parameter $\varPi$ to the turbulent Prandtl number model of Venayagamoorthy & Stretch (Reference Venayagamoorthy and Stretch2010) implies that their key mixing parameter $\gamma$ may be expressed as
where $\gamma$ too has been empirically observed to assume relatively insensitive values at moderate Reynolds numbers in HSST (Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2006). Finally, if the objective is to construct a mixing length scale for applications in Prandtl mixing length models, such as the framework of Odier et al. (Reference Odier, Chen, Rivera and Ecke2009), where $\kappa _b = l_m^{2} S$, (2.18) may be equivalently expressed as
3. Direct numerical simulations
In order to test the underlying hypotheses described in § 2.2 leading to the irreversible mixing model embodied in (2.18), we must consider a stratified model flow wherein the behaviour of $\varPi$ may be evaluated as a function of the dynamic range associated with isotropy. The consideration of S-HSST is justified by the emergence of a Reynolds number as the only independent non-dimensional descriptive flow parameter due to two distinct phenomenologies: (i) the coupling of the outer turbulence scales to that of the shear production mechanism such that $ {{S_\ast }} \approx 5$ and (ii) the gradient Richardson number assuming a balanced apparently critical value that maintains a stationary balance of energy production and dissipation mechanisms. Both phenomenologies were discovered in direct numerical simulations for this flow configuration by the investigations of Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000, Reference Shih, Koseff, Ivey and Ferziger2005) and subsequently verified more recently in a broader Reynolds number space in simulations where the Richardson number is controlled to maintain constant turbulent kinetic energy (Portwood et al. Reference Portwood, de Bruyn Kops and Caulfield2019). We therefore utilise the flows presented in Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019) to verify the modelling framework presented in § 2.2 and relevant associated consequences highlighted later in this manuscript.
3.1. Numerical method
The homogeneous stratified shear system defined by (2.1a) and (2.1b) is solved in a triply periodic domain with a Fourier pseudo-spectral method. Timestepping is performed with a fractional step method where nonlinear terms are advanced with a third-order Adams–Bashforth scheme. Aliasing is treated by a sharp-spectral filter at $15/16 k_{max}$, which proved to be sufficient in a posteriori analysis. Temporal integration of the inhomogeneous shear term is handled by the integrating factor approach of Sekimoto, Dong & Jiménez (Reference Sekimoto, Dong and Jiménez2016). The rest of the solver is derived from the method used in Hebert & de Bruyn Kops (Reference Hebert and de Bruyn Kops2006), with further details discussed therein.
The suite of simulations are designed to resolve sufficiently both the smallest and largest spatio-temporal turbulent flow scales. The streamwise extent, $L_x$, is chosen to be 40 times larger than the large-eddy length scale. Because of the anisotropy of integral length scales that was observed in a limited parameter space study, the cross-stream domain length $L_y$ and the vertical domain length $L_z$ are chosen such that $L_x/L_y=2$ and $L_x/L_z=4$. A visualisation of the computational domain is featured in figure 1.
The smallest spatial scales are resolved by enforcing $k_{max}L_{K} \approx 2$, where $k_\text {max}$ indicates the maximum wavenumber supported by the domain discretisation when fully resolved.
Finally, we have adopted a mass–spring–damper-like system (as similarly adopted in Overholt & Pope Reference Overholt and Pope1998; Rao & de Bruyn Kops Reference Rao and de Bruyn Kops2011), to tune the Richardson number to its stationary value by setting a constant kinetic energy target $E_t$. We note that this is similar to the constant power criteria as used in stationary HSST by Chung & Matheou (Reference Chung and Matheou2012) but here we include damping. The system is controlled by time-dependent variation of $Ri$ being required to satisfy the equation
where the prime notation denotes a temporal derivative, $\tilde {E}_k(t)\equiv E_k/E_t$ is the normalised turbulent kinetic energy, $\omega$ is the characteristic frequency of oscillation, $\alpha$ is a dimensionless damping factor and $c_0$ is a dimensionless parameter. This control system has been derived by assuming that the kinetic energy follows a second-order linear system (e.g. Rao & de Bruyn Kops Reference Rao and de Bruyn Kops2011), and then by applying a first-order approximation to the temporal evolution of kinetic energy (cf. Jacobitz et al. Reference Jacobitz, Sarkar and Van Atta1997)
The choice for the parameter $c_0\approx -1$ is supported by Jacobitz et al. (Reference Jacobitz, Sarkar and Van Atta1997), the characteristic frequency $\omega$ is determined by the large-eddy time scale which itself is proportional to $1/S$ (Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000) and finally a damping coefficient $\alpha =1.5$ was found to work well.
3.2. Summary of experiments
The stationarity constraint and the emergent empirical observation that $S_\ast \approx 5$ lead to the natural consideration of sets of simulations which can test the effects of variation of the Reynolds number, essentially independently of the other parameters. The simulations presented here are largely equivalent to those reported in Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019), although some cases have been run for longer times to ensure convergence of statistics. The various parameters of the simulations are summarised in table 1, and some flow visualisations of the density field are shown in figure 1. Two key aspects are apparent, as noted in the caption. First, in the horizontal plane visualisations, as the Reynolds number increases, the magnitude of the density fluctuations tends to decrease, while, unsurprisingly the dynamic range of scales increases, with enhanced smaller-scale structures. Second, in the vertical plane visualisation, inclined large-scale structures are apparent, consistently with the results of Chung & Matheou (Reference Chung and Matheou2012) and Jacobitz & Moreau (Reference Jacobitz and Moreau2016).
Time series are shown in figure 2 to illustrate convergence. Richardson numbers tend to fluctuate within approximately 10 % of their means, at a confidence interval at least as precise as critical Richardson numbers reported by Jacobitz et al. (Reference Jacobitz, Sarkar and Van Atta1997) and Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000). The fluctuations in the Richardson number occur at large integral time scales such that effects of $Ri'(t)$ are small. Energetic fluctuations are more substantial than higher-dimensional forcing schemes (i.e. Overholt & Pope Reference Overholt and Pope1998; Rao & de Bruyn Kops Reference Rao and de Bruyn Kops2011) that dynamically control the flow with more granularity by accessing a broader range of turbulent length scales compared with a single mean-scale parameter, as is done here. Nonetheless, energy remains approximately stationary over the large time scales of the simulations.
Recalling that the critical Richardson number is an emergent quantity from the control system defined by (3.1), and is not determined a priori, an increase of the critical Richardson number is not observed in these simulations as suggested by Holt et al. (Reference Holt, Koseff and Ferziger1992), Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000), although the emergent critical values of the Richardson numbers are broadly consistent with those reported by Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000). Specifically, we observe the critical gradient Richardson number which induces statistical stationarity in these flows to be approximately 0.16. Similarly, we observe $ {{Fr}}\approx 0.5$ in all simulations, which implies $ {{S_\ast }} \approx 5$, as observed in other studies (Jacobitz et al. Reference Jacobitz, Sarkar and Van Atta1997; Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000).
Curiously, the lowest Reynolds number to maintain stationarity robustly according to our stringent criteria corresponds approximately to case R1. Case R1 corresponds to $ {{Re_b}} \approx 36$, which is consistent with critical values observed and estimated for sustained three-dimensional turbulence previously (Gibson Reference Gibson1980; Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016).
For further comparison with the flows considered by Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000) and Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005), we plot their relevant cases against the solutions presented for this research in a $Fr-Re_b$ parameter space, as shown in figure 3. We remark that the configurations reported by Shih and co-authors in those publications feature transitional and, crucially, variable $Fr$ when $Re_b$ is $O(100)$. As is apparent in the figure, $Fr$ and $Re_b$ appear to be coupled, with an approximate scaling relationship $Fr \propto Re_b^{1/2}$, as shown with a dashed line. Therefore, the parameterisation of mixing by $Re_b$ cannot be made independent of $Fr$ with the data presented in Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000, Reference Shih, Koseff, Ivey and Ferziger2005). We note that this explanation does not rely on any argument that the simulations are inadequately resolved, either at small or large scales. Our perspective is qualitatively different from the perspective presented in Kunze et al. (Reference Kunze, MacKay, McPhee-Shaw, Morrice, Girton and Terker2012) and Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018), which suggests the $Re_b$-sensitive regime with $Re_b>100$ is due to inadequate small-scale spatial resolution in the simulations. However, it is possible that the solutions analysed in Shih et al. suffer from inadequate large-scale resolution in the domain as suggested by other reports (see Kunze Reference Kunze2011; Kunze et al. Reference Kunze, MacKay, McPhee-Shaw, Morrice, Girton and Terker2012; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018).
Therefore, disentangling the behaviour of their simulations as $Re_b$ varies from the effects of variation in $Fr$ is challenging, particularly when $Re_b$ (or $Re_s$) is large. It is important to remember that it was precisely for $Re_b \gtrsim 100$ that Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) argued that the turbulent flux coefficient $\varGamma \propto Re_b^{-1/2}$, whereas it is entirely plausible that variable $Fr$ is playing a significant role. This figure highlights the novelty of our simulations in this particular flow configuration, which allow the isolated study of the effects of variations in Reynolds number from the influence of other non-dimensional turbulent parameters, and in particular at constant $Fr$.
4. Emergent phenomena
4.1. Energetics
Just as non-dimensional parameters are emergent quantities in our simulations, so too are the partitionings of potential energy and the various components of kinetic energy. Such partitionings are significant, not least because the ratio of potential energy to kinetic energy,
is a critical component of mixing models wherein the Reynolds number, or $ {{Re_b}}$, dependence is often omitted and the mixing is assumed to be a function of $ {{Ri}}$ (Osborn & Cox Reference Osborn and Cox1972; Schumann & Gerz Reference Schumann and Gerz1995; Pouquet et al. Reference Pouquet, Rosenberg, Marino and Herbert2018). Furthermore, in ‘strongly’ stratified turbulent flow, Billant & Chomaz (Reference Billant and Chomaz2001) suggests that there should be approximate equipartition between potential and kinetic energy, i.e. $R_{PK}\approx 1$, an assumption also used by Lindborg (Reference Lindborg2006). It is important to remember that these models do not account for Reynolds number effects. Therefore, parameterisation of such energetic partitioning is an important a priori validation necessary for model evaluation.
Recalling that the kinetic energy is the same for all cases, by construction from our definition of stationarity, we observe a significant decrease in potential energy with increasing dynamic range, as shown in figure 4(b). In comparison, the results of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), which span up to $ {{Re_b}} \sim O(10)$ at much lower Froude number in homogeneous stratified turbulence, report $R_{PK} \approx 0.15$ at $ {{Re_b}} = {{Ri}} {{Re_s}} \approx 16$, noting an asymptotic trend to this value from $R_{PK}\approx 0.05$ at $ {{Re_b}}\approx 0.1$. Although we observe similar values at $ {{Re_s}}\approx 5$ (corresponding to $ {{Re_b}} \approx 30$), higher values of $ {{Re_b}}$ reveal a subsequent decline in $R_{PK}$, as suggested by Remmler & Hickel (Reference Remmler and Hickel2012). For Reynolds numbers larger than $ {{Re_s}} \approx 40$ ($ {{Re_b}} \approx 300$), $R_{PK}$ assumes an apparently asymptotic value of $0.08$ in these simulations. This relative decrease of potential energy as a function of $ {{Re_s}}$ can be observed in the domain slice visualisations of $\rho$ featured in figure 1, as there is a clear reduction in the variance of the density fluctuations. In terms of the important question concerning the modelling of mixing, and in particular parameterising the turbulent flux coefficient $\varGamma$, (2.21) shows that, if $R_{PK}$ approaches an asymptotic value for large $Re_s$, $\varGamma$ can exhibit parameter dependence only if the key parameter $\varPi$ (defined in (2.17)) does.
Furthermore, a distinguishing characteristic of sheared stratified turbulence, relative to other stratified flows, is the absence of axisymmetry normal to gravity with the result that anisotropy must be parameterised in each dimension. Inertially relevant anisotropy can be characterised by the variance of individual velocity vector components as they contribute to the mean turbulent kinetic energy. These velocity variances resulting from an anisotropic dynamics are shown in figure 4(a). We observe turbulent kinetic energy is dominated by streamwise velocity fluctuations, with the smallest contributions coming from vertical velocity fluctuations. Notably distinct from axisymmetric flows, the horizontal components of velocity variance represent approximately 80 per cent of the total contributions to kinetic energy. The cross-stream velocity variance accounts for only 30 per cent of the total energy. The vertical variance, subject to exchanges to potential energy, is typically suggested to vary as a function of Froude number (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). Here, the vertical variance accounts for approximately 20 per cent of the total energy and we observe that it remains largely constant as a function of $ {{Re_b}}$.
Perhaps surprisingly, anisotropy of velocity fluctuations increases with increasing dynamic range. That is, the streamwise velocity component becomes increasingly dominant with increasing Reynolds number, apparently at the expense of the cross-stream velocity variance. One possible interpretation is that reducing viscous effects at inertial scales as the dynamic range increases allows the flow to evolve towards an inertial equilibrium state.
4.2. Dynamics
It might be assumed that the significant energetic transitions at sufficiently high Reynolds number would be explained by accompanying clear transitions in the dynamics. Here, we analyse the stationary behaviour of the dynamics relevant to the kinetic and potential energy balances. The ensemble-averaged energy equations, as derived from (2.1), obey
where the turbulent production from mean shear $P=-\langle u_x u_z \, {{\rm d}}\bar {u}_x/{{{\rm d}}z}\rangle$ and the ‘buoyancy’ flux $B=\langle ({g}/{\rho _0}) {u_z} \rho \rangle$ have anisotropic effects on the evolution of the streamwise and vertical components of the kinetic energy $\langle u_x^{2}\rangle$ and $\langle u_z^{2}\rangle$, respectively.
Note that the left-hand sides of (4.3), (4.2) reduce to approximately zero due to the stationarity constraint. Therefore, the dynamics should, at least in principle, be describable by the turbulent flux coefficient $\varGamma$ defined in terms of $\chi$ in (2.21), and a flux Richardson number, defined as
since, if the flow actually enters an emergent stationary state the flux $B \simeq \chi$, and remembering the underlying original definition (1.3) $\varGamma _b \simeq \varGamma$. Furthermore, in this circumstance, the turbulent viscosity $\nu _b$ can be naturally defined as
and so we can obtain an expression consistent with (2.22). Henceforth, we continue to describe the evolution of the energetics by $\varGamma$ defined as in (2.21) in recognition of the long-standing application of different definitions of the turbulent flux coefficient to estimating turbulent diffusivity and turbulent viscosity (Osborn & Cox Reference Osborn and Cox1972; Crawford Reference Crawford1982). Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018) discuss in detail the subtleties and potential for confusion concerning different definitions for turbulent flux coefficients, and so it must always be remembered that here we concentrate on using the definition in (2.21), i.e. $\varGamma \equiv \chi /\epsilon$.
The ratio of terms in the energetic balance is shown in figure 5. We observe a slight decline of $\varGamma$ from $0.2$ to approximately $0.175$, as noted in Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019), until $ {{Re_s}} \approx 40$. However, $\varGamma$ remains approximately constant as $Re_s$ increases further, and in particular we do not observe the proposed $ {{Re_b}}^{-1/2}$ scaling reported in a similar flow by Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) for their energetic regime of $ {{Re_b}} > 100$. Remembering the data presented in figure 3, this is perhaps unsurprising due to the apparent correlated variation of $Fr$ and $Re_b$ in those simulations. Due to the condition of stationarity being imposed on the kinetic energy in our simulations, accompanying induced stationarity of the potential energy is also guaranteed, although still plotted for reference in figure 5, demonstrating that $\chi \approx B$. It is also apparent that $ {{Rf}} \equiv B/P$ is less than $\varGamma \equiv \chi /\epsilon$, as expected.
4.3. Length scales
Since we have observed $ {{Fr}}$ and $ {{Ri}}$ to be essentially independent of the Reynolds number, the relatively decreasing potential energy as $ {{Re_s}}$ increases for $Re_s \lessapprox 50$ implies that scales characteristic of the scalar should also exhibit a similar asymptotic trend. The scalar outer scale associated with Obukhov–Corrsin similarity, $L_{OC}$ is defined in (2.12). In figure 6(a), we plot the ratio of this outer scale to the Ozmidov length scale $L_O$, as defined in (2.4). As $Re_s$ increases, this ratio decreases for $Re_s \lessapprox 40.$ For larger $Re_s$, there is some evidence that the ratio remains constant near unity, suggesting of course that $L_{OC}$ anomalously scales with $L_O$. This implies that the scaling $\chi \approx E_p N$ is valid in this high Reynolds number regime and that the outer scales of the scalar adjust to the regime associated with local isotropy, which is indeed a condition of a Obukhov–Corrsin subrange in anisotropic flows.
Length scales associated with the mixing length models suggested by Odier et al. (Reference Odier, Chen, Rivera and Ecke2009) can also be shown to remain constant in our large $Re_s$ regime, and indeed even for smaller $Re_s$. The scalar mixing length, $L_\rho$, and the momentum mixing length $L_m$, both defined by Odier et al. (Reference Odier, Chen, Rivera and Ecke2009) as
are coupled via the turbulent Prandtl number $ {{Pr_t}}$, which here is very close to one, since
Furthermore, they coincide with the Corrsin length scale by the stationarity constraint from (4.3). The expected scalings $L_\rho \approx L_C$ and $L_\rho \approx L_m$ are both clearly apparent in figure 6(a) across a wide range of Reynolds numbers.
The decrease of $L_{OC}$ implies that the total dynamic range associated with the scalar anomalously scales with the shear Reynolds number $ {{Re_s}}$. We define the total range of dynamic scales for turbulence and the scalar as
respectively, and show them as a function of Reynolds number in figure 6(b). With respect to the turbulence scales, after an initial small $ {{Re_s}}$ transient, the range of total dynamic scales appears to increase with $ {{Re_s}}^{4/3}$, as predicted by our definitions in § 2.1.1. For the scalar, on the other hand, there is clearly a flatter-than-predicted scaling regime at high Reynolds number which is characteristic of the transient outer scales observed in figure 6(a). Furthermore, the range of dynamic scales associated with the scalar is significantly smaller than that of the turbulence, suggesting that, if a dynamic-range threshold exists, it will not occur simultaneously for the scalar and the turbulence.
5. Implications of inertial subrange scaling
Before turning to model verification, an intermediate phenomenological test of the universality hypothesis outlined in § 3.1 may be performed by analysis of the one-dimensional correlations. The one-dimensional spectrum of the potential energy is expected to be determined by the energetic dissipation rates and a wavenumber, i.e. for the cross-stream spectra
In this expression, $\beta _1$ is the one-dimensional Obukhov–Corrsin constant. Furthermore, the functions denoted by $f^{OC}_v$, $f^{OC}_L$ are the viscous and outer-scale correction functions which are necessary if no implicit assumptions of local isotropy are made. The notation $L^{O}_{ij}$ indicates the outer scale of locally isotropic turbulence associated with the velocity component $i$ with respect to the direction $j$. This is ostensibly the Corrsin length scale but is inevitably anisotropic, as discussed in Kaimal (Reference Kaimal1973) for the specific case of the stably stratified boundary layer.
Equivalently, according to Kolmogorov (Reference Kolmogorov1941) for an isotropic inertial subregime, the one-dimensional cross-stream longitudinal and streamwise transverse energy spectra are
where $C_1$ and $C_1'$ are universal constants, and $f^{K}_L$ and $f^{K}_v$ are the correction functions. The ratio of (5.1) and (5.2) yields
From empirical studies of turbulence, even in the absence of the effects of stratification and shear, we expect the correction functions $f_L$, $f_v$ to be non-trivial (e.g. Muschinski & de Bruyn Kops Reference Muschinski and de Bruyn Kops2015). However, it is at least plausible that a subregime exists wherein
either due to sufficiently high scale separation of isotropic scales, i.e. $ {{Re_s}} \gg 1$, or due to shared functional forms such that $f_L^{OC}=f_L^{K}$. Subject to either condition, it would be expected that (5.4) would reduce to
Indeed, similar arguments would motivate the $4/3$ relation between transverse and longitudinal velocity spectra in the presence of anisotropic outer scales such that $F_L^{K}(k_yL^{O}_{xy})\neq F_L^{K}(k_yL^{O}_{yy})$ in order to obtain
The validity of the relation (5.6) is tested with appropriately scaled one-dimensional spectra in figure 7(a), where we would expect the relation to hold in the locally isotropic wavenumber regime where $L_C k_y \sim O(1)$ (Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994). We observe a tendency for the local minima of the spectra to decrease with increasing Reynolds number. The minima for cases R8, R9 and R10 are all close to, and bounded below by 0.72, within the range of values reported for the universal constants (Sreenivasan Reference Sreenivasan1995; Sreenivasan & Kailasnath Reference Sreenivasan and Kailasnath1996). Furthermore, in the stratified boundary layer, Wyngaard & Coté (Reference Wyngaard and Coté1971) report measurements of each constant independently, in a flow configuration very similar to the one studied here, which imply $0.76 \pm 0.11$, also in agreement with the results observed here.
The wavenumbers associated with this regime also tend toward smaller scales as $ {{Re_s}}$ increases. In cases R8, R9 and R10, the lower bound for this regime appears to remain fixed at approximately $k_yL_C=0.4$ as the right bound begins to increase such that the flat region of the spectra broadens. We note that this apparently asymptotic behaviour in cases R8, R9 and R10 occurs while $ {{Re_s}}$ increases by a factor of two. The ratio of transverse to longitudinal one-dimensional spectra, as expressed in (5.7), is shown in figure 7(b). We observe a similar downward trend of the locally isotropic wavenumber regime with $ {{Re_s}}$, leading eventually to good agreement with the $4/3$ law for case R10.
Similarly to the data shown in figure 7(a), the locally isotropic regime at higher wavenumbers becomes more consistent with the Corrsin length scale at approximately $k_yL_C\approx 2$. The misalignment of the wavenumber regimes consistent with (5.6) might be seen to be a consequence of the anisotropy of the outer scales, as accounted for in our definitions of anisotropic length scales in $f_L^{K}$ and also in the observation that the outer scale of the scalar largely coincides with the Ozmidov length. Therefore, it is expected that $L_{OC}/L_C \sim {{Ri}}^{-3/4}$ and hence the wavenumber regime associated with the four-thirds relation is expected to be approximately a factor of four larger than (5.6), essentially as observed here.
We compare measured three-dimensional potential and kinetic energy spectra with the model spectra, (2.16) and (2.15), in figure 8(a,b). In both panels, a $|\boldsymbol {k}|^{-5/3}$ curve is plotted for reference. Our intention is to verify that the proposed spectra are approximated by our analytic model spectra at high $Re_s$, in the sense that we do not observe substantial systematic deviations which would otherwise indicate alternative two-point parameterisations would be more appropriate. Here, we observe a trend wherein the range of scales associated with three-dimensional turbulence, $|\boldsymbol {k}|L_C>1$, is approximated by the $-5/3$ power law as $ {{Re_s}}$ increases. Whereas deviations from the power-law behaviour are observed at all $ {{Re_s}}$, the quantitative sensitivity of the parameter $\varPi$ to deviations from the model spectra is unclear at this stage. Therefore, we next turn to evaluating the impact of Reynolds number on our key parameter $\varPi$.
5.1. Model verification
In the previous section, we have presented evidence that the ratio of the potential energy spectrum to a longitudinal energy spectrum appears to be consistent with the underlying scaling arguments central to coexistent Kolmogorov and Obukhov–Corrsin regimes. The value of the key parameter $\varPi$ can be expressed in terms of the one-dimensional universal constants, at high Reynolds number, by
For these flows, the parameter $\varPi$ should be assumed to be universal to the extent that $C$ and $\beta$ are. For reasons outlined in Sreenivasan (Reference Sreenivasan1991) and also observed here, measuring individual constants from a sheared flow is fundamentally problematic and, furthermore, measuring $\beta$ requires an accurate measurement of $\chi$, which is inherently difficult in experimental flows.
Nonetheless, estimates of the (three-dimensional) Kolmogorov and Obukhov–Corrsin constants in the literature would imply that $\varPi \approx 0.42$ (Wyngaard & Coté Reference Wyngaard and Coté1971; Sreenivasan Reference Sreenivasan1991; Sreenivasan & Kailasnath Reference Sreenivasan and Kailasnath1996; Yeung Reference Yeung2002; de Bruyn Kops Reference de Bruyn Kops2015) with estimates of the individual constants typically being reported within approximately 15 % of each other in the literature cited. Robust estimates of $\varPi$ from the relation (5.8) and empirical observations in figure 7 indicate $\varPi \approx 0.39$, which is certainly within the range of values reported in the literature. We demonstrate that the important expression (5.8) appears to be valid within the range of uncertainty for reported values of $\beta /C$ in figure 9. We also remark that the range of $\varPi$ is small with respect to values of $\beta /C$ for the Reynolds number space considered. Such an apparently limited range of $\varPi$ has also been reported by Venayagamoorthy & Stretch (Reference Venayagamoorthy and Stretch2006).
Furthermore, given the empirical observation that $L_{OC}=L_O$ (figure 6a) at sufficiently high Reynolds number,
Therefore, the stationary Froude number can be thought of as a simple consequence of the two universal constants at high Reynolds numbers, as validated in figure 9, which further implies $N_\ast \approx 1$, remembering (2.19). This finally suggests that the behaviour at lower Reynolds number may be parameterised by a function $f_\varPi$, i.e. at unity Prandtl number
6. Concluding remarks
Inspired by increasing evidence that high dynamic-range stratified turbulence exhibits scale similarity as proposed by Kolmogorov (Reference Kolmogorov1941), Oboukhov (Reference Oboukhov1949) and Corrsin (Reference Corrsin1951), we have considered the formal adaptation of a passive scalar mixing model (Beguier et al. Reference Beguier, Dekeyser and Launder1978) to stratified scalar mixing. By deriving the model from scaling theories, we demonstrate how model parameters have explicit functional dependence and are Reynolds number independent, at sufficiently high Reynolds number. In order to verify our hypothesis, we have employed a model flow (i.e. S-HSST) which equilibrates such that the effects of decreasing dissipation scales, here parameterised by increasing shear Reynolds number $ {{Re_s}}$ (as defined in (2.7a–c)), may be investigated independently of other flow parameters.
By direct numerical simulations of S-HSST, we have shown that the outer length scales of turbulence and the active scalar become parametrically coupled such that the relationship between energetics and dissipation rates is predicted by the scaling theories of Kolmogorov and Obukhov–Corrsin similarity. Therefore, the fundamental parameter associated with this relationship, $\varPi$, as defined in (2.17), is the ratio (5.8) of the Obukhov–Corrsin constant $\beta$ and the Kolmogorov constant $C$, at sufficiently high $ {{Re_s}}$, and should be considered universal to the extent that its constitutive constants are. As detailed in the context of various mixing models in § 2.2, robust characterisations of $\varPi$ have strong implications for the calibration of a number of turbulence and mixing models.
This universality has profound, but still not fully explained, implications for the turbulent flux coefficient $\varGamma$ of such steady stratified and sheared turbulence. From the expression (2.21), $\varGamma$ tending to a universal asymptotic value can now be understood as being due to $\varPi$ tending to a constant, which follows from the scaling theory and similarity arguments presented here, that in turn rely upon classical, and relatively well-established, turbulence modelling approaches. Therefore, to ‘understand’ physically and theoretically why the empirical modelling of Osborn & Cox (Reference Osborn and Cox1972) and Osborn (Reference Osborn1980) works so well (with $\varGamma \sim 0.2$), the remaining interesting open question to address is why the ratio of energies $R_{PK}=E_p/E_k \simeq 0.08$ for $Re_s \gtrapprox 50$, as shown in figure 4.
The parameter $\varPi$ may be stated exactly as a function of $N_\ast, {{S_\ast }} , {{Ri}}$. However, we demonstrate that it may be effectively reduced to a function of $ {{Fr}}$ at high Reynolds number, or a function of $ {{Fr}}$ and $ {{Re_s}}$ at low Reynolds numbers for the flows considered here. The transition of the high and low Reynolds number regimes is observed approximately at $ {{Re_s}}\approx 40$, or $ {{Re_b}}\approx 300$ at stationary Richardson number, for unity Prandtl number. We demonstrate that large-scale characteristics of the active scalar, as parameterised by $N_\ast$, exhibit significant dependence on the Reynolds number for $ {{Re_s}} \lessapprox 40$. The correlation between the transition point, in Reynolds number space, of the parameter $\varPi$ and $N_\ast$ supports the underlying assumptions of the proposed model.
A stated and significant limitation of the present study is the restriction of the model to $ {{Fr}} \sim O(1)$ and $Pr=1$. An objective of future work is to perform a similar procedure to § 2.2, except using scaling laws which account for expanded parameter regimes (for instance Batchelor Reference Batchelor1959; Lindborg Reference Lindborg2006; Kunze Reference Kunze2019).
Acknowledgements
We acknowledge the guidance of Professor J.J. Riley, the input of Dr F.A.V. de Bragança Alves and helpful discussions with Dr R. Ristorcelli and Dr J.A. Saenz.
Funding
This work was funded by the U.S. Office of Naval Research via grant N00014-15-1-2248. High performance computing resources were provided through the U.S. Department of Defense High Performance Computing Modernization Program by the Army Engineer Research and Development Center and the Army Research Laboratory under Frontier Project FP-CFD-FY14-007. This manuscript is approved for public release from Los Alamos National Laboratory as LA-UR-20-28841.
Declaration of interests
The authors report no conflict of interest.