Baropycnal Work: A Mechanism for Energy Transfer Across Scales

The role of baroclinicity, which arises from the misalignment of pressure and density gradients, is well-known in the vorticity equation, yet its role in the kinetic energy budget has never been obvious. Here, we show that baroclinicity appears naturally in the kinetic energy budget after carrying out the appropriate scale decomposition. Strain generation by pressure and density gradients, both barotropic and baroclinic, also results from our analysis. These two processes underlie the recently identified mechanism of"baropycnal work,"which can transfer energy across scales in variable density flows. As such, baropycnal work is markedly distinct from pressure-dilatation into which the former is implicitly lumped in Large Eddy Simulations. We provide numerical evidence from 1,024^3 direct numerical simulations of compressible turbulence. The data shows excellent pointwise agreement between baropycnal work and the nonlinear model we derive, supporting our interpretation of how it operates.


I. INTRODUCTION
Energy transfer across length scales is one of the defining characteristics of turbulent flows, the subject of which fits well under the "Multiscale Turbulent Transport" theme of this special issue in Fluids. In constant density turbulence, the only pathway for transferring energy across scales is deformation work [1], which we represent below by Π. This is often referred to as the turbulence production term in the turbulent kinetic energy (TKE) budget within the Reynolds averaging decomposition [2] or the spectral flux within the Fourier decomposition [3][4][5]. Deformation work gives rise to the cascade in incompressible turbulence, which is largely believed to operate by vortex stretching in 3-dimensions [1,[6][7][8][9], an idea which may be traced back to G.I Taylor [10,11].
Recent studies [12][13][14][15][16][17] have shown that in the presence of density variations, such as in compressible turbulence, there exists another pathway across scales called "baropycnal work," represented by Λ below. In the traditional formulation of the compressible Large Eddy Simulation (LES) equations [18][19][20][21], which are essentially a coarse-graining decomposition of scales, Λ is almost always implicitly lumped with P ∇·u, where P is pressure and u is velocity, and treated as a large scale (resolved) pressure-dilatation which does not require modeling. Ref. [14] argued that baropycnal work, Λ, is more similar in nature to deformation work, Π, in that it involves large-scales interacting with small-scales thereby allowing it to transfer energy across scales. As such, it is fundamentally distinct from pressure dilatation which involves only large-scales and cannot transfer energy directly across scales.
In this work, we shall investigate the mechanisms by which baropycnal work transfers energy across scales. The main result is embodied in eq. (16) below, which shows that Λ transfers energy by two processes: I) Barotropic and baroclinic generation of strain, S, from gradients of pressure and density, ρ: where the dyadic product ∇ρ (∇P ) T is a tensor. Length scale is that at which density and pressure gradients are evaluated as we make clear in eqs. (16)- (18) below. To our knowledge, these results are the first to show how baroclinicity, ∇ρ×∇P , appears in the kinetic energy budget. Baroclinicity is often analyzed within the vorticity budget but its role in the kinetic energy budget has never been obvious. We shall show here that baroclinicity appears naturally in the kinetic energy budget after performing the appropriate scale decomposition.
Strain generation by pressure and density gradients (both barotropic and baroclinic) also results from our analysis, highlighting its potential significance which is often overlooked in the literature. The processes of strain and vortex generation show baropycnal work Λ to be markedly distinct from the process of pressure dilatation, further supporting the argument against lumping the two terms.
There is a diverse array of applications in this research subject. Density variability can arise in high Mach number flows, but it is also pertinent in the limit of low Mach numbers along contact discontinuities such as in multi-species or multi-phase flows [22]. Variable density (VD) flows are relevant in a wide range of systems, such as in molecular clouds in the interstellar medium [23][24][25]), in inertial confinement fusion [26][27][28], in high-speed flight and combustion [29,30], and in air-sea interaction in geophysical flows [31][32][33][34].
The outline of the paper is as follows. In section II, we shall use coarse-graining to decompose scales and identify baropycnal work, Λ. In section III, we use scale locality to approximate Λ with a nonlinear model, which in turn shows how Λ is due to a combination of strain generation and baroclinic vorticity generation. In section IV, we describe our direct numerical simulations (DNS) used in section V to present evidence that indeed Λ and its nonlinear model exhibit excellent agreement. This justifies our analysis of Λ via its nonlinear model, similar to what was done by [6] in their analysis of Π. We conclude with section VI.

