1D sediment column

[stanleesocca]

Hi @jagoosw and the rest. Thanks for the package. I have been testing it lately and assessing its capabilities. I see there is a sediment model included in the model. However, it seems to be only a reflective boundary or 0D box model. Is it possible for OceanBioME to have a 1D sediment model (some like the classic OMEXDIA model)?

How difficult would it be to get it integrated in the oceananigans/OceanBioMe framework?

Thanks

[jagoosw]

Hi @stanleesocca, thank you for your interest in OceanBioME!

We don't currently have any 1D sediment models but I think it wouldn't be too hard to implement one with the sediment infrastructure we've built.

I'm currently rewriting the sediment module here, which is nearly done. After that, I think you would be able to make a 1D model by defining a sediment biogeochemistry and writing new methods for these functions:

OceanBioME.jl/src/Sediments/compute_tendencies.jl

Lines 31 to 50 in b92d540

	@kernel function compute_sediment_tendency!(G, grid, args)
	i, j = @index(Global, NTuple)
	# tendencies in X/m²/s
	@inbounds G[i, j, 1] = sediment_tendencies(i, j, grid, args...)
	end
	@inline sediment_tendencies(i, j, grid, val_name, biogeochemistry, fields, tracked_fields, clock) =
	biogeochemistry(i, j, grid, val_name, clock, fields, tracked_fields)
	@inline function sediment_tendencies(i, j, grid, val_name, biogeochemistry::ACFSBGC, fields, tracked_fields, clock)
	model_fields = merge(fields, tracked_fields)
	field_values = extract_biogeochemical_fields(i, j, 1, grid, model_fields, keys(model_fields))
	x = xnode(i, j, 1, grid, Center(), Center(), Center())
	y = ynode(i, j, 1, grid, Center(), Center(), Center())
	return biogeochemistry(val_name, x, y, clock.time, field_values...)
	end

and

OceanBioME.jl/src/Sediments/timesteppers.jl

Lines 15 to 28 in b92d540

	function ab2_step!(model::BiogeochemicalSediment, Δt)
	grid = model.grid
	arch = architecture(grid)
	model_fields = prognostic_fields(model)
	χ = model.timestepper.χ
	for (i, field) in enumerate(model_fields)
	kernel_args = (field, Δt, χ, model.timestepper.Gⁿ[i], model.timestepper.G⁻[i])
	launch!(arch, grid, :xy, ab2_step_flat_field!, kernel_args...; exclude_periphery=true)
	end
	return nothing
	end

If this was something you were interested in doing, I would be happy to help and we could think about moving something like AbstractOneDSediment into the source code.

[stanleesocca]

Hi @jagoosw Thanks for the reply. yes, I would be interested in having a 1D sediment model here. I think we can start by first testing with a rather simple 1D model with one or two species before hitting the multi-component diagenetic models as in the literature. Bear in the mind, these sediment biogeochemical models are highly stiff and from the way I see it, Oceananigans/OceanBioME seems to have made deliberate design choice in not using the suite of solvers from SciML. I already have a test-like OMEXDIA julia code integrated with Chris's SciML. If needed I can share that with you.

In any case, I will need to study the OceanBioME/oceananigans code a bit to understand how things are called. Perhaps you can help with guiding me through the maze here.

[jagoosw]

Yeah, that sounds ideal. I'm quite busy at the moment but am happy to help get a PR going to add the missing pieces in the source code.

Timestepping bgc models is definitely a bit of a problem here. I thought about this a bit and think the best solution would be to substep bgc models with SciML or DifferentialEquations timesteppers, but there is some work that needs to be done to get that working in Oceananigans. I made a PR here but tabled it for now. Maybe a bit of help is that sediments do get sub stepped at each fluid time step, so you can RK3 step the sediment for each RK3 substep of the fluid solver (https://github.com/OceanBioME/OceanBioME.jl/blob/jsw/refactor-sediment/src/Sediments/update_state.jl, https://github.com/OceanBioME/OceanBioME.jl/blob/jsw/refactor-sediment/src/Sediments/timesteppers.jl).

As for getting into the code: in Oceananigans BGC models provide tendancies ($\partial X/\partial t += ...$) for tracer fields, but also can have an update_biogeochemical_state! function which is called between every substep. In OceanBioME, we call this on various things (particles, sediments, modifiers), and that is when the sediment gets integrated, and the tendency contributions to the tracers get added.

When this happens the method in https://github.com/OceanBioME/OceanBioME.jl/blob/jsw/refactor-sediment/src/Sediments/update_state.jl gets called and if you look through the code in that folder hopefully you can see how sediments currently work.

Let me know if it would be helpful for me to tell you more, or if it would be useful to talk.