API Reference

Diffusivity(case, t, var)

Compute the effective numerical diffusivity from variance dissipation diagnostics.

Diffusivity(case, t, var)

Arguments

• case: Case dictionary containing dissipation and gradient fields
• t: Time index
• var: Variable symbol (e.g., :b for buoyancy)

Returns

• KernelFunctionOperation representing the numerical diffusivity field

The diffusivity is computed as κ = -0.5 * (Ax + Az) / Gz, where Ax and Az are the advective dissipation rates and Gz is the vertical gradient squared. This provides a local measure of numerical mixing, as described in the paper.

KineticEnergy(case, i)

Compute the volume-averaged kinetic energy at time index i.

KineticEnergy(case, i)

Arguments

• case: Case dictionary containing velocity and volume fields
• i: Time index

Returns

• Volume-averaged kinetic energy [m²/s²]: (u² + v² + w²) / V

Kinetic energy is computed using volume-weighted velocity components at their respective grid locations (u at Face-Center-Center, v at Center-Face-Center, w at Center-Center-Face).

MeanKineticEnergy(case, i)

Compute the volume-averaged mean kinetic energy (zonally averaged) at time index i.

MeanKineticEnergy(case, i)

Arguments

• case: Case dictionary containing velocity and volume fields
• i: Time index

Returns

• Volume-averaged mean kinetic energy [m²/s²]

Mean kinetic energy is computed from zonally averaged velocity components, representing the energy in the mean flow rather than eddy kinetic energy.

buoyancy_flux(i, j, grid, clock, model_fields, p)

Compute the surface buoyancy flux forcing.

buoyancy_flux(i, j, grid, clock, model_fields, p)

Arguments

• i, j: Grid indices
• grid: Grid object
• clock: Simulation clock
• model_fields: Model fields
• p: Parameters named tuple containing Qᵇ (flux magnitude), y_shutoff (shutoff location), and Ly (domain length)

Returns

• Surface buoyancy flux [m²/s³] following a cosine profile, zero beyond y_shutoff

The flux profile is Qᵇ * cos(3π * y / Ly) for y < y_shutoff, providing differential heating/cooling that drives meridional circulation.

buoyancy_relaxation(i, j, k, grid, clock, model_fields, p)

Compute the buoyancy relaxation forcing in the sponge layer.

buoyancy_relaxation(i, j, k, grid, clock, model_fields, p)

Arguments

• i, j, k: Grid indices
• grid: Grid object
• clock: Simulation clock
• model_fields: Model fields
• p: Parameters named tuple containing relaxation parameters

Returns

• Buoyancy relaxation rate [m/s³] that restores buoyancy toward the initial profile in the sponge layer

This forcing relaxes buoyancy toward the initial exponential profile in the northern sponge region, with strength controlled by mask(y, p) / λt where λt is the relaxation timescale.

calculate_z★_diagnostics(b::Field, vol)

Calculate the re-sorted height field z★ from buoyancy and volume fields.

calculate_z★_diagnostics(b, vol)

Arguments

• b: Buoyancy field
• vol: Volume field

Returns

• z★ field representing the height each fluid parcel would occupy after adiabatic re-sorting

The re-sorted height is computed by sorting the buoyancy field and computing cumulative volume distribution, providing a reference state for potential energy calculations.

channel_simulation(; momentum_advection, tracer_advection, closure, zstar, restart_file, arch, bottom_height, timestepper, grid, initial_file, testcase)

Set up and run the idealized re-entrant channel flow simulation.

channel_simulation(
;
    momentum_advection,
    tracer_advection,
    closure,
    zstar,
    restart_file,
    arch,
    bottom_height,
    timestepper,
    grid,
    initial_file,
    testcase
)

Keyword Arguments

• momentum_advection: Momentum advection scheme (default: WENOVectorInvariant())
• tracer_advection: Tracer advection scheme (default: 7th-order WENO)
• closure: Turbulence closure (default: CATKEVerticalDiffusivity)
• zstar: Whether to use z-star vertical coordinates (default: true)
• restart_file: Optional restart file path (default: nothing)
• arch: Architecture to run on (default: CPU())
• bottom_height: Optional bottom height function (default: nothing for flat bottom)
• timestepper: Timestepper symbol (default: :SplitRungeKutta)
• grid: Grid to use (default: default_grid(...))
• initial_file: Optional initial condition file (default: "tIni_80y_90L.bin")
• testcase: Test case identifier string (default: "0")

Returns

• Simulation object after running to completion

