2 Theory and background

We consider an incompressible fluid with a velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ and a density field determined by a perturbation $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)$ to a constant background linear stratification. We apply the Boussinesq approximation that density changes are negligible compared to the mean density and furthermore assume that the density field has a linear equation of state and hence satisfies an advection–diffusion equation. The flow is thus governed by the Navier–Stokes equations in the form

(2.1)$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-Ri_{0}\unicode[STIX]{x1D70C}\hat{\boldsymbol{z}}+\boldsymbol{F}_{\boldsymbol{u}}, & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D70C}=\frac{1}{RePr}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D70C}+w+F_{\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$

where we have non-dimensionalised the equations using length and velocity scales $L_{0}$ and $U_{0}$, and $\hat{\boldsymbol{z}}$ is the unit vector in the vertical direction. The density perturbation has also been non-dimensionalised by $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$. External forcing acting on the velocity and density fields are denoted by $\boldsymbol{F}_{\boldsymbol{u}}$ and $F_{\unicode[STIX]{x1D70C}}$, the precise forms of which are detailed in the next section. The three dimensionless parameters in the equations are the Reynolds number, Prandtl number and bulk Richardson number

(2.4a-c)$$\begin{eqnarray}Re=\frac{L_{0}U_{0}}{\unicode[STIX]{x1D708}},\quad Pr=\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D705}},\quad Ri_{0}=\frac{g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}L_{0}}{\unicode[STIX]{x1D70C}_{0}{U_{0}}^{2}}=\frac{{N_{0}}^{2}{L_{0}}^{2}}{{U_{0}}^{2}},\end{eqnarray}$$

where $\unicode[STIX]{x1D708}$ is the kinematic viscosity, $\unicode[STIX]{x1D705}$ is the density diffusivity and $\unicode[STIX]{x1D70C}_{0}$ is the mean density. The constant background density gradient is $-\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/L_{0}$, which can be used to define $N_{0}=\sqrt{g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}L_{0}}$ as a background buoyancy frequency. Since $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ also acts as the dimensional scale for the density perturbation $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)$, the Boussinesq approximation requires that $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\ll \unicode[STIX]{x1D70C}_{0}$. For clarity, the full dimensional density field will be written as $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}[-z+\unicode[STIX]{x1D70C}(\boldsymbol{x},t)]$ in this formalism, where $z=z^{\ast }/L_{0}$ and $z^{\ast }$ are respectively the dimensionless and dimensional vertical coordinates.

Stably stratified flows are commonly anisotropic, with horizontal length scales much larger than vertical length scales. When analysing these flows, it is therefore natural to consider the decomposition of the velocity and density fields into horizontally averaged mean quantities and perturbations from them. We will use the following notation, denoting mean quantities with an overbar and perturbations with a prime, i.e.

(2.5)$$\begin{eqnarray}f(\boldsymbol{x},t)=\overline{f}(z,t)+f^{\prime }(\boldsymbol{x},t),\quad \overline{f}(z,t)=\frac{1}{L_{x}L_{y}}\int _{0}^{L_{x}}\int _{0}^{L_{y}}f(\boldsymbol{x},t)\,\text{d}x\,\text{d}y.\end{eqnarray}$$

Taking a horizontal mean of the incompressibility condition (2.1) implies that $\unicode[STIX]{x2202}\overline{w}/\unicode[STIX]{x2202}z=0$, and if $\overline{w}=0$ initially then it remains zero for all time. From now on, we will assume that this is the case and hence that the mean velocity $\overline{\boldsymbol{u}}(z,t)=(u,v,0)$ is purely horizontal.

In this paper we consider flow in a triply periodic domain, which allows us to construct simple equations for the energy of the system from (2.2) and (2.3). Implementing the decomposition (2.5) yields four energy quantities of interest: the kinetic and potential energies (per unit mass) associated with both the mean and perturbation fields, namely

(2.6a,b)$$\begin{eqnarray}\overline{E_{K}}={\textstyle \frac{1}{2}}\langle |\overline{\boldsymbol{u}}|^{2}\rangle ,\quad E_{K}={\textstyle \frac{1}{2}}\langle |\boldsymbol{u}^{\prime }|^{2}\rangle ,\end{eqnarray}$$

(2.7a,b)$$\begin{eqnarray}\overline{E_{P}}=\frac{Ri_{0}}{2}\langle \overline{\unicode[STIX]{x1D70C}}^{2}\rangle ,\quad E_{P}=\frac{Ri_{0}}{2}\langle \unicode[STIX]{x1D70C}^{\prime 2}\rangle ,\end{eqnarray}$$

where $\langle \cdot \rangle$ denotes a volume average. This positive semi-definite form of potential energy is valid since $\unicode[STIX]{x1D70C}$ is a departure from a linear background profile. Motivated by the ‘pancake vortices’ description of stratified turbulence, we refer to $E_{K}$ as the turbulent kinetic energy and to $E_{P}$ as the turbulent potential energy. Since $\overline{w}=0$ and $\langle \boldsymbol{u}\rangle =\mathbf{0}$, the mean kinetic energy $\overline{E_{K}}$ is often associated with what are conventionally referred to as ‘shear modes’, and there is a body of literature investigating its development in stratified turbulence (e.g. Smith & Waleffe Reference Smith and Waleffe2002; Augier, Billant & Chomaz Reference Augier, Billant and Chomaz2015).

Multiplying (2.2) and (2.3) by the velocity and density fields respectively leads to the following evolution equations for the energy:

