Buoyancy as sum of tracers

[raphaelouillon]

Hi all,
I am trying to setup a simulation of a gravity-driven flow where different scalar fields contribute to the buoyancy. From what I understand, this is done by setting buoyancy=BuoyancyTracer(), in which case :b is the buoyancy tracer. What I am hoping to do is have :b1,:b2,:b3... all contribute to the buoyancy. Is this currently possible? Alternatively, can one define a custom "equation of state"? Thanks!

[simone-silvestri]

Hi Raphael, an "easy" way to do this is to define your own BuoyancyModel following the BuoyancyTracer struct.
To do that you have to import and extend the functions that are found in src/BuoyancyModels/buoyancy_tracer.jl.
For example, if you want your buoyancy to be the sum of 3 tracers

using Oceananigans.Buoyancy: AbstractBuoyancyModel
import Oceananigans.Buoyancy: required_tracers, buoyancy_perturbation,  ∂x_b,  ∂y_b,  ∂z_b, top_buoyancy_flux, bottom_buoyancy_flux

struct MyBuoyancyTracer <: AbstractBuoyancyModel{Nothing} end

required_tracers(::MyBuoyancyTracer) = (:b1, :b2, :b3)

@inline buoyancy_perturbation(i, j, k, grid, ::MyBuoyancyTracer, C) = @inbounds C.b1[i, j, k] + C.b2[i, j, k] + C.b3[i, j, k]

@inline ∂x_b(i, j, k, grid, ::MyBuoyancyTracer, C) = ∂xᶠᶜᶜ(i, j, k, grid, C.b1) + ∂xᶠᶜᶜ(i, j, k, grid, C.b2) + ∂xᶠᶜᶜ(i, j, k, grid, C.b3)
@inline ∂y_b(i, j, k, grid, ::MyBuoyancyTracer, C) = ∂yᶜᶠᶜ(i, j, k, grid, C.b1) + ∂yᶜᶠᶜ(i, j, k, grid, C.b2) + ∂yᶜᶠᶜ(i, j, k, grid, C.b3)
@inline ∂z_b(i, j, k, grid, ::MyBuoyancyTracer, C) = ∂zᶜᶜᶠ(i, j, k, grid, C.b1) + ∂zᶜᶜᶠ(i, j, k, grid, C.b2) + ∂zᶜᶜᶠ(i, j, k, grid, C.b3)

@inline    top_buoyancy_flux(i, j, grid, ::MyBuoyancyTracer, top_tracer_bcs, clock, fields)    = getbc(top_tracer_bcs.b1, i, j, grid, clock, fields) + getbc(top_tracer_bcs.b2, i, j, grid, clock, fields) + getbc(top_tracer_bcs.b3, i, j, grid, clock, fields)
@inline bottom_buoyancy_flux(i, j, grid, ::MyBuoyancyTracer, bottom_tracer_bcs, clock, fields) = getbc(bottom_tracer_bcs.b1, i, j, grid, clock, fields) + getbc(bottom_tracer_bcs.b2, i, j, grid, clock, fields) + getbc(bottom_tracer_bcs.b3, i, j, grid, clock, fields)

Note that the only function required for buoyancy calculation in the Navier-Stokes equations is buoyancy_perturbation.
Buoyancy derivatives (∂x_b, ∂y_b, ∂z_b) and buoyancy boundary conditions (top_buoyancy_flux and bottom_buoyancy_flux) are only used for parameterizations. If you don't have particular parameterizations you can avoid extending those functions.

You can then use this buoyancy in the model constructor

model = NonhydrostaticModel(; grid, buoyancy = MyBuoyancyTracer(), tracers = (:b1, :b2, :b3))

This might be a better solution than defining your own equation of state. This is because the machinery in place assumes that an equation of state requires a :T and an :S, so you would have to carry those additional fields over

[glwagner]

Just want to add a bit about terminology to help --- in terms of the way we normally think of "equations of state", this approach is defining a custom equation of state.

The concept of "buoyancy" has two components which is encapsulated in this API: required_tracers (the tracers that contribute to buoyancy), and buoyancy_perturbation (the relationship between those tracers and buoyancy). So the function buoyancy_perturbation implements what is often called the "equation of state".

We have provided a special buoyancy type SeawaterBuoyancy, which is designed for equations of state of two tracers called T and S, representing temperature and salinity.

But there are plenty of other equations of state depending on other tracers (concentration of sediments of various kinds, temperature, salinity, microorganisms...) that are not what we have chosen to call "seawater".

[raphaelouillon]

Got it, thanks for the additional information @glwagner!

[mkrapp]

Hey @glwagner. Thanks for the context. I recently came across a related problem.

In my model, I want to add a concentration ($C$) as a tracer. For example, some sediment with density, say $\rho_s$, suspended in seawater. What is the best way to implement this? Do I have to build it using AbstractBuoyancyModel, as suggested, or is there a quick hack for the SeawaterBuoyancy type?

The modified density for sediment with concentration $C$ would be

$$ \rho = \rho_{sw}(1 - C)[1 + \beta S - \alpha T] + \rho_s C $$

[simone-silvestri]

Extending AbstractBuoyancyModel with your own buoyancy model would be the fastest and easiest way. SeawaterBuoyancy assumes that buoyancy depends only on the two tracers T and S. It is possible to modify it to include another tracer but that would require importing and changing more functions than building your own custom buoyancy model.

Note that the name of function to extend was refactored to buoyancy_perturbationᶜᶜᶜ. The function signature is

@inline function buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, b::MyBuoyancyModel, tracers)

where tracers is a named tuple containing all the tracers passed to the model (in your case (T, S, C)).

[glwagner]

Also worth noting that you can implement this kind of buoyancy model as a Forcing function if you're using NonhydrostaticModel. (For HydrostaticFreeSurfaceModel, you can't use forcing for a buoyancy model.)

[raphaelouillon]

Thanks @simone-silvestri, that's great!