GitHub

Flow diagnostics

The FlowDiagnostics module provides a collection of commonly-used diagnostics for characterizing the state of ocean and turbulent flows. These are not budget terms from a specific governing equation but rather derived quantities that summarize important flow properties at a glance.

The module includes:

  • Stability parameters: the gradient Richardson number ($Ri$) comparing stratification to shear, useful for predicting shear instabilities ($Ri < 0.25$).
  • Rotation diagnostics: the Rossby number ($Ro$), which measures relative vorticity against planetary vorticity and indicates whether flow dynamics are geostrophically balanced or ageostrophic.
  • Potential vorticity: the Ertel PV ($\boldsymbol{\omega}_{\text{tot}} \cdot \nabla b$), a materially conserved quantity under adiabatic, inviscid conditions that is fundamental for understanding large-scale ocean dynamics. A thermal-wind-balance approximation and a directional decomposition are also provided.
  • Velocity gradient tensor invariants: the strain rate tensor modulus ($\|S_{ij}\|$), vorticity tensor modulus ($\|\Omega_{ij}\|$), and the $Q$-criterion for vortex identification.
  • Mixed layer depth: computed by scanning downward from the surface to find where buoyancy or density departs from the surface value by more than a user-specified threshold.
  • Bottom cell value: extracts the value of any diagnostic at the bottommost active cell (respecting immersed boundaries).

Example

julia> using Oceananigans, Oceanostics

julia> using Oceanostics: FlowDiagnostics

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> model = NonhydrostaticModel(grid; buoyancy=BuoyancyTracer(), tracers=:b, coriolis=FPlane(f=1e-4));

