Size structure example

This example changes the number of plankton and the size structure of the Agate.jl-NiPiZD model.

Loading dependencies

The example uses Agate.jl, Oceananigans.jl, and OceanBioME.jl for the ecosystem simulation. CairoMakie.jl is used for plotting.

using Agate
using Agate.Introspection: tracer_names, plankton_tracers, plankton_diameters
using Agate.Library.Light
using OceanBioME
using OceanBioME: Biogeochemistry
using Oceananigans
using Oceananigans.Units
using CairoMakie

Ecosystem model

The plankton size structure is configured separately for phytoplankton and zooplankton. Here, phytoplankton are split into three classes spanning 1 to 10 μm equivalent spherical diameter (ESD), using logarithmic spacing.

phyto_size_structure = (n=3, min_esd=1, max_esd=10, splitting=:log_splitting)

Alternatively, an explicit array can be used. This zooplankton size structure also defines three classes, specifying the ESD of each class directly.

zoo_size_structure = [10.0, 32.0, 100.0]

3-element Vector{Float64}:
  10.0
  32.0
 100.0

The model is constructed following the default model setup but with the custom size structures. This creates three phytoplankton tracers and three zooplankton tracers.

bgc = Agate.Models.NiPiZD.construct(; phyto_size_structure, zoo_size_structure)

We can inspect the tracer names required by the model.

println(tracer_names(bgc))

[:N, :D, :Z1, :Z2, :Z3, :P1, :P2, :P3]

As well as the plankton sizes:

for (name, diameter) in zip(plankton_tracers(bgc), plankton_diameters(bgc))
    println(name, ": ", diameter, " μm")
end

Z1: 10.0 μm
Z2: 32.0 μm
Z3: 100.0 μm
P1: 1.0 μm
P2: 3.1622776601683795 μm
P3: 10.000000000000002 μm

Physical model

Next, we define the light model and combine it with the Agate.jl ecosystem model using OceanBioME.jl. The BoxModelGrid represents a well-mixed zero-dimensional water column.

light_attenuation = FunctionFieldPAR(; grid=BoxModelGrid())

bgc_model = Biogeochemistry(bgc; light_attenuation)

full_model = BoxModel(; biogeochemistry=bgc_model)

Initial conditions

set!(full_model; N=8.0, P1=0.01, P2=0.05, P3=0.1, Z1=0.01, Z2=0.01, Z3=0.01, D=0.01)

Simulation

filename = "size_structure.jld2"

simulation = Simulation(full_model; Δt=240minutes, stop_time=1095days)

simulation.output_writers[:fields] = JLD2Writer(
    full_model,
    full_model.fields;
    filename,
    schedule=TimeInterval(1day),
    overwrite_existing=true,
)

run!(simulation)

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (2.418 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (822.959 ms).
[ Info: Simulation is stopping after running for 8.637 seconds.
[ Info: Simulation time 1095 days equals or exceeds stop time 1095 days.

Plotting

Load the simulated tracer time series from disk.

tracer_syms = tracer_names(bgc)

timeseries = (;
    (s => FieldTimeSeries(filename, string(s))[1, 1, 1, :] for s in tracer_syms)...
)

Plot each tracer concentration through time.

times = FieldTimeSeries(filename, string(first(tracer_syms))).times

fig = Figure(; size=(650, 520), fontsize=12)

for (idx, sym) in enumerate(tracer_syms)
    ax = Axis(
        fig[floor(Int, (idx - 1) / 2), Int((idx - 1) % 2)];
        ylabel=string(sym),
        xlabel="Days",
        title="$(sym) concentration (mmol N / m³)",
        titlesize=16,
        xlabelsize=12,
        ylabelsize=12,
        xticklabelsize=10,
        yticklabelsize=10,
    )
    lines!(ax, times / day, getproperty(timeseries, sym); linewidth=1.5)
end

fig

This page was generated using Literate.jl.