Available potential energy equation

AvailablePotentialEnergy(
    model;
    method,
    geopotential_height,
    location
)

Return a KernelFunctionOperation computing the local available potential energy density,

    eₐ(b, z) = ∫_{z✶}^{z} [b✶(z̃) - b] dz̃ ,   equivalently   (g/ρ₀) ∫_{z✶}^{z} [ρ - ρ✶(z̃)] dz̃

the work needed to bring a parcel from the reference height z✶ it would occupy in the adiabatically resorted state to the height z where it actually sits. The parcel's own buoyancy b is held fixed along the path; only the reference profile b✶ varies with z̃.

This is the spatially local APE density of Holliday & McIntyre (1981) and it is also used in Wenegrat, Chor & Barkan (2026) as a basis for a filtered APE framework. It is non-negative everywhere whenever the reference profile is one-dimensional and gravitationally stable.

z✶ is the reference height computed by reference_height; pass one explicitly to share a single sort with BackgroundPotentialEnergy, or pass method through to choose how it is built. All four give the same eₐ volume integral when each sorts the field itself.

Integral(eₐ) recovers the global APE Integral(PotentialEnergy(model)) - Integral(BackgroundPotentialEnergy(model)) only in the continuum limit: the local density samples the reference profile at the model's cell centers while the global split effectively samples it at the sorted column's, and the two midpoint quadratures differ at finite Δz. The gap is second order in the vertical spacing, so it is a fraction of a percent on a well resolved grid but a few percent on a coarse one. eₐ does vanish, cell by cell and exactly, for a statically stable and horizontally uniform stratification.

The result lives at (Center, Center, Center), per unit mass (units m² s⁻²), and is defined for the same buoyancy formulations as PotentialEnergy. Under VerticalSort it lands on the sorted column, indexed by rank rather than by position in the flow.

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)AvailablePotentialEnergy(model)# outputAvailablePotentialEnergy KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: local_ape_ccc (generic function with 2 methods)└── arguments: ("Field", "Field", "Field")└── computes: local available potential energy density  eₐ = ∫[b✶(z̃) - b]dz̃ ≥ 0

BuoyancyDisplacementPotential(
    model;
    method,
    geopotential_height,
    location
)

Return a KernelFunctionOperation computing the buoyancy displacement potential

    Υ = z✶ - z

how far below its actual height a parcel's reference height sits, and so how far it would have to travel to reach the adiabatically resorted state. It is the derivative of the local available potential energy with respect to buoyancy, Υ = ∂eₐ/∂b, which is what makes it the natural conjugate of b: contracting it with a buoyancy gradient gives an APE dissipation rate (AvailablePotentialEnergyDissipationRate), and contracting it with a sub-filter buoyancy flux gives a cross-scale APE flux.

This is the buoyancy form of Υ(ρ, z) = g(z - z✶(ρ))/ρ₀ as Wenegrat, Chor & Barkan (2026) write it in their Eq. (7) for density. The two differ by the factor -g/ρ₀ that converts between buoyancy and density, which cancels wherever Υ is contracted with a buoyancy gradient. The result lives at (Center, Center, Center) and is a length (units m).

z✶ is the reference height computed by reference_height; pass one explicitly to share a single sort with the other reference-state diagnostics, or pass method through to choose how it is built. It has to be one that lives on the model grid, so VerticalSort is rejected.

HeavisideIntegral is the default here rather than the package-wide ThreeDimensionalSort because Υ is a map, and every use of it differentiates that map. Only Eq. (11) of Winters et al. makes z✶ a function of buoyancy alone, so tied cells share one reference height; with ThreeDimensionalSort a run of equal buoyancy takes consecutive slots and spreads z✶ over the depth it fills, which is harmless in a volume integral but shows up in ∇Υ as grid-scale noise.

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)BuoyancyDisplacementPotential(model)# outputBuoyancyDisplacementPotential KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: upsilon_ccc (generic function with 1 method)└── arguments: ("Field",)└── computes: buoyancy displacement potential  Υ = z✶ - z

AvailablePotentialEnergyDissipationRate(
    model;
    method,
    geopotential_height,
    location
)

Return a KernelFunctionOperation computing the rate at which diffusion destroys available potential energy,

    ε_A = κ ∂ᵢb ∂ᵢΥ = κ [(∂z✶/∂b)|∇b|² - ∂b/∂z] ,

the sink of the local available potential energy equation of Wenegrat, Chor & Barkan (2026) (their Eqs. 11 and 14, where it appears as -ε_A), with Υ the BuoyancyDisplacementPotential. It follows from ∂eₐ/∂b = Υ, which makes the diffusive part of Deₐ/Dt equal to Υκ∇²b = ∇·(κΥ∇b) - κ∇Υ·∇b: once the flux divergence is set aside, ε_A = κ∇Υ·∇b is what remains.

Written out, the first part is the diapycnal mixing rate of Winters et al. (1995), the work done rearranging the reference state, and the second is PotentialEnergyDiffusiveVerticalBuoyancyFlux, the diffusion that state undergoes on its own, which carries no APE with it. The two cancel exactly for a statically stable, horizontally uniform stratification, where z✶ = z and there is no available energy to destroy, so ε_A measures only the APE actually lost — it is not the sign-definite κ|∇b|²-like quantity the name might suggest.

κ ∂ᵢb is taken from the closure's own diffusive flux rather than from a diffusivity supplied here, so this follows whatever closure the model runs with, and is written in the same conservative form TracerVarianceDissipationRate uses. The result lives at (Center, Center, Center), per unit mass (units m² s⁻³).

The buoyancy has to be a tracer the closure diffuses, so this is defined for BuoyancyTracer models only — SeawaterBuoyancy would need the diffusive fluxes of temperature and salinity combined through the equation of state.

upsilon is a keyword of the two-argument form only. Anyone holding a Υ also holds the z✶ it was built from, and passing that is both cheaper and unambiguous: this form would have to sort the domain to build a z✶ it then uses for nothing but a grid check.

z✶ is the reference height computed by reference_height, and has to be one that lives on the model grid, since ∇b is taken there; HeavisideIntegral is the default for the reason BuoyancyDisplacementPotential gives. upsilon takes a Υ you already have, so that writing both out costs one sort and one Υ rather than two of each:

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(κ=1e-4))z✶ = reference_height(model, method=HeavisideIntegral())Υ = Field(BuoyancyDisplacementPotential(model, z✶))AvailablePotentialEnergyDissipationRate(model, z✶; upsilon=Υ)# outputAvailablePotentialEnergyDissipationRate KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ape_dissipation_rate_ccc (generic function with 1 method)└── arguments: ("Field", "ScalarDiffusivity", "Nothing", "Val", "Field", "Clock", "NamedTuple", "BuoyancyForce")└── computes: available potential energy dissipation rate  ε_A = κ ∂ᵢb ∂ᵢΥ