Intermittency Scaling for Mixing and Dissipation in Rotating Stratiﬁed Turbulence at the Edge of Instability

: Many issues pioneered by Jackson Herring deal with how nonlinear interactions shape atmospheric dynamics. In this context, we analyze new direct numerical simulations of rotating stratiﬁed ﬂows with a large-scale forcing, which is either random or quasi-geostrophic (QG). Runs were performed at a moderate Reynolds number Re and up to 1646 turn-over times in one case. We found intermittent ﬂuctuations of the vertical velocity w and temperature θ in a narrow domain of parameters as for decaying ﬂows. Preliminary results indicate that parabolic relations between normalized third-and fourth-order moments of the buoyancy ﬂux ∝ (cid:104) w θ (cid:105) and of the energy dissipation emerge in this domain, including for passive and active scalars, with or without rotation. These are reminiscent of (but not identical to) previous ﬁndings for other variables and systems such as oceanic and atmospheric ﬂows, climate re-analysis data, fusion plasmas, the Solar Wind, or galaxies. For QG forcing, sharp scaling transitions take place once the Ozmidov length scale (cid:96) Oz is resolved— (cid:96) Oz being the scale after which a turbulent Kolmogorov energy spectrum likely recovers at high Re .


Introduction
Much progress has been made in our understanding of turbulence, developing and analyzing experiments, observations, direct numerical simulations (DNSs), as well as numerous theoretical and modeling techniques, among them the direct interaction approximation (DIA) stemming from a renormalization-group (RG) procedure [1], the test-field model (TFM), and related second-order closures, often Markovianized such as the Eddy-Damped Quasi-Normal Markovian (EDQNM), or the quasi-normal scale elimination [2] (see [3][4][5] for recent reviews, and this volume).These flows are complex, due to the large number of eddies interacting non-linearly, as in homogeneous isotropic fully developed turbulence (FDT), and due to their inherent lack of predictability [6][7][8].Flow structures change in the presence of waves, stemming for example from imposed gravity, rotation, or a strong uniform magnetic field and, specifically, their behaviors differ according to how fast the waves are relative to the eddy turn-over time τ NL = L int /U rms with U rms , L int the r.m.s.velocity and the integral length scale.One of the governing parameters of the resulting dynamics, beyond the Reynolds number, is the ratio G = τ P /τ NL with τ P the characteristic period of the waves; G is the Rossby number in the presence of an imposed rotation, or it is the Froude number in the presence of stratification.
Accurate direct numerical simulations with pseudo-spectral methods were developed in particular at MIT and NCAR in the early 70s and were immediately put to use by Jack Herring to study turbulence, as e.g., in the context of the growth of errors in numerical modeling of the atmosphere [8], of two-dimensional (2D), axisymmetric or QG turbulence [9][10][11], for comparing closures to DNS [12,13], including for non-Gaussian statistics and intermittency in the far dissipation range [14,15].Jack Herring was central to this introduction of models and DNS to complement experimental and theoretical studies of FDT, and he then moved on to study stratified turbulence [16,17].
Characteristic of turbulence is the presence of localized extreme dissipation in many settings [18,19].In the ocean, satellite altimetry allows to study the interactions of largescale planetary or (sub-)synoptic-scale waves with shear structures, as observed for example in the Gulf Stream and the Kuroshio Current [20,21].Turbulence leads to filamentation, as seen in phytoplankton dynamics forced by hurricanes [22], as well as to energy transfer and dissipative events (e.g., [23,24]).There is also recent radar data indicative of extreme vertical drafts in the (strongly stratified) polar summer mesosphere [25].Furthermore, waves can change the energy spectra, as demonstrated in weak turbulence [3], detected in the atmosphere [26] and in DNS [27], and as shown in magnetohydrodynamic (MHD) as well [28,29].Atmospheric data can be related to numerical results for strongly stratified flows showing quantitatively that turbulence interacting with fast waves can be more intermittent than FDT [30,31] in a narrow domain of parameters (see [32] for Solar Wind observations, and [33] for energy dissipation in reconnection events leading also to particle acceleration in the magnetotail).
In this context, we present preliminary results using DNS on relationships between normalized third and fourth order moments taken as proxies to characterize the intermittency and dissipative properties of rotating strongly stratified turbulence, contrasting behaviors for high vs. low kurtosis of the vertical velocity with varying Froude and Reynolds numbers (see next Section for definitions).Similar relationships have been found in the turbulent atmosphere and ocean [34][35][36][37][38], climate data re-analysis and glaciology [39,40], and in fusion plasmas [41][42][43], the magnetosheath, the interplanetary magnetic field or the cosmos [44][45][46][47][48]. Equations, definitions and our numerical methodology are given in the next section.Some of these properties for several field variables for the new runs performed herein are presented in Section 3 for Navier-Stokes and stratified flows.Rotating stratified turbulence is analyzed in somewhat more detail in Section 4.An overview is given in Section 5, and Section 6 presents a discussion and conclusions.

