3.1.1 Power-law scaling and the virtual origins

We begin by examining the decay of non-stratified homogeneous, isotropic turbulence for the three cases used as the initial conditions for the stratified cases, and refer the reader to George (Reference George1992), Skrbek & Stalp (Reference Skrbek and Stalp2000) and Davidson (Reference Davidson2010), and the references therein, for theoretical and experimental studies of the decay of homogeneous, isotropic turbulence. A primary interest here is in the power-law decay relationships

where $u^{\prime }$ is the volume-averaged r.m.s. value of any of the velocity components, $t$ is measured from a virtual origin $t_{0}$ and $L^{\prime }$ is the corresponding integral length scale, either longitudinal or transverse. If $L^{\prime }$ is proportional to ${u^{\prime }}^{3}/\unicode[STIX]{x1D716}$ , where $\unicode[STIX]{x1D716}$ is the volume-averaged dissipation rate of turbulence kinetic energy, then $m=1-n$ , although with laboratory or numerical simulation data for homogeneous, isotropic turbulence this relationship may not be observed exactly and many studies report only $n$ (cf. de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2000). From the plots by Perot & de Bruyn Kops (Reference Perot and de Bruyn Kops2006) of $n$ as a function of Reynolds number based on data from numerous sources, and from Saffman’s theory (Saffman Reference Saffman1967; Davidson Reference Davidson2010), we expect $n\approx 0.6$ for the current simulations.

The interdependence of $n$ and the virtual origin $t_{0}$ makes the calculation of these parameters and their corresponding uncertainties dependent on how the curve fitting is carried out, which is beyond the scope of this paper. Here our immediate objective is to establish that the non-stratified simulations exhibit power-law decay with a value of $n$ consistent with that in the literature. Figure 2 indicates that this is the case. In the plots $u^{\prime }$ , $v^{\prime }$ and $w^{\prime }$ denote the r.m.s. velocities in the $x$ , $y$ and $z$ directions, respectively, with $z$ in the vertical direction. In the figure, the symbols representing the non-stratified flow data lie on curves corresponding to $n=0.56$ , 0.56 and 0.55 for Reynolds numbers 160 (Case I), 600 (Case II) and 2325 (Case III), respectively. These are consistent with the expectation that $n\approx 0.6$ .

Figure 2.

Root-mean-square velocities versus time $t$ . Symbols are for non-stratified cases (only for Case I in d). The reference slopes are determined by least-squares fits to the logarithms of the data in the adjacent lines or symbols.

Figure 3.

The longitudinal and transverse integral lengths of $u$ versus time. The longitudinal scales $L_{ux}$ are plotted true and the others are offset vertically in increments of half-decades. The symbols are for non-stratified cases. The dashed-line reference slopes are $1-n$ from the non-stratified simulations. The solid-line reference slopes are as derived in Davidson (Reference Davidson2010).

3.1.2 Power-law scaling in the stratified flows

Also plotted in figure 2 are the r.m.s. velocities and in figure 3 the integral length scales for the three stratified cases, the latter which are computed as recommended in appendix E of Comte-Bellot & Corrsin (Reference Comte-Bellot and Corrsin1971). The notation $L_{ux}$ denotes the integral length scale of $u$ in the $x$ -direction, defined in terms of the autocorrelation of $u$ , and similarly for the other length scales.

Consider first the r.m.s. velocities, which are plotted in figure 2. There is an adjustment immediately following the onset of stratification at $t=1$ , the abruptness of which can be observed in the evolution of the r.m.s. density in figure 4. By approximately $t=3$ ( $T=1$ , i.e. one buoyancy period), the decay rates for $u^{\prime }$ and $v^{\prime }$ begin to deviate from their non-stratified values, and settle into power-law decay rates which are less than the non-stratified values; they are: $n=0.31$ for Case I, and $n=0.25$ for Cases II and III. From these results, and from figure 2(d), it is seen that the decay rates are very weakly dependent on the Reynolds number, with the fastest decay rate observed for Case I, the case with the lowest Reynolds number. This observed decrease in the decay rate due to stable stratification is consistent with the prediction of Davidson (Reference Davidson2010), although the decay rates are slightly less (in absolute value) than his estimate of $t^{-0.4}$ . The reduction in decay rate is also qualitatively consistent with the experiments of Lin & Pao (Reference Lin and Pao1979), Stillinger et al. (Reference Stillinger, Helland and Van Atta1983), Yoon & Warhaft (Reference Yoon and Warhaft1990) and Spedding (Reference Spedding1997), and with the observations of Almalkie & de Bruyn Kops (Reference Almalkie and de Bruyn Kops2012b ) who, in simulations forced to maintain constant energy in a narrow band of small wavenumbers, observed that stronger stratification slightly reduces the dissipation rate. These results are inconsistent with the laboratory experiments of Liu (Reference Liu1995), however, who finds a decay rate of the horizontal velocities of approximately $n=0.825$ with stratification, faster than the rate of 0.65 for his non-stratified case. They are also inconsistent with the results of Praud et al. (Reference Praud, Fincham and Sommeria2005) who found the decay rate of approximately $n=0.65$ for all their non-stratified and stratified cases. Note that, in our simulations with stable stratification, after some initial adjustment, $w^{\prime }$ decays at almost the same rate as with no stratification for cases II and III. When the Reynolds number is low (Case I at later times) however, $w^{\prime }$ decays faster with stratification than without. This is most easily observed in figure 2(a).

