2 Formulation

We begin by considering the following version of the non-dimensionalised, monochromatic body-forced, incompressible, Boussinesq equations

where we define the (external) Reynolds number $Re$ , the bulk stratification parameter $B$ and the Prandtl number $Pr$ as

Here, $\boldsymbol{u}(x,y,z,t)=u\hat{\boldsymbol{x}}+v\hat{\boldsymbol{y}}+w\hat{\boldsymbol{z}}$ is the three-dimensional velocity field, $p$ is the pressure and the density is decomposed into $\unicode[STIX]{x1D70C}_{tot}=\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D70C}_{B}(z)+\unicode[STIX]{x1D70C}(\boldsymbol{x},t)$ , i.e as the sum of a Boussinesq reference density, a constant linear background stratification and a fully varying disturbance density. We have non-dimensionalised using the characteristic length scale $L_{y}/2\unicode[STIX]{x03C0}$ , characteristic time scale $\sqrt{L_{y}/2\unicode[STIX]{x03C0}\unicode[STIX]{x1D706}}$ and density gradient scale $\unicode[STIX]{x1D6FD}=\text{d}\unicode[STIX]{x1D70C}_{b}/\text{d}z$ . (see LCK17 for details). Furthermore, $n$ is the forcing wavenumber, $\unicode[STIX]{x1D706}$ is the forcing amplitude, $\unicode[STIX]{x1D708}$ is the kinematic viscosity, and $\unicode[STIX]{x1D705}$ is the molecular diffusivity. We impose periodic boundary conditions in all directions and solve over the cuboid $[0,2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}]\times [0,2\unicode[STIX]{x03C0}]^{2}$ , where $\unicode[STIX]{x1D6FC}=L_{y}/L_{x}$ defines the horizontal aspect ratio of the domain. Vorticity $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$ is used as the prognostic variable and DNS are performed using the fully dealiased (two-thirds rule) pseudospectral method with mixed fourth-order Runge–Kutta and Crank–Nicolson timestepping implemented in CUDA to run on GPU cards. We initialise the velocity field’s Fourier components with uniform amplitudes and randomised phases in the range $2.5\leqslant |\boldsymbol{k}|\leqslant 9.5$ such that the total enstrophy $\langle |\unicode[STIX]{x1D74E}|^{2}\rangle _{V}=1$ and $\unicode[STIX]{x1D70C}=0$ initially. Throughout, $\langle (\cdot )\rangle _{V}:=\unicode[STIX]{x1D6FC}\iiint (\cdot )\,\text{d}x\,\text{d}y\,\text{d}z/(2\unicode[STIX]{x03C0})^{3}$ denotes a volume average, $\langle (\cdot )\rangle _{h}:=\unicode[STIX]{x1D6FC}\iint (\cdot )\,\text{d}x\,\text{d}y/(2\unicode[STIX]{x03C0})^{2}$ denotes a horizontal average and $\langle (\cdot )\rangle _{v}:=\int (\cdot )\,\text{d}z/(2\unicode[STIX]{x03C0})$ denotes a vertical average.

We define the diagnostics involved in the energetic budgets as

where $\text{d}{\mathcal{K}}/\text{d}t={\mathcal{I}}-{\mathcal{B}}-{\mathcal{D}}$ , $\text{d}{\mathcal{P}}/\text{d}t={\mathcal{B}}-\unicode[STIX]{x1D712}$ and ${\mathcal{K}}$ is the total kinetic energy density, ${\mathcal{P}}$ the density variance, ${\mathcal{I}}$ is the energy input by the forcing, ${\mathcal{B}}$ is the buoyancy flux, ${\mathcal{D}}$ is the viscous dissipation rate, and $\unicode[STIX]{x1D712}$ is the density variance dissipation rate. We fix $\unicode[STIX]{x1D6FC}=0.5$ to avoid subcritical transition, $Pr=1$ for numerical efficiency and $n=1$ , i.e. the flow is forced with $\sin (y)\hat{\boldsymbol{x}}$ , to mimic closely other large-scale shear profiles previously studied (e.g. stratified Taylor–Couette flow (Woods et al. Reference Woods, Caulfield, Landel and Kuesters2010) and vertically sheared Kolmogorov flow (Garaud, Gallet & Bischoff Reference Garaud, Gallet and Bischoff2015)) and denote time averages with overbars, i.e. $\overline{(\cdot )}=[\int _{0}^{T}(\cdot )\text{d}t]/T$ , where $T$ is normally the full simulation time after discarding the first five time units to allow for flow spin-up.

