Prandtl number effects on extreme mixing events in forced stratified turbulence

Abstract Relatively strongly stratified turbulent flows tend to self-organise into a ‘layered anisotropic stratified turbulence’ (LAST) regime, characterised by relatively deep and well-mixed density ‘layers’ separated by relatively thin ‘interfaces’ of enhanced density gradient. Understanding the associated mixing dynamics is a central problem in geophysical fluid dynamics. It is challenging to study LAST mixing, as it is associated with Reynolds numbers $Re := UL/\nu \gg 1$ and Froude numbers $Fr :=(2{\rm \pi} U)/(L N) \ll 1$ ($U$ and $L$ being characteristic velocity and length scales, $\nu$ the kinematic viscosity and $N$ the buoyancy frequency). Since a sufficiently large dynamic range (largely) unaffected by stratification and viscosity is required, it is also necessary for the buoyancy Reynolds number $Re_{b} := \epsilon /(\nu N^{2}) \gg 1$, where $\epsilon$ is the (appropriately volume-averaged) turbulent kinetic energy dissipation rate. This requirement is exacerbated for oceanically relevant flows, as the Prandtl number $Pr := \nu /\kappa = {O}(10)$ in thermally stratified water (where $\kappa$ is the thermal diffusivity), thus leading (potentially) to even finer density field structures. We report here on four forced fully resolved direct numerical simulations of stratified turbulence at various Froude ($Fr=0.5, 2$) and Prandtl ($Pr=1, 7$) numbers forced so that $Re_{b}=50$, with resolutions up to $30\,240 \times 30\,240 \times 3780$. We find that, as $Pr$ increases, emergent ‘interfaces’ become finer and their contribution to bulk mixing characteristics decreases at the expense of the small-scale density structures populating the well-mixed ‘layers’. However, extreme mixing events (as quantified by significantly elevated local destruction rates of buoyancy variance $\chi _0$) are always preferentially found in the (statically stable) interfaces, irrespective of the value of $Pr$.


