2.1. Mathematical model

In order to assess the nature of turbulent mixing induced by a sheared and density stratified parallel flow, the initial height-dependent profiles for horizontal velocity and density are fixed to the following idealized and dimensional profiles:

(2.1a,b ) $$\begin{eqnarray}U(z)=U_{0}\tanh \left(\frac{z}{d}\right),\quad {\it\rho}(z)={\it\rho}_{0}-{\rm\Delta}{\it\rho}\tanh \left(\frac{z}{d}\right),\end{eqnarray}$$

in which the shear and density layers are both assumed to have thickness $2d$ . Further, $U_{0}$ and ${\it\rho}_{0}$ denote reference velocity and density while ${\rm\Delta}{\it\rho}$ fixes the density difference across the layer. The coordinate system is chosen such that $x,y$ and $z$ are aligned with the streamwise, spanwise and vertical directions respectively.

Subject to the Boussinesq approximation and in dimensionless form, the conservation equations become

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}u_{i}}{\text{D}t}=-\frac{\partial p}{\partial x_{i}}-\mathit{Ri}_{0}{\it\rho}{\it\delta}_{i3}+\frac{1}{\mathit{Re}}\frac{\partial ^{2}u_{i}}{\partial x_{j}^{2}}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}}{\partial x_{i}}=0, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}{\it\rho}}{\text{D}t}=\frac{1}{\mathit{Re}\,\mathit{Pr}}\frac{\partial ^{2}{\it\rho}}{\partial x_{j}^{2}}, & \displaystyle\end{eqnarray}$$

$i,j=1,2,3$

(2.5a−c ) $$\begin{eqnarray}\mathit{Re}=\frac{U_{0}d}{{\it\nu}},\quad \mathit{Ri}=\frac{N^{2}}{S^{2}}=-\frac{g}{{\it\rho}_{0}}\frac{\partial {\it\rho}}{\partial z}\left(\frac{\partial u}{\partial z}\right)^{-2},\quad \mathit{Pr}=\frac{{\it\nu}}{{\it\kappa}}.\end{eqnarray}$$

The minimum value in the initial profile of $\mathit{Ri}$ , obtained from the initial density and velocity profiles in (2.1), occurs at the midpoint of the shear layer (at $z=0$ ). We denote this specific choice of gradient Richardson number by $\mathit{Ri}_{0}$ , which also appears as the coefficient of the non-dimensional buoyancy term in (2.2). It should be noted that while the value of $\mathit{Ri}_{0}$ is important as a necessary condition for the growth of the primary KH instability, bulk measures of the Richardson number (e.g. by vertically averaging $N^{2}$ and $S^{2}$ ) are more relevant insofar as a comparison with oceanic observations is concerned and are usually much higher than $\mathit{Ri}_{0}$ . Also, in the above non-dimensionalization ${\it\rho}_{0}$ , $U_{0}$ and $d$ (half of the shear layer thickness) have been chosen as the characteristic density, velocity and length scales. In contrast to theoretical studies, observational oceanographers commonly employ the full shear layer thickness as the proper choice of length scale as well as $2{\it\rho}_{0}$ and $2U_{0}$ as the characteristic density and velocity, resulting in Reynolds numbers that are four times higher than the values employed in this study (e.g. Smyth et al. Reference Smyth, Moum and Caldwell2001).

2.2. Direct numerical simulations

The three-dimensional DNS of stratified turbulent flows to be discussed in the main body of this paper were performed using the modern spectral element code, Nek5000, originally developed at Argonne National Laboratory (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008), which is especially suited to simulations of transitional and turbulent flows, including those that develop in complex geometries. The methodology has been shown to be characterized by minimal numerical dissipation and dispersion, as described in Fischer, Kruse & Loth (Reference Fischer, Kruse and Loth2002). As is most common in the application of Galerkin methods, a weak formulation of the governing equations is employed, and these are discretized in space using $N$ th-order Lagrange polynomial interpolants following the $P_{N}{-}P_{N-2}$ formulation of Maday, Patera & Rønquist (Reference Maday, Patera and Rønquist1990) for the velocity ( $N$ ) and pressure ( $N-2$ ) spaces respectively. The high-order polynomials employed are based on the tensor product of Gauss–Lobatto–Legendre (GLL) quadrature points in each of $K$ hexahedral elements, a procedure that produces a number of grid points equal to  $KN^{3}$ . The decoupled velocity and pressure fields are solved iteratively using conjugate gradients and GMRES methods, which in turn employ as preconditioners the scalable Jacobi and additive overlapping Schwarz methods, respectively (Fischer Reference Fischer1997). The latter sparse linear solver exploits the tensor product structure with highly optimized matrix–matrix ( $\text{m}\times \text{m}$ ) routines as well as a parallel coarse-grid solver, scaling to thousands of processors (Tufo & Fischer Reference Tufo and Fischer2001). Such scalability is critical in the present context, in which the interest is in stratified turbulent flows in the high- $\mathit{Re}$ , high- $\mathit{Pr}$ regime. The time advance of the Boussinesq equations is based on use of third-order semi-implicit operator-splitting methods (Maday et al. Reference Maday, Patera and Rønquist1990; Fischer Reference Fischer1997). At high $\mathit{Re}$ , in order to achieve numerical stability while maintaining solution accuracy, Nek5000 employs dealiasing through application of the $3/2$ rule (Deville, Fischer & Mund Reference Deville, Fischer and Mund2002), as well as polynomial filtering (Fischer & Mullen Reference Fischer and Mullen2001). In the latter approach, a fraction ${\it\beta}$ of the accumulated energy associated with the highest wavenumbers is transferred to the third-highest order polynomial through which the high-frequency oscillations are damped. In the present study, all computations to be reported have been carried out using both dealiasing and 5 % filtering ( ${\it\beta}=0.05$ ). Benchmark simulations with no filtering resulted in negligible change to the mean turbulent quantities. The Nek5000 code has been employed previously in a limited number of problems in environmental hydrodynamics of stratified flows, and has been shown to perform extremely well, for example in the reproduction of a benchmark result for a lock-exchange problem (Özgökmen, Iliescu & Fischer Reference Özgökmen, Iliescu and Fischer2009) and in the representation of entrainment processes in density currents and overflows (Özgökmen et al. Reference Özgökmen, Fischer, Duan and Iliescu2004a ,Reference Özgökmen, Fischer, Duan and Iliescu b ).

For the purposes of the current work the computational domain is taken to consist of a parallelepiped of dimensions $(L_{x},L_{y},L_{z})$ which is discretized using uniformly spaced hexahedral elements in the streamwise ( $x$ ) and spanwise ( $y$ ) directions. The grid elements in the vertical direction are taken to be more finely spaced (uniformly) in the vicinity of the inflection point in the vertical profile of horizontal velocity over the non-dimensional range $-5\leqslant z\leqslant 5$ . In the region above and below this innermost region of enhanced resolution the vertical spacing between adjacent elements in the grid is gradually stretched by 1.25 % between successive elements. The innermost region in the vertical generously accommodates the entire height of the evolving Kelvin–Helmholtz billow (mainly controlled by the value of $\mathit{Ri}_{0}$ ) and the surrounding fluid that it entrains, while providing the necessary resolution where required for accurate representation of the evolving turbulence. The smallest energized scales of the turbulent flows are determined based on the choices of the non-dimensional control parameters $\mathit{Re}$ and $\mathit{Pr}$ . At the high Reynolds and Prandtl numbers that are of interest to us in this paper, structures persist all the way down to the dissipation scale (and whose contributions during the transition phase into turbulence can be significant), so compromise is required through the introduction of an appropriate cut-off scale in order to numerically simulate such flows. Nonetheless, we expect that as long as the required resolutions are constructed consistently for different values of the non-dimensional control parameters, the results obtained will be of high quality. In this regard it is now generally accepted that a reliable DNS simulation should resolve down to 3–6 times the dissipation scale (see e.g. Moin & Mahesh Reference Moin and Mahesh1998). When the momentum and scalar diffusivities are the same (i.e. $\mathit{Pr}=1$ ), the dissipation scale is characterized by the Kolmogorov length scale where $L_{K}=({\it\nu}^{3}/{\it\varepsilon})^{1/4}$ , in which ${\it\nu}$ is the kinematic viscosity ( ${\it\nu}=U_{0}d/\mathit{Re}$ in our non-dimensionalization) and ${\it\varepsilon}$ denotes the dimensionless viscous dissipation. For higher Prandtl numbers, however, the resolution requirement becomes significantly more stringent, as the smallest scales involved in the diffusion of density are then on the order of the Batchelor length scale  $L_{B}=(L_{K}/\mathit{Pr}^{1/2})$ (see e.g. Smyth & Moum Reference Smyth and Moum2000b for a recent discussion). On these grounds we have elected to design all of our numerical simulations consistently, based on an a priori estimation of $L_{B}$ using a series of benchmark simulations. All of the results to be reported in this paper were verified a posteriori as having resolved length scales extending down to approximately $4.5L_{B}$ , and at least down to $2.5L_{K}$ during the most turbulently active stages of flow evolution, and scales as small as approximately $2L_{B}$ outside these most active stages of flow evolution. Note that for the purpose of these calculations of the resolution requirements, the total viscous dissipation ${\it\varepsilon}$ (see (2.14)) is averaged over the range $-1\leqslant z\leqslant 1$ (rather than the whole domain), where presumably the most energetic turbulent activity occurs. For further discussion of this topic please refer to the supplementary data accompanying this paper, available at http://dx.doi.org/10.1017/jfm.2015.225. Table 1 summarizes the details of the set of different numerical experiments for which results will be presented below. The required extent of the computational domain has been chosen based on the previously described numerical analyses of Mashayek et al. (Reference Mashayek, Caulfield and Peltier2013) and Mashayek & Peltier (Reference Mashayek and Peltier2013). All of the analyses to be reported in this paper have been conducted at a fixed nominal initial Reynolds number of $\mathit{Re}=6000$ , which leads to a Taylor microscale Reynolds number of $\mathit{Re}_{{\it\lambda}}=O(10^{2})$ (see appendix A for further discussion) during the most intense phase of turbulence. At $\mathit{Pr}=1$ , an initial Reynolds number of $\mathit{Re}=6000$ has been shown to be high enough for mixing efficiency to remain unchanged for higher $\mathit{Re}$ and also to prohibit the pairing instability and the upscale cascade of energy which this vortex merger instability supports, as demonstrated in Mashayek & Peltier (Reference Mashayek and Peltier2013). This latter conclusion will be re-established in this paper to be as appropriate at high Prandtl number as it has been shown to be at $\mathit{Pr}=1$ . Thus for the primary mixing investigations, it will suffice to set the streamwise extent of our domain to include only a single wavelength ( $L_{x}=14.27$ ) of the fastest growing mode (FGM) of inviscid instability. In order to investigate the pairing instability, of course, we will set the streamwise extent to $2L_{x}$ . The most demanding simulation in the sequence that we will describe, that for $\mathit{Pr}=16$ , required approximately 5 701 632 core-hours on a BG/Q machine, which is approximately equivalent to a DNS study of a case with $\mathit{Pr}=1$ and $\mathit{Re}=40\,000$ , or $\mathit{Pr}=700$ and $\mathit{Re}=500$ (assuming the same numerical requirements concerning spatial resolution), all of which are unprecedented regimes that have not previously been investigated in the literature but which are highly relevant to the understanding of diapycnal mixing in the oceans.

