Nutrients, plankton, bacteria, detritus

This example illustrates how to use ClimaOceanBiogeochemistry's CarbonAlkalinityNutrients model in a single column context.

using ClimaOceanBiogeochemistry: CarbonAlkalinityNutrients

using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity

using CUDA

using Printf
using CairoMakie

A single column grid

We set up a single column grid with 4 m grid spacing that's 256 m deep:

grid = RectilinearGrid(#GPU(),
                       size = 64,
                       z = (-256meters, 0),
                       topology = (Flat, Flat, Bounded))

1×1×64 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 0×0×3 halo
├── Flat x
├── Flat y
└── Bounded  z ∈ [-256.0, 0.0]    regularly spaced with Δz=4.0

Convection then quiet

To illustrate the dynamics of CarbonAlkalinityNutrients, we set up a physical scenario in which strong convection drives turbulent mixing for 4 days, and then abruptly shuts off. Once the convective turbulence dies down, plankton start to grow.

Qᵇ(t) = ifelse(t < 4days, 1e-7, 0.0)
b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ))

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: ContinuousBoundaryFunction Qᵇ at (Nothing, Nothing, Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

We put the pieces together. The important line here is biogeochemistry = CarbonAlkalinityNutrients(). We use all default parameters.

model = HydrostaticFreeSurfaceModel(; grid,
                                    biogeochemistry = CarbonAlkalinityNutrients(),
                                    closure = CATKEVerticalDiffusivity(),
                                    tracers = (:b, :e),
                                    buoyancy = BuoyancyTracer(),
                                    boundary_conditions = (; b=b_bcs))

HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×64 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 0×0×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: (b, e, DIC, Alk, PO₄, NO₃, DOP, Fe)
├── closure: CATKEVerticalDiffusivity{VerticallyImplicitTimeDiscretization}
├── buoyancy: Oceananigans.BuoyancyModels.BuoyancyTracer with ĝ = NegativeZDirection()
├── advection scheme: 
│   ├── momentum: Centered reconstruction order 2
│   ├── b: Centered reconstruction order 2
│   ├── e: Centered reconstruction order 2
│   ├── DIC: Centered reconstruction order 2
│   ├── Alk: Centered reconstruction order 2
│   ├── PO₄: Centered reconstruction order 2
│   ├── NO₃: Centered reconstruction order 2
│   ├── DOP: Centered reconstruction order 2
│   └── Fe: Centered reconstruction order 2
└── coriolis: Nothing

Initial conditions

We initialize the model with reasonable nutrients, detritus, and a nutrient mixed layer.

N₀ = 1e-1 # Surface nutrient concentration
D₀ = 1e-1 # Surface detritus concentration
dᴺ = 50.0 # Nutrient mixed layer depth
N² = 1e-5 # Buoyancy gradient, s⁻²

bᵢ(z) = N² * z
Nᵢ(z) = N₀ * max(1, exp(-(z + dᴺ) / 100))
Dᵢ(z) = D₀ * exp(z / 10)

Dᵢ (generic function with 1 method)

set!(model, b=bᵢ, P=1e-1, B=1e-1, D=Dᵢ, N=Nᵢ, e=1e-6)

set!(model, b=bᵢ, e=1e-6)

A simulation of physical-biological interaction

We construct a simple simulation that emits a message every 10 iterations and outputs tracer fields.

simulation = Simulation(model, Δt=10minutes, stop_time=12days)

progress(sim) = @printf("Iteration: %d, time: %s\n", iteration(sim), prettytime(sim))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

filename = "single_column_carbon_alkalinity_nutrients.jld2"

simulation.output_writers[:fields] = JLD2OutputWriter(model, model.tracers;
                                                      filename,
                                                      schedule = TimeInterval(20minutes),
                                                      overwrite_existing = true)

run!(simulation)

[ Info: Initializing simulation...
Iteration: 0, time: 0 seconds
[ Info:     ... simulation initialization complete (1.957 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (30.875 seconds).
Iteration: 10, time: 1.667 hours
Iteration: 20, time: 3.333 hours
Iteration: 30, time: 5 hours
Iteration: 40, time: 6.667 hours
Iteration: 50, time: 8.333 hours
Iteration: 60, time: 10 hours
Iteration: 70, time: 11.667 hours
Iteration: 80, time: 13.333 hours
Iteration: 90, time: 15 hours
Iteration: 100, time: 16.667 hours
Iteration: 110, time: 18.333 hours
Iteration: 120, time: 20 hours
Iteration: 130, time: 21.667 hours
Iteration: 140, time: 23.333 hours
Iteration: 150, time: 1.042 days
Iteration: 160, time: 1.111 days
Iteration: 170, time: 1.181 days
Iteration: 180, time: 1.250 days
Iteration: 190, time: 1.319 days
Iteration: 200, time: 1.389 days
Iteration: 210, time: 1.458 days
Iteration: 220, time: 1.528 days
Iteration: 230, time: 1.597 days
Iteration: 240, time: 1.667 days
Iteration: 250, time: 1.736 days
Iteration: 260, time: 1.806 days
Iteration: 270, time: 1.875 days
Iteration: 280, time: 1.944 days
Iteration: 290, time: 2.014 days
Iteration: 300, time: 2.083 days
Iteration: 310, time: 2.153 days
Iteration: 320, time: 2.222 days
Iteration: 330, time: 2.292 days
Iteration: 340, time: 2.361 days
Iteration: 350, time: 2.431 days
Iteration: 360, time: 2.500 days
Iteration: 370, time: 2.569 days
Iteration: 380, time: 2.639 days
Iteration: 390, time: 2.708 days
Iteration: 400, time: 2.778 days
Iteration: 410, time: 2.847 days
Iteration: 420, time: 2.917 days
Iteration: 430, time: 2.986 days
Iteration: 440, time: 3.056 days
Iteration: 450, time: 3.125 days
Iteration: 460, time: 3.194 days
Iteration: 470, time: 3.264 days
Iteration: 480, time: 3.333 days
Iteration: 490, time: 3.403 days
Iteration: 500, time: 3.472 days
Iteration: 510, time: 3.542 days
Iteration: 520, time: 3.611 days
Iteration: 530, time: 3.681 days
Iteration: 540, time: 3.750 days
Iteration: 550, time: 3.819 days
Iteration: 560, time: 3.889 days
Iteration: 570, time: 3.958 days
Iteration: 580, time: 4.028 days
Iteration: 590, time: 4.097 days
Iteration: 600, time: 4.167 days
Iteration: 610, time: 4.236 days
Iteration: 620, time: 4.306 days
Iteration: 630, time: 4.375 days
Iteration: 640, time: 4.444 days
Iteration: 650, time: 4.514 days
Iteration: 660, time: 4.583 days
Iteration: 670, time: 4.653 days
Iteration: 680, time: 4.722 days
Iteration: 690, time: 4.792 days
Iteration: 700, time: 4.861 days
Iteration: 710, time: 4.931 days
Iteration: 720, time: 5 days
Iteration: 730, time: 5.069 days
Iteration: 740, time: 5.139 days
Iteration: 750, time: 5.208 days
Iteration: 760, time: 5.278 days
Iteration: 770, time: 5.347 days
Iteration: 780, time: 5.417 days
Iteration: 790, time: 5.486 days
Iteration: 800, time: 5.556 days
Iteration: 810, time: 5.625 days
Iteration: 820, time: 5.694 days
Iteration: 830, time: 5.764 days
Iteration: 840, time: 5.833 days
Iteration: 850, time: 5.903 days
Iteration: 860, time: 5.972 days
Iteration: 870, time: 6.042 days
Iteration: 880, time: 6.111 days
Iteration: 890, time: 6.181 days
Iteration: 900, time: 6.250 days
Iteration: 910, time: 6.319 days
Iteration: 920, time: 6.389 days
Iteration: 930, time: 6.458 days
Iteration: 940, time: 6.528 days
Iteration: 950, time: 6.597 days
Iteration: 960, time: 6.667 days
Iteration: 970, time: 6.736 days
Iteration: 980, time: 6.806 days
Iteration: 990, time: 6.875 days
Iteration: 1000, time: 6.944 days
Iteration: 1010, time: 7.014 days
Iteration: 1020, time: 7.083 days
Iteration: 1030, time: 7.153 days
Iteration: 1040, time: 7.222 days
Iteration: 1050, time: 7.292 days
Iteration: 1060, time: 7.361 days
Iteration: 1070, time: 7.431 days
Iteration: 1080, time: 7.500 days
Iteration: 1090, time: 7.569 days
Iteration: 1100, time: 7.639 days
Iteration: 1110, time: 7.708 days
Iteration: 1120, time: 7.778 days
Iteration: 1130, time: 7.847 days
Iteration: 1140, time: 7.917 days
Iteration: 1150, time: 7.986 days
Iteration: 1160, time: 8.056 days
Iteration: 1170, time: 8.125 days
Iteration: 1180, time: 8.194 days
Iteration: 1190, time: 8.264 days
Iteration: 1200, time: 8.333 days
Iteration: 1210, time: 8.403 days
Iteration: 1220, time: 8.472 days
Iteration: 1230, time: 8.542 days
Iteration: 1240, time: 8.611 days
Iteration: 1250, time: 8.681 days
Iteration: 1260, time: 8.750 days
Iteration: 1270, time: 8.819 days
Iteration: 1280, time: 8.889 days
Iteration: 1290, time: 8.958 days
Iteration: 1300, time: 9.028 days
Iteration: 1310, time: 9.097 days
Iteration: 1320, time: 9.167 days
Iteration: 1330, time: 9.236 days
Iteration: 1340, time: 9.306 days
Iteration: 1350, time: 9.375 days
Iteration: 1360, time: 9.444 days
Iteration: 1370, time: 9.514 days
Iteration: 1380, time: 9.583 days
Iteration: 1390, time: 9.653 days
Iteration: 1400, time: 9.722 days
Iteration: 1410, time: 9.792 days
Iteration: 1420, time: 9.861 days
Iteration: 1430, time: 9.931 days
Iteration: 1440, time: 10 days
Iteration: 1450, time: 10.069 days
Iteration: 1460, time: 10.139 days
Iteration: 1470, time: 10.208 days
Iteration: 1480, time: 10.278 days
Iteration: 1490, time: 10.347 days
Iteration: 1500, time: 10.417 days
Iteration: 1510, time: 10.486 days
Iteration: 1520, time: 10.556 days
Iteration: 1530, time: 10.625 days
Iteration: 1540, time: 10.694 days
Iteration: 1550, time: 10.764 days
Iteration: 1560, time: 10.833 days
Iteration: 1570, time: 10.903 days
Iteration: 1580, time: 10.972 days
Iteration: 1590, time: 11.042 days
Iteration: 1600, time: 11.111 days
Iteration: 1610, time: 11.181 days
Iteration: 1620, time: 11.250 days
Iteration: 1630, time: 11.319 days
Iteration: 1640, time: 11.389 days
Iteration: 1650, time: 11.458 days
Iteration: 1660, time: 11.528 days
Iteration: 1670, time: 11.597 days
Iteration: 1680, time: 11.667 days
Iteration: 1690, time: 11.736 days
Iteration: 1700, time: 11.806 days
Iteration: 1710, time: 11.875 days
Iteration: 1720, time: 11.944 days
[ Info: Simulation is stopping after running for 39.149 seconds.
[ Info: Simulation time 12 days equals or exceeds stop time 12 days.

Visualization

All that's left is to visualize the results.

bt   = FieldTimeSeries(filename, "b")
et   = FieldTimeSeries(filename, "e")
DICt = FieldTimeSeries(filename, "DIC")
Alkt = FieldTimeSeries(filename, "Alk")
PO₄t = FieldTimeSeries(filename, "PO₄")
NO₃t = FieldTimeSeries(filename, "NO₃")
DOPt = FieldTimeSeries(filename, "DOP")
Fet  = FieldTimeSeries(filename, "Fe")

t = bt.times
nt = length(t)
z = znodes(bt)

fig = Figure(size=(1200, 600))

axb = Axis(fig[1, 1], xlabel="Buoyancy (m² s⁻³)", ylabel="z (m)")
axe = Axis(fig[1, 2], xlabel="Turbulent kinetic energy (m² s²)")
axP = Axis(fig[1, 3], xlabel="Concentration (mmol)")
axN = Axis(fig[1, 4], xlabel="Nutrient concentration (mmol)")

xlims!(axe, -1e-5, 1e-3)
xlims!(axP, 0, 0.2)

slider = Slider(fig[2, 1:4], range=1:nt, startvalue=1)
n = slider.value

title = @lift @sprintf("Convecting plankton at t = %d hours", t[$n] / hour)
Label(fig[0, 1:4], title)

bn = @lift interior(bt[$n], 1, 1, :)
en = @lift interior(et[$n], 1, 1, :)
DICn = @lift interior(DICt[$n], 1, 1, :)
Alkn = @lift interior(Alkt[$n], 1, 1, :)
PO₄n = @lift interior(PO₄t[$n], 1, 1, :)
NO₃n = @lift interior(NO₃t[$n], 1, 1, :)
DOPn = @lift interior(DOPt[$n], 1, 1, :)
Fen  = @lift interior(Fet[$n], 1, 1, :)

lines!(axb, bn, z)
lines!(axe, en, z)

lines!(axP, DICn, z, label="DIC")
lines!(axP, PO₄n, z, label="Phosphate")
lines!(axP, Fen, z, label="Iron")
axislegend(axP)

lines!(axN, DOPn, z)

record(fig, "carbon_alkalinity_nutrients.mp4", 1:nt, framerate=24) do nn
    n[] = nn
end

This page was generated using Literate.jl.