Introduction
In density-stratified flows, small-scale turbulence is known to enhance the rate at which density gradients are irreversibly smoothed by diffusive processes (i.e.mixed) and quantifying † Email address for correspondence: np546@cam.ac.uk Abstract must not spill onto p.2 arXiv:2310.15365v1[physics.flu-dyn]23 Oct 2023 2 this rate in geophysical, environmental and industrial settings is of crucial importance.For instance, turbulent mixing is known to have a leading order impact on global oceanic circulations (Wunsch & Ferrari 2004).Stratified turbulent flows are characterised by a variety of length scales associated with the structure of the turbulent velocity and density fields.Velocity fluctuations are smoothed by molecular viscosity and dissipated into heat below the Kolmogorov scale   := ( 3 /) 1/4 (where  is the molecular viscosity and  is the appropriately volume-averaged dissipation rate of turbulent kinetic energy).Similarly, density fluctuations are dissipated by molecular diffusion below the Batchelor scale   := ( 2 /) 1/4 ,  being the molecular diffusivity.The (molecular) Prandtl number  := /, that effectively quantifies the relative strength of viscous dissipation to diffusion, is the ratio of the square of these two length scales.Unlike the atmosphere, in which  ≃ 0.7 (  ≈   ), in the ocean  ∼ O (10) when thermally stratified and the equivalent Schmidt number  := / ∼ O (1000) ( being the salt diffusivity) when dominantly salt-stratified.
For turbulent flows with  ≫ 1, the inevitable scale separation with   ≪   is such that fully resolved numerical simulations are highly challenging.However, there is increasing evidence indicating that the value of  has (perhaps unsurprisingly) a leading order impact on irreversible scalar-mixing properties of stratified turbulent flows.For instance, the effects of  variations on the properties of secondary instabilities arising from the breakdown of Kelvin-Helmholtz billows were reported by Mashayek & Peltier (2011) and Salehipour et al. (2015).Similarly, Zhou et al. (2017) studied the influence of  variations on fully developed turbulence in stratified plane Couette flows using direct numerical simulations (DNS) and reported significant effects on density and momentum fluxes.Also using DNS, Legaspi & Waite (2020) analysed the effects of  variations on homogeneous forced stratified turbulence and showed, at least for the range of parameters they considered, that the  = 1 simulations are able to describe  > 1 dynamics at large scales and also that kinetic energy spectra (which importantly do not directly contain information from the density field), are largely unaffected by variations in .However, they did observe  effects at small scales below   and in the density flux spectra.Therefore, here we focus on the influence of  variations on small-scale structures of the density field and their associated effects on the mixing properties of forced stratified turbulent flows.
Another important length scale associated with stratified turbulent flows is the buoyancy length scale   := 2/ 0 where  is a characteristic velocity of the flow and  0 is a characteristic (background) value of the buoyancy frequency.Scaled using a characteristic horizontal length scale  ℎ , this parameter then defines a horizontal Froude number  ℎ :=   / ℎ = 2/ 0  ℎ that quantifies the relative strength of the stratification of the flow.For small horizontal Froude numbers  ℎ (i.e.relatively strongly stratified flows), the density field is known to self-organise into relatively well-mixed 'layers' (whose size scales as   ; the buoyancy scale can therefore be thought of as the largest energetically possible overturn) separated by relatively thin 'interfaces' of enhanced density gradient (Billant & Chomaz 2001;Waite 2011).Such density 'staircase' structure has important implications for irreversible scalar mixing in density-stratified turbulent flows.Couchman et al. (2023) showed that (in the  = 1 case) while static instabilities are most prevalent within the wellmixed layers (that are hence characterised by relatively high values of ), much of the scalar mixing, as described by the (appropriately volume-averaged) buoyancy variance destruction rate , is located in the relatively strongly stratified interfaces, a phenomenon that is not apparent when considering  only.Hence, we focus here on the contribution of the different structures described above to the overall mixing properties of the flow, as described by .
Fundamentally, we aim to understand how the density field's small-scale organisation shapes the bulk mixing properties of forced stratified turbulent flows at different Prandtl numbers.To this end, we analyse fully resolved DNS data of forced stratified turbulent flows at  = 1 and  = 7, each with  = 0.5 and  = 2, for values of  such that the emergent   = 50.The rest of this paper is organised as follows.We describe the DNS data in section 2 and then present the methodology used to extract distinct classes of density field structures in section 3.In section 4, we apply this methodology to the DNS data to segment the density field into weakly and strongly stratified interfaces, relatively small-scale 'lamella-like' structures and larger scale density inversions and analyse the associated mixing properties of the extracted density structures.Brief conclusions are drawn in section 5.

