Exercise 12: Seasonal diffusivity

This exercise uses a seasonal mixed-layer-depth diffusivity profile in a simple 1D water-column model. The physical model setup is based on an example provided in the OceanBioME.jl documentation and represents an idealized 200m deep North Atlantic time series.

Loading dependencies

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

using Agate
using Agate.Introspection: tracer_groups
using Agate.Library.Light
using OceanBioME
using OceanBioME: Biogeochemistry
using Oceananigans
using Oceananigans.Units
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)

const year = years = 365day

Ecosystem model

First, we construct our ecosystem model. Here, we use a default 2 phytoplankton, 2 zooplankton Agate.jl-NiPiZD ecosystem model.

bgc = Agate.Models.NiPiZD.construct()
groups = tracer_groups(bgc)

Forcings

Second, we define the model physical forcings. Mixed layer depth (MLD) forces the physical mixing (diffusivity), while PAR is held at a fixed surface value with depth-dependent attenuation.

#diffusivity
@inline function diffusivity(x, y, z, t)
    H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0)
    function fmld1(t)
        return H(t, 50days, year) *
               (1 / (1 + exp(-(t - 100days) / 5days))) *
               (1 / (1 + exp((t - 330days) / 25days)))
    end
    function MLD(t)
        return -(
            10 +
            340 * (
                1 - fmld1(year - eps(year)) * exp(-mod(t, year) / 25days) -
                fmld1(mod(t, year))
            )
        )
    end
    return 1e-2 * (1 + tanh((z - MLD(t)) / 10)) / 2 + 1e-4
end

#irradiance
@inline function constant_PAR(x, y, z, t)
    PAR⁰ = 80
    attenuation = 0.04
    return PAR⁰ * exp(attenuation * z)
end

#plots
t_range = 0.0:days:(365.0 * days)  # Time range from 0 to 365 days
z_range = -200.0:10.0:0.0  # Depth range from -200m to 0m
x, y, z = 0.0, 0.0, 0.0
κₜ_values = [diffusivity(x, y, z, t) for t in t_range, z in z_range]
PAR_values = [constant_PAR(x, y, z, t) for t in t_range, z in z_range]

fig_forcing = Figure(; size=(800, 600), fontsize=14)
ax1 = Axis(fig_forcing[1, 1]; xlabel="Time (days)", ylabel="Depth (m)", title="irradiance")
hm1 = CairoMakie.heatmap!(ax1, t_range ./ days, z_range, PAR_values; colormap=:viridis)
Colorbar(fig_forcing[1, 2], hm1)

ax2 = Axis(fig_forcing[2, 1]; xlabel="Time (days)", ylabel="Depth (m)", title="diffusivity")
hm2 = CairoMakie.heatmap!(ax2, t_range ./ days, z_range, κₜ_values; colormap=:viridis)
Colorbar(fig_forcing[2, 2], hm2)

display(fig_forcing)
fig_forcing

Physical model

grid = RectilinearGrid(; size=(1, 1, 20), extent=(20meters, 20meters, 200meters))

bgc_model = Biogeochemistry(
    bgc; light_attenuation=FunctionFieldPAR(; grid, PAR_f=constant_PAR)
)

full_model = NonhydrostaticModel(;
    grid,
    clock=Clock(; time=0.0),
    timestepper=:QuasiAdamsBashforth2,
    closure=ScalarDiffusivity(
        VerticallyImplicitTimeDiscretization(); ν=diffusivity, κ=diffusivity
    ),
    biogeochemistry=bgc_model,
)

Initial conditions

set!(full_model; default_initial_conditions(bgc; detritus = 0.0, total_plankton_biomass = 0.12)...) # mmol N / m³

Simulation

filename = "12_diffusivity_seasonal.jld2"

simulation = Simulation(full_model; Δt=1hours, stop_time=1year)

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

run!(simulation)

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (11.123 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (17.410 seconds).
[ Info: Simulation is stopping after running for 45.064 seconds.
[ Info: Simulation time 365 days equals or exceeds stop time 365 days.

Plotting

#Load time series data
timeseries = NamedTuple{keys(full_model.tracers)}(
    FieldTimeSeries(filename, "$field") for field in keys(full_model.tracers)
)

#Use Agate's introspection helpers to recover the structural tracer layout
all_keys = [groups.plankton..., groups.nonplankton...]

#Create figure with appropriate size
fig = Figure(; size=(800, 1200), fontsize=16)

#Plot all fields
for (i, key) in enumerate(all_keys)
    x_nodes, y_nodes, z_nodes = nodes(timeseries[key])
    z_vals = collect(z_nodes)
    times = collect(timeseries[key].times / days)

    ax = Axis(
        fig[i, 1];
        title="$(key) concentration (mmol N / m³)",
        xlabel="Time (days)",
        ylabel="z (m)",
        limits=((0, 365), (-200, 0)),
    )
    hm = heatmap!(
        ax,
        times,
        z_vals,
        Float32.(interior(timeseries[key],1,1,:,:)');
        colormap=:viridis,
        rasterize=true,
    )  # Rasterize for smaller output
    Colorbar(fig[i, 2], hm)
end

#Save figure
save("12_diffusivity_seasonal.png", fig)

display(fig)
fig  # Display the figure