GitHub

Buoyancy models and equations of state

The buoyancy option selects how buoyancy is treated in NonhydrostaticModels and HydrostaticFreeSurfaceModels (ShallowWaterModels do not have that option given the physics of the model). There are currently three alternatives:

  1. No buoyancy (and no gravity).
  2. Evolve buoyancy as a tracer.
  3. Seawater buoyancy: evolve temperature $T$ and salinity $S$ as tracers with a value for the gravitational acceleration $g$ and an equation of state of your choosing.

No buoyancy

To turn off buoyancy (and gravity) you can simply pass buoyancy = nothing to the model constructor. For example to create a NonhydrostaticModel:

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

julia> model = NonhydrostaticModel(; grid, buoyancy=nothing)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
└── coriolis: Nothing

The option buoyancy = nothing is the default for NonhydrostaticModel, so omitting the buoyancy keyword argument from the NonhydrostaticModel constructor yields the same:

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

The same is true for HydrostaticFreeSurfaceModel,

julia> model = HydrostaticFreeSurfaceModel(; grid)
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
├── advection scheme: 
│   └── momentum: Vector Invariant, Dimension-by-dimension reconstruction
└── coriolis: Nothing

Buoyancy as a tracer

Both NonhydrostaticModel and HydrostaticFreeSurfaceModel support evolving a buoyancy tracer by including :b in tracers and specifying buoyancy = BuoyancyTracer():

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

Similarly for a HydrostaticFreeSurfaceModel with buoyancy as a tracer:

julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b)
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
├── advection scheme:
│   ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction
│   └── b: Centered(order=2)
└── coriolis: Nothing

Seawater buoyancy

NonhydrostaticModel and HydrostaticFreeSurfaceModel support modeling the buoyancy of seawater as a function of the gravitational acceleration, the conservative temperature $T$, and the absolute salinity $S$. The relationship between $T$, $S$, the geopotential height, and the density perturbation from a reference value is called the equation_of_state.

Specifying buoyancy = SeawaterBuoyancy() returns a buoyancy model with a linear equation of state, Earth standard gravitational_acceleration = 9.80665 (in S.I. units $\text{m}\,\text{s}^{-2}$) and requires to add :T and :S as tracers:

julia> model = NonhydrostaticModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection()
└── coriolis: Nothing

and the same is true for HydrostaticFreeSurfaceModel,

julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection()
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
├── advection scheme:
│   ├── momentum: Vector Invariant, Dimension-by-dimension reconstruction
│   ├── T: Centered(order=2)
│   └── S: Centered(order=2)
└── coriolis: Nothing

To model flows near the surface of Europa where gravitational_acceleration = 1.3 $\text{m}\,\text{s}^{-2}$, we might alternatively specify

julia> buoyancy = SeawaterBuoyancy(gravitational_acceleration=1.3)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 1.3
└── equation_of_state: LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078)

julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=(:T, :S))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=1.3 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection()
└── coriolis: Nothing

for example.

Linear equation of state

To specify the thermal expansion and haline contraction coefficients $\alpha = 2 \times 10^{-3} \; \text{K}^{-1}$ and $\beta = 5 \times 10^{-4} \text{psu}^{-1}$,

julia> buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-3, haline_contraction=5e-4))
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: LinearEquationOfState(thermal_expansion=0.002, haline_contraction=0.0005)

Idealized nonlinear equations of state

Instead of a linear equation of state, six idealized (second-order) nonlinear equations of state as described by Roquet et al. (2015) may be used. These equations of state are provided via the SeawaterPolynomials.jl package.

julia> using SeawaterPolynomials.SecondOrderSeawaterPolynomials

julia> eos = RoquetEquationOfState(:Freezing)
BoussinesqEquationOfState{Float64}:
    ├── seawater_polynomial: SecondOrderSeawaterPolynomial{Float64}
    └── reference_density: 1024.6
    
julia> eos.seawater_polynomial # the density anomaly
ρ' = 0.7718 Sᴬ - 0.0491 Θ - 0.005027 Θ² - 2.5681e-5 Θ Z + 0.0 Sᴬ² + 0.0 Sᴬ Z + 0.0 Sᴬ Θ

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

TEOS-10 equation of state

A high-accuracy 55-term polynomial approximation to the TEOS-10 equation of state suitable for use in Boussinesq models as described by Roquet et al. (2015) is implemented in the SeawaterPolynomials.jl package and may be used.

julia> using SeawaterPolynomials.TEOS10

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

The direction of gravitational acceleration

To simulate gravitational accelerations that don't align with the vertical (z) coordinate, we use BuoyancyForce(formulation; gravity_unit_vector), wherein the buoyancy formulation can be BuoyancyTracer, SeawaterBuoyancy, etc, in addition to the gravity_unit_vector. For example,

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

julia> θ = 45; # degrees

julia> g̃ = (0, sind(θ), cosd(θ));

julia> buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector=g̃)
BuoyancyForce:
├── formulation: BuoyancyTracer
└── gravity_unit_vector: (0.0, 0.707107, 0.707107)

julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=:b)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered(order=2)
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = (0.0, 0.707107, 0.707107)
└── coriolis: Nothing