Table 1.

Imposed parameters and diagnostic outputs for the three groups of simulations. The non-dimensional Kolmogorov microscale is given by $\unicode[STIX]{x1D702}=Re^{-3/4}\overline{{\mathcal{D}}}^{-1/4}$ and is scaled by the maximum wavenumber allowable $k_{max}=N/3$ . Groups A and B have resolution $128\times 64^{2}$ and group C has $256\times 128^{2}$ to retain spectral convergence. The non-dimensional Ozmidov length scale $l_{O}=\overline{{\mathcal{D}}}^{1/2}/B^{3/4}$ , the average buoyancy Reynolds number $\overline{Re}_{B}=\overline{{\mathcal{D}}}Re/B$ , the average bulk Richardson number $\overline{Ri}_{B}$ and the Taylor microscale Reynolds number $Re_{\unicode[STIX]{x1D706}}={\mathcal{K}}_{turb}\sqrt{10Re/\overline{{\mathcal{D}}}}$ are also listed. Here we define ${\mathcal{K}}_{turb}=(1/2)\langle (\boldsymbol{u}-\overline{\boldsymbol{u}})^{2}\rangle _{V}$ and overbars are always time averages over the full $T=1000$ window after discarding the first five time units for spin-up.

3 Irreversible mixing in the direct numerical simulations

We begin by characterising the mixing in this system by conducting a set of simulations across a range of $Re$ and $B$ . The pertinent single-point diagnostics are presented in table 1. In particular we examine the irreversible mixing efficiency which, following Salehipour & Peltier (Reference Salehipour and Peltier2015) and Maffioli, Brethouwer & Lindborg (Reference Maffioli, Brethouwer and Lindborg2016), we define in terms of the dissipation rates:

The dependence of $\mathscr{E}$ on typical flow parameters is a crucial ongoing problem in stratified turbulence, and a source of some controversy. As discussed by Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014), the flow of a turbulent, shear-driven stratified fluid has several competing time scales, and it is natural to attempt to parameterise $\mathscr{E}$ in terms of parameters quantifying the relative importance of these time scales. Using shear-instability simulations and comparison with observations, Salehipour et al. (Reference Salehipour, Peltier, Whalen and MacKinnon2016) developed a parameterisation using two parameters: an appropriately defined ‘buoyancy Reynolds number’, effectively a ratio of the time scale of the stratification to the time scale of the turbulence; and an appropriately defined Richardson number, effectively a ratio of the square of the time scale of the vertical shear (assumed to be the dominant driver of the mixing) to the time scale of the stratification. In order to determine how the mixing efficiency varies with these dimensionless numbers in this flow, we define the buoyancy Reynolds number $Re_{B}(t)$ and the (local) gradient Richardson number $Ri_{G}(\boldsymbol{x},t)$ as (with the usual caveat that comparing different studies is difficult if the key parameters are defined differently)

where $(\unicode[STIX]{x2202}u_{h}/\unicode[STIX]{x2202}z)^{2}=(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z)^{2}+(\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z)^{2}$ . $Ri_{G}$ is a pointwise quantity which characterises the relative stability of the flow to overturning shear instabilities, and calculating an appropriate overall average for this quantity to characterise a particular flow should be performed with care, not least because locations of low vertical shear significantly skew the distributions of $Ri_{G}$ . The appropriate averaging used here is thus

There are two important points to appreciate about these definitions. First, we have chosen to define $Ri_{B}$ with the total gradients of density while $Re_{B}$ only depends explicitly on the background $N_{B}.$ Second, although in principle these are independent non-dimensional parameters, it remains to be established that they are independent in this flow, as there exist situations where they are strongly correlated (see e.g. Zhou, Taylor & Caulfield Reference Zhou, Taylor and Caulfield2017).

Figure 1.