Table 1.

Details of the three-dimensional numerical experiments in which the total grid points are about $N^{3}K_{x}K_{y}K_{z}$ . Also, $K_{z}^{c}$ represents the number of uniform elements resolving the innermost region of the domain in the vertical direction ( $-5\leqslant z\leqslant 5$ ).

In these simulations, flow evolution is initialized by adding a small two-dimensional perturbation in the form of the FGM of the Taylor–Goldstein equation (see e.g. Klaassen & Peltier Reference Klaassen and Peltier1985c ; Kundu & Cohen Reference Kundu and Cohen2008). The amplitude of this initial normal mode perturbation is set to 0.01, which translates into a kinetic energy of the seed two-dimensional disturbance that is $10^{-6}$ times smaller than the kinetic energy of the background flow. We also experimented with an alternative form for the initial perturbation with a shape function $\exp (-z^{2})$ , and discovered that while the results (including the shape of the billow and the sequence of secondary instabilities that arise) differ significantly at a low Reynolds number of $\mathit{Re}=750$ , they become noticeably insensitive to such choice of initial perturbation at sufficiently high $\mathit{Re}$ , which is the most relevant for the geophysical applications of interest to us. Thus at high $\mathit{Re}$ , the only difference incurred as a consequence of any such difference in the initialization procedure involves the computational time required to reach the saturation state for the initial KH billows, where our choice of FGM as the initial perturbation tends to minimize the required computational resources. In addition to the two-dimensional perturbation in the streamwise vertical plane ( $x{-}z$ ), it is absolutely critical to initially perturb the spanwise velocity with small-amplitude white noise in order for the three-dimensional instabilities to grow without bias from the initial conditions. The amplitude of such spanwise perturbations is set to 0.001, which translates into a level of three-dimensional disturbance that is initially one order of magnitude smaller than the two-dimensional disturbance. In the absence of such spanwise perturbation, the essential growth of three-dimensional instabilities would eventually be influenced by truncation errors at the grid scale and thus the flow field would become grid-dependent and hence physically implausible. Insofar as boundary conditions are concerned, the velocity and density fields are assumed to be periodic in both streamwise and spanwise directions, while the top and bottom boundaries are assumed to be free-slip and impermeable for velocity and insulated (zero density flux) for the density field.

2.3. Energetics and mixing analysis of three-dimensional DNS data

In order to analyse the flow energetics during the course of its evolution, we employ a Reynolds decomposition of the flow field into a mean (associated with the background flow) and perturbation components. Furthermore, the perturbation components of the velocity and density fields are themselves further subdivided into a spanwise averaged component associated with two-dimensional eddies (mainly associated with the initial KH wave) as well as an inherently three-dimensional component (Klaassen & Peltier Reference Klaassen and Peltier1991; Caulfield & Peltier Reference Caulfield and Peltier2000), expansions that we may write explicitly as

(2.6) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}(x,y,z,t)=\overline{U}(z,t)+\boldsymbol{u}^{\prime }(x,y,z,t), & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}^{\prime }(x,y,z,t)=[u_{2d}+u_{3d},v_{3d},w_{2d}+w_{3d}], & \displaystyle\end{eqnarray}$$

$2d$

$3d$

(2.8) $$\begin{eqnarray}\displaystyle & \overline{U}(z,t)=\langle \boldsymbol{u}\rangle _{xy}, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & (u_{2d},0,w_{2d})(x,z,t)=\langle (u-\overline{U},v,w)\rangle _{y}, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & (u_{3d},v_{3d},w_{3d})(x,y,z,t)=(u-\overline{U}-u_{2d},v,w-w_{2d}) & \displaystyle\end{eqnarray}$$

$\langle .\rangle _{p}$

$p$

An appropriate kinetic energy budget can be obtained by taking the inner product of (2.2) with the total velocity vector. The derivation is straightforward, and common in the derivation of turbulent balance equations (see e.g. Kundu & Cohen Reference Kundu and Cohen2008). The resulting evolution equation for the total kinetic energy $\mathscr{K}=\langle u^{2}+v^{2}+w^{2}\rangle _{xyz}/2$ becomes (Peltier & Caulfield Reference Peltier and Caulfield2003)

(2.11) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\mathscr{K}=-\mathscr{H}+\mathscr{D}.\end{eqnarray}$$

The total kinetic energy in (2.11) is either transferred into buoyancy flux $\mathscr{H}$ , thereby increasing the potential energy of the system (i.e. it appears with a reversed sign in (2.15) below), or dissipated through viscous stresses at a rate $\mathscr{D}$ (which is negative definite), and these have the explicit definitions

(2.12) $$\begin{eqnarray}\displaystyle & \mathscr{H}=\mathit{Ri}_{0}\langle {\it\rho}w\rangle _{xyz}, & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \mathscr{D}=\langle {\it\varepsilon}\rangle _{xyz}, & \displaystyle\end{eqnarray}$$

${\it\varepsilon}$

(2.14) $$\begin{eqnarray}{\it\varepsilon}=-\frac{2}{\mathit{Re}}(\unicode[STIX]{x1D634}_{ij}\unicode[STIX]{x1D634}_{ij}),\end{eqnarray}$$

in which $\unicode[STIX]{x1D634}_{ij}=(\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i})/2$ is the strain rate tensor.

The total averaged potential energy of the system, on the other hand, is defined as $\mathscr{P}=\mathit{Ri}_{0}\langle {\it\rho}z\rangle _{xyz}=\mathit{Ri}_{0}\langle z\overline{{\it\rho}}\rangle _{z}$ (Caulfield & Peltier Reference Caulfield and Peltier2000). Thus the time evolution of $\mathscr{P}$ is governed by the following evolution equation:

(2.15) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\mathscr{P}=\mathscr{H}+\mathscr{D}_{p},\end{eqnarray}$$

where $\mathscr{H}$ is as defined in (2.12), and $\mathscr{D}_{p}$ represents the irreversible conversion of fluid internal energy into potential energy through the action of diffusive processes acting preferentially upon the smallest scales of the density heterogeneity in a statically stable, motionless environment, and is hence strictly positive. This contribution to the growth of potential energy may be approximated as (Caulfield & Peltier Reference Caulfield and Peltier2000)

(2.16) $$\begin{eqnarray}\mathscr{D}_{p}=\frac{\mathit{Ri}_{0}{\rm\Delta}{\it\rho}}{L_{z}\mathit{Re}\mathit{Pr}},\end{eqnarray}$$

where ${\rm\Delta}{\it\rho}$ is the density difference across the computational domain of height $L_{z}$ .

It is important to note that not all the exchange of energy between the kinetic and potential energy ‘reservoirs’ that occurs as a consequence of the action of the buoyancy flux $\mathscr{H}$ can be interpreted as irreversible mixing. In fact, this buoyancy flux is composed of two distinct components, the first of which we refer to as the stirring rate, $\mathscr{S}$ , and the second we identify as the mixing rate, $\mathscr{M}$ , which are associated with reversible and irreversible processes respectively. Such partitioning may be explicitly expressed as

(2.17) $$\begin{eqnarray}\mathscr{H}=\mathscr{M}+\mathscr{S}.\end{eqnarray}$$

This distinction between reversible stirring and irreversible mixing is fundamental to the original definition of available potential energy (APE) by Lorenz (Reference Lorenz1955), whose interest was in the application of this concept to the understanding of the general circulation of the atmosphere (see Peltier & Caulfield Reference Peltier and Caulfield2003 for a review). Following Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) and Caulfield & Peltier (Reference Caulfield and Peltier2000), who applied this concept to the understanding of small-scale turbulent flows, one can subdivide the total averaged potential energy of the flow into APE, ( $\mathscr{P}_{A}$ ), and background potential energy or BPE, ( $\mathscr{P}_{B}$ ), by writing

(2.18) $$\begin{eqnarray}\mathscr{P}=\mathscr{P}_{A}+\mathscr{P}_{B}.\end{eqnarray}$$

In other words, at any instant of time in the course of flow evolution, the APE involved in reversible stirring of the turbulent flow is the difference between the actual potential energy of the flow field ( $\mathscr{P}$ ) and the notional statically stable background state ( $\mathscr{P}_{B}$ ). It follows from (2.15), (2.17) and (2.18) that the APE for ‘stirring’, which is available to be converted back into kinetic energy, evolves as $(\text{d}/\text{d}t)\mathscr{P}_{A}=\mathscr{S}=\mathscr{H}-\mathscr{M}$ and all the irreversible processes contribute to increase the background potential energy as $(\text{d}/\text{d}t)\mathscr{P}_{B}=\mathscr{M}+\mathscr{D}_{p}$ . So the amount of irreversible mixing is obtained from the background potential energy as

(2.19) $$\begin{eqnarray}\mathscr{M}=\frac{\text{d}}{\text{d}t}\mathscr{P}_{B}-\mathscr{D}_{p}.\end{eqnarray}$$

Here $\mathscr{P}_{B}=\mathit{Ri}_{0}\langle {\it\rho}_{B}z\rangle$ represents the minimum (background) average potential energy that would be obtained if the fluid flow were strictly statically stable with a mean background density profile of ${\it\rho}_{B}(z,t)$ . Finding ${\it\rho}_{B}(z,t)$ involves two steps: (i) adiabatic rearrangement of the fluid parcels according to their individual densities into a monotonically decreasing upwards assembly, and (ii) sequentially deforming the shape of individual fluid parcels and stacking them as individual isopycnal sheets in order to conserve mass (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995; Caulfield & Peltier Reference Caulfield and Peltier2000). Such distortion, in our discretized space, converts a fluid parcel with mass ${\it\delta}m={\it\rho}{\it\delta}x{\it\delta}y{\it\delta}z$ into a sheet with thickness ${\it\delta}z_{B}={\it\delta}m/({\it\rho}_{B}L_{x}L_{y})$ .

Now the usual approach whereby this adiabatic rearrangement of the density field may be produced has been to numerically sort ${\it\rho}(x,y,z,t)$ at every time step in the evolution of the flow (see for example Caulfield & Peltier Reference Caulfield and Peltier2000; Peltier & Caulfield Reference Peltier and Caulfield2003 and Mashayek & Peltier Reference Mashayek and Peltier2013). However, for large simulations such as those required for the purpose of the current study, such serial sorting algorithms become prohibitively expensive (or even impossible) to employ because of the computational time and memory requirements. An alternative and superior approach is one based on the use of probability density functions (PDFs), proposed by Tseng & Ferziger (Reference Tseng and Ferziger2001), which we have implemented for present purposes (for a simple MATLAB code illustrating this approach for a one-dimensional profile, see the supplementary data). The idea is to discretize the probability sample space of density $\tilde{{\it\rho}}$ into an appropriate number of a priori selected bins. The probability that ${\it\rho}(x,y,z,t)\in [\tilde{{\it\rho}},\tilde{{\it\rho}}+\text{d}\tilde{{\it\rho}}]$ can be defined as a volume integral of the delta function as