Figure 4.

Root-mean-square of the fluctuating density $\unicode[STIX]{x1D70C}^{\prime }$ versus time $t$ . The dashed reference slopes are from the total r.m.s. velocities in the stratified cases.

In figure 3 it is seen that the growth rates of all the horizontal integral length scales are enhanced by the density stratification. The growths of the vertical integral scales of $u^{\prime }$ , however, are so strongly affected by the stratification as to become negative at later times for cases Case II and Case III. Both of these results are qualitatively consistent with the predictions of Davidson (Reference Davidson2010). The result for the vertical scales is consistent with the suggestion by Lilly (Reference Lilly1983) that, as the effects of stratification become important, the vertical shearing of the horizontal velocities increases. It is also consistent with the scaling analysis of Billant & Chomaz (Reference Billant and Chomaz2001). The fact that all but the vertical lengths are larger in the stratified cases than in the non-stratified cases does not necessarily mean that coherent structures are larger, but rather that the small scales are suppressed. This will become more apparent when we discuss energy spectra.

The r.m.s. of the density fluctuations provides yet another measure of the decay rate of turbulence with stable density stratification. For a passive scalar $\unicode[STIX]{x1D719}$ in isotropic turbulence with $u^{\prime }\propto t^{-n}$ , similarity analysis and data indicate that $\unicode[STIX]{x1D719}^{\prime }\propto t^{-n}$ as well (cf. de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2000). Since the stratified flows exhibit different decay rates for the horizontal and vertical velocities, we might not expect the decay rate of the scalar fluctuations to agree with that of any single velocity component. From figure 4, however, it is observed that for cases II and III the decay rate of the density fluctuations is almost the same as that of the total r.m.s. velocity. Interpreting ${\unicode[STIX]{x1D70C}^{\prime }}^{2}$ in terms of potential energy (see § 3.1.3), this implies that the potential energy decays at almost the same rate as the horizontal components of the kinetic energy for these cases. Significant Reynolds number effects are observed in the figure. For the case with lowest Reynolds number, the decay rate of the scalar after $T=1$ is faster than for the total r.m.s. velocity. At later times in this case (Case I), the decay rate is even faster than for the non-stratified simulation, consistent with it following decay rate of $w^{\prime }$ . Recall, for example, that Liu (Reference Liu1995) observed that stratification increased the decay rate of velocity, especially the vertical velocity, and so Case I at later times may provide some insight into the dynamics of those laboratory experiments. This issue will be addressed further in § 3.2.

Noticeable in the plots of $w$ and $\unicode[STIX]{x1D70C}$ is oscillatory behaviour, indicating a wave-like behaviour for these flows. The period of oscillation is approximately $0.4\,T_{B}$ for Case I, and appears to lengthen slightly with Reynolds number. The early time dynamics of these flows is strongly nonlinear, however, so that one would not expect internal waves with distinct dispersion relationships to exist at early times. Note that the Froude number $Fr_{h}$ can be considered as an estimate of the nonlinear term compared to the time derivative term in the momentum equation. Therefore, as the flows evolve and $Fr_{h}$ drops below order 1, one would expect internal waves to exist in the flows, but strongly modified by other motions in the flows. It is possible that the reduced decay rates in the flows is partially due to the appearance of internal wave-like motions.

3.1.3 Flow energetics

Information regarding the energy decay rates is given in figure 5. Here $E_{h}$ is the kinetic energy associated with $u$ and $v$ , and $E_{v}$ is that associated with $w$ . We define $E_{p}$ , the potential energy of the density fluctuations relative to the ambient density, as

