Correlation between Buoyancy Flux, Dissipation and Potential Vorticity in Rotating Stratified Turbulence

We study in this paper the correlation between the buoyancy flux, the efficiency of energy dissipation and the linear and nonlinear components of potential vorticity, PV, a point-wise invariant of the Boussinesq equations, contrasting the three identified regimes of rotating stratified turbulence, namely wave-dominated, wave–eddy interactions and eddy-dominated. After recalling some of the main novel features of these flows compared to homogeneous isotropic turbulence, we specifically analyze three direct numerical simulations in the absence of forcing and performed on grids of 10243 points, one in each of these physical regimes. We focus in particular on the link between the point-wise buoyancy flux and the amount of kinetic energy dissipation and of linear and nonlinear PV. For flows dominated by waves, we find that the highest joint probability is for minimal kinetic energy dissipation (compared to the buoyancy flux), low dissipation efficiency and low nonlinear PV, whereas for flows dominated by nonlinear eddies, the highest correlation between dissipation and buoyancy flux occurs for weak flux and high localized nonlinear PV. We also show that the nonlinear potential vorticity is strongly correlated with high dissipation efficiency in the turbulent regime, corresponding to intermittent events, as observed in the atmosphere and oceans.


Introduction
Atmospheric and oceanic flows are complex. The rotation of the Earth together with density stratification at large scale give rise to inertia-gravity waves, and such flows also support turbulent eddies that interact nonlinearly with the waves. The energy input, mostly at large scale, comes from a variety of sources such as solar radiation, surface winds and tides, and this energy has to be dissipated at small scales; but how does it get there? In fact, topography plays an essential role; for example, in the deep ocean, topography has been shown to increase the level of turbulence through breaking of internal waves [1][2][3][4][5]. Moreover, in strongly stratified and inhomogeneous turbulence with low buoyancy Reynolds numbers and in the absence of rotation and forcing, the flow structures consist of localized regions of turbulence, and gravity wave packets travel horizontally from these patches into surrounding quiescent regions, leading to mixing in the flow that is not related only to potential vorticity [6].
Shear instabilities lead to turbulent mixing as well [7][8][9]), but it is not clear what happens when shear layers are created self-consistently by the nonlinear dynamics of turbulent flows. The role of an externally imposed shear S has been demonstrated explicitly through the analysis of numerical simulations and observational data [10]; this was done in terms of the shear Richardson number Ri S = S 2 /N 2 , with N the Brunt-Väisälä frequency, Ri S being a central parameter in classical models of geophysical flows [11,12]. In the presence of rotation, an inverse cascade of energy to the large scales can take place, with the upscale transfer of energy possibly enhanced by shear [13]. The competition between rotation leading to the energy moving to large scales, at a rate UP and stratification transferring energy to the (small) dissipative scales at a rate down , has been analyzed recently [14,15]. It can be shown that the ratio of these two flux rates depends on the governing parameters of the problem, namely the Froude and Rossby numbers [15].
Similarly, oceanographic measurements in the vicinity of the Hawaiian ridge indicate the presence of patches of strong dissipation: with N = 10 −3 s −1 and u rms ≈ 0.1ms −1 the rms velocity, the measured kinetic energy dissipation is V ≈ 10 −6 W [16]. These data correspond to the dimensional evaluation of dissipation for a length scale of 1km, which is typical of the tidal excitation due to the bathymetry. With such strong dissipation, the flow recovers isotropy at small scale, at least locally, and indeed the enhanced dissipation is observed on vertical length scales on the order of a kilometer.
In this paper, we examine the link between kinetic energy dissipation, buoyancy flux and the magnitude of both the linear and nonlinear part of potential vorticity (see next section for definitions) through the analysis of joint Probability Distribution Functions (PDFs). We focus on three (fiducial) runs that are characteristic of the three regimes encountered in rotating stratified turbulence (RST), namely regime I dominated by quasilinear inertia-gravity waves (run rI ), regime II in which there is a competition between these waves and nonlinear eddies due to the advection terms of the hydrodynamics (run rII) and regime III dominated by the nonlinear eddies which can at times be faster than the waves (run rIII). Earlier global studies of these flows within a large parametric study have been presented elsewhere [17][18][19]. The main characteristics of the runs analyzed herein are given in Table 1 in the next section, together with the governing equations and dimensionless parameters. In Section 3 we recall global scaling properties of rotating stratified turbulence, and in Section 4, we analyze joint PDFs of several relevant variables for these three regimes. A discussion and concluding remarks are presented in Section 5. Table 1. Characteristics of the decay runs analyzed here: Froude and Rossby numbers, Fr, Ro and their ratio Ro/Fr = N/ f , as well as Reynolds and buoyancy Reynolds numbers, Re, R B (see text for definitions, and see [18][19][20] for more data on the runs). The last line refers to a run on a grid of 4096 3 points analyzed in [17], run r4k. Like run rII, it belongs to the intermediate (second) regime of parameters as defined in [18] for rotating stratified turbulence. The four runs have isotropic random initial conditions centered in the large scales, with zero initial potential energy.