This function sets up the idealized re-entrant channel test case described in the paper. The configuration consists of a 1000 km × 2000 km × 3 km periodic channel forced by sinusoidal wind stress and variable surface heat flux. Buoyancy is restored at the northern boundary to maintain an equilibrated state.

The simulation includes a spin-up phase followed by a long integration (40 years) with time-averaged outputs. This test case demonstrates how numerical mixing interacts with explicit physical mixing in an equilibrated configuration, as shown in the paper.

compute_rpe_density(case::Dict)

Compute reference potential energy (RPE) and available potential energy (APE) time series from a case dictionary.

compute_rpe_density(case)

Arguments

• case: Dictionary containing :b (buoyancy FieldTimeSeries) and :VCCC (volume FieldTimeSeries)

Returns

• Named tuple with rpe and ape arrays containing volume-averaged RPE and APE at each time step

RPE represents the minimum potential energy achievable by adiabatic re-sorting of the density field, while APE is the difference between actual and reference potential energy. These diagnostics are used to quantify irreversible mixing, as described in the paper.

compute_rpe_density(b::Field, vol)

Compute RPE and APE density fields from buoyancy and volume fields.

compute_rpe_density(b, vol)

Arguments

• b: Buoyancy field
• vol: Volume field

Returns

• Named tuple with ze (re-sorted height), εe (RPE density), and αe (APE density) fields

This function computes the re-sorted height z★ by sorting buoyancy and computing the cumulative volume distribution. RPE density is ρ * z★ and APE density is ρ * (z - z★), where ρ is computed from buoyancy using the linear equation of state.

compute_rpe_density_two(case::Dict)

Compute RPE and APE time series using buoyancy directly (without volume weighting).

compute_rpe_density_two(case)

Arguments

• case: Dictionary containing :b (buoyancy FieldTimeSeries) and :VCCC (volume FieldTimeSeries)

Returns

• Named tuple with rpe and ape arrays

This variant computes RPE/APE using the buoyancy field directly, suitable for cases where buoyancy is already volume-weighted or when using different volume conventions.

compute_tracer_dissipation!(sim)

Compute the time rate of change of tracer variance (dissipation) for tracers b and c.

compute_tracer_dissipation!(sim)

Arguments

• sim: The Simulation object containing the model with tracers b and c

Returns

• nothing (modifies sim.model.auxiliary_fields in place)

This function computes the dissipation rate of tracer variance, defined as Δtc² = (c² - c⁻²) / Δt where c⁻ is the tracer value from the previous time step. The results are stored in sim.model.auxiliary_fields.Δtc² and sim.model.auxiliary_fields.Δtb².

This diagnostic is used to quantify numerical mixing introduced by the time discretization, as described in the paper's appendix.

default_closure()

Return the default turbulence closure for the channel simulation.

default_closure()

Returns

• CATKEVerticalDiffusivity closure with default parameters

This closure uses the CATKE (Convective Adjustment and TKE) parameterization for vertical mixing, which is appropriate for equilibrated channel flows with mesoscale forcing.

default_grid(arch, zstar, bottom_height)

Construct the default grid for the channel simulation.

default_grid(arch, zstar, bottom_height)

Arguments

• arch: Architecture to use (CPU() or GPU())
• zstar: Whether to use z-star vertical coordinates
• bottom_height: Optional bottom height function (x, y) -> z, or nothing for flat bottom

Returns

• RectilinearGrid or ImmersedBoundaryGrid with 200×400×90 grid points

The grid spans 1000 km zonally (periodic), 2000 km meridionally (bounded), and uses 90 non-uniform vertical layers with finer spacing near the surface and bottom.

hasclosure(closure, ClosureType)

Check if a closure or tuple of closures contains a specific closure type.

hasclosure(closure, ClosureType)

Arguments

• closure: A single closure or tuple of closures
• ClosureType: The closure type to check for

Returns

• true if closure is of type ClosureType or contains it in a tuple, false otherwise

This helper function is used to determine if CATKE or other specific closures are present in the closure configuration.

idealized_coast(timestepper::Symbol; arch, forced, lowres, free_surface)

Set up and run the idealized coastal baroclinic adjustment test case simulation.

idealized_coast(
    timestepper;
    arch,
    forced,
    lowres,
    free_surface
)

Arguments

• timestepper: Symbol indicating the timestepper (:QuasiAdamsBashforth2 or :SplitRungeKutta3)

Keyword Arguments