The variation of average mixing efficiency $\overline{\mathscr{E}}$ with: (a) average bulk Richardson number $\overline{Ri}_{B}$ ; and (b) average buoyancy Reynolds number $\overline{Re}_{B}$ . The horizontal line is $\overline{\mathscr{E}}=1/6$ , consistent with the upper bound of Osborn (Reference Osborn1980). The closed loops mark the projections onto the time-dependent $(Ri_{B},\mathscr{E})$ and $(Re_{B},\mathscr{E})$ planes of the two unstable relative periodic orbits UPO-o1 (solid line) and UPO-l1 (dashed line) discussed in § 4. (c) The variation of $\overline{Re}_{B}$ with $\overline{Ri}_{B}$ . Different symbols mark different values of $Re$ . The red curve in (b) corresponds to $\overline{\mathscr{E}}=\overline{Re}_{B}^{-1/2}$ , the large $Re_{B}$ scaling behaviour suggested by Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005).

In figure 1(a,b) we plot the time-averaged mixing efficiency $\overline{\mathscr{E}}$ against the time-averaged parameters $\overline{Ri}_{B}$ and $\overline{Re}_{B}$ . We observe non-monotonic dependence of $\overline{\mathscr{E}}$ on both parameters, in at least qualitative agreement with the experimental analysis of Linden (Reference Linden1979) and numerical simulations of Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005). Indeed, $\mathscr{E}$ appears to saturate near the critical value of $1/6$ (marked with a horizontal line), equivalent to the upper bound on the turbulent flux coefficient $\unicode[STIX]{x1D6E4}\simeq \mathscr{E}/(1-\mathscr{E})\leqslant 0.2$ proposed by Osborn (Reference Osborn1980), and then decreases for $\overline{Re}_{B}\gtrsim 30$ , with a dependence not entirely unlike the $\overline{\mathscr{E}}\propto \overline{Re}_{B}^{-1/2}$ observed by Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) in what they referred to as the ‘energetic’ regime (see Ivey, Winters & Koseff (Reference Ivey, Winters and Koseff2008) for a more detailed discussion).

However, although these similarities are intriguing, there are still significant differences. Firstly, we do not observe a ‘transitional’ plateau of approximately constant mixing efficiency at intermediate $\overline{Re}_{B}$ . Indeed, the qualitative structure of the mixing efficiency curve is more reminiscent of the Padé approximant proposed by Mashayek et al. (Reference Mashayek, Salehipour, Bouffard, Caulfield, Ferrari, Nikurashin, Peltier and Smyth2017) using (vertical) shear-instability numerical data as well as observations, where $\unicode[STIX]{x1D6E4}\propto Re_{B}^{1/2}$ and $Re_{B}^{-1/2}$ for small and large $Re_{B}$ respectively. Second, and once again reminiscent of the approach proposed by Mashayek et al. (Reference Mashayek, Salehipour, Bouffard, Caulfield, Ferrari, Nikurashin, Peltier and Smyth2017) to parameterise mixing in terms of $\overline{Re}_{B}$ alone, the two parameters $\overline{Ri}_{B}$ and $\overline{Re}_{B}$ are actually closely correlated, and our simulations actually trace out a (monotonic) curve in $\overline{Ri}_{B}{-}\overline{Re}_{B}$ space, as shown in figure 1(c), independently of the externally imposed parameters $B$ and $Re$ . The low mixing efficiency observed at high $\overline{Ri}_{B}$ can thus be understood as being entirely associated with weaker turbulence (smaller $\overline{Re}_{B}$ ) and the maximum mixing efficiency occurs at a sweet spot of sufficiently vigorous turbulence at sufficiently high stratification, analogously to the results of Zhou et al. (Reference Zhou, Taylor and Caulfield2017). We stress that we are not claiming that this clear correlation between $\overline{Ri}_{B}$ and $\overline{Re}_{B}$ is generic (see for example the discussion in Scotti & White (Reference Scotti and White2016)), just that it occurs for this flow.

