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 the quantity the filtered APE framework is built on. 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  ∫[b✶(z̃) - b]dz̃ ≥ 0

AvailablePotentialEnergyDisplacementPotential(
    model;
    method,
    geopotential_height,
    location
)

Return a KernelFunctionOperation computing the 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 subfilter buoyancy flux gives a cross-scale APE flux.

This is the buoyancy form of the displacement potential written for density as Υ(ρ, z) = g(z - z✶(ρ))/ρ₀. 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.

using Oceanostics.AvailablePotentialEnergyEquation additionally brings in DisplacementPotential, short enough to read beside the other terms of an eₐ budget. That alias is scoped to this module, as DissipationRate is: using Oceanostics does not bring it in, since unprefixed it says nothing about which budget's displacement it names.

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)AvailablePotentialEnergyDisplacementPotential(model)# outputAvailablePotentialEnergyDisplacementPotential 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: displacement potential  z✶ - z

ReferenceBuoyancyAnomaly(
    model;
    method,
    geopotential_height,
    location
)

Return a KernelFunctionOperation computing the buoyancy anomaly a parcel carries relative to the adiabatically sorted reference profile taken at the parcel's own height,

    bᵣ = b - b✶(z) ,

the buoyancy form of the anomaly written for density as b_r(ρ, z) = -g(ρ - ρ✶(z))/ρ₀. It is what the available potential energy exchanges with the kinetic energy (AvailablePotentialToKineticEnergyConversion), and the buoyancy-space counterpart of AvailablePotentialEnergyDisplacementPotential: Υ = z✶ - z measures a parcel's displacement from the reference state as a height, bᵣ measures it as a buoyancy, and both vanish exactly where the fluid is already sorted.

Note that b✶(z) is the reference profile at the height the parcel actually occupies, which is not reference_buoyancy: that pairs the profile with z✶ instead, and b✶(z✶) is the parcel's own buoyancy, so the anomaly built from it would be zero everywhere. The profile is sampled by reference_buoyancy_at_height, which reads it off the sort rather than repeating it.

The result lives at (Center, Center, Center) and is a buoyancy (units m s⁻²). 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.

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

AvailablePotentialToKineticEnergyConversion(
    model;
    method,
    geopotential_height,
    location
)

Return a KernelFunctionOperation computing the rate at which available potential energy is converted into kinetic energy,

    w bᵣ = w [b - b✶(z)] ,

the exchange term of the local available potential energy equation, which the eₐ budget takes with a minus sign and the kinetic energy budget with a plus.

This is not PotentialToKineticEnergyConversion, which computes uᵢbᵢ over the total buoyancy (wb under the vertical gravity these diagnostics require) and belongs to the eₚ budget. The two differ by w b✶(z), the exchange between the kinetic energy and the background state, and that difference is why the pressure that goes with this budget is the deviation from the hydrostatic pressure of the reference profile. As fields they are different maps; w b✶(z) is a flux divergence, so the two only agree once integrated over a periodic or closed domain:

    ∫ w bᵣ dV = ∫ w b dV .

The anomaly bᵣ is ReferenceBuoyancyAnomaly, built here unless one is passed as anomaly, which is worth doing when both are wanted, since a Field of it is then computed once. z✶ is the reference height computed by reference_height; it has to be one that lives on the model grid, so VerticalSort is rejected. The result lives at (Center, Center, Center) and is a conversion rate per unit mass (units m² s⁻³).

using Oceananigans, Oceanosticsgrid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)z✶ = reference_height(model, method=HeavisideIntegral())bᵣ = Field(ReferenceBuoyancyAnomaly(model, z✶))   # share the anomaly between the twoAvailablePotentialToKineticEnergyConversion(model, z✶; anomaly=bᵣ)# outputAvailablePotentialToKineticEnergyConversion KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: ape_to_ke_conversion_ccc (generic function with 1 method)└── arguments: ("Field", "Field")└── computes: available potential to kinetic energy conversion  wbᵣ

AvailablePotentialEnergyDissipationRate(
    model;
    method,
    geopotential_height,
    location
)

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

    εₐ = -qᵢ ∂ᵢΥ = -(∂z✶/∂b) qᵢ ∂ᵢb + q₃ ,

the sink of the local available potential energy equation, where it appears as -εₐ, with Υ the AvailablePotentialEnergyDisplacementPotential and qᵢ the diffusive buoyancy flux the closure supplies. It follows from ∂eₐ/∂b = Υ, which makes the diffusive part of Deₐ/Dt equal to -Υ ∂ᵢqᵢ = -∂ᵢ(Υqᵢ) + qᵢ∂ᵢΥ: once the flux divergence is set aside, εₐ = -qᵢ∂ᵢΥ is what remains. Nothing here assumes a form for qᵢ: it is -κ∂ᵢb for Fickian diffusion, and whatever an LES closure returns otherwise.

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 εₐ measures only the APE actually lost — it is not the sign-definite buoyancy-variance-like quantity the name might suggest.

qᵢ 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 AvailablePotentialEnergyDisplacementPotential 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(AvailablePotentialEnergyDisplacementPotential(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  -qᵢ∂ᵢΥ