Potential energy equation

The PotentialEnergyEquation module provides a diagnostic for the specific gravitational potential energy (per unit mass). In a Boussinesq fluid, the specific potential energy is defined as

\[e_p = -bz = \frac{g\rho}{\rho_0} z\]

where $b = -g\rho/\rho_0$ is buoyancy, $z$ is the vertical coordinate, $g$ is gravitational acceleration, $\rho$ is density, and $\rho_0$ is a reference density. The quantity $e_p$ has units of m² s⁻² (energy per unit mass).

Lower case for densities, upper case for their integrals

Throughout Oceanostics a lower-case $e$ is an energy density, the pointwise quantity a diagnostic returns, and the matching upper-case $E$ is its volume integral,

\[E_p = \int e_p \, \mathrm{d}V ,\]

and likewise for the kinetic, background, and available potential energies.

The potential energy equation

Throughout what follows we assume that gravity_unit_vector points towards the NegativeZDirection, as is the default in Oceananigans, but tilted domains where this is not true are possible. Since $e_p$ is just $-z$ times the buoyancy, and $z$ does not change in time, the equation for $e_p$ follows directly from the tracer equation applied to $b$. In the convention Oceananigans uses (and ignoring background fields), that equation reads

\[\partial_t b = -\partial_j (u_j b) - \partial_j q_j + F_b ,\]

where $q_j$ is the diffusive flux of buoyancy supplied by the closure ($q_j = -\kappa\,\partial_j b$ for Fickian diffusion with diffusivity $\kappa$) and $F_b$ is any applied forcing. Multiplying through by $-z$,

\[\partial_t e_p = -z\,\partial_t b = z\,\partial_j(u_j b) + z\,\partial_j q_j - z F_b .\]

Pulling $z$ inside the derivatives separates them, and since $\partial_j z$ is non-zero only in the vertical, the product rule gives

\[\partial_t e_p = \underbrace{-\partial_j(u_j e_p)}_{\text{advection}} \underbrace{- wb}_{\text{PE to KE conversion}} + \underbrace{\partial_j(z q_j)}_{\text{diffusive transport}} \underbrace{- q_3}_{\text{diffusive buoyancy flux}} \underbrace{- z F_b}_{\text{forcing}}.\]

The two terms written as divergences transport $e_p$ and vanish when integrated over a periodic or closed domain (with impermeable, insulating walls).

Terms and what is implemented

Every term above is implemented, two of them elsewhere in Oceanostics.

Quantity	Expression	Diagnostic
Potential energy	$e_p = -bz$	`PotentialEnergy`
Tendency	$\partial_t e_p = -z\,\partial_t b$	`PotentialEnergyTendency`
Buoyancy advection	$z\,\partial_j(u_j b) = -\partial_j(u_j e_p) - wb$	`PotentialEnergyBuoyancyAdvection`
Advection	$\partial_j(u_j e_p)$	`PotentialEnergyAdvection`
PE to KE conversion	$wb$	`PotentialToKineticEnergyConversion`
Buoyancy diffusion	$z\,\partial_j q_j = \partial_j(z q_j) - q_3$	`PotentialEnergyBuoyancyDiffusion`
Diffusive transport	$\partial_j(z q_j)$	`PotentialEnergyDiffusion`
Diffusive vertical buoyancy flux	$\Phi = -q_3$	`DiffusiveVerticalBuoyancyFlux`
Forcing	$-z F_b$	`PotentialEnergyForcing`

Each of the two flux terms appears twice, once in each of the forms the derivation above relates:

\[\underbrace{z\,\partial_j(u_j b)}_{\texttt{BuoyancyAdvection}} = -\underbrace{\partial_j(u_j e_p)}_{\texttt{Advection}} - wb , \qquad \underbrace{z\,\partial_j q_j}_{\texttt{BuoyancyDiffusion}} = \underbrace{\partial_j(z q_j)}_{\texttt{Diffusion}} + \Phi .\]

The Buoyancy* pair are the $-z\,\times$ forms, which is the convention the kinetic energy equation follows as well ($u_i\partial_j(u_ju_i)$ rather than $\partial_j(u_jK)$). Taken that way they are the model's own buoyancy tendency split apart, so

\[\texttt{Tendency} = \texttt{BuoyancyAdvection} + \texttt{BuoyancyDiffusion} + \texttt{Forcing}\]

holds cell by cell rather than to within a truncation error.

