Potential energy equation

The PotentialEnergyEquation module provides a diagnostic for computing 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). Potential energy is a key quantity in ocean energetics: its conversion to/from kinetic energy (via the buoyancy production term $wb$) drives ocean circulation and mixing.

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.

The diagnostic requires gravity to be aligned with the negative $z$-direction (NegativeZDirection).

Example

julia> using Oceananigans, Oceanostics

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

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

julia> Ep = PotentialEnergyEquation.PotentialEnergy(model)
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",)

Potential energy

Oceanostics.PotentialEnergyEquation.PotentialEnergyType
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 ones 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 ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergyEquation.PotentialEnergy(model)
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",)

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 ĝ = NegativeZDirection()
└── coriolis: Nothing

julia> PotentialEnergyEquation.PotentialEnergy(model)
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")

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)
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")
source