• arch: Architecture to run on (default: CPU())
• forced: Whether to apply wind forcing (default: true)
• lowres: Whether to use low resolution (default: false)
• free_surface: Free surface formulation (default: SplitExplicitFreeSurface with CFL=0.7)

Returns

• Simulation object configured but not yet run

This function sets up the idealized coastal baroclinic adjustment test case described in the paper. The configuration features a rectangular domain with a linearly sloping bathymetry and initial meridional salinity gradient that drives baroclinic instabilities. The test case demonstrates how numerical mixing can suppress submesoscale variability, as shown in the paper.

The simulation runs for 40 days and outputs velocity, temperature, salinity, buoyancy, and variance dissipation diagnostics.

idealized_coast_timestep(::Val{timestepper})

Return the recommended time step for the idealized coast test case given a timestepper.

idealized_coast_timestep(_)

Arguments

• timestepper: Symbol indicating the timestepper (:QuasiAdamsBashforth2 or :SplitRungeKutta3)

Returns

• Recommended time step [s] for the given timestepper

Time steps are chosen to match computational cost between AB2 and RK schemes.

initial_buoyancy(z, p)

Compute the initial buoyancy profile as a function of depth.

initial_buoyancy(z, p)

Arguments

• z: Vertical position [m]
• p: Parameters named tuple containing ΔB (surface buoyancy gradient), h (decay scale), and Lz (domain depth)

Returns

• Initial buoyancy [m/s²] following an exponential profile

The profile follows ΔB * (exp(z/h) - exp(-Lz/h)) / (1 - exp(-Lz/h)), which provides a stable stratification that decays exponentially with depth.

internal_tide(timestepper::Symbol; free_surface, tracer_advection)

Set up and run the internal tide test case simulation.

internal_tide(timestepper; free_surface, tracer_advection)

Arguments

• timestepper: Symbol indicating the timestepper (:QuasiAdamsBashforth2, :SplitRungeKutta3, or :SplitRungeKutta4)

Keyword Arguments

• free_surface: Free surface formulation (default: SplitExplicitFreeSurface with 60 substeps)
• tracer_advection: Tracer advection scheme (default: 7th-order WENO)

Returns

• Simulation object after running to completion

This function sets up the internal tide test case described in the paper, which simulates tidal flow over a Gaussian seamount. The domain is initially stratified and forced by an oscillatory tidal forcing. The simulation runs for 40 days and outputs velocity, buoyancy, tracer fields, and dissipation diagnostics.

The test case isolates the role of time discretization in numerical mixing, as spatial advection plays a secondary role in this mostly linear configuration.

internal_tide_grid()

Construct the grid for the internal tide test case with a Gaussian seamount.

internal_tide_grid()

Returns

• ImmersedBoundaryGrid with periodic horizontal boundaries and a Gaussian seamount bottom

The grid spans from -L to L horizontally (periodic) and from -H to 0 vertically (bounded). A Gaussian seamount of height h₀ and width width is imposed using an immersed boundary method. The grid parameters are taken from internal_tide_parameters().

internal_tide_parameters()

Return a named tuple containing the default parameters for the internal tide test case.

internal_tide_parameters()

Returns

• Nx: Number of grid points in the horizontal (x) direction
• Nz: Number of grid points in the vertical (z) direction
• H: Domain depth [m]
• L: Domain half-length [m]
• h₀: Seamount height [m]
• width: Seamount width [m]
• T₂: Tidal period [s]
• ω₂: Tidal frequency [rad/s]
• ϵ: Excursion parameter (dimensionless)
• U₂: Characteristic tidal velocity [m/s]
• f: Coriolis parameter [s⁻¹]
• A₂: Tidal forcing amplitude [m/s²]
• Nᵢ²: Initial stratification [s⁻²]

These parameters define a 2-kilometer deep domain with a Gaussian seamount that interacts with an oscillatory tidal forcing, as described in the paper.

internal_tide_timestep(::Val{timestepper})

Return the recommended time step for the internal tide test case given a timestepper.

internal_tide_timestep(_)

Arguments

• timestepper: Symbol indicating the timestepper (:QuasiAdamsBashforth2, :SplitRungeKutta3, or :SplitRungeKutta4)

Returns

• Recommended time step [s] for the given timestepper

The time steps are chosen to match computational cost between AB2 and RK schemes while maintaining stability, as described in the paper.

load_channel(folder, case_number; arch)

Load and process channel simulation output data.

load_channel(folder, case_number; arch)

Arguments