(2.20) $$\begin{eqnarray}\text{PDF}(\tilde{{\it\rho}})=\frac{1}{V}\int _{V}{\it\delta}(\tilde{{\it\rho}}-{\it\rho})\,\text{d}V.\end{eqnarray}$$

Once this PDF reference state has been constructed by scanning the entire domain, a corresponding reference distribution of height levels ( $z_{B}$ ) that is associated with a monotonically decreasing upwards density profile ( ${\it\rho}_{B}$ ) is obtained as (note that $L_{x}L_{y}{\it\delta}z_{B}=V\,\text{PDF}(\tilde{{\it\rho}})\,\text{d}\tilde{{\it\rho}}$ )

(2.21) $$\begin{eqnarray}z_{B}({\it\rho}_{B})=L_{z}\int _{{\it\rho}_{B}}^{{\it\rho}_{M}}\text{PDF}(\tilde{{\it\rho}})\,\text{d}\tilde{{\it\rho}}.\end{eqnarray}$$

Hence the background density profile ( ${\it\rho}_{B}$ ) is assumed constant between adjacent height levels (here ${\it\rho}_{M}$ denotes the maximum density). In order for this method to work properly, especially for high-gradient profiles such as $\tanh (z)$ , it is critical to employ a non-uniform mesh spacing in the sample space so as to concentrate the bins in regions with high density gradient. A Chebyshev transformation is therefore used in the present study, where the non-uniform spacing is calculated as $\text{d}\tilde{{\it\rho}}_{j}=\text{d}\tilde{{\it\rho}}_{u}\cos {\it\theta}_{j}$ , with $\text{d}\tilde{{\it\rho}}_{u}$ being the presumed uniform spacing and ${\it\theta}=0,\dots ,{\rm\pi}/2$ near both tails of $\tanh$ and ${\it\theta}=0,\dots ,{\rm\pi}$ in the central smoother region.

Figure 1 provides a schematic that illustrates the process of finding PDFs in a parallel computation. Having constructed a user-defined non-uniform distribution of bins (using a Chebyshev transformation) associated with a discretized probability sample space characterized by a monotonically decreasing upward assembly (figure 1 b), every computational core (depicted by gridded cubes in figure 1 a) searches among the grid elements under its control and assigns the density of each element to the appropriate bin of an a priori discretization in order to calculate the complete PDF using (2.20). Finally a global sum among all processors determines the total PDF in every bin, which is then interpolated to form the desired background density profile (figure 1 c).

Figure 1.

Schematic illustration of a parallel computation of the sorted background density profile ${\it\rho}_{B}(z_{B})$ using the probability density functions as defined in (2.20) and (2.21). The numbers inside the bins in panel (b) are associated with each individual CPU of a cluster. The heights of the bins indicate the PDF (%) and their widths show the total thickness that is associated with the stacked layers of density that belong to that bin, i.e.  ${\it\rho}(x,y,z,t)\in [\tilde{{\it\rho}},\tilde{{\it\rho}}+\text{d}\tilde{{\it\rho}}]$ .

Having isolated the amount of change in $\mathscr{P}$ that is associated with mixing through (2.18), one can calculate the mixing efficiency as proposed in Caulfield & Peltier (Reference Caulfield and Peltier2000). Combining (2.11) and (2.17), the rate of change of kinetic energy of the flow, $(\text{d}/\text{d}t)\mathscr{K}$ is controlled by stirring, mixing and viscous dissipation processes ( $(\text{d}/\text{d}t)\mathscr{K}(t)=-\mathscr{S}-\mathscr{M}+\mathscr{D}$ ), where the last two terms are irreversible. So the instantaneous efficiency of mixing $\mathscr{E}_{i}$ may be defined as the ratio of ‘useful’ expenditure of $\mathscr{K}$ to the total irreversible loss of $\mathscr{K}$ , i.e.

(2.22) $$\begin{eqnarray}\mathscr{E}_{i}=\frac{\mathscr{M}}{\mathscr{M}-\mathscr{D}}.\end{eqnarray}$$

3.2. Secondary instabilities

We have conducted our secondary stability analysis on a series of background non-parallel shear flows obtained using two-dimensional DNS simulations at $\mathit{Re}=6000$ , $\mathit{Ri}_{0}=0.12$ and $\mathit{Pr}=1,8,16$ . These background states included not only the initial parallel velocity shear flow, but also the rolled-up ‘spaghetti-on-a-fork-like’ density and vorticity fields of the KH wave component of the flow. Since the analysis assumes that the background state is frozen in time (i.e. quasi-steady assumption), we have performed this analysis for several snapshots in time about the saturation point of the primary KH wave. The selected times are ideally suited to the purpose of our stability analysis since the growth rate of three-dimensional perturbations during the saturation of the two-dimensional KH wave is much faster than that associated with the background state (i.e.  ${\it\sigma}_{3d}=O(-1)\gg {\it\sigma}_{kh}=O(-4)$ ), thus satisfying our fundamental assumption of time scale separation. Furthermore, inasmuch as the two-dimensional flow at this high $\mathit{Re}$ becomes turbulent even without adding the spanwise three-dimensional perturbations (unlike lower $\mathit{Re}$ experiments, which stay laminar in a two-dimensional setting and only become turbulent due to the added three-dimensional perturbations, as discussed in Klaassen & Peltier Reference Klaassen and Peltier1985b ), it is important to choose a limited range of time slices at which the background state has not undergone such unrealistic two-dimensional turbulence. The selected background states are proximate in time to the saturation time of the KH wave and so also satisfy the latter consideration.

This study extends the previous work of Mashayek & Peltier (Reference Mashayek and Peltier2012a ,Reference Mashayek and Peltier b ) by exploring the stability of secondary modes at a higher Reynolds number for the background two-dimensional flows whose stability characteristics are to be investigated (their highest $\mathit{Re}$ was 2000) and, more significantly, by investigating the influence of elevating molecular Prandtl number above the value of unity. As a result of this shift in the region of parameter space of interest to us, and in particular due to a higher $\mathit{Re}$ , the conventional secondary instabilities studied so far will be seen to play more complicated roles in mediating the transition to turbulence, and to be modified importantly in terms of the physical mechanisms they represent. Understanding these differences in the operative secondary instabilities in the high- $\mathit{Re}$ , high- $\mathit{Pr}$ regime forms the main focus of this subsection. Besides the more conventional instabilities, we have also identified further inhabitants of the ‘zoo’ of such modes, but these will not be described here.

The conventionally understood secondary instabilities in the context of growing three-dimensional perturbations on a stratified KH wave include (i) the secondary shear-aligned convective instability (Klaassen & Peltier Reference Klaassen and Peltier1985c ) which produces shear-aligned streamwise vortex tubes of alternating sign in the spanwise direction in the billow eyelids, (ii) the secondary shear instability of the braid (manifested by inclined mini-replicas of the primary KH wave on the braid) originally discussed by Corcos & Sherman (Reference Corcos and Sherman1976), and (iii) the SPI (Mashayek & Peltier Reference Mashayek and Peltier2012a ) which is located at the stagnation point of the braid. In all the cases for which we have conducted non-separable stability analysis (for various selected times as well as the full bandwidth of the spanwise wavenumber), the SPI appeared to be the FGM with the secondary shear-aligned convective mode being next in the list, arranged in decreasing order of growth rate. Interestingly, the eigenfunction associated with the secondary shear instability on the braid does rise in prominence briefly as the FGM prior to its being replaced in that role by the SPI, and was never found to persist significantly in time. In the following subsections we will focus on the characteristics of these important instabilities in the parameter space of interest to us, by categorizing them into ‘braid instability’ and ‘eyelid instability’ types.

3.2.1. Braid instability

The SPI dominated the braid region in all of our two-dimensional simulations, while the shear instability on the braid did not emerge at $\mathit{Re}=6000$ until the braid was perturbed by the surrounding instabilities at later times during flow evolution. Our secondary instability analysis also established similar results regarding the dominance of SPI, while the other braid instability could be most prominently realized at times prior to the actual emergence of SPI. Figure 2(a,b) illustrates the evolving two-dimensional background states before and after the emergence of SPI for a case with $\mathit{Pr}=1$ . Figure 2(c,d) depicts the eigenfunction associated with the perturbation kinetic energy of the FGM at these times. The eigenfunction at $t=64$ suggests KH-like vortices on the braid with strong perturbations along the braid representing the secondary shear instability. Nevertheless, at later times in the flow evolution this eigenfunction transitions into a more localized mode which is focused upon the stagnation point, and characterized by weaker perturbations (albeit non-zero) along the braid. As has been previously discussed in Mashayek & Peltier (Reference Mashayek and Peltier2012a ), the emergence of these two modes depends on the relative strength of the strain and shear fields when mapped into a braid-aligned coordinate system. For our $\mathit{Re}=6000$ analysis, the strength of the braid strain field must have initially overwhelmed velocity shear normal to the braid such that the secondary shear instability was inhibited and the SPI was dominant. With this view in mind, these two instabilities may be viewed simply as comprising two branches of the ‘braid instability’.

Figure 2.

Demonstration of the braid instability (FGM) for $\mathit{Pr}=1$ at $t=64$ (a,c) and $t=70$ (b,d) corresponding to times before and after the emergence of SPI. (a,b) The density contour plots and (c,d) the perturbation kinetic energy eigenfunction that is associated with this mode of instability.

The competition between the shear and strain deformation fields on the braid is subtle. Based on an equilibrium similarity theory, Smyth (Reference Smyth2003) suggested that increasing $\mathit{Pr}$ directly increases the shear across the braid-aligned coordinate system by further sharpening the density gradients. As will be discussed in the following section, the molecular Prandtl number also indirectly influences the strain field such that it suppresses the growth rate of SPI, resulting in a postponed emergence of this mode. Figure 3(a) plots the variation of growth rate ${\it\sigma}_{3d}$ with respect to spanwise wavenumber $d$ for $\mathit{Pr}=1,8,16$ obtained from a background state at the time of saturation of KH wave amplitude, $t_{2d}^{s}=74$ . As discussed by Mashayek & Peltier (Reference Mashayek and Peltier2012a ) and also shown in this figure, this mode is highly three-dimensional as it exists at a wide range of cross-stream wavelengths, unlike the shear instability on the braid which is primarily two-dimensional. Moreover, at this time, the growth rate associated with SPI decreases at almost all wavelengths as $\mathit{Pr}$ increases. In order to fully assess the suppression of SPI due to higher $\mathit{Pr}$ , figure 3(b) compares the time variation of ${\it\sigma}_{FGM}$ associated with the wavenumber of maximum growth rate for $\mathit{Pr}=1$ and $\mathit{Pr}=16$ . This plot suggests that for $\mathit{Pr}=1$ the highest growth rate is sustained for a longer period of time, whereas for $\mathit{Pr}=16$ , ${\it\sigma}_{FGM}$ reaches similar rates only briefly. Since SPI extracts its energy from the strain deformation field of the braid, its suppression due to higher $\mathit{Pr}$ must have been caused by weakening of this field. As will be suggested in the next section, the strain and shear fields are indirectly influenced by $\mathit{Pr}$ due to the promotion of the ‘eyelid instability’, which modifies these fields in such a way that the emergence of SPI is postponed. We will also confirm the suppression of SPI due to an increase in the molecular Prandtl number based on our three-dimensional DNS results below.