julia> Ri = FlowDiagnostics.RichardsonNumber(model)
KernelFunctionOperation at (Center, Center, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: richardson_number_ccf (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Field", "Tuple")

julia> Ro = FlowDiagnostics.RossbyNumber(model)
KernelFunctionOperation at (Face, Face, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: rossby_number_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "NamedTuple")

julia> EPV = FlowDiagnostics.ErtelPotentialVorticity(model)
KernelFunctionOperation at (Face, Face, Face)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: ertel_potential_vorticity_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Field", "Int64", "Int64", "Float64")

julia> S = FlowDiagnostics.StrainRateTensorModulus(model)
KernelFunctionOperation at (Center, Center, Center)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── kernel_function: strain_rate_tensor_modulus_ccc (generic function with 1 method)
└── arguments: ("Field", "Field", "Field")

Richardson number

The gradient Richardson number measures the ratio of stratification to shear:

\[Ri = \frac{\partial b / \partial z}{|\partial \mathbf{u}_h / \partial z|^2}\]

where $z$ is the true vertical (anti-parallel to gravity). Computed at (Center, Center, Face).

Rossby number

The Rossby number measures the ratio of relative vorticity to planetary vorticity in the direction of the rotation axis:

\[Ro = \frac{\omega^z}{f}\]

Computed at (Face, Face, Face).

Oceanostics.FlowDiagnostics.RossbyNumber — Type
RossbyNumber(
    model;
    loc,
    add_background,
    dWdy_bg,
    dVdz_bg,
    dUdz_bg,
    dWdx_bg,
    dUdy_bg,
    dVdx_bg
)

Calculate the Rossby number using the vorticity in the rotation axis direction according to model.coriolis. Rossby number is defined as

    Ro = ωᶻ / f

where ωᶻ is the vorticity in the Coriolis axis of rotation and f is the Coriolis rotation frequency.

source

Potential vorticity

Ertel potential vorticity

The full Ertel potential vorticity:

\[\text{EPV} = \boldsymbol{\omega}_{\text{tot}} \cdot \nabla b\]

where $\boldsymbol{\omega}_{\text{tot}}$ is the absolute (relative + planetary) vorticity. Computed at (Face, Face, Face).

Oceanostics.FlowDiagnostics.ErtelPotentialVorticity — Type
ErtelPotentialVorticity(model; tracer_name, loc)

Calculate the Ertel Potential Vorticty for model, where the characteristics of the Coriolis rotation are taken from model.coriolis. The Ertel Potential Vorticity is defined as

EPV = ωₜₒₜ ⋅ ∇b

where ωₜₒₜ is the total (relative + planetary) vorticity vector, b is the buoyancy and ∇ is the gradient operator.

julia> using Oceananigans

julia> grid = RectilinearGrid(topology = (Flat, Flat, Bounded), size = 4, extent = 1);

julia> N² = 1e-6;

julia> b_bcs = FieldBoundaryConditions(top=GradientBoundaryCondition(N²));

julia> model = NonhydrostaticModel(grid; coriolis=FPlane(f=1e-4), buoyancy=BuoyancyTracer(), tracers=:b, boundary_conditions=(; b=b_bcs));

julia> stratification(z) = N² * z;

julia> set!(model, b=stratification)

julia> using Oceanostics: ErtelPotentialVorticity

julia> EPV = ErtelPotentialVorticity(model)
KernelFunctionOperation at (Face, Face, Face)
├── grid: 1×1×4 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── kernel_function: ertel_potential_vorticity_fff (generic function with 1 method)
└── arguments: ("Field", "Field", "Field", "Field", "Int64", "Int64", "Float64")

julia> interior(Field(EPV))
1×1×5 view(::Array{Float64, 3}, 1:1, 1:1, 4:8) with eltype Float64:
[:, :, 1] =
 0.0

[:, :, 2] =
 1.0000000000000002e-10

[:, :, 3] =
 9.999999999999998e-11

[:, :, 4] =
 1.0000000000000002e-10

[:, :, 5] =
 1.0e-10

Note that EPV values are correctly calculated both in the interior and the boundaries. In the interior and top boundary, EPV = f×N² = 10⁻¹⁰, while EPV = 0 at the bottom boundary since ∂b/∂z is zero there.

source

Thermal wind potential vorticity

A simplified PV assuming thermal wind balance:

\[\text{TWPV} = (f + \omega^z)\,\partial_z b - f\left[(\partial_z U)^2 + (\partial_z V)^2\right]\]

Oceanostics.FlowDiagnostics.ThermalWindPotentialVorticity — Type
ThermalWindPotentialVorticity(model; tracer_name, loc)

Calculate the Potential Vorticty assuming thermal wind balance for model, where the characteristics of the Coriolis rotation are taken from model.coriolis. The Potential Vorticity in this case is defined as

    TWPV = (f + ωᶻ) ∂b/∂z - f ((∂U/∂z)² + (∂V/∂z)²)

where f is the Coriolis frequency, ωᶻ is the relative vorticity in the z direction, b is the buoyancy, and ∂U/∂z and ∂V/∂z comprise the thermal wind shear.

source

Directional Ertel potential vorticity

The contribution to the Ertel PV from a single user-specified direction.

Oceanostics.FlowDiagnostics.DirectionalErtelPotentialVorticity — Type
DirectionalErtelPotentialVorticity(
    model,
    direction;
    tracer_name,
    loc
)

Calculate the contribution from a given direction to the Ertel Potential Vorticity basde on a model and a direction. The Ertel Potential Vorticity is defined as

EPV = ωₜₒₜ ⋅ ∇b

where ωₜₒₜ is the total (relative + planetary) vorticity vector, b is the buoyancy and ∇ is the gradient operator.

source

Velocity gradient tensor invariants

Strain rate tensor modulus

\[\|S_{ij}\| = \sqrt{S_{ij} S_{ij}}, \qquad S_{ij} = \tfrac{1}{2}(\partial_j u_i + \partial_i u_j)\]

Oceanostics.FlowDiagnostics.StrainRateTensorModulus — Type
StrainRateTensorModulus(model; loc)

Calculate the modulus (absolute value) of the strain rate tensor S, which is defined as the symmetric part of the velocity gradient tensor:

    Sᵢⱼ = ½(∂ⱼuᵢ + ∂ᵢuⱼ)

Its modulus is then defined (using Einstein summation notation) as

    || Sᵢⱼ || = √(Sᵢⱼ Sᵢⱼ)
source

Vorticity tensor modulus

\[\|\Omega_{ij}\| = \sqrt{\Omega_{ij} \Omega_{ij}}, \qquad \Omega_{ij} = \tfrac{1}{2}(\partial_j u_i - \partial_i u_j)\]

Oceanostics.FlowDiagnostics.VorticityTensorModulus — Type
VorticityTensorModulus(model; loc)

Calculate the modulus (absolute value) of the vorticity tensor Ω, which is defined as the antisymmetric part of the velocity gradient tensor:

    Ωᵢⱼ = ½(∂ⱼuᵢ - ∂ᵢuⱼ)

Its modulus is then defined (using Einstein summation notation) as

    || Ωᵢⱼ || = √(Ωᵢⱼ Ωᵢⱼ)
source

Q-criterion

\[Q = \tfrac{1}{2}(\Omega_{ij}\Omega_{ij} - S_{ij}S_{ij})\]

Used to identify vortices in fluid flow; positive values indicate vortex-dominated regions.

Oceanostics.FlowDiagnostics.QVelocityGradientTensorInvariant — Type
QVelocityGradientTensorInvariant(model; loc)

Calculate the value of the Q velocity gradient tensor invariant. This is usually just called Q and it is generally used for identifying and visualizing vortices in fluid flow.

The definition and nomenclature comes from the equation for the eigenvalues λ of the velocity gradient tensor ∂ⱼuᵢ:

    λ³ + P λ² + Q λ + T = 0

from where Q is defined as

    Q = ½ (ΩᵢⱼΩᵢⱼ - SᵢⱼSᵢⱼ)

and where Sᵢⱼ= ½(∂ⱼuᵢ + ∂ᵢuⱼ) and Ωᵢⱼ= ½(∂ⱼuᵢ - ∂ᵢuⱼ). More info about it can be found in doi:10.1063/1.5124245.

source

Mixed layer depth

Oceanostics.FlowDiagnostics.MixedLayerDepth — Function
MixedLayerDepth(grid, args; criterion)

Returns the mixed layer depth defined as the depth at which criterion is true.

Defaults to DensityAnomalyCriterion where the depth is that at which the density is some threshold (defaults to 0.125kg/m³) higher than the surface density.

When DensityAnomalyCriterion is used, the arguments buoyancy_formulation and C should be supplied where buoyancy_formulation should be the buoyancy model, and C should be a named tuple of (; T, S), (; T) or (; S) (the latter two if the buoyancy model specifies a constant salinity or temperature).

source
Oceanostics.FlowDiagnostics.BuoyancyAnomalyCriterion — Type
struct BuoyancyAnomalyCriterion{FT} <: Oceanostics.FlowDiagnostics.AbstractAnomalyCriterion

Defines the mixed layer to be the depth at which the buoyancy is more than threshold greater than the surface buoyancy (but the pertubaton is usually negative).

When this model is used, the arguments buoyancy_formulation and C should be supplied where C should be the named tuple (; b), with b the buoyancy tracer.

source
Oceanostics.FlowDiagnostics.DensityAnomalyCriterion — Type
struct DensityAnomalyCriterion{FT} <: Oceanostics.FlowDiagnostics.AbstractAnomalyCriterion

Defines the mixed layer to be the depth at which the density is more than threshold greater than the surface density.

When this model is used, the arguments buoyancy_formulation and C should be supplied where buoyancy_formulation should be the buoyancy model, and C should be a named tuple of (; T, S), (; T) or (; S) (the latter two if the buoyancy model specifies a constant salinity or temperature).

source

Bottom cell value