A local diagnostic energy to study diabatic effects on a class of degenerate Hamiltonian systems with application to mixing in stratified flows
Alberto Scotti Pierre-Yves Passaggia








^1 Dept. of Marine Sciences, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599


Oceanography, Fluid dynamics





A. Scotti

ascotti@unc.edu



In Hamiltonian systems characterized by a degenerate Poisson algebra, we show how to construct a local energy-like quantity that can be used to study diabatic effects on the evolution of the Available Energy of the system, the latter concept formalizing the original idea of Margules'.  
We calculate the local diagnostic energy for geophysically relevant flows. For the particular case of stratified Boussinesq flows, we show that under moderately general conditions, in inertial frames where the initial distribution of potential vorticity is even around the origin, our framework recovers the Available Potential Energy introduced by Holliday and McIntyre <cit.>, and as such depends only on the mass distribution of the flow. 
In non-inertial frames, we show that the local diagnostic energy of flows which are, in an appropriate sense, characterized by a low-Rossby number Ro ground state, has to lowest order in Ro, a universal character. 









§ INTRODUCTION


A major  challenge in contemporary oceanography is to understand the role of small scale turbulence in regulating the oceanic Meridional Overturning Circulation (MOC), the process that over
millennial time scales exchanges surface with deep water <cit.>. Given the large storage capacity of the deep ocean for heat and greenhouse gases, 
understanding the drivers of the MOC is essential for climate prediction <cit.>. 
<cit.> used energetics arguments in an attempt to quantify the amount of energy required to sustain the MOC,
but significant gaps remain in quantifying the pathways that energy injected
 
A local diagnostic energy to study diabatic effects on a  class of degenerate Hamiltonian systems with application to mixing in stratified flows
    
A. Scotti^1 and P.-Y. Passaggia^1
    December 30, 2023
================================================================================================================================================

at large scales 
by winds and tides take to reach the small 
scales at which mixing occurs <cit.>.
Particularly vexing is the problem of estimating the energetic cost incurred when  turbulent processes irreversibly mix the stratifying agents 
<cit.>, that is how to estimate the energy input necessary to sustain a given rate of dissipation of the variance of the stratifying agents. 




For a fluid in the Boussinesq approximation, <cit.> used the concept of Available Potential Energy (APE) as a diagnostic tool to achieve such a relation. APE has a long history, starting with the pioneering work of <cit.> and <cit.>, later developed by <cit.>  into a working tool, the so-called Lorenz Energy Cycle, which is still used today <cit.>. Essentially, Margules' idea was to calculate a minimum-energy state, compatible with certain constraints, and use the energy of such state to "gauge" the amount of potential energy effectively available to produce mechanical work, yielding a definition of the APE of the system.

Winters et al. considered the global effects of mixing (i.e., diabatic effects) on a closed system. For this purpose, they only needed a definition of APE for the system as a whole. There are however situations that call for a definition of APE that apply to localized regions of the domain, e.g. when measuring the  energy carried by nonlinear internal waves <cit.>, when partitioning energy between  mean and fluctuating components in a cyclone <cit.>, in determining  the efficiency of different mixing systems <cit.>, and in studying mixing and turbulence in spatially inhomogeneous systems <cit.>. In all these studies, the starting point was  the local definitions of APE developed in the early 1980's by <cit.> for incompressible flows and by <cit.> for compressible flows, based on a reference state that depends only on the mass distribution, though <cit.> already suggested that more general reference states can be considered. Indeed, in the atmospheric literature, more general reference states have been considered <cit.>, whereas in the oceanographic applications the Lorenz paradigm still dominates <cit.>. 
For a recent review see <cit.>. 

The inadequacy of the standard definition of energy as  the sum of kinetic, potential and, for compressible fluids, internal energy to quantify the capacity of the system to do actual work, is connected to the degeneracy of the Poisson algebra in  the Hamiltonian formulation  of the problem <cit.>. By degeneracy, we mean that the center of the Poisson algebra contains non-constant functionals of the phase space, the so-called Casimir functionals, or Casimirs for short. The set of Casimirs is an ideal of the algebra, and thus can be used to define a notion of equivalence
on the set of Hamiltonians:  Two Hamiltonians are equivalent, in the sense that they give the same dynamics, if they differ by a Casimir.  In other words, the Hamiltonian possesses a local (in phase space) gauge symmetry. 
From this point of view,   
the local APE formulations can be seen as selecting, out of a specific equivalence class, an Hamiltonian that satisfies one or more additional conditions <cit.>, given by a gauge-fixing condition.  As we shall clarify later, the Casimir includes the effects of constraints on the system. 

In this paper, we aim to revisit the issue of constructing a local diagnostic energy that can be used to diagnose the effect of diabatic processes on the energetics of fluid systems that, in the adiabatic limit, are described by a Hamiltonian with a degenerate Poisson algebra. In particular  we are seeking a quantity with the following properties:


  * Accounts for relevant  constraints.

  * Locality, and such that it satisfies (up to diabatic effects) suitable conservation laws.

  * Can be connected to Margules' intuitive notion of Available Energy when the latter is properly formalized.

  * Convexity in phase space, so that it can be meaningfully partitioned into a mean and eddy (or turbulent) component.

  * Its evolution under diabatic conditions reflects the loss of Available Energy. 

Following <cit.>, we first lay down in sec. <ref> the general theoretical framework that applies to generic systems that in the adiabatic limit have a Hamiltonian description characterized by a degenerate Poisson algebra. 
From there, we introduce the specific gauge-fixing condition that, when applied to the equivalence class of Hamiltonians, identifies the one whose density is the local diagnostic energy that we seek. 
In practice, the gauge-fixing condition identifies which Casimir needs to the added to the energy. At the same time, we also obtain the equations that specify the appropriate reference state. 

We then consider  two models which are commonly used in oceanography and which have an adiabatic limit described by a degenerate Hamiltonian structure: the incompressible shallow water equations and the incompressible Euler equations in the Boussinesq approximation for a continuously stratified flow. Both models are considered in inertial and non-inertial (i.e., rotating) frames.  For these systems, we give general properties for both the gauge-fixed Casimir and the reference state. In simple geometric configurations, we calculate analytically the solution or a suitable approximation. An interesting result that applies to both shallow water and Boussinesq equations is that, in rotating frames, the local diagnostic energy associated to low Rossby number reference states has a universal character. 



§ THEORETICAL BACKGROUND


The phase space which describes a Hamiltonian mechanical system characterized by a finite number of degrees of freedom typically, but not always <cit.>, has a non-degenerate Poisson algebra. In this case, if no other constraints exist, two Hamiltonians are dynamically indistinguishable, and thus belong to the same equivalence class, if they differ by at most a time-dependent constant. We can think of this as a global symmetry. 


If, however, the Poisson algebra is degenerate, then the equivalence class, as pointed out in the introduction, is wider. If that is the case,  the global symmetry, manifested by the invariance of the dynamics under the addition of  a (time varying) constant to the energy at every point  in phase space, widens to a local (in phase space) symmetry, i.e. the energy can be altered by the addition of a Casimir that can depend on the location in phase space, and still yield the same dynamical equations. From a field-theoretical point of view, the Hamiltonian in this case possesses a local gauge symmetry <cit.>. 






























It is important to remark that the existence of non-trivial Casimirs depends on the degeneracy of the Poisson algebra, rather than on the particular choice of the Hamiltonian. They represent conserved quantities whose existence does not depend via  Noether's theorem on symmetries of the Hamiltonian. If the algebra is degenerate, equilibrium solutions need not be (and in general, are not)  extrema of the Hamiltonian. 

 
 
 
For notational purposes, here and thereafter we will denote with H[q] the "naive" Hamiltonian, i.e. the sum of kinetic, potential and, when needed, internal energy, while H[q] will denote a  member of the class of equivalent Hamiltonians, which can be written as H[q]=H[q]+C[q], the sum of the naive Hamiltonian and a Casimir. 
 It is also convenient to introduce the following notation: if F[q] is a functional over the phase space which can be represented as the integral of an appropriate density n-form 𝔉(q) over the n-dimensional oriented manifold 𝔻, then F(q)≡⋆𝔉(q) is the corresponding scalar density (⋆ is the Hodge-star operator), i.e.

    F[q]=∫_𝔻𝔉(q)=∫_𝔻 F(q)𝔙,

where 𝔙 is the volume n-form. To summarize, square brackets [ ] denote a volume integrated quantity, with the corresponding local density indicated by the use of (). 
In fluid systems <cit.>, we find that 
 the Casimirs are generally underdetermined.  A set of Casimirs can be constructed as follows <cit.>: Let s_1(𝐱),⋯,s_p(𝐱) be p Lagrangian invariants, by which we mean they are  intensive quantities (0-forms) that are constant along Lagrangian trajectories, i.e. 

    ∂ s_i/∂ t+ L_ v(s_i)=0, i=1,…,p,