It is apparent that at small $\overline{Ri}_{B}$ , $\overline{Ri}_{B}\propto \overline{Re}_{B}^{-1}$ . This is unsurprising, as the turbulence and shear are largely unaffected by such weak stratification, and the variation of both parameters with $B$ completely dominates. However, as $\overline{Ri}_{B}$ increases beyond the value associated with the most efficient mixing, the turbulence and shear do indeed become suppressed by the strengthening stratification, and the power-law dependence steepens so $\overline{Ri}_{B}\propto \overline{Re}_{B}^{-3/2}$ . There is undoubtedly also an increasing significance of viscosity, which becomes dominant for the largest values of $\overline{Ri}_{B}>1$ , where $\overline{Re}_{B}$ drops below one, and the turbulence is essentially completely suppressed, with very inefficient mixing. However, it is always important to remember that both $\overline{Ri}_{B}$ and $\overline{Re}_{B}$ are global measures averaged in both space and time. Both $Ri_{B}$ and $Re_{B}$ are functions of time, and there could also be spatial variation of $Ri_{G}$ within the flow domain. Such variation proves to be a key part of the behaviour of the periodic orbits which we identify, and hence of the mixing occurring in these flows.

4 Recurrent flow analysis

In order to determine what processes underpin the turbulence in our simulations, we use ‘recurrent flow analysis’. The general approach is as follows. First a simulation is conducted during which near recurrences are located in the chaotic trajectory. This is achieved by storing a historical record of state vectors and periodically checking the history against the current state vector. When an appropriately ‘near repeat’ is identified, it is stored for a later convergence attempt with a Newton–GMRES–hookstep algorithm. For more details on this approach see Chandler & Kerswell (Reference Chandler and Kerswell2013), Lucas & Kerswell (Reference Lucas and Kerswell2015). Here, we are interested in finding representative unstable periodic orbits (UPOs) with qualitatively different mixing properties in a broad range of flows, rather than identifying a large number of UPOs for a specific set of parameters. Therefore, we have identified five orbits in two broad classes: one class of two orbits found in the relatively weakly stratified simulation A1 with $Re=20$ and $B=1$ , and the other class of three orbits in the more strongly stratified simulation B3 with $Re=50$ and $B=50$ . The first class has $\overline{Re}_{B}\sim O(10)$ , and we label it as class ‘o’, for ‘overturning’, while the second class has $\overline{Re}_{B}\sim O(1)$ , and we label it as class ‘l’ for layered. We plot projections onto the time-dependent $(Ri_{B},\mathscr{E})$ and $(Re_{B},\mathscr{E})$ planes of the trajectories of characteristic relative periodic orbits UPO-o1 (solid line) and UPO-l1 (dashed line) on figures 1(a) and 1(b), demonstrating that these two classes do indeed have qualitatively different mixing properties, and also that these time-dependent properties are consistent with the time-averaged properties of the ‘full’ DNS.

Table 2.

Table cataloguing the unstable recurrent flows showing the DNS from which they come (table 1): period, $T_{p}$ ; relative shifts in $x$ and $z$ directions, $s_{x}$ , $s_{z}$ ; discrete shift/reflect in $y$ , $m_{y}$ ; inverse stability coefficient, $\unicode[STIX]{x1D6EC}^{-1}$ (see Chandler & Kerswell (Reference Chandler and Kerswell2013)); the number of unstable directions, $N_{\unicode[STIX]{x1D706}}$ ; buoyancy Reynolds number, $\overline{Re}_{B}$ ; maximum and time-averaged mixing efficiency, $\mathscr{E}$ ; and time-averaged bulk Richardson number, $\overline{Ri}_{B}$ .

Table 2 lists properties of the converged orbits, including period, relative shifts due to the continuous symmetries in $x$ and $z$ , stability and some diagnostic averages. By comparison with the averages of the DNS in table 1, it is clear the UPOs are reproducing the bulk time averages quite well. The projection of the trajectories onto the $(Re_{B},{\mathcal{E}})$ plane of the five converged orbits are plotted in figure 2, with greyscale colours representing the probability density function (p.d.f.) of the turbulent DNS. Notice that the darker colours represent regions where the turbulent trajectory spends more time, and lighter colours represent less frequent excursions to that part of phase space. For simulation B3, the members of the ‘l’ class of converged UPOs sit across the darkest region of the p.d.f. where the turbulence spends most of its time but miss the higher-efficiency, intermittent turbulent bursts that the DNS exhibits. Missing extreme events is a known failing of the recurrent flow analysis in some circumstances and is a topic of ongoing research, see Lucas & Kerswell (Reference Lucas and Kerswell2015). However, for the ‘o’ class in simulation A1, the orbit UPO-o1 (in particular) appears to span the turbulent attractor well in this projection, missing only some very rare excursions to high $Re_{B}$ . As shown in table 2, the ‘o’ class orbits are considerably more unstable than the orbits in the ‘l’ class.