$E_{p}$ can be considered as a surrogate for available potential energy (Lorenz Reference Lorenz1955). In panel (a) of the figure, all the energies have been normalised by the local (in time) total energy, i.e. by $E_{h}+E_{v}+E_{p}$ . From this panel it is observed that the horizontal components of the kinetic energy become increasingly dominant over the sum of the vertical component and the potential energy. This could be due to one or more of several effects. Possibly internal waves are being absorbed by the more horizontal motions, as in the laboratory experiments of Godoy-Diana et al. (Reference Godoy-Diana, Chomaz and Billant2004), leading to a decrease in both $E_{v}$ and $E_{p}$ while enhancing $E_{h}$ . Furthermore, as kinetic energy is exchanged with potential energy through the action of the buoyancy flux $-(2\unicode[STIX]{x03C0}/Fr)^{2}\langle \unicode[STIX]{x1D70C}w\rangle$ , with $\langle \cdot \rangle$ denoting a spatial average over the entire domain, the dissipation rate of potential energy $\unicode[STIX]{x1D712}$ (see § 3.1.4) provides another avenue for dissipation in addition to the kinetic energy dissipation rate $\unicode[STIX]{x1D716}$ . From this panel, it is also observed that the partition between the contributions to total energy depends on the Reynolds number, since Reynolds number is the only parameter that is varied between the simulations. This Reynolds number dependence is not surprising. For example, Perot & de Bruyn Kops (Reference Perot and de Bruyn Kops2006) reported, from examining many sources, that the decay exponents for non-stratified turbulence are larger in magnitude for lower Reynolds numbers. Finally, note from panel (b), and the fact that the non-stratified energy decay rate is approximately $t^{-1.1}$ , that the decay rate of the total energy is inhibited by the stratification, consistent with the inhibition of the decay rates of the horizontal components of the velocity.

Figure 5.

(a) Ratio of energy contributions to total energy versus time. (b) Total energy versus time. In places where it is difficult to see the curves for Case II, they almost coincide with those for Case III. In (b), the two reference slopes correspond to $n=-0.55$ and $n=-0.28$ in figure 2(c).

The behaviour of the buoyancy flux for all three stratified cases is seen in figure 6. The oscillations about 0 in the buoyancy flux in figure 6(a) represent net exchanges from the vertical component of the kinetic energy to potential energy (positive values) and back from potential into kinetic energy (negative values), consistent with the plots of these components of the energy, and corresponding to the noted oscillations in $w^{\prime }$ and $\unicode[STIX]{x1D70C}^{\prime }$ . These oscillatory results are qualitatively consistent with the numerical simulations of Riley et al. (Reference Riley, Metcalfe, Weissman and West1981) and the laboratory experiments of Liu (Reference Liu1995). Figure 6(b) is a plot of the cumulative (integrated) value of the buoyancy flux, divided by the initial total energy, giving the net transfer of kinetic into potential energy. It is seen that, for cases II and III, the cumulative flux asymptotes to approximately 30 % of the initial total energy. Also note that the cumulative flux is somewhat smaller for Case I, the lowest Reynolds number case. These facts are closely related to the behaviour of the mixing efficiency $\unicode[STIX]{x1D702}$ discussed in the next section.

Figure 6.

(a) Buoyancy flux, normalised by $w^{\prime }\unicode[STIX]{x1D70C}^{\prime }$ , versus time for all three stratified cases. (b) Cumulative buoyancy flux normalised by initial total energy.

3.1.4 Dissipation rates and mixing efficiency

The dissipation rates of potential and kinetic energies are

and

respectively, where $\unicode[STIX]{x1D64E}$ is the non-dimensional strain rate tensor. The time evolutions of their volume averages are shown in figure 7. Noticeable is the decrease in $\unicode[STIX]{x1D716}$ due to the stratification. Furthermore, there is a slight decrease and then increase in $\unicode[STIX]{x1D716}$ as the Reynolds number is increased, whereas there is a slight increase in $\unicode[STIX]{x1D712}$ as it is increased. Also evident is that the non-stratified cases are not quite in power-law decay at $t=1$ , which is also apparent in figure 2 for all the non-stratified cases; however, power-law decay does occur in the non-stratified simulations well before $T=1$ . In addition note that, after the flow adjustment, the ratio of $\unicode[STIX]{x1D712}$ to $\unicode[STIX]{x1D716}$ is approximately 0.54 except for the lowest Reynolds number case (Case I), where the ratio continues to decrease in time. We will find below in § 3.2 that, for this lowest Reynolds number case, the flow might not be fully turbulent beyond approximately 2 buoyancy periods.

Figure 7.

Values of $\unicode[STIX]{x1D716}$ (upper curves), $\unicode[STIX]{x1D712}$ (lower curves). Symbols are for Case III with no stratification.

Of special interest is the ratio $\unicode[STIX]{x1D702}\equiv \unicode[STIX]{x1D712}/(\unicode[STIX]{x1D712}+\unicode[STIX]{x1D716})$ , which is sometimes taken to be the definition of the mixing efficiency (Osborn Reference Osborn1980; Smyth, Moum & Caldwell Reference Smyth, Moum and Caldwell2001). From figure 8 it is seen that, for the higher Reynolds number cases, $\unicode[STIX]{x1D702}$ asymptotes to a value of approximately 0.35, compared to the value of 0.3 obtained by Riley & de Bruyn Kops (Reference Riley and de Bruyn Kops2003) in a similar regime, to 0.25 obtained by Maffioli, Brethouwer & Lindborg (Reference Maffioli, Brethouwer and Lindborg2016) at very low Froude number, to 0.41 inferred from the data of Liu (Reference Liu1995), and 0.17 assumed by Osborn (Reference Osborn1980) based upon ocean data. In the forced homogeneous simulations reported by Almalkie & de Bruyn Kops (Reference Almalkie and de Bruyn Kops2012b ), $\unicode[STIX]{x1D702}$ is approximately 0.36 for $Fr_{h}=$ between 0.258 to 1.47. Note that the value for $\unicode[STIX]{x1D702}$ of 0.35 is consistent with the asymptotic ratio of $\unicode[STIX]{x1D712}/\unicode[STIX]{x1D716}$ of approximately 0.54 indicated by figure 7.