(2.8a,b)$$\begin{eqnarray}\frac{\text{d}\overline{E_{K}}}{\text{d}t}=-S_{p}-\overline{\unicode[STIX]{x1D700}},\quad \frac{\text{d}E_{K}}{\text{d}t}=S_{p}-\unicode[STIX]{x1D700}+J_{b}+P_{K},\end{eqnarray}$$

(2.9a,b)$$\begin{eqnarray}\frac{\text{d}\overline{E_{P}}}{\text{d}t}=-N_{p}-\overline{\unicode[STIX]{x1D712}},\quad \frac{\text{d}E_{P}}{\text{d}t}=N_{p}-\unicode[STIX]{x1D712}-J_{b}+P_{P}.\end{eqnarray}$$

The terms on the right-hand side of these equations act as inputs, exchanges and outputs of energy for the system as sketched in figure 1 and detailed below.

Figure 1. Schematic detailing the energy pathways.

The rate at which energy is dissipated by the flow is quantified by the expressions

(2.10a,b)$$\begin{eqnarray}\overline{\unicode[STIX]{x1D700}}=\frac{1}{Re}\left\langle \left|\frac{\unicode[STIX]{x2202}\overline{\boldsymbol{u}}}{\unicode[STIX]{x2202}z}\right|^{2}\right\rangle ,\quad \unicode[STIX]{x1D700}=\frac{1}{Re}\left\langle \frac{\unicode[STIX]{x2202}{u_{i}}^{\prime }}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}{u_{i}}^{\prime }}{\unicode[STIX]{x2202}x_{j}}\right\rangle ,\end{eqnarray}$$

(2.11a,b)$$\begin{eqnarray}\overline{\unicode[STIX]{x1D712}}=\frac{Ri_{0}}{RePr}\left\langle \left|\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}z}\right|^{2}\right\rangle ,\quad \unicode[STIX]{x1D712}=\frac{Ri_{0}}{RePr}\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}x_{j}}\right\rangle ,\end{eqnarray}$$

where $\unicode[STIX]{x1D700}$ is commonly known as the turbulent kinetic energy (TKE) dissipation rate. It is important to appreciate that these various rates are defined in terms of volume averages over the whole computational domain.

Energy can be exchanged between kinetic and potential energy through the buoyancy flux $J_{b}$. Since we have assumed that the mean flow is purely horizontal, this exchange can only take place between the turbulent energies $E_{K}$ and $E_{P}$. The small-scale turbulence instead interacts with the mean flow via the turbulent shear production $S_{p}$, and by an analogous term that appears in the potential energy equations which we refer to as buoyancy production $N_{p}$. The three energy exchange terms are defined

(2.12a-c)$$\begin{eqnarray}J_{b}=-Ri_{0}\langle w^{\prime }\unicode[STIX]{x1D70C}^{\prime }\rangle ,\quad S_{p}=-\left\langle \overline{w^{\prime }\boldsymbol{u}^{\prime }}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\overline{\boldsymbol{u}}}{\unicode[STIX]{x2202}z}\right\rangle ,\quad N_{p}=-Ri_{0}\left\langle \overline{w^{\prime }\unicode[STIX]{x1D70C}^{\prime }}\,\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}z}\right\rangle .\end{eqnarray}$$

The energy input provided by the forcing is prescribed to not act directly on the mean flow, so the energy input rates only appear in the perturbation energy equations and are defined as

(2.13a,b)$$\begin{eqnarray}P_{K}=\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{F}_{\boldsymbol{ u}}\rangle ,\quad P_{P}=Ri_{0}\langle \unicode[STIX]{x1D70C}^{\prime }F_{\unicode[STIX]{x1D70C}}\rangle .\end{eqnarray}$$

This is consistent with the forcing terms used in our numerical simulations.

In the formulation, we have chosen to retain flexibility to force both the velocity and density fields. In particular, whether or not the density field has explicit forcing has important implications for the energy budget of a quasi-steady turbulent state when $\text{d}E_{P}/\text{d}t\approx 0$. It is worth noting that in the flows considered by this study the buoyancy production $N_{p}$ is typically much smaller than the other terms on the right-hand side of (2.9), so if there is no density forcing then the turbulent potential energy budget reads

(2.14)$$\begin{eqnarray}0\approx -\unicode[STIX]{x1D712}-J_{b},\end{eqnarray}$$

in a steady state. Since $\unicode[STIX]{x1D712}$ is positive semi-definite this implies that $J_{b}\leqslant 0$ and so the buoyancy flux acts to transfer energy from kinetic to potential. When $F_{\unicode[STIX]{x1D70C}}\neq 0$ this restriction is not enforced, and we will investigate the effect of introducing density forcing on the energy pathways later on. This flexibility is also physically motivated since nonlinear interactions between internal gravity waves can transfer kinetic and potential energy across spatial scales.

As noted in the introduction, we are interested in defining an appropriate measure of the ‘efficiency’ of mixing. We would wish to define an instantaneous mixing efficiency $\unicode[STIX]{x1D702}$ as the proportion of energy lost by turbulence that leads to irreversible mixing. In our triply periodic domain it is unfortunately technically impossible to unambiguously define a background potential energy quantity as is used to quantify irreversible mixing in, for example, Peltier & Caulfield (Reference Peltier and Caulfield2003). We therefore treat $E_{P}$ as a proxy for available potential energy and use $\unicode[STIX]{x1D712}$ (as defined in (2.11)) to quantify the irreversible loss of $E_{P}$ that leads to mixing, yielding

(2.15)$$\begin{eqnarray}\unicode[STIX]{x1D702}:=\frac{\unicode[STIX]{x1D712}}{\unicode[STIX]{x1D712}+\unicode[STIX]{x1D700}},\end{eqnarray}$$