• folder: Directory containing the simulation output files
• case_number: Case identifier string
• arch: Architecture to load data on (default: CPU())

Returns

• Dictionary containing velocity, buoyancy, volume, and energy diagnostics

This function loads channel simulation snapshots and computes kinetic energy, mean kinetic energy, and potential energy diagnostics for equilibrated channel flow analysis.

load_idealized_coast(folder, closure, suffix, timestepper)

Load and process idealized coast simulation output data.

load_idealized_coast(folder, closure, suffix, timestepper)

Arguments

• folder: Directory containing the simulation output files
• closure: Closure name string (e.g., "CATKE")
• suffix: Filename suffix string (e.g., "split_free_surface_")
• timestepper: Timestepper name string (e.g., "SplitRungeKutta3")

Returns

• Dictionary containing velocity, tracer, dissipation, and energy diagnostics

Similar to load_internal_tide, but includes additional fields for temperature and salinity dissipation diagnostics, as well as three-dimensional volume fields.

load_internal_tide(folder, timestepper, free_surface)

Load and process internal tide simulation output data.

load_internal_tide(folder, timestepper, free_surface)

Arguments

• folder: Directory containing the simulation output files
• timestepper: Timestepper name string (e.g., "SplitRungeKutta3")
• free_surface: Free surface type string (e.g., "split_free_surface_DFB")

Returns

• Dictionary containing:
  • :u, :w, :b, :η: Velocity, buoyancy, and free surface FieldTimeSeries
  • :VCCC, :VFCC, :VCCF: Volume FieldTimeSeries at different grid locations
  • :Abx, :Abz: Advective buoyancy dissipation in x and z directions
  • :Gbx, :Gbz: Buoyancy gradient squared in x and z directions
  • :abx, :abz, :abt: Volume-averaged dissipation rates
  • :gbx, :gbz, :gbt: Volume-averaged gradient squared
  • :κb: Effective numerical diffusivity
  • :KE, :MKE: Kinetic energy and mean kinetic energy time series
  • :η2: Free surface variance time series
  • :RPE, :APE: Reference and available potential energy time series

This function loads simulation output, computes volume fields accounting for free surface variations, and calculates various diagnostics used for analyzing numerical mixing.

mask(y, p)

Compute the sponge layer mask function for buoyancy relaxation.

mask(y, p)

Arguments

• y: Meridional position [m]
• p: Parameters named tuple containing Ly (domain length) and Lsponge (sponge layer width)

Returns

• Mask value between 0 and 1, increasing linearly from 0 at y = Ly - Lsponge to 1 at y = Ly

This function creates a sponge layer mask used for buoyancy restoring at the northern boundary.

print_progress(sim)

Print simulation progress information including completion percentage, iteration number, simulation time, wall clock time, maximum velocity components, and next time step.

print_progress(sim)

Arguments

• sim: The Simulation object to print progress for.

Returns

• nothing

This callback function is typically added to simulations using add_callback! to monitor progress during long-running simulations.

simulation_Δt(::Val{timestepper})

Return the recommended time step for the channel simulation given a timestepper.

simulation_Δt(_)

Arguments

• timestepper: Symbol indicating the timestepper (:QuasiAdamsBashforth2, :SplitRungeKutta3, or :SplitRungeKutta6)

Returns

• Recommended time step [s] for the given timestepper

Time steps are chosen to match computational cost between different timesteppers.

tidal_forcing(x, z, t, p)

Compute the tidal forcing amplitude at position (x, z) and time t.

tidal_forcing(x, z, t, p)

Arguments

• x: Horizontal position [m]
• z: Vertical position [m]
• t: Time [s]
• p: Parameters named tuple containing A₂ (forcing amplitude) and ω₂ (tidal frequency)

Returns

• Forcing amplitude [m/s²] as a function of time only: A₂ * sin(ω₂ * t)

This function implements a simple oscillatory tidal forcing that varies sinusoidally in time with amplitude A₂ and frequency ω₂.

u_stress(i, j, grid, clock, model_fields, p)

Compute the zonal wind stress at the surface.

u_stress(i, j, grid, clock, model_fields, p)

Arguments

• i, j: Grid indices
• grid: Grid object
• clock: Simulation clock
• model_fields: Model fields
• p: Parameters named tuple containing τ (stress magnitude) and Ly (domain length)

Returns

• Zonal wind stress [m²/s²] following a sine profile: -τ * sin(π * y / Ly)

This implements a sinusoidal wind stress profile that drives zonal flow, with maximum stress at the domain center and zero at the boundaries.