Equations and Numerical Settings
The Boussinesq equations are written below for an incompressible fluid with velocity u = [u ⊥ , w] and temperature fluctuations θ about a mean that is linear in z, in the presence of rotation and gravity but in the absence of forcing: with ω = ∇ × u the vorticity, P the total pressure, N the Brunt-Väisälä frequency, f = 2Ω and Ω the imposed rotation frequency; ν and κ are the kinematic viscosity and scalar diffusivity;ê z is the unit vector in the vertical direction of the imposed rotation, Ω = Ωẑ, with gravity, g, in the opposite direction. These equations are written in units such that both u and θ have physical dimensions [L][T −1 ] per unit mass, whereas the buoyancy b = Nθ has, of course, the units of acceleration, [L][T −2 ]. The periodic domain has length L max = 2π. Wavenumbers thus vary from k min = 1 to k max = n p /3 with n p the number of points in all three directions, and using a standard 2/3 de-aliasing method. The boundary conditions are periodic in all three directions. Finally, we take a Prandtl number Pr = ν/κ equal to unity (see [21] for an evaluation of the turbulent Prandtl number in a plume, where it is indeed found to be of order unity). The dimensionless parameters governing the evolution of Rotating Stratified Turbulence (RST) as governed by Equations (1)-(3) are the Reynolds, Froude and Rossby numbers defined as usual as They all correspond to ratios of characteristic times, namely the gravity and inertial wave periods τ N = 1/N, τ f = 1/ f , the nonlinear turn-over time τ NL = L int /u rms and the dissipation time τ diss = L 2 int /ν. They are all defined and evaluated using large-scale characteristic quantities, with L int the integral scale and u rms the rms velocity. The Froude number can also be seen as the ratio of the typical turbulent velocity to the phase velocity of the gravity waves. Of course, in RST, we actually deal with inertia-gravity waves with a dispersion relation that combines rotation and stratification, and the perpendicular and parallel directions (with respect to the common vertical direction of gravity and the rotation axis). One also defines where, in terms of the strain-rate tensor, s ij = 1 2 (∂ i v j + ∂ j v i ), V = 2νs ij s ij is the pointwise kinetic energy dissipation rate, and β measures its efficiency against a dimensional evaluation, D . Finally, the point-wise gradient Richardson number, in the absence of imposed vertical shear, is written as The coupling between the velocity and the temperature fields occurs through the buoyancy flux B f defined as B N f represents a normalization of the buoyancy flux when compared to kinetic energy dissipation, and both can be defined as well in a point-wise manner.
The invariant of the Boussinesq equations, in the absence of dissipation (ν = 0, κ = 0), is the total energy E T = E V + E P , with E V = 1 2 |u| 2 the kinetic energy, and E P = 1 2 θ 2 the potential energy. Furthermore, the point-wise potential vorticity defined (neglecting the constant, N f ) as with ω = ∇ × u the vorticity, is also conserved in a Lagrangian sense (D t PV = 0 with D t = ∂ t + u · ∇). The code we have used for the simulations is pseudo-spectral and periodic in all three directions. It parallelizes efficiently using a combination of MPI, Open-MP as well as CUDA [22]. It also has a version with non-periodic boundary conditions in the vertical [23]. We take initial conditions that are typical for homogeneous isotropic turbulence (HIT), with zero temperature fluctuations (thus, zero initial potential energy), and velocity modes centered on the large scales, 2 ≤ |k 0 | ≤ 3. The total kinetic energy is normalized so that, at t = 0, u rms ≈ 1; isotropy is assumed and the phases are taken randomly. A second set of initial conditions is taken to be in quasi-geostrophic (QG) equilibrium which reflects a balance between pressure gradient, Coriolis and gravity force [18,19]. The QG regime is studied in detail, for example, in [24]; here we rather concentrate on the intermediate regime between the wave-dominated and the nonlinear eddy regimes [20].
The three runs analyzed in this paper are listed in Table 1 with their principal characteristics. They correspond, respectively, to run 5, 32 and 58 of Table 1 in [20] (see also Table I in [18]). All runs are performed for a few turnover times in the absence of external forcing, and grids of 1024 3 points are used. Results from each run presented herein draw mainly from a single snapshot of all fields on the computational domain taken at the peak of dissipation (see [20]). A brief analysis of an additional run is performed here in order to illustrate small-scale behavior of these flows. This stems from a large direct numerical simulation (DNS) on a grid of 4096 3 points in the absence of forcing (see [17] for a rather complete description of the properties of that large run).

