Subfilter available potential energy equation

SubFilterAvailablePotentialEnergy(
    model,
    filter;
    method,
    geopotential_height
)

Return the subfilter-scale (SFS) available potential energy eₐˢ, the available potential energy carried by the scales that a low-pass filter removes from the buoyancy field — the filtered full APE minus the APE of the filtered buoyancy b̄ = filter(b):

    eₐˢ = filter(eₐ(b, z)) - eₐ(b̄, z) ,   eₐ(b, z) = ∫_{z✶(b)}^{z} [b✶(z̃) - b] dz̃

where eₐ is the local available potential energy density (AvailablePotentialEnergy), eₐ(b̄, z) is the APE of the filtered buoyancy (FilteredAvailablePotentialEnergy), and both terms are measured against one shared reference profile (b✶, z✶): the filtered buoyancy is looked up in the same profile the full field is measured against, which is exactly what ProfileLookup was built for. It is the potential-energy counterpart of the subfilter kinetic energy eₖˢ (SubFilterKineticEnergy), just as FilteredAvailablePotentialEnergy is that of the filtered kinetic energy eₖˡ.

Because the two states have to share one profile, method must be a ProfileLookup:

• ProfileLookup() (the default) sorts the profile from the model's own buoyancy — a VerticalSort column built internally and re-sorted on every compute!, so the reference state tracks the flow.
• ProfileLookup(z✶_column) borrows a column you already built with reference_height(model, method=VerticalSort()), sharing its sort across several diagnostics.
• ProfileLookup(b✶, z✶) with plain arrays holds the reference profile fixed in time, which also makes the diagnostic sort-free: each compute! is then a filter plus two binary-search lookups.

eₐ is convex in buoyancy (∂²eₐ/∂b² = ∂z✶/∂b ≥ 0 on a stable profile), so for a filter that averages at fixed z Jensen's inequality gives eₐˢ ≥ 0 pointwise, up to a marginal error where a filtered buoyancy falls between the profile's entries (the lookup then takes the nearest class). A filter with a vertical component mixes heights as well as buoyancies and needs extra steps, not currently implemented, to guarantee a non-negative eₐˢ, so filter has to act only in the horizontal and a filter with a vertical component throws an ArgumentError.

filter is any callable mapping a field to its low-pass-filtered counterpart, e.g. a reusable GaussianFilter or BoxFilter over dims = (1, 2). The filtered buoyancy and the filtered APE are materialized as Fields internally (so the separable filter takes its fast staged path), and the returned object is a lazy operation over them, ready for Field, Integral and OutputWriters. It lives at (Center, Center, Center), 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)filter = GaussianFilter(; dims=(1, 2), σ=0.1)SubFilterAvailablePotentialEnergy(model, filter)# outputSubFilterAvailablePotentialEnergy KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: subfilter_ape_ccc (generic function with 1 method)└── arguments: ("Oceananigans.AbstractOperations.BinaryOperation",)└── computes: subfilter available potential energy  filter(eₐ) - eₐˡ

A convenience method SubFilterAvailablePotentialEnergy(model; σ, dims, boundary, N) builds the Gaussian filter for you along dims, (1, 2) by default, from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ). geopotential_height enters the buoyancy construction exactly as in reference_height.

SubFilterAvailablePotentialEnergyDissipationRate(
    model,
    filter;
    method,
    geopotential_height
)

Return the subfilter-scale (SFS) available potential energy dissipation rate εₐˢ, the APE destruction by diffusion carried by the scales that a low-pass filter removes:

    εₐˢ = filter(εₐ) - εₐˡ ,   εₐˡ = -q̄ᵢ ∂ᵢΥˡ ,   q̄ᵢ = filter(qᵢ) ,   Υˡ = z✶(b̄) - z

where εₐ = -qᵢ∂ᵢΥ is the dissipation rate of the full field (AvailablePotentialEnergyDissipationRate), qᵢ is the closure's own diffusive buoyancy flux, and εₐˡ is the same contraction evaluated on the filtered state (FilteredAvailablePotentialEnergyDissipationRate): the filtered flux q̄ᵢ against the displacement potential Υˡ of the filtered buoyancy b̄ = filter(b). Filtering the flux rather than recomputing it from b̄ is what makes εₐˡ the sink of the filtered-state budget when κ varies in space; that docstring has the details. Both states are measured against one shared reference profile, exactly as in SubFilterAvailablePotentialEnergy, whose budget this is the diffusive sink of; it mirrors what SubFilterKineticEnergyDissipationRate is to the subfilter kinetic energy.