Figure 8.

Mixing efficiency $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D712}/(\unicode[STIX]{x1D712}+\unicode[STIX]{x1D716})$ .

There is some evidence that $\unicode[STIX]{x1D702}$ depends on the Prandtl (or Schmidt) number. For example, from numerical simulations of stratified shear flows, Smyth et al. (Reference Smyth, Moum and Caldwell2001) found that $\unicode[STIX]{x1D702}$ decreased from 0.35 to 0.22 as the Prandtl number increased from 1 to 7. More recently Salehipour & Peltier (Reference Salehipour and Peltier2015) found in numerical simulations of stratified shear flows that $\unicode[STIX]{x1D702}$ decreased from 0.24 to 0.18 as the Prandtl number increased from 1 to 16. Salehipour & Peltier (Reference Salehipour and Peltier2015) also gave some physical and mathematical reasoning for this decrease with Prandtl number. For a subset of the simulations of de Bruyn Kops, Riley & Portwood (Reference de Bruyn Kops, Riley and Portwood2016) run with $Pr=7$ , $\unicode[STIX]{x1D702}$ ranges from 0.17 at $Fr_{h}=0.4$ to 0.28 at $Fr_{h}=2.8$ .

3.1.5 The buoyancy Reynolds number and activity parameter

Dimensional reasoning suggests that three initial parameters are required to describe decaying, homogeneous, stratified turbulence. For equations (2.2), these are chosen to be the Prandtl number, $Pr$ , in addition to the Reynolds and Froude numbers, $Re$ and $Fr$ . The latter are defined here in terms of the integral length scale $\hat{L}$ and r.m.s. velocity $\hat{U}$ of the simulated flow at $t=1$ , i.e. before there is any effect of stratification. The conditions at $t=1$ provide a useful basis for non-dimensionalising the flow statistics. To better understand the effects of stable stratification, however, it is useful to consider the local in time values of $Re$ and $Fr$ . Because of the non-isotropy which develops in the flows, the corresponding vertical and horizontal scales in the flow differ significantly, as seen in figure 3. Therefore, considering the local r.m.s. horizontal velocity $u_{h}$ , one can define two Froude numbers of potential importance, $Fr_{h}=2\unicode[STIX]{x03C0}u_{h}/NL_{ux}$ and $Fr_{v}=2\unicode[STIX]{x03C0}u_{h}/NL_{uz}$ based upon the horizontal longitudinal integral scale $L_{ux}$ and the vertical lateral integral scale $L_{uz}$ , respectively.

Figure 9.

Horizontal and vertical Froude numbers, $Fr_{h}$ and $Fr_{v}$ respectively, versus time $t$ .

The temporal decay of the horizontal Froude number, $Fr_{h}$ , and the vertical Froude number, $Fr_{v}$ , are seen in figure 9 for all three cases. The decay curves for both Froude numbers are very similar for all three different Reynolds numbers, although at later times the Froude numbers for Case I are decaying slightly faster. Furthermore, as suggested by Riley et al. (Reference Riley, Metcalfe, Weissman and West1981), $Fr_{h}<O(1)$ for about $T>1$ , indicating that stratification begins to dominate in this temporal region. This is consistent with the results in figures 2 and 3, where the velocity decay rates and integral scale growth rates change somewhat abruptly near $T=1$ , as the flows enter the strongly stratified regime. On the other hand, the scaling arguments of Billant & Chomaz (Reference Billant and Chomaz2001) suggest that when $Fr_{h}<O(1)$ , at least for large $Re$ , then $Fr_{v}$ should become and remain $O(1)$ . This is also seen to occur in figure 9. More discussion of the late time behaviour of both $Fr_{h}$ and $Fr_{v}$ is presented in § 3.2.

Since $Fr_{v}$ becomes approximately constant for $T>1$ , its local in time behaviour need not be considered. Therefore, it is important to locally consider only $Re_{h}$ and $Fr_{h}$ , where $Re_{h}$ is also based upon $u_{h}$ and $L_{ux}$ . As reviewed in § 1, these two parameters combine to form ${\mathcal{R}}=Fr_{h}^{2}Re_{h}$ , whose potential relevance is based on the shearing dynamics observed to occur when strong stratification leads to the tendency for vertical decoupling of horizontal motions (Riley & de Bruyn Kops Reference Riley and de Bruyn Kops2003).

In the stratified turbulence literature, ${\mathcal{R}}$ is often referred to as the ‘buoyancy Reynolds number’, a terminology we have adopted in this paper. In the geophysical literature, this terminology usually denotes