Equations and Definitions
We perform several sets of numerical simulations for an incompressible (∇ • u = 0) velocity field u = (u, v, w) = (u ⊥ , w) with rotation and stratification.The active scalar θ (called in the following temperature) is in velocity units; strictly speaking, θ and the temperature have opposite signs, but they are linearly related through a thermal expansion coefficient.With ν, κ the kinematic viscosity and thermal diffusivity, taken equal, p the pressure, N = [−g∂ z θ/θ 0 ] 1/2 the Brunt-Väisälä frequency, and f 0 = 2Ω, Ω = Ωz an imposed rotation in the vertical direction of unit vector z , the Boussinesq equations are: Two types of forcing for F u , F θ have been implemented, both set in the large scales at a wavenumber k F ≈ 2.5.The first one has random phases and is white noise in time.The second forcing, for the rotating stratified turbulence (RST) case, ensures a quasi-geostrophic balance at all times in the forcing modes (see [11]), following the procedure described in detail in [49] for QG initial conditions with zero vertical velocity.Pressure balance with the Coriolis force in the horizontal is enforced, as well as hydrostatic momentum balance due to gravity (see [50] for a different implementation in a study of balanced dynamics).One might note that the time-scale associated with the forcing, in particular in the case of white noise, may matter in the interactions between waves and eddies, with more waves for short correlation times, in particular at low Froude number for which N becomes comparable to 1/∆t, where ∆t is the time-step [27].The buoyancy flux B f , the kinetic energy dissipation du and the potential vorticity P V (point-wise invariant in the absence of dissipation) are: with ω = ∇ × u the vorticity.B f is an energy exchange term with the same physical dimension as du .It leads to local changes in density in the presence of gravity waves which need to be modeled in large-eddy simulations.Characteristic scales can be defined as the Taylor micro-scale λ T that factors in the dissipation in the inertial range, the dissipation scale η K based on a Kolmogorov spectrum, the buoyancy scale L B measuring the effect of stratification in the large scales, and the Ozmidov scale Oz : The last expression indicates the multi-scale link in the dynamics of the flow between the development of turbulence and its stratification.The dimensionless parameters governing the fluid behavior are the Reynolds, Rossby and Froude numbers, Re, Ro, Fr, with the Prandtl number Pr ≡ ν/κ = 1 here: with R λ , R B the Taylor and buoyancy Reynolds numbers.For the purely stratified runs, the forcing scale L f = 2π/k F is used in the evaluation of the parameters instead of the integral scale L int .One also defines the interaction parameter R IB and, when involving gravity through N, the Richardson number Ri based on an overall vertical gradient of u ⊥ , as well as a local gradient Richardson number Ri g : , I Ri = Ri g .(6) The phenomenological evaluation of kinetic energy dissipation, D = U 3 rms /L int , is based on expressing that it occurs in an eddy turn-over time, irrespective of the (small) viscosity.When du becomes comparable to D , as in FDT, one has R B = R IB .Note also that Ri, Ri g , in the absence of uniform shear, are based on large-scale shear resulting from the overall nonlinear dynamics.One can also define a Taylor Froude number F λ and relate it to the other parameters.After simple manipulations, one has: The RST runs cover the three regimes identified in [51]: wave and viscous-dominated at small [Fr, R B ], eddy-wave interactions at intermediate values, and the eddy-dominated regime with a resolved turbulent inertial range beyond the Ozmidov scale ( Oz > η K ).These three regimes can likely be related to the classification proposed in [52] for the nocturnal planetary boundary layer (PBL) into weakly stable, transitional, and very stable regimes: when increasing the Froude number, there is a drop in the Richardson number, particularly marked in the transition regime, and a sharp increase in eddy diffusivities.
Another classification of stratified flows can be found in [53] (see also [54,55]) in which the role of viscous dissipation, governed by a small Reynolds number, is particularly emphasized.This is in contrast to the role of fast waves compared to nonlinear eddies, as governed by a small Froude number; the competition between these two effects can be encompassed in the parameter R B = ReFr 2 , or in R IB , but not necessarily uniquely [56].At low R B , anisotropy comes strongly into play.Previous studies (see e.g., [57]) have shown that when the system is forced three-dimensionally, this results in a significant excitation of wave modes which, as a result of viscous effects, can excite less vortical modes.However, their interactions with waves are still present [58], and we call this regime wave dominated since no turbulence range is excited.We also note that, at a given (high) buoyancy Reynolds number, all three regimes can be identified in the energy spectra in a given flow, provided L B , Oz and η K are resolved (see [59] for an analysis of these regimes using a turbulence closure).Different values for R IB can be found when using the energy injection rate F u • u instead of du , the two equilibrating only on average.This is particularly relevant when there are strong intermittent fluctuations in the system [60,61], and this balance can also be broken in the presence of an inverse energy cascade.Finally, the skewness and excess kurtosis of a random variable V, with averages taken over three-dimensional (3D) space and with S G V = 0, K G V = 0 for a Gaussian, are defined as usual, with the following parabolic relation to be tested in this paper, e.g., for V = B f :