wind_stress(i, j, grid, clock, fields, p)

Compute the wind stress forcing at grid point (i, j).

wind_stress(i, j, grid, clock, fields, p)

Arguments

• i, j: Grid indices
• grid: Grid object
• clock: Simulation clock
• fields: Model fields
• p: Parameters named tuple containing τ₀ (stress amplitude) and f (frequency)

Returns

• Wind stress [m²/s²] that activates after 4 days: τ₀ * sin(f * t) if t > 4 days, else 0

This function implements a time-dependent wind stress that begins after a 4-day spin-up period and oscillates with frequency f and amplitude τ₀.

(ϵ::BuoyancyVarianceDissipation)(model)

Compute buoyancy variance dissipation and update flux caches for the next time step.

BuoyancyVarianceDissipation(grid; Uⁿ⁻¹, Uⁿ)

Arguments

• model: The HydrostaticFreeSurfaceModel containing temperature and salinity tracers

Returns

• nothing (modifies ϵ in place)

This function is called as a simulation callback to compute the advective dissipation of buoyancy variance at each time step. It first computes dissipation from previously cached fluxes, then updates the flux caches for the next time step.

The model must have :T and :S tracers to compute buoyancy.

c_grid_vector(grid)

Create a named tuple of C-grid vector fields (x, y, z components at face locations).

c_grid_vector(grid)

Arguments

• grid: Grid on which to create the fields

Returns

• Named tuple with x, y, z fields at XFace, YFace, and ZFace locations respectively

This helper function creates the vector fields needed to store advective fluxes and dissipation rates at their proper staggered grid locations.

cache_fluxes!(dissipation, model)

Cache advective fluxes and update velocity fields for dissipation computation.

cache_fluxes!(dissipation, model)

Arguments

• dissipation: BuoyancyVarianceDissipation object
• model: The model containing velocities and tracers

Returns

• nothing (modifies dissipation in place)

This function updates the velocity fields (Uⁿ, Uⁿ⁻¹) and caches the advective buoyancy fluxes (Fⁿ, Fⁿ⁻¹) needed for computing dissipation at the next time step. The caching procedure depends on the timestepper type and stage.

compute_dissipation!(dissipation, model)

Compute the numerical dissipation of buoyancy variance from previously calculated advective fluxes.

compute_dissipation!(dissipation, model)

Arguments

• dissipation: BuoyancyVarianceDissipation object containing cached fluxes
• model: The model containing tracers and timestepper information

Returns

• nothing (modifies dissipation.advective_production in place)

The dissipation is computed using the formulation:

A = 2 * δc★ * F - U δc²  # For advective dissipation

where F is the flux associated with advection, U is the advecting velocity, c★ and c² are functions of the buoyancy field.

For multi-stage timesteppers, the fluxes and velocities must account for the substepping procedure. For AB2: F = 1.5 Fⁿ - 0.5 Fⁿ⁻¹ and U = 1.5 Uⁿ - 0.5 Uⁿ⁻¹. For RK schemes, the fluxes from the appropriate substep are used.

This method is adapted from the exact variance budget methodology described in the paper's appendix.

finally_cache_fluxes!(dissipation, model)

Cache advective buoyancy fluxes from temperature and salinity tracers.

finally_cache_fluxes!(dissipation, model)

Arguments

• dissipation: BuoyancyVarianceDissipation object
• model: The model containing tracers and advection scheme

Returns

• nothing (modifies dissipation.advective_fluxes in place)

This function computes the advective fluxes of buoyancy by combining temperature and salinity fluxes using the linear equation of state. The fluxes are computed at face locations and cached for use in dissipation computation.

flatten_dissipation_fields(t::BuoyancyVarianceDissipation)

Flatten the dissipation fields into a named tuple for output.

flatten_dissipation_fields(t)

Arguments

• t: BuoyancyVarianceDissipation object

Returns

• Named tuple containing Abx, Aby, Abz fields representing advective buoyancy dissipation in the x, y, and z directions respectively

This function extracts the dissipation fields from the internal C-grid vector structure into a flat named tuple suitable for use as simulation outputs.

flux_parameters(grid)

Compute kernel parameters for flux computation based on grid topology.

flux_parameters(grid)

Arguments

• grid: Grid object

Returns

• KernelParameters object with appropriate index ranges for flux computation

This function determines the correct index ranges for computing fluxes at face locations, accounting for flat (1D/2D) dimensions in the grid topology.