Here the turbulence Froude and Reynolds numbers, $Fr_{t}$ and $Re_{t}$ , are defined in terms of the turbulence length scale $\hat{\ell }=\hat{U} ^{3}/\hat{\unicode[STIX]{x1D716}}$ , with $\hat{U}$ being the total r.m.s. velocity. These are the definitions used in Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007); however, they refer to $Gn$ as the buoyancy Reynolds number. (Note then that ${\mathcal{R}}$ and $Gn$ differ in the choice of their characteristic velocity and characteristic length scales.) We use the symbol $Gn$ in deference to the introduction of this quantity by Gibson (Reference Gibson and Nihous1980), and to its identification by Gargett, Osborn & Nasmyth (Reference Gargett, Osborn and Nasmyth1984) as a measure of the scale separation available for turbulence between the Ozmidov, $\hat{L}_{o}=(\hat{\unicode[STIX]{x1D716}}/\hat{N}^{3})^{1/2}$ , and Kolmogorov, $\hat{L}_{k}=(\unicode[STIX]{x1D708}^{3}/\hat{\unicode[STIX]{x1D716}})^{1/4}$ , length scales; we call $Gn$ the ‘activity parameter’ following some of the early literature on the subject (e.g. Dillon & Caldwell Reference Dillon and Caldwell1980).

Figure 10.

Values of ${\mathcal{R}}$ and $Gn$ , normalised by the initial Reynolds number, versus  $t$ .

In figure 10 the evolution of ${\mathcal{R}}$ and $Gn$ in time are shown. Note that, as seen in terms of their decay rates in figure 10, ${\mathcal{R}}$ and $Gn$ are not proportional in the present simulations, so that they cannot be used interchangeably. Also note that there is a change in slope in ${\mathcal{R}}$ at about $T=1$ , consistent with the results for $u$ and $v$ seen in figure 2. It is also apparent that ${\mathcal{R}}$ scaled by the initial Reynolds number evolves differently for Case I than for Cases II and III. All three cases evolve in a similar manner up to about $T=2$ , at which point Cases II and III appear to diverge from Case I.

Our hypothesis is that the differences in the evolution of the three cases, evident in figure 10, can be explained by when $Gn$ drops below $O(1)$ , as will be discussed more in § 3.2 on late time results. Gibson (Reference Gibson and Nihous1980) predicted that $Gn>30$ is necessary for sustained turbulence, a value supported experimentally by Stillinger et al. (Reference Stillinger, Helland and Van Atta1983). In simulations of forced homogeneous stratified turbulence with $Gn=13,48$ and 220 (Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), only in the case with $Gn=220$ is there a range of length scales in which the turbulence is approximately isotropic throughout much of the simulation domain, although all the cases exhibit at least some patches of approximately isotropic turbulence. For the same data, de Bruyn Kops (Reference de Bruyn Kops2015) shows that, even in the case with $Gn=220$ , there is no range of scales consistent with the assumptions of Kolmogorov, Oboukhov and Corrsin regarding inertial and convective subranges (Kolmogorov Reference Kolmogorov1941; Oboukhov Reference Oboukhov1941a ,Reference Oboukhov b , Reference Oboukhov1949; Corrsin Reference Corrsin1951; Kolmogorov Reference Kolmogorov1962; Oboukhov Reference Oboukhov1962). For the same three data sets, Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) find that at $Gn=220$ the domain is dominated by space-filling patches of turbulence, while when $Gn=13$ and 48 the simulated flows consist of a mixture of space-filling patches and horizontally layered turbulence. As seen in figures 1 and 11, the flows studied here become very layered beyond approximately $T=2$ .

Figure 11.

Value of $\log _{10}(Gn)$ on a vertical slice of Case III at various times. Only the upper $1/4$ and left $1/2$ of the domain is plotted.

3.2.1 Energy spectra

The horizontal velocity spectra $E_{xx}$ versus horizontal wavenumber $\unicode[STIX]{x1D705}_{x}$ are plotted in figure 13(a), where the spectra of $u$ have been averaged with the corresponding spectra of $v$ . The corresponding compensated horizontal spectra versus horizontal wavenumber times the Kolmogorov length scale $L_{k}$ are plotted in figure 13(b). The horizontal velocity spectra $E_{xx}$ versus vertical wavenumber $\unicode[STIX]{x1D705}_{z}$ are plotted in figure 13(c), again where the spectra of $u$ have been averaged with the corresponding spectra of $v$ . It is seen that, at $T=0$ , when the stratification is applied, the non-stratified spectra have moderate inertial subranges. The compensated spectrum in figure 13(b) is consistent with the Kolmogorov constant of approximately $C_{1}=0.53$ , which is slightly higher than the value of $0.49$ in Pope (Reference Pope2000) but well within the range for the many data sets reviewed by Sreenivasan (Reference Sreenivasan1995).

Figure 13.

Uncompensated (a) and compensated (b) longitudinal spectra and uncompensated vertical spectrum (c) of $u$ averaged with that of $v$ for Case III.

