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=(64, 64, 64), extent=(1, 1, 1));

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

buoyancy=nothing is the default option forNonhydrostaticModel, so omitting buoyancy from the NonhydrostaticModel constructor yields an identical result:

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

To create a HydrostaticFreeSurfaceModel without a buoyancy term we explicitly specify buoyancy=nothing flag. The default tracers T and S for HydrostaticFreeSurfaceModel may be eliminated when buoyancy=nothing by specifying tracers=():

julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=nothing, tracers=())
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×64 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
└── 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: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with -ĝ = ZDirection
└── coriolis: Nothing

We follow the same pattern to create a HydrostaticFreeSurfaceModel with buoyancy as a tracer:

julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b)
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with -ĝ = ZDirection
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
└── coriolis: Nothing

Seawater buoyancy

NonhydrostaticModel and HydrostaticFreeSurfaceModel support modeling the buoyancy of seawater as a function of gravitational acceleration, conservative temperature $T$ and 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 ($\text{m}\,\text{s}^{-2}$) and requires the tracers :T and :S:

julia> model = NonhydrostaticModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×64 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 -ĝ = ZDirection
└── coriolis: Nothing

With HydrostaticFreeSurfaceModel,

julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×64 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 -ĝ = ZDirection
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
└── coriolis: Nothing

is identical to the default,

julia> model = HydrostaticFreeSurfaceModel(; grid)
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×64 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 -ĝ = ZDirection
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
└── 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=buoyancy, tracers=(:T, :S))
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (T, S)
├── closure: Nothing
├── buoyancy: SeawaterBuoyancy with g=1.3 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with -ĝ = ZDirection
└── 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, five idealized (second-order) nonlinear equation of state as described by Fabien Roquet, Gurvan Madec, Laurent Brodeau, J. Nycander (2015) may be used. These equations of state are provided via the SeawaterPolynomials.jl package.

julia> using SeawaterPolynomials.SecondOrderSeawaterPolynomials

julia> eos = RoquetSeawaterPolynomial(:Freezing)
0.7718 Sᴬ - 0.0491 Θ + 0.0 Θ² - 2.5681e-5 Θ Z + 0.0 Sᴬ² - 0.005027 Sᴬ Z + 0.0 Sᴬ Θ

julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation of state: SecondOrderSeawaterPolynomial{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 F. Roquet, G. Madec, Trevor J. McDougall, Paul M. Barker (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 wrap the buoyancy model in Buoyancy() function call, which takes the keyword arguments model and gravity_unit_vector,

julia> θ = 45; # degrees

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

julia> model = NonhydrostaticModel(; grid, 
                                   buoyancy=Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃), 
                                   tracers=:b)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×64 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with -ĝ = Tuple{Int64, Float64, Float64}
└── coriolis: Nothing