Irreversible mixing by unstable periodic orbits in buoyancy dominated stratified turbulence
TThis draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Irreversible mixing by unstable periodicorbits in buoyancy dominated stratifiedturbulence
Dan Lucas † & C. P. Caulfield School of Computing and Mathematics, Keele University, Staffordshire, ST5 5BG Department of Applied Mathematics & Theoretical Physics, University of Cambridge, Centrefor Mathematical Sciences, Wilberforce Road, Cambridge, CB3 0WA, UK BP Institute, University of Cambridge, Madingley Rise, Madingley Road, Cambridge, CB30EZ, UK(Received xx; revised xx; accepted xx)
We consider turbulence driven by a large-scale horizontal shear in Kolmogorov flow(i.e. with sinusoidal body forcing) and a background linear stable stratification withbuoyancy frequency N B imposed in the third, vertical direction in a fluid with kinematicviscosity ν . This flow is known to be organised into layers by nonlinear unstable steadystates, which incline the background shear in the vertical and can be demonstratedto be the finite-amplitude saturation of a sequence of instabilities, originally from thelaminar state. Here, we investigate the next order of motions in this system, i.e. the time-dependent mechanisms by which the density field is irreversibly mixed. This investigationis achieved using ‘recurrent flow analysis’. We identify (unstable) periodic orbits, whichare embedded in the turbulent attractor, and use these orbits as proxies for the chaoticflow. We find that the time average of an appropriate measure of the ‘mixing efficiency’ ofthe flow E = χ/ ( χ + D ) ( D is the volume-averaged kinetic energy dissipation rate and χ is the volume-averaged density variance dissipation rate) varies non-monotonically withthe time-averaged buoyancy Reynolds numbers Re B = D / ( νN B ), and is bounded aboveby 1 /
6, consistently with the classical model of Osborn (1980). There are qualitativelydifferent physical properties between the unstable orbits that have lower irreversiblemixing efficiency at low Re B ∼ O (1) and those with nearly optimal E (cid:46) / Re B ∼
10. The weaker orbits, inevitably embedded in more stronglystratified flow, are characterised by straining or ‘scouring’ motions, while the moreefficient orbits have clear overturning dynamics in more weakly stratified, and apparentlyshear-unstable flow.
1. Introduction
A fundamental outstanding problem in stratified turbulence relates to the mixingof the stratifying agent, e.g. heat or salinity. Moreover such flows typically exhibitstrong anisotropy with vertical motions being suppressed due to gravitational forces,and layerwise motions commonly being observed. The route by which the flow fieldrearranges the density field into well-defined mixed regions or ‘layers’ separated by sharpgradients or ‘interfaces’ as well as the layer/interface structure’s robustness are subjectsof continuing research. A widely utilised characterisation of mixing comes in the formof some appropriate measure of the ‘mixing efficiency’, often defined as the ratio ofirreversible potential energy increase relative to the irreversible kinetic energy loss (Peltier † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] S e p & Caulfield 2003). One reason for this interest in mixing efficiency stems from the effortto devise accurate parameterisations of diapycnal diffusivity for use in ocean models, andthere is growing evidence that such parameterisations should take into account variationof the efficiency with control parameters (Salehipour et al. et al. unstable solutionsembedded within the stratified turbulence. Such an approach has seen great success inclarifying the behaviour of transitional wall-bounded shear flows (Kawahara & Kida 2001;Cvitanovi´c & Gibson 2010; Kawahara et al. et al. et al. (2017), (henceforth LCK17) have shownhow such an approach can advance our understanding of layer formation by locatingnonlinear layered steady states about which the turbulence organises. Given this success,we are motivated to investigate whether the next order of motions can be identified,those time-dependent simple invariant manifolds, i.e. periodic orbits, which capture somesalient signature of the processes by which buoyancy is mixed. In particular, we areinterested in how ‘efficient’ (defined in a precise fashion below) such mixing is, and howsuch processes vary with control parameters.To address these issues, the paper is organised as follows. Section 2 contains theformulation and discussion of the methods employed, while section 3 presents a set ofpreliminary and motivational direct numerical simulations (DNS) with a discussion oftheir mixing properties. Section 4 shows results from the recurrent flow analysis togetherwith a discussion of the processes exhibited by the periodic orbits discovered, and finallysection 5 presents our conclusions.
2. Formulation
We begin by considering the following version of the non-dimensionalised, monochro-matic body-forced, incompressible, Boussinesq equations ∂ u ∂t + u · ∇ u + ∇ p = 1 Re ∆ u + sin( ny ) ˆx − Bρ ˆz , (2.1) ∂ρ∂t + u · ∇ ρ = w + 1 ReP r ∆ρ, ∇ · u = 0 (2.2)where we define the (external) Reynolds number Re , the bulk stratification parameter B and the Prandtl number P r as Re := √ λν (cid:18) L y π (cid:19) / , B := gβL y ρ λ π , P r = νκ . (2.3)Here, u ( x, y, z, t ) = u ˆx + v ˆy + w ˆz is the three-dimensional velocity field, p is the pressureand the density is decomposed into ρ tot = ρ + ρ B ( z ) + ρ ( x , t ), i.e as the sum of aBoussinesq reference density, a constant linear background stratification and a fullyvarying disturbance density. We have non-dimensionalised using the characteristic lengthscale L y / π , characteristic time scale (cid:112) L y / πλ and density gradient scale β = d ρ ∗ B / d z. (see LCK17 for details). Furthermore, n is the forcing wavenumber, λ is the forcingamplitude, ν is the kinematic viscosity, κ is the molecular diffusivity. We impose periodicboundary conditions in all directions and solve over the cuboid [0 , π/α ] × [0 , π ] where α = L y /L x defines the horizontal aspect ratio of the domain. Vorticity ω = ∇ × u isused as the prognostic variable and DNS are performed using the fully dealiased (two-thirds rule) pseudospectral method with mixed fourth order Runge-Kutta and Crank-Nicolson timestepping implemented in CUDA to run on GPU cards. We initialise thevelocity field’s Fourier components with uniform amplitudes and randomised phasesin the range 2 . (cid:54) | k | (cid:54) . (cid:104)| ω | (cid:105) V = 1 and ρ (cid:48) =0 initially. Throughout, (cid:104) ( · ) (cid:105) V := α (cid:82)(cid:82)(cid:82) ( · ) dxdydz/ (2 π ) denotes a volume average, (cid:104) ( · ) (cid:105) h := α (cid:82)(cid:82) ( · ) dxdy/ (2 π ) denotes a horizontal average and (cid:104) ( · ) (cid:105) v := (cid:82) ( · ) dz/ (2 π )denotes a vertical average.We define the diagnostics involved in the energetic budgets as K = 12 (cid:104)| u | (cid:105) V , P = B (cid:104) ρ (cid:105) V , I = (cid:104) u · f (cid:105) V = (cid:104) u sin( ny ) (cid:105) V , (2.4) B = (cid:104) u · Bρ ˆz (cid:105) V = B (cid:104) wρ (cid:105) V , D = 1 Re (cid:104)|∇ u | (cid:105) V , χ = BP rRe (cid:104)|∇ ρ | (cid:105) V , (2.5)where d K d t = I − B − D , d P d t = B − χ and K is the total kinetic energy density, P the density variance, I is the energy input by the forcing, B is the buoyancy flux, D isthe viscous dissipation rate, and χ the density variance dissipation rate. We fix α =0 . P r = 1 for numerical efficiency and n = 1, i.e. theflow is forced with sin( y ) ˆx , to mimic closely other large scale shear profiles previouslystudied (e.g. stratified Taylor-Couette flow (Woods et al. et al. · ) = [ (cid:82) T ( · )d t ] /T where T is normally the full simulation time.
3. Irreversible mixing in the direct numerical simulations
We begin by characterising the mixing in this system by conducting a set of simulationsacross a range of Re and B. The pertinent single point diagnostics are presented in table 1.In particular we examine the irreversible mixing efficiency which we, following Salehipour& Peltier (2015) and Maffioli et al. (2016), define in terms of the dissipation rates: E ( t ) = χχ + D . (3.1)The dependence of E on typical flow parameters is a crucial ongoing problem instratified turbulence, and a source of some controversy. As discussed by Mater & Ve-nayagamoorthy (2014), the flow of a turbulent, shear-driven stratified fluid has severalcompeting time scales, and it is natural to attempt to parameterise E in terms of param-eters quantifying the relative importance of these time scales. Using shear-instabilitysimulations and comparison with observations, Salehipour et al. (2016) developed aparameterisation using two parameters: an appropriately defined ‘buoyancy Reynoldsnumber’, effectively a ratio of the time scale of the stratification to the time scale ofthe turbulence; and an appropriately defined Richardson number, effectively a ratio ofthe square of the time scale of the vertical shear (assumed to be the dominant driver ofthe mixing) to the time scale of the stratification. In order to determine how the mixingefficiency varies with these dimensionless numbers in this flow, we define the buoyancyReynolds number Re B ( t ) and the (local) gradient Richardson number Ri G ( x , t ) as (withthe usual caveat that comparing different studies is difficult if the key parameters aredefined differently) Re B ( t ) = D ReB , Ri G ( x , t ) = − B ∂ρ tot ∂z (cid:0) ∂u h ∂z (cid:1) , (3.2)where ( ∂u h /∂z ) = ( ∂u/∂z ) +( ∂v/∂z ) . Ri G is a pointwise quantity which characterisesthe relative stability of the flow to overturning shear instabilities, and developing anappropriate overall average for this quantity to characterise a particular flow should beperformed with care, not least because locations of low vertical shear significantly skewthe distributions of Ri G . The appropriate averaging used here is thus Ri B ( t ) = (cid:42) (cid:68) − B ∂ρ tot ∂z (cid:69) h (cid:68)(cid:0) ∂u h ∂z (cid:1) (cid:69) h (cid:43) z . (3.3)There are two important points to appreciate about these definitions. First, we havechosen to define Ri B with the total gradients of density while Re B only dependsexplicitly on the background N B . Second, although in principle these are independentnondimensional parameters, it remains to be established that they are in this flow, asthere exist situations where they are strongly correlated (see e.g. Zhou et al. (2017)).In figure 1(a) and (b) we plot the time-averaged mixing efficiency E against the time-averaged parameters Ri B and Re B . We observe non-monotonic dependence of E on bothparameters, in at least qualitative agreement with the experimental analysis of Linden(1979) and numerical simulations of Shih et al. (2005). Indeed, E appears to saturatenear the critical value of 1 /
6, (marked with a horizontal line) equivalent to the upperbound on the turbulent flux coefficient Γ (cid:39) E / (1 − E ) (cid:54) . Re B (cid:38) , with a dependence not entirely unlike the E ∝ Re − / B observed by Shih et al. (2005) in what they referred to as the ‘energetic’ regime (see Ivey et al. (2008) for a more detailed discussion).However, although these similarities are intriguing, there are still significant differences.Firstly, we do not observe a ‘transitional’ plateau of approximately constant mixingefficiency at intermediate Re B . Indeed, the qualitative structure of the mixing efficiencycurve is more reminiscent of the Pad´e approximant proposed by Mashayek et al. (2017)using (vertical) shear-instability numerical data as well as observations, where Γ ∝ Re / B and Re − / B for small and large Re B respectively. Second, and once again reminiscentof the approach proposed by Mashayek et al. (2017) to parameterise mixing in termsof Re B alone, the two parameters Ri B and Re B are actually closely correlated, and oursimulations actually trace out a (monotonic) curve in Ri B − Re B space, as shown in figure1(c), independently of the externally imposed parameters B and Re . The low mixingefficiency observed at high Ri B can thus be understood as being entirely associated withweaker turbulence (smaller Re B ) and the maximum mixing efficiency occurs at a sweetspot of sufficiently vigorous turbulence at sufficiently high stratification, analogously tothe results of Zhou et al. (2017). We stress that we are not claiming that this clearcorrelation between Ri B and Re B is generic, (see for example the discussion in Scotti &White (2016)) just that it occurs for this flow.It is apparent that at small Ri B , Ri B ∝ Re − B . This is unsurprising, as the turbulenceand shear are largely unaffected by such weak stratification, and the variation of bothparameters with B completely dominates. As Ri B increases beyond the value associatedwith the most efficient mixing however, the turbulence and shear do indeed becomesuppressed by the strengthening stratification, and the power law dependence steepensso Ri B ∝ Re − / B . There is undoubtedly also an increasing significance of viscosity, whichbecomes dominant for the largest values of Ri B >
1, where Re B drops below one, and theturbulence is essentially completely suppressed, with very inefficient mixing. However, itis always important to remember that both Ri B and Re B are global measures averaged . . . . . . . . . . .
01 0 . . . . . . . . . . . . . .
01 0 . E Ri B UPO-o1UPO-l1 / E Re B Re < Re = 50 Re = 100 Re B Ri B Re < Re = 50 Re = 100 Re − B Re − / B Re − / B r: (a)for:(b). (c) Figure 1.
The variation of average mixing efficiency E with: (a) average bulk Richardsonnumber Ri B and (b) average buoyancy Reynolds number Re B . The horizontal line is E = 1 / Ri B , E ) and ( Re B , E ) planes of the two unstable relative periodic orbitsUPO-o1 (solid line) and UPO-l1 (dashed line) discussed in section 4. (c) shows the variation of Re B with Ri B . Different symbols mark different values of Re . The red curve in (b) correspondsto E = Re − / B , the large Re B behaviour suggested by Shih et al. (2005). in both space and time. Both Ri B and Re B are functions of time, and there could alsobe spatial variation of Ri G within the flow domain. Such variation proves to be a keypart of the behaviour of the periodic orbits which we identify, and hence of the mixingoccuring in these flows.
4. Recurrent flow analysis
In order to determine what processes underpin the turbulence in our simulations,we use ‘recurrent flow analysis’. The general approach is as follows. First a simulationis conducted during which near recurrences are located in the chaotic trajectory. Thisis achieved by storing a historical record of state vectors and periodically checking thehistory against the current state vector. When an appropriately ‘near repeat’ is identified,it is stored for a later convergence attempt with a Newton-GMRES-hookstep algorithm.For more details on this approach see Chandler & Kerswell (2013); Lucas & Kerswell(2015). Here, we are interested in finding representative unstable periodic orbits (UPOs)with qualitatively different mixing properties in a broad range of flows, rather thanidentifying a large number of UPOs for a specific set of parameters. Therefore, we haveidentified five orbits in two broad classes: one class of two orbits found in the relativelyweakly stratified simulation A1 with Re = 20 and B = 1, and the other class of three Re B Ri B ηk max E Re B l o Re λ A1 20 1 0.325 2.7 0.0859 9.8 0.695 31.3A2 30 1 0.235 2.0 0.1026 14.9 0.698 51.2A3 30 5 0.525 1.8 0.0892 4.04 0.244 25.5A4 40 5 0.489 1.5 0.114 4.64 0.226 52.1A5 40 7.5 0.58 1.5 0.099 3.45 0.177 31.6A6 40 10 0.835 1.6 0.096 1.94 0.124 23.6B1 50 5 0.398 1.3 0.121 5.89 0.228 44.7B2 50 10 0.608 1.3 0.112 3.21 0.141 44.1B3 50 50 3.18 1.2 0.054 0.83 0.048 31.3B4 50 100 14.95 1.4 0.038 0.25 0.022 45.9C1 100 0.5 0.046 1.6 0.109 94.82 1.15 98.5C2 100 0.75 0.065 1.6 0.129 65.19 0.86 92.5C3 100 1.0 0.089 1.6 0.144 48.03 0.69 94.1C4 100 1.5 0.126 1.6 0.156 31.41 0.50 91.2C5 100 2.0 0.156 1.6 0.160 24.3 0.41 88.6C6 100 5.0 0.278 1.5 0.153 10.8 0.22 76.2C7 100 10.0 0.423 1.5 0.137 5.91 0.14 62.9
Table 1.
Imposed parameters and diagnostic outputs for the three groups of simulations.The nondimensional Kolmogorov microscale is given by η = Re − / D − / and is scaled by themaximum wavenumber allowable k max = N/
3. Groups A and B have resolution 128 × andgroup C has 256 × to retain spectral convergence. The nondimensional Ozmidov lengthscale l O = D / /B , the average buoyancy Reynolds number Re B = D Re/B, the average bulkRichardson number Ri B and the Taylor microscale Reynolds number Re λ = K turb (cid:113) Re/ D are also listed. Here we define K turb = (1 / (cid:104) ( u − u ) (cid:105) V and overbars are always time averagesover the full T = 1000 window. orbits in the more strongly stratified simulation B3 with Re = 50 and B = 50. The firstclass has Re B ∼ O (10), and we label it as class ‘o’, for ‘overturning’, while the secondclass has Re B ∼ O (1), and we label it as class ‘l’ for layered. We plot projections ontothe time-dependent ( Ri B , E ) and ( Re B , E ) planes of the trajectories of characteristicrelative periodic orbits UPO-o1 (solid line) and UPO-l1 (dashed line) on figures 1(a) and(b), demonstrating that these two classes do indeed have qualitatively different mixingproperties, and also that these time-dependent properties are consistent with the time-averaged properties of the ‘full’ DNS.Table 2 lists properties of the converged orbits, including period, relative shifts dueto the continuous symmetries in x and z, stability and some diagnostic averages. Bycomparison with the averages of the DNS in table 1, it is clear the UPOs are repro-ducing the bulk time-averages quite well. The projection of the trajectories onto the( Re B , E ) plane of the five converged orbits are plotted in figure 2, with greyscale coloursrepresenting the probability density function (p.d.f.) of the turbulent DNS. Notice thatthe darker colours represent regions where the turbulent trajectory spends more time,and lighter colours represent less frequent excursions to that part of phase space. Forsimulation B3, the ‘l’ class of converged UPOs sit across the darkest region of the p.d.f.where the turbulence spends most of its time but miss the higher efficiency, intermittentturbulent bursts that the DNS exhibits. Missing extreme events is a known failing ofthe recurrent flow analysis in some circumstances and is a topic of ongoing research, seeLucas & Kerswell (2015). However, for the ‘o’ class in simulation A1, the orbit UPO-o1 No. DNS
Re B T p s x s z m y Λ − N λ Re B max( E ) Ri B E UPO-o1 A1 20 1 21.28 5.98 0.35 0 1.65 25 8.5 0.13 0.333 0.078UPO-o2 A1 20 1 7.51 7.81 -0.06 0 2.30 27 8.9 0.097 0.265 0.067UPO-l1 B3 50 50 8.75 0.016 0 0 0.211 7 0.62 0.07 3.17 0.052UPO-l2 B3 50 50 8.74 0 0 0 0.207 6 0.63 0.058 3.13 0.048UPO-l3 B3 50 50 9.21 0 0.007 0 0.116 6 0.57 0.051 4.64 0.045
Table 2.
Table cataloguing the unstable recurrent flows showing the DNS from which theycome (table 1), period T p , relative shifts in x and z directions- s x , s z discrete shift/reflect in ym y , inverse stability coefficient ( Λ − , see Chandler & Kerswell (2013)), the number of unstabledirections N λ , Re B , maximum and time averaged mixing efficiency E and time averaged bulkRichardson number Ri B . (in particular) appears to span the turbulent attractor well in this projection, missingonly some very rare excursions to high Re B . As shown in table 2, the ‘o’ class orbits areconsiderably more unstable than the orbits in ‘l’ class.In order to compare how the timescales of the UPOs relate to those exhibited bythe turbulence, we plot in figure 3 the power spectra of the kinetic energy for the twoDNS signals from which the UPOs have been extracted (i.e. A1 and B3). Vertical linesshow that the the periods of the orbits are approximately coincident with observedfrequencies. For simulation A1 the largest non-zero peak in the frequency spectrumoccurs at approximately half the fundamental frequency observed in the UPOs, i.e. f ≈ π/
15 = 0 .
42 with the fundamental period around 7 .
5, UPO-o1 having just underthree such periods. For simulation B3 the periods of the orbits are on the edge of themain distribution of frequencies near f = 0 . . This is reasonable evidence that thetimescales of the UPOs are representative of the typical timescales of the turbulence.We can also consider the buoyancy frequency and its influence on these cases. Forsimulation A1 where B = 1, the background buoyancy frequency N B = (cid:112) ( B ) = 1so the timescale or ‘buoyancy period’ is T N = 2 π which is not far from the “fundamentalperiod” we have established from the UPOs and spectrum. For B3, the buoyancy period T N = 2 π/ (cid:112) (50) ≈ . Ri B and E for UPO-o1 and UPO-l1. The dynamicsof UPO-o1 show relatively low Ri B , strongly correlated to the mixing efficiency E . Thereare three distinct sub-periods, each of which shows the bulk Ri B decrease to small values(indeed below the canonical value of 1 / ρ tot andstreamwise velocity u are chosen before and after the final minimum of Ri B and nearthe subsequent maximum of Ri B following the peak E as marked on figure 4(a). (Ananimation of the time evolution of these fields is available as supplementary material.)The density field shows very distinct Kelvin-Helmholtz-like billows forming at t = 16 . Re B E UPO-l1UPO-l2UPO-l3 Re B E UPO-o1UPO-o2
Figure 2.
Projection onto the ( Re B , E ) plane of the UPOs and the p.d.f. of the direct numericalsimulations rendered in greyscale, with darker shades denoting that the turbulent trajectoryspends more time there for: (a) class ‘o’ UPOs from simulation A1; and (b) class ‘l’ UPOs fromsimulation B3. ˆ K ( f ) fT = 21 T = 7 . Re = 20 , B = 1 ˆ K ( f ) fT = 9 Re = 50 , B = 50 Figure 3.
Power spectra for the kinetic energy ˆ K ( f ) = (cid:82) T K ( t )e i ft d t for the DNS signals fromsimulation A1 (left) and simulation B3 (right) with vertical lines showing approximately theperiods of the UPOs extracted. R i B E t (a) Ri B E R i B E t (b) Ri B E Figure 4.
Time evolution of Ri B and E for: (a) UPO-o1; and (b) UPO-l1. Notice the correlationbetween Ri B and E for (a) and not for (b). The times highlighted in figures 5 and 6 are markedwith stars. where E and Ri B are beginning their growth phase, justifying the labelling of this UPOas being in class ‘o’. Of course, a bulk measure of the Richardson number does not capturethe stability properties of the flow, and so in figure 5 we also plot the p.d.f. of Ri G ( x ) overthe spatial domain at the same times, with a vertical line indicating Ri G = 1 /
4. Thesedistributions show a striking increase of the proportion of the domain with Ri G < / bulk Ri B some regionsof the domain remain with Ri G < /
4. Apparently, it is necessary for an appreciableproportion of the domain to have a local gradient Richardson number Ri G < / Ri B ,seemingly uncorrelated to the behaviour of the markedly lower mixing efficiency E . Thisis also apparent in the snapshots shown in figure 6 at the characteristic times markedon figure 4(b) (and the associated animations available as supplementary materials). Nooverturns occur, and the cycle of mixing behaviour is associated with straining or scouringmotions drawing the perturbation density ρ into thin layers, justifying the labelling asclass ‘l’. The instantaneous p.d.f.s of the gradient Ri G plotted in figure 6 now show thatnowhere in the domain has Ri G < / , even near t = 5 . E , t ≈ . , the snapshot of ρ shows increased layers and arguablylarger ∂ρ/∂z. In order to characterise the straining/scouring motions exhibited by the ‘l’ class, weexamine the properites of the strain tensor S ij = 12 (cid:18) ∂u i ∂x j + ∂u j ∂x i (cid:19) . This tensor has three real eigenvalues, which we denote by ˜ α, ˜ β, ˜ γ, referred to as theprincipal strains, which sum to zero in an incompressible flow. In general ˜ α > γ < β will control the three dimensionality, i.e. ˜ β < β > β = 0 a purely two-dimensional straining field.Figure 7 shows time series of the volume-averaged strains for the two characteristic UPOs,UPO-o1 and UPO-l1. In the overturning case, UPO-o1 shows that (cid:104) ˜ β (cid:105) V is well correlatedto the mixing efficiency, being positive and maximum when E is largest. This is in goodagreement with Smyth (1999) who examined these strains for a freely decaying Kelvin-Helmholtz unstable shear layer; (cid:104) ˜ β (cid:105) V is positive and has a distinct maximum during thesecond phase of the overturn where the flow becomes more isotropic (isotropic turbulenceis observed to have the strains approximately in the ratio 3:1:-4 (Smyth 1999; Ashurst et al. (cid:104) ˜ β (cid:105) V are an order of magnitude smaller and are actuallynegative for much of the cycle. At the point where the mixing efficiency is maximal( t ≈ . (cid:104) ˜ β (cid:105) V also has a maximum, but here this amounts to crossing the axes and (cid:104) ˜ α (cid:105) V ≈ −(cid:104) ˜ γ (cid:105) V . Our interpretation is then that the flow approaches two-dimensionalitywhere the mixing is largest and increased vertical gradients of ρ are formed by the changein sign of (cid:104) ˜ β (cid:105) V . This is in stark contrast to the overturning case where mixing has a moreisotropic straining signature with large (cid:104) ˜ β (cid:105) V . As a secondary observation, the higherfrequency oscillations of (cid:104) ˜ β (cid:105) V for UPO-l1 are of the order of N = √ B the buoyancyfrequency, i.e. the buoyancy period here is T N = 2 π/ √ ≈ . . This suggests thatinternal waves are also playing a role in this dynamical cycle.
5. Discussion and Conclusions
In this paper we have successfully applied recurrent flow analysis to stratified flows forthe first time. As expected, success is restricted to relatively weak turbulence at modestReynolds numbers. Nevertheless, the approach has still provided detailed insight into0 ˆ x ˆ y ˆ z t = 15 t = 16 . t = 20 ρ tot u . . . . . .
01 0 . Ri G . . . − ( Ri G ) . . . . . .
01 0 . Ri G . . . − ( Ri G ) . . . . . .
01 0 . Ri G . . . − ( Ri G ) t = 15 t = 16 . t = 20 t = 15 t = 16 . t = 20 Figure 5.
Three dimensional rendering of: (top) ρ tot = ρ − z ; (middle) u ; and (bottom) p.d.f.sof Ri G in the entire computational domain at t = 15 , . Ri G , while the inset shows the result whendistributing the intervals on a log scale. Vertical lines mark Ri G = 1 / the sustaining mechanisms and mixing processes at work in these flows. In particular wehave established that mixing can occur from two distinct mechanisms; overturning andscouring, consistently with the classification presented by Woods et al. (2010). Scouringis observed when stratification is strong and overturns when the stratification is weakenough to allow small gradient Richardson numbers. By examining the cycles in spaceand time we have demonstrated that the spatial distribution of local gradient Ri G ismarkedly skewed before high overturning mixing, but remains bounded from below by1 / Ri G . Crucially, weaker average mixing may not be simply controlled byspatiotemporally intermittent shear instability, but by other scouring, straining processes,which may require other predictive diagnostics. Furthermore, since the bulk measures ofRichardson number and buoyancy Reynolds number are functionally related, the resultsof our motivational DNS for horizontally forced, sustained turbulence suggest that an( Ri B , Re B ) parameterisation of mixing should be treated with caution, particularlysince it proved difficult to access a strongly stratified, yet strongly turbulent regime,analogously to the situation arising in stratified plane Couette flow (Zhou et al. . . . . . Ri G . . . . . Ri G . . . . . Ri G ˆ x ˆ y ˆ z uρ t = 4 t = 4 t = 4 t = 5 . t = 5 . t = 5 . t = 7 t = 7 t = 7 Figure 6.
Three dimensional rendering of: (top) ρ tot = ρ − z ; (middle) u ; and (bottom) p.d.f.sof Ri G in the entire computational domain at t = 4 , . Ri G , while the inset shows the result whendistributing the intervals on a log scale. Vertical lines mark Ri G = 1 / h ˜ β i V h ˜ α i V −h ˜ γ i V t -0.02-0.015-0.01-0.00500.0050.01 0 1 2 3 4 5 6 7 8 9 1.41.451.51.551.61.651.71.751.81.851.91.95 h ˜ β i V h ˜ α i V −h ˜ γ i V t h ˜ β i V h ˜ α i V −h ˜ γ i V h ˜ β i V h ˜ α i V −h ˜ γ i V Figure 7.
Time series of the volume averaged principal strains for UPO-o1 (left) and UPO-l1(right). For the overturning case, peaks of E are coincident with large positive intermediatestrain (cid:104) ˜ β (cid:105) V , nearer to the values of isotropic turbulence, where as the more layered case (right)is very anisotropic with small (cid:104) ˜ β (cid:105) V , approximately zero when E is largest. One open question is the interplay between scouring and overturning. Clearly over-turning shear instability will overwhelm any background straining or scouring motions,however it remains a challenge of rationalising spatiotemporal chaos to predict theswitching between these processes. For instance simulation B3 shows that intermittent2bursts can raise the local (in time) mixing efficiency, despite the majority of the mixingin this case being controlled by the straining and scouring layered motions characteristicof class ‘l’. Examination of a segment of the trajectory reaching large E in this casesuggests straining and shear instability, saturating at various finite amplitudes, combineto populate the distribution of figure 1(b). Furthermore, as already widely discussed,research into exact coherent structures embedded in turbulent flows faces a seriouschallenge in the extension of methods to handle higher Reynolds numbers and spatiallylocalised dynamics. It remains unclear whether such UPOs can be identified in ‘energetic’stratified turbulence with Re B (cid:38) O (100), yet it is of particular interest to understandwhether the observed reduction of mixing efficiency at such high Re B is associated witha qualitative change in mixing properties. Acknowledgements . We extend our thanks, for many helpful and enlightening discussions,to Paul Linden, John Taylor, Stuart Dalziel and the rest of the ‘MUST’ team inCambridge and Bristol. We also thank the three anonymous referees whose construc-tive comments have significantly improved the clarity of the manuscript. The sourcecode used in this work is provided at https://bitbucket.org/dan_lucas/psgpu andthe associated data including initialisation files and converged states can be foundat
This work is supported by EPSRC Programme GrantEP/K034529/1 entitled ‘Mathematical Underpinnings of Stratified Turbulence’. Themajority of the research presented here was conducted when DL was a postdoctoralresearcher in DAMTP as part of the MUST programme grant.
REFERENCESAshurst, Wm. T., Kerstein, A. R., Kerr, R. M. & Gibson, C. H.
Phys.Fluids (8), 2343–2353. Chandler, G. J. & Kerswell, R. R.
J. Fluid Mech. , 554–595.
Cvitanovi´c, P. & Gibson, J. F.
Physica Scripta , 4007.
Garaud, Pascale, Gallet, Basile & Bischoff, Tobias
Physics of Fluids (8), 084104–22. Ivey, G. N., Winters, K. B. & Koseff, J. R.
Annu. Rev. Fluid Mech. , , vol. 40, pp. 169–184.
Kawahara, G. & Kida, S.
J. Fluid Mech. , 291.
Kawahara, G., Uhlmann, M. & van Veen, L.
Annu. Rev. Fluid Mech. (1), 203–225. Linden, P. F.
Geophys. Astrophys. Fluid Dyn. (1), 3–23. Lucas, D., Caulfield, C. P. & Kerswell, R. R.
Lucas, D. & Kerswell, R. R.
Phys. Fluids (4), 045106–27. Lucas, D. & Kerswell, R. R.
Journal of Fluid Mechanics , R3–11. Maffioli, A., Brethouwer, G. & Lindborg, E.
J. Fluid Mech. , R3–12.
Mashayek, A., Salehipour, H., Bouffard, D., Caulfield, C. P., Ferrari, R.,Nikurashin, M., Peltier, W. R. & Smyth, W. D.
Geophys. Res. Lett.
Mater, B. D. & Venayagamoorthy, S. K.
Geophys. Res.Lett. (13), 4646–4653. Osborn, T. R.
Journal of Physical Oceanography (1), 83–89. Peltier, W. R. & Caulfield, C. P.
Annu.Rev. Fluid Mech. , 135–167. Salehipour, H. & Peltier, W. R.
J. Fluid Mech. , 464–500.
Salehipour, H., Peltier, W. R., Whalen, C. B. & MacKinnon, J. A.
Geophys. Res. Lett. . Scotti, A. & White, B.
J.Phys. Oceanogr. (10), 3181–3191. Shih, L. H., Koseff, J. R., Ivey, G. N. & Ferziger, J. H.
J.Fluid Mech. , 193–214.
Smyth, W. D.
J. Fluid Mech. , 209–242. van Veen, L., Kida, S. & Kawahara, G.
Fluid Dyn. Res. (1), 19–46. Woods, A. W., Caulfield, C. P., Landel, J. R. & A., Kuesters
J. FluidMech. , 347–357.
Zhou, Q., Taylor, J. R. & Caulfield, C. P.
J. Fluid Mech.820