In isotropic turbulence at high enough Reynolds numbers, there is found a bump in the compensated spectrum to the right of $\unicode[STIX]{x1D705}L_{k}\approx 0.02$ , which is sometimes attributed to a ‘bottleneck’ of energy transferring to smaller scales. Meyers & Meneveau (Reference Meyers and Meneveau2008) review efforts to understand and model this phenomenon. de Bruyn Kops (Reference de Bruyn Kops2015) shows that stratification suppresses this bump and that the $-5/3$ slope can occur at the top of the dissipation range even though it does not occur in the inertial range. The suppression of the bump by stratification is consistent with the hypothesis that the bump is due to the bottlenecking of the energy transfer rate combined with the observation that, for the current simulations, stratification reduces the energy decay rate, and hence the energy transfer rate.

Note that in figure 13(b) there is a bottleneck bump to the right of $\unicode[STIX]{x1D705}L_{k}\approx 0.02$ at $T=0$ . By $T=0.5$ the bump has been suppressed and a peak in the compensated spectrum has begun to develop in the range where the inertial-range plateau exists at $T=0$ . By $T=12$ , the compensated spectrum is very similar to those in figure 4(a) of de Bruyn Kops (Reference de Bruyn Kops2015). This is somewhat surprising because it has been assumed that the peak in the compensated spectra is an artefact of the energy forcing used to maintain the statistical stationarity of those simulations, yet in the current simulations there is no forcing. The plots of the compensated horizontal spectra are also qualitatively similar to those obtained by Maffioli & Davidson (Reference Maffioli and Davidson2016), and in addition the growth of the compensated spectra in the low wavenumber region, and the decay in the high wavenumber region, are also consistent.

Figure 14.

Horizontal spectra of various shear terms.

The question arises regarding whether a stratified turbulence inertial range, as proposed by Lindborg (Reference Lindborg2006), should be expected for our simulations. The possibility for this to occur will depend upon (i) the ratio of the energy-containing range scale, characterised by $L_{h}$ , to the Ozmidov scale $L_{o}$ ; and (ii) the ratio of the Ozmidov scale to the Kolmogorov scale, given by $Gn^{4/3}$ . A stratified turbulence inertial range would be expected to exist between the energy-containing and Ozmidov wavenumbers, so that $L_{h}/L_{o}$ should be at least 100 or more. Furthermore, for turbulent behaviour to exist, we find in § 3.2 that $Gn$ must be at least of order 1. Considering the highest Reynolds number case (Case III), we find that at $T=1$ , where the flow becomes strongly affected by stratification, $L_{h}/L_{o}=32.3$ , while $L_{o}/L_{k}=10.1$ . Clearly the ratio of $L_{h}/L_{o}$ is not large enough for an inertial range to exist. When the ratio of $L_{h}/L_{o}$ becomes large enough for the inertial range to exist, at $T=6$ , where $L_{h}/L_{o}=216$ , then $L_{o}/L_{k}$ has dropped to 1.4, which is too small to expect active turbulence. Therefore we conclude that there is not enough separation of scales for a stratified flow inertial range to occur. The flow would have to be started at a substantially higher Reynolds number for the possibility to exist.

A comparison of the horizontal and vertical wavenumber spectra in figures 13(a) and 13(c) give some insight into the turbulence dynamics. The high wavenumber energy in the horizontal decays much more rapidly in time than does the vertical wavenumber energy. On the other hand, the low wavenumber energy in the horizontal decays much more slowly. These facts together indicate a much weaker spectral energy transfer in the horizontal than in the vertical, a fact qualitatively consistent with the arguments of Lilly (Reference Lilly1983) and Billant & Chomaz (Reference Billant and Chomaz2001). Note that at higher wavenumbers the vertical spectra steepen considerably at later times, consistent with the flow layering. A clear $k_{z}^{-3}$ spectrum is not obtained, however.

Some additional insight into the dynamics of strongly stratified flows can be seen from examining the horizontal spectra of various components of the velocity gradient tensor. Figure 14 contains plots of the $k_{x}$ spectra of $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z$ , $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}y$ , $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}y$ and $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ for four points in time for Case III, the highest Reynolds number case. Note that, for isotropic turbulence, it is expected that the spectra of $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z$ and $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}y$ will be identical, as will the spectra of $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}y$ , and $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ . It is clear that at $T=0$ this is the case. As the flow evolves, the spectra of $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}y$ , and $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ remain approximately the same, although the spectrum of $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}y$ decreases slightly faster than that of $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ at higher wavenumbers. This supports the scaling arguments of Billant & Chomaz (Reference Billant and Chomaz2001) that all three terms in the continuity equation, $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}x$ , $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}y$ and $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ , should be of the same order. Somewhat surprisingly, the spectra of $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}y$ and $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ become almost identical at the latest time. What is noteworthy, however, is the rapid decrease with time of the spectra for $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}y$ , $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}y$ and $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}z$ compared to that for $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z$ , with the spectrum for $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z$ evolving well above the rest. Furthermore, the peak in the spectrum of $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z$ remains into smaller wavenumbers. Both of these facts are indications of the development of strong vertical shearing of the horizontal motions as suggested by Lilly (Reference Lilly1983).