Note that the KE to PE conversion term can go by three names since it is shared by different budgets. It stays defined in the kinetic energy equation, where PotentialEnergyConversion names what it does to the kinetic energy. using Oceanostics brings in PotentialToKineticEnergyConversion, which names the exchange and says which way it runs. This module and the available potential energy equation both re-export that and additionally give it KineticEnergyConversion, short enough to read well beside the other terms of a potential energy budget. That last alias is scoped to those two modules: using Oceanostics does not bring it in, because unprefixed it says nothing about which budget it belongs to.

$e_p$ also splits into two parts that do have their own modules. The background potential energy $e_b$ is the portion the flow could never release, and the available potential energy $e_a = e_p - e_b$ is the remainder.

Buoyancy formulations

PotentialEnergy is implemented for three buoyancy model types:

BuoyancyTracer: uses the buoyancy field $b$ directly as $e_p = -bz$.
SeawaterBuoyancy with LinearEquationOfState: computes buoyancy from a linear equation of state applied to temperature and/or salinity tracers.
SeawaterBuoyancy with BoussinesqEquationOfState (from SeawaterPolynomials.jl): computes density from a nonlinear equation of state. An optional geopotential_height keyword argument allows using a potential density referenced to a fixed depth instead of in-situ density.

For now the full budget requires gravity to be aligned with the negative $z$-direction (NegativeZDirection), and every term checks it at construction. Tilting gravity does not make these terms approximate, it makes them a different quantity: the height that works against gravity becomes $z\cos\theta - y\sin\theta$, so $-bz$ is wrong both by a factor and by a cross-slope term, and nothing in the output would reveal it.

Two diagnostics are exempt because they do not depend on that alignment. DiffusiveVerticalBuoyancyFlux never touches $z$; it returns the vertical component of the closure's diffusive flux, which is what it promises whatever gravity does. PotentialToKineticEnergyConversion is the full $u_i b_i$ contraction over the gravity-projected buoyancy components, so it is correct under any tilt and only reads as $wb$ when $\hat{g} = -\hat{z}$. Neither is the term the budget above needs once gravity tilts, so the split still does not close there — only the individual quantities remain meaningful.

Summary of $e_p$ equation terms

Oceanostics.PotentialEnergyEquation.PotentialEnergy — Type

PotentialEnergy(model; location, geopotential_height)

Return a KernelFunctionOperation to compute the PotentialEnergy per unit volume,

\[eₚ = \frac{gρ}{ρ₀}z = -bz\]

at each grid location in model. PotentialEnergy is defined for both BuoyancyTracer and SeawaterBuoyancy. See the relevant Oceananigans.jl documentation on buoyancy models for more information about available options.

The optional keyword argument geopotential_height is only used if one wishes to calculate eₚ with a potential density referenced to geopotential_height, rather than in-situ density, when using a BoussinesqEquationOfState.

Example

Usage with a BuoyancyTracer buoyacny model

julia> using Oceananigans, Oceanostics

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded))
1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-1000.0, 0.0] regularly spaced with Δz=10.0

julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=(:b,))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergyEquation.PotentialEnergy(model)
PotentialEnergy KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: minus_bz_ccc (generic function with 3 methods)
└── arguments: ("Field",)
└── computes: potential energy per unit volume  eₚ = -bz

The default behaviour of PotentialEnergy uses the in-situ density in the calculation when the equation of state is a BoussinesqEquationOfState:

julia> using Oceananigans, SeawaterPolynomials.TEOS10, Oceanostics

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded))
1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-1000.0, 0.0] regularly spaced with Δz=10.0

julia> tracers = (:T, :S)
(:T, :S)

julia> eos = TEOS10EquationOfState()
BoussinesqEquationOfState{Float64}:
├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64}
└── reference_density: 1020.0

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: BoussinesqEquationOfState{Float64}

julia> model = NonhydrostaticModel(grid; buoyancy, tracers)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergyEquation.PotentialEnergy(model)
PotentialEnergy KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: minus_bz_ccc (generic function with 3 methods)
└── arguments: ("KernelFunctionOperation", "NamedTuple")
└── computes: potential energy per unit volume  eₚ = -bz

To use a reference density set a constant value for the keyword argument geopotential_height and pass this the function. For example,

julia> using Oceananigans, SeawaterPolynomials.TEOS10, Oceanostics

julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded));

julia> tracers = (:T, :S);

julia> eos = TEOS10EquationOfState();

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos);

julia> model = NonhydrostaticModel(grid; buoyancy, tracers);

julia> geopotential_height = 0; # density variable will be σ₀