II. MULTI-SCALE DYNAMICS
To analyze the dynamics of different scales in a compressible flow, we use the coarsegraining approach, which has proven to be a natural and versatile framework to understand and model scale interactions (e.g. [35,36]). The approach is standard in partial differential equations and distribution theory (e.g., see Refs. [37,38]). It became common in Large Eddy Simulation (LES) modeling of turbulence thanks to the original work of Leonard [39] and the later work of Germano [40]. Eyink [36,[41][42][43] subsequently developed the formalism mathematically to analyze the fundamental physics of scale coupling in turbulence.
For any field a(x), a coarse-grained or (low-pass) filtered field, which contains modes at scales > , is defined in n-dimensions as where G(r) is a normalized convolution kernel and G (r) = −n G(r/ ) is a dilated version of the kernel having its main support over a region of diameter . The scale decomposition in (1) is essentially a partitioning of scales in the system into large ( ), captured by a , and small ( ), captured by the residual a = a − a . In the remainder of this paper, we shall omit subscript from variables if there is no risk for confusion.

A. Variable Density Flows
In incompressible turbulence, our understanding of the scale dynamics of kinetic energy centers on analyzing |u | 2 /2. In variable density turbulence, scale decomposition is not as straightforward due to the density field ρ(x). Several definitions of "large-scale" kinetic energy have been used in the literature, corresponding to different scale-decompositions as discussed in [14]. These include ρ |u | 2 /2, which has been used in several studies (e.g. [67][68][69][70]), and |( √ ρu) | 2 /2, which has also been used extensively in compressible turbulence studies (e.g. [16,[71][72][73] [18,76], |ρu | 2 /2ρ , satisfies the inviscid criterion, which allows for properly disentangling the dynamical ranges of scales. We will use the common notation which yields ρ | u | 2 /2 for kinetic energy at scales larger than . The budget for the large-scale KE can be easily derived [14] from the momentum equation (20) below: where J (x) is space transport of large-scale kinetic energy, −P ∇·u is large-scale pressure dilatation, D (x) is viscous dissipation acting on scales > , and inj (x) is the energy injected due to external stirring. These terms are defined in eqs. (16)-(18) of Ref. [14]. The Π (x) and Λ (x) terms account for the transfer of energy across scale , and are defined as where is a 2 nd -order generalized central moment of fields f (x), g(x) (see [40]).
The first flux term, Π , is similar to its incompressible counterpart and is often called deformation work. The second flux term, Λ , was identified in [12,14] and called "baropycnal work." It is inherently due to the presence of a variable density and vanishes in the incompressible limit. Recent work by Eyink and Drivas [17,77] identified a third possible pathway for energy transfer, which they called "pressure-dilatation defect" and arises when the joint limits of κ, µ → 0 and → 0 do not commute. Eyink and Drivas [17] showed that the pressure-dilatation defect mechanism transfers energy downscale in 1D normal shocks.
In this paper, we shall focus on understanding the mechanisms by which baropycnal work transfers energy across scales.

III. THE MECHANISM OF BAROPYCNAL WORK
Our investigation of the mechanism behind baropycnal work is inspired by the work of Borue & Orszag [6], where they used the nonlinear model of the energy flux Π to show that it operates, on average, by vortex stretching (see also eq. (32) in [7]).
Our derivation of the nonlinear model of Λ will follow that in [78], which is somewhat different from the standard derivation of nonlinear models [6,55,79]. We utilize the property of scale-locality [36] (specifically, ultraviolet locality) of the subscale mass flux which was proved to hold in variable density flows by [12] under weak assumptions. Specifically, for our present purposes, we require that the spectra of density and velocity decay faster than k −1 in wavenumber. In other words, density and velocity should have finite second-order moments, ρ 2 < ∞ and |u| 2 < ∞, in the limit of infinite Reynolds number. Here, . . .
is a space average. Ultraviolet scale locality implies that contributions to the subscale mass flux τ (ρ, u) at scale from smaller scales δ are negligible [12]: By assuming the validity of eq.(7) in the limit δ → , we can justify the approximation which neglects any contribution from scales < to the subscale mass flux. Using the usual definition of an increment: the subscale mass flux term can be rewritten exactly in terms of δρ and δu [12,36]: Equation (10) is exact, where is a local average around x over all separations r weighted by the kernel G . A spatially localized kernel effectively limits the average to separations |r| /2. Since a filtered field f (x) is smooth, we can Taylor expand its increments around x where we neglect higher order terms.
Substituting the first term in the Taylor expansion of each of δρ and δu into eq. (10) gives This is the nonlinear model of the subscale mass flux τ (ρ, u i ). In deriving the second line, we used the symmtery of the kernel such that r k = 0. In the final expression, C 2 = d 3 r G(r) |r| 2 depends solely on the shape of kernel G and, in particular, is independent of scale .
We now have a nonlinear model of Λ that is only a function of filtered fields, which are resolved in LES simulations. We can use this model to gain insight into the mechanism by which Λ transfers energy across scales. The velocity gradient tensor ∂ k u j can be decomposed into symmetric and antisymmetric parts with where ω = ∇×u is vorticity and ijk is the Levi-Civita symbol. Therefore, Λ at any scale can be approximated by a nonlinear model, Λ m , everywhere in space: where is the strain generation process of baropycnal work and is its baroclinic vorticity generation process. Equation (16) is the main result of this paper.
In the following sections, we will provide numerical support showing excellent pointwise agreement between Λ and its nonlinear model Λ m .
To illustrate how strain generation by baropycnal work takes place, consider an unstably stratified flow configuration in which ∇P and ∇ρ in Λ SR are anti-aligned (∇ρ·∇P < 0) as illustrated in Fig. 1 of [14], and both are parallel to a contracting eigenvector of S (associated with a negative eigenvalue of S). Remember that the strain S in our nonlinear model arises from τ (ρ, u) which represents scales smaller than . In such a configuration, the contraction (and therefore strain) is enhanced leading to the generation of kinetic energy in the form of straining motion at scales smaller than (Λ > 0 in eq. (3)). The ultimate source of kinetic energy being transferred by Λ to motions at scales < is potential energy due to the large-scale pressure gradient, ∇P .
The baroclinic component, Λ BC , demonstrates how baroclinicity, ∇ρ×∇P , plays a role in the energetics across scales. The importance of baroclinicity is well known [80,81] but it has always been analyzed within the vorticity budget. Its contribution to the energy budget has never been clear. Baroclinicity in the kinetic energy budget arises naturally from our scale decomposition and the identification of Λ as a scale-transfer mechanism. The need for a scale decomposition in order for Λ and, as a result, baroclinic energy transfer, to appear in the kinetic energy budget should not be surprising. This is similar to the scale transfer term Π, which does not appear in the budget without disentangling scales due to energy conservation. In the same vein, the appearance of baroclinicity in the vorticity equation can be interpreted as being a consequence of an effective scale decomposition performed by the curl operator ∇×, which a high-pass filter.
As mentioned in the introduction, in the compressible LES literature, Λ is almost always lumped with pressure-dilatation, P ∇·u in the form of P ∇· u [18][19][20][21], thereby completely missing the physical processes inherent in baropycnal work. Our analysis here supports the argument in [12,14] to separate Λ from pressure-dilatation. In those studies, it was reasoned that Λ and P ∇·u are fundamentally different; the former involves interactions between the large scale pressure gradient with subscale fluctuations, allowing the transfer energy across scales, whereas the latter is solely due to large-scale fields and cannot participate in the transfer of energy across scales.

IV. SIMULATIONS
To provide empirical support to our nonlinear model of Λ, we carry out a suite of DNS of forced compressible turbulence in a periodic box of size 2π on which we perform a priori tests of our derived model against simulation data. The DNS solve the fully compressible Navier Stokes equations: Here, u is velocity, ρ is density, E = |u| 2 /2 + e is total energy per unit mass, where e is specific internal energy, P is thermodynamic pressure, µ is dynamic viscosity, q = −κ∇T is the heat flux with a thermal conductivity κ and temperature T . Both dynamic viscosity and thermal conductivity are spatially variable, where µ(x) = µ 0 (T (x)/T 0 ) 0.76 . Thermal conductivity is set to satisfy a Prandtl number P r = c p µ/κ = 0.7, where c p = R γ/(γ − 1) is the specific heat with specific gas constant R and γ = 5/3. We use the ideal gas equation of state, P = ρRT . We stir the flow using an external acceleration field F i , and RL represents radiation losses from internal energy. S ij = (∂ j u i + ∂ i u j )/2 is the symmetric strain tensor and σ ij is the the deviatoric (traceless) viscous stress We solve the above equations using the pseudo-spectral method with 2/3rd dealiasing.
We advance in time using the 4 th -order Runge-Kutta scheme with a variable time step.
The acceleration F we use is similar to that in [82]. In Fourier space, the acceleration is defined as where the complex vector f is constructed from independent Ornstein-Uhlenbeck stochastic processes [83] and the projection operator P ζ ij (k) = ζδ ij + (1 − 2ζ) k i k j |k| 2 allows to control the ratio of solenoidal (∇ · F = 0) and dilatational (∇ × F = 0) components of the forcing using the parameter ζ. When ζ = 0, the forcing is purely dilatational and when ζ = 1, the forcing is purely solenoidal. The acceleration is constrained to low wavenumbers |k| < k F . For the internal energy loss term RL, we have tested two schemes: a spatially varying radiative loss, RL = ρu·F, and another that is independent of space, RL = ρu·F . The two schemes yield indistinguishable results. Including an internal energy loss term is similar to what was done in other studies [84,85] to allow total energy to remain stationary. After an initial transient, mean kinetic energy reaches a statistically stationary state as shown in The compressibility metrics in Table I show the relative importance of dilatational versus solenoidal velocity modes. We use the Helmholtz decomposition, u = u d + u s to obtain the dilatational (∇ × u d = 0) and solenoidal (∇ · u s = 0) components of the velocity field.
The dilatational kinetic energy is K d = ρu d i u d i /2 and the solenoidal kinetic energy is K s = ρu s i u s i /2 . Their ratio K d /K s yields a measure of compressibility at large scales. It is well-known [82,86] that the dilatational kinetic energy K d becomes significant when forced directly using F with a small ζ. This holds, even though the low-ζ runs have a lower Mach number than high-ζ runs. The ratio (∇·u) rms /(∇×u) rms yields a measure of compressibility at small scales.
The last two columns in table I summarize the effect of ζ on the relative importance of Π Low-ζ corresponds to high compressibility in the external forcing. The spatially averaged cascade terms Λ and Π are calculated using the sharp-spectral cutoff filter with k = 2π/ = 6. ∆x η is the ratio of grid size to the Kolmogorov length. and Λ. While the deformation work Π is significant in all cases, baropycnal work Λ, which arises only in variable density flows, is greatly affected by the type of forcing used. At the Reynolds numbers we simulate, we find that Λ becomes important only for low ζ when the dilatational modes are directly forced. Even at relatively high Mach numbers, Λ remains small for high ζ. We caution, however, that these observations might be Reynolds number dependent. Moreover, Λ has been shown to dominate in non-dilatational low Mach number variable density flows [87,88]. The stark qualitative difference shows the significance of dilatational forcing on the flow [82,89], at least in limited resolution simulations. It has been argued [84,90,91] that at sufficiently high Reynolds numbers, flows forced dilatationally will produce sufficient vortical motion to resemble those forced solenoidally. Since we are primarily interested in the Λ term here, unless stated otherwise, plots that follow will be from low-ζ simulations where Λ is significant. Figure 3 shows the velocity spectra in the case of highly compressive (low-ζ) forcing. At intermediate scales, the spectrum of u d seems to follow a power law close to k −2 , which is expected for the dilatational component [16,24,86]. It is well-known [e.g. 16,24] that obtaining a clear power-law scaling of the solenoidal velocity is challenging when forcing dilatationally, even at our 1,024 3 resolution. Figure 4 shows the cascade terms Π and Λ averaged over the domain as a function of the filter wavenumber k. As is the case in 3D isotropic incompressible turbulence, Π is positive for all wavenumbers, transferring kinetic energy from large to small scales. On the other hand, Λ is negative, effectively reducing the total amount of energy transferred across scales. This is consistent with previous studies which measured Λ in homogeneous isotropic compressible turbulence [16]. Across a shock, the pressure and density gradients have the same direction and are aligned with the contracting strain eigenvector, leading to negative baropycnal work [14,17], thereby reducing the intensity of the cascade. The situation is different in buoyancy driven (unstably stratified) flows, where pressure and density gradients are in opposite directions leading to positive baropycnal work [14]. Within the framework of Reynolds-averaged Navier-Stokes (RANS), it has been shown that for variable density flows, such as turbulence generated by the Rayleigh-Taylor instability, the RANS equivalent of Λ is the largest contributor to the kinetic energy cascade [87]. We shall present our results on Λ in buoyancy driven flows in forthcoming work [88]. for x ≥ 0 and H(x) = 0 for x < 0. The correlation coefficient R c is shown at two scales k = 2π/ .
Filter type Kernel

V. NUMERICAL RESULTS
To quantify the pointwise agreement between baropycnal work Λ and its nonlinear model Λ m in our DNS, we measure the correlation coefficient: We also analyze the joint probability density function (PDF) in Figs. 5-7, and visualize Λ and Λ m in x-space in Fig. 8.
In our study, we use the filters defined in table II. Both the box and Gaussian filters are positive in physical space, which is important to guarantee physical realizability of filtered quantities [45]. On the other hand, the sharp spectral filter is not sign definite in x-space, which limits its utility in analyzing scale process in physical space.
Our results indicate an excellent agreement between baropycnal work and its nonlinear model. Using either the Gaussian or Box filters, the correlation coefficients are very high, R c > 0.9, for all the length scales we analyzed. The sharp spectral filter, on the other hand, yields poor agreement. This is not surprising since the sharp spectral filter can yield negative filtered densities [14,74] and physically unrealizable subscale stresses [45] due to its non-positivity in x-space. showing an almost perfect pointwise correlation. We note that in our dilatationally forced flows, most of the contribution to Λ m is from its straining component, Λ SR , with a negligible contribution from Λ BC (see eq. (16)). This is due to the shocks which contribute mostly to Λ SR . In flows dominated by baroclinicity, such as in the Rayleigh-Taylor instability, a significant contribution to Λ comes from Λ BC , as will be shown in forthcoming work [88].

VI. SUMMARY
Past work [12,14,17] has identified baropycnal work, Λ, as a process capable of transferring kinetic energy across scales in addition to deformation work, Π. This paper aimed at elucidating the physical mechanism by which Λ operates.
Using scale-locality [12] and a multiscale gradient expansion [78], we derived a nonlinear model, Λ m , of baropycnal work. Using DNS, we showed excellent agreement between Λ and Λ m everywhere in space and at any time, giving further empirical justification for our analysis of Λ via its model Λ m .
We found that baropycnal work operates by the baroclinic generation of vorticity, and also by strain generation due to pressure and density gradients, both barotropic and baroclinic.
While the role of pressure and density gradients in generating vorticity is well recognized, their role in strain generation has been less emphasized in the literature.
As far as we know, this is the first direct demonstration of how baroclinicity enters the kinetic energy budget, which arises naturally from our scale decomposition and the identification of Λ as a scale-transfer mechanism. Baroclinicity is often analyzed within the vorticity budget but its role in the energetics has never been obvious. The need for a scale decomposition in order for Λ and, as a result, baroclinic energy transfer, to appear in the kinetic energy budget is similar to the scale transfer term Π, which only appears in the budget after decomposing scales due to energy conservation. In the same vein, the appearance of baroclinicity in the vorticity equation can be interpreted as being a consequence of an effective scale decomposition performed by the curl operator ∇×, which is a high-pass filter.
Our findings here support the argument in [12,14] to separate Λ from pressure-dilatation, P ∇·u in compressible LES, where the two terms are often lumped together in the form of P ∇· u.
In forthcoming work, we shall present further evidence of the excellent agreement between Λ and Λ m using low Mach number buoyancy driven flows with significant density variability [88]. their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, complete-ness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof.
The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.