method has to be a ProfileLookup, for the reason SubFilterAvailablePotentialEnergy gives, and the lookup also makes each z✶ a function of buoyancy alone — the property that differentiating Υ and Υˡ needs (see AvailablePotentialEnergyDisplacementPotential). filter has to act only in the horizontal, since the current implementation cannot guarantee a non-negative APE with a filter that has a vertical component; such a filter throws an ArgumentError. Like AvailablePotentialEnergyDissipationRate, this diagnostic needs the buoyancy to be a tracer the closure diffuses (BuoyancyTracer only) and a closure that supplies a diffusive flux.

filter is any callable mapping a field to its low-pass-filtered counterpart, e.g. a reusable GaussianFilter or BoxFilter. The filtered fluxes, the filtered buoyancy, Υˡ and the filtered full-field dissipation are materialized as Fields internally, so the returned object is a lazy operation ready for Field, Integral and OutputWriters, recomputing (re-sorting included) as the simulation evolves. It lives at (Center, Center, Center), 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, closure=ScalarDiffusivity(κ=1e-4))filter = GaussianFilter(; dims=(1, 2), σ=0.1)SubFilterAvailablePotentialEnergyDissipationRate(model, filter)# outputSubFilterAvailablePotentialEnergyDissipationRate KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: subfilter_ape_dissipation_rate_ccc (generic function with 1 method)└── arguments: ("Oceananigans.AbstractOperations.BinaryOperation",)└── computes: subfilter available potential energy dissipation rate  filter(εₐ) - εₐˡ

A convenience method SubFilterAvailablePotentialEnergyDissipationRate(model; σ, dims, boundary, N) builds the Gaussian filter for you along dims, (1, 2) by default, from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ).

SubFilterAvailablePotentialToKineticEnergyConversion(
    model,
    filter;
    method,
    geopotential_height
)

Return the subfilter-scale (SFS) conversion of available potential energy into kinetic energy τˡ(w, bᵣ), the rate at which the scales a low-pass filter removes release their APE to the subfilter flow:

    τˡ(w, bᵣ) = filter(w bᵣ) - w̄ b_rˡ ,   bᵣ = b - b✶(z) ,   b_rˡ = b̄ - b✶(z)

It is the subfilter half of the split whose filtered half is FilteredAvailablePotentialToKineticEnergyConversion w̄b_rˡ: the two sum to filter(w bᵣ), so the subfilter and filtered budgets between them exchange exactly what the full field converts. It enters this budget as -τˡ(w, bᵣ) and the subfilter kinetic energy budget (SubFilterKineticEnergy) as +τˡ(w, bᵣ), so it is a reversible exchange rather than a source or a sink.

The reference profile is not filtered in either half — b_rˡ is b̄ - b✶(z), not filter(bᵣ), which would filter the reference along with the buoyancy. That is what makes the two halves an exact decomposition; FilteredAvailablePotentialToKineticEnergyConversion gives the reason. The two choices coincide for the horizontal filters this diagnostic accepts, since b✶ is a function of z alone.

method has to be a ProfileLookup, for the reason SubFilterAvailablePotentialEnergy gives, and both halves are built on the one profile it supplies. filter has to act only in the horizontal, since the current implementation cannot guarantee a non-negative APE with a filter that has a vertical component; such a filter throws an ArgumentError. filter is any callable mapping a field to its low-pass-filtered counterpart, e.g. a reusable GaussianFilter or BoxFilter over dims = (1, 2). The result lives at (Center, Center, Center), 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)filter = GaussianFilter(; dims=(1, 2), σ=0.1)SubFilterAvailablePotentialToKineticEnergyConversion(model, filter)# outputSubFilterAvailablePotentialToKineticEnergyConversion KernelFunctionOperation at (Center, Center, Center)├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo├── kernel_function: subfilter_ape_to_ke_conversion_ccc (generic function with 1 method)└── arguments: ("Oceananigans.AbstractOperations.BinaryOperation",)└── computes: subfilter APE to KE conversion  filter(wbᵣ) - w̄b_rˡ

A convenience method SubFilterAvailablePotentialToKineticEnergyConversion(model; σ, dims, boundary, N) builds the Gaussian filter for you along dims, (1, 2) by default, from a standard deviation σ (with σ = ℓ / (2√(2 ln 2)) for a FWHM ℓ).