as the expression for mixing efficiency. The denominator $\unicode[STIX]{x1D712}+\unicode[STIX]{x1D700}$ represents the total instantaneous energy lost due to turbulence, and specifically excludes any laminar diffusion of the mean flow through $\overline{\unicode[STIX]{x1D712}}$ or $\overline{\unicode[STIX]{x1D700}}$. We focus on mixing by turbulence because geophysical flows will typically occur at much larger Reynolds numbers than we can accurately simulate, leading to negligible laminar diffusion. Even at our modest $Re$ the results are qualitatively unchanged by including the dissipation of the mean flow, with the average mixing efficiency only decreasing by between 4 % and 8 % across the simulations.

The mixing efficiency $\unicode[STIX]{x1D702}$ is closely related to $\unicode[STIX]{x1D6E4}$, the turbulent flux coefficient commonly used in oceanography to infer measures of mixing from observations. The original definition of Osborn (Reference Osborn1980) postulates a linear relationship between the buoyancy flux and the turbulent dissipation rate in a quasi-steady state of fully developed turbulent flow, with $\unicode[STIX]{x1D6E4}$ as the constant of proportionality. However, since we are considering flows where the buoyancy flux may well have a significant reversible component associated with internal waves, we believe it is more appropriate to define $\unicode[STIX]{x1D6E4}$ in terms of the ratio between $\unicode[STIX]{x1D712}$ and $\unicode[STIX]{x1D700}$, so that

(2.16)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}:=\frac{\unicode[STIX]{x1D712}}{\unicode[STIX]{x1D700}}=\frac{\unicode[STIX]{x1D702}}{1-\unicode[STIX]{x1D702}}.\end{eqnarray}$$

Recent stratified turbulence studies by Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) and Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) have derived scalings for such a defined $\unicode[STIX]{x1D6E4}$ in terms of the turbulent Froude number

(2.17)$$\begin{eqnarray}Fr:=\frac{\unicode[STIX]{x1D700}}{{Ri_{0}}^{1/2}E_{K}}.\end{eqnarray}$$

In the strongly stratified regime associated with $Fr\ll O(1)$ the scalings and associated simulations suggest that $\unicode[STIX]{x1D6E4}$ is independent of $Fr$. Although this might be thought to provide some justification for the use of a constant $\unicode[STIX]{x1D6E4}$ to infer mixing in geophysical flows, there is most definitely no consensus in the fluid dynamical community as to what value $\unicode[STIX]{x1D6E4}$ takes in this regime, the region of validity of the regime or indeed the variability of $\unicode[STIX]{x1D6E4}$ outside of this regime.

In particular, there is an alternative approach to parameterisation based around the argument that the appropriate parameter to use is the buoyancy Reynolds number

(2.18)$$\begin{eqnarray}Re_{b}:=\frac{\unicode[STIX]{x1D700}Re}{Ri_{0}}\end{eqnarray}$$

(see for example Monismith, Koseff & White Reference Monismith, Koseff and White2018), which can be considered to quantify how ‘energetic’ the turbulence is. Monismith et al. (Reference Monismith, Koseff and White2018) present data from numerical simulations and energetic near-shore flow observations that suggest the mixing efficiency scales as $\unicode[STIX]{x1D702}\sim Re_{b}^{-1/2}$ when the flow is ‘energetic’, which may also loosely be thought of as being weakly stratified, defined as $Re_{b}>O(100)$.

Monismith et al. (Reference Monismith, Koseff and White2018) still support the hypothesis that $\unicode[STIX]{x1D6E4}$ is constant in strongly stratified flows with $Re_{b}<O(100)$, taking an approximate value of 0.2, although we caution associating smaller values of $Re_{b}$ with ‘strong’ stratification, as smaller values of $Re_{b}$ should really be considered to be associated with flows which are viscously dominated, or at least viscously affected. Gregg et al. (Reference Gregg, D’Asaro, Riley and Kunze2018) also argue in favour of using the estimate $\unicode[STIX]{x1D6E4}\simeq 0.2$ which dates back to the early parameterisation of Osborn (Reference Osborn1980). They, however, caution the use of $Re_{b}$ as a sole parameter for the functional dependence of $\unicode[STIX]{x1D6E4}$, and note that the turbulence produced by internal waves typically has $Re_{b}\lesssim 200$, where $\unicode[STIX]{x1D6E4}$ is thought to be constant. Indeed, it is not at all clear at the moment whether the independence of $\unicode[STIX]{x1D6E4}$ with respect to $Fr$ when $Fr\ll 1$ is in any way associated with the classic empirically useful estimate that $\unicode[STIX]{x1D6E4}\simeq 0.2$, not least because it is exceptionally computationally challenging to consider flows with $Fr\ll 1$ and larger values of $Re_{b}$.

Since we are primarily concerned with internal wave-driven mixing in the open ocean, we choose to focus on ‘strongly stratified’ flows associated with $Fr\ll O(1)$, rather than classifying flows in terms of $Re_{b}$. Even in this regime, there are still discrepancies in the value of $\unicode[STIX]{x1D6E4}$ between observations and numerical simulations. Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) find a trend towards the constant value of $\unicode[STIX]{x1D6E4}=0.33$ from their simulations of forced stratified turbulence, and simulations of decaying stratified turbulence by de Bruyn Kops & Riley (Reference de Bruyn Kops and Riley2019) have shown sustained values of $\unicode[STIX]{x1D6E4}$ as large as 0.54. Understanding how these discrepancies may arise is vital if we hope to relate these numerical studies to observed mixing in the ocean.