Scaling Properties of Rotating Stratified Turbulent Flows
How do eddies and waves, at both large and small scales, interact in rotating stratified turbulence? We first recall a few recent results obtained using a large data base of more than 70 runs mainly on grids of 1024 3 points, together with a few other runs at lower resolution, and one data point for a grid of 4096 3 points (run r4k). This set of runs has already been analyzed in [17][18][19][20], and we now summarize some of the major results obtained so far. A table with the main characteristics of all these runs can be found in [18] (see also [17] for an earlier analysis focused on the temporal evolution and on the variation of characteristic time scales of kinetic and potential energy decay with governing parameters). The domain of controlling parameters is large and encompasses values (nearly) adequate for the atmosphere and oceans, except for lower Reynolds numbers as usual for DNS studies.
Turbulence is what couples together the nonlinear eddies and the inertia-gravity waves, leading to a remarkable balance between strong vertical gradients (giving rise to the observed layers in the atmosphere) and strong horizontal eddies or micro-vortices. A striking property of RST is the fact that the effective flow dissipation β (see Equation (5)) scales linearly with the controlling parameter, i.e., the Froude number, measuring how fast gravity waves are, when compared to nonlinear eddies. This is amply discussed in [18]. The behavior of β with Fr delineates three regimes: in regime I, the waves dominate and the dissipation that occurs depends on the Reynolds number. By contrast, in regime III, eddies can be faster than waves, and the flow behaves in ways similar to homogeneous isotropic turbulence, at least at small scales. Note however that there are strong signatures of persistent anisotropy even in that regime, see e.g., [19]. Between regime I and regime III, for an intermediate set of parameters (roughly, 0.01 ≤ Fr ≤ 1), the eddies become progressively faster, sharper gradients form allowing for more energy dissipation and a simple phenomenological argument compatible with weak wave turbulence leads to a linear relationship between β and Fr in what defines regime II.
We have also shown that there is a narrow interval of parameters in which rotating stratified flows display strong non-Gaussian wings in the vertical velocity w and in the temperature fluctuations, as illustrated in Figure 1 by the kurtosis K W =< w 4 > / < w 2 > 2 of w, as a function of Froude number, Fr; with this definition, a normally distributed vertical velocity will have K W = 3. The binning is done here in Rossby number, and the colors/symbols are given in the caption. Smaller symbols are for the 1024 3 runs, and larger symbols are for runs with larger viscosity and lower resolution. Finally, the stars indicate quasi-geostrophic initial conditions, in that case all with N/ f = 5; they are runs (Q12 through Q15, of Table 2 in [18])), with 0.067 ≤ Fr ≤ 0.111. Thus, there is intermittency of the vertical velocity for rotating stratified flows, as is already the case for non-rotating stratified turbulence, found in the presence of forcing [25][26][27].  [3,6), red squares for [6,10) and inverted magenta triangles for [10, ∞); stars are for runs with QG initial conditions.
The peak observed in Figure 1 occurs for Fr ≈ 0.07, a value comparable to that found in [26] in the purely stratified forced case. In the presence of rotation, the effect is particularly notable for the runs (indicated by stars), with initial conditions being in quasigeostrophic balance. It may be related to the fact that, for these flows, there is no vertical velocity initially; thus, w must be produced by the waves and the nonlinear dynamics occur in abrupt events leading to strong localized dissipation [27]. Note that this intermittency can be modeled in a simple way in the stratified case [25], a result that persists in the presence of rotation as we find here. By contrast, recall that, for HIT, the three components of the velocity are Gaussian [27,28]. Such intense and concentrated dissipative events, much stronger than their surroundings, are observed in the ocean, as in the Hawaiian Ridge or the Puerto-Rico Trench [16,29]. Note that extreme events in the distribution of aerosols in the atmosphere have been observed as well (see for example [30]).
There is of course intermittency also present in the small scales, as measured for example by large skewness and kurtosis of velocity and temperature gradients, with, again, a peak for the same intermediate regime of parameters. This is illustrated in Figure 2  In order to take a closer look at these PDFs, we give them in Figure 3 specifically for the three fiducial runs of Table 1 in terms of the ω 2 (horizontal) component of the vorticity; both Gaussian profiles (red circles; refer to the caption in Figure 3) and exponential fits of the form P(ω 2 ) = B exp (−Cω S 2 ) are indicated, with blue triangles giving the domain in which the fit is performed. This preliminary analysis indicates a (somewhat) lesser range for ω 2 in regime I, which is mostly dominated by waves. Such stretched-exponential behavior is suggested in many studies as being due to localized turbulence structures, as for shocks occurring in the Burgers equation. We note that the intermediate regime scales differently from the other two, but none of them scales as for the shocks in onedimensional pressure-less Burgers turbulence, or for models of drift-wave turbulence, for both of which S = 3/2 for quadratic nonlinearities and first-order moments [31,32]. Finally, in all cases, we postulate that these exponential wings in the PDFs are directly linked to the presence of strong, vertically sheared structures. Our study differs in several ways from preceding analysie of instantons in various turbulent flows [31][32][33][34][35][36][37]. Namely, it includes stratification and rotation, it includes only decaying flows, and the turbulent domain that appears at scales smaller than the Ozmidov scale (or, in other words, for high buoyancy Reynolds numbers) is not sufficiently resolved (see also Section 5). However, this preliminary analysis does indicate the presence of stretched exponentials, the more detailed analysis of which is left for future work.  ; the fit domain is given by blue triangles, and the red circles correspond to Gaussian profiles computed from the distribution mean and standard deviation.