Overview of the Direct Numerical Simulations
We use a variety of DNS (see Table 1).The runs labeled H are for the hydrodynamic case (Navier-Stokes, N ≡ 0, Ω ≡ 0); the S i=1,6 runs deal with the Boussinesq equations with Ω ≡ 0. The R, Q runs are RST with random or QG forcing; all these runs have N/ f 0 ≈ 5, leading to 0.28 ≤ Ro ≤ 3.7.Smaller Rossby numbers would imply the development of a strong inverse energy cascade (see e.g., [62] and references therein), a cascade that would be altered by the large-scale forcing when there is insufficient scale separation as here.The dependence of kurtosis on skewness for various fields in the presence of both inverse and direct energy cascades is left for future work.
The linear dimension of the cubic grid is n p , varying from 128 to 512.The code is a versatile pseudo-spectral formulation of the time-integration of a large set of partial differential equations for fluid and plasma turbulence, with efficient hybrid parallelism implemented in a periodic box [63].There is a version of the code in non-cubic geometry, as well as with non-periodic boundary conditions in one direction for the incompressible case [64].T M is the extension in terms of turn-over time τ NL of the segment of the simulation used to build the statistics, with 41 ≤ T M ≤ 1646 across all runs.Data is gathered roughly four times every τ NL , but for purely stratified runs, statistics of the dimensionless parameters are taken around the first peak of dissipation in order to avoid the effect of energy accumulation in the slow-mode leading to vertically sheared horizontal winds [65,66].An estimate of the error on the skewness, and thus of its departure from Gaussianity, is √ 6/N S , and twice that for the kurtosis [38,67], N S being the number of independent samples.For our runs, the worst-case gives S > 0.38 for a flow to be determined to be non-Gaussian, and ≈0.06 in the best case.All Taylor Reynolds numbers are above ≈27 and should thus suffice to determine small-scale intermittent properties of the flow [68,69] when properly averaged over long times.In fact, local quasi-singular behavior in the intensity and curvature (morphology) of the vorticity and strain tensor can already be discerned at moderate R λ when accurately evaluating the geometrical properties and anisotropy of the flow close to such structures through a so-called sparseness scale [70][71][72].
Finally, in Table 1, K M w is the maximum over time of the kurtosis of the vertical velocity, a signature of large-scale intermittency and strong dissipation in the transitional regime [30,31,73,74].Most of the RST runs are at low to intermediate R B with high K M w , as for the stratified runs.How well the runs are resolved is measured through Res N = k max η K , with k max = n p /3 using a standard dealiasing method.Res N ≥ 1 implies that the dissipation scale is resolved by the numerical grid; a, b are the parameters when performing aS 2 + b fits analyzing K(S) for the buoyancy flux; r.m.s.errors are between 1.8 and 6.3% for all runs for a, and between 0.8 and 2.3% for b, with respective averages for the eight QG runs of 4.27% and 1.44%.The extent of the turbulent range beyond the Ozmidov scale is Res T = Oz /η K = F 3/2 λ ; when greater than unity, a sharp transition occurs in the K(S) scaling parameters for both B f and the kinetic energy dissipation du (see Section 4 below).
Table 1.Characteristics of the hydrodynamics (H), stratified (S) and RST runs with random (R) or QG (Q) forcing.K M w is the temporal maximum kurtosis of vertical velocity, T M the final time in units of τ NL , and is the extent of the fully turbulent range, and I Ri = Ri g is the averaged gradient Richardson number.Res T and I Ri mark a clear scaling transition.Many studies of stratified flows have been performed [27,53,54,[75][76][77][78][79][80].Small-scale intermittency is present once Oz /η K > 1, whereas at large-scale it is observed in a restricted range of values of Ri g strongly peaked for shear instabilities (as well as in a narrow range of Fr, R B , Res T ), without rotation [30], or with it as shown in this paper.An increase in the amplitude of nonlinear P V when going to higher buoyancy Reynolds number is also observed [74,75].Note also that, in the presence of a strong Coriolis force, a quasigeostrophic balance is obtained that can lead not only to an inverse energy cascade to large scales (Ref.[62] and references therein), but in fact to a dual constant-flux energy cascade to both large and small scales [81][82][83].

Large-Scale and Small-Scale Intermittency in Turbulence
Intermittency, measured by non-Gaussian statistics [14], can modify passive scalar advection such as pollutants, and thus the chemistry within turbulent flows, e.g., through segregation as in the convective boundary layer, in particular when strong shear is present [84] (see [48,85] in the interstellar medium and galactic contexts).These baroclinic instabilities can be abrupt and result in turbulent bursts, reminiscent of large-scale on-off intermittency (LSI).There are in fact several such findings of LSI, starting with the probability distribution function (PDF) of the temperature which can be non-Gaussian (e.g., Métais and Herring [16], Rorai et al. [73]).One of the striking new results in [31] is that dissipation can be stronger locally than for FDT, as quantified by measuring the minimum flow volume needed to reach some level of energy dissipation, say 50%.Close to the peak of intermittency of w, θ, this volume can be as small as 12%, compared to 15% for FDT at a similar Reynolds number, and less than 4% for the active scalar, possibly indicating than a log-normal distribution is not adequate for the statistics of these flows [86] (see [87] for the FDT case).The waves can enhance dissipation in hot spots even when, globally, the flow is constrained by the large-scale, smooth dynamics.Extreme events are observed frequently in the atmosphere, e.g., in terms of w, θ [88], of the rate of energy dissipation [89] as well as for extreme precipitation [90,91].A possible role in the formation of such structures is through the precession resonance of Rossby waves which was analyzed e.g., in [92,93], leading to their breaking as well as to their interactions with potential vorticity [94,95].Similarly, in the ocean, wave breaking leads to efficient mixing processes [77-79], and observations of double-signed skewness can be the signature of a front, modeled in a dichotomic way for relative vorticity [67] as for the Agulhas current.
A simple non-Gaussian model is to express the PDF for field V as the sum of a normal distribution and its square, the relative weight of the two being governed by an open parameter.When the quadratic term dominates, one finds . In fusion plasmas, similar relationships are obtained [41,42,97,98], with different quadratic laws followed by different experiments [43].The properties of the high-variability parts of turbulent flows, such as shear layers and their relaxation, seem universal and can be described by the model in [34], such shear layers developing at the boundaries [97], or as internal structures due to the turbulence.In fact, high kurtosis for the velocity and temperature at the edge of a shear layer have been found in a DNS of a Kelvin-Helmholtz instability [99], and singular shear layers can lead to intrinsic stochasticity [100].Moreover, interactions with large scales in fluid turbulence represent roughly 20% of the kinetic energy flux at moderate Re.For structure functions beyond the energy level, including in the absence of waves, the nonlinear interactions are found to be mostly non-local, involving the integral scale [101], but with a tendency towards locality as Re → ∞.Large-scale extreme events can arise as well in climatological studies, for example from quasi-resonant amplifications of planetary waves [92].One also recovers a quadratic K(S) law by adding three Gaussian PDFs, a model suitable for climate temperature variability [102], with the result that knowledge of skewness becomes a predictor of kurtosis.In fact, (S, K) mapping can be taken as an indicator of the turbulent dynamics in all these flows [44,46,67,96].

Is There a Skewness-Kurtosis Relationship for Fluid Turbulence with a Passive Scalar?
Curiously perhaps, FDT (with a passive scalar) has barely been analyzed in the context of K(S) laws.This problem was studied in Herring and Kerr [12] in the decaying context and for Pr = 0.5.The skewness of the velocity was compared with experiments and with predictions of turbulence closures at R λ comparable to what we compute in this paper, but here in the forced case.The flatness evaluated from the DNS was also contrasted with that given by closures (DIA and TFM).Good agreement was observed up to R λ ≈ 30.At high R λ (above 800 or so), the experimental and numerical data reviewed in [103] for fluid turbulence does lend some credence to a relationship K ∼ S j , j 2.5 for the velocity derivative, as noted in [43] (see also [104], Figure A6, for a recent estimate of K, S for numerical and experimental fluid turbulence with R λ in excess of 6000).Furthermore, a statistical study of passive scalars in the atmosphere under many physical conditions gives a parabolic relationship, with in some cases rather high values of both S P and K P [36,37].Several PDFs have been invoked to justify a quadratic relationship between K and skewness S, depending on the tails of such distributions [105], and exponential wings in PDFs of the horizontal vorticity for RST flows have been found [74].Thus, we briefly discuss first the kurtosis and skewness data for FDT and a passive scalar with Pr = 1.
Figure 1 (top) displays K(S) for run H2 for the vertical velocity, temperature and ω z (a-c).At the bottom, K(S) is shown for the buoyancy flux B f (we are not considering here the joint w, θ distribution), for runs H1 (d) and H2 (e), with least-square fits for [a, b] given in Table 1.Note that, for a unimodal and symmetric distribution, one can derive estimates that improve on the Cauchy-Schwarz inequality, namely K ≥ S 2 − 6/5 [106], thus limiting potential closures [107].In all the K(S) figures, two parabolae are drawn, K(S) = 1.5S 2 − 1.1 and K(S) = 3(S 2 + 2) (using n = 5 in [108], Equation (9b)).The first formulation appears in a variety of data, as described earlier.It should also be mentioned that it is quite close to the approximation developed in [109] (who also considers joint distributions), namely K ≈ 16S 2 /9, emanating from an expansion of a PDF assumed to be quasi-Gaussian and first analyzed on data of coastal waves [110] (see e.g., [111] in the context of the Weibull distribution for a model of sea surface wind velocity, and [112] for recent data analysis of strong waves on a beach with a sloping bottom).The second K(S) formulation arises from the analysis of stochastically-generated skewed distributions (SGS) for which general K(S) relations are obtained.These SGS distributions can be related to Markov processes, thereby shedding some light on the underlying stochastic systems [108].As expected, no measurable intermittency is seen in Figure 1 in the large-scale field w, with both S and K close to zero, but θ, ω z display measurable non-zero values that are at times higher than those given in Herring and Kerr [12], the forcing likely allowing for wider excursions.In plasmas, the K ∼ S 2 law obtains for both quiet and active regions and, in that sense, the relationship itself may not be related to the development of strong turbulence, but rather to intrinsic multi-scale properties of the fields themselves [41]; du (not shown) is also close to a quadratic law, but with higher coefficients.Taking ω z as a (skewed) velocity gradient, the values of S, K are in fact consistent with data in [113], but higher Reynolds numbers will have to be explored to see whether a parabolic law actually appears.Many models describe the interactions between widely separated scales from energycontaining to dissipative eddies.The combination of (weak) fast additive noise and largescale multiplicative noise is known to produce PDFs with power-law tails [114].Further examples are on-off intermittency in FDT close to instability [115], the scale-dependent creation of intermittency [116], or the magnetic dynamo problem [117].Chaotic dynamics can also be seen as quasi-linear behavior forced by non-Gaussian processes whose events in the tail can trigger bifurcations in the system [118].In this context, let us mention the climate model of Hasselmann [119] which begins with a linear stochastic (Langevin) equation for the temperature anomaly Θ, namely ∂ t Θ = −λΘ + η; here, λ represents the heat flux and η is a noise term modeling small-scale fast fluctuations.One can in fact write λ = λ + λ with λ an average, making the noise multiplicative [38] (see [120] for a recent review).Indeed, small-scale eddies are in large numbers, have low energy compared to that of the large scales and are intermittent, leading to complex behavior.This makes it plausible to model them through a stochastic linear noise acting on the slower large scale eddies (see also [121] for another Langevin approach, and [122] for stochastic turbulence models).One can then derive the parabolic relationship in Equation ( 8) with the expectation that a ≥ 3/2, b ≥ 0 [38], a finding corroborated by several analyses of climate data.The long-term (monthly onward) climate statistics are Gaussian, but the daily averages reflect the stochasticity of small-scale eddies [38] which also interact with waves.Therefore, on some intermediate temporal scale for the large structures, these small eddies can be viewed as noise, although there exist nonlocal (in scale) interactions between small and large scales.This leads to enhanced dissipation at the level that is needed in order that, in its adimensionalized form and at least in the absence of waves, it be of order unity, both for fluid turbulence [19,123] and for fusion and space plasmas [124,125].

K(S) Relationships for the Buoyancy Flux in the Purely Stratified Case
For strong stratification, the vertical velocity intermittency can be significant, like in the nocturnal PBL [18,34,52], or for DNS in a narrow range of Froude number [30].Does a K(S) relationship ensue in that regime?To answer this question, we examine the results of a series of purely stratified runs with 0.032 ≤ Fr ≤ 0.145, and with high Re ≈ 3000.At left in Figure 2 is K(S) for w for the run with the highest kurtosis of w [30,31].We observe that there is a large number of data points with large excess kurtosis, but the skewness of w remains close to zero.On the other hand, such is not the case for the buoyancy flux B f shown at right in Figure 2 for all the S runs.In the inset, three parabolae are given (dash; dots; and dash-dots) corresponding to different [a,b] fits.The first remark is that, for stratified turbulence, a parabolic K(S) law seems to emerge as well for B f (and for du , not shown).Note that practically all runs lie on the lower parabola marked with dash, and showing as an orange cloud of points in the plot.Superimposed on this, for two runs close to the maximum of K w , there is an explosion in the values taken by the kurtosis (blue triangles and black squares, runs S3 and S4).This confirms the analysis performed in [31] that dissipation in this Fr range takes high local values; it corresponds to a transitional phenomenon with high intermittency of vertical velocity, flux and dissipation itself, and thus leading to enhanced mixing as we move from regime I (wave dominated) to regime II (eddy-wave interactions), likely similar to what happens in the nocturnal PBL [52].The a coefficient in front of the quadratic term in the K(S) law for the vertical velocity in the PBL is thought to depend on the stability of the layer through Ri g [54,126].Indeed, strong vertical drafts are clearly correlated with Ri g ≈ 1/4 (see [30], Figure 3).It may be that K(S) laws are emerging here from the flow being at the edge of instability [127,128].An analysis similar to that of Vieillefosse [129] yields nontrivial relations between high-order moments of the field gradients near the onset of convective instabilities [130].An inverse scaling, on more than two orders of magnitude, between two a-priori unrelated parameters, namely one governing stratification (the gradient Richardson number) and one governing the turbulence in the presence of stratification (the buoyancy Reynolds number)-was observed in DNS of rotating stratified flows [51].One could thus infer that the a parameter just mentioned would depend in a similar fashion on R B .This inverse Ri g ∼ R −1 B relationship is significant because it implies a strong dynamical effect of vertical shear, and a balance between vertical gradients (of the turbulent horizontal velocity field) and the turbulent isotropic eddies beyond the Ozmidov scale [51], as hinted in fact by several authors over the years in varied contexts.

The Case of Rotating Stratified Turbulence (RST)
Intermittent vertical velocities are also found in RST, for decaying flows [74] and, as shown here, in forced runs.We look first at global dynamics.In Figure 3a is given as a function of time the total kinetic energy for run H1 as a base, and for three runs with varying Fr and 0.36 ≤ Ro ≤ 3.7.The higher level of saturation for the RST runs is due to the fact that on the one hand, there is a lesser dissipation due to the strong effect of waves, and on the other hand, Ro ≈ 0.36 for run Q1, which is low enough that an inverse transfer of energy to the large scales is plausible, although it cannot be distinguished directly on energy spectra which are still peaked at k ≈ k F (see Figure 3d).Higher Ro may correspond to the environment of tropical storms in which turbulence also plays a dual role [131,132].In Figure 3b we plot K w (dotted magenta line), the kinetic (solid blue), potential (black dashed) and total energy dissipation rates for run R1, with the time expressed in units of the turnover time τ NL .We see that K w can take large values, as in the purely stratified case [30].There are also discernible excursions in the dissipation rates.The high-K values are associated with local instabilities and overturning regions; these excursions in K w are also present in several other runs.We recover here the behavior analyzed in detail in [30,31] for purely stratified flows: narrow peaks in K w as a function of Fr and R B occur, compatible with observed strong updrafts in the DNS.The temporal evolution of the kinetic energy E V and enstrophy |ω| 2 are given for the total length of that run in (c), showing a persistence in the amplitude of fluctuations.The sporadic energetization of the flow has its signature in energy spectra shown in Figure 3d for run R1; indeed, the kinetic and potential spectra gain close to a factor of ten at the smallest scales for an excited state (with a peak in K w ) when compared to quiet times (with a trough in K w ), whereas the spectra have comparable amplitudes at the largest scales.Similar differences in small-scale spectra occur for the QG runs.Note also that, the forcing being centered in the large scales, no inverse cascade of energy that would occur in the presence of strong rotation can be observed in these runs (see [81,82] for DNS resolving both the inverse and direct energy cascades).In Figure 4, we give K(S) for w for runs R1-R3 (a-c) , and for B f for runs Q1, Q3, and Q8 (d-f), with buoyancy Reynolds numbers varying from 1.4 to 385.At top, the skewness is typically bounded by ±1 or less, and there is a clear decrease in the kurtosis of the vertical velocity as R B increases from left to right.We note that, for all these runs, and similarly to the FDT and purely stratified cases, the vertical velocity has basically no skewness but its kurtosis can take moderate values, especially so at intermediate R B (here, the first plot), although no parabolic law emerges from the data.For the buoyancy flux, there is a marked transition towards a standard scaling (≈ 1.5S 2 + b), once the Kolmogorov regime is sufficiently resolved, i.e., for high enough R B , Ri g , and with a classical FDT dissipation taking place (Figure 4e,f).Note that the intermittency uncovered in this paper likely deals more with intrinsic properties of dissipation than local accelerations that might be observed in a boundary layer.A different approach, following [109], is to expand the PDF in terms of its cumulant under the assumption of closeness to a Gaussian distribution.In the simplest case, a relationship is obtained, namely K = [4S/3] 2 , close to the previous expression (a ≈ 1.5 for small b).The K(S) plots for the temperature (Figure 5a,b) and the vertical component of the vorticity (c,d plots in Figure 5) are similar as well, as shown for runs Q2 (a,c) and Q7 (b,d).In all cases, the skewness stays close to zero and there is a decrease of kurtosis as the flows become more dynamical.
A last issue is how much mixing occurs, estimated through B f , du taken as proxies [79,86].Thus, we give in Figure 6, K(S) for du for runs Q1, Q2 (a,b) with Res T ≈ 1, and Q5, Q6 (d,e, Res T > 10).Again, a transition is observed, with steeper scaling in the latter case; similar results obtain for FDT, with a best fit for du (run H2) of K = 2.6S 2 + 7.2 (not shown).Note that, for the components of the strain tensor itself, S ij = [∂ i u j + ∂ j u i ]/2, we have K ≥ 21S 2 /4 in the incompressible case [133].Figure 6 (right) gives the PDFs of du for run Q3 for times with a high (c) and low (f) K w ; a ≈ −3 power law is plausible, with a larger extension when the intermittency of w is stronger; note that a shallower scaling for du has been observed in oceanic data, together with substantially higher values of K, S [134].In the fits, α = 1 + S 2 du /4 (see [108]).

Possible Frameworks for K(S) Relationships
It is not clear what the link of the above results is with the analysis of complex turbulent flows using the concept of avalanches and of self-organized criticality [127,135].In that latter case, the dissipative structures at small scale propagate in some fashion their intermittency to the large scales, just as small avalanches become large through merging.Saying it another way, dissipative structures near strain and vortices, or in current sheets and filaments, exist at the large scale characteristic of the flow and at the small scale corresponding to their thickness, of the order of the dissipation length.Therefore, the governing parameter may not be Re or R λ , but could be related to the ability of the small-scale structures to organize into larger clusters, or in other words to a transition at the local (in scale) Reynolds number of unity in the vicinity of the Kolmogorov scale in fluid turbulence.The occurrence of a quadratic K(S) law, albeit with varying coefficients, has been noted before.One should analyze large-enough samples [108,136]), but the non-Gaussianity of turbulence is ascertained, for non-vanishing dissipation, by the exact laws stemming from the conservation of quadratic invariants, as for kinetic energy in fluid turbulence [137,138], under several conditions: with u L the longitudinal velocity along the distance |r|, and δ expressing velocity differences.In MHD, cross-field correlations emerge [139] (see [140] for a review), which are central to the structure of reconnection events [33,47,141].Inertia-gravity waves, which exchange energy between the kinetic and potential modes, will not alter the balance at the level of total energy dissipation.However, when writing an equation for the temporal evolution of B f , it occurs at the scale of the turn-over time since it deals with large-scale variables, whereas the equipartition deficiency N [θ 2 − w 2 ] due to waves evolves on a fast time for small Fr and can be averaged to be modeled as a stochastic forcing for B f , du ; this makes their temporal evolutions akin to systems already discussed herein (see e.g., [38,119]).Another remark is that, obtaining quadratic K(S) laws for both weak and strong turbulence, diffusive or not, and dispersive or not, could lead one to conclude that the issue is more of a mathematical constraint than a physical one.However, the K(S) relationship is important in limiting possible closures applied to turbulence, as for example in the context of solar convection (see [107] for recent developments of closures beyond the quasi-normal framework).It has also been remarked that these systems have in common a convective nonlinearity, however weak or intense [97].The specific transition in the scaling of the K(S) laws obtained here is linked to the fact that waves become progressively more strongly coupled to nonlinear eddies, or in other words steepen significantly, thus leading to enhanced dissipation and to the onset of a Kolmogorov energy spectrum at small scale.The existence of K ∼ S 2 laws at high and low R λ leads to another possible origin for it, namely the presence of intermittency in the dissipation range common to all such flows.As pointed out in [15,142], this range can be prone to extreme intermittency because of a faster than algebraic decrease of the intensity of turbulence at high wavenumber, a spectral domain that exists even in the near absence of a nonlinear inertial range at low R λ .This likely leaves only a few shear structures dominating the overall statistics.In this context, studies of Lagrangian turbulence have been able to show the link between the dissipation range and large scales, with both short and long-range correlations [143].
Computations performed at moderate R λ identify small-scale intermittency with numerous large-scale structures [15], as opposed to the more familiar vortex filaments [144][145][146].Such nonlocal structures could be detected in rotating stratified flows using specific algorithms (see [135] for the identification of current sheets at the onset of the dissipation range of plasmas in the MHD regime).In that spirit, it would be of interest to consider as well the case of quasi-geostrophy, that of cubic nonlinearities like for compressible flows, or for solitons (such as for the Korteweg-deVries or the cubic nonlinear Schrödinger (NLS) equations and their multi-dimensional extensions), all in the presence of forcing and dissipation (see also [147] where the K(S) behavior for the NLS equation is associated with on-off intermittency).The large-scale behavior linked with the presence of Rossby waves and their possible nonlinear coupling and breaking may also reinforce or even alter large-scale intermittency.
Finally, a link between inertial and dissipative ranges is consistent with a linear stochastic model for the fast small scales which is successful in giving K(S) ∼ S 2 for climate data [38,39].In decaying experiments at both low and high R λ using sulfur hexafluoride, the inertial range does not follow a strict power-law, but rather may display a logarithmic correction that is independent of Reynolds number and that may arise from the dissipative range, putting into question their independence [148].Nonlinearities do not seem to be central to the establishment of non-Gaussian statistics beyond their (necessary) mode-coupling role [149], and it is rather the susceptibility of the system being close to criticality, linking all scales through dissipation that is central to turbulence dynamics and intermittency.This can manifest itself as avalanches in granular media, or in the development of shear instabilities in (e.g., RST) flows close to the critical Richardson number ( [74,75,89,127,150] and references therein).