julia> PotentialEnergyEquation.PotentialEnergy(model)
PotentialEnergy KernelFunctionOperation at (Center, Center, Center)
├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: minus_bz_ccc (generic function with 3 methods)
└── arguments: ("KernelFunctionOperation", "NamedTuple")
└── computes: potential energy per unit volume  eₚ = -bz

source

Oceanostics.PotentialEnergyEquation.PotentialEnergyTendency — Type

PotentialEnergyTendency(model; location)

Return a KernelFunctionOperation computing the tendency of the potential energy eₚ = -bz,

    ∂ₜeₚ = -z ∂ₜb ,

where ∂ₜb is Oceananigans' own tracer tendency for the buoyancy. Since that kernel is the one the model steps, this is the whole right-hand side of the eₚ equation in one term: advection by the total (perturbation plus background) flow, advection of a background buoyancy field, diffusion, and forcing. The individual terms are PotentialEnergyBuoyancyAdvection, PotentialEnergyBuoyancyDiffusion and PotentialEnergyForcing, and they sum to this one cell by cell — but only when the model has none of the features listed below. Note it is those two Buoyancy* terms that sum, not PotentialEnergyAdvection and PotentialEnergyDiffusion — the latter are the transports alone, and differ from them by the two conversions wb and Φ.

Because this term is tracer_tendency itself, it carries everything the model steps, including four things the other three diagnostics do not. With any of them present the split falls short by exactly that term, silently:

a BackgroundField buoyancy B, contributing -z ∂ⱼ(uⱼB), which has no diagnostic yet;
an AdvectiveForcing or a biogeochemical drift velocity, which tracer_tendency folds into the advecting flow but PotentialEnergyBuoyancyAdvection does not pick up by default;
an ImmersedBoundaryGrid, whose immersed_∇_dot_qᶜ no diagnostic here mirrors;
a biogeochemical transition term for b.

Defined for BuoyancyTracer models, where b is one of the model's tracers.

using Oceananigans, Oceanostics

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(κ=1e-4))

PotentialEnergyTendency(model)   # or `Tendency` inside the module

# output

PotentialEnergyTendency KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: minus_z_∂ₜb_ccc (generic function with 1 method)
└── arguments: ("Val", "Val", "Centered", "ScalarDiffusivity", "Nothing", "BuoyancyForce", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Clock", "Returns")
└── computes: potential energy tendency  ∂ₜeₚ = -z ∂ₜb

source

Oceanostics.PotentialEnergyEquation.PotentialEnergyAdvection — Type

PotentialEnergyAdvection(model; velocities, location)

Return a KernelFunctionOperation computing the advection of the potential energy itself,

    ADV = ∂ⱼ(uⱼeₚ) ,

with eₚ = -bz handed to the model's own advection scheme the way a tracer would be. uⱼ defaults to the total velocity, perturbation plus background, which is what the model advects with; pass velocities to override it. It enters the eₚ equation with a minus sign, as TracerEquation.Advection does for a tracer.

This is a transport and nothing else: it is a flux divergence, so over a periodic or closed domain it integrates to zero to roundoff. What it does not include is the conversion wb that arises from weighting the buoyancy equation by -z; the two together make PotentialEnergyBuoyancyAdvection, which is the term that sums with the others to PotentialEnergyTendency:

    z ∂ⱼ(uⱼb) = -∂ⱼ(uⱼeₚ) - wb .

using Oceananigans, Oceanostics

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)

PotentialEnergyAdvection(model)   # or `Advection` inside the module

# output

PotentialEnergyAdvection KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: div_U_eₚ_ccc (generic function with 1 method)
└── arguments: ("Centered", "NamedTuple", "KernelFunctionOperation")
└── computes: potential energy advection  ∂ⱼ(uⱼeₚ)

source

Oceanostics.PotentialEnergyEquation.PotentialEnergyBuoyancyAdvection — Type

PotentialEnergyBuoyancyAdvection(
    model;
    velocities,
    location
)

Return a KernelFunctionOperation computing the advective term of the eₚ = -bz equation,

    ADV = z ∂ⱼ(uⱼb) ,

the advection of the buoyancy weighted by -z. uⱼ defaults to the total velocity, perturbation plus background, which is what the model advects with; pass velocities to override it.

Pulling z inside the derivative writes this as a transport of eₚ plus the conversion term, ADV = -∂ⱼ(uⱼeₚ) - wb, so over a periodic or closed domain its volume integral is -∫wb dV, the (negated) PotentialToKineticEnergyConversion.

A background buoyancy field is advected by a term of its own, which this one does not include and which has no diagnostic here yet.