One issue in comparing observations with these scaling arguments is that $Fr$ is rarely recorded and requires simultaneous measurement of multiple turbulent quantities. Observationally it is easier to obtain the fundamental length scales named after Ellison and Ozmidov, defined as

(2.19a,b)$$\begin{eqnarray}L_{E}:=\frac{{\unicode[STIX]{x1D70C}^{\ast }}_{rms}}{|\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}^{\ast }}/\unicode[STIX]{x2202}z|},\quad L_{O}:=\left(\frac{\unicode[STIX]{x1D700}}{N^{3}}\right)^{1/2}.\end{eqnarray}$$

For example, Ivey, Bluteau & Jones (Reference Ivey, Bluteau and Jones2018) use a mixing length model to infer diapycnal diffusivity from $L_{E}$ and the mean shear measured by moorings. They find good agreement with microstructure measurements, but it is unclear whether this estimate would work well in regions where the background state is dominated by the internal wave spectrum rather than a mean shear. Ivey & Imberger (Reference Ivey and Imberger1991) use $L_{E}/L_{O}$ more generally to infer $Fr$, and therefore determine whether a flow is strongly turbulent or significantly affected by stratification. Many observational studies, however, assume that these length scales are approximately equal, as originally postulated by Dillon (Reference Dillon1982), based on limited experimental data and due to restrictions in measurement equipment. A detailed comparison of these length scales in the thermocline can be found in Moum (Reference Moum1996).

3 Numerical simulations

We use the Diablo software (Taylor Reference Taylor2008) to perform three-dimensional numerical simulations of (2.1)–(2.3). The software implements pseudo-spectral methods to calculate spatial derivatives and a third-order Runge–Kutta scheme for time stepping. The equations are solved in a cubic domain of length $2\unicode[STIX]{x03C0}$ represented by a uniformly spaced grid of $1024^{3}$ points. A 2/3 rule is applied for dealiasing the calculation of nonlinear terms. Periodic boundary conditions are applied in every direction to the velocity and density fields $\boldsymbol{u}$ and $\unicode[STIX]{x1D70C}$. Recall that $\unicode[STIX]{x1D70C}$ represents the density perturbation so periodicity in the vertical does not contradict our use of a stable background density gradient. Table 1 summarises the input parameters used across all simulations.

Table 1. Input parameters for the numerical simulations.

Motivated by the existence of a background internal wave field in the ocean, we construct the initial condition for the simulations as follows. Computational constraints mean that we cannot resolve the range of scales required to represent a full Garrett–Munk spectrum in our domain. We therefore take an approach similar to that of Furue (Reference Furue2003) to construct an initial state where the large scales of the flow field are representative of the small-scale portion of the GM spectrum as defined by Munk (Reference Munk, Warren and Wunsch1981). To obtain the desired vertical energy spectrum of $E\sim m^{-2}$ we need to account for waves with horizontal wavelengths larger than the domain. Furue (Reference Furue2003) achieves this by integrating the GM spectrum over small horizontal wavenumbers to obtain a shear flow containing all of the ‘missed’ energy. We simply define the initial shear as a sum of shear modes $\boldsymbol{u}_{0}\sim (A/m)\text{e}^{\text{i}mz}$ that give an energy spectrum of $m^{-2}$. The shear modes are large scale in the domain and are thus limited to $m\leqslant m_{c}=7$. Each shear mode is randomly phased, and the total energy in this component is normalised such that the mean gradient Richardson number $Ri_{g}=Ri_{0}/\langle S^{2}\rangle$ is equal to $Ri_{0}$. Here $S^{2}=(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z)^{2}+(\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}z)^{2}$ is the dimensionless squared shear, and $\langle \cdot \rangle$ denotes a volume average (simply equivalent to a vertical average in this case). The shear component is complemented by a collection of randomly phased internal waves that satisfy the three-dimensional GM energy spectrum $E(\boldsymbol{k})$ defined in Furue (Reference Furue2003). These waves contribute 10 % of the initial energy and are non-zero for $|\boldsymbol{k}|\leqslant 7$.

We numerically integrate the system for approximately 20 time units without body forcing to allow the initial transient dynamics to dissipate, and for the associated dissipation rate to reach its maximum value. From this state we perform three simulations, each with a different form of body forcing applied. All three types of forcing can be expressed as

(3.1a,b)$$\begin{eqnarray}\boldsymbol{F}_{\boldsymbol{u}}=\mathop{\sum }_{\substack{ 2.5\leqslant |\boldsymbol{k}|\leqslant 3.5 \\ \unicode[STIX]{x1D705}\neq 0}}\widetilde{\boldsymbol{F}_{\boldsymbol{u}}}(\boldsymbol{k})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}},\quad F_{\unicode[STIX]{x1D70C}}=\mathop{\sum }_{\substack{ 2.5\leqslant |\boldsymbol{k}|\leqslant 3.5 \\ \unicode[STIX]{x1D705}\neq 0}}\widetilde{F_{\unicode[STIX]{x1D70C}}}(\boldsymbol{k})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}},\end{eqnarray}$$

where $\boldsymbol{k}=(k,l,m)$ is the wave vector and $\unicode[STIX]{x1D705}=\sqrt{k^{2}+l^{2}}$ is the horizontal wavenumber.

The first type of forcing we consider is that used by Maffioli (Reference Maffioli2017). We refer to this forcing as case H since the forcing acts purely on the horizontal components of velocity and therefore $F_{w}=F_{\unicode[STIX]{x1D70C}}=0$. Forcing H is representative of vertically uniform ‘vortical modes’ with $(\widetilde{F_{u}},\widetilde{F_{v}})\propto (l,-k)$ and the modes being non-zero only when $m=0$. Each mode is randomly phased at every time step.