Global Properties
We now concentrate on properties of potential vorticity, as it relates to other flow diagnostics. PV is an ideal conserved quantity that can be seen, in the framework of Noether's theorem, as stemming from the symmetry due to the invariance of the Boussinesq equations under the relabeling of Lagrangian particles [38,39]. PV can be decomposed into its linear (L) and nonlinear (NL) parts: Note that a time-scale can be associated with PV, namely τ PV = [PV] −1/2 . We then define the normalized versions of Π L and Π NL as Π 1 and Π 2 : (note that these notations differ from what is used in [40,41]); ω is sometimes called the relative vorticity, fê z the planetary vorticity and ω + fê z the absolute vorticity (see e.g., [42]). PV can actually be measured in the stratosphere [43,44]. It is also shown in [45] that PV conservation plays an important role in the dynamical evolution of stratified turbulence, as for example in tropical cyclones. Another remark is that the point-wise conservation of PV is not the only strong constraint exerted by the nonlinear dynamics on turbulent flows. There is in fact detailed conservation of total (kinetic + potential) energy for each allowable triadic interaction between (three) Fourier modes. Correspondingly, one can establish a flux conservation in configuration space for each distance |r| in the inertial range, assuming isotropy and homogeneity. This is expressed in terms of exact laws (under several hypotheses) for each point separation |r| entering the argument of velocity correlation functions, such as in the so-called 4/5th law of Kolmogorov [46] for HIT.
Moreover, it was shown using DNS with Taylor-Green initial conditions-corresponding to two counter-rotating vortices-that even though PV is an invariant in stratified flows, it does not undergo a classical cascade ( [47]). In particular, its dissipation is not concentrated around the classical dissipative scale (that is, the Kolmogorov scale η = [ V /ν 3 ] −1/4 ) and beyond, but in fact pervades many scales.
Stratified flows are known to support shear layers that go unstable, for example forming Kelvin-Helmoltz (KH) rolls, as observed in the atmosphere, the ocean or in the magnetopause [48][49][50]. We show in Figure 4a thin slice of potential vorticity PV, zooming on a Kelvin-Helmoltz event for the run r4k (see Table 1), at a resolution of 4096 3 points [17]. The zoom encompasses roughly a tenth of the compute box size on the side, and the KH event is barely discernible in terms of PV in the middle top right of the picture (see [17], Figure 10, for the visualization of the vorticity, temperature and gradient Richardson number in the same local view). This is due to the fact that PV emphasizes the small scales, more so than vorticity itself as it combines as well fronts developing in the scalar field. Indeed, PV can reach high values and is typical of small-scale dynamics, involving two derivatives of the basic fields, through the vorticity and the temperature gradient. Also, there are sharp contrasts of signs (purple and red for positive and negative respectively) with an imbrication of structures in this small-scale flow, potentially leading to different types of instabilities.
The relative importance of the different terms in the expression of PV has been measured for the Boussinesq equations [51][52][53]; as expected, the nonlinear part of PV is found by these authors to be negligible when waves are strong enough. On the other hand, when Π NL is large compared to the linear part, one can show that potential vorticity may become important as shear instabilities take place (see, e.g., [53]). PV is conserved point-wise, but its different elements (separating horizontal and vertical components) can undergo huge variations when following a single particle trajectory [53,54].
Using two-point closures of turbulence, which allow for numerical experiments at a high Reynolds number, it was shown in [55] that β (defined in Equation (5)) tends to constant values at high Taylor Reynolds numbers [56], values that may differ for different flows [57], by up to a factor of 5 [58], as in the presence or absence of shear (see also [59,60]). In particular, a relationship is found in [59] between the numerical value of β and the presence or not of stagnation points in the flow. This immediately implies that the efficiency of small-scale dissipation is related to the structure of the large-scale flow that determines its overall topology. In this light, the advantage of the present data base, which is consistent throughout as to the large-scale initial conditions, leaving a clear imprint of the influence of stratification on the dissipation, is also a drawback since the actual values of the dissipation efficiency depend on the overall geometry of the flow set by the large scales.
For example, we already know that, in the presence of a strong shear, dissipation is much more efficient, in particular through the formation of localized fronts and filaments [61]. In fact, one commanding question regarding turbulent flows is whether they are best described through a state of maximal energy dissipation, or of maximum entropy production (see, e.g., [62] and references therein). Rotating stratified flows in general develop sharp gradients with localized dissipative events. For N = 0, vertical Taylor columns form, but even weak stratification can alter them substantially. Such flows were studied experimentally in [63] using salt for 0 ≤ N/ f ≤ 0.24 with a Rossby number of ≈ 9 × 10 −3 . The effect of stratification was noticeable with a substantial shortening of the column: already for N/ f = 0.07 (Fr ≈ 0.13), it is reduced by a factor of three. In RST in general, structures are slanted, with an angle depending on the relative value of the stratification and rotation periods, as analyzed in [64] (see also Figure 4). For example, for flows with N/ f = L ⊥ /L (unit Burger number), the layers that form make an angle δ with the horizontal, with tan δ = L /L ⊥ . Moreover, the layers broaden when rotation increases [65], switching progressively from a vertical buoyancy scale u rms /N to a vertical QG scale [ f /N]L int . In addition, it can be shown analytically that frontogenesis disappears below a small critical Fr, for high rotation or stratification [66], when, on the other hand, strong, anomalous, kinetic and potential energy dissipation takes place in such flows through the formation of intermittent events at small scale.
Version January 14, 2021 submitted to Atmosphere 9 of 18