We now consider the late time behaviour of these flows, when $Gn$ falls to $O(1)$ and below, and so viscosity becomes a dominating influence. Figure 11 contains grey scale plots, for the same conditions as in figure 1, of the turbulence activity parameter $Gn$ , defined in (3.5), which gives an indication of the strength of the turbulence. Although at $T=0.5$ the effects of stable stratification are little noticeable, by two and especially four buoyancy periods a quasi-horizontal nature and a patchiness of the turbulence are clearly evident. By $T=30$ , at which point $Gn$ is well below 1 (see figure 16 below) the behaviour of the flow has substantially changed.

Figure 15 contains plots of $u^{\prime }$ , $w^{\prime }$ and $E_{p}^{1/2}$ versus time carried out to 150 buoyancy periods or more for each case, and where $Gn$ has dropped to values of 0.01 or less at very late times. The $Gn$ scale at the bottom of each plot gives the local values of $Gn$ at different times for each case. There are two things to note about these plots. First, for each case, for values of $Gn$ between 1 and 0.1, $u^{\prime }$ begins to decay at a slightly steeper rate. Second, for values of $Gn$ in this same range, the vertical velocity and $E_{p}^{1/2}$ (proportional to $\unicode[STIX]{x1D70C}^{\prime }$ ) begin to decay precipitously, indicating that the flows are becoming quasi-horizontal.

Figure 15.

Evolution of $u^{\prime }$ , $v^{\prime }$ , and $E_{p}^{1/2}$ into late time. Reference slopes are shown for $Gn>1$ and $Gn\ll 1$ except for Case I $Gn>1$ which is not in the range of times on this figure. For velocity data at early times see figure 2.

These plots can be better understood by examining the trajectory of each flow in the $Fr_{h}$ – $Gn$ plane, as given in figure 16. Each flow is seen to go through three regimes. In the first, the flows decay at the non-stratified rate, which is seen in figure 2 to occur up to approximately 1 buoyancy period after flow initiation, when the values of the local Froude numbers have fallen to approximately 0.6 (see figure 9). We take this to be the value of $Fr_{h}$ when a flow has entered the strongly stratified regime. The decay rate in this regime diminishes due to the effects of stable stratification, approximately as predicted by Davidson (Reference Davidson2010). Finally, in the third regime, when $Gn$ has decayed to between 1 and 0.1 for each case, the flows become viscously dominated and begin to decay at a faster rate. Note that this viscous domination had already been observed for Case I at earlier times in the behaviour of $\unicode[STIX]{x1D70C}^{\prime }$ in figure 4, in the behaviour of the various components of energy in figure 5(a), and in the behaviour of the mixing efficiency in figure 8.

Figure 16.

The trajectory of each flow in the $Fr_{h}$ – $Gn$ plane. The symbols mark the time in buoyancy periods, with the corresponding times indicated on the curve for Case III.

These results are analogous to those obtained in the experiments of Spedding (Reference Spedding1997) discussed in § 1. He studied the decay of turbulence in the wake of a towed sphere in a salt-stratified towing tank, and found in measuring the mean axial velocity that, initially, the flow decayed at the non-stratified rate. But at a time of the order of a buoyancy period, the decay rate decreased appreciably in what he called the non-equilibrium regime. Further in time, in the region he called the late wake, or quasi-two-dimensional regime, the flow again decayed at a faster rate. We speculate that in this last regime $Gn$ had become small in his flows, which would therefore be viscously dominated. Furthermore the vertical velocity had been suppressed compared to the horizontal velocity, as the flows exhibited quasi-horizontal vortex structure.

Figure 17.

Froude numbers versus time.

Additional insight can be gained by examining the behaviours of the Froude numbers, based upon the horizontal and vertical integral scales, throughout the flow decay, as shown in figure 17. The Froude number based upon the horizontal integral scale, $Fr_{h}$ , is seen to continually decay, although the decay is somewhat diminished when the effects of stratification become strong, and then the decay becomes stronger as viscous effects begin to dominate. At first the Froude number based upon the vertical scale, $Fr_{v}$ , begins to decay at a rate similar to that of $Fr_{h}$ . As the effects of stratification become strong, however, $Fr_{v}$ becomes approximately constant, as predicted by the scaling arguments of Billant & Chomaz (Reference Billant and Chomaz2001), and as seen in the simulations of Maffioli & Davidson (Reference Maffioli and Davidson2016). Later in time, as the effects of viscosity begin to dominate, then $Fr_{v}$ begins to decay.