Figure 3.

The suppression of SPI due to higher $\mathit{Pr}$ at a fixed $\mathit{Re}$ and $\mathit{Ri}$ . (a) The variation of growth rate ${\it\sigma}_{3d}$ associated with FGM at $t_{2d}^{s}$ for the full bandwidth of spanwise wavenumbers. (b) Time variation of ${\it\sigma}_{3d}$ at the wavenumber associated with its maximum rate.

Table 2 summarizes the growth rates associated with SPI at the spanwise wavenumber of maximum growth rate for $\mathit{Pr}=1$ and $\mathit{Pr}=16$ . Moreover, the energetic contribution of different mechanisms to the growth of perturbation kinetic energy as introduced in (3.6) is also reported as a percentage in this table. Notice that the share of $\mathscr{H}_{3d}$ increases monotonically in time for all $\mathit{Pr}$ , highlighting the three-dimensional nature of this mode of instability, which is enhanced in the course of flow evolution.

Table 2.

The growth rate ${\it\sigma}_{3d}$ and energetic properties (in percentage and as introduced in (3.6)) of the most unstable three-dimensional mode associated with the SPI.

3.2.2. Eyelid instability

In all our two-dimensional simulations the secondary shear-aligned convective instability dominated the billow eyelid (as opposed to the braid). Besides the conventional picture for this mode, which consists of a periodic array of streamwise vortices in the spanwise direction and are convectively driven, we have discovered a new characteristic associated with this mode of instability at high Reynolds and Prandtl number, consisting of spanwise vortices that are counter-rotating in the vertical direction and are shear-driven. The latter property for this mode is new and has not been realized in the previously published literature on the stability analysis of these secondary modes, primarily due to the much lower values of $\mathit{Re}$ at which these previous analyses have been performed. This new characteristic prompts us to refer to this mode simply as an ‘eyelid instability’. In this section, we will first analyse the impact of $\mathit{Pr}$ on the convectively driven component of the instability by investigating its influence on the Rayleigh number of the most statically unstable region of the billow. We will then explain the shear-driven counter-rotating vortex component of the mode that seems to follow the development of the convectively driven component and eventually appears most prominently around the diagonal corners of the eyelid(s). Finally in this section we will address the critical role of $\mathit{Pr}$ in transforming this mode into a short-wave instability, capable of directly injecting turbulence into small scales.

Figure 4.

The contour plots of buoyancy frequency $N^{2}$ for (a)  $\mathit{Pr}=1$ and (b)  $\mathit{Pr}=16$ at $t_{2d}^{s}$ in which the shaded regions illustrate the statically unstable SARs in which $N^{2}<0$ .

Figure 4 illustrates the buoyancy frequency contours associated with $\mathit{Pr}=1$ and $\mathit{Pr}=16$ at $t_{2d}^{s}$ . The super adiabatic regions (SARs) in which the density distribution is statically unstable are shaded in this figure. As shown in this plot, the degree of static instability associated with these regions is enhanced as $\mathit{Pr}$ increases, which is mainly a consequence of a higher density gradient across the SARs. This is because the decreased density diffusion inside the core at higher values of $\mathit{Pr}$ tends to only weakly thin the SAR but to significantly increase the density contrast across it (Klaassen & Peltier Reference Klaassen and Peltier1985a ). In analogy with the Rayleigh–Bénard problem, we define a Rayleigh number representative of the outermost statically unstable region that is most likely to exceed the critical value (which is usually of order $10^{3}$ ) for the onset of convective instability (Klaassen & Peltier Reference Klaassen and Peltier1985c ) as

(3.11) $$\begin{eqnarray}\mathit{Ra}=\mathit{Re}^{2}\,\mathit{Ri}\,\mathit{Pr}\,{\rm\Delta}{\it\rho}^{\star }{\it\delta}^{3},\end{eqnarray}$$

where ${\rm\Delta}{\it\rho}^{\star }$ represents the dimensionless density contrast across the outermost SAR and ${\it\delta}$ is its dimensionless thickness. This expression might suggest that $\mathit{Ra}\propto \mathit{Pr}$ . However, as mentioned earlier, ${\rm\Delta}{\it\rho}^{\star }$ and ${\it\delta}$ are both functions of $\mathit{Pr}$ , and for a fluid with higher molecular Prandtl number they have opposing effects on the Rayleigh number. In order to understand these intricate connections, table 3 contains the measured values of ${\rm\Delta}{\it\rho}^{\star }$ and ${\it\delta}$ for a vertical profile in the middle of the billow as well as the calculated Rayleigh number using (3.11). As expected ${\rm\Delta}{\it\rho}^{\star }$ increases for higher $\mathit{Pr}$ , primarily due to reduced density diffusion inside the billow core at the bottom edge of the SAR, while the density at the upper edge remains unchanged as $\mathit{Pr}$ increases. Moreover, ${\it\delta}$ decreases for higher $\mathit{Pr}$ . Overall the Rayleigh number increases monotonically but not linearly with the Prandtl number, and we can conclude that the indirect impact of $\mathit{Pr}$ on $\mathit{Ra}$ via both ${\rm\Delta}{\it\rho}^{\star }$ and ${\it\delta}$ is of second order while its direct impact on $\mathit{Ra}$ through (3.11) (i.e.  $\mathit{Ra}\propto \mathit{Pr}$ or lower molecular diffusivity is associated with higher Prandtl numbers and thus higher Rayleigh numbers) is of first order.

Figure 5.

Illustrations of the different components of the eigenfunction associated with the eyelid instability. (a) Perturbation streamwise vorticity, ${\it\omega}_{3d}^{\dagger }|_{x}$ , associated with the convectively driven counter-rotating streamwise vorticity component of the mode. (b) Perturbation spanwise vorticity, ${\it\omega}_{3d}^{\dagger }|_{y}$ , associated with the shear-driven counter-rotating spanwise vorticity component of the ‘eyelid’ mode.

Table 3.

The measured values of density contrast across the SAR ( ${\rm\Delta}{\it\rho}^{\star }$ ) and the thickness of the SAR ( ${\it\delta}$ ) as well as the calculated Rayleigh number $\mathit{Ra}$ (3.11) for a vertical profile in the middle of the billow.

A new property of the eyelid instability at high $\mathit{Re}$ which was discovered through our non-separable stability analysis is that the usual low- $\mathit{Re}$ convectively driven instability is transformed into a mixed shear/convection-driven mode. Figure 5 illustrates different characteristics of the eigenfunction associated with this mixed mode. Figure 5(a) depicts the perturbation streamwise vorticity component thus corresponding to the convectively driven aspect of the eyelid instability. This perturbation, if it were to grow in isolation, would have the form of counter-rotating vortex tubes that are periodic in the cross-stream spanwise direction and aligned by the shear. Notice in this figure that the convective component of this instability contains noticeable perturbations that are inclined towards the diagonal edges of the billow in the shear direction, as indicated by a dashed line.

In addition, figure 5(b) illustrates the perturbation spanwise vorticity component of the same eigenfunction, thus corresponding with the shear component of the eyelid instability. If it were to develop in isolation, this perturbation would be manifested as counter-rotating vortices in the streamwise vertical $x{-}z$ plane (notice the positive and negative vorticity regions in the eyelid) wrapped around the core. The required shear for the existence of these counter-rotating spanwise vortices is created due to the KH roll-up process, during which fluid regions characterized by different streamwise velocities become superimposed and produce alternately retrograde and prograde ‘jets’. At high $\mathit{Re}$ the shear becomes sufficiently strong that this component of the vorticity eigenfunction strongly emerges. In fact, the convective and shear components of the eyelid instability are dynamically coupled. The convectively driven streamwise vortex tubes modulate the jet, which is formed by the fluid entering the eyelid, in such a way that the jet becomes alternately compressed (hence strengthened) and expanded (hence weakened). The strengthening and weakening of the jet in turn modulates the shear, which induces counter-rotating spanwise vortices associated with the shear-driven component of the eyelid instability. It is also important to note that the shear-driven component of the eyelid instability was also apparent in the high- $\mathit{Re}$ DNS results of Mashayek & Peltier (Reference Mashayek and Peltier2013); however, in Mashayek & Peltier (Reference Mashayek and Peltier2012a ) (see their figure 16) it was plausibly but incorrectly attributed to an extension of the secondary shear instability on the braid into the eyelid region. We will confirm the emergence of both properties of the eyelid instability in our three-dimensional DNS results, which will be discussed below.

The fluid’s Prandtl number plays a primary role in the growth of the eyelid instability, while its influence on the braid instabilities is only secondary through possible modification of the shear (Smyth Reference Smyth2003) and strain fields around the braid. The indirect impact of $\mathit{Pr}$ on the SPI, which was noted earlier in § 3.2.1, can be recognized by inspection of figure 5(a), where the convectively driven perturbations extend to regions surrounding the stagnation point (shown by dashed circles in this figure), which might modify the shear and strain fields of the braid such that at higher $\mathit{Pr}$ the onset of SPI is postponed. This expected impact will be fully verified in our three-dimensional DNS results, to be discussed in the next section. On the other hand, changing $\mathit{Pr}$ affects the eyelid instability dramatically. We discussed earlier the role of increasing $\mathit{Pr}$ on $\mathit{Ra}$ and hence intensification of the convective forcing of this mode of secondary instability. Furthermore, the shear-driven component of this mode is also promoted as $\mathit{Pr}$ increases (due to the dynamic coupling explained earlier) such that the counter-rotating vortices grow in abundance, appearing first as KH-like vortices on the top and bottom peripheries which propagate along the shear direction and later populate the diagonal corners (in the shear direction) of the eyelid.

Similar to table 2, table 4 includes the growth rate as well as the detailed energetic information associated with the eyelid instability during different times for $\mathit{Pr}=1$ and $\mathit{Pr}=16$ . Notice that irrespective of $\mathit{Pr}$ , the shear contribution ( $\mathscr{S}h$ ) to the growth rate of perturbation kinetic energy is always approximately 70–80 %, while convection ( $\mathscr{H}_{3d}$ ) has a 20–30 % share in the growth rate of the eyelid instability at our choice of $\mathit{Re}$ and $\mathit{Ri}_{0}$ . We also expect the convectively driven component of the eyelid instability to precede the shear-driven component and therefore be manifested earlier, simply because for earlier times during the growth of three-dimensional perturbations (not shown in this table) the relative contribution of convection to shear ( $\mathscr{H}_{3d}/\mathscr{S}h$ ) increases.