Discussion, Conclusions, and Perspectives
In the nocturnal planetary boundary layer, " ... turbulence is expected to become more intermittent in very stable conditions" (Mahrt et al. [52]).Relating to these observations, we are likely dealing in this paper with two phenomena, possibly superposed and intertwined.On the one hand, we have the increased intermittency, as diagnosed by a high kurtosis of vertical velocity and temperature, if in a narrow domain of Froude and buoyancy Reynolds number [30,74], and signaling the onset of small-scale turbulence interacting with the larger-scale (inertia)-gravity waves.The signature in physical space consists of a high local vertical velocity and hence a high local Reynolds number, allowing for instabilities to develop until the system evolves (in a turn-over time) back to a state dominated by the fast waves [31,127,130,151], and the cycle can restart.Thus, it leads to strong local energy dissipation in an otherwise wave-dominated rather placid system, and as such these isolated hot spots are a critical factor in controlling the overall dissipation.This dynamics can be modeled e.g., through the temporal evolution of the velocity and temperature gradient matrices, taking into account the anisotropy of the system and using a closure for the pressure Hessian [130].High kurtosis values for the vertical velocity are narrowly centered around a peak in Froude number (and in R IB ), although it is not known at this point how much the peak depends on the Reynolds number itself.They can be associated for example with jet meandering in the ocean and modeled through stochastic parametrization.In fact, it has been hypothesized -and studied in depth over the yearsthat shear flows are central to the critical instability of turbulence, a point further developed recently in a global fluid context [152][153][154][155], and applied to stratified fluids as well [127] (see [128] for RST, and [156] for plasmas).This transition could also be linked to a rolechange in potential vorticity P V at small scales due to the breaking, at the Ozmidov scale onward, of gravity waves [80].
On the other hand, different K(S) laws emerge with different coefficients for the buoyancy flux, the kinetic energy dissipation as well as for mixing as measured through the ratio B f / du .This second phenomenon of a transition between regimes for K(S) is new to our knowledge, certainly in the specific context of rotating stratified flows, but of course long computations at higher Re, R IB with a resolved range beyond Oz are needed.The scaling corresponds to a shift, in the small scales and possibly as well in the early dissipation range, in the nature of an intermittency which becomes stronger as the nonlinear eddies get stronger, and with a K(S) scaling for B f which seems to approach that found for other turbulent systems [38,41,46].This is linked to the well-known transition at the Ozmidov scale of an energy spectrum becoming progressively more Kolmogorov-like.Many models have been proposed for these K(S) parabolic laws, as partially reviewed all along in the paper.More work will be needed to sort out the different possibilities, and to reach higher buoyancy Reynolds numbers corresponding to conditions pertaining to those in atmospheric and oceanic flows.In order to obtain dissipative structures at higher R B for long-time runs, a computation in a non-cubic geometry will help tremendously.For example, it was shown in [157] that in a flat box with a 1:8 aspect ratio, the process of front and filament formation can be very efficient in sheared flows [130].The evolution of extreme events in such flows with comparable geometries was also considered in [158] for their effect on the mixing of neutrally buoyant particles (see also [159] for a substantially higher aspect ratio together with using hyper-viscosity).Similarly, peak vertical vorticity can be close to two orders of magnitude larger than the imposed rotation, as shown in [160] within a Large Eddy Simulation (LES) framework.Another feature of these flows is their persistent anisotropy [59,128,161,162], associated, for example, with seasonal oceanic variability [163].Thus, the precise formulation of the small-scale parametrization and of the dissipation of resolved scales in numerical codes may affect the mesoscale dynamical properties such as energy spectra [164].Related evaluations for plasma disruptions and fusion, from MHD to Hall-MHD and beyond, might also be of interest to understand betterand possibly control to some extent-these strong small-scale disruptive features.
One of the several models of dissipative events leading to a quadratic K(S) relationship not yet mentioned here is in the context of so-called 1/ f n noise as a signature of longtime correlations in turbulent hydrodynamic and MHD flows [165] where it is shown that the transition from Gaussian statistics in the inviscid case to a non-Gaussian one in the presence of dissipation can be associated with varying properties of energy transfer through scales.Are there several different classes of universality, though?In this context, we have already noted that, for the Navier-Stokes equations, there is a clear departure from standard coefficients for the law K(S) = 3S 2 /2 − 1.1, as well as for low Froude numbers in the case of stratified turbulence and RST.
Another potential issue is the form of the dissipation operator.When using either eddy viscosity or hyper-viscosity as a model , or LES, a sub-grid scale algorithm, or the α-model Lagrangian description for fluids and MHD, do similar quadratic K(S) relations arise?What about for the cyclostrophic balance in helical stratified flows [166]?Or in the ideal (ν = 0) case in its pseudo-dissipative range stemming from the eddy viscosity due to the small-scale equilibrium range [167], or in 2D MHD or even in 1D (allowing for the long temporal integrations which are needed in this intermediate inviscid but turbulence-like phase)?In other words, does modifying the way energy is dissipated in turbulence alter the K(S) relation for velocity or dissipation, therefore suggesting that it is a crucial feature for dissipative events?We already know that the PDFs of velocity gradients in the NS case are affected [168], but it is not clear whether the K(S) law would be.Such renormalized dissipation arise at high Re, whereas a K ∼ S 2 relation is present as well at moderate Re, for which the viscosity is directly active.If a change in the K(S) law occurs, beyond the one identified in this paper and related to the onset of a Kolmogorov range beyond Oz , what functional form does it take, and are there consequences for accuracy in an evaluation of the large-scale dynamics that weather and climate codes seek to reproduce faithfully while using techniques such as LES (see, e.g." the discussion in [108])?
We have presented preliminary results, limited in Reynolds number, and much remains to be done.One issue is to examine the dual constant-flux energy cascade when forcing is at smaller scale [81,82].Other questions arise, such as the persistence of large-scale intermittency of vertical velocity in RST in the presence of a resolved turbulence range, or of strong helicity (see e.g., [169,170] on intermittency and helicity in FDT).Nonlinear interactions ∝ u × ω are weakened in vortex tubes, potentially allowing waves to play a more prominent role in driving the flow closer to the critical state in which strong vertical drafts and ensuing dissipation occur.A wave-vortex decomposition, such as that studied e.g., in [171,172], could be of help to distinguish between these different regimes.One can also ask how anisotropy in RST flows is affected by these strong updrafts.Large-scale intermittency has also been linked to the role of helicity (velocity-vorticity correlations) [173], and the eigenvectors of strain can play an important role as well in intermittent energy transfer, as shown in a simple model in which alignment of vorticity with strain eigenvector corresponding at first to its intermediate eigenvalue, and featuring the combination of a shear flow with a horizontal straining field [174] (see also [175]), and yet anisotropy is essential in possibly stopping the formation of singular structures [70].
We wish to end this text with a direct quote from a review article on waves and coherent flows in the tropics (Stephan et al. [176], see page 2619): "Theory provides the qualitative framework-the causal narrative-within which the quantitative information from models and observations can be meaningfully compared and usefully combined."The research of Jack Herring was on the idealized, theoretical side of the study of the dynamics of the atmosphere, pioneering the inclusive use of theoretical and numerical approaches.As such, it has provided invaluable information, and causal accounts, for the idealized behavior of turbulence and the central role played by nonlinear effects in such-and other geophysical, astrophysical as well as engineering, chemical, biological and climatological-complex systems.

Figure 2 .
Figure 2. (a): Kurtosis vs. skewness of the vertical velocity for run S3 at the peak of K w [30,31].Note the narrow range of skewness allowing to focus on small S w , high K w values.(b): K B (S B ) for the buoyancy flux B f for all stratified runs; parabolic fits and color symbols are given at right.

Figure 3 .
Figure 3.Time variation of (a) kinetic energy for runs (see inset) up to t/τ NL ≈ 100; (b) kinetic, potential and total dissipation (scale at right), and K w , run R1; (c): kinetic energy (black) and enstrophy (red, scale at right) for the total length of Q3. (d): Kinetic (solid black/dash blue) and potential (solid red/dash green) energy spectra at resp.maximum and minimum kurtosis of w, run R1.

Figure 4 .Figure 5 .
Figure 4. K w (S w ) of vertical velocity for RST runs (a) R1, (b) R2 and (c) R3.K B f (S B f ) of buoyancy flux for QG runs (d) Q1, (e) Q3 and (f) Q8.Insets give best fits indicated by dashed red lines.