In order to compare how the time scales of the UPOs relate to those exhibited by the turbulence, we plot in figure 3 the power spectra of the kinetic energy for the two DNS signals from which the UPOs have been extracted (i.e. A1 and B3). Vertical lines show that the periods of the orbits are approximately coincident with observed frequencies. For simulation A1 the largest non-zero peak in the frequency spectrum occurs at approximately half the fundamental frequency observed in the UPOs, i.e. $f\approx 2\unicode[STIX]{x03C0}/15=0.42$ , with the fundamental period around $7.5$ , UPO-o1 having just under three such periods. For simulation B3 the periods of the orbits are on the edge of the main distribution of frequencies near $f=0.6$ . This is reasonable evidence that the time scales of the UPOs are representative of the typical time scales of the turbulence. We can also consider the buoyancy frequency and its influence on these cases. For simulation A1 where $B=1$ , the background buoyancy frequency $N_{B}=\sqrt{B}=1$ so the time scale or ‘buoyancy period’ is $T_{N}=2\unicode[STIX]{x03C0}$ , which is not far from the ‘fundamental period’ we have established from the UPOs and spectrum. For B3, the buoyancy period $T_{N}=2\unicode[STIX]{x03C0}/\sqrt{50}\approx 0.9$ which is much faster than the time scales found in the spectrum. We find some signature of this frequency later in this section. Henceforth, we will focus attention on UPO-o1 and UPO-l1 since, having larger amplitude, they represent better proxies for the turbulence, although the dynamics within each class is broadly similar.

Figure 2.

Projection onto the ( $Re_{B},\mathscr{E}$ ) plane of the UPOs and the p.d.f. of the direct numerical simulations rendered in greyscale, with darker shades denoting that the turbulent trajectory spends more time there for: (a) class ‘o’ UPOs from simulation A1; and (b) class ‘l’ UPOs from simulation B3.

Figure 3.

Power spectra for the kinetic energy $\hat{{\mathcal{K}}}(f)=\int _{T}{\mathcal{K}}(t)\text{e}^{\text{i}ft}\text{d}t$ for the DNS signals from simulation A1 (a) and simulation B3 (b) with vertical lines showing approximately the periods of the UPOs extracted.

Figure 4.

Time evolution of $Ri_{B}$ and $\mathscr{E}$ for: (a) UPO-o1; and (b) UPO-l1. Notice the correlation between $Ri_{B}$ and $\mathscr{E}$ for (a) and not for (b). The times highlighted in figures 5 and 6 are marked with stars.

Figure 5.

Three-dimensional rendering of: (a) $\unicode[STIX]{x1D70C}_{tot}=\unicode[STIX]{x1D70C}-z$ ; (b) $u$ ; and (c) p.d.f.s of $Ri_{G}$ in the entire computational domain at $t=15,\,16.5$ and $20$ for UPO-o1. The p.d.f.s are computed with linear intervals of $Ri_{G}$ , while the inset shows the result when distributing the intervals on a log scale. Vertical lines mark $Ri_{G}=1/4$ .

Figure 6.

Three-dimensional rendering of: (a) $\unicode[STIX]{x1D70C}_{tot}=\unicode[STIX]{x1D70C}-z$ ; (b) $u$ ; and (c) p.d.f.s of $Ri_{G}$ in the entire computational domain at $t=4,\,5.5$ and $7$ for UPO-l1. The p.d.f.s are computed with linear intervals of $Ri_{G}$ , while the inset shows the result when distributing the intervals on a log scale. Vertical lines mark $Ri_{G}=1/4$ .