Similar to figure 3(a), which showed the persistence of SPI throughout the entire bandwidth of the spanwise wavenumber, figure 6 plots the growth rate of the eyelid instability (second FGM) at $t_{2d}^{s}$ for various spanwise wavenumbers and for $\mathit{Pr}=1,8,16$ . In general due to a higher $\mathit{Re}$ , the eyelid instability is active at shorter wavelengths as noted by Mashayek & Peltier (Reference Mashayek and Peltier2012a ). Nevertheless, as shown in figure 6, by increasing $\mathit{Pr}$ the wavenumber associated with the highest growth rate also increases, which implies a length scale contraction associated with this secondary mode of instability. In particular, note the significant increase of the growth rates of the shortest wavelengths which tends to flatten the growth rate versus cross-stream wavenumber curve. Unlike lower- $\mathit{Re}$ simulations (at $\mathit{Pr}=1$ ) in which there exists a distinct wavenumber at which the growth rate becomes maximum ( $d\approx 3{-}4$ ) (see e.g. Caulfield & Peltier Reference Caulfield and Peltier2000), in our simulations with a much higher $\mathit{Re}$ and $\mathit{Pr}$ , a preferred wavelength ceases to exist as the curve becomes increasingly flat at short wavelengths as $\mathit{Pr}$ increases, suggesting that in the limit of extremely high $\mathit{Pr}$ there would be no wavenumber selectivity at all. This fundamental modification of length scales at higher $\mathit{Pr}$ suggests that the ‘eyelid instability’ is capable of directly injecting turbulence into all scales of motion. We will further discuss this topic in the next sections using three-dimensional DNS results.

Figure 6.

The variation of perturbation growth rate ${\it\sigma}_{3d}$ at $t_{2d}^{s}$ associated with the eyelid instability (second FGM) for $\mathit{Pr}=1,8,16$ for the full bandwidth of spanwise wavenumbers  $d$ .

Table 4.

Similar to table 2 but for the second most unstable three-dimensional mode associated with the eyelid instability.

4.1. The phases of turbulent flow evolution

By employing the velocity decomposition in (2.6) and (2.7), the total kinetic energy budget of the flow can also be decomposed as $\mathscr{K}=\overline{\mathscr{K}}+\mathscr{K}_{2d}+\mathscr{K}_{3d}$ . In figure 7, we plot the time variation of the reservoir of total kinetic energy $\mathscr{K}$ as well as that of two of its individual components associated respectively with the kinetic energy of the spanwise averaged motions $\mathscr{K}_{2d}$ and the kinetic energy of the three-dimensional ‘turbulence’ $\mathscr{K}_{3d}$ . It proves useful at this point to characterize the life cycle of shear layer evolution as involving three primary phases, namely (I) the initial growth phase, (II) the fully turbulent regime, and (III) the final decay regime.

(I) The initial growth phase coincides with the growth of the primary two-dimensional KH instability as well as the secondary (or indeed tertiary) instabilities, catalysed from the initial laminar shear flow with initial minimum gradient Richardson number $\mathit{Ri}_{0}$ below the Miles–Howard critical value of $\mathit{Ri}=0.25$ . For all of the analyses to be discussed in this paper we have fixed the properties of the initial parallel flow to be characterized by $\mathit{Ri}_{0}=0.12$ , a value that is slightly below the Richardson number found to be maximally efficient for mixing at $\mathit{Pr}=1$ in the analyses of Mashayek et al. (Reference Mashayek, Caulfield and Peltier2013). In this initial phase of flow evolution the amplitude of the two-dimensional finite amplitude billow saturates at a time denoted by $t_{2d}^{s}$ in the previous section (the time of the maximum in the $\mathscr{K}_{2d}$ time series). Following the saturation of $\mathscr{K}_{2d}$ , as is evident from figure 7, the growth of intrinsically three-dimensional fluctuations commences, characterized by an abrupt rise in $\mathscr{K}_{3d}$ , leading eventually to complete turbulent collapse of the density stratified shear flow. Such fast growth of $\mathscr{K}_{3d}$ confirms the separation of time scales between the growth rate of primary two-dimensional KH instability and the secondary three-dimensional instabilities, which is an essential condition underlying the validity of our non-separable linear stability analysis as discussed in § 3.1.

A second phase of flow evolution is the fully turbulent regime (II), during which the most energetic small-scale turbulent motions exist and the associated three-dimensional perturbations themselves saturate. The onset of this phase, identified as $t_{f}^{o}$ , can be characterized by the time at which $\mathscr{K}_{3d}$ overwhelms $\mathscr{K}_{2d}$ , which occurs after the full collapse of the KH wave into turbulence. The sustained turbulent phase of flow evolution continues until a time denoted by $t_{f}^{e}$ at which the buoyancy Reynolds number $\mathit{Re}_{b}=|\mathscr{D}|/({\it\nu}\langle N^{2}\rangle _{z})$ (which is a dimensionless measure of turbulent intensity) reaches a critical value, usually taken to be approximately $\mathit{Re}_{b}=20$ (Stillinger, Helland & van Atta Reference Stillinger, Helland and van Atta1983; Smyth & Moum Reference Smyth and Moum2000b ) representing a lower limit below which stratification would be expected to suppress turbulent activity. Our DNS results during this phase correspond to $20<\mathit{Re}_{b}<150{-}250$ , which reaches the respective range of $\mathit{Re}_{b}=O(10{-}10^{4})$ that is commonly associated with oceanographic observations (e.g. Smyth et al. Reference Smyth, Moum and Caldwell2001; van Haren & Gostiaux Reference van Haren and Gostiaux2010).

The final phase of flow evolution is associated with the final decay regime (III), during which the most intense turbulence of the second phase gradually decays, eventually leading to a complete re-laminarization of the flow and the return to a new state of stratified parallel flow, which again might be unstable to small-amplitude perturbations according to the Miles–Howard criterion, provided that further sufficiently strong shear were applied. As a consequence of the mixing that has occurred during the evolution of the turbulent flow, the initial shear layer has thickened and a substantial conversion of initial kinetic energy into potential energy has occurred.

We next focus upon the effects of changing molecular Prandtl number while keeping other non-dimensional control parameters fixed at $\mathit{Re}=6000$ and $\mathit{Ri}_{0}=0.12$ . Figure 7 demonstrates that the two-dimensional KH billow saturates at a time that is essentially independent of $\mathit{Pr}$ , and hence $t_{2d}^{s}$ is almost identical in all cases. In other words, the differences between the molecular diffusion rates for momentum and density are overwhelmed by the faster growth rates of the primary KH instability, and thus $\mathit{Pr}$ has a negligible impact on $\mathscr{K}_{2d}$ during the initial growth phase prior to $t_{2d}^{s}$ . In contrast, insofar as the abrupt growth of the three-dimensional secondary instabilities is concerned, $\mathit{Pr}$ plays a key role due to its influence on emergence, abundance and growth rate of the secondary instabilities that were discussed in depth in § 3.2, and hence the time series associated with the total kinetic energy and its individual components undergo complex variations. We will reassess the impact of $\mathit{Pr}$ on these instabilities and the associated mixing they engender in the following subsections, but here we will first focus on a number of additional characteristics of these phases of flow development.

First, prior to $t_{2d}^{s}$ , two-dimensional KH billows are energized by the background shear flow, but at a time subsequent to $t_{2d}^{s}~(t\approx 75{-}100)$ , $\mathscr{K}_{2d}$ drops, which corresponds to a transfer of energy from the two-dimensional KH wave to the background mean flow and hence to the development of a local maximum in all the time series describing the evolution of $\mathscr{K}$ . In addition, this loss in $\mathscr{K}_{2d}$ also provides the energy for the growth of secondary instabilities (see the abrupt rise in $\mathscr{K}_{3d}$ ). Such extraction of kinetic energy from the two-dimensional KH billow becomes less intense as $\mathit{Pr}$ increases (notice the suppression of the local maximum in $\mathscr{K}$ that occurs in the range $t\approx 90{-}100$ ). The subsequent drop in $\mathscr{K}$ , following the collapse of the original two-dimensional KH billow, is associated with the most intense mixing activity that develops during flow evolution (as will be discussed in detail below) and which results in the transfer of kinetic energy into background potential energy more efficiently than during the sustained turbulent phase.

Second, at low $\mathit{Pr}$ and during the transition to the fully turbulent regime, $\mathscr{K}_{3d}$ has a pronounced double-peak structure, indicating the occurrence of successive periods of gain and loss of turbulent kinetic energy which is associated with concurrent oscillations of $\mathscr{K}_{2d}$ , suggesting that there is a tendency for three-dimensional flow to retain some portion of its initial two-dimensionality immediately after the KH billow collapses, especially at $\mathit{Pr}=1$ and $\mathit{Pr}=2$ (i.e. the turbulent flow is more anisotropic in these cases). However, this tendency appears to be markedly reduced at higher $\mathit{Pr}$ , suggesting that the fully turbulent phase becomes more isotropic as $\mathit{Pr}$ increases. The first of these two successive extrema (more noticeable for the $\mathit{Pr}=1$ and $\mathit{Pr}=2$ cases) occurs near $t\approx 100$ subsequent to the collapse of the KH billow, indicating that the three-dimensional turbulence benefits from such transfer of the APE (that had built up during the KH roll-up) into the background potential energy reservoir. However, such transfer from $\mathscr{K}_{2d}$ to $\mathscr{K}_{3d}$ diminishes at higher $\mathit{Pr}$ , as shown by a marked reduction in $\mathscr{K}_{2d}$ and similarly the reduced gain in $\mathscr{K}_{3d}$ , which is eventually manifested in the disappearance of the first local maximum peak in $\mathscr{K}_{3d}$ for the higher $\mathit{Pr}$ cases. This suggests that the secondary instabilities, which extract their energy for growth from the two-dimensional KH wave, become less ‘effective’ as $\mathit{Pr}$ increases. However, it is interesting to note that the apparent ineffectiveness of the secondary instabilities for higher $\mathit{Pr}$ leads to a more isotropic distribution of energy, one in which $\mathscr{K}_{2d}$ does not regain its energy content (i.e. no rise in $\mathscr{K}_{2d}$ following the initial increase of three-dimensionality) after the initial transfer into $\mathscr{K}_{3d}$ . The fact that an increase in molecular $\mathit{Pr}$ supports a more isotropic turbulent state is in accord with previous findings in the literature (see e.g. Smyth & Moum Reference Smyth and Moum2000a ). We will return to the underlying physical reason behind this in § 4.3 and discuss its connection with the apparent role of $\mathit{Pr}$ in rendering the secondary instabilities less ‘effective’.

A third observation that follows from the results of the analyses documented in figure 7 focuses once more on the properties of the fully turbulent regime, during which the turbulence seems to become increasingly stationary as $\mathit{Pr}$ increases. This is clearly evident on the basis of the almost uniform, time-independent variation of $\mathscr{K}_{3d}$ at high $\mathit{Pr}$ (note especially the results for $\mathit{Pr}=8$ and $\mathit{Pr}=16$ ), during which even the second peak of the former double-peak structure (which is still apparent at $\mathit{Pr}=4$ ) is essentially eliminated. As a consequence of this sustained period of stationary turbulent flow at high $\mathit{Pr}$ , the fully turbulent regime reaches the saturation state much earlier as turbulent evolution of the flow continues, although the maximum intensity of turbulent activity is diminished.

