Exercise 04: Number size classes

This exercise changes the number of phytoplankton and zooplankton size classes in the Agate.jl NiPiZD model. The total initial plankton biomass is held fixed and split evenly across the available plankton tracers, so the comparison isolates the effect of resolving more size classes.

Loading dependencies

using Agate
using Agate.Introspection: tracer_names
using CairoMakie

workshop_script = let dir = @__DIR__
    while !isfile(joinpath(dir, "src", "AgateWorkshop.jl"))
        parent = dirname(dir)
        parent == dir && error("Could not find src/AgateWorkshop.jl")
        dir = parent
    end
    joinpath(dir, "src", "AgateWorkshop.jl")
end
include(workshop_script)

mkpath("outputs")
mkpath("figures")

Number-of-classes cases

The default model has two phytoplankton classes and two zooplankton classes. The alternative cases keep the same broad phytoplankton and zooplankton size ranges and logarithmic spacing, but increase both groups to five and ten classes.

default_phyto_size_structure = (n = 2, min_esd = 2, max_esd = 10, splitting = :log_splitting)
default_zoo_size_structure = (n = 2, min_esd = 20, max_esd = 100, splitting = :linear_splitting)

function construct_size_class_bgc(n)
    return Agate.Models.NiPiZD.construct(;
        phyto_size_structure = (; default_phyto_size_structure..., n),
        zoo_size_structure = (; default_zoo_size_structure..., n),
    )
end

bgc_default = Agate.Models.NiPiZD.construct(;
    phyto_size_structure = default_phyto_size_structure,
    zoo_size_structure = default_zoo_size_structure,
)
bgc_5_each = construct_size_class_bgc(5)
bgc_10_each = construct_size_class_bgc(10)

bgcs = [bgc_default, bgc_5_each, bgc_10_each]
case_labels = ["default: 2 P, 2 Z", "5 P, 5 Z", "10 P, 10 Z"]

for (label, bgc) in zip(case_labels, bgcs)
    println(label)
    println(tracer_names(bgc))
end

default: 2 P, 2 Z
[:N, :D, :Z1, :Z2, :P1, :P2]
5 P, 5 Z
[:N, :D, :Z1, :Z2, :Z3, :Z4, :Z5, :P1, :P2, :P3, :P4, :P5]
10 P, 10 Z
[:N, :D, :Z1, :Z2, :Z3, :Z4, :Z5, :Z6, :Z7, :Z8, :Z9, :Z10, :P1, :P2, :P3, :P4, :P5, :P6, :P7, :P8, :P9, :P10]

Run zero-dimensional ecosystem simulations

default_initial_conditions distributes the same total plankton biomass evenly across whichever plankton tracers are present in each configuration.

runs = [
    run_box_model(
        bgc;
        filename = joinpath("outputs", "04_number_size_classes_$(i).jld2"),
        initial_conditions = default_initial_conditions(bgc; nutrient = 8.0, total_plankton_biomass = 0.15),
    )
    for (i, bgc) in enumerate(bgcs)
]

timeseries = [read_box_tracer_timeseries(run.filename, run.tracer_syms) for run in runs]

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (1.774 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (162.155 μs).
[ Info: Simulation is stopping after running for 4.403 seconds.
[ Info: Simulation time 1095 days equals or exceeds stop time 1095 days.
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (1.219 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.454 seconds).
[ Info: Simulation is stopping after running for 9.749 seconds.
[ Info: Simulation time 1095 days equals or exceeds stop time 1095 days.
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (2.102 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.357 seconds).
[ Info: Simulation is stopping after running for 19.299 seconds.
[ Info: Simulation time 1095 days equals or exceeds stop time 1095 days.

Tracer concentration comparison

This diagnostic compares the tracer set across the three configurations. Tracers that exist only in the higher-resolution cases appear only for those cases.

comparison_variables = sort!(unique(vcat([collect(keys(ts.data)) for ts in timeseries]...)); by = string)
comparison_figure_path = joinpath("figures", "04_number_size_classes_timeseries_comparison.png")
fig_comparison = plot_box_timeseries(timeseries; labels = case_labels, variables = comparison_variables)
save(comparison_figure_path, fig_comparison; px_per_unit = 1)
display(fig_comparison)
fig_comparison

Relative nitrogen contributions: default size classes

nitrogen_default_figure_path = joinpath("figures", "04_number_size_classes_relative_nitrogen_default.png")
fig_nitrogen_default = plot_contributions(timeseries[1].times, timeseries[1].data)
save(nitrogen_default_figure_path, fig_nitrogen_default; px_per_unit = 1)
display(fig_nitrogen_default)
fig_nitrogen_default

Relative nitrogen contributions: five classes each

nitrogen_5_figure_path = joinpath("figures", "04_number_size_classes_relative_nitrogen_5_each.png")
fig_nitrogen_5 = plot_contributions(timeseries[2].times, timeseries[2].data)
save(nitrogen_5_figure_path, fig_nitrogen_5; px_per_unit = 1)
display(fig_nitrogen_5)
fig_nitrogen_5

Relative nitrogen contributions: ten classes each

nitrogen_10_figure_path = joinpath("figures", "04_number_size_classes_relative_nitrogen_10_each.png")
fig_nitrogen_10 = plot_contributions(timeseries[3].times, timeseries[3].data)
save(nitrogen_10_figure_path, fig_nitrogen_10; px_per_unit = 1)
display(fig_nitrogen_10)
fig_nitrogen_10

Community-weighted mean size comparison

size_comparison_figure_path = joinpath("figures", "04_number_size_classes_cwm_size_comparison.png")
fig_size_comparison = plot_cwm_size(timeseries, bgcs; labels = case_labels)
save(size_comparison_figure_path, fig_size_comparison; px_per_unit = 1)
display(fig_size_comparison)
fig_size_comparison