GitHub

One-dimensional column example

In this example we setup a simple 1D column with the LOBSTER biogeochemical model and observe its evolution. The example demonstrates:

  • How to setup OceanBioME's biogeochemical models
  • How to visualise results

This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first.

Install dependencies

First we check we have the dependencies installed

using Pkg
pkg"add OceanBioME, Oceananigans, CairoMakie"

Model setup

We load the packages and choose the default LOBSTER parameter set

using OceanBioME, Oceananigans, Printf
using Oceananigans.Fields: FunctionField, ConstantField
using Oceananigans.Units

const year = years = 365days

Surface PAR and turbulent vertical diffusivity based on idealised mixed layer depth

Setting up idealised functions for PAR and diffusivity (details here can be ignored but these are typical of the North Atlantic)

@inline PAR⁰(x, y, t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2

@inline H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0)

@inline fmld1(t) = H(t, 50days, year) * (1 / (1 + exp(-(t - 100days) / 5days))) * (1 / (1 + exp((t - 330days) / 25days)))

@inline MLD(t) = - (10 + 340 * (1 - fmld1(year - eps(year)) * exp(-mod(t, year) / 25days) - fmld1(mod(t, year))))

@inline κₜ(x, y, z, t) = 1e-2 * (1 + tanh((z - MLD(t)) / 10)) / 2 + 1e-4

@inline temp(x, y, z, t) = 2.4 * cos(t * 2π / year + 50days) + 10

Grid

Define the grid.

grid = RectilinearGrid(size = (1, 1, 50), extent = (20meters, 20meters, 200meters))
1×1×50 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 1×1×3 halo
├── Periodic x ∈ [0.0, 20.0)   regularly spaced with Δx=20.0
├── Periodic y ∈ [0.0, 20.0)   regularly spaced with Δy=20.0
└── Bounded  z ∈ [-200.0, 0.0] regularly spaced with Δz=4.0

Model

First we define the biogeochemical model including carbonate chemistry (for which we also define temperature ($T$) and salinity ($S$) fields) and scaling of negative tracers(see discussion in the positivity preservation) and then setup the Oceananigans model with the boundary condition for the DIC based on the air-sea CO₂ flux.

biogeochemistry = LOBSTER(; grid,
                            surface_photosynthetically_active_radiation = PAR⁰,
                            carbonates = true,
                            scale_negatives = true)

CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition()

clock = Clock(; time = 0.0)
T = FunctionField{Center, Center, Center}(temp, grid; clock)
S = ConstantField(35.0)

model = NonhydrostaticModel(; grid,
                              clock,
                              closure = ScalarDiffusivity(ν = κₜ, κ = κₜ),
                              biogeochemistry,
                              boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),),
                              auxiliary_fields = (; T, S))

set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2409.0)

Simulation

Next we setup a simulation and add some callbacks that:

  • Show the progress of the simulation
  • Store the model and particles output
simulation = Simulation(model, Δt = 3minutes, stop_time = 100days)

progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n",
                                iteration(sim),
                                prettytime(sim),
                                prettytime(sim.Δt),
                                prettytime(sim.run_wall_time))

simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(10days))

filename = "column"
simulation.output_writers[:profiles] = JLD2Writer(model, model.tracers,
                                                  filename = "$filename.jld2",
                                                  schedule = TimeInterval(1day),
                                                  overwrite_existing = true)
JLD2Writer scheduled on TimeInterval(1 day):
├── filepath: column.jld2
├── 9 outputs: (NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk)
├── array_type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 60.8 KiB

Also output the surface CO₂ flux

@inline qCO₂_kernel(i, j, k, grid, clock, DIC, Alk, T, S, carbon_boundary_condition) =
    carbon_boundary_condition.condition.func(i, j, grid, clock, (; DIC, Alk, T, S))

qCO₂ = KernelFunctionOperation{Center, Center, Face}(qCO₂_kernel,
                                                     grid,
                                                     clock,
                                                     model.tracers.DIC,
                                                     model.tracers.Alk,
                                                     model.auxiliary_fields.T,
                                                     model.auxiliary_fields.S,
                                                     CO₂_flux)

simulation.output_writers[:carbon_flux] = JLD2Writer(model, (; qCO₂),
                                                     indices = (:, :, grid.Nz),
                                                     filename = filename * "_carbon.jld2",
                                                     schedule = TimeInterval(1day),
                                                     overwrite_existing = true)

┌ Warning: Attempting to store typeof(Main.var"##225".κₜ).
│ JLD2 only stores functions by name.
│  This may not be useful for anonymous functions.
└ @ JLD2 ~/.julia-oceananigans/packages/JLD2/WDhXU/src/data/writing_datatypes.jl:447