Our final observation in this section focuses on the total loss in kinetic energy as can be observed from the time series of $\mathscr{K}$ . During the course of turbulent evolution from the initial stratified shear flow, the molecular $\mathit{Pr}$ does impact the total loss in the initial kinetic energy and hence is expected to show a similar trend in the increase of the background potential energy, associated with the total entrainment of the ambient fluid into the mixing layer. More precisely, as $\mathit{Pr}$ increases, the total loss of kinetic energy and the associated rise in the background potential energy decrease, implying a less vigorous entrainment of fluid into the mixing layer from the flanks of the shear flow where the ambient speed is highest. We will return to the above preliminary observations after analysing in detail the impacts of the molecular Prandtl number on secondary instabilities and the mixing that ensues. These will be the foci of the following sections.

4.2. Visualizations of secondary instabilities growing on the primary Kelvin–Helmholtz billows

Figure 8 illustrates the three-dimensional structures associated with the important eyelid instability for various values of $\mathit{Pr}$ and at several time instants during the initial growth phase. In this figure, at the time when the two-dimensional flow associated with the main KH wave saturates (i.e. at $t_{2d}^{s}=72{-}74$ ), the convectively driven component of this instability (the so-called shear-aligned convective instability) emerges as streamwise vorticity tubes (e.g. shown for ${\it\omega}_{x}=\pm 0.01$ ) with negative and positive signs. These vortex tubes form clear convective cells which locally overturn the statically unstable fluid that forms in the eyelids of the billows during roll-up. Comparing the structure of these vortex tubes for different $\mathit{Pr}$ at $t_{2d}^{s}$ demonstrates the development of an increasing number of shear-aligned convective rolls in the cross-stream direction as $\mathit{Pr}$ increases and thus a decrease in their characteristic cross-stream wavelength. This is completely in accord with our analytical findings in § 3.2 regarding the effect of higher $\mathit{Pr}$ in shifting the maximum growth rates towards higher wavenumbers (see figure 6) leading to length scale contraction. Furthermore, note that at any spanwise cross-section, the convective rolls possess opposite signs, again in accord with our prediction in figure 5(a).

Figure 8.

Iso-surfaces of ${\it\omega}_{x}=\pm 0.01$ at various $\mathit{Pr}$ showing the growth of shear-aligned secondary convective instability, manifested as vortex tubes with alternating signs peripheral to the billow core.

Figure 9.

Contour plots of spanwise vorticity ${\it\omega}_{y}$ at the mid-point in the cross-stream $y$ direction, demonstrating the growth of both types of braid instability: shear instability on the braid and SPI for various $\mathit{Pr}$ .

At a later time ( $t=80$ ), the other component of the eyelid instability is seen to develop on the initially dominant shear-aligned convective rolls, which is manifested in the form of small-scale spanwise vortices, which were also predicted for this region based upon our linear secondary instability theory of the previous section – a characteristic of this instability that is inherently shear-driven, not convectively driven. These vortices then travel in the shear direction and are eventually obscured by the complete turbulent collapse of the flow, during which evidence also appears of the action of tertiary Rayleigh–Taylor instability, which can be recognized at $t=82$ in the region surrounding the eyelids, for example. Meanwhile, as the secondary convective instability perturbs the eyelid region, the SPI dominates the action in the central portion of the braid with a spanwise wavenumber similar to that of the eyelid instability (also shown in figure 8).

The growth of the secondary shear-driven instabilities (either on the braid or the eyelid region of the billow) can be better captured in a streamwise vertical ( $x{-}z$ plane) illustration of the spanwise vorticity. Figure 9, illustrates contours of spanwise vorticity at the mid-point in the cross-stream $y$ direction, during the same phase of the KH cycle, a time period after the saturation of the KH billow and before its full collapse into turbulence. Notice in this figure that the stagnant core associated with the SPI mode reduces in size, implying its suppression at higher $\mathit{Pr}$ in agreement with our analytical findings in the previous section. Moreover, the spanwise structure associated with this unstable mode resembles that of the shear-aligned convective rolls (as shown in figure 8); this may be understood after recalling the contribution of $\mathscr{H}_{3d}$ to the energetics of the SPI mode as noted in table 2. As the sequence of plots at different $\mathit{Pr}$ in figure 9 shows, secondary shear instabilities of the braid vorticity emerge after both shear-aligned convective instability and SPI and are excited by perturbations induced by the stagnant vortex core or extensions of the convective instability near the braid (see for example figure 9, $t=82$ at $\mathit{Pr}=1$ ). Nevertheless, comparing the number of such vortex cores along the braid at $t=82$ for all $\mathit{Pr}$ , suggests that the induced shear instability on the braid tends to be promoted as $\mathit{Pr}$ increases. This can be explained by recalling the direct influence of $\mathit{Pr}$ in sharpening the density gradients across the braid (Smyth Reference Smyth2003) as well as the indirect influence of $\mathit{Pr}$ through the suppression of SPI at higher $\mathit{Pr}$ , which has to do with the weakening of the braid strain field (see figure 3). As a consequence, braid shear dominates, leading to the emergence of secondary shear instability on the braid. Note that the eigenfunction associated with the braid instability already supports the potential growth of these perturbations (albeit weakly), as shown in figure 2.

A critical observation from figure 9 involves identification of the shear-driven component of the eyelid instability. Notice in this figure the counter-rotating shear instabilities that first appear on the top and bottom of the eyelid and then grow both in size and abundance as they travel in the shear direction. This instability leads to a core periphery that is most intensely populated with small-scale coherent structures along the top-right and bottom-left diagonal of the billow (when the initial parallel shear flow is to the right (left) in the upper (lower) region), consistent with the eigenfunction associated with the spanwise vorticity of the eyelid instability (see figure 5 b). In addition, as $\mathit{Pr}$ increases, the billow becomes populated even more densely with these shear-driven instabilities.

As far as the pairing instability is concerned, we have performed similar simulations (to be discussed in the next section) but with double the streamwise wavelength at $\mathit{Pr}=1,2,4$ (at $\mathit{Re}=6000$ and $\mathit{Ri}_{0}=0.12$ with $2L_{x}$ ) in order to allow the emergence of the pairing instability. Nevertheless, none of these new three-dimensional DNS experiments showed any visually identifiable signature of the pairing mode, verifying our initial conjecture that a value of $\mathit{Re}=6000$ would continue to be large enough, even at $\mathit{Pr}>1$ , to rule out the possibility of an upscale cascade of energy to a wavelength longer than that which is occurring at the injection scale. It is also shown in the next section that the mode associated with the pairing instability is in fact unchanged or further suppressed with an increase in $\mathit{Pr}$ .

The sequence of plots in figures 8 and 9 also demonstrates that, in all cases, the KH billow includes a rather ‘quiet’ core, one which is less amenable to the loss of its structural integrity and hence more difficult to mix, as opposed to the intense activity that occurs in both the eyelids and the braid. In other words, while an increase in $\mathit{Pr}$ leads to a proliferation of coherent small-scale features in both the eyelid and on the braid, the core remains the last portion of the KH wave anatomy to be engulfed in the turbulent collapse process in all of the cases investigated here. We believe that this is also a characteristic of any such high- $\mathit{Re}$ flow, and is plausibly relevant to the understanding of the results reported for the oceans in Geyer et al. (Reference Geyer, Lavery, Scully and Trowbridge2010). Based on their observations of KH billows in the high-Reynolds-number flow of a salt-stratified estuary, characterized by extremely high $\mathit{Pr}$ , Geyer et al. (Reference Geyer, Lavery, Scully and Trowbridge2010) inferred ‘quiet conditions within the cores’, on the basis of the weak acoustic backscattered signals corresponding with the in situ measurements of low conductivity variance within the KH cores observed (as opposed to the intense signals and high-conductivity variance they identified to be originating from the braids). Their picture is completely consistent with our interpretation of figure 9 concerning the quietness of the core region as a focus of irreversible mixing. In order to further clarify the abundance of secondary instabilities on the KH braid and in the eyelid(s), as opposed to its core at relatively high $\mathit{Re}$ and high $\mathit{Pr}$ , figure 10 illustrates contours of total viscous dissipation ( ${\it\varepsilon}$ in (2.14)) at several instants of time during the transition period, for different $\mathit{Pr}$ , and provides a similar picture to those of the eco-sounding images usually obtained from oceanographic observations. Similar contours for the scalar variance field ${\it\chi}=|{\rm\nabla}{\it\rho}|^{2}$ (which might be more suitable for comparison with their in situ measurements of conductivity variance or density structure) were also investigated, which provided visually indistinguishable results (i.e. no density variance within the core together with intensified density variance on its periphery and the braid) compared to figure 10, and thus are not presented here. These plots further illustrate the lack of viscous dissipation (i.e. turbulent activity) as well as the lack of density variance in the core region, while all the activity occurs on its periphery and on the braid. Although high $\mathit{Re}$ primarily controls such quietness of the core (as shown in figure 10 for $\mathit{Pr}=1$ ), increasing $\mathit{Pr}$ further reduces the diffusion inside the core and localizes the billow instabilities to smaller scales around the periphery with marginal interaction with the central core; thus it might be expected that this influence would also reduce the core activity.

Figure 10.

Contour plots at the mid-point in the cross-stream $y$ direction, associated with viscous dissipation ${\it\varepsilon}$ as defined in (2.14), demonstrating a ‘quiet core’ with almost no turbulent activity (i.e. zero dissipation associated with white regions) during the growth of three-dimensional perturbations and before being engulfed by the surrounding turbulence.

Our findings here concerning the quietness of KH cores for high- $\mathit{Re}$ flows does not support one of the earlier conjectures of Mashayek & Peltier (Reference Mashayek and Peltier2011), in which it was plausibly conjectured that the absence of turbulent activity in the billow cores could be due to an early homogenization of the core by the rapid emergence and vigour of secondary instabilities in this region at very high $\mathit{Re}$ and high $\mathit{Pr}$ . As shown in figure 10, the core region, once formed, stays almost inactive until it is eventually engulfed by the surrounding turbulent activity (at $t\simeq 110$ for example); thus a ‘quiet core’ is a direct consequence of high- $\mathit{Re}$ flow (as hypothesized by Geyer et al. Reference Geyer, Lavery, Scully and Trowbridge2010), which leads to an almost passive core, and is not a consequence of the growth of instabilities internal to the core leading to its early homogenization. This must also be seen to follow from our non-separable linear stability analyses, which have never revealed the existence of core-centred modes of secondary instability.

Figure 11.

The evolution of PSD (4.3) associated with (a) ${\it\xi}_{1}$ pairing mode and (b) ${\it\xi}_{2}$ KH mode.

4.3. Spectral representation