The other two types of forcing are intended to be representative of flows induced by internal gravity waves, with the forcing components satisfying the internal wave polarisation relations

(3.2a-c)$$\begin{eqnarray}(\widetilde{F_{u}},\widetilde{F_{v}})={\mathcal{A}}\frac{(k,l)m}{\unicode[STIX]{x1D705}|\boldsymbol{k}|},\quad \widetilde{F_{w}}=-{\mathcal{A}}\frac{\unicode[STIX]{x1D705}}{|\boldsymbol{k}|},\quad \widetilde{F_{\unicode[STIX]{x1D70C}}}={\mathcal{A}}\frac{i}{{Ri_{0}}^{1/2}}.\end{eqnarray}$$

We denote one variant of this forcing as case R where the phase of the complex amplitude ${\mathcal{A}}$ for each mode is chosen randomly at every time step. The final type of forcing represents energy input from a propagating wave field and we refer to it as case P. In this case the phase of each ${\mathcal{A}}$ is shifted at time $t$ by $-\unicode[STIX]{x1D714}t$ where the frequency $\unicode[STIX]{x1D714}$ is determined by the linear internal gravity wave dispersion relation, which in our non-dimensionalisation is given by

(3.3)$$\begin{eqnarray}\unicode[STIX]{x1D714}=\frac{Ri_{0}^{1/2}\unicode[STIX]{x1D705}}{|\boldsymbol{k}|}.\end{eqnarray}$$

To ensure that the dissipation rates are comparable across the simulations, we enforce the total energy input rate $P_{K}+P_{P}$ to be constant. We normalise the amplitude of the forcing at each time step to achieve the constant energy input rate shown in table 1. We also use the ‘constant power minimal forcing’ method from Maffioli (Reference Maffioli2017) to avoid large artificial energy inputs arising from discrete time stepping. Each simulation is run for a total of 150 time units.

5 Discussion and conclusions

We have performed direct numerical simulations of stratified turbulence sustained by different forms of large-scale body forcing. The simulation of case H implements the horizontal vortical mode forcing prescribed by Maffioli (Reference Maffioli2017) that injects horizontal kinetic energy into randomly phased columnar vortex modes. The other two simulations force a field of resonant internal gravity waves at large scales of our computational domain as motivated by Furue (Reference Furue2003). The simulation of case R randomly phases each wave at every time step, whereas the simulation of case P shifts the phase of each wave according to the dispersion relation for linear internal gravity waves. Each simulation is initialised with a flow dominated by vertically sheared horizontal modes that is motivated by ambient internal waves with large horizontal wavelengths. The forcing is activated after some initial transient behaviour, and after a further adjustment time the turbulence characterised by the perturbations from the horizontal mean reaches a statistically quasi-steady state.

We find that the quasi-steady states in the wave-forced simulations (cases R and P) have significantly more potential energy than the state achieved by the vortical mode forcing in the simulation of case H. This increased potential energy is provided by the direct forcing of the density field in cases R and P. In the simulation of case H, the irreversible mixing, quantified by the dissipation of the perturbation potential energy, must come via a transfer of energy from the TKE to the potential energy through the buoyancy flux. The density forcing in cases R and P allows this energy transfer to reverse, with mixing made possible without a mean transfer from kinetic to potential energy. This reversal in buoyancy flux can be associated with larger overturning, as seen in figure 2, and thus more convective mixing. The wave-forced simulations also exhibit a higher mixing efficiency than the horizontally forced simulation (case H), which is consistent for flows where mixing occurs through convective mechanisms. The vortical mode forcing used in the simulation of case H forces neither the vertical velocity nor the density field, and therefore cannot produce such large-scale convective overturns.

The qualitatively different energy pathways for mixing in each case also coincide with a change in the interaction between the mean shear flow and the perturbations. Whereas the turbulence extracts energy from the mean flow in the simulation for case H, the wave-forced simulations (cases R and P) show a small transfer in energy from the perturbations to the mean flow on average. This small change partially contributes to the reduced value of the TKE dissipation rate $\unicode[STIX]{x1D700}$ in the wave-forced cases.

The kinetic energy associated with the mean shear is not constant in each case, but varies slowly over time due to the exchange with the perturbation field. In our simulations, the mean shear modes are intended to represent waves at horizontal scales larger than our computational domain. Since the perturbation energy and dissipation rates are statistically quasi-steady, we believe that the small-scale turbulence is still representative of the geophysical flows which motivate us to conduct these idealised simulations. Background shear modes appear in many studies of forced stratified turbulence (Smith & Waleffe Reference Smith and Waleffe2002; Lindborg Reference Lindborg2006; Waite & Bartello Reference Waite and Bartello2006; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007) and are shown not to impact the turbulent dynamics significantly. Furthermore, a recent study by Fitzgerald & Farrell (Reference Fitzgerald and Farrell2018) shows that the emergence of vertically sheared horizontal flows also occurs in a forced 2-D Boussinesq system. This result suggests that energy increase in the shear modes is due to a wave–mean flow interaction, which may explain why we observe the energy increase in the wave-forced cases but not from the horizontal vortical mode forcing utilised in case H.