Joint Probability Distribution Functions
In order to examine the relationship between dissipation and mixing, we show now joint PDFs of various fields. We define as usual the joint distribution function of two random variables (u, v) as the probability of u having value u 0 and v taking simultaneously the value v 0 , that is P J (u, v) = P (u = u 0 , v = v 0 ). Such joint PDFs are plotted in the next few figures for a variety of variables, with a color bar to the right of each plot indicating the amplitude of the correlation between the two fields, defined as the log of the number of samples in the two-dimensional bin divided by the maximum number in any bin.
In Figure 5, we display the buoyancy flux density wθ against the point-wise Richardson number Ri g for the three fiducial runs chosen in this paper; note the change of range for Ri g between rI and rII, rIII. For all runs, the highest probability occurs for a zero buoyancy flux, but as the Froude number increases, by roughly a factor 10 between rI and rII, and rII and rIII, the centroid of the joint distribution moves to Ri g ≈ 0, with the range of values attained by Ri g narrowing: whereas in regime I, the flow is mostly stable, and in regime III the flow is close to instability (be it convective or KH), almost everywhere. One can envision the development of turbulence in these systems as being a spatial propagation of quasi-linear instabilities progressively over the entire flow, leading to secondary instabilities and so forth. In this manner, one could say that the localized turbulent hot-spots at low Fr, R IB are increasing in size as Fr, R IB increases and that regime III first obtains when the direct surroundings of these different hot-spots become contiguous, and the flow becomes dominated by nonlinear eddies with Kolmogorov spectra and energy dissipation efficiency of the order of that found in HIT.  Table 1). Note the extended range of values for Ri g in the left panel. AP: Normalized or not, see Equation (7). Figure 6 shows how the efficiency of kinetic energy dissipation β (defined in Equation (5)) varies with point-wise buoyancy flux, in two ways. At the top, β is correlated with the normalized buoyancy flux B N f (see Equation (7)), whereas at the bottom β is given using logarithms in order to display its wide range of values, and it is plotted now against the raw buoyancy flux (divided by N).
Local values of β can be much larger than its average, which should remain close to 1 (from below), by construction. However, note that the dimensional evaluation of dissipation, namely D , is constructed on the rms value of the overall velocity (and on the overall integral scale). However, in fact, dissipation does occur in highly localized patches of a typical dimension that is smaller than L int , and with higher velocities since it is in these patches that the vertical velocity becomes significant (see for example [26] for the purely stratified case). Thus the local dimensional dissipation may be larger, leading to higher local βs. The way higher probabilities populate the space from rI to rII is also significant: it indicates that in the wave regime, most data points in the flow have barely any dissipation, whereas in the regime dominated by turbulent eddies, the highest correlation is for significant efficiency for all values of B N f , but with more at small B f , and with a significant growth of dissipation efficiency for low buoyancy fluxes. In that regard, run rII is clearly transitional. Finally, we note in the bottom plot that the centroid of the distribution moves upward, to higher local (and global) values of the dissipation as the turbulence strengthens from rI to rIII. When we plot the joint PDFs of the normalized linear and nonlinear parts of the potential vorticity, Π 1 , Π 2 for the same three runs (see Figure 7, top and bottom respectively) against normalized buoyancy flux, we observe again a drastic change between the wavedominated and eddy-dominated regimes. Note also that the two plots (for Π 1,2 ) give complementary information. For Π 2 , at low Froude number (for rI, left), the highest joint probability corresponds, as expected, to very low kinetic energy dissipation (normalized buoyancy flux close to unity) and zero nonlinear PV. Conversely, for Π 1 , most of the points with high probability lie for B N f ≈ 1, Π 1 ≈ ±1 (in normalized terms). By contrast, for run III for which β has reached a value close to the efficiency of HIT [18], the highest probability for Π 2 is for PV dominated by its nonlinear part and for B N f ≈ 0, although in fact the buoyancy flux takes a large range of values. Run rII, in the intermediate regime, is a combination of these two behaviors, although it does display its highest probability for V still quite small, with B N f close to unity. The wavy appearance of the data for rII for PV at high B N f is likely linked to a symmetry between the ± parts of Π 1,2 , with extrema at ±0.5 corresponding to the linear and nonlinear parts of PV balancing each other. This occurs for B N f ≈ 1, that is for low kinetic energy dissipation. For rI, the highest joint probability is for V ≈ 0, which of course makes sense; for rIII, one finds high joint probability for Π 2 ≈ ±1, which again is to be expected since Π 1 ≈ 0. rII shows characteristics of both the rI and rIII regimes. We can say that whatever the value of the buoyancy flux B f , the highest joint probability is found for Π 1 ≈ 0 (and consequently Π 2 ≈ 1) for the run in the strongly turbulent regime, and the converse is true in regime I in which waves are fast. In Figure 8, we present the joint PDFs of the linear (top) and nonlinear (bottom) parts of potential vorticity, each normalized by the square of the turn-over time τ NL built on large-scale fields (rms velocity and integral scale). The highest joint probability, in all cases, takes place for β ≈ 0 and PV ≈ 0. As the Froude number progressively grows from run rI (left) to run rIII (right), there are higher values of the efficiency of kinetic energy dissipation that are populated, i.e., more local patches of highly dissipative flow, with the highest probability for rIII as we get closer to the fully turbulent regime. This can be verified by remarking that, taking some fixed value of joint probability, like 10 −6 , the level of Π NL increases with Froude number, going from rI to rIII, at the same time that the level of Π L decreases.
Finally, note in Figure 8 that, while the range for β remains the same for all six plots, the scale for the nonlinear PV is a factor of >5 times that of the linear PV. Again, these high values for the nonlinear component correspond to high vorticity and high temperature gradients, with no preferred direction of orientation. It is also noteworthy that nonlinear PV does not reach the whole range of values for the run dominated by the waves (rI, left). In fact, in general, the functional support for the linear PV varies in the opposite sense from the nonlinear PV as we traverse regimes, with that for linear PV decreasing from rI to rIII, while that for nonlinear PV increases. We also note that PV takes strong values for both ± signs, a probable signature of different instabilities such as the symmetric and baroclinic instabilities that can occur within geophysical flows.
Plots are for the three fiducial runs with, as before, rI at left, rII in the middle and rIII at right (see Table 1). Note the different vertical scales between the linear and nonlinear PV.