Before concluding this subsection, it will prove useful to attempt to confirm our analytical findings in § 3 and gain further quantitative insight concerning the impacts of molecular Prandtl number on both length scales and growth rates of the secondary instabilities which facilitate the transition of the stratified shear flow into turbulence.

First, insofar as the pairing instability is concerned, employing the additional simulations at $\mathit{Pr}=1,2,4$ (with fixed $\mathit{Ri}_{0}=0.12$ and $\mathit{Re}=6000$ ) in which the computational domain includes two wavelengths of the primary KH wave in the streamwise direction (i.e. $2L_{x}$ ), we have spectrally decomposed the spanwise averaged density field ${\it\rho}_{2d}(x,z,t)$ in the streamwise direction, for all the times throughout flow evolution such that ${\it\rho}_{2d}$ can be represented as

(4.1) $$\begin{eqnarray}{\it\rho}_{2d}(x,z,t)=\mathop{\sum }_{m=-M}^{M}c_{m}({\it\xi}_{m},z,t)\text{e}^{\text{i}{\it\xi}_{m}x},\end{eqnarray}$$

in which $c_{m}({\it\xi}_{m},z,t)$ is the complex Fourier coefficient at a discrete wavelength ${\it\xi}_{m}=2{\rm\pi}m/(2L_{x})$ with streamwise wavenumber $m$ . Note that with this definition for ${\it\xi}_{m}$ , the first mode ( $m=1$ ) represents the pairing mode while the second mode ( $m=2$ ) represents the primary KH wave. Furthermore, $c_{m}$ can be expressed as

(4.2) $$\begin{eqnarray}c_{m}({\it\xi}_{m},z,t)=\langle {\it\rho}_{2d}\text{e}^{-\text{i}{\it\xi}_{m}x}\rangle _{x}.\end{eqnarray}$$

Figure 11 illustrates the variation in the power spectral density (PSD), defined as

(4.3) $$\begin{eqnarray}P({\it\xi}_{m},t)=\langle |c_{m}({\it\xi}_{m},z,t)|^{2}\rangle _{z}\end{eqnarray}$$

for the first and second modes corresponding respectively with the pairing and original KH structures. The first mode (when normalized by the amplitude of the second) associated with the pairing instability becomes evident only during the fully turbulent phase in figure 11(a), and interestingly its maximum PSD decreases or remains unchanged with $\mathit{Pr}$ , indicating that an increase in $\mathit{Pr}$ still suppresses the pairing mode of instability. In the laboratory experiments, $\mathit{Re}$ is low but also the working fluid is usually salt-stratified water, which is characterized by a high value of $\mathit{Pr}$ . These factors, low $\mathit{Re}$ and high $\mathit{Pr}$ , seem to support the occurrence of an inverse cascade of energy through merging of adjacent KH billows, because $\mathit{Re}$ is the primary factor inhibiting this mode of instability as shown by Mashayek & Peltier (Reference Mashayek and Peltier2013). If $\mathit{Re}$ is sufficiently large (which is the case here), the pairing instability is not sufficiently strong to survive the competition with other modes of instability that are all contributing to the transition of the flow into turbulence. As a result, the amalgamation of KH billows, which is the manifestation of the inverse cascade of energy from the injection scale to longer wavelengths, does not occur. However, energy at twice the wavelength of the FGM would always be weakly evident in the fully turbulent flow, simply due to the imposed boundary conditions.

We next focus on the other modes of instability (braid and eyelid instability) to further investigate the impact of molecular Prandtl number on their growth rates and length scales. For this purpose, at every time step during the transition period we have spectrally decomposed the streamwise vorticity field ${\it\omega}_{x}=\partial w/\partial y-\partial v/\partial z$ , using the fast Fourier transform. This spectral representation is assumed to take the form

(4.4) $$\begin{eqnarray}{\it\omega}_{x}(x,y,z,t)=\mathop{\sum }_{n=-N}^{N}c_{n}({\it\gamma}_{n},x,z,t)\text{e}^{\text{i}{\it\gamma}_{n}y},\end{eqnarray}$$

in which $c_{n}({\it\gamma}_{n},x,z,t)$ is the complex Fourier coefficient at a discrete wavelength ${\it\gamma}_{n}=2{\rm\pi}n/L_{y}$ with wavenumber $n$ , and can be expressed as

(4.5) $$\begin{eqnarray}c_{n}({\it\gamma}_{n},x,z,t)=\langle {\it\omega}_{x}\text{e}^{-\text{i}{\it\gamma}_{n}y}\rangle _{y}.\end{eqnarray}$$

In order to avoid aliasing errors due to under-sampling the spanwise direction where we have approximately $KN_{y}$ data points (see table 1), we set the maximum wavenumber $N$ for our calculation to be below the Nyquist frequency. Of particular interest is the relative significance of various components of the spectrum, as represented by the PSD $P({\it\gamma}_{n})$ defined as

(4.6) $$\begin{eqnarray}P({\it\gamma}_{n},t)=\langle |c_{n}({\it\gamma}_{n},x,z,t)|^{2}\rangle _{xz}.\end{eqnarray}$$

We also compartmentalize the flow domain into (a) the core and its periphery, (b) the core region only, (c) the braid including the stagnation point, and (d) the stagnation point only, in order to most accurately estimate the characteristics of the instabilities that occupy these different regions of flow evolution (for an illustration of these compartments refer to figure 1 in Caulfield & Peltier Reference Caulfield and Peltier2000). In order to compare the results obtained at different $\mathit{Pr}$ , we normalize the above function by the total spectral power of the entire computational domain at $\mathit{Pr}=1$ , (i.e.  $P_{tot}|_{Pr=1}=\sum _{n=0}^{N}\langle |c_{n}({\it\gamma}_{n},x,z,t)|^{2}\rangle _{xz}$ ). Figures 12 and 13 compare the power spectrum associated with the first 25 wavenumbers at $t_{2d}^{s}$ and $t=82$ respectively, for each of these compartments, as well as the range of $\mathit{Pr}$ investigated in this work. The panels in each figure correspond to compartments (a)–(d) discussed above. Also note that as we have verified, replacing ${\it\omega}_{x}$ with ${\it\omega}_{y}$ (which might seem more relevant to the stagnation point and secondary shear instability on the braid) in the above formalism, does not change the conclusions to be discussed below.

Figure 12.

Wavenumber ( ${\it\gamma}$ ) dependence of the normalized PSD, $\hat{P}({\it\gamma})$ , of streamwise vorticity ${\it\omega}_{x}$ (4.6) at time $t=t_{2d}^{s}$ (refer to the text for normalization details). Different panels are associated with four regions enclosing different parts of KH anatomy, namely (a) eyelid, (b) core, (c) braid and (d) stagnation point. Refer to the legend in panel (a) for different values of $\mathit{Pr}$ .

Figure 13.

Similar to figure 12 but at time $t=82$ .

At $t_{2d}^{s}$ and for the eyelid (figure 12 a), there is a clear shift towards higher wavenumbers for the FGM associated with the highest PSD. This agrees completely with our earlier analytical predictions in figure 6 and also observations of numerical results in figure 8, in which more rolls of convective instability appeared in the spanwise direction as $\mathit{Pr}$ increased. Thus the eyelid instability must be characterized by a shorter wavelength with increasing $\mathit{Pr}$ . Not only does the wavenumber increase at higher $\mathit{Pr}$ , but the power contained in the whole spectrum is also enhanced as $\mathit{Pr}$ increases. This is an indication of a more vigorous eyelid instability in agreement with our discussion in § 3.2.2. The core region shows a similar trend but contains about two orders of magnitude less power at the same wavenumbers, highlighting again the fact that most of the activity occurs on its periphery (see figure 10). At this time during the course of transition to turbulence, the braid and the region enclosing the stagnation point (figures 12 c and 12 d respectively) are characterized by an almost identical spectrum, which is a manifestation of the dominance of SPI over the shear instability on the braid and is in accord with our observation of no KH-like shear instability on the braid at $t_{2d}^{s}$ in all cases illustrated in figure 9. However, at the higher end of the wavenumber spectrum, an increase in $\mathit{Pr}$ directly promotes more energetic motions at smaller scales on the braid, a feature that was also evident in figure 12(a) for the eyelid, as discussed above. As we showed in figure 5(a), the perturbation eigenfunction associated with the convectively driven characteristic of the eyelid instability extends to regions near the braid and is thus responsible for respective modification of the PSD at higher wavenumbers, similar to the eyelid region. Having related the PSD at $t_{2d}^{s}$ to the stagnation point only, figure 12(d) demonstrates that higher $\mathit{Pr}$ has minimal or even a suppressing influence on the PSD associated with this instability, which agrees with our discussion of figure 3 and the fact that the stagnation point vortex diminishes in size as $\mathit{Pr}$ increases, as shown in figure 9.

Figures 13(a)–13(d) focus on a later time (e.g.  $t=82$ ), at which secondary shear instability of the braid has had enough time to emerge. The first two panels further reinforce the points made earlier concerning the role of increasing $\mathit{Pr}$ on the growth of secondary convective instability, as well as the striking shift towards higher wavenumber for the FGM, and the fact that the core remains ‘quiet’ compared to its periphery. In addition, and as a consequence of increasing molecular $\mathit{Pr}$ and the associated changes in the secondary instabilities, the distinction between the power spectrum associated with the braid region and that of a zone enclosing the stagnation point has become evident at this time which demonstrates (i) the delayed dominance of secondary shear instability on the braid after the earlier dominance of the SPI and (ii) promotion of secondary KH-like shear instability on the braid, as was suggested on the basis of the flow visualizations of figure 9.

6.1. Mixing efficiency

The instantaneous efficiency of mixing for the DNS analyses we have performed is calculated using (2.22) and is plotted in figure 14(a) for all values of the Prandtl number investigated herein. In analysing the behaviour of mixing efficiency, we have to give separate treatment to different phases of the flow as defined in § 4.1. Insofar as the influence of $\mathit{Pr}$ during the initial period is concerned, the irreversible process of molecular diffusion of density inside the shear layer diminishes as $\mathit{Pr}$ increases and therefore the instantaneous mixing efficiency decreases (i.e. the growth of $\mathscr{P}_{B}$ is delayed). During this initial phase the viscous dissipation is minimal, leading to high mixing efficiency, which is irrelevant to the turbulent regime. Furthermore, prior to the complete collapse of the KH billows (i.e. before $t\approx 85{-}90$ ) when the maximum efficiency occurs, the impact of $\mathit{Pr}$ on mixing efficiency continues to have a clear trend such that as $\mathit{Pr}$ increases the efficiency of mixing decreases. The high efficiency during this period is entirely due to the absence of intense viscous dissipation, which is accompanied by a rapid transfer of APE $\mathscr{P}_{A}$ into background potential energy $\mathscr{P}_{B}$ (see figure 14 b). Furthermore, the decrease in the efficiency maximum is associated with the increased viscous dissipation associated with higher $\mathit{Pr}$ .