Summary of the DNS datasets
We consider statistically steady, forced, fully-resolved DNS of stratified turbulence from the simulation campaign originally reported by Almalkie & de Bruyn Kops (2012).The nonhydrostatic dimensionless Navier-Stokes equations under the Boussinesq approximation, (2.1) are numerically integrated using a pseudospectral code in a triply periodic domain, where ẑ is the (upward) vertical unit vector.The dimensionless parameters are: the above-defined Prandtl number ; the Froude number  := 2/( 0 ) (where  and  are characteristic velocity and length scales associated with the forcing and  0 is the background buoyancy frequency); and the Reynolds number  :=  /.The forcing term F corresponds to the deterministic 'Rf' scheme described in Rao & de Bruyn Kops (2011) and is designed to match a target low wavenumber kinetic energy spectrum at steady-state.This ensures a buoyancy Reynolds number   := /( 2 0 ) ≃ 50 for these simulations, where  is the volume-average of the point-wise local dissipation rate of turbulent kinetic energy  0 , defined in terms of the symmetric part of the strain-rate tensor    as: The density field is a superposition of a background linear density field (), characterised by a reference density  0 and reference density gradient − 2 0  0 /, and a perturbation  ′ (, ) so that, in dimensional form, We can use  ′ to compute the (dimensional) local destruction rate of buoyancy variance: (2.4)This quantity, appropriately scaled in this way by − 0  2 0 / (i.e. the density gradient against which the turbulence acts), corresponds to the local destruction rate of available potential energy (density) and can therefore be used as a proxy of local irreversible mixing (Howland et al. 2021) which can still be calculated pointwise-locally in the flow domain.Henceforth, we consider dimensionless quantities (the prefactor being 1/( ) in our system) and denote by  the volume-average of local  0 across the entire computational domain.
Here, we consider two values of the Prandtl number ( = {1, 7}) as well as two values of the Froude number ( = {0.5, 2}).The simulation parameters are summarised in table 1.We choose the grid spacing Δ ≲   , and a vertical domain dimension of approximately 1.5  .We consider a single snapshot in time (at statistical steady-state) of the various simulations.
For the relatively weakly stratified  = 2 simulations, a patch of elevated turbulence, described by relatively high values of the local dissipation rate  0 , is generated by a relatively large-scale vertically aligned vortex such as the one reported in (Couchman et al. 2023) for the P1F200 dataset.Outside the vortex, the density field self-organises into relatively well-mixed layers separated by interfaces characterised by enhanced density gradients (see figures 1 and 2).In the relatively strongly stratified  = 0.5 case, the vortex does not appear.Couchman et al. (2023) showed (at least for the P1F200 simulation) that the interfaces account for relatively large values of  0 and low values of  0 and hence 'extreme' values of the (local) flux coefficient Γ 0 :=  0 / 0 , whereas the well-mixed layers were potentially more turbulent (relatively large values of  0 ) but characterised by relatively low values of  0 (owing to relatively low values of the density gradient) and hence relatively low values of Γ 0 .Here we aim to understand whether such a picture holds for flows with different values of  and more importantly of , and understand how different structures in the density field (defined more precisely in the next section) influence the mixing within the flow, quantified by values of (local)  0 .

Focus on Fluids articles must not exceed this page length 3. Density field segmentation methodology
Following Couchman et al. (2023), we are interested in the contribution to mixing of the emergent stably stratified interfaces appearing in the density field and how this contribution varies as a function of  and .As  increases, we expect finer structures to arise in the intervening well-mixed regions and we also are interested in understanding how these structures shape the mixing properties of the flow.Therefore, prior to any analysis, the flow is segmented into four different categories based on the local (vertical) density gradient    as well as the neighbouring structure of the flow.We apply the following methodology, as summarized in figure 1: • Stably stratified interfaces: We sort vertical density profiles by values of the density  in order to create statically stable density profiles  * ().We identify points of the dataset whose values of the density remain unaltered by the sorting procedure, i.e. points that lie in statically stable regions of the dataset.These points have relatively high values of the background buoyancy gradient  2 * := −/( 0  2 0 )   * (see figure 1(b-c)) and correspond to stably stratified interfaces lying between relatively well-mixed regions, such as those reported by Couchman et al. (2023) for the P1F200 dataset.We thus refer to such points as belonging to the INT ('interface') cluster, which is further subdivided into relatively strongly stratified (SINT) and relatively weakly stratified (WINT) interfaces, for  2 * > 1 and  2 * ⩽ 1 respectively.
• Relatively well-mixed layers: We subdivide points that are altered by the sorting procedure (i.e. are not within an interface INT) into two categories depending on the local density gradient    and neighbouring structure of the (unsorted) density field , in particular relative to   , the Taylor microscale of the flow.  describes the scale at which viscosity starts to affect the development of turbulent eddies significantly, defined (dimensionally) as   := √︁ 10K/, where K is the volume averaged turbulent kinetic energy.As illustrated in figure 1(e), by considering a point in the dataset at vertical coordinate   and density   satisfying   (  ) > 0, we define  + as the closest point moving upwards for which ( + ) =   .As equality is difficult to ensure numerically, we define depending on the local value of the density gradient), we have identified a relatively small scale structure in the density field, such as a blob or stretched 'lamella' as discussed in detail by Villermaux (2019), which is strongly affected by viscosity.Therefore we classify these structures as belonging to the LAM ('lamella') cluster.Conversely, if we have identified a point belonging to a relatively large scale density inversion largely unaffected by viscosity, and so we classify it as being in the INV ('inversion') cluster.This segmentation methodology, summarised in figure 1, is applied to the four datasets considered in this work, as shown in figure 2.