Further evidence that the two classes are qualitatively different is shown in figure 4, where we plot the time dependence of $Ri_{B}$ and $\mathscr{E}$ for UPO-o1 and UPO-l1. The dynamics of UPO-o1 show relatively low $Ri_{B}$ , strongly correlated to the mixing efficiency $\mathscr{E}$ . There are three distinct subperiods, each of which shows the bulk $Ri_{B}$ decrease to small values (indeed below the canonical value of $1/4$ for the linear instability of a stratified (vertical) shear layer given by the Miles–Howard theorem and frequently invoked as a diagnostic for overturning), followed by a maximum in mixing efficiency as the flow overturns. Such overturnings are visualised in figure 5, where snapshots of total density $\unicode[STIX]{x1D70C}_{tot}$ and streamwise velocity $u$ are chosen before and after the final minimum of $Ri_{B}$ and near the subsequent maximum of $Ri_{B}$ following the peak $\mathscr{E}$ as marked on figure 4(a). Supplementary movies of the time evolution of these fields (movies 1 and 2) are available at https://doi.org/10.1017/jfm.2017.695. The density field shows distinctive Kelvin–Helmholtz-like billows forming at $t=16.5$ where $\mathscr{E}$ and $Ri_{B}$ are beginning their growth phase, justifying the labelling of this UPO as being in class ‘o’. Of course, a bulk measure of the Richardson number does not capture the stability properties of the flow, and so in figure 5 we also plot the p.d.f. of $Ri_{G}(\boldsymbol{x})$ over the spatial domain at the same times, with a vertical line indicating $Ri_{G}=1/4$ . These distributions show a striking increase of the proportion of the domain with $Ri_{G}<1/4$ prior to the overturning mixing event. This proportion then decreases with time, but does not vanish completely. Even in the periods of maximum bulk $Ri_{B}$ some regions of the domain remain with $Ri_{G}<1/4$ . Apparently, it is necessary for an appreciable proportion of the domain to have a local gradient Richardson number $Ri_{G}<1/4$ before overturning is observed.

By contrast, as is apparent in figure 4(b), UPO-l1 has much larger overall $Ri_{B}$ , seemingly uncorrelated with the behaviour of the markedly lower mixing efficiency $\mathscr{E}$ . This is also apparent in the snapshots shown in figure 6 at the characteristic times marked on figure 4(b) (and in supplementary movies 3, 4 and 5). No overturns occur, and the cycle of mixing behaviour is associated with straining or scouring motions drawing the perturbation density $\unicode[STIX]{x1D70C}$ into thin layers, justifying the labelling as class ‘l’. The instantaneous p.d.f.s of the gradient $Ri_{G}$ plotted in figure 6 now show that nowhere in the domain has $Ri_{G}<1/4$ , even near $t=5.5$ when mixing efficiency is maximised. At peak $\mathscr{E}$ , $t\approx 5.5$ , the snapshot of $\unicode[STIX]{x1D70C}$ shows increased layers and arguably larger localised magnitudes of density gradient.

In order to characterise the straining/scouring motions exhibited by the ‘l’ class, we examine the properties of the strain tensor

This tensor has three real eigenvalues, which we denote by $\tilde{\unicode[STIX]{x1D6FC}},\,\tilde{\unicode[STIX]{x1D6FD}},\,\tilde{\unicode[STIX]{x1D6FE}}$ , referred to as the principal strains, which sum to zero in an incompressible flow. In general $\tilde{\unicode[STIX]{x1D6FC}}>0$ and is therefore the stretching strain, $\tilde{\unicode[STIX]{x1D6FE}}<0$ and is the compressional strain and the intermediate strain $\tilde{\unicode[STIX]{x1D6FD}}$ will control the three-dimensionality, i.e. $\tilde{\unicode[STIX]{x1D6FD}}<0$ represents compression in two directions and fluid elements drawn into filaments, $\tilde{\unicode[STIX]{x1D6FD}}>0$ represents compression in one direction and fluid elements are drawn into sheets and $\tilde{\unicode[STIX]{x1D6FD}}=0$ represents a purely two-dimensional straining field. Figure 7 shows time series of the volume-averaged strains for the two characteristic UPOs, UPO-o1 and UPO-l1. In the overturning case, UPO-o1 shows that $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}$ is well correlated with the mixing efficiency, being positive and maximum when $\mathscr{E}$ is largest. This is in good agreement with Smyth (Reference Smyth1999), who examined these strains for a freely decaying Kelvin–Helmholtz unstable shear layer; $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}$ is positive and has a distinct maximum during the second phase of the overturn where the flow becomes more isotropic (isotropic turbulence is observed to have the strains approximately in the ratio 3 : 1 :  $-4$ (Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987; Smyth Reference Smyth1999)) and approaches zero in the quiescent periods as the flow relaxes back to the parallel shear mean flow. For the layered case, UPO-l1 is distinctly more anisotropic; the overall absolute values for $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}$ are an order of magnitude smaller and are actually negative for much of the cycle. At the point where the mixing efficiency is maximal ( $t\approx 5.5$ ) $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}$ also has a maximum, but here this amounts to crossing the axes and $\langle \tilde{\unicode[STIX]{x1D6FC}}\rangle _{V}\approx -\langle \tilde{\unicode[STIX]{x1D6FE}}\rangle _{V}.$ Our interpretation is then that the flow approaches two-dimensionality where the mixing is largest and increased (in magnitude) vertical gradients of $\unicode[STIX]{x1D70C}$ are formed by the change in sign of $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}.$ This is in stark contrast to the overturning case, where mixing has a more isotropic straining signature with large $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}.$ As a secondary observation, the higher-frequency oscillations of $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}$ for UPO-l1 are of the order of $N=\sqrt{B}$ the buoyancy frequency, i.e. the buoyancy period here is $T_{N}=2\unicode[STIX]{x03C0}/\sqrt{50}\approx 0.9.$ This suggests that internal waves are also playing a role in this dynamical cycle.