Subsequent to the emergence of secondary instabilities and the corresponding growth of inherently three-dimensional perturbations which eventually leads to turbulent KH collapse, an abrupt drop in efficiency, independent of $\mathit{Pr}$ , occurs. This can be explained by a dramatic drop in the magnitude of APE, $\mathscr{P}_{A}$ (see e.g. figure 14 b). The main mechanism responsible for this drop in $\mathscr{E}_{i}$ or APE is driven by the growth of secondary instabilities which leads to a gradual destruction of the KH billow. From an energetics perspective, this mechanism is associated with the transfer of kinetic energy from $\mathscr{K}_{2d}$ , that had been stored in the form of APE within the KH billows, to $\mathscr{K}_{3d}$ (see figure 7 and the associated discussion in § 4.1). It is important to note that this drop in APE and the corresponding rise in $\mathscr{P}_{B}$ should be attributed to the growth of all secondary instabilities, which are themselves growing during this period and which drive the increase in $\mathscr{P}_{B}$ . Simultaneously within this period, the viscous dissipation is also intensified due to the increased small-scale activity. The molecular diffusion, once more, becomes a dominant mechanism for controlling mixing and its efficiency during the decay phase until the flow re-laminarizes. Consequently, during the phase of turbulence decay, higher $\mathit{Pr}$ results consistently in a reduced efficiency of mixing.

By the time three-dimensional turbulence is sufficiently energetic (i.e. during the onset of fully developed flow at $t_{f}^{o}\approx 125{-}150$ ), the mixing efficiency appears to be uninfluenced by $\mathit{Pr}$ . The instantaneous mixing efficiency undergoes complicated but moderate shifts about the canonical value of 0.2, and it is unclear how secondary instabilities could conceivably impact irreversible mixing during this phase. Figure 14(c) illustrates the $\mathit{Pr}$ -dependence associated with the mean value of $\mathscr{E}_{i}$ averaged over the fully turbulent phase (precisely during $t\in [t_{3d}^{s}~t_{f}^{e}]$ , where $t_{3d}^{s}$ denotes the saturation point of three-dimensional kinetic energy) and the error bars represent one root-mean-square deviation from the mean values. It seems that the mean value for instantaneous efficiency of mixing decreases monotonically for higher Prandtl numbers. However, given the uncertainties associated with this averaging (indicated by error bars), the statistical significance ascribed to this trend is questionable. To further investigate this, we instead employ the cumulative mixing efficiency defined by Caulfield & Peltier (Reference Caulfield and Peltier2000), which is consistent with the definition of mixing efficiency employed in experiments (e.g. Linden Reference Linden1979) based on energy amounts rather than their rates of change (cf. (2.22)), i.e.

(6.1) $$\begin{eqnarray}\mathscr{E}_{c}=\frac{\langle \mathscr{M}\rangle _{f}}{\langle \mathscr{M}-\mathscr{D}\rangle _{f}},\end{eqnarray}$$

where $\langle F(t)\rangle _{f}=\int _{t_{3d}^{s}}^{t_{f}^{e}}F(t)\,\text{d}t$ represents the time integral over the period of sustained turbulence (i.e. $t_{3d}^{s}\leqslant t\leqslant t_{f}^{o}$ ).

Figure 14(d) plots the $\mathit{Pr}$ -dependence associated with $\mathscr{E}_{c}$ which avoids (by definition) the uncertainties involved in the former approach. As this figure shows, the efficiency of mixing does in fact decrease monotonically for the range of $\mathit{Pr}$ investigated herein. This interesting result may be construed to support the hypothesis of Mashayek & Peltier (Reference Mashayek and Peltier2013) that the route to turbulence (i.e. the detailed nature of the ‘zoo’ of secondary instabilities that populate the transition phase of flow evolution) matters insofar as mixing efficiency is concerned (Mashayek & Peltier Reference Mashayek and Peltier2013). We can argue that with an increase in $\mathit{Pr}$ , due to the contraction of energy-containing length scales, the turbulent eddies responsible for enhancing the background potential energy by extracting it from the APE (or equivalently the turbulent eddies that energize APE by extracting it from the kinetic energy reservoir) become less effective. Alternatively, due to an increased level of isotropy associated with the turbulent flow as $\mathit{Pr}$ increases, the strain field associated with turbulent eddies is weakened and thus it supports weaker scalar gradients, which in turn hinders turbulent mixing (Batchelor Reference Batchelor1959). As a consequence of this ineffectiveness associated with the energy transfer into the background potential energy, the impact of molecular $\mathit{Pr}$ on the irreversible mixing efficiency becomes one in which $\mathscr{E}_{c}$ decreases slightly at higher $\mathit{Pr}$ .

The moderate decrease in the mixing efficiency associated with the fully turbulent flow in our high- $\mathit{Pr}$ , high- $\mathit{Re}$ DNS analyses must be distinguished from the relatively substantial drop in $\mathscr{E}_{i}$ that was shown to be the case by Smyth et al. (Reference Smyth, Moum and Caldwell2001) for their low- $\mathit{Re}$ simulations. This new finding points once more to the importance of high $\mathit{Re}$ , which is the geophysically relevant regime in which the mixing efficiency remains essentially fixed to the near vicinity of the canonical value of 0.2 (as shown by the dashed lines in figure 14 d) for the initial minimum Richardson number of $\mathit{Ri}_{0}=0.12$ investigated here. Also, the upward concavity associated with the variations of $\mathscr{E}_{c}$ at higher $\mathit{Pr}$ suggests that mixing efficiency could asymptote to values closer to 0.1 for much higher $\mathit{Pr}$ , assuming that our $\mathit{Re}$ has reached its asymptotic value above which the mixing efficiency does not change markedly. Besides $\mathit{Pr}$ , the other non-dimensional number that dominantly controls the efficiency of mixing during this phase is $\mathit{Ri}_{0}$ , as is already evident in the analyses of Caulfield & Peltier (Reference Caulfield and Peltier2000) and most recently investigated by Mashayek et al. (Reference Mashayek, Caulfield and Peltier2013), in which they have established the existence of a non-monotonic dependence of $\mathscr{E}_{c}$ on $\mathit{Ri}_{0}$ having a maximum efficiency near $\mathit{Ri}_{0}=0.16$ . However, the question of whether this non-monotonicity shifts towards higher/lower $\mathit{Ri}$ or remains unchanged at higher $\mathit{Pr}$ (which is more relevant to the ocean) remains an open question.

Note that the nominal values of $\mathit{Re}$ and $\mathit{Ri}_{0}$ are simply the initial values of the investigated simulations (as listed in table 1) and do not characterize the fully developed turbulent flows generated. An alternative parameter that may be employed to characterize the stratified turbulent flows is the buoyancy Reynolds number, $\mathit{Re}_{b}$ as introduced in § 4.1. Figure 15 illustrates the time variations of $\mathit{Re}_{b}$ for all of the values of $\mathit{Pr}$ that have been investigated, showing that they all reach values as high as $\mathit{Re}_{b}\approx 150{-}250$ at the onset of fully developed flow, as pointed out earlier in § 4.1. Insofar as the impact of $\mathit{Pr}$ on the mixing efficiency is concerned, in this paper we have developed a physical understanding as to why mixing efficiency decreases at higher $\mathit{Pr}$ . Ascribing such behaviour to $\mathit{Re}_{b}$ or $\mathit{Re}_{{\it\lambda}}$ (see the appendix A for a discussion of $\mathit{Re}_{{\it\lambda}}$ ) or other variations of parameters would simply be to rephrase the question as to why these additional parameters change in the way they do with $\mathit{Pr}$ . Furthermore, since mixing efficiency varies non-monotonically with $\mathit{Re}_{b}$ , it is not a trivial matter as to how one might simply characterize the decrease in mixing efficiency with $\mathit{Re}_{b}$ .

Figure 15.

Time variations of the buoyancy Reynolds number $\mathit{Re}_{b}=$ $|\mathscr{D}|/({\it\nu}\langle N^{2}\rangle _{z})$ for all of the three-dimensional DNS analyses with different values of $\mathit{Pr}$ .

Figure 16.

(a) The time series associated with background potential energy $\mathscr{P}_{B}$ for various $\mathit{Pr}$ , demonstrating a decrease in entrainment of ambient flow (reduced rise in $\mathscr{P}_{B}$ ). (b) Vertical profile of density $\overline{{\it\rho}}=\langle {\it\rho}\rangle _{xy}$ showing a reduction in stratified shear layer expansion with increasing $\mathit{Pr}$ after the flow is re-laminarized.

6.2. The effect of $\mathit{Pr}$ on entrainment

Figure 16(a) illustrates the time variation of the background potential energy, $\mathscr{P}_{B}$ , a notional state of static stability which has been obtained after an adiabatic rearrangement of the density field into a monotonically decreasing upward assembly following the method depicted in figure 1. The total rise in $\mathscr{P}_{B}$ is directly linked to the total entrainment of the ambient flow on the flanks of the shear layer into the central mixing region. As figure 16(a) demonstrates, such a rise is suppressed as $\mathit{Pr}$ increases, which is in accord with what we had predicted on the basis of our observation from figure 7 regarding the total loss in kinetic energy of the flow. Another indication of the occurrence of less effective entrainment at higher $\mathit{Pr}$ can be seen on the basis of the ultimate stratified shear layer thickness after the flow is re-laminarized. Figure 16(b) illustrates the density profiles associated with the re-laminarized shear layer for the range of $\mathit{Pr}$ we have investigated. This plot shows that the eventual expansion of the stratified shear layer is reduced following the end of the turbulent life cycle in fluids with higher molecular Prandtl number, which also suggests that at higher $\mathit{Pr}$ , the processes involved during the course of turbulence evolution would become less effective insofar as entrainment is concerned.

Here it might be puzzling why different rates of density and momentum diffusion would physically lead to such variation in the re-laminarized shear layer thickness, or the total gain in the background potential energy, given the impact of $\mathit{Pr}$ on promoting the strength of the eyelid instability (see § 3), as the primary catalyst for driving the irreversible mixing process. A plausible argument resolving this conundrum can be stated in two equivalent ways. (I) Increasing $\mathit{Pr}$ not only strengthens the major secondary instabilities but also reduces the associated energy-containing length scales (for example, see the reduced horizontal scale of the shear-aligned convective rolls in figure 8), which overall reduces their effectiveness in entraining fast moving fluid on the flanks of the shear layer into the well-mixed region near the original inflection point of the shear profile. Alternatively (II), due to the presence of the ultraviolet catastrophe at higher $\mathit{Pr}$ (actually plausible only in the limit of extremely high $\mathit{Pr}$ ), kinetic energy of the flow is directly injected into the dissipation subrange (i.e. the Batchelor scale) and hence a quasi-steady turbulent phase occurs. However, at low $\mathit{Pr}$ , interactions between filaments of considerably larger length scales are precursors to the full downscale cascade of energy. As a consequence, shear layers with lower $\mathit{Pr}$ contain significantly larger coherent structures than their higher- $\mathit{Pr}$ counterparts during the same phase and the entrainment of fast-moving fluid from the flanks of the shear layer becomes more efficient, thereby leading to the increase of the background potential energy and the associated entrainment.