Cluster properties
4.1.Contributions to  We consider the relative importance of each cluster defined in section 3 in terms of their contribution to the local and bulk mixing properties of the flow, as described by local  0 and volume-averaged .We particularly focus on the role of 'extreme' mixing events and therefore sort the data by decreasing values of  0 , defining a sorted vector where  is the number of points in the dataset and  0 * 0 and   * 0 are the greatest and lowest values of  0 found in the dataset, respectively.This vector can be used to construct the normalized The cumulative contribution   is plotted in figure 3 for each dataset (solid black line), highlighting the fact that the dominant contribution to  is generated by a relatively small set of localised 'extreme' events, regardless of  and .More specifically, for the four sets of parameters considered in this work, approximately 80% of the contribution to  is contained within the first 10% of the points sorted by decreasing values of  0 (and hence 10% of the box volume).Similar statistics are observed in oceanographic data (Couchman et al. 2021).
The total contribution of each cluster to bulk  is shown in figure 3(e).We also assign the data points sorted by values of  0 into  = 20 equal volume bins, in order to calculate a relative contribution of points in each cluster to the bin-volume-averaged ⟨⟩ bin := bin  0 /(/), where  is the total number of points in the domain, as illustrated by colored shading in figures 3(a-d).Focusing on the  = 1 case, the contribution (figure 3(e)) from the strongly stratified interfaces (SINT cluster) to bulk  is roughly 25 − 35%, whereas the contribution from small-scale structures (LAM) reaches ∼ 60%.As  decreases and stratification strengthens, the total contribution from large-scale inversions (INV) shrinks from ∼ 10% for  = 2 to less than 3% for  = 0.5, possibly due to the suppression of large-scale overturnings by relatively strong stratification.We observe similar trends for the relative contributions in each bin.In increasing the Prandtl number from  = 1 to 7, strong interfaces become (perhaps unsurprisingly) finer (figure 2(f,h)) and their total contribution shrinks to Data points are assigned to 20 equal volume bins, sorted by decreasing  0 and clustered using the method presented in section 3.For each bin, we compute the relative contribution of each cluster to ⟨⟩ bin , as shown by the heights of the colored regions.For each bin, we compute the fraction of points within each cluster to the bin's total number of points (red circled dots), along with the (arithmetic) mean value of Γ 0 :=  0 / 0 for each cluster (dashed lines).(e) Total contributions to  from each cluster for the four simulations.about 10%, whereas the total contribution from small-scale structures in 'lamella' (LAM) increases to about 80 − 90%.The total contribution from large-scale inversions does not exceed 10%.Again, we observe similar trends in the relative contributions as  0 decreases.