Figure 7.

Time series of the volume-averaged principal strains for UPO-o1 (a) and UPO-l1 (b). For the overturning case (a), peaks of $\mathscr{E}$ are coincident with large positive intermediate strain $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}$ , nearer to the values of isotropic turbulence, where as the more layered case (b) is very anisotropic with small $\langle \tilde{\unicode[STIX]{x1D6FD}}\rangle _{V}$ , approximately zero when $\mathscr{E}$ is largest around $t\simeq 5.5$ .

5 Discussion and conclusions

In this paper we have successfully applied recurrent flow analysis to stratified flows for the first time. As expected, success is restricted to relatively weak turbulence at modest Reynolds numbers. Nevertheless, the approach has still provided detailed insight into the sustaining mechanisms and mixing processes at work in these flows. In particular we have established that mixing can occur from two distinct mechanisms; overturning and scouring, consistently with the classification presented by Woods et al. (Reference Woods, Caulfield, Landel and Kuesters2010). Scouring is observed when stratification is strong and overturning is observed when the stratification is weak enough to allow small gradient Richardson numbers. By examining the cycles in space and time we have demonstrated that the spatial distribution of local gradient $Ri_{G}$ is markedly skewed before high overturning mixing, but remains bounded from below by $1/4$ when scouring mixing dominates.

From this analysis it is clear that characterising mixing by appropriate measures of Richardson number is important from both a statistical and dynamical point of view, although we have also shown that weak mixing can be controlled by processes uncorrelated to $Ri_{G}$ . Crucially, weaker average mixing may not be simply controlled by spatiotemporally intermittent shear instability, but by other scouring, straining processes, which may require other predictive diagnostics. Furthermore, since the bulk measures of Richardson number and buoyancy Reynolds number are functionally related, the results of our motivational DNS for horizontally forced, sustained turbulence suggest that an $(\overline{Ri}_{B},\overline{Re}_{B})$ parameterisation of mixing should be treated with caution, particularly since it proved difficult to access a strongly stratified, yet strongly turbulent regime, analogously to the situation arising in stratified plane Couette flow (Zhou et al. Reference Zhou, Taylor and Caulfield2017).

One open question is the interplay between scouring and overturning. Clearly overturning shear instability will overwhelm any background straining or scouring motions; however, it remains a challenge of rationalising spatiotemporal chaos to predict the switching between these processes. For instance, simulation B3 shows that intermittent bursts can raise the local (in time) mixing efficiency, despite the majority of the mixing in this case being controlled by the straining and scouring layered motions characteristic of class ‘l’. Examination of a segment of the trajectory reaching large $\mathscr{E}$ in this case suggests straining and shear instability, saturating at various finite amplitudes, combine to populate the distribution of figure 1(b). Furthermore, as already widely discussed, research into exact coherent structures embedded in turbulent flows faces a serious challenge in the extension of methods to handle higher Reynolds numbers and spatially localised dynamics. It remains unclear whether such UPOs can be identified in ‘energetic’ stratified turbulence with $\overline{Re}_{B}\gtrsim O(100)$ , yet it is of particular interest to understand whether the observed reduction of mixing efficiency at such high $\overline{Re}_{B}$ is associated with a qualitative change in mixing properties.