How to have vertical drift in only one speices of a biogeochemistry model?

[francispoulin]

Below is a simple minimum working example that shows how you can add vertical advection to nutrients using a built in advection scheme. This uses first order upwinding to make sure that it preserves positivity and does not go beyond the maximum.

What I actually want to do is have an NPZ model say, and only have vertical drift of one species on top of the regular advection.

I think to do this i need to specify forcing to one of the species. Can someone show me how to do this using a built in advection scheme, rather than having to rebuild my own?

using Oceananigans
using Oceananigans.Advection
using  Printf

Δt      = 0.1         
T_end   = 20.0
Δt_save = 1.0

Nz = 100
grid = RectilinearGrid(size=(Nz,), z=(-1.0, 0.0), topology=(Flat, Flat, Periodic))

const w0 = -0.005
wfield = ZFaceField(grid);  set!(wfield, w0)

N_initial(z) = (tanh(10*(z + 0.5)) + 1) / 2

model = NonhydrostaticModel(;
    grid,
    tracers     = (:N,),
    advection   = WENO(bounds=(0.0, 1.0)), 
    timestepper = :RungeKutta3              
)

set!(model; w = wfield, N = N_initial)

simulation = Simulation(model; Δt=Δt, stop_time=T_end)

progress(sim) = @printf("t=%s  min(N)=%.3e  max(N)=%.3e\n",
                        prettytime(time(sim)), minimum(sim.model.tracers.N), maximum(sim.model.tracers.N))
simulation.callbacks[:progress] = Callback(progress, TimeInterval(Δt_save))

run!(simulation)

[jagoosw]

I think the easiest way to do this would be with AdvectiveForcing where you specify u, v, and/or w as fields and then put it in like a normal forcing.

If you wanted to use the built in biogeochemistry infrastructure you can define your model like:

using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry
struct NPZModel <: AbstractContinuousFormBiogeochemistry
    parameters...
end

import Oceananigans.Biogeochemistry: biogeochemical_drift_velocity

biogeochemical_drift_velocity(bgc::NPZModel, ::Val{:D}) = 
    (u = ZeroField(), v = ZeroField(), w = ConstantField(-0.005))

[jagoosw]

I think AdevectiveForcing is combined into the inbuilt advection scheme, and the biogeochemistry way definitely does.

[francispoulin]

Hello Jago, and thanks for your reply! In this example here there isn't any motion so it's pretty easy. But very soon we want to consider a moving fluid, and then add a sinking velocity to one of the speices.

What you suggested looks attractive but I don't know how this would then be added to the fluid velocity, which must be done somewhere?

Also, if we want the velocity only affecting one species, is that a problem?

Our model is a new model and not quite an NPZ model so we do have to create it from scratch. But I am happy to use the framekwork that exists, as its easier for us and others hopefully will benefit from it in the future.

[jagoosw]

No problem!

The fields all get added together in the tendency kernel and the total used in the advection:

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl

Lines 274 to 292 in 4e9d709

	biogeochemical_velocities = biogeochemical_drift_velocity(biogeochemistry, val_tracer_name)
	total_velocities = sum_of_velocities(velocities, background_fields.velocities, biogeochemical_velocities)
	total_velocities = with_advective_forcing(forcing, total_velocities)
	@inbounds c = tracers[tracer_index]
	@inbounds background_fields_c = background_fields.tracers[tracer_index]
	closure_c = if background_fields isa BackgroundFieldsWithClosureFluxes
	sum_fields(c, background_fields_c)
	else
	c
	end
	closure_velocities = assemble_closure_velocities(velocities, background_fields)
	closure_model_fields = merge(closure_velocities, tracers, auxiliary_fields)
	model_fields = merge(velocities, tracers, auxiliary_fields)
	return ( - div_Uc(i, j, k, grid, advection, total_velocities, c)

So they can be added to any model velocity field

[francispoulin]

Oh very nice! So now I see very clearly that if we define biogeochemical_drift_velocity, it will get added to the velocity when needed. Great!

I also see that this is all explained in the docs, https://oceanbiome.github.io/OceanBioME.jl/stable/model_implementation/

Sorry that I didn't find this before, but thanks for pointing me in the right direction!