The default velocities is the perturbation plus background flow. That is what the model advects with in the ordinary case, but not when an AdvectiveForcing or a biogeochemical drift velocity is in play: tracer_tendency folds those in too, through with_advective_forcing and biogeochemical_drift_velocity. Pass velocities explicitly to match the model in that case.

using Oceananigans, Oceanostics

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b)

PotentialEnergyBuoyancyAdvection(model)   # or `BuoyancyAdvection` inside the module

# output

PotentialEnergyBuoyancyAdvection KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: z_div_Uc_ccc (generic function with 1 method)
└── arguments: ("Centered", "NamedTuple", "Field")
└── computes: potential energy buoyancy advection  z ∂ⱼ(uⱼb)

source

Oceanostics.PotentialEnergyEquation.PotentialEnergyDiffusion — Type

PotentialEnergyDiffusion(model; location)

Return a KernelFunctionOperation computing the diffusive transport of the potential energy,

    DIFF = ∂ⱼ(z qⱼ) ,

where qⱼ is the closure's own diffusive flux of buoyancy (qⱼ = -κ ∂ⱼb for Fickian diffusion). It is the diffusive counterpart of PotentialEnergyAdvection: a flux divergence, so over a periodic or closed domain it integrates to zero to roundoff.

What it does not include is the vertical flux -q₃ that arises from weighting the buoyancy equation by -z; the two together make PotentialEnergyBuoyancyDiffusion, which is the term that sums with the others to PotentialEnergyTendency:

    z ∂ⱼqⱼ = ∂ⱼ(z qⱼ) - q₃ = ∂ⱼ(z qⱼ) + Φ ,

with Φ the DiffusiveVerticalBuoyancyFlux.

Like the other diffusive terms this reads κ∇b off the closure, so it needs the buoyancy to be a tracer the closure diffuses.

Two limits are worth knowing. The divergence is assembled from the unconditional diffusive_flux_*, while PotentialEnergyBuoyancyDiffusion reaches the closure through Oceananigans' ∇_dot_qᶜ, which uses the conditional _diffusive_flux_* that ImmersedBoundaries overrides to zero across immersed faces. Off an immersed grid the two are the same flux; on one this term includes fluxes its partner zeroes, and neither the telescoping nor the identity above survives. And the horizontal fluxes are weighted by z at the cell centre rather than at the x- and y-faces they live on, which is exact wherever znode does not vary with i and j — every RectilinearGrid and LatitudeLongitudeGrid with an immutable vertical coordinate — but not under a MutableVerticalDiscretization.

using Oceananigans, Oceanostics

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(κ=1e-4))

PotentialEnergyDiffusion(model)   # or `Diffusion` inside the module

# output

PotentialEnergyDiffusion KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: div_z_q_ccc (generic function with 1 method)
└── arguments: ("ScalarDiffusivity", "Nothing", "Val", "Field", "Clock", "NamedTuple", "BuoyancyForce")
└── computes: potential energy diffusive transport  ∂ⱼ(z qⱼ)

source

Oceanostics.PotentialEnergyEquation.PotentialEnergyBuoyancyDiffusion — Type

PotentialEnergyBuoyancyDiffusion(model; location)

Return a KernelFunctionOperation computing the diffusive term of the eₚ = -bz equation,

    DIFF = z ∂ⱼqⱼ ,

where qⱼ is the closure's own diffusive flux of buoyancy (qⱼ = -κ ∂ⱼb for Fickian diffusion).

Pulling z inside the derivative writes this as a transport plus the vertical flux, DIFF = ∂ⱼ(zqⱼ) - q₃, so over a periodic or closed domain its volume integral is ∫Φ dV, the DiffusiveVerticalBuoyancyFlux. Diffusion in the horizontal drops out of that integral entirely, since z does not vary along it.

Like Φ, this reads κ∇b off the closure, so it needs the buoyancy to be a tracer the closure diffuses.

using Oceananigans, Oceanostics

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(κ=1e-4))

PotentialEnergyBuoyancyDiffusion(model)   # or `BuoyancyDiffusion` inside the module

# output

PotentialEnergyBuoyancyDiffusion KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: z_∇_dot_qᶜ_ccc (generic function with 1 method)
└── arguments: ("ScalarDiffusivity", "Nothing", "Val", "Field", "Clock", "NamedTuple", "BuoyancyForce")
└── computes: potential energy buoyancy diffusion  z ∂ⱼqⱼ

source

Oceanostics.PotentialEnergyEquation.DiffusiveVerticalBuoyancyFlux — Type

DiffusiveVerticalBuoyancyFlux(model; location)