Based upon our simulation results, we find that it is most useful and instructive to divide up the flows into three different regimes, depending upon the values of $Fr_{h}$ and $Gn$ . When $Fr_{h}>O(1)$ the flows behave almost as non-stratified flows and we term this regime the weakly stratified regime. When $Fr_{h}$ becomes less than $O(1)$ , in our case less than approximately 0.6 which occurs at approximately one buoyancy period, stratification begins to strongly affect the flows, and we term this the strongly stratified regime, as discussed in § 1. Depending on the local value of $Gn$ , the strongly stratified regime splits into two parts. When $Gn$ is greater than $O(1)$ , the flow obeys the decay laws predicted by Davidson (Reference Davidson2010), and the scaling arguments predicted by Billant & Chomaz (Reference Billant and Chomaz2001), and we call this the stratified turbulence regime, in deference to the terminology introduced by Lilly (Reference Lilly1983). Finally, in the strongly stratified regime, when $Gn$ has dropped somewhat below 1, the flows become quasi-horizontal, and the dissipation rate is approximately governed by the vertical shearing of the horizontal velocity, and we refer to this as the viscous-dominated regime.

The behaviour of the flow in this late time, viscous-dominated regime can be seen from the scaling arguments of Riley et al. (Reference Riley, Metcalfe, Weissman and West1981), as modified by Godoy-Diana et al. (Reference Godoy-Diana, Chomaz and Billant2004) to include the effects of viscosity. The arguments suggest that the ratio of the r.m.s. vertical velocity to the r.m.s. horizontal velocity should be proportional to $Fr_{h}/Fr_{v}$ , so that the flow should be quasi-horizontal in this regime. Godoy-Diana et al. (Reference Godoy-Diana, Chomaz and Billant2004) also argued that as $Gn$ becomes small, the vertical viscous term in the horizontal momentum equation becomes large compared to the vertical transport term, suggesting that the vertical integral scale behaves as $L_{v}=L_{h}Re_{h}^{-1/2}$ , where $Re_{h}$ is defined in terms of $L_{h}$ . Figure 18 is a plot of the ratio $Re_{h}^{1/2}L_{v}/L_{h}$ versus $Gn$ for Cases I, II and III. It is seen that for approximately $Gn>O(1)$ , this ratio is not constant, but increases as $Gn$ increases. As $Gn$ decreases below $O(1)$ , however, the ratio becomes approximately constant with a value slightly below 5, and in fact the curves for Cases I and II collapse together, although the value for Case I is somewhat lower. This gives some confidence in the arguments of Godoy-Diana et al. (Reference Godoy-Diana, Chomaz and Billant2004) applies to these simulated flows.

Figure 18.

A test of the relationship $L_{v}=L_{h}Re_{h}^{-1/2}$ for small $Gn$ .

Examining these arguments further, neglecting the vertical transport term and assuming that the aspect ratio $L_{v}/L_{h}$ is small as suggested by the scaling arguments and by figure 17, so that the horizontal diffusion terms can be neglected, the result is that (2.2) reduce to

There are several indications of the appropriateness of these equations to describe the late time behaviour of flows with small $Gn$ . One is the fact that the flows are approximately horizontal, i.e.  $u^{\prime }$ , $v^{\prime }\gg w^{\prime }$ , which in the present case is clearly indicated in figure 15. Another indication can be found by examining the ratios For (3.7) to have near validity, $Rt_{G}$ must be approximately $15/4$ . In contrast, in isotropic turbulence, all of the preceding ratios should be unity. Both of these conditions are observed in figure 19, which shows the ratios for case II. At early times, each of the ratios is approximately unity whereas at late time $Rt_{G}$ appears to by asymptoting to $15/4$ .

Figure 19.

Plot of the five statistically independent contributions to $\unicode[STIX]{x1D716}$ , normalised as in (3.8), for Case II.

The results from other studies exhibit this behaviour of the dissipation rate as well. Both Fincham et al. (Reference Fincham, Maxworthy and Spedding1996) and Praud et al. (Reference Praud, Fincham and Sommeria2005) performed experiments in salt-stratified tanks, and both generated turbulence by towing a vertical grid through the tank. At the later stages of decay in both experiments the value of $Gn$ had dropped to values well below one. In addition, the flows had become approximately horizontal, and in both experiments it was reported that the ratio $Rt_{G}/(15/4)$ was well above 0.9. These features for both experiments all indicate the approximate validity of (3.7) at later times. Hebert & de Bruyn Kops (Reference Hebert and de Bruyn Kops2006b ) performed direct numerical simulations in a strongly stratified fluid of the decay of flows initially consisting of Taylor–Green vortices, similar to the work of Riley & de Bruyn Kops (Reference Riley and de Bruyn Kops2003), but over a very broad range of Reynolds number, Froude number and $Gn$ . They found that the ratio $Rt_{G}$ is a very strong function of $Gn$ , and that $Rt_{G}/(15/4)$ being near 0.9 only occurred when $Gn$ was order one or less. From the present results, and from these experimental and numerical results, there is reason to speculate that (3.7) hold in many other flows that are strongly stratified and $Gn$ is small.