Statistics of Γ
For each bin in figure 3, we compute a mean value of the (local) flux coefficient Γ 0 :=  0 / 0 associated with each cluster (see dashed lines), although such mean values of ratios of local quantities should always be treated with caution.The leftmost bins (corresponding to the most 'extreme' values of  0 ) have the largest values of Γ 0 (of order unity, significantly above the canonical value Γ = 0.2) suggesting strongly that the extreme mixing events (large  0 ) are not necessarily correlated with extreme turbulent events (large  0 ).Among the extreme events in  0 , those in the strongly stratified interfaces (SINT) cluster correspond to the largest mean values of Γ 0 (as suggested by Couchman et al. (2023) for the  = 1,  = 2 case), independently of variations in  and .This point is further emphasized in figure 4 where we plot the probability density functions (PDFs) of log 10 (  0 ), log 10 ( 0 ) and log 10 (Γ 0 ) for the different clusters and simulations.We also present the statistical mean of  0 ,  0 and Γ 0 .For each simulation, the statistical mean of the flux coefficient Γ 0 is maximized for the strongly stratified interfaces (almost twice as large as the statistical mean for the entire dataset), because of relatively low values of  0 and large values of  0 .As the Prandtl number Figure 4: Probability density functions (PDF) for log 10 (  0 ) (panel (a), (e)), log 10 ( 0 ) (panel (b), (f )) and log 10 (Γ 0 ) (panel (c), (g)) for the different clusters defined in section 3 and for the four simulations.The statistical mean of each field (without the logarithm) is represented by a colour-coded circle ( = 1) or triangle ( = 7).The green dotted vertical line corresponds to the canonical value Γ = 0.2.increases, the statistical mean of Γ 0 decreases (a result also observed by Salehipour et al. (2015)) for both the weakly and strongly stratified cases, for the following (different) reasons.
In the weakly stratified simulations with  = 2 (figures 4(a-c)), density effectively acts as a passive scalar (as can be seen from the weakly stratified scaling of Riley et al. (1981), for instance) and hence the statistics of  0 do not depend on , as can be seen in figure 4(b).However, as  increases, the left tails of the  0 distributions become more significant, as might be expected for the mixing of an effectively passive scalar, as  0 is multiplied by ( ) −1 where  is fixed.We note that a widening of the tails of    ′ (not shown here, see for instance Riley et al. (2023)) effectively rebuilds the right tail of  0 , explaining why the PDF of  0 is not just shifted toward lower values.A more in-depth discussion of this phenomenon is provided by Bragg & de Bruyn Kops (2023).As a result of the statistics of  0 remaining roughly independent of  but those of  0 decreasing with increasing , the statistics of Γ 0 decrease as  increases for the weakly stratified ( = 2) case.
Conversely, in the strongly stratified simulations with  = 0.5, buoyancy acts 'actively' on the momentum field, as demonstrated for instance by the strongly stratified scaling analysis of Billant & Chomaz (2001).Moreover, as  increases, the volume contribution of the smallscale structures (LAM cluster) increases at the expense of the strongly stratified interfaces (SINT) (see figure 3).LAM structures, which locally are only weakly affected by stratification because they populate the relatively well-mixed regions of the flow, are viscously affected (as their vertical extent is, by definition, smaller than   ) but still significantly disordered.As a result, their local static instability inevitably encourages increased viscous dissipation, and so their enhanced prevalence actually shifts the statistics of  0 towards higher values for the  = 7 flows as compared to  = 1 (yellow curves, figure 4(f )).In this sense, enhanced stratification at higher  actually makes the flow 'more turbulent', reminiscent of the prediction by Pearson & Linden (1983) of relatively long-lived 'approximately horizontal striations' in high- decaying stratified turbulence.Also, though this is a second order effect, for the strongly stratified simulations the left tail of the  0 distribution once again becomes somewhat more significant as  increases.We note that while the forcing scheme endeavors to maintain a constant   , bulk  is slightly mismatched between the two  simulations at  = 0.5 (see Table 1).Therefore, while increased  dramatically changes the relative prevalence of LAM versus SINT structures, part of the observed increase in local  0 statistics at  = 7,  = 0.5 is due to the slight mismatch in targeted bulk  rather than being a purely  effect.Importantly, it is apparent in figure 4(e) that the total amount of irreversible mixing for the flow with  = 7 actually slightly increases compared to  = 1 in the more strongly stratified  = 0.5 simulations, essentially because the  = 7 flow is more vigorous, despite the fact that Γ 0 appears to decrease with increasing .Therefore, consideration of Γ 0 , or even Γ := / (constructed from volume-averaged dissipation rates), in isolation should be treated with caution.