Return a KernelFunctionOperation computing the vertical diffusive buoyancy flux, taken upward,

    Φ = κ ∂b/∂z = -q₃ ,

which is the work diffusion does against gravity as it smooths the stratification. It is the diffusive conversion term of the e_p equation, and the only way diffusion can change the potential energy of a closed domain: the rest of the diffusive contribution is a flux divergence that integrates to zero. The result lives at (Center, Center, Center), per unit mass (units m² s⁻³).

Φ is also the second of the two parts AvailablePotentialEnergyDissipationRate is written out of,

    ε_A = κ (∂z✶/∂b) |∇b|² - Φ ,

and it is the part that carries no available energy with it: Φ enters the E_p and E_b budgets identically, so it cancels from their difference, which is E_a. A statically stable, horizontally uniform stratification is its own reference state and has ε_A = 0 cell by cell, so there Φ accounts for the whole of the diapycnal mixing rate. Adding it back to ε_A recovers that rate in general, and its volume integral is the rate BackgroundPotentialEnergy grows:

    d/dt ∫e_b dV = ∫(ε_A + Φ) dV .

For a constant κ, Φ volume-integrates to a boundary term, since no buoyancy crosses the top or the bottom of a closed domain:

    ∫Φ dV = κ A [b(z_top) - b(z_bottom)]        (constant κ only)

with A the domain's horizontal area, so the flow enters only through the buoyancy difference across it. That collapse is exact rather than approximate: the discrete sum telescopes, which is also what makes the cells against a wall report half of what the interior does, since a no-flux wall zeroes the flux on the outer face and the cell center averages the two faces bounding it.

A κ that varies in space breaks the telescoping, and the boundary form is then simply wrong. What the integral always is, is the flux summed over the interior z faces,

    ∫Φ dV = A Σ κ [b(above) - b(below)] ,

which collapses only when κ factors out of the sum. Under a depth-dependent diffusivity or any closure that computes κ from the flow, the integral depends on the interior arrangement and can come out with the opposite sign to κ A [b(z_top) - b(z_bottom)], so do not reach for the boundary form as a check there. Φ is returned as a field rather than as either number, since it is the pointwise partner of ε_A.

κ ∂b/∂z 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, including ones that compute κ from the flow. That needs the buoyancy to be a tracer the closure diffuses, so this is defined for BuoyancyTracer models only.

Not to be confused with TracerEquation.ZDiffusiveFlux(model, :b), which is the same closure call read raw: that one is q₃ itself, down-gradient and on the z face. This is -q₃, interpolated to the cell center, which are the sign and the location the energy budgets want. Reach for the tracer diagnostic when you want the flux as a flux, and this one when you want it as a term of the e_p equation.

using Oceananigans, Oceanostics

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=(Periodic, Periodic, Bounded))
model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, closure=ScalarDiffusivity(κ=1e-4))

PotentialEnergyDiffusiveVerticalBuoyancyFlux(model)   # or `DiffusiveVerticalBuoyancyFlux` inside the module

# output

PotentialEnergyDiffusiveVerticalBuoyancyFlux KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: diffusive_buoyancy_flux_ccc (generic function with 1 method)
└── arguments: ("ScalarDiffusivity", "Nothing", "Val", "Field", "Clock", "NamedTuple", "BuoyancyForce")
└── computes: diffusive vertical buoyancy flux  Φ = κ ∂b/∂z = -q₃

source

Oceanostics.PotentialEnergyEquation.PotentialEnergyForcing — Type

PotentialEnergyForcing(model; location)

Return a KernelFunctionOperation computing the forcing term of the eₚ = -bz equation,

    FORC = -z Fᵇ ,

where Fᵇ is whatever forcing is applied to the buoyancy. Unlike the transport terms this does not drop out of a volume integral, so a forced run needs it in the budget.

This calls the forcing pointwise, as Fᵇ(i, j, k, grid, clock, fields). An AdvectiveForcing returns zero there by construction — it acts through the advecting velocity instead — so this term reports nothing for one, and PotentialEnergyBuoyancyAdvection will not see it either unless it is handed matching velocities.

using Oceananigans, Oceanostics

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b,
                            forcing = (; b = Forcing((x, y, z, t) -> 1e-8)))

PotentialEnergyForcing(model)

# output

PotentialEnergyForcing KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: minus_z_Fᵇ_ccc (generic function with 1 method)
└── arguments: ("Oceananigans.Forcings.ContinuousForcing", "Clock", "NamedTuple")
└── computes: potential energy forcing  -z Fᵇ

source