1. Introduction
In a density-stratified fluid, turbulence enhances the rate at which scalars are irreversibly diffused throughout the flow, a process of great interest in a variety of geophysical, environmental and industrial settings (e.g. Fernando Reference Fernando1991). Of particular importance is characterizing the role of turbulence in the vertical transport of heat within the ocean, a crucial mechanism for driving the required upwelling of cold bottom waters to maintain the ocean's vertical stratification profile and to complete global circulation currents (Wunsch & Ferrari Reference Wunsch and Ferrari2004). Turbulence in the ocean generates dynamically relevant motions of the order of millimetres, which cannot currently be resolved in numerical circulation models and must therefore be parametrized, with the choice of parametrization found to influence future climate projections strongly (Whalen et al. Reference Whalen, de Lavergne, Garabato, Klymak, Mackinnon and Sheen2020). Considerable observational, numerical and theoretical work has thus been focused on developing more accurate and universal mixing models which account for the wide range of turbulent processes observed in different flow regimes within the ocean (Caulfield Reference Caulfield2020).
The rate at which turbulence mixes a non-uniform density field is often defined in terms of an appropriately averaged vertical density flux 
$B \equiv \langle \rho 'w'\rangle$, where 
$\rho '$ and 
$w'$ denote fluctuations in density and vertical velocity away from the mean flow, respectively. If 
$B$ is to be used as a robust indicator of irreversible mixing, it is critical that measurements of 
$B$ are averaged over sufficiently large spatial volumes or time intervals to isolate irreversible diffusive processes from reversible stirring motions (Villermaux Reference Villermaux2019). Stirring, occurring on relatively large scales, may be thought of as the adiabatic rearrangement of fluid parcels of different density induced by the underlying turbulence, which in principle is reversible. Hence, a pointwise measurement of 
$B$ would not be a sufficient indicator that irreversible mixing had occurred, as the sign of 
$B$ could subsequently switch direction yielding a net flux of zero. Thus, we here use the term mixing to refer specifically to the diffusive transport of density across gradients that have been enhanced by such macroscopic stirring motions, irreversibly leading the system towards a state of greater homogenization. To isolate only irreversible contributions to mixing, Lorenz (Reference Lorenz1955) introduced the concept of an available potential energy (APE). The APE quantifies the difference between a system's current potential energy and its minimum background potential energy (BPE) that could be achieved if fluid parcels were adiabatically sorted into their most stable configuration. For a Boussinesq fluid, Winters et al. (Reference Winters, Lombard, Riley and D'Asaro1995) demonstrated that irreversible mixing may be described as the conversion of APE into BPE, with a system's BPE increasing in time as it homogenizes. Generalizing this mixing framework to compressible flows, Tailleux (Reference Tailleux2009) argued that the mixing of a thermally stratified fluid should most rigorously be defined as the conversion of APE into internal energy, which in the Boussinesq limit then exactly matches the generation of BPE.
Given a variety of sampling limitations involved with collecting turbulence data within the ocean, it is exceedingly difficult to perform the averaging required to extract the irreversible component of density fluxes from direct observational measurements of 
$B$. Therefore, a number of indirect methods have been proposed that infer such fluxes from more readily available quantities, which may be computed locally (Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018). Two such quantities, associated with what is conventionally referred to as the turbulent microstructure, are the dissipation rates of kinetic energy 
$\epsilon$ and scalar variance 
$\chi$,
\begin{equation} \epsilon=\frac{\nu}{2}\left(\frac{\partial u'_{i}}{\partial x_{j}}+\frac{\partial u'_{j}}{\partial x_{i}}\right)^{2}, \quad \chi=\frac{g^{2}\kappa}{\rho_{0}^{2}N^{2}}\left(\frac{\partial\rho'}{\partial x_{i}}\frac{\partial\rho'}{\partial x_{i}}\right), \end{equation}
representing the rates at which viscosity 
$\nu$ and molecular diffusivity 
$\kappa$ smooth gradients in the turbulent velocity 
$\boldsymbol {u'}$ and density 
$\rho '$ fields, respectively. In (1.1a,b), 
$g$ denotes the gravitational acceleration, 
$\rho _0$ a reference background density and 
$N=\sqrt {(-g/\rho _{0})\partial \bar {\rho }/\partial z}$ the buoyancy frequency, defined by an appropriately averaged ambient density gradient 
$\partial \bar {\rho } / \partial z$ against which the turbulence acts. The quantities 
$\epsilon$ and 
$\chi$ are intimately related to the irreversible processes associated with the conversion of kinetic energy and available potential energy into internal energy, respectively, as is further described by Caulfield (Reference Caulfield2021). In particular, for the class of direct numerical simulation considered here, characterized by an imposed uniform background stratification 
$N^2_0$, Howland, Taylor & Caulfield (Reference Howland, Taylor and Caulfield2021) demonstrated that 
$\chi$ computed using 
$N^2=N_0^2$ in (1.1a,b) provides an excellent approximation for the destruction rate of APE and is therefore a good measure of local irreversible mixing.
As discussed by Ivey, Bluteau & Jones (Reference Ivey, Bluteau and Jones2018), 
$\chi$ also arguably provides the most robust method for estimating irreversible mixing from oceanographic measurements, since 
$\langle B\rangle \simeq \langle \chi \rangle$ in the steady state provided that averaging is performed over sufficiently long times and large volumes so that reversible processes and transport terms are negligible (Osborn & Cox Reference Osborn and Cox1972). Importantly, 
$\chi$ is both directly proportional to the scalar diffusivity 
$\kappa$ and sign-definite, providing a robust local measure of the irreversible fluxes associated with molecular diffusion, which does not require the averaging of the density flux 
$B$ needed to filter our reversible local stirring motions in the turbulent flow. Due to a scarcity of 
$\chi$ measurements, however, 
$\epsilon$ is more commonly used to infer mixing following the method of Osborn (Reference Osborn1980), which requires the introduction of a flux coefficient 
$\varGamma \equiv \chi / \epsilon$ to prescribe the fraction of turbulent kinetic energy that leads to irreversible mixing, as opposed to being directly dissipated by viscosity. A constant value 
$\varGamma = 0.2$ is commonly assumed when estimating global patterns of oceanic mixing (MacKinnon et al. Reference MacKinnon2017), which has been found to be in agreement with tracer release experiments (Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018). However, there is significant evidence suggesting that 
$\varGamma$ varies appreciably in different flow regimes (Caulfield Reference Caulfield2021) and so a clear physical picture has not yet emerged explaining why 
$\varGamma =0.2$ is a reasonable assumption.
In the absence of measurements of 
$\epsilon$ or 
$\chi$, mixing locations are primarily inferred from the presence of unstable overturns in vertical density profiles, as proposed by Thorpe (Reference Thorpe1977). Assuming that the vertical extent of an overturn is correlated with the Ozmidov length 
$L_{O}=\sqrt {\epsilon /N^{3}}$, 
$\epsilon$ may be inferred from the measurement of overturns which can then be converted into a flux via 
$\varGamma$. However, this assumed correlation between the vertical overturning scale and 
$L_O$ is not always robust, as has recently been discussed, for example, by Ivey et al. (Reference Ivey, Bluteau and Jones2018), Ijichi et al. (Reference Ijichi, St. Laurent, Polzin and Toole2020) and Mashayek, Caulfield & Alford (Reference Mashayek, Caulfield and Alford2021). Using a forced direct numerical simulation (DNS) similar to that considered here, Taylor et al. (Reference Taylor, de Bruyn Kops, Caulfield and Linden2019) quantified the errors associated with the indirect flux estimates of Osborn & Cox (Reference Osborn and Cox1972), Osborn (Reference Osborn1980) and Thorpe (Reference Thorpe1977) by sparsely sampling vertical profiles of the computational domain to mimic oceanographic measurements.
Spatio-temporal intermittency in stratified turbulence greatly reduces the applicability of classical turbulence modelling assumptions, including the common assumption of log-normal distributions for 
$\epsilon$ and 
$\chi$ (de Bruyn Kops Reference de Bruyn Kops2015). Cael & Mashayek (Reference Cael and Mashayek2021) found that global ocean measurements of 
$\epsilon$ were not well approximated by an assumed log-normal distribution but instead had a skewed right tail, indicating that a small number of extreme events dominated the bulk statistics. By considering local correlations between direct ocean measurements of 
$\epsilon$ and 
$\chi$, Couchman et al. (Reference Couchman, Wynne-Cattanach, Alford, Caulfield, Kerswell, MacKinnon and Voet2021) further emphasized the importance of extreme events, finding that while the majority of the sampled domain was characterized by the canonical flux coefficient 
$\varGamma = 0.2$, isolated mixing events containing the largest 
$\chi$ were not reflected by a corresponding local increase in 
$\epsilon$, yielding a dramatic increase in 
$\varGamma$.
Vertical layering is also known to be a canonical feature of stratified turbulent flows, with the density field often forming ‘staircases’ of deep, relatively well-mixed layers separated by thin interfaces with strong gradients (Caulfield Reference Caulfield2021). For sufficiently stratified environments, vertical shearing induced by the decoupling of horizontal and vertical motions in such a layered structure becomes an important mechanism for triggering instability and the ensuing generation of turbulence (Lilly Reference Lilly1983; Billant & Chomaz Reference Billant and Chomaz2001). Parametrizations of mixing based on simple domain averages are thus unlikely to be accurate as rare extreme events and spatial heterogeneity within the flow will be missed, a potential cause of the highly scattered mixing statistics currently reported throughout the literature (Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018).
In an attempt to classify such intermittency in an automatic, yet robust and interpretable manner, Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) devised an algorithm for splitting a snapshot from a forced DNS into three dynamically distinct regions: quiescent regions, intermittent layers and turbulent patches. These regions were distinguished by an increasing degree of local overturning, as determined by computing the fraction of unstable points 
$\partial \rho / \partial z > 0$ within an extended neighbourhood. Local overturning fractions and dissipation rates 
$\epsilon$ were found to be strongly correlated, in agreement with the arguments of Thorpe (Reference Thorpe1977). For the relatively large filter sizes used to segment the domain, of the order of a buoyancy length 
$L_B = 2{\rm \pi} u_h /N$ (where 
$u_h$ denotes a characteristic horizontal velocity scale), distributions of 
$\chi$ associated with each region were also found to be correlated with 
$\epsilon$, although the finer spatial distributions of 
$\epsilon$ and 
$\chi$ within each region, and the resulting flux coefficient 
$\varGamma$, were not considered.
Motivated by the automated flow segmentation of Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) in terms of unstable local density gradients 
$\partial \rho / \partial z > 0$, and the observation of Couchman et al. (Reference Couchman, Wynne-Cattanach, Alford, Caulfield, Kerswell, MacKinnon and Voet2021) that, within the ocean, extreme events in 
$\chi$ are not necessarily correlated with those in 
$\epsilon$, we here analyse spatial mixing distributions within a computational domain by considering local correlations among 
$\epsilon$, 
$\chi$ and 
$\partial \rho / \partial z$. In particular, we wish to probe whether overturning alone provides a robust indicator for local mixing, or if significant mixing as revealed by 
$\chi$ might occur in other regions that would seem inconspicuous based on consideration of only 
$\epsilon$ or 
$\partial \rho / \partial z$.
In line with the previous investigations of Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) and Taylor et al. (Reference Taylor, de Bruyn Kops, Caulfield and Linden2019), we consider a forced DNS of stratified turbulence using the methodologies presented by Almalkie & de Bruyn Kops (Reference Almalkie and de Bruyn Kops2012). In § 2, we summarize the DNS dataset considered here and highlight the presence of a (previously unreported) robust vertically aligned vortex generated by the forcing scheme, that injects energy into the domain at large scales and induces vertical layering in the surrounding flow. In § 3, we then consider pointwise correlations among 
$\partial \rho / \partial z$, 
$\epsilon$, 
$\chi$ and the flux coefficient 
$\varGamma$, which suggests that mixing occurs not only in overturning regions, but also in areas of local static stability. In § 4, we move beyond pointwise statistics to consider extended mixing structures within the flow, highlighting two ways in which local static instability in the density gradient fails to be a sufficient indicator of mixing: within the vortex, a lateral density gradient is correlated with the majority of 
$\chi$, and outside the vortex, extreme values of 
$\chi$ are localized to relatively ‘sharp’ stable density interfaces at the bounding edges between overturning layers. Finally, in § 5, we summarize our results and discuss implications for parametrizing turbulent mixing within the ocean.
2. Summary of DNS dataset
We consider a statistically steady, forced DNS of stratified turbulence from the simulation campaign originally reported by Almalkie & de Bruyn Kops (Reference Almalkie and de Bruyn Kops2012), and subsequently analysed by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) and Taylor et al. (Reference Taylor, de Bruyn Kops, Caulfield and Linden2019). Using a characteristic root-mean-square horizontal velocity scale 
$u_h$, length scale 
$L$ and background buoyancy frequency 
$N_0$, the non-hydrostatic Boussinesq approximation of the Navier–Stokes equations may be written in the following dimensionless form:
\begin{gather}\begin{aligned}
\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0,
\quad \frac{\partial\boldsymbol{u}}{\partial
t}&+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}
=-\left(\frac{2{\rm \pi}}{Fr}\right)^{2}\rho\hat{\boldsymbol{z}}-\boldsymbol{\nabla}
p+\frac{\nabla^{2}\boldsymbol{u}}{Re}+\mathcal{F}, \\
\frac{\partial\rho}{\partial
t}&+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho-w =\frac{\nabla^{2}\rho}{RePr}.\end{aligned} \end{gather} The governing equations (2.1a–c) are numerically integrated using a pseudospectral technique in a triply periodic domain, as detailed by Almalkie & de Bruyn Kops (Reference Almalkie and de Bruyn Kops2012). The dimensionless parameters governing the flow are the Prandtl number 
$Pr = \nu / \kappa$, Froude number 
$Fr_h = 2 {\rm \pi}u_h / (N_0 L)$ and Reynolds number 
$Re_h = u_h L / \nu$. The density field satisfies
\begin{equation} {\rho(\boldsymbol{x},t)=\rho_{0}(1 -N_{0}^{2}z/g)+\rho' (\boldsymbol{x},t)}, \end{equation}
where 
$\rho _{0}(1 -N_{0}^{2}z/g)$ defines a time-independent, linear background density gradient characterized by a reference density 
$\rho _0$ and an imposed constant background buoyancy frequency 
$N_0$. Density perturbations 
$\rho '$ away from this linear background state satisfy the periodic boundary conditions and are used to compute 
$\chi$ in (1.1a,b). The imposed constant background buoyancy frequency 
$N_0$ is used as the characteristic ‘appropriately averaged’ buoyancy frequency 
$N$ required to compute 
$\chi$ in (1.1a,b), as is widely considered the natural choice when quantifying irreversible mixing in numerical simulations with an imposed background stratification (see e.g. Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Maffioli, Brethouwer & Lindborg Reference Maffioli, Brethouwer and Lindborg2016; Garanaik & Venayagamoorthy Reference Garanaik and Venayagamoorthy2019; Portwood, de Bruyn Kops & Caulfield Reference Portwood, de Bruyn Kops and Caulfield2019). By explicitly computing the APE of a triply periodic domain with an imposed uniform background stratification 
$N_0$, Howland et al. (Reference Howland, Taylor and Caulfield2021) confirmed that normalizing 
$\chi$ by 
$N_0$ indeed provides an excellent approximation to the true irreversible mixing rate as computed through changes in the system's APE. The forcing term 
$\mathcal {F}$ is governed by the deterministic scheme denoted ‘Rf’ by Rao & de Bruyn Kops (Reference Rao and de Bruyn Kops2011), which adds energy to horizontal motions larger than 1/8th of the horizontal box size so as to match a target kinetic energy spectrum at small wavenumbers. We consider a simulation characterized by 
$Pr = 1$, 
$Fr_h = 2.23$ and 
$Re_h = 1271$, in a domain of size 
$2 {\rm \pi}\times 2 {\rm \pi}\times {\rm \pi}$ with 
$4096 \times 4096 \times 2048$ grid points, resulting in a grid spacing of 
$\varDelta \approx L_K/2$, with 
$L_K$ denoting the Kolmogorov length scale. For reference, the characteristic buoyancy Reynolds number of the simulation is 
$Re_{b}=\langle \epsilon \rangle /\nu N_{0}^{2}=50$. We consider a single snapshot of the flow in time and all figures are displayed on grids that have been sparsed by a factor of eight in each dimension.
The main characteristics of the dataset are summarized in figure 1. The vertically averaged dissipation rate 
$\epsilon$ (figure 1a) reveals a dominant patch of elevated turbulence that is generated by a large-scale vertically aligned vortex in the velocity field, rotating counterclockwise (figure 1b). Radial averages, centred on the vortex, of the angular velocity 
$u_\theta$ and vertical component of vorticity 
$\omega _z = \partial v/\partial x-\partial u/\partial y$ are plotted in figure 1(c), indicating a Rankine-type vortex that is approximately characterized by rigid-body rotation at small radii 
$r$ from the vortex core, followed by a transition to roughly irrotational flow at larger 
$r$. It is important to note that such a description characterizes the radially averaged flow, and that smaller-scale vortical motions will still certainly be present in the turbulent patch surrounding the vortex. A series of horizontal currents travelling in alternating directions are found to emanate from the vortex, characterized by a vertical scale of the order of a buoyancy length 
$L_B$. In figure 1(d), we plot the vertical slice of the density field that corresponds to the velocity field shown in figure 1(b), highlighting an analogous vertically layered structure outside of the vortex, with relatively well-mixed density layers separated by sharp, stable interfaces. The approximate locations of these density interfaces (delineated by green contours) correspond to minima in the histogram of 
$\rho$ (figure 1e), which highlights a strong perturbation of the density field away from its uniform background gradient. Superimposing these density contours on the velocity field in figure 1(b) highlights that the sheared interfaces in the velocity field are strongly correlated with the stable interfaces in the density field. This correlation is further demonstrated in supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.253, where rotations of the slices in figure 1(b,d) around the centre of the vortex are shown. In § 4, we demonstrate that these interfaces, characterized by both high shear and a strong statically stable density gradient, are critically important for the mixing generated outside of the vortex.
Figure 1. Large-scale characteristics of the velocity and density fields within the computational domain. (a) Vertical average of the dissipation rate 
$\epsilon$ at gridpoint 
$(x,y)$, normalized by the domain average. The whole domain is shown, with gridpoints sparsed by a factor of eight in each dimension. The dashed circle highlights a region of elevated 
$\epsilon$, coinciding with a vortex in the velocity field, as is further examined in figure 3. A vertical slice of the domain at 
$y=200$ is considered in figure 4. (b) Horizontal velocity normal to the plane for a vertical slice at constant 
$y$ passing through the centre of the dashed circle in panel (a), revealing a vertically aligned vortex rotating counterclockwise. The grid has been shifted in 
$x$ relative to panel (a) so as to centre the vortex. The buoyancy length 
$L_B$ and Taylor length 
$L_T \approx 25 L_K$ (where 
$L_K$ is the Kolmogorov length) are marked for reference, as were the filter sizes used in the segmentation analysis of Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016). Green contours mark stable interfaces in the density field, as shown in panel (d). (c) Radially averaged angular velocity 
$u_\theta$ and vertical component of vorticity 
$\omega _z$ as a function of the distance 
$r$ from the centre of the vortex. The dashed line at 
$r\approx 80$ corresponds to the dashed circle plotted in panel (a). (d) Density field corresponding to the vertical slice of velocity plotted in panel (b). The green lines in panels (b) and (d) illustrate contours at the minimum values of the histogram of the density field in panel (e), as are marked with vertical green lines, delineating the interfaces between relatively well-mixed density layers in the flow surrounding the vortex.
The spontaneous formation of a persistent vortex is a key, yet previously unreported feature of the forcing scheme of Rao & de Bruyn Kops (Reference Rao and de Bruyn Kops2011) used to generate statistically steady turbulence. In particular, the identification of the vortex in figure 1 provides insight into how the segmentation results of Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) (see their figure 2c), who used an identical forcing scheme, are related to the background flow field. Specifically, the roughly cylindrical patch of most vigorous turbulence detected by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), using a filter of size 
$L_B$, extends across the entire vertical domain and almost certainly corresponds to an analogous vortical structure in their DNS. Similarly, their ‘intermittent layers’ are primarily composed of horizontal offshoots from the central vertically aligned turbulent patch, and are shaped by a similar pattern to the sheared velocity interfaces observed in figure 1(b). As it is now evident that Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) have broadly identified such a vortex to be a turbulent hotspot, a goal of this study is to perform a finer analysis of mixing patterns both within and outside of the vortex to determine how patterns in the small-scale turbulent microstructure, as described by 
$\epsilon$ and 
$\chi$, are related to the larger-scale layered structure of the flow.
3. Pointwise statistics conditioned on local density gradient
Motivated by the flow segmentation of Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) in terms of the local fraction of overturning 
$\partial \rho / \partial z > 0$, we first consider how the pointwise distributions of 
$\epsilon$, 
$\chi$ and 
$\varGamma = \chi / \epsilon$ depend on the magnitude of 
$\partial \rho / \partial z$, for both statically stable and unstable points, as shown in figure 2. For illustration, in figure 2(a), we split the distribution of 
$\partial \rho / \partial z$ into three regions: two tails containing 
$10\,\%$ by volume of the most stable and unstable points (coloured blue and red, respectively), and the remaining 
$80\,\%$ of the intermediate values (green). For such a division, we then consider the distributions of 
$\epsilon$, 
$\chi$ and 
$\varGamma$ within each region, as shown in figure 2(b–d). Although the distribution characterizing the bulk of the domain (green) is centred around the canonical flux coefficient 
$\varGamma = 0.2$ (see figure 2d), such points contain the lowest 
$\chi$ (figure 2c) and are thus not of primary importance for the total mixing arising within the computational domain. Instead, it is the extreme tails of the 
$\partial \rho / \partial z$ distribution that must be considered, containing the most significant values of 
$\chi$. While both the blue and red tails contain elevated but similar distributions of 
$\epsilon$, they may be distinguished by their asymmetry in 
$\chi$; the stable tail (blue) contains disproportionately elevated 
$\chi$ as compared to 
$\epsilon$, and therefore some of the highest values of 
$\varGamma$ within the domain.
Figure 2. Pointwise mixing statistics of the DNS data, conditioned by the local vertical density gradient. (a) Histogram of the perturbed density gradient 
$\partial \rho ' / \partial z$ normalized by the magnitude of the imposed uniform background gradient, with values greater than one indicating local overturning. The distribution is split into three regions, by assigning a fixed volume (here 
$10\,\%$) to each tail. Panels (b–d) illustrate the distributions of 
$\epsilon$, 
$\chi$ and 
$\varGamma$, respectively, for the whole domain (black) and the subdomains encompassed by the coloured regions in panel (a). Circles mark the median values of each distribution, and the dashed line in panel (d) indicates the canonical flux coefficient 
$\varGamma = 0.2$ for reference. Panels (e–g) illustrate how the medians of the respective distributions in panels (b–d) vary with the tail volume selected in panel (a). The circles mark the medians for the segmented distributions shown in panels (b–d) for the case of 
$10\,\%$ tail volume. For panels (e) and (f), the fraction of each quantity (
$\epsilon$, 
$\chi$) contained within each tail relative to the entire domain is also indicated. The dashed diagonal lines mark what would be expected for a uniform distribution of each quantity throughout the domain.
In figure 2(e–g), we then analyse how the medians of the 
$\epsilon$, 
$\chi$ and 
$\varGamma$ distributions change as a function of the volume contained within the blue and red tails, and additionally plot the relative contributions of these tails to the domain total. Comparing the right panels of figures 2(e) and 2(f) reveals that 
$\chi$ is far more dominated by extreme events than 
$\epsilon$, in agreement with the analysis of oceanographic data by Couchman et al. (Reference Couchman, Wynne-Cattanach, Alford, Caulfield, Kerswell, MacKinnon and Voet2021). For instance, when each tail contains 
$10\,\%$ volume, the stable (blue) and unstable (red) tails each contain approximately 
$20\,\%$ of the total 
$\epsilon$ in the domain, but 
$45\,\%$ and 
$30\,\%$ of the total 
$\chi$, respectively. Furthermore, while the contributions to 
$\epsilon$ from both tails are roughly equal, the contribution to 
$\chi$ from the stable tail is always roughly 
$50\,\%$ greater than for the unstable tail. While 
$\varGamma =0.2$ may thus be a suitable approximation for the bulk of the domain, it may here not be relied upon for capturing the most extreme events in 
$\chi$, which dominate the bulk mixing statistics.
Additionally, the statistics in figure 2 suggest that local instability may not be a sufficient indicator for mixing, given the significance of the blue stable tail. However, we note that such a conclusion cannot definitively be drawn from the pointwise distributions of 
$\partial \rho / \partial z$, as such a distribution provides no information about the extended spatial environment around each point. For example, in regions of fully developed turbulence that might emerge after the collapse of a shear-induced billow, there is likely a random mixture of neighbouring unstable and stable points in close proximity (roughly a 
$50\,\%$ mixture as identified by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) in their most turbulent patches), and so points within the red and blue tails of figure 2(a) could be direct neighbours in space. Therefore, in § 4, we extend our pointwise analysis by identifying spatially extended and coherent stable regions, which appear to take the form of ‘interfaces’ with enhanced density gradients. We then assess the significance of these non-overturning structures to the overall mixing statistics.
4. Extended mixing structures
We now consider coherent spatial distributions of the microstructure quantities 
$\epsilon$ and 
$\chi$, and their relation to the large-scale flow patterns observed in figure 1, by focusing on mixing structures arising both within and outside of the vortex. We first perform a closer examination of mixing within the vortex, as shown in figure 3. Vertical averages of 
$\epsilon$, 
$\chi$, 
$w$ and 
$\rho '$ are plotted in figures 3(a)–3(d), respectively, highlighting clear differences in the spatial distributions of 
$\epsilon$ and 
$\chi$. Such differences are further illustrated in figures 3(e)–3(h), which show the respective radial distributions of the azimuthal velocity 
$u_\theta$, 
$\epsilon$, 
$\chi$ and 
$\varGamma$ with respect to the vortex core. These radial distributions illustrate that the inner section of the vortex, characterized by roughly rigid-body rotation, is well mixed and contains the largest values of 
$\epsilon$ despite having minimal scalar diffusion rates 
$\chi$. This observation is consistent with the density field shown in figure 1(d), where initially horizontal contours of constant density (green lines) are strongly deflected towards the vertical before reaching the centre of the vortex, resulting in a vertically extended core of roughly constant density (seen predominantly in the vertical interval 
$25\lesssim z\lesssim 175$). Conversely, the majority of 
$\chi$ is found outside the core at radii where the angular velocity begins to decay, and is distributed in a roughly spiral pattern (figure 3b). Examination of the vertically averaged perturbed density field 
$\rho '$ (figure 3d) reveals the presence of a strong lateral density gradient, induced by the alternating upwelling of dense fluid and downwelling of lighter fluid within the vortex as a function of 
$r$. Superimposing the position of this lateral gradient (grey) onto the distribution of 
$\chi$ in figure 3(b) reveals that this gradient is strongly correlated with the spiral distribution of the most intense 
$\chi$. While the vortex was identified by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) to be a patch of vigorous turbulence with elevated 
$\epsilon$ due to its generation of significant local vertical overturning, our analysis suggests that much of the mixing within the vortex, as quantified by 
$\chi$, instead results from diffusion across a strong lateral gradient in the perturbed density field.
Figure 3. Mixing patterns within the vortex. Vertical averages of the (a) dissipation rate of kinetic energy 
$\epsilon$, (b) dissipation rate of scalar variance 
$\chi$, (c) vertical velocity 
$w$ and (d) density perturbations 
$\rho '$, for the vortex region delineated in figure 1(a). The outer green circle in panels (a–d) coincides with the black circle in figure 1(a). The grey curve in panel (d) delineates the sharp lateral gradient separating regions of positive and negative 
$\rho '$ and is found to be correlated with the spiral distribution of 
$\chi$ observed in panel (b). Radial averages of (e) the azimuthal velocity 
$u_\theta$, (f) 
$\epsilon$, (g) 
$\chi$ and (h) 
$\varGamma$, as a function of the distance 
$r$ from the centre of the vortex, with shading denoting the standard deviation around the radial mean. In panel (h), 
$\langle \varGamma \rangle _{r}=\langle \chi \rangle _{r}/\langle \epsilon \rangle _{r}$ is the ratio of the red lines in panels (f) and (g). The vertical green dashed lines in panels (e–h) mark the radial locations of the maximum and first zero of the radially averaged azimuthal velocity, and correspond to the radii of the green dashed circles in panels (a–d).
Outside of the vortex, the vertical homogeneity of the velocity and density fields collapses, forming a vertically layered structure. In figure 4, we consider a vertical 
$(x,z)$ slice of the domain at position 
$y=200$ in figure 1(a), to understand how this large-scale layering pattern gives rise to mixing at the microscale. Motivated by the significance of the stable tail (blue) in the pointwise distribution of 
$\partial \rho / \partial z$ in figure 2(a), and the observation of horizontally extended stable filaments of 
$\partial \rho / \partial z$ in figure 4(a), we examine whether such structures contribute substantially to mixing in the layered flow surrounding the vortex. To isolate these stable filaments, we apply a Gaussian filter to the density field with standard deviation 
$\sigma \approx 6 L_K$ (corresponding to two grid points in figure 4), where 
$L_K$ denotes the Kolmogorov length scale. The intent of such a filter is to isolate spatially coherent stable structures from patchy overturning regions that would contain a random assortment of stable and unstable neighbouring points. We note that our filter length is of the order of 
$10 L_K$ as suggested by Kuo & Corrsin (Reference Kuo and Corrsin1971) for removing internal intermittency. Further, it is significantly finer than the Taylor length 
$L_T \approx 25 L_K$, which was the smallest filter size considered by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) in their identification of ‘intermittent layers’, allowing us to examine the importance of finer-scale structures within the flow.
Figure 4. Characteristics of the layered flow outside the vortex. Vertical slices at 
$y=200$ in figure 1(a) of: (a) the density field 
$\partial \rho / \partial z$ normalized by the magnitude of the imposed background gradient; (b) the result of applying a Gaussian filter with a standard deviation of two gridpoints (
$\approx 6L_K$) to the density field in panel (a); (c) extracted stable filaments from panel (b) obtained by retaining points in the bottom (most stable) 
$15\,\%$ of the filtered density distribution; (d) the local fraction of unstable overturned points computed with a filter size corresponding to the Taylor length (
$L_T \approx 25L_K$), following the method of Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016); (e,f) the 
$x$ and 
$y$ horizontal components of velocity, respectively; and (g,h) the dissipation rates of kinetic energy 
$\epsilon$ and scalar variance 
$\chi$, respectively. The stable filaments from panel (c) are overlaid on panels (d–f) for reference. The contours coloured green in panels (c–f) and white in panels (g,h) correspond to those plotted in figure 1(b,d), marking the stable interfaces separating relatively well-mixed density layers. Green crosses in panels (g,h) indicate examples of locations where 
$\chi$ is locally high, due to the presence of a stable density interface, without a corresponding increase in local 
$\epsilon$.
Having filtered the density field (figure 4b), we then extract the most stable density structures by considering points in the bottom (most stable) 
$q$ percent of the filtered 
$\partial \rho / \partial z$ distribution. For illustration, we here extract structures comprising the most stable 
$q=15\,\%$ of points (figure 4c), and in the Appendix demonstrate the effect of changing this percentage. The green contours from figure 1(b,d) are overlaid on figure 4(c), demonstrating that the extracted filaments correspond to segments of the sharp interfaces separating relatively well-mixed layers in the density field. Importantly, figure 4(d) highlights that the concentration of locally overturned points (the segmentation indicator used by Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) is greatest in the regions between these stable interfaces. In figure 4(e,f), we again highlight that these stable interfaces are also roughly correlated with regions of high vertical shear in the layered velocity field, as are generated by the vortex. Corresponding slices of the dissipation rates of kinetic energy 
$\epsilon$ and scalar variance 
$\chi$ are shown in figure 4(g,h), respectively. The spatial distribution of 
$\epsilon$ is seen to be much more diffuse than that of 
$\chi$, with extreme values of 
$\chi$ being primarily concentrated within thin filamentary structures such as those identified in figure 4(c). Crucially, there are many examples of locations in the flow (see green crosses, figure 4g,h) where the stable interfaces identified in figure 4(c) contain highly elevated local signatures of 
$\chi$ without a proportional local increase in 
$\epsilon$. Such an observation thus raises the question as to whether these stable interfaces contribute significantly to the total mixing within the domain, in addition to the mixing occurring in more conventionally studied isotropic overturning regions (such as the large overturn located in the vicinity of 
$(x,y)=(450,125)$ in figure 4).
We address the question of whether the identified stable filaments contribute substantially to the total mixing occurring within the domain in figure 5, where we consider the relative contributions of both the vortex and the isolated stable interfaces to the domain totals of 
$\epsilon$ and 
$\chi$. In agreement with Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), despite occupying less than 
$10\,\%$ of the domain volume, the vortex contributes approximately a third of the entire domain's 
$\epsilon$ and 
$\chi$ (red bars, figure 5a). Outside of the vortex, however, it is the stable interfaces that play a key role in the overall mixing, contributing
\begin{equation} \frac{\text{(% in interface)} \cap \text{(% outside vortex)}}{(\text{% outside vortex})} =\frac{26\,\%}{26\,\%+40\,\%}=39\,\% \end{equation}
of the total 
$\chi$ outside the vortex, despite appearing unremarkable based on their much smaller contribution to 
$\epsilon$ (
$11\,\%/(11\,\%+55\,\%)=17\,\%$). Figure 4(d) thus highlights a key conclusion of this study: while the concentration of overturned points is most prevalent within the well-mixed density layers, relatively thin stable interfaces between such relatively deep layers, which are also correlated with high vertical shear, yield a crucial component of the bulk scalar mixing rate 
$\chi$. In particular, figure 5(b,c) highlights that while these interfaces may be strongly distinguished by their distributions of 
$\chi$, where the median values differ by almost an order of magnitude, they are virtually indistinguishable based on their distributions of 
$\epsilon$. This mismatch between the spatial distributions of 
$\epsilon$ and 
$\chi$ results in significantly elevated 
$\varGamma$ within the interfaces, well above the canonical value 
$\varGamma = 0.2$ (figure 5d). It thus appears crucial to consider the independent information provided by the distributions of 
$\epsilon$ and 
$\chi$ within a domain when quantifying mixing, particularly for identifying the locations of the most extreme scalar mixing events.
Figure 5. Mixing contributions of the vortex and stable interfaces. (a) Fractional contributions to the domain total of the dissipation rates of kinetic energy 
$\epsilon$ and scalar variance 
$\chi$, from within the vortex (red), stable interfaces outside the vortex (dark blue), as illustrated in figure 4(c), and the rest of the domain outside both the vortex and interfaces (light blue). Histograms of (b) 
$\epsilon$, (c) 
$\chi$ and (d) 
$\varGamma$ outside the vortex (black), further split according to whether points are contained within a stable interface (dark blue) or not (light blue). Dashed lines indicate median values of the respective distributions.
5. Discussion
We have considered local correlations between the vertical density gradient 
$\partial \rho / \partial z$ and the dissipation rates of kinetic energy 
$\epsilon$ and scalar variance 
$\chi$ to characterize the spatial distributions of mixing within a forced direct numerical simulation of density-stratified turbulence. The forcing scheme is found to generate a vertically aligned vortex within the domain, largely explaining the concentrated ‘patch’ region of vigorous turbulence reported by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016). Outside of the vortex, the flow is characterized by a layered density profile, with thin, highly stable interfaces separating relatively well-mixed layers. While a mixing analysis based solely on the identification of local overturning would deem the well-mixed layers to be of primary importance, as in the identification of ‘intermittent layers’ with elevated 
$\epsilon$ by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), we have demonstrated that a significant fraction of 
$\chi$ is localized to the edges of such layers, within the stable intervening interfaces. Notably, these interfaces appear unremarkable if looking at 
$\epsilon$ alone (see figure 5b), emphasizing the importance of 
$\chi$ as an independent indicator of local mixing. A number of other studies have also highlighted that significant mixing rates may be found in regions devoid of local overturning, emphasizing the importance of considering other mixing mechanisms present within stratified flows. For instance, by considering a different class of forced direct numerical simulations to those analysed here, Basak & Sarkar (Reference Basak and Sarkar2006) demonstrated that horizontal shear is able to generate a complex pattern of vorticies which efficiently mix the density field without local overturning. A striking experimental demonstration of dye being transported across stationary, highly stable density interfaces has been demonstrated by Oglethorpe, Caulfield & Woods (Reference Oglethorpe, Caulfield and Woods2013), where a ‘scouring’ rather than overturning dynamic generates the mixing. As the flux coefficient 
$\varGamma = \chi / \epsilon$ has been found to strongly depend on the time history of a turbulent event (Mashayek et al. Reference Mashayek, Caulfield and Alford2021), it would be instructive to now consider the time evolution and formation of the stable interfaces identified in our study, characterized by strongly elevated 
$\varGamma$. For instance, as the density interfaces are correlated with regions of high vertical shear, it is conceivable that they might be remnants of the previous collapse of shear-induced billows that are now only visible in signals of 
$\chi$ but not 
$\epsilon$, as coined ‘fossil turbulence’ by Nasmyth (Reference Nasmyth1970).
Our findings have two potential implications for the parametrization of ocean mixing. First, our analysis highlights the importance of adequately sampling rare, yet extreme mixing events in a turbulent flow, as was also recently discussed by Cael & Mashayek (Reference Cael and Mashayek2021). In agreement with the analysis of oceanographic data by Couchman et al. (Reference Couchman, Wynne-Cattanach, Alford, Caulfield, Kerswell, MacKinnon and Voet2021), figure 2(d) demonstrates that although the majority of the domain indeed appears to be well characterized by the canonical flux coefficient 
$\varGamma = 0.2$, significantly elevated 
$\varGamma$ is associated with the most extreme events in 
$\chi$, events that are not reflected by a corresponding local increase in 
$\epsilon$. Given the current relative sparsity of measurements within the ocean, mixing parametrizations may thus be biased towards the most commonly measured events, which are not necessarily the most significant. Second, even with perfect sampling, different proxies for mixing are likely to yield contrasting predictions for the amount and spatial distribution of mixing within the highly anisotropic layered flow considered here. For example, if measurements of 
$\chi$ were not available, the stable filaments at the edges of the overturning layers (figure 4d) would appear unremarkable, as they appear locally quiescent based on their density gradient and are not correlated with any discernible increase in 
$\epsilon$. Further, given the strong spatial variability of 
$\varGamma$ within the vertically layered flow (see figure 5d), it is unclear what value of 
$\varGamma$ should be used in the method of Osborn (Reference Osborn1980) if trying to infer a flux from values of 
$\epsilon$ measured directly by a microstructure profiler or derived from a Thorpe overturning analysis.
As discussed by Caulfield (Reference Caulfield2021), an accurate parametrization of the flux coefficient 
$\varGamma$ is likely to depend on multiple dimensionless groups characterizing the underlying flow, such as the buoyancy Reynolds number 
$Re_b$, Froude number 
$Fr$ and Prandtl number 
$Pr$. For instance, DNS studies have demonstrated that bulk-averages of 
$\varGamma$ decrease with increasing 
$Pr$ (Salehipour, Peltier & Mashayek Reference Salehipour, Peltier and Mashayek2015) and decreasing 
$Fr$ (Maffioli et al. Reference Maffioli, Brethouwer and Lindborg2016). A promising future direction of inquiry would be to try and rationalize such variations in 
$\varGamma$ in terms of differences in the prevalence and structure of smaller-scale extreme events within the flow, such as analysing changes in the morphology of the stable filaments considered here. It would also be instructive to extend our analysis to simulations of decaying turbulence which also develop layered structures (de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019) to establish whether the spatial distribution of mixing events observed here changes significantly in forced versus unforced scenarios.
Finally, following Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016) and typical oceanographic measurements, we have here primarily relied upon the local density gradient 
$\partial \rho / \partial z$ to inform our analysis of spatial mixing patterns. However, there are likely more optimal flow variables, or linear combinations thereof, that could lead to a more robust segmentation of the turbulent domain into distinct regimes. For example, one could imagine constructing more insightful indicator functions of mixing from components of the velocity gradient tensor 
$\partial u_i / \partial x_j$, as suggested by de Bruyn Kops et al. (Reference de Bruyn Kops, Saunders, Rietman and Portwood2019). Applying data-driven techniques, such as unsupervised clustering or dimensionality reduction, to the wealth of observational, experimental and numerical stratified turbulence data currently available has the potential to discover automatically optimal mixing indicators free of human bias. Such an analysis would hopefully further our understanding of the dominant mixing mechanisms present in different flow regimes, along with their prevalence, guiding the search for a more universal and accurate mixing parametrization.
Supplementary movie
A supplementary movie is available at https://doi.org/10.1017/jfm.2023.253.
Acknowledgements
For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission. We are grateful for the detailed feedback received from anonymous referees during the development of this article.
Funding
This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. S.deB.K. was supported under U.S. Office of Naval Research Grant number N00014-19-1-2152.
Declaration of interests
The authors report no conflict of interest.
Appendix. Thresholding of stable filaments
The stable filaments (black) plotted in figure 4(c) were extracted by identifying points within the bottom (most stable) 
$q=15\,\%$, by volume, of the Gaussian-filtered distribution of the density gradient 
$\partial \rho / \partial z$ (figure 4b). We here briefly consider how changing this thresholding percentage 
$q$ influences the characteristics of the extracted stable structures.
In figure 6, we plot the stable structures that are identified by varying the percentage 
$q$ from 
$5\,\%$ to 
$30\,\%$. As 
$q$ is increased, meaning that more points in the stable tail of the filtered 
$\partial \rho / \partial z$ distribution are considered, the identified stable structures are found to grow primarily in the horizontal direction, tracing out more of the stable interfaces identified by the green contours from figure 1(b,d). Figure 6 thus highlights that the magnitude of the vertical density gradient along such stable contours is not uniform, with certain segments having stronger gradients (as identified by using a smaller 
$q$) and thus being characterized by a larger local 
$\chi$.
Figure 6. Extracted stable filaments as a function of the thresholding percentage 
$q$. Filaments (black) are identified as the union of points within the bottom (most stable) 
$q\,\%$ of the filtered density distribution (see discussion in § 4). The filaments detected using values of 
$q$ between 
$5\,\%$ and 
$30\,\%$ are shown in panels (a–f), respectively. Panel (c) corresponds to figure 4(c).
It is also natural to consider how the mixing statistics presented in figure 5(a) depend on the thresholding percentage 
$q$. In figure 7, considering only the computational domain outside of the vortex, we plot the percent contribution of the extracted interfaces to 
$\epsilon$ and 
$\chi$, as a function of 
$q$. The points at 
$q=15\,\%$ correspond to the statistics presented in figure 5(a), noting that in figure 7, the percent contributions are normalized by the domain total outside the vortex, and not the entire domain including the vortex as in figure 5(a). Figure 7 demonstrates that over a wide range of threshold percentages 
$q$, the identified stable filaments always contribute over twice the amount of 
$\chi$ as compared to 
$\epsilon$.
Figure 7. Percent contributions of 
$\epsilon$ and 
$\chi$ contained within the stable interfaces identified in figure 6, normalized by the domain totals outside of the vortex, as a function of the thresholding percentage 
$q$. The statistics at 
$q=15\,\%$ correspond to those presented in figure 5(a), noting that in figure 5(a) the contributions are normalized by the total domain including the vortex.
Open access