Discussion
We have analysed the influence of the Prandtl number  and Froude number  on density structures in forced turbulent stratified flows and their contribution to irreversible scalar mixing properties.Using fully resolved DNS data and a flow segmentation algorithm based on the local value of the (vertical) density gradient and the local (vertical) structure of the density field, we have extracted distinct regions of the turbulent density field -interfaces (both relatively strong and relatively weak) separating well-mixed density layers made up of both small-scale lamellar structures and larger scale density inversions -and analysed their contribution to the bulk value of the destruction rate of buoyancy variance  and to the statistics of the (locally evaluated) flux coefficient Γ 0 .
As  increases, the strongly stratified density interfaces become finer and their contribution to average values of  decreases at the expense of the small-scale lamellar structures in the relatively well-mixed density layers.However, similarly to the flow with  = 1 (Couchman et al. 2023), these structures are 'quiescent', yet mixing hotspots.The points in these structures associated with the largest values of (local)  0 are associated with extreme values of Γ 0 and hence with relatively low values of  0 .More generally, the flux coefficient associated with strongly stratified structures is (in an averaged sense) essentially twice as large as its value for the three other segmented regions, regardless of the values of  and .All in all, strongly stratified interfaces are therefore characterised by relatively large values of  0 and relatively weak values of  0 compared to the other density structures considered in this work and are therefore likely to be overlooked if considering  0 (or indeed the volume average ) alone as a proxy for significant irreversible mixing.
As is becoming increasingly well-appreciated, a universal value of Γ is not able to capture the complex and inhomogeneous structures of the density and density gradient turbulent fields.Moreover, our description of the density fields as well as consideration of the considered mixing regime (i.e.'passive' mixing at higher  as compared to 'active' mixing at lower ) offers a potential explanation for the empirical observation that the 'mixing efficiency' Γ/(1 + Γ) decreases with .On the one hand, mixing in the relatively weakly stratified case can be described as 'passive': the velocity field and hence  0 are largely unaffected by  but the small-scale structure of the density field evolves as  increases from 1 to 7 in a way that results in a decrease in  0 and thus the flux coefficient Γ 0 .On the other hand, mixing in the relatively strongly stratified case can be thought of as 'active', with the density field and its small-scale lamellar structure having a leading order impact on the rate at which the turbulent kinetic energy is dissipated ( 0 ), ultimately leading to a decrease in Γ 0 as  increases.This shows another way in which a focus on the flux coefficient can be misleading.Since the decrease in Γ 0 is actually principally related to an increase in  0 , the total amount of mixing (quantified by ) actually increases in the higher , strongly stratified flow considered here.Another key result of our analysis is that the Prandtl number  has a leading order impact (at least compared to the Froude number ) on the density field structures and their mixing properties, further emphasizing the need to take this parameter into account when 'measuring mixing', a process that is inherently a diffusive one (Villermaux 2019;Caulfield 2021) and that affects the density field.
We here considered steady-state forced stratified turbulent flows, and so it would now be interesting to apply our segmentation and associated mixing contribution analysis to timedependent decaying stratified turbulent flows.For example, what is the fate of the different emerging structures?Can a separate analysis of the temporal evolution of each structure help us better understand the mixing history of a stratified turbulent flow as a whole and unveil the fundamental rules of stratified mixing?These are questions left for future work.

Figure 1 :
Figure1: Flow segmentation methodology.For the P1F200 simulation, the density field, represented here by its vertical gradient as shown in panel (), is vertically sorted into a statically stable field, whose vertical gradient is denoted  2 * (up to a multiplicative constant), as shown in panel ().The sorting algorithm highlights the stable interfaces of the density field, i.e. the points of the density field unaltered by the sorting procedure.The strongly stratified interfaces (SINT, in blue with  2 * > 1) are then extracted in panel ().The entire segmented field is shown in panel (), including the weakly stratified interfaces with  2 * ⩽ 1 (WINT, in purple), and the relatively well-mixed regions between interfaces, further sub-divided into small-scale 'lamella' structures (LAM, in orange) and larger-scale density inversions (INV, in red) using the procedure described in section 3 and shown schematically in panel (), based on the Taylor microscale   .

Figure 3 :
Figure 3: (-) Normalized cumulative contribution to  (black line; see equation (4.1)) for each simulation.Data points are assigned to 20 equal volume bins, sorted by decreasing  0 and clustered using the method presented in section 3.For each bin, we compute the relative contribution of each cluster to ⟨⟩ bin , as shown by the heights of the colored regions.For each bin, we compute the fraction of points within each cluster to the bin's total number of points (red circled dots), along with the (arithmetic) mean value of Γ 0 :=  0 / 0 for each cluster (dashed lines).(e) Total contributions to  from each cluster for the four simulations.

Funding.
This project received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 956457 and used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, 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.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.

Table 1 :
Summary of the DNS data.  ,   and   denote the Kolmogorov, Taylor and buoyancy scales.