Run!

We are ready to run the simulation

run!(simulation)
[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 3 minutes, wall time: 0 seconds
[ Info:     ... simulation initialization complete (21.042 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (13.736 seconds).
Iteration: 4800, time: 10 days, Δt: 3 minutes, wall time: 56.375 seconds
Iteration: 9600, time: 20 days, Δt: 3 minutes, wall time: 1.289 minutes
Iteration: 14400, time: 30 days, Δt: 3 minutes, wall time: 1.633 minutes
Iteration: 19200, time: 40 days, Δt: 3 minutes, wall time: 1.983 minutes
Iteration: 24000, time: 50 days, Δt: 3 minutes, wall time: 2.336 minutes
Iteration: 28800, time: 60 days, Δt: 3 minutes, wall time: 2.676 minutes
Iteration: 33600, time: 70 days, Δt: 3 minutes, wall time: 2.989 minutes
Iteration: 38400, time: 80 days, Δt: 3 minutes, wall time: 3.329 minutes
Iteration: 43200, time: 90 days, Δt: 3 minutes, wall time: 3.659 minutes
[ Info: Simulation is stopping after running for 3.992 minutes.
[ Info: Simulation time 100 days equals or exceeds stop time 100 days.
Iteration: 48000, time: 100 days, Δt: 3 minutes, wall time: 3.992 minutes

Load saved output

Now we can load the results and do some post processing to diagnose the air-sea CO₂ flux. Hopefully, this looks different to the example without kelp!

   P = FieldTimeSeries(filename * ".jld2", "P")
 NO₃ = FieldTimeSeries(filename * ".jld2", "NO₃")
   Z = FieldTimeSeries(filename * ".jld2", "Z")
sPOM = FieldTimeSeries(filename * ".jld2", "sPOM")
bPOM = FieldTimeSeries(filename * ".jld2", "bPOM")
 DIC = FieldTimeSeries(filename * ".jld2", "DIC")
 Alk = FieldTimeSeries(filename * ".jld2", "Alk")

air_sea_CO₂_flux = FieldTimeSeries(filename * "_carbon.jld2", "qCO₂")

x, y, z = nodes(P)
times = P.times

We compute the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth that corresponds to k = grid.Nz - 20.

carbon_export = zeros(length(times))

using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity

for (n, t) in enumerate(times)
    clock.time = t

    carbon_export[n] = (sPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] +
                        bPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry)
end

Both air_sea_CO₂_flux and carbon_export are in units mmol CO₂ / (m² s).

Plot

Finally, we plot!

using CairoMakie

fig = Figure(size = (1000, 1500), fontsize = 20)

axis_kwargs = (xlabel = "Time (days)", ylabel = "z (m)", limits = ((0, times[end] / days), (-150meters, 0)))

axP = Axis(fig[1, 1]; title = "Phytoplankton concentration (mmol N / m³)", axis_kwargs...)
hmP = heatmap!(times / days, z, interior(P, 1, 1, :, :)', colormap = :batlow)
Colorbar(fig[1, 2], hmP)

axNO₃ = Axis(fig[2, 1]; title = "Nitrate concentration (mmol N / m³)", axis_kwargs...)
hmNO₃ = heatmap!(times / days, z, interior(NO₃, 1, 1, :, :)', colormap = :batlow)
Colorbar(fig[2, 2], hmNO₃)

axZ = Axis(fig[3, 1]; title = "Zooplankton concentration (mmol N / m³)", axis_kwargs...)
hmZ = heatmap!(times / days, z, interior(Z, 1, 1, :, :)', colormap = :batlow)
Colorbar(fig[3, 2], hmZ)

axD = Axis(fig[4, 1]; title = "Detritus concentration (mmol N / m³)", axis_kwargs...)
hmD = heatmap!(times / days, z, interior(sPOM, 1, 1, :, :)' .+ interior(bPOM, 1, 1, :, :)', colormap = :batlow)
Colorbar(fig[4, 2], hmD)

CO₂_molar_mass = (12 + 2 * 16) * 1e-3 # kg / mol

axfDIC = Axis(fig[5, 1], xlabel = "Time (days)", ylabel = "Flux (kgCO₂/m²/year)",
                         title = "Air-sea CO₂ flux and Sinking", limits = ((0, times[end] / days), nothing))
lines!(axfDIC, times / days, air_sea_CO₂_flux[1, 1, grid.Nz, :] / 1e3 * CO₂_molar_mass * year * 10, linewidth = 3, label = "Air-sea flux x10")
lines!(axfDIC, times / days, carbon_export / 1e3    * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export")
Legend(fig[5, 2], axfDIC, framevisible = false)

fig

This page was generated using Literate.jl.