with L_𝐯 being the Lie derivative w.r.t. the flow induced by 𝐯. 
We also need a conserved density 𝔇, that is an n-form 
  representing a materially conserved extensive property, i.e.   such that 

    ∂𝔇/∂ t+ L_𝐯(𝔇)=0.

 For our purposes, we will use the mass n-form or, for incompressible flows, the volume n-form.   
Then 

    C[q]=∫_𝔻 C(s_1(𝐱),⋯,s_p(𝐱))𝔇,

where C(s_1,⋯,s_p) is any real-valued function which depends on its arguments algebraically (i.e., we are not allowing C to depend on derivatives of the s_i's).
It is trivial to verify that dC/dt≅ 0 (here and thereafter ≅ means equality up to boundary terms, or, in the case of n-forms, equality up to an exact form). Note that there may be more general classes of Casimirs (e.g., which may depend on derivatives of the conserved quantities). Here we limit to Casimirs that can be written as (<ref>).



 §.§ Holonomic brakes and Available Energy

It is straightforward to reformulate Margules' intuitive idea of "Available Energy" within the framework of a Hamiltonian system with a degenerate Poisson algebra. 

Consider a system described at time t=0 by a point q_0 in phase space. Under adiabatic conditions, the system is described by a trajectory q(t) in phase space such that along the trajectory H[q(t)]= H[q_0], where H[q] is any member of the class of equivalent Hamiltonian functionals. 

While <cit.> modified the Hamiltonian dynamics to dissipate the Casimirs, here we introduce into the system what we may call a holonomic brake, whose purpose is to extract energy from the system (hence a brake) in such a way that the s_i's
remain Lagrangian invariants (holonomic), and 𝔇 remains an integral invariant, i.e. (<ref>) and (<ref>) are still satisfied. Thus, the holonomic brake leaves the Casimirs invariant. Later, we will consider how holonomic brakes can be realized for particular classes of flows.






With the brake turned on, d H[q(t)]/dt≤ 0   and, since the brake is holonomic, any decrease in H[q(t)] is due, entirely,   to a decrease in the "naive" Hamiltonian H[q(t)]. We assume that, in the limit t→∞, the application of the brake causes the system to relax to a stationary "ground" state q_*. If such a ground state exists and is unique,  the energy extracted by the holonomic brake from the system which is initially at q_0, is the difference E_ AE[q_0]≡ H[q_0]-H[q_*]. We identify E_ AE with the Available Energy envisioned by Margules. We have thus the following definition:

The Available Energy E_ AE is the energy of the system referenced to the energy of the ground state, the latter being the state that the system settles to when subject to a holonomic brake.  

We shall see presently that, in the process of determining q_*, a gauge-fixing condition will naturally arise which selects an element H_* of  the equivalence class of Hamiltonians. The rate of change in time of H_*[q(t)] measures the loss of Available Energy when q(t) evolves under full diabatic conditions. 



 §.§ A special family of Casimirs: the Mass Distribution Function

The main postulate of our approach is that the q_* to which the system relaxes after the holonomic brake is turned on, is the point that minimizes the naive Hamiltonian subject to a set of constraints required by  the conservation of the s_i's along trajectories, which are not affected by the brake.  Thus, the extremal point of H[q] must be sought among the points that satisfy an extra set of conditions, the definition of which is the subject of this section. 


Consider a Lagrangian set of coordinates α=(α^1,…,α^n).
A Lagrangian particle characterized by a p-tuple s=(s_1,…,s_p) of Lagrangian Invariants retain its identity. The mass of all Lagrangian particles whose p-tuple falls between s and s+d s is given by

    Δ M( s)=(∫_𝔻[∏_i=1^p(s_i(α)-s_i)]𝔇(α))d s,

where 𝔇(α) is the mass  n-form expressed in Lagrangian coordinates, and  is the Dirac's delta (not to be confused with the δ used to denote variations). Integrating Δ M over the p-tuples, from s_i to its maximum value S_i≡max{s_i(α), α∈𝔻} for each i, we obtain the  Mass Distribution Function  (MDF)

    M( s)=∫_s_1^S_1⋯∫_s_p^S_pΔ M( s')=∫_𝔻[∏_i=1^pθ(s_i(α)-s_i)]𝔇(α),

where  

    θ(s)=∫_∞^s(t)dt,

is the Heaviside function. It represents the mass of all Lagrangian particles with values of the p-tuple greater than (s_1,⋯,s_n).[ 
<cit.> used a similar approach to construct a background state in an atmospheric context. 

]  
If we wish to express the theory in Eulerian coordinates, we just need to replace Eulerian coordinates x^i=x^i(α,t) in (<ref>) and integrate over the corresponding Eulerian mass density form





    𝔇( x)=(ρ(α)∂(α^1,⋯,α^n)/∂(x^1,⋯,x^n)(t))𝔙( x),

where 𝔙( x) is the volume form in Eulerian coordinates and ρ(α) the scalar mass density in Lagrangian coordinates. 

By construction, the MDF  is a p-parameter family of Casimirs, one for each p-tuple.  MDFs introduce a phase space "foliation": with the holonomic brake turned off, q(t) moves on a leaf of the foliation along isolines H[q]; when the brake is turned on, q(t) still moves on the same leaf, but crosses the isolines of H[q] on its way to q_*, which, by our hypothesis, is where the naive Hamiltonian attains its minimum on the leaf (see figure <ref>).

Two remarks are in order: While it may be possible to use (<ref>) to quantify the constraints, the integral version (<ref>) has the advantage of including bulk constraints, such as conservation of total mass; also, for incompressible flows in the Boussinesq approximation, we can replace the mass n-form with the volume n-form, in which case we should properly speak of a Volume Distribution Function. For simplicity, we will still refer to M( s) as the MDF regardless of what n-form is used. 

To solve the minimization problem, we introduce  a suitable Lagrange multiplier ψ( s), defined on the support of the MDF M^0( s) that specifies the leaf. That is, 
we seek 
to render the Lagrange functional

    L[q,ψ]≡ H[q]+∫ψ( s)[∫_𝔻∏_i=1^pθ(s_i( x)-s_i)𝔇( x)-M^0( s)]d s

stationary, 
by varying both q and ψ (recall that s_i(𝐱) is the value of the i-th Lagrangian invariant at 𝐱 calculated from the fields in q), that is we seek the solution[Though we should admit the possibility that more than one extremal (minimal) solution exists, for simplicity we assume that the minimal solution is unique.] q_*,ψ_*  of

    δ L/δ q≅ 0, δ L/δψ=0.

 If all the s_i's depend algebraically on the field in q, then the ≅ 0 can be replaced by equality stricto sensu.  
From now on, q_* will be referred to as the "ground state". 



 §.§ The gauge-fixing condition

Let 

    Ψ[q]≡∫_𝔻[∫ψ( s)(∏_i=1^pθ(s_i( x)-s_i))d s]𝔇( x),

be the Casimir associated to the Lagrange multiplier and let

    H[q]=H[q]+Ψ[q].

By switching the integral over the p-tuples s with the integral over the manifold 𝔻, the Lagrange functional can be rewritten as 

    L[q,ψ]=
     H[q]-∫ψ( s)M^0( s)d s.


Since the last term in (<ref>) does not depend on q, the Lagrangian multiplier ψ_* associated to the ground state q_* 
can be used to define the following  gauge-fixing condition to the equivalence class of Hamiltonians 

    .δ H/δ q|_q_*≅  0,

which is satisfied when H[q]= H_*[q]≡ H[q]+Ψ_*[q], where the Casimir Ψ_*[q] is related to the Lagrange multiplier ψ_* by (<ref>). 







We leave to the reader to verify that if two states q_1 and q_2 belong to the same leaf (i.e., they have the same MDF), then

    Ψ_*[q_1]=Ψ_*[q_2].

We are now in a position to define the local diagnostic energy as follows: 

The local diagnostic energy E is the scalar density of the Hamiltonian which satisfies  the gauge-fixing condition (<ref>), with the ground state q_* obtained from (<ref>) ,  i.e.

    E(q)=H(q)+Ψ_*(q).


 



In terms of the local diagnostic energy, 
the Available Energy  is 

    E_AE[q]≡ E[q]- E[q_*]=H[q]-H[q_*],

where the last equality follows from (<ref>), since q and q_* belong to the same leaf. Thus, the gauge-fixing condition selects, out of the class of dynamically equivalent Hamiltonians, the one that returns the same Available Energy as the naive Hamiltonian, when referenced to the ground state. Of course, the naive Hamiltonian is physically interpreted as  being the energy of the system.     
We also make the following conjecture: There exists a neighborhood of q_* over which  E[q] is a convex function (we implicitly assume that the phase space has a vector space structure). 
Note that (<ref>) is a necessary condition for the conjecture to hold. 
What makes the local diagnostic energy interesting is how the evolution of the Available Energy under full diabatic conditions is related to the local diagnostic energy.



 §.§ Evolution of the Available Energy under diabatic conditions

The Hamiltonian selected by the gauge-fixing condition (<ref>) acquires a special role when we apply  non-holonomic forces  which break the  invariance of the s_i's along Lagrangian trajectories. Without brakes, the motion occurs on a leaf along lines of constant energy. With a holonomic brake, the motion still occurs within the leaf, crossing the isolines if the naive Hamiltonian. Finally with the full non-holonomic brake on, the MDF becomes time dependent and  q(t) drifts across leaves. At each time, we can still calculate a ground state solving (<ref>-<ref>) with M^0( s) replaced by its instantaneous value M^t( s)  calculated from  q(t). Thus, starting from a time series of q(t) states, we obtain a sequence of ground states q_*(t) with the corresponding Lagrange multipliers ψ_*( s,t). Note that (<ref>) still holds on the pairs (q(t),q_*(t)).

Based on our definition of local diagnostic energy, 

    E[q_*(t)]=H[q_*(t)]+Ψ_*[q_*(t),t],

where the second argument in Ψ_*[·,t] is a reminder that the time dependence of the Casimir is due both to the time dependence of the ground state as well as to the intrinsic time dependence of ψ_*. The rate of change of the integrated local diagnostic energy of the ground state is 

    d E[q_*(t)]/dt=. (δ H/δ q+δΨ_*/δ q)|_(q_*(t),t)[d q_*/dt]+.∂Ψ_*/∂ t|_(q_*(t),t)≅.∂Ψ_*/∂ t|_(q_*(t),t),

where the last equality is made possible by the gauge-fixing condition.
Note that the last partial derivative only applies to the explicit dependence on time in the gauge-fixed Casimir[ In the following, when we speak of the Casimir, we will intend the gauge-fixed Casimir, as opposed to a generic Casimir. The context will usually suffice to make the distinction clear.  ]. Since q(t) and q_*(t) are on the same leaf, we have 

    d E[q_*(t)]/dt≅.∂Ψ_*/∂ t|_(q(t),t),

where now the Casimir is evaluated on  q(t), rather than on the ground state q_*(t). 
Therefore, the rate of change of the Available Energy 

    dE_AE[q(t)]/dt≅.δ E/δ q|_(q(t),t)[dq/dt]≅. (δ H/δ q+δΨ_*/δ q)|_(q(t),t)[d q/dt],

can be expressed solely in terms of the rate of change of the local diagnostic energy evaluated on q. 
The rate of change of the naive Hamiltonian reflects the action of the holonomic and non-holonomic components of the diabatic processes, while the change in the Casimir is due solely to the non-holonomic component. Interestingly, the non-holonomic component may actually increase the Available Energy, at least temporarily and/or locally. 



 The geometric structure of the leaf provides the "maximal" constraint:  the evolution of the system in the presence of non-holonomic forces can relax the constraint, and in so doing  frees up energy that would have been frozen in the ground state. 



§ SPECIFIC FLUID MODELS

 The theory developed in the preceding sections is quite general, and applies to any type of  flow regime that can be described in the adiabatic limit with a Hamiltonian structure having a degenerate Poisson algebra. 
 
In this section, we apply the theory to two types of flows: Shallow water flows, and 
continuously stratified flows that can be studied in the Boussinesq approximation. 
The first is chosen because it is often used as a first approximation to describe large-scale geophysical flows. Since it has only one Lagrangian invariant, it provides a good pedagogical introduction.  The second is commonly used to study low-Mach number flows in which vertical oscillations are confined to vertical scales much smaller than a suitably defined compressibility scale.  
While there are issues with the Boussinesq approximation, especially regarding the proper relation to  thermodynamics <cit.>, it nonetheless provides a good model to study flows in which the internal energy reservoir, in an approximate sense, merely acts as a sink of energy dissipated via diabatic effects. 
Most oceanographic models are based on the Boussinesq approximation, and therefore it is worthwhile to provide energy-based diagnostic tools applicable to Boussinesq flows.

The main result is that in the flows considered here,  the Casimir is related via a simple transformation to the Bernoulli invariant  of the ground state, and that in strongly rotating systems, to leading order, the Casimir is independent on the details of the MDF. These results can be obtained relatively easily in the shallow water case. The generalization to Boussinesq flows is more involved. For this reason, we will isolate the main results as theorems and corollaries to allow readers who are not interested in the derivation  to skip over the mathematical details.   



 §.§ Shallow water flows

Shallow water flows provide a relatively simple system where we can apply  the ideas outlined in the previous section. While relevant in its own from a Geophysical Fluid Dynamics point of view, it is mostly used as a  propaedeutic tool, having the advantage of possessing a single Lagrangian invariant, the potential vorticity. We will keep the discussion informal. 
A more formal approach will be reserved for the more interesting case of continuously stratified flows in the Boussinesq approximation. 

Shallow water flows are described by the pair (h,λ), where h is the local surface elevation, measured relative to a surface of constant geopotential height[In a rotating system, the geopotential height includes the "potential" of the centrifugal force.],  and,   λ= v^♭=v_idx^i the  1-form associated via the ♭ musical isomorphism to the velocity vector v=v^i∂/∂ x^i.  Here and thereafter, the Einstein convention of summation over repeated indexes is applied to Latin indexes.  Together, the pair (h,λ)  satisfies 


    ∂λ/∂ t=-i_ vΩ-d(gh+E_k),
    ∂ h𝔄/∂ t+ L_ v(h𝔄)=0,


where  𝔄  is the area element (a 2-form) of the manifold which contains the flow, E_k is the kinetic energy,  and
 the total vorticity Ω=dλ+Ξ is  the sum of the relative vorticity dλ and the frame vorticity Ξ. In an inertial frame Ξ=0.  . Because of (<ref>), the role of 𝔇 is played by h𝔄. 
For a flow in the shallow water approximation, the only Lagrangian invariant (up to a constant multiplicative factor) is the potential vorticity

    p≡ω/h,

where ω=⋆Ω is the Hodge-dual of the  vorticity Ω.  

The naive Hamiltonian  is

    H[q]=∫_𝔻 h(E_k+gh/2)𝔄.


The Casimir that satisfies the gauge fixing condition takes the form

    Ψ[q]=∫_𝔻 hΨ(p)𝔄.

A holonomic brake can be obtained adding to the r.h.s. of (<ref>) a term proportional to d(⋆ d(hΦ)). Physically, it represents a conservative (so as not to affect the vorticity) force field that extracts energy from surface waves.  


  §.§.§ Ground state and gauge-fixed Casimir

The ground state on a leaf characterized by the MDF M(s), and the Casimir that specifies the gauge-fixed Hamiltonian are found solving


    δ H/δ h=(E_k+g h+Ψ-pΨ_,p)𝔄=0,
    δ H/δλ≅ h(⋆λ)+dΨ_,p=0,
    ∫_𝔻 hθ(p-s)𝔄=M(s), p=⋆
    (dλ+Ξ)/h.


for the three unknown  h_*, λ_* and Ψ_*(p).    To avoid notational clutter, we denote with (·)_,s_1… s_p the p-th derivative w.r.t. s_1,…,s_p (note that the indexes are to the right of the comma). We assume that M is differentiable, and that M_,s<0. In other words, we assume that there is a continuous distribution of potential vorticity, so that with the introduction of a suitable second coordinate t, (p,t) can be used as local[It is important to realize that more than one chart may be needed to cover the entire manifold 𝔻.] coordinates over the manifold. 


Before considering actual solutions to (<ref>–<ref>), 
we point out a few general features. 


  * Eq. (<ref>-<ref>) is a system of three equations relating 4 quantities (λ_*,h_*,Ψ_*,M). The first two equations define a steady state solution of the shallow water equations in terms of Ψ_*. 
Indeed, from (<ref>) we have, with Φ_*≡-⋆λ_* being the flux form, 

    d(h_*Φ_*)=d(d(Ψ_*,p))=0,

by Poincare's lemma, thus (<ref>) is satisfied. 
Substituting dλ_*=-⋆ ph_*=-ph_*𝔄 in the l.h.s. of (<ref>) and using (<ref>) we have

    -i_ v(ph_*𝔄)+d(pΨ_*,p-Ψ_*)=-ph_*Φ_*+pΨ_*,ppdp=-pd(Ψ_*,p)+pΨ_*,ppdp=0,

thus satisfying (<ref>). 




  * If Ψ_* is a solution, so is Ψ_*+ap, where a is any constant. This indeterminacy simply reflects the fact that to fully determine the solution we need to specify the circulation on the boundary of the manifold. Indeed, the Casimir associated to ap is

    a∫_𝔻 hp𝔄=a∫_𝔻 dλ=a∫_∂𝔻λ=aΓ,

where Γ is the circulation along the boundary of the manifold. In particular, if the manifold is closed (i.e., has no boundary) Γ=0. 

  * The first equation tells us that the Casimir is related to the quantity  E_*k+gh_*  familiar from Bernoulli's theorem. Though it does not appear to have an agreed-upon name, for convenience we will refer to it as the Bernoulli invariant. 
Thus, if we know h_* and λ_*, say by integrating (<ref>-<ref>) with a holonomic brake until a steady state is reached, we obtain  the Casimir via 

    Ψ_*=p∫_p_ min^pE_*k(s)+gh_*(s)/s^2ds.




  §.§.§ The natural coordinates of the ground state

The form of the equations, and the considerations just made, suggest  that the ground state is more conveniently expressed with a "natural" set of coordinates, which in the present case are the potential vorticity and a second, time-like coordinate t. The latter is defined so that the relationship between "geometric" (x^1,x^2) and natural coordinates (p,t) is locally given integrating along the trajectories

    dx^i=v_*^idt,

to which we associate the space-time 1-forms ϕ^i≡ dx^i-v_*^idt (the subscript _* reminds us that these are ground state fields). 
To calculate the volume element in natural coordinates, consider the 2-form defined in space-time 

    𝔔=√(g)ϕ^1∧ϕ^2=𝔄-Φ∧ dt.

Here √(g) is the square root of the determinant of the metric in geometric coordinates. Let (p,t)ζ→ (x^1,x^2) be the (local) map from natural to geometric coordinates, and let (p,t)ζ̂→(ζ(p,t),t) be the map that embeds space into space-time. Under the pullback ζ̂^*𝔔=0 or
 
    𝔄=Φ∧ dt=Ψ_,pp/hdp∧ dt,

which gives the volume form in natural coordinates. Alternatively, we can define t as the second coordinate such that the r.h.s. of (<ref>) is the volume form of the manifold  (from now on, when multiplying simple differentials, we will follow standard practice and omit the exterior product symbol ∧). 
Incidentally, note that (<ref>) is integrable since 

    dϕ^i∧𝔔=0, i=1,2.


In natural coordinates, v=∂/∂ t, thus if Ω=ω dpdt is the vorticity, from (<ref>) we have the following integrability condition

    0=d(i_ vΩ_*)=-d(ω_* dp)=ω_*,t dpdt.


We want to consider under what conditions the ground state coincides with a constant geopotential surface, i.e. when dh=0 (recall that h is measured relative to the geopotential). Since the total volume is conserved, such a ground state has the minimum attainable potential energy consistent with conservation of volume. Setting h_*=const.  in (<ref>) we obtain

    0=ω_* dp-dEk_*.

The MDF is frame-invariant, by which we mean that observers in different frames of reference (inertial or otherwise) will measure the same MDF. 
Thus, we can always take the point of view of an  
inertial observer, for which the vorticity is given by 

    dλ=d(⋆Φ)=d(g_ptdp+g_ttdt)=(g_tt,p-g_pt,t)dpdt,

where the g_XY=g_ijx^i_,Xx^j_,Y  are the covariant components of the metric in natural coordinates.
Given that g_tt=g_ijx^i_,tx^j_,t=2E_k*, we have 

    ω_*=Ek_*,p-x^i_,pg_ij(x^j_,tt+Γ^j_knx^k_,tx^n_,t),

where the Γ^j_km's are the Christoffel symbols describing the Levi-Civita connection of the manifold 𝔻. Substituted in (<ref>) we obtain the following two conditions


    E_k*,t=0,
    
    x^i_,pg_ij(x^j_,tt+Γ^j_knx^k_,tx^n_,t)=0.


The first condition states that the kinetic energy must be constant along streamlines. A sufficient condition for the second condition to be satisfied is that the streamlines are geodesics. 



  §.§.§ The low-Rossby number ground state of shallows flows in rotating channels

A particularly interesting application of the observations made so far is to the case of a shallow water flow contained in a rotating channel, with Coriolis frequency f.  We seek a geostrophic ground state where the contribution of the relative vorticity dλ to the total vorticity is O( Ro) relative the frame vorticity Ξ, the latter in Cartesian coordinates being fdx^1dx^2, where Ro is the Rossby number.
In a geostrophic ground state, the kinetic energy must be O( Ro^2), and thus to leading order the Bernoulli invariant  is simply gh_*(p)=gfp^-1+O( Ro). 
The MDF of such a ground state is narrowly distributed around a non-zero p, i.e. 

    M(p)=hL_1L_2m(p-p/Δ p).

with Ro≡Δ p/p and h an average depth. Vice versa, it is possible to show that if the MDF is narrowly distributed, there is a rotating frame in which  the ground state is geostrophic. The interesting aspect of geostrophic ground states is that to leading order the Casimir is given by 

    Ψ_*=-fg/2p+O( Ro),

and thus to O( Ro)  the local diagnostic energy does not depend on the details of the MDF. 


  §.§.§ Non-holonomic dynamics on geostrophic leaves

The evolution of the Casimir of shallow-water geostrophic leaves under the influence of non-holonomic forces provides an interesting perspective on how the erosion of the potential vorticity constraint affects the Available Energy of a system.  For simplicity, we choose a simple model for the non-holonomic dynamics, which includes the effects of the bottom frictional Ekman layers on the vorticity dynamics <cit.>. Denoting with T the spin-down time of an Ekman layer, the potential vorticity satisfies 

    ∂ p/∂ t+ L_ v(p)=-⋆ dλ/Th,

i.e., the potential vorticity changes along trajectories by an amount proportional to the relative vorticity scaled with the local depth (Note that in a geostrophic flow, the r.h.s. is O( Ro pT^-1).). Thus, the change of the Casimir density along trajectories is to leading order

    ∂ hΨ_*/∂ t𝔄+ L_ v(hΨ_*𝔄)=Ψ_*,p(∂ p/∂ t+ L_ v(p))=-fg/2p^2⋆ dλ/T𝔄.

Therefore, Ekman layers with cyclonic relative vorticity  act as a sink of Available Energy along trajectories, whereas layers with anticyclonic relative vorticity are a source of Available  Energy. Of course, both cyclonic and anticyclonic layers act as a net sink on the naive energy.  


 §.§ Stratified Boussinesq flows

A continuously stratified flow in the Boussinesq approximation is described by the pair (b,λ), where b≡ g(ρ_0-ρ)/ρ_0 is the buoyancy,  and, as before,  λ≡ v^♭=v_idx^i satisfying 


    ∂λ/∂ t=-i_ vΩ-d(P+E_k)+bdZ,
    ∂ b/∂ t+i_ vdb=0,
    
    dΦ=0,


where P is the pressure, Φ=⋆λ=i_ v𝔙 is the flux form, and Z the geopotential height <cit.>. 
As before, the total vorticity Ω≡ dλ+Ξ is the sum of the relative vorticity dλ and the frame vorticity Ξ. Also, we assume that the geometry of the manifold 𝔻 that contains the fluid is Euclidean. Thus, there exists geometric coordinates (x^1,x^2,x^3) where the metric tensor assumes the simple form g_ij=1 if i=j and g_ij=0 if i≠ j.   
Also note that (<ref>) implicitly assumes that the volume form 𝔙 is constant in time. While the extension to non-inertial frames whose volume form is not constant in time is possible, in the interest of simplicity we do not pursue here. 

In a continuously stratified Boussinesq fluid there are two independent Lagrangian invariants, the buoyancy b itself and the potential vorticity p≡⋆ (db∧Ω).
Note that if f(b) is any (smooth) function of b, then f(b) is also a Lagrangian invariant, and so is ⋆(df(b)∧Ω). This will be reflected in certain degrees of freedom available in the definition of the gauge-fixed Casimir.  
The mass form 𝔇 can be replaced by the volume form 𝔙 of the manifold that contains the flow.  
 A holonomic brake  is provided by adding a term proportional to (i_ vdb)db to the r.h.s. of (<ref>), which extracts energy from the internal wave field without affecting the component of the vorticity aligned along the buoyancy gradient. 
 
 The naive Hamiltonian is 
 
    H[q]=∫_𝔻(E_k-bZ)𝔙,

 where E_k is the kinetic energy per unit mass,
 while Casimirs  have the general form 
 
    Ψ[q]=∫_𝔻Ψ(b,p)𝔙.

 In terms of the Lagrange multiplier, the density of the gauge-fixed Casimir is 
 
    Ψ(b,p)=∫_b_ min^b(∫_p_ min^pψ(q,s)dq)ds.

 
 

  §.§.§ The gauge-fixing condition

 The Frèchet derivatives of the naive energy w.r.t. to λ and b are trivial. We concentrate here on the variations of the Casimir written as in (<ref>). 
 The variation w.r.t. to b is given by
 
    δΨ/δ bδ b=(Ψ_,bδ b+Ψ_,pδ p)𝔙=(Ψ_,bδ b+Ψ_,p(⋆(d(δ b)∧Ω)))𝔙=
    Ψ_,bδ b𝔙+Ψ_,pd(δ b)∧Ω≅Ψ_,bδ b𝔙-d(Ψ_,pΩ)δ b=
     [(Ψ-pΨ_,p)_,b𝔙-Ψ_,ppdp∧Ω]δ b

  while the variation w.r.t. λ is given by 
  
    δΨ/δλ∧δλ=Ψ_,pδ p𝔙=Ψ_,pdb∧ d(δλ)≅ d(Ψ_,pdb)∧δλ.

  
 Given the MDF V(p,b), the ground state and the gauge-fixed Casimir are found solving
 
 
    δ H/δ b≅[-Z+(Ψ-pΨ_,p)_,b]𝔙-Ψ_,ppdp∧Ω= 0,
    δ H/δλ≅Φ+d(Ψ_,pdb)= 0,
    δ H/δψ=∫_𝔻[θ(p( x)-p)θ(b( x)-b)]𝔙=V(p,b)

 
 An attentive reader may have wondered why we did not enforce the incompressibility condition (<ref>) via an explicit Lagrangian multiplier. This is not necessary, as (<ref>) guarantees that the flux form of the ground state is closed. Also, as before, we assume that the MDF is a continuous and differentiable function, i.e. buoyancy and potential vorticity are continuously distributed. More complicated cases consisting of layers where buoyancy and/or potential vorticity are interleaved with layers where one or the other are constant will not be considered here. 
 
 

  §.§.§ The gauge-fixed Casimir
 
 The purpose of this section is to prove the following
 
Up to an inconsequential constant, the gauge-fixed Casimir is given by

    Ψ_*=p∫_ p_ min^pB_*(s,b)/s^2ds,

where B_*=Ek_*(p,b)-bZ_*(p,b)+P_*(p,b) is the Bernoulli invariant 
of the ground state expressed in natural coordinates. 

 
As we did for the shallow water case, we introduce the natural set of coordinates, which now include  the buoyancy,  the potential vorticity and the time-like coordinate, in terms of which the volume 3-form is 𝔙=Ψ_,pp  dbdpdt.  
In natural coordinates, we write the total vorticity as 

    Ω_*=Ω^bdpdt+Ω^pdtdb+Ω^tdbdp.

Since Ω is closed by definition

    (Ω^b)_,b+(Ω^p)_,p+(Ω^t)_,t=0.

From the definition of potential vorticity, 

    p=⋆(db∧Ω)=Ω^b⋆(dbdpdt)=Ω^b(Ψ_,pp)^-1.


We introduce Ψ_*=pΨ_*,p-Ψ_*. 
Taking the Hodge star of (<ref>) and using  (<ref>) we get


    Ψ_*,b=-(Z_*+Ω^p).
    Ψ_*,p=Ω^b,


where we have used the fact that  Ψ_*,p=pΨ_*,pp and (<ref>). These two equations can be combined into a single equation 

    d(Ψ_*+Z_*b)=bdZ-Ω^pdb+Ω^bdp=bdZ_*-i_ v_*Ω_*.

We have derived the above expression in natural coordinates, where v_* has the simple expression  v_*=∂/∂ t, but 
(<ref>) is coordinate-free. Thus,  comparing (<ref>) with (<ref>), we see immediately that Ψ_* is, up to a constant, the Bernoulli invariant of the ground state, and thus we have proven the theorem. <cit.>

 

  §.§.§ The apedic theorem and its corollaries

 Let us introduce a general 
  
  An apedic[From the Greek  a)/pedos, meaning level, flat.] ground state is a ground state such that b_*=b_*(Z). By extension, the leaf in phase space to which an apedic ground state belongs is called an apedic leaf.  
   
  To characterize under what conditions a ground state is apedic we have the following 
  
(Apedic theorem) Let the absolute vorticity of a ground state be 

    Ω_*=Ω^tdbdp+Ω^bdpdt+Ω^pdtdb.
 The ground state is  apedic if and only if


    Ω^p_,t=0,
    Ω^t_,t=0,


that is, the absolute circulation over loops embedded in surfaces of constant potential vorticity and the circulation over loops embedded in surfaces of constant "time" must be independent on "time" for a ground state to be apedic.

This follows from the integrability condition for (<ref>). Indeed, taking the external derivative of (<ref>) and leveraging the closedness of Ω, we obtain after some simple algebra 




    dZdb+Ω^t_,tdbdp+Ω^p_,tdtdb=0.

 The ground state must satisfy (<ref>). Since dpdt and dtdb are independent bases of the module of 2-forms over 𝔻, and  on an apedic manifold dbdZ=0 by definition, we have the proof of the apedic theorem.  

Let us now explore some consequences of the  apedic theorem. From (<ref>) we know that in an inertial system, the vorticity of the ground state

    Ω_*=dλ_*=d⋆(Ψ_*,ppdbdp)=d(g_btdb+g_ptdp+g_ttdt).

  For a manifold to be apedic, (<ref>) substituted in (<ref>-<ref>) requires that the metric satisfy  


    g_pt,bt-g_bt,pt=0,
    
    g_bt,tt-g_tt,bt=0.


In terms of Cartesian coordinates, the apedic condition becomes


    x^i_,px^i_,btt-x^i_,bx^i_,ptt=0,
    
    x^i_,bx^i_,ttt-x^i_,tx^i_,btt=0.


which provides the proof to  the first apedic corollary

A sufficient condition for a ground state to be apedic is that streamlines are geodesics and the velocity is uniform along streamlines.

Note that the Cartesian nature of the coordinates is of the essence. 
Consider, for example,  x^i's that are cylindrical coordinates (z,r,ϕ) and a ground state flow such that in its natural coordinates the x^i's at most depend linearly on t. Then, for example, we have from (<ref>)  that 

    Ω^t_,t=2ϕ_,tr(r_,bϕ_,pt-r_,pϕ_,bt),

which in general will not be equal to zero (and indeed, the trajectories are not geodesics!). In other words, the acceleration is relative to an inertial frame. Consequently, we can expect that when the system is non-inertial, a generic ground state will not be apedic. 
Also, it is worth remarking that the notion of apedicity, depends on the local direction of the plumb line, and thus it is not frame-independent. Thus, a fluid in solid body rotation will appear apedic to an observer in solid body rotation with the flow, but not to an inertial observer.









  §.§.§ The Casimir of inertial apedic ground states

Under the assumption that the ground state is apedic in an inertial system, we can make further analytic inroads. Let E_k=x^i_,tx^i_,t/2=g_tt/2 be the kinetic energy of the ground state. Assuming that apedicity is due to the trajectories being inertial (recall that the first apedic corollary proves only the sufficiency condition), we have 

    Ω^b=g_tt,p-g_pt,t=2Ek_,p-Ek_,p=Ek_,p,
    Ω^p=-g_tt,b+g_bt,t=-2Ek_,b+Ek_,b=-Ek_,b,

Then (<ref>-<ref>) can be combined

    d(Ψ-Ek+Zb)=bdZ=d(∫^Zb dZ),

and since we know that Ψ is the Bernoulli function of the ground state, we immediately deduce that the pressure distribution in apedic ground states is hydrostatic. 

The definitions of local Available Potential Energy proposed in the past that rely on a Casimir <cit.>, usually, <cit.> include only the potential and pressure term, combined in
∫^bZdb. This is 
 not surprising, since 
the traditional Available Energy approach  is to define the Background Potential Energy as the energy of the isochorically restratified flow. In our formalism, it amounts to enforcing conservation of buoyancy but not potential vorticity, which is to say we are minimizing the energy on a wider leaf. If we want to enforce both buoyancy and potential vorticity conservation, the corresponding ground state in general must have kinetic energy, and thus the system will have a lower Available Energy. 
This, however, does not exclude the existence  of interesting classes of flows for which the traditional approach is valid. They will be discussed in the next section,   
where we apply the theory to simple channel geometries,  both inertial and rotating. 




  §.§.§ Boussinesq Channel flows

For  arbitrary MDFs and geometries, (<ref>-<ref>) constitute a formidable nonlinear problem for the ground state. In these situations, to calculate the Casimir it is better to solve by other means (more likely, numerically) the Euler equations with a holonomic brake to find the ground state. Once the pair b_*,λ_* is known, the Casimir can be calculated from its Bernoulli invariant. 

Analytic progress can be made if we are willing to sacrifice geometric complexity. Thus, here we consider simple "open" channel geometries, which geometrically can be described by Cartesian coordinates. For definiteness, let 0≤ x^1≡ Z≤ L_z be the vertical direction, 0≤ x^2≡ y≤ L_y the spanwise direction and 0≤ x^3≡ x≤ L_x the streamwise direction.  The bottom, top and side walls are free-slip boundaries, while periodicity in streamwise direction.  
We seek what we may call 2-1/2 dimensional ground states: these are states such that x^i=x^i(b,p), i=1,2, and such that the velocity, which depends on two dimensions (b,p) is normal to the x^1,x^2 plane (the 1/2 dimension).





  
Inertial channels and the second apedic corollary. 
In inertial channels, we have the following

Let V be the MDF of a flow contained in an inertial channel which can be described by a 2-1/2 dimensional ground state as defined above. Then the ground state is apedic. 
Further,  if, for any value of b, V_,pb is an even function of p, the
kinetic energy of the ground state is zero.

The method to construct an apedic ground state for simple flow geometries extends the technique introduced by <cit.> and consists of a series of steps:


  – Determine the resorted height Z_*(b) solving (<ref>) in the limit p→ p_ min. This coincides with the isochoric restratification of the fluid. 

  – Compute the resorted distribution of potential vorticity along the spanwise direction y_*(p).

  – Finally, determine the kinetic energy E_k* of the ground state.


Later, we will modify this recipe to calculate the ground state and associated Casimir of rotating flows characterized by a low-Rossby number ground state, with the Rossby number being defined as Ro=O(ω_*/f) where ω_* is the relative vorticity and f is the Coriolis frequency.

 
It is straightforward to verify that 2-1/2 dimensional ground states in inertial channels are apedic.  Indeed, in geometric coordinates, the components of the velocity vector are [0,0,v_*(b(Z,y),p(Z,y))], thus dx=v(b,p)dt and dy=dZ=0, hence the velocity is uniform along geodesics.  Rather than working in "natural" coordinates, we work with "hybrid" coordinates (b,p,x). 
Being apedic Z_*=Z(b), (here, as we have done before, quantities used as a coordinates will not carry the _* subscript) which can be calculated from the limit p→ p_ min of (<ref>)

    ∫θ(b_*(Z)-s)dZdydx=-L_yL_x∫_b_ min^b_ Maxθ(b-s)dZ_*/dbdb=V(p_ min,s),

whose solution is 

    Z_*(b)=L_z(1-V(p_ min,b)/V(p_ min,b_ min)).

Its functional inverse b_*(Z) is the buoyancy distribution of the ground state. From this, we can calculate the hydrostatic pressure of the ground state.  
Also, for later convenience, we define N^-2_*≡ Z_*,b(b). 

Next, consider the volume element

    𝔙=dZdydx=|∂(Z,y)/∂(b,p)|dbdpdx.

Using again (<ref>) 

    |∂(Z,y)/∂(b,p)|=V_,pb/L_x,

from which y=y(p,b) is obtained solving 

    V_,pb/L_x=|∂(Z,y)/∂(b,p)|=|V_,b(p_ min,b)/L_xL_y∂ y/∂ p|.

The solution that spans the domain is

    y_*=L_yV_,b(p_ min,b)-V_,b(p,b)/V_,b(p_ min,b).

Note that in this apedic ground state, isopycnals are flat, whereas surfaces of constant potential vorticity are 
not and that the relationship between natural and geometric coordinates is invariant under the dilation
V→ e^s V.

At this point, all we have left is to calculate the kinetic energy of the ground state. By construction, we must have

    λ_*= Fdx,

and since Ω=dλ=F_,bdbdx+F_,pdpdx, using the definition of potential vorticity (<ref>) and (<ref>) we may try as a first Ansatz a solution to 

    F_,p=pV_,pb/L_x,

 with the constant of integration set along planes of constant buoyancy so that the total momentum along the x direction is conserved.

In the special case

    V_,pb=δ(p)Ṽ_,b,

which corresponds to a leaf characterized by a  uniformly zero distribution of potential vorticity,
F and the kinetic energy of the ground state are zero.  The Casimir then reduces (up to a factor pf(b) which does not concern us) to 

    Ψ_*=∫_b_ min^bZ(s)ds,

and the ground state to a quiescent fluid isochorically restratified according to (<ref>).
Two-dimensional flows contained in a vertical plane have zero potential vorticity, and thus  satisfy (<ref>).  In this case, we recover the same definition of Available Energy used by Winters et al.. 


Consider a flow at time t which evolved from a purely two-dimensional flow at t=0, i.e. at t=0 the MDF was V_,p=V_,bδ(p). As instabilities develop, the delta-like distribution in the MDF develops tails. Let us assume that in the absence of  any mechanism that breaks the symmetry, the potential vorticity that the diabatic processes generate ex nihilo  is equally distributed around the origin. In other words, we assume that whatever mechanism generates potential vorticity, it is not biased against either positive or negative potential vorticity. 
Then, without loss of generality, we can write for the MDF at time t as

    V_,pb≃L_yL_xN_*^-2(b)/σ(b)w(p/σ(b)),

where w(x) is a suitable even function. Analysis of numerical experiments indicate that a stretched exponential appears to give a good fit in different configurations (figure <ref>). As long as it decays sufficiently rapidly, we can replace p_ min with -∞ and p_ Max with ∞.  σ(b) measures the spread (i.e. the variance if the distribution was Gaussian) in potential vorticity near a given buoyancy level. Denoting with W(x)=∫_-∞^x tw(t)dt and with W_n=∫_-∞^∞W^n(t)dt, integration of (<ref>) yields for the velocity the following expression (setting the integration constant to zero for the moment) 

    F(p,b)=∫_-∞^ppV_,pb/L_xdp=L_y W(p/σ(b))N_*^-2(b)σ(b),

from which the kinetic energy per unit volume  of the ground state can be easily calculated 

    Ek_*=∫_-∞^+∞(∫_b_ min^b_ MaxF^2/2V_,bpdb )dp/∫_-∞^+∞(∫_b_ min^b_ MaxV_,bpdb )dp.=
    L_y^2/2(W_2/L_z∫_0^L_z(N_*^-2(b_*(Z)σ(b_*(Z)))^2dZ),

the inescapable consequence of which is that the kinetic energy per unit volume of this point on the leaf   grows with the width L_y of the channel squared. While we have derived this result assuming that the MDF has the particular form given by (<ref>), it fundamentally rests on the fact that under the dilation V→ e^sV we have F→ e^s F and Ek→ e^2s Ek. Therefore, what we have derived cannot be the ground state that we are seeking. To understand what went wrong, and to find a cure for it, in figure <ref> we plot F calculated using (<ref>) where the MDF was obtained from a DNS of stratified Couette flow <cit.> (indicated with the label 1).     

In this set up, the Lagrangian parcels with high negative potential vorticity are shunted on one side of the domain, most of the domain is filled by the (many) particles with little potential vorticity, while the parcels with large positive vorticity are pushed to the other side of the physical domain. We can see immediately two problems with this approach: First, having separated the large-negative   from the large-positive potential vorticity particles, we have created a large, sustained shear that forces the velocity to reach a large value in the center of the domain; and second, while this setup does not violate a free-slip boundary condition at the side walls, it cannot accommodate periodic boundary conditions. 


Indeed,  no matter what the details of the MDF are, a single (b,p,x)ζ→(z,y,x) chart cannot cover the manifold if the latter is periodic in the y direction. 
Let us explore then what happens if we allow two charts to cover the manifold. The first chart maps the (b,p,x) space to  the half of the channel such that 0≤ y≤ L_y/2, while the second covers the other half. Further, since the MDF is an extensive quantity by definition, we require that each half of the physical domain contributes exactly one half to the total MDF. In each half, we proceed exactly as we did before, and since F(p,b)p→±∞⟶0 it is possible to stitch together the velocity at the edge of every subdomain to obtain the profile shown on the right panel of figure <ref> labelled (2). Because in each half the MDF is half of the original one, the energy per unit volume in each half (an intensive quantity) is now 25% of the energy per unit volume obtained using only one chart. As a bonus, the velocity field is now periodic. Again, what makes this possible is the fact that V_,pb is an even function of p. 

Of course, there is no reason to stop at two charts. The procedure can be repeated with 4,8,…, charts, each covering 1/4,1/8,… of the channel, and each time obtaining a field with a kinetic energy per unit volume 1/16,1/64,… of the previous (figure <ref> shows the profiles obtained using one, two and four charts). Thus, we generate a sequence of velocity fields in which potential vorticity is interleaved at finer and finer scales, while at the same time reducing the overall kinetic energy by a factor 4  with each iteration. Barring the existence of a cut-off scale, we have proven the second apedic corollary.





Naturally, it is possible to consider initial conditions in which the symmetry in p is broken ab initio, even in the absence of rotation. For example, consider   filling a channel, which, as before, has dimensions (L_z,L_y,L_x) with a fluid whose  buoyancy 
has a uniform gradient along the spanwise direction db=N_0^2dy, and flowing such that λ=Δ v_0 (Z/L_z)dx, so that Ω=-Δ v_0/L_zdxdz. Thus, the potential vorticity is uniform and equal to -N^2_0Δ v_0/L_z. We can handle this singular limit by "smearing" the potential vorticity distribution around its constant value by an amount O(ϵ), carry out the calculation, and finally taking the ϵ→0 limit. In this case, the buoyancy of the ground state is such that   db_*=N^2_0(L_y/L_z)dZ and λ_*=Δ v_0(y/L_y)dx. Interestingly, the ground state has as much kinetic energy as the initial state, and  thus the Available Energy per unit volume is simply the difference between the potential energy of the initial and ground state (Δ b L_z/12). 





  
Rotating channels: geostrophic ground states. 

We want to consider here the ground states of flows contained in rotating channels where the frame vorticity is much larger than the relative vorticity.   In other words, we look for ground states characterized by a small value of the Rossby number Ro≡ O(ω/f), where ω is the order of magnitude of the relative vorticity and f the Coriolis frequency.  As we did in the shallow-water case, we  call low-Ro ground states geostrophic ground states. We work with natural (b,p,t) coordinates assuming that the ground state is 2-1/2 dimensional. 


Under these assumptions, the main source of non-inertiality along trajectories in the ground state is due to the frame vorticity, which we write as 

    Ξ=fdydx=Ξ^tdbdp+Ξ^bdpdt+Ξ^pdtdb,

with


    Ξ^t=f(y_,bx_,p-y_,px_,b),
    Ξ^b=fy_,px_,t,
    Ξ^p=-fy_,bx_,t,


where f is the Coriolis parameter.

Under the 2-1/2 assumption,  x_,t=√(2Ek), and (<ref>) becomes

    dΨ=f√(2Ek) dy-Zdb+O( Ro),

where we made use of (<ref>-<ref>). 

To O( Ro) the integrability condition 

    f√(2Ek)_,ZdZdy=b_,ydZdy,

shows that the ground state is in the so-called thermal wind balance <cit.>. As we did for flows in inertial channels, we endeavor to construct 
the ground state explicitly by seeking first  the coordinate maps Z_*=Z(b,p) and y_*=y(b,p). 

Along a streamline extending from one end of the domain to the other, ∫ dt=L_x/√(2Ek). From the definition of MDF

    L_x/√(2Ek)Ψ_,pp= V_,pb,

which, combined with (<ref>) gives

    fy_*,p=pV_,pb/L_x.

A second equation involving the coordinates is (<ref>),

    |∂(Z_*,y_*)/∂(b,p)|=Z_*,by_*,p-Z_*,py_*,b=V_,pb/L_x.

In choosing the sign, we assume that f>0. 
  We seek solutions of (<ref>-<ref>) in power series of (b-b_ min) as follows
  
  
    Z_*=∑_n=1^∞ Z_n(p)(b-b_ min)^n,
     
      y_*=∑_n=0^∞ y_n(p)(b-b_ min)^n,

  
 where we assume that the level of zero motion is the flat isopycnal b=b_ min. We also write
 
    V_,pb/L_x=∑_n=0^∞F_n(p)(b-b_ min)^n.

 as power series.[From a computational point of view, it may be more advantageous to use orthogonal polynomials. Here, we use a simple power series for convenience.] Inserting (<ref>) into (<ref>) we obtain the coefficients of the power series for y as 
 
    y_n=f^-1∫^p_p_ minsF_n(s)ds,

 while combining (<ref>) and (<ref>) 
 
    ∑_l=0^∞∑_n=1^l+1(nZ_ny_l+1-n,p-(l+1-n)Z_n,py_l+1-n)(b-b_ min)^l=∑_l=0^∞ F_l(b-b_ min)^l,

 which yields the following upper triangular set of algebraic equations for the coefficients Z_n of the series for Z
 
    ∑_n=1^l+1(nZ_ny_l+1-n,p-(l+1-n)Z_n,py_l+1-n)=F_l.

The equation for the first term  

    Z_1y_0,p=F_0

 can be readily solved using (<ref>)
 
    Z_1=f/p,

 while the other terms satisfy the following recurrence relation
 
    Z_l+1=[(l+1)y_0,p]^-1[lZ_1,py_l-∑_n=2^l(nZ_ny_l+1-n,p-(l+1-n)Z_n,py_l+1-n) ], l≥ 1.

 Note how Z_l+1 depends on the y_n's with n≥ 1.  
 
Having determined the map from natural to physical coordinates, we now calculate the kinetic energy of the ground state, by solving  
 the thermal wind equation (<ref>), which expressed in (b,p) coordinates reads
 
    (√(2Ek))_,py_,b-(√(2Ek))_,by_,p=f^-1Z_,p.

 Once again, we look for a power series solution
 
    √(2Ek)=∑_l=0^∞v_l(p)(b-b_ min)^l,

and, since we have assumed that b=b_ min is the level of no-motion, v_0=0. 

Substituting into (<ref>) and rearranging we again obtain a triangular algebraic set of equations for the v_l's (no surprise here, as we are dealing with the same nonlinear operator)


    v_1=0,
    
    v_l+1=[ -Z_l,p+f∑_n=1^l (nv_l+1-n,py_n-(l+1-n)v_l+1-ny_n,p)]/[f(l+1)y_0,p],l≥ 1.

 
As before, the v_l's with l≥ 3 depend on the y_n's with n≥ 1. 
In particular,  
 
    fy_0,pv_2=f/2p^2.

 An  special case is when  F_n=0, n>0, i.e. the VDF is linear in the buoyancy. In this case, the solution has the very simple form
 
    y_*=f^-1∫^p_p_ minsF_0(s)ds, Z_*=fb-b_ min/p, √(2Ek_*)=f(b-b_ min)^2/2p^2[pF_0(p)]^-1.











 
 With the complete solution in hand, it is now time to revisit the small-Ro assumption made at the beginning. 
 The goal is to characterize the MDFs that lead to small-Ro ground states.  
 It is straightforward to calculate the relative vorticity of the ground state 
 
    Ω_*=-(2Ek)_,bdtdb+(2Ek)_,pdpdt.

 Let us start from the special case  that leads to (<ref>). The most stringent condition is   
 
    Ro=O((√(2Ek_*))_,Z/f)=L_z/f^2pF_0(p)≪ 1

 We can estimate 
 
    O(pF_0)=pL_yL_z/Δ pΔ b=pL_y/Δ pN^2,

 where N^2=Δ b/L_z. Here Δ p is a measure of the spread of potential vorticity, e.g. the rms of the fluctuations,  Δ b is a measure of the spread of buoyancy and p the mean potential vorticity (which, when needed, can be estimated as fN^2).   
 Thus, the ground state satisfies the low-Ro condition if the following condition on the spread of PV 
 
    Δ p/p≪L_yL_z/L_R^2,

 where L_R≡ L_zN/f is the internal Rossby radius of deformation, is satisfied. 
 
 Let us now consider the more general case of a MDF which is not linear in the buoyancy, but which still satisfies (<ref>). We have
 
    O(pF_n)= Ro^-1L_R^2/Δ b(Δ b)^-n,

 whence
 
 
    O(fy_n)=O(Δ p  pF_n)=fL_y(Δ b)^-n,
    
     O(fy_n,p)=O(pF_n)= Ro^-1L_R^2/Δ b(Δ b)^-n.

 
Using (<ref>), we have 

    O(Z_2(b-b_ min)^2)= RoL_yL_z/L_R^2L_z.

Given the nature of the recurrence relation (<ref>), it  follows that all the other terms in the series 
for Z are O( Ro),  and thus 

    Z_*=f/p(b-b_ min)+O( Ro).

Along essentially the same lines we have that  √(2Ek_*)=O( Ro), and thus to O( Ro) the pressure in the ground state is hydrostatic. 

Hence

    Ψ=-f/2(b-b_ min)^2/p+O( Ro),
























which shows that the Casimir of geostrophic leaves, defined as leaves characterized by a geostrophic ground state (i.e., a low-Ro state)  does not depend, to leading order, on the details of the MDF that defines the leaf to which the ground state belongs. 
A field on a geostrophic leaf does not need to be in geostrophic equilibrium, but the ageostrophic component must be such that the overall spread in potential vorticity is small when looked over the entire field.    
 As long as the diabatic dynamics does not cause a significant spread in potential vorticity, the flow remains on a geostrophic leaf and the gauge-fixed Hamiltonian is to leading order independent on the MDF. Thus, the local diagnostic energy has, to O( Ro) a universal character. 
 

 §.§ The non-holonomic dynamics of Casimirs in Boussinesq flows


Once the gauge-fixed Casimir is known, non-holonomic effects on the Available Energy can be analyzed by considering how the Casimir evolves along Lagrangian trajectories. In apedic manifolds characterized by an even distribution of potential vorticity, our theory recovers the standard SBB local APE formulation (here, we reinstate the qualified Potential, since the ground state has no kinetic energy). In particular, under Boussinesq approximation, the sink of APE is quantified by -κψ_,bb|∇ b|^2, where κ is the diffusivity of the stratifying agent. Archetypal flows for small-scale mixing (e.g., shear layers, Couette flows,…) fall into this category. On geostrophic manifolds, the evolution of the Casimir depends on the details of the non-holonomic brake. Assuming a standard diffusion term in the buoyancy equation, the sink term due to mixing of the stratifying agent -κψ_,bb|∇ b|^2 is equal to κ f/p|∇ b|^2+O( Ro). However, turbulent momentum fluxes will change the potential vorticity along trajectories and thus modify the Available Energy in either direction, just as we saw for the shallow-water case. The specific way depends in general on how the turbulent fluxes are modeled, and details are left as future work. 



§ RELATION TO OTHER DEFINITIONS OF LOCAL AND GLOBAL AVAILABLE ENERGY

 Examples of global <cit.>  as well as local Available Energy definitions <cit.> applied to study the effects of mixing on the energetics of three-dimensional flows have been previously considered in the literature. 

At the level of a global definition of Available Energy, the novelty of the approach followed in sec. <ref> lies in setting a framework that allows, at least in principle, to calculate the ground state from the MDF that describes the state of a generic state, though it can also be used to calculate the MDF associated to a prescribed ground state. While the framework is amenable to analytic treatment in some special cases, it suggests a way to calculate numerically the ground state via the application of suitably defined holonomic brakes. Once the ground state is known, the gauge-fixed Casimir can be calculated, and thus the local diagnostic energy is obtained. 
Out of the cases that can be treated analytically, we show that for an important class of flows, the  isochoric restratification proposed by Winters et al. coincides with the  ground state obtained with our framework. 



The major departure from previous local formulations of APE 
<cit.>, aside from the way the Casimir is calculated, resides in not subtracting the contribution of the ground state at the local level. If we were to follow Scotti and White's "recipe", we would define the local diagnostic energy as E^ SW(q)= H_*(q)- H_*(q_*), with H_* given by the gauge-fixing condition (<ref>). 


Subtracting the gauge-fixed Hamiltonian density of the ground state  in the local definition  has an intuitive appeal, because the global Available Energy would coincide with E^ SW[q], whereas in the approach followed here the integral of the local diagnostic energy does not have an immediate physical interpretation. However, the quantity of interested is the rate of change of the Available Energy under diabatic conditions, which coincides with the rate of change of  E[q]. Thus, the local diagnostic energy allows to explore at the local level the processes that are responsible for energy degradation, and at the local level there are reasons to question the inclusion of the ground state energy in the local definition. 

From a functional point of view, the inclusion of H_*(q_*) amounts to no more than a "constant" term, in the sense that its variation w.r.t. q is zero, and it is not necessary to impose either the gauge condition or to ensure   the convexity property.

There is also a subtler reason why it is not advisable to include q_* in the definition at the local level.  Consider the local evolution of H_*(q_*) under the flow v associated to the point q in phase space. For simplicity, we consider the stratified  Boussinesq case, where the volume form is an integral invariant. Then H_*(q_*)= H_*(λ_*,b_*,Z) (we do not need to worry about the explicit time dependence of H_*) and 

    ∂( H_*𝔙)/∂ t+ L_ v( H_*𝔙)=
     δ ( H_*𝔙)/δλ_*(∂λ_*/∂ t+ L_ v(λ_*)) 
    
    +δ ( H_*𝔙)/δ b_*(∂b_*/∂ t+ L_ v(b_*))
    +δ ( H_*𝔙)/δ Z L_ v(Z)
    ≅δ ( H_*𝔙)/δ Z L_ v(Z),

where, as before, ≅ means equal up to an exact form. For an apedic manifold, 

    δ ( H_*𝔙)/δ Z L_ v(Z)=-b_*(Z)d(ZΦ)=-d(P_*(Z)Φ)≅ 0, Φ=i_ v𝔙, P_*=∫^Z b_*(s)ds,

so that in the adiabatic limit, E^ SW is, up to boundary terms, an integral invariant.  

However, under more general ground states (e.g., geostrophic), E^ SW is not an integral invariant, whereas E is. Therefore, referencing the local diagnostic energy to the local background state does not contribute to the desirable properties of the local diagnostic energy, while, aside from special circumstances, adding terms to the local evolution equation of dubious physical interpretation.








§ CONCLUSION




In this paper, we have developed a framework to diagnose  diabatic effects on systems that in the adiabatic limit are described by a degenerate Hamiltonian structure. The degeneracy has two important consequences. The first is that the phase space foliates into a collection of leaves, each labeled by a particular volumetric distribution of the Lagrange invariants associated to the Casimirs. Trajectories in phase space are contained within a leaf. The Available Energy on a leaf is thus the energy referenced to the minimum energy attainable on the leaf. In other words, it is the energy that can be extracted from the system without disturbing the dynamics of  the Lagrangian invariants. 

The second consequence is that the Hamiltonian that describes the flow in the adiabatic limit possesses a local (in phase space) gauge symmetry. Thus, there exists a class of dynamically equivalent Hamiltonians. 

We show that a specific gauge-fixing condition can be imposed to select a specific Hamiltonian that can be used to define a local diagnostic energy with the following property: the temporal rate of change of the Available Energy  under diabatic conditions is given by the diabatic evolution of the gauge-fixed Hamiltonian.


The present study provides expressions of the gauge-fixed Hamiltonians for a number of geophysically relevant flows, together with general properties of the reference state. 
For non-rotating flows in the Boussinesq limit, we have shown that the local diagnostic energy based on the gauge-fixed Hamiltonian recovers the Available Potential Energy of <cit.>, which form the basis of <cit.> and <cit.> application of APE to mixing, provided the distribution of potential vorticity in the flow is even around the origin. Most archetypal flows used to study small-scale mixing start as two-dimensional flows contained in a vertical plane. For such flows, under the assumption that the diabatic dynamics does not favor negative over positive potential vorticity,  the standard  approach remains correct. 


Our definition of local diagnostic energy extends naturally to rotating flows: under the effect of rotation,  the conservation of potential vorticity imposes a strong constraint, which is not accounted for in the Winter et al.'s definition of APE. In particular, when the flow is obtained from a general-form perturbation  which preserves both the potential vorticity and the buoyancy distribution of a low-Rossby number state in near-geostrophic equilibrium,
to lowest order in the Rossby number, we show that the expression for the local diagnostic energy has a universal character.

Once the functional form of the local diagnostic energy is known, it can be used to study the effect of diabatic processes. The effect of these processes, which in general  break the conservation of the Lagrangian invariants on the change of the Available Energy, can be studied  following the evolution of the local diagnostic energy along Lagrangian trajectories.
Using a simple model for diabatic processes in  shallow-water flows, we have shown that such processes can locally increase the Available Energy. 
This is not too surprising, since diabatic processes erode the potential vorticity constraint, keeping energy locked up in the ground state field. 


The results were derived under the assumption that the ground state is unique. Of course, we can envision situations where the naive Hamiltonian has multiple minima on a given leaf. In this case, it cannot be ruled out that the application of different holonomic brakes may lead the different ground states, especially if the system is close to the separatrix between the basin of attractions. On the other hand, if the system is not too far from a minimum, then it is not unreasonable to assume that most holonomic brakes will lead to the same extremal point, even though there may exists other ground states with overall less energy.    


We considered applications of the general formalism developed in section  <ref> to flows in which thermodynamic effects are neglected. 
However, the formalism is quite general, and can be applied to any flow regime described by a Hamiltonian on a phase space with a degenerate Poisson algebra. In particular, it could be applied to flows where the internal energy plays a more active role, such as low Mach number applications <cit.>.

In conclusion, the formalism  presented in this paper extends Margules' original idea to highlight the role that multiple constraints play in the available energetics framework.

Since the Mass Distribution Function, or, for incompressible flows, the Volume Distribution Function plays a central role in defining the gauge-fixing condition, our approach applies to flows contained on a finite volume. 
 However, the possibility that other ways to enforce the constraints, not dependent on the system having a finite volume, exist should not be discarded. 
Since  the applicability of any local diagnostic tool is constrained by the data available for the analysis, it can rarely be considered independent from the flow data itself 
6pt

20pt

The research described in this paper did not involve any living or deceased organism.



AS conceived the study and drafted the manuscript. AS and PYP worked together on the calculations, PYP heavily edited the manuscript and contributed some of the figures. 

We have no competing interests.

The authors acknowledge the support by the National Science Foundation Grant N^o OCE-1155558.

We would like to thank Dr. K. Lamb, Dr. A. Hoggs and Dr. R. Tailleux for numerous engaging discussions on the subject. 
AS would like to thank Dr. E. Santilli for opening his eyes to the beauty of exterior calculus. 

Insert disclaimer text here.


rspublicnat