Discussion And Conclusions
The results presented in this paper show in particular that the nonlinear part of potential vorticity, Π NL = ω · ∇θ, which is weak initially for initial conditions centered on large scales, comes to dominate small-scale dynamics and is strongly correlated with highly efficient local dissipation involving local instabilities. As the Froude number increases, the turbulence becomes stronger, and small-scale gradients form. The point-wise conservation of potential vorticity implies that, when strong gradients form reinforcing Π NL , the linear part of PV has to decrease in the same proportion, explaining the large ratio between the two as observed in Figures 7 and 8.
One of the many remaining issues is whether and how Π NL can be weakened through the orthogonalization of the vorticity and the temperature gradient. Note that in fact there is another invariant of a helical type, I H = A · B with B = ∇ × A, where I H involves, this time, potential vorticity gradients [67]. The dynamics of I H were considered in the presence of an imposed vertical shear in [68], where it was shown that it does become weaker through alignment of A and B. The quasi-suppression, in strong local structures, of nonlinearities in turbulent flows is in fact a common feature where, for example, in HIT small-scale filaments of vorticity form in which ω is close to being parallel to u; thus, the Lamb vector u × ω that represents the nonlinearity of the Navier-Stokes equations to within a pressure term ∇|u| 2 /2 is negligible, making these filaments long-lived. Other examples of systematic weakening of nonlinear interactions can be found in MHD fluids and plasmas. The issue of the relative orientation of vorticity and temperature gradients is left for future work.
In the presence of a strong imposed vertical shear S and of strong rotation, as controlled by the Richardson and Rossby numbers, it was found in [69] that the baroclinic instability can develop efficiently. Shear and rotation were essential to determine the unstable state that leads to nonlinear mode coupling, and to the formation of strong small-scale gradients, in turn giving rise to localized enhanced dissipation as observed in the atmosphere and in the ocean. PDFs of the linear and nonlinear parts of PV (with a normalization different from what is done here, using τ N instead of τ NL ) are also given in that same paper, where it is observed that they are highly non-Gaussian. Indeed, strong (order unity) negative skewness develops in these flows, significant of strong localized turbulence [69]. The localization of turbulence, both in space and in time, has been studied for many years as it leads, for example, to anomalous scaling of structure functions. There is renewed interest in these issues within the context of instanton theory for both fluids and plasmas [32][33][34]37], and its link to the so-called PQR analysis of the invariants of the velocity gradient matrix [36], following the Vieillefosse model which, in the context of strongly stratified flows, indicates the essential role of convective instabilities [26,54,70].
Furthermore, it was shown in [71] that for HIT, but in the presence of a large-scale imposed anisotropy, the ratio of the dissipation rates estimated with either vertical (w 3 /L z ) or horizontal (u 3 ⊥ /L ⊥ ) variables are in fact equal, showing that there is a disappearance of the imposed anisotropy when one reaches the small scales, provided the Reynolds number is large enough. The link between anisotropy and the structure of potential vorticity will be the topic of a future paper.
Many other points require further discussion. For example, homogeneous turbulence in the presence of shear can be studied with a combination of rapid distortion theory and DNS. In that case, it can be shown that rotation has a huge impact on the dynamics of a passive scalar such as a chemical tracer, on the velocity-scalar cross-correlation (related to the buoyancy flux wθ ), as well as on the ratio of the turn-over time to the characteristic time associated with the dynamics of the scalar [72]. In our case, the Froude number is quite small for most of the runs, making the scalar interacting in a significant way with the velocity field. On the other hand, when the Froude number becomes close to unity (from below), atmospheric and oceanic constraints (that is, a ratio N/ f larger than 5 for the atmosphere and more like 100 for the oceans) lead to a Rossby number that is is too large for rotation to be significant, except in the largest scales, and the scalar becomes passive.
More to the point, in purely stratified flows, for scales smaller than the Ozmidov scale [73], temperature fluctuations become passive, and isotropy and classical (Kolmogorov) scaling in terms of isotropic wavenumber for the kinetic and potential energy spectra are conjectured to recover for a large enough Reynolds number. This was shown clearly in the purely rotating case beyond the so-called Zeman scale (where f replaces N in the above expression), analyzing a DNS run on a grid of 3072 3 points with Re = 27,000, Ro = 0.07) [74]. However, in the presence of both rotation and stratification, the recovery of anisotropy is slow, in particular for the vorticity, a fact which could be attributed to the invariance of PV [19]. There are, of course, several potential limitations to our work. One concerns the Prandtl number Pr taken equal to one. However, one can argue that the turbulent Prandtl number is of order unity, as demonstrated analytically using the Renormalization Group formalism for homogeneous isotropic turbulence in [75]. For stratified flows, it was also shown analytically for Ri < 0.2, using the Quasi-Normal Scale Elimination model [76], and it was demonstrated numerically for R B ≥ 50 in [77]. Furthermore, one can find a compilation of various atmospheric data sources in [78], where Pr is shown to have a limit of one for a Richardson number Ri ≤ 0.1, that is when the turbulence becomes energetic. Such a result, in terms of magnetic Prandtl number close to unity, is also present when coupling the velocity to a magnetic field for anisotropic flows [79,80].
We also want to remark that the runs that are the least well resolved numerically (comparing the cut-off frequency and the Kolmogorov dissipation wavenumber) are not those with the highest R B , as could have been expected a priori, since as R B grows, one enters the fully turbulent regime with more efficient and vigorous dissipation at small scales. However, it is shown in [25] for the purely stratified case (see Figure 1 in this paper for RST) that at intermediate values of R B in the so-called saturation regime in which there is a balance between nonlinear advection and buoyancy flux, intermittency, as measured by Probability Density Functions of the fields themselves, have non-Gaussian wings contrary to the HIT case (see also [26,27,81]), thus requiring more numerical resolution to accurately be able to follow the small scales of the intermediate regime.
Analyzing forced flows will give access to longer temporal averages to compute accurately the vertical buoyancy and momentum fluxes. However, one will then have to cope with issues of non-stationarity at small Rossby number. Indeed, it is known that below a critical Rossby number Ro c ≈ 0.1, the energy cascade goes predominantly to large scales [82,83] with a linear temporal growth of the kinetic energy, unless a large-scale friction term is present.
Finally, novel numerical methods of decimation, either progressive [84] or fractal [85], as well as new data analysis algorithms in the spirit of artificial intelligence algorithms [86,87], may be found to be of some use in this context. They could help tackle these important problems at the core of the understanding of weather and climate systems.
Author Contributions: All authors contributed equally. All authors have read and agreed to the published version of the manuscript.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data in whole or in part may be requested from authors.