In every simulation, the turbulent dissipation organises into quasi-horizontal layers. The vertical location of these layers varies depending on the forcing type, but it is currently unclear what determines the change in vertical structure between the simulations. Initial analysis shows that regions of high dissipation do not simply correlate with local changes in the background shear and stratification, and thus further research is needed to investigate the mechanisms by which these layers form and are sustained. This predominantly vertical variation in the dissipation rate allows us to investigate correlations between the turbulent dissipation rate and the density variance destruction rate over orders of magnitude. Intriguingly, we find that the ratio between the two, a local measure of the coefficient $\unicode[STIX]{x1D6E4}$, remains close to the volume-averaged ratio in both regions of high and low dissipation. We deduce that the local mixing efficiency is independent of turbulent intermittency below the scale of the forcing, and instead depends predominantly on the type of large-scale forcing implemented. This further supports the notion that the larger overturnings in the wave-forced cases correspond to a fundamentally different (and in some sense ‘convective’) mixing mechanism from that observed in the simulation of case H.

In the wave-forced simulations the dissipation rate $\unicode[STIX]{x1D700}$ is not correlated with the background stratification, so we also obtain a wide range of values for the horizontally averaged turbulent Froude number $Fr_{H}$ defined in (4.2). This confirms the lack of dependence of $\unicode[STIX]{x1D6E4}$ on $Fr$ in the low-$Fr$ regime for these quasi-steady states. The vortical mode-forced simulation (case H) does not, however, produce such a wide range in $Fr_{H}$. This may be a consequence of the purely horizontal forcing allowing vertical scales to adjust locally, as in the scalings of Billant & Chomaz (Reference Billant and Chomaz2001) where the vertical Froude number $F_{v}$ adjusts to a constant of $O(1)$. Despite the reduced $Fr$ range in case H, there is still at least a hint of $Fr$-dependence for $\unicode[STIX]{x1D6E4}$ in figure 8(a) at the largest values of $Fr_{H}\approx 0.1$. Previous studies (e.g. Lindborg Reference Lindborg2006) have shown that the strongly stratified limit of low $Fr$ requires $Fr\leqslant O(10^{-2})$, suggesting that the observed dependence in case H is outside this regime. The variation in $\unicode[STIX]{x1D6E4}$ is more clearly displayed in figure 12 where, following Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019), henceforth denoted GV19, we plot $\unicode[STIX]{x1D6E4}$ against the length scale ratio $L_{E}/L_{O}$. In particular, we plot these quantities defined in terms of horizontal averages as $\unicode[STIX]{x1D6E4}_{H}=\unicode[STIX]{x1D712}_{H}/\unicode[STIX]{x1D700}_{H}$ and

(5.1)$$\begin{eqnarray}\left(\frac{L_{E}}{L_{O}}\right)_{H}=\frac{{Ri_{0}}^{3/4}\overline{\unicode[STIX]{x1D70C}^{\prime 2}}^{1/2}}{{\unicode[STIX]{x1D700}_{H}}^{1/2}(1-\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}/\unicode[STIX]{x2202}z)^{1/4}}.\end{eqnarray}$$

At first the results shown in figure 12(a) for case H appear inconsistent with any of the scalings proposed by GV19, with $\unicode[STIX]{x1D6E4}\sim (L_{E}/L_{O})^{2}$ for values of $L_{E}/L_{O}$ being $O(1)$. However, this scaling can be reproduced by combining various self-consistent assumptions used in that paper, as follows.

Figure 12. Two-dimensional p.d.f. of $\unicode[STIX]{x1D6E4}$ and $L_{E}/L_{O}$ calculated from horizontally averaged quantities. Blue dashed lines show the mean value of $\unicode[STIX]{x1D6E4}$ in each case and the black dotted line in (a) plots the $(L_{E}/L_{O})^{2}$ scaling found in (5.5). As before, (a) plots data from case H, (b) plots data from case R, (c) plots data from case P.

Firstly, taking $Fr=O(1)$, we assume that the turbulent kinetic energy $E_{K}\sim w^{\prime 2}$ and that its dissipation rate is governed by the turbulent time scale $T_{L}=E_{K}/\unicode[STIX]{x1D700}$, such that $\unicode[STIX]{x1D700}\sim w^{\prime 2}/T_{L}$. We then scale the density field by taking $g\unicode[STIX]{x1D70C}_{rms}^{\prime }/\unicode[STIX]{x1D70C}_{0}$ to be a ‘reduced gravity’ acceleration with velocity scale $w^{\prime }$ and time scale $N^{-1}$. Typical vertical displacements are taken to be $z\sim w^{\prime }/N$ so that the potential energy density $E_{P}=g\unicode[STIX]{x1D70C}^{\prime }z/\unicode[STIX]{x1D70C}_{0}\sim w^{\prime 2}$, and its dissipation rate is governed by the buoyancy time scale $N^{-1}$, giving $\unicode[STIX]{x1D712}\sim E_{P}/T\sim w^{\prime 2}N$. These results provide the scalings

(5.2)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}=\frac{\unicode[STIX]{x1D712}}{\unicode[STIX]{x1D700}}\sim \frac{w^{\prime 2}N}{w^{\prime 2}/T_{L}}=NT_{L}, & \displaystyle\end{eqnarray}$$

(5.3)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{L_{E}}{L_{O}}=\left(\frac{g\unicode[STIX]{x1D70C}_{rms}^{\prime }}{\unicode[STIX]{x1D70C}_{0}}\right)\unicode[STIX]{x1D700}^{-1/2}N^{-1/2}\sim (w^{\prime }N)\left(\frac{w^{\prime 2}}{T_{L}}\right)^{-1/2}N^{-1/2}=(NT_{L})^{1/2}. & \displaystyle\end{eqnarray}$$

Recalling that $NT_{L}=Fr^{-1}$, these scalings become

(5.4a,b)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}\sim Fr^{-1},\quad Fr\sim \left(\frac{L_{E}}{L_{O}}\right)^{-2},\end{eqnarray}$$

and hence at $Fr=O(1)$, we have

(5.5)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}\sim \left(\frac{L_{E}}{L_{O}}\right)^{2}.\end{eqnarray}$$

We therefore recover the behaviour shown in figure 12(a). GV19 instead find that $\unicode[STIX]{x1D6E4}\sim Fr^{-1}$ for $Fr=O(1)$, but $Fr\sim (L_{E}/L_{O})^{-2}$ for $Fr<O(1)$. Indeed, the derivation of (5.3) relies only on assumptions that are also applicable in the low $Fr$ regime. We do not believe that combining these scalings is inconsistent, since both rely on the assumptions that the dominant time scale for density-related terms is $N^{-1}$ and the dominant time scale for kinetic energy dissipation is $T_{L}=E_{K}/\unicode[STIX]{x1D700}$, and at $Fr=O(1)$ both of these time scales may affect the dynamics. On the other hand in a flow dominated by internal gravity waves, it is likely that $N^{-1}$ is the dominant time scale for both the velocity and density fields. This is evident in figures 12(b) and 12(c), where the wave-forced cases show $\unicode[STIX]{x1D6E4}_{H}$ to be less dependent on $(L_{E}/L_{O})_{H}$, consistent with the low $Fr$ and high $L_{E}/L_{O}$ regime, and at least conceivably consistent with a flow dominated by internal waves, thus still strongly affected by stratification. The larger mean value of $L_{E}/L_{O}$ in cases R and P can be associated with larger buoyancy excursions, providing further evidence that the wave-forced cases exhibit more convective behaviour. It appears that forcing with vortical modes at the same rate of energy input produces turbulence that can (at least locally) access a higher Froude number regime than the internal wave forcing. The fact that this $Fr=O(1)$ regime does not show the same $L_{E}/L_{O}$ scaling as in GV19 may hinder the ability to infer Froude numbers from observations in this intermediate range. The frequent appearance of $O(1)$ values of $L_{E}/L_{O}$ in observations (e.g. Moum Reference Moum1996) motivates the need for further research into mixing for flows in which $Fr=O(1)$.

The appearance of a $Fr=O(1)$ scaling in case H could suggest that the difference in volume-averaged $\unicode[STIX]{x1D6E4}$ between the simulations is purely related to a difference in the average Froude number. However, figure 8(a) also shows a region at low values of $Fr_{H}$ where $\unicode[STIX]{x1D6E4}_{H}$ appears independent of $Fr_{H}$ and importantly takes a much lower value than in the wave-forced cases R and P. The results of Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) also show $\unicode[STIX]{x1D6E4}\approx 0.35$ for values of $Fr$ similar to cases R and P but using a forcing scheme with a greater similarity to the forcing scheme used in case H. We are therefore confident that the difference in the large-scale forcing is the primary contributor to the changes in $\unicode[STIX]{x1D6E4}$ between the simulations, rather than a simple dependence on $Fr$. Although all the simulations here exhibit $\unicode[STIX]{x1D6E4}_{H}{-}Fr_{H}$ scalings consistent with the regimes predicted by Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) and GV19, we believe that it is important to understand how different flows lead to scatter around these regimes. To that end, a better understanding of how the local correlations in figures 9 and 10 are distributed would be invaluable. It is not currently clear how the $Fr_{L}$ dependence of $\unicode[STIX]{x1D6E4}_{L}$ shown in figure 10 leads to a global $\unicode[STIX]{x1D6E4}$ that is independent of $Fr$ for $Fr\ll 1$.

Despite the significant differences in the mixing properties between the vortical mode and wave-forced simulations, spectral analysis reveals remarkable similarity between the energy spectra associated with each flow. At moderate wavenumbers of $O(10)$, each component of the energy spectra appears to follow universal scalings, with wavelet analysis revealing distinct vertical wavenumber scalings between the turbulent and quiescent regions. The emergence of an $m^{-3}$ scaling in the low-$Re_{b}$ portions of the domain is particularly of note, since it is consistent with the energy spectra being determined exclusively by the buoyancy frequency (as discussed for example in § 14.3 of Davidson Reference Davidson2013). Differences in the energy spectra between the three simulations are only noticeable at low wavenumbers associated with the large-scale flow, despite the contrasting mixing efficiencies for each simulation persisting throughout the domains. This emphasises the importance of understanding the larger-scale flow dynamics when inferring small-scale mixing from spectra. With regions of the flows producing an $m^{-3}$ scaling consistent with Billant & Chomaz (Reference Billant and Chomaz2001), but variations in $\unicode[STIX]{x1D6E4}$ associated with the various larger-scale forcing strategies, it appears that the particular form of the larger-scale forcing retains a fundamental imprint on the properties of the small-scale mixing. Therefore, it is at least plausible that mixing events in the ocean could be sensitive to the particular form of the large-scale energy injection, suggesting that generic, ‘unified’ arguments such as those presented by Kunze (Reference Kunze2018) should be treated with caution.

Our results highlight a significant challenge in the measurement and parameterisation of turbulent mixing in the ocean. The turbulent flux coefficient $\unicode[STIX]{x1D6E4}$, commonly used to infer mixing rates from the TKE dissipation rate, varies by over 30 % depending on the nature of the larger-scale flow, although local estimates of $\unicode[STIX]{x1D6E4}$ can of course vary by much more. Mixing generated by internal gravity waves results in $\unicode[STIX]{x1D6E4}\approx 0.5$, consistent with recent studies of decaying stratified turbulence (de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019) but above results from numerical studies forcing turbulence through purely horizontal motion (Maffioli et al. Reference Maffioli, Brethouwer and Lindborg2016). This is also far higher than the value 0.2 from Osborn (Reference Osborn1980) typically used in observational studies, and also observed in simulations of statistically steady forced linearly sheared stratified turbulence with high dynamic range (Portwood et al. Reference Portwood, de Bruyn Kops and Caulfield2019). We conjecture that this difference may be associated with the mixing being more appropriately characterised as being ‘convectively driven’ rather than ‘shear instability driven’. This is consistent with previous studies (e.g. Davies Wykes, Hughes & Dalziel Reference Davies Wykes, Hughes and Dalziel2015) that show purely buoyancy-driven flows with non-monotonic stratification can achieve very high values of mixing efficiency. Mixing via shear instabilities often also occurs through a secondary convective instability arising due to the roll-up of the density field in a Kelvin–Helmholtz billow. The larger values of $L_{E}/L_{O}$ in our wave-forced simulations suggest the overturns are larger than from such shear-driven flows, consistent with the idea that the flow is ‘convectively driven’ by ‘breaking’ waves.

When interpreting our results in the context of ocean mixing, some caveats must be addressed. As mentioned in § 1, a significant fraction of mixing in the ocean occurs in surface and bottom boundary layers. The physics determining mixing efficiency in these regions is rather different from the ocean interior, with wind-driven shear and tidal flows as examples of important drivers of turbulence and mixing (Thorpe Reference Thorpe2005).

Furthermore, our simulations are performed with a molecular Prandtl number of 1 for computational efficiency, rather than a typical oceanic value of 7 for a thermally stratified region. Numerical studies of shear instabilities (Smyth, Moum & Caldwell Reference Smyth, Moum and Caldwell2001; Salehipour, Peltier & Mashayek Reference Salehipour, Peltier and Mashayek2015) have shown that higher Prandtl numbers can lead to a significant decrease in the mixing efficiency. This factor may bring our results closer to the value used in oceanographic estimates, but it is unclear how the differences in mixing efficiency between the simulations would change at higher $Pr$.

Another issue which needs to be considered is the potential ‘patchiness’ of mixing, with the turbulent mixing in the ocean exhibiting significant spatio-temporal variability. Observational studies that focus on small-scale mixing frequently isolate patches of turbulence for their measurements in intermittent oceanic flows (e.g. Moum Reference Moum1996). Both Smyth et al. (Reference Smyth, Moum and Caldwell2001) and Ijichi & Hibiya (Reference Ijichi and Hibiya2018) produce a $\unicode[STIX]{x1D6E4}\sim (L_{T}/L_{O})^{4/3}$ scaling from such patches, where $L_{T}$ is the Thorpe scale. The construction of this scaling (see Ijichi & Hibiya (Reference Ijichi and Hibiya2018) for further details) is fundamentally associated with two assumptions consistent with high $Fr$ dynamics, in particular that the characteristic time and length scales are determined from the turbulence properties alone, largely unaffected by the ambient stratification. The first assumption is that the turbulence is largely unaffected by stratification, and using a classical mixing length argument then leads to

(5.6a,b)$$\begin{eqnarray}L\sim \frac{E_{K}^{3/2}}{\unicode[STIX]{x1D700}};\quad \unicode[STIX]{x1D705}_{T}\equiv \frac{\unicode[STIX]{x1D6E4}\unicode[STIX]{x1D700}}{N^{2}}\sim E_{K}^{1/2}L.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D705}_{T}$ is the turbulent ‘eddy’ diffusivity of density, defined in terms of the ‘mixing length’ $L$ and the characteristic turbulent velocity scale $E_{K}^{1/2}$. The second assumption is that this mixing length $L\sim L_{T}$, leading to $\unicode[STIX]{x1D6E4}\propto (L_{T}/L_{O})^{4/3}$ in patches where the ambient stratification has little effect on the turbulent mixing properties, although it is of course possible that this scaling can be observed in situations where the underlying assumptions are no longer completely justified.

Indeed, the results presented here pose an interesting question of how best to model mixing in a spatio-temporally intermittent flow. The combination of our wavelet analysis and the work of Maffioli (Reference Maffioli2017) suggests that the $m^{-3}$ portion of the energy spectrum may be associated with larger scales and regions with smaller turbulent dissipation rates. The associated nonlinear buoyancy-dominated flow could be thought to act as a background from which the turbulent patches develop intermittently. Since high $Fr$ flows are associated with lower values of mixing efficiency, it is therefore important to quantify the relative contributions to mixing of these highly energetic isolated patches compared to the total background.

The main differences between our simulations, coinciding with the change in $\unicode[STIX]{x1D6E4}$, are an increase in the energy component of the vertical velocity and the available potential energy in our simulations. Despite these changes being most significant at large scales in our domain, validation of our results in the field would be difficult. Scales we refer to as large require high resolution equipment to resolve in the ocean: for example if we take a velocity scale of $U_{0}=10^{-2}~\text{m}~\text{s}^{-1}$ and a buoyancy frequency $N_{0}=10^{-2}~\text{s}^{-1}$, then our domain length will be less than 10 m given the $Re$, $Ri_{0}$ values we have chosen. An investigation of high-resolution measurements of vertical velocity and density fluctuations that coincide with the appearance of an $E\sim m^{-3}$ vertical wavenumber spectrum would be invaluable for determining the nature of flows at these scales. A better understanding of which flows are dominated by convective wave breaking and which mix through primarily horizontal motion or even shear instabilities, would then allow us to constrain mixing estimates better and improve our predictions of spatial patterns in diapycnal transport.