OceanBioME was developed with generous support from the Cambridge Centre for Climate Repair CCRC and the Gordon and Betty Moore Foundation as a tool to study the effectiveness and impacts of ocean carbon dioxide removal (CDR) strategies.
OceanBioME is a flexible modelling environment written in Julia for modelling the coupled interactions between ocean biogeochemistry, carbonate chemistry, and physics. OceanBioME can be run as a stand-alone box model, or coupled with Oceananigans.jl to run as a 1D column model or with 2 and 3D physics.
First, download and install Julia
From the Julia prompt (REPL), type:
julia> using Pkg
julia> Pkg.add("OceanBioME")
As a simple example lets run a Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model in a two-dimensional simulation of a buoyancy front:
using OceanBioME, Oceananigans
using Oceananigans.Units
grid = RectilinearGrid(CPU(), size=(256, 32), extent=(500, 100), topology=(Bounded, Flat, Bounded))
biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid, open_bottom=true)
model = NonhydrostaticModel(; grid, biogeochemistry, buoyancy=BuoyancyTracer(), tracers=:b, advection=WENO(; grid), closure = AnisotropicMinimumDissipation())
báµ¢(x, y, z) = ifelse(x < 250, 1e-4, 1e-3)
set!(model, b = báµ¢, N = 5.0, P = 0.1, Z = 0.1, T = 18.0)
simulation = Simulation(model; Δt=1.0, stop_time=3hours)
wizard = TimeStepWizard(cfl=0.3, max_change=1.5)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))
simulation.output_writers[:tracers] = JLD2OutputWriter(model, model.tracers, filename = "buoyancy_front.jld2", schedule = TimeInterval(1minute), overwrite_existing=true)
run!(simulation)
We can then visualise this:
using CairoMakie
b = FieldTimeSeries("buoyancy_front.jld2", "b")
P = FieldTimeSeries("buoyancy_front.jld2", "P")
n = Observable(1)
b_lims = (minimum(b), maximum(b))
P_lims = (minimum(P), maximum(P))
b_plt = @lift b[1:grid.Nx, 1, 1:grid.Nz, $n]
P_plt = @lift P[1:grid.Nx, 1, 1:grid.Nz, $n]
fig = Figure(resolution = (1600, 160 * 4))
supertitle = Label(fig[0, :], "t = 0.0")
ax1 = Axis(fig[1, 1], xlabel = "x (m)", ylabel = "z (m)", title = "Buouyancy pertubation (m / s)", width = 1400)
ax2 = Axis(fig[2, 1], xlabel = "x (m)", ylabel = "z (m)", title = "Phytoplankton concentration (mmol N / m³)", width = 1400)
hm1 = heatmap!(ax1, xnodes(grid, Center(), Center(), Center())[1:grid.Nx], znodes(grid, Center(), Center(), Center())[1:grid.Nz], b_plt, colorrange = b_lims, colormap = :batlow, interpolate=true)
hm2 = heatmap!(ax2, xnodes(grid, Center(), Center(), Center())[1:grid.Nx], znodes(grid, Center(), Center(), Center())[1:grid.Nz], P_plt, colorrange = P_lims, colormap = Reverse(:bamako), interpolate=true)
Colorbar(fig[1, 2], hm1)
Colorbar(fig[2, 2], hm2)
record(fig, "buoyancy_front.gif", 1:length(b.times)) do i
n[] = i
msg = string("Plotting frame ", i, " of ", length(b.times))
print(msg * " \r")
supertitle.text = "t=$(prettytime(b.times[i]))"
end
In this example OceanBioME
is providing the biogeochemistry
and the remainder is taken care of by Oceananigans
. For comprehensive documentation of the physics modelling see Oceananigans' Documentation, and for biogeochemistry and other features we provide read below.
See the documentation for full description and examples.