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.Fields: FunctionField
using Oceananigans.Units
using CairoMakie

Ecosystem model

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

phyto_size_structure = (n=3, min_esd=1, max_esd=10, spacing=:log)

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

zoo_size_structure = [10.0, 32.0, 100.0]

size_structure = (;
    phytoplankton=(P=phyto_size_structure,),
    zooplankton=(Z=zoo_size_structure,),
)

bgc = Agate.Models.NiPiZD.construct(; size_structure)

We can inspect the tracer names required by the model.

println(tracer_names(bgc))

[:N, :D, :P_1, :P_2, :P_3, :Z_1, :Z_2, :Z_3]

As well as the plankton sizes:

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

P_1: 1.0 μm
P_2: 3.1622776601683795 μm
P_3: 10.000000000000002 μm
Z_1: 10.0 μm
Z_2: 32.0 μm
Z_3: 100.0 μ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.

grid = BoxModelGrid()
clock = Clock(; time=zero(grid))
PAR = FunctionField{Center,Center,Center}(CyclicalPAR(-10), grid; clock)
light_attenuation = PrescribedPhotosyntheticallyActiveRadiation(PAR)

bgc_model = Biogeochemistry(bgc; light_attenuation)

full_model = BoxModel(; grid, clock, biogeochemistry=bgc_model)

Initial conditions

set!(full_model; N=8.0, P_1=0.01, P_2=0.05, P_3=0.1, Z_1=0.01, Z_2=0.01, Z_3=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...
┌ Warning: error ArgumentError("a group or dataset named grid is already present within this group") thrown when trying to serialize the grid at serialized/grid
└ @ Oceananigans.OutputWriters ~/.julia/packages/Oceananigans/9oDU5/src/OutputWriters/jld2_writer.jl:251
[ Info:     ... simulation initialization complete (6.639 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (881.174 ms).
[ Info: Simulation is stopping after running for 11.905 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.