Plankton mixing and blooming

In this example, we simulate the mixing of phytoplankton by convection that decreases in time and eventually shuts off, thereby precipitating a phytoplankton bloom. A similar scenario was simulated by Taylor and Ferrari (2011), providing evidence that the "critical turbulence hypothesis" explains the explosive bloom of oceanic phytoplankton observed in spring.

The phytoplankton in our model are advected, diffuse, grow, and die according to

\[∂_t P + \bm{u ⋅ ∇} P - κ ∇²P = (μ₀ \exp(z / λ) - m) \, P \, ,\]

where $\bm{u}$ is the turbulent velocity field, $κ$ is an isotropic diffusivity, $μ₀$ is the phytoplankton growth rate at the surface, $λ$ is the scale over which sunlight attenuates away from the surface, and $m$ is the mortality rate of phytoplankton due to viruses and grazing by zooplankton. We use Oceananigans' Forcing abstraction to implement the phytoplankton dynamics described by the right side of the phytoplankton equation above.

This example demonstrates

  • How to use a user-defined forcing function to simulate the dynamics of phytoplankton growth in sunlight and grazing by zooplankton.
  • How to set time-dependent boundary conditions.
  • How to use the TimeStepWizard to adapt the simulation time-step.
  • How to use AveragedField to diagnose spatial averages of model fields.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, Plots, JLD2, Measures"

The grid

We use a two-dimensional grid with 64² points and 1 m grid spacing:

using Oceananigans

grid = RegularRectilinearGrid(size=(64, 1, 64), extent=(64, 1, 64))
RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 64.0], y ∈ [0.0, 1.0], z ∈ [-64.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (64, 1, 64)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (1.0, 1.0, 1.0)

Boundary conditions

We impose a surface buoyancy flux that's initially constant and then decays to zero,

using Oceananigans.Utils

buoyancy_flux(x, y, t, p) = p.initial_buoyancy_flux * exp(-t^4 / (24 * p.shut_off_time^4))

buoyancy_flux_parameters = (initial_buoyancy_flux = 1e-8, # m² s⁻³
                                    shut_off_time = 2hours)

buoyancy_flux_bc = BoundaryCondition(Flux, buoyancy_flux, parameters = buoyancy_flux_parameters)
BoundaryCondition: type=Flux, condition=buoyancy_flux(x, y, t, p) in Main.ex-convecting_plankton at none:2

The fourth power in the argument of exp above helps keep the buoyancy flux relatively constant during the first phase of the simulation. We produce a plot of this time-dependent buoyancy flux for the visually-oriented,

using Plots, Measures

time = range(0, 12hours, length=100)

flux_plot = plot(time ./ hour, [buoyancy_flux(0, 0, t, buoyancy_flux_parameters) for t in time],
                 linewidth = 2, xlabel = "Time (hours)", ylabel = "Surface buoyancy flux (m² s⁻³)",
                 size = (800, 300), margin = 5mm, label = nothing)

The buoyancy flux effectively shuts off after 6 hours of simulation time.

The flux convention in Oceananigans.jl

Fluxes are defined by the direction a quantity is carried: positive velocities produce positive fluxes, while negative velocities produce negative fluxes. Diffusive fluxes are defined with the same convention. A positive flux at the top boundary transports buoyancy upwards, out of the domain. This means that a positive flux of buoyancy at the top boundary reduces the buoyancy of near-surface fluid, causing convection.

The initial condition and bottom boundary condition impose the constant buoyancy gradient

N² = 1e-4 # s⁻²

buoyancy_gradient_bc = BoundaryCondition(Gradient, N²)
BoundaryCondition: type=Gradient, condition=0.0001

In summary, the buoyancy boundary conditions impose a destabilizing flux at the top and a stable buoyancy gradient at the bottom:

buoyancy_bcs = TracerBoundaryConditions(grid, top = buoyancy_flux_bc, bottom = buoyancy_gradient_bc)
Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions
├── x: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
├── y: CoordinateBoundaryConditions{BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing},BoundaryCondition{Oceananigans.BoundaryConditions.Periodic,Nothing}}
└── z: CoordinateBoundaryConditions{BoundaryCondition{Gradient,Float64},BoundaryCondition{Flux,Oceananigans.BoundaryConditions.ContinuousBoundaryFunction{Nothing,Nothing,Nothing,Nothing,typeof(Main.ex-convecting_plankton.buoyancy_flux),NamedTuple{(:initial_buoyancy_flux, :shut_off_time),Tuple{Float64,Float64}},Tuple{},Nothing,Nothing}}}

Phytoplankton dynamics: light-dependent growth and uniform mortality

We use a simple model for the growth of phytoplankton in sunlight and decay due to viruses and grazing by zooplankton,

growing_and_grazing(x, y, z, t, P, p) = (p.μ₀ * exp(z / p.λ) - p.m) * P

with parameters

plankton_dynamics_parameters = (μ₀ = 1/day,   # surface growth rate
                                 λ = 5,       # sunlight attenuation length scale (m)
                                 m = 0.1/day) # mortality rate due to virus and zooplankton grazing
(μ₀ = 1.1574074074074073e-5, λ = 5, m = 1.1574074074074074e-6)

We tell Forcing that our plankton model depends on the plankton concentration P and the chosen parameters,

plankton_dynamics = Forcing(growing_and_grazing, field_dependencies = :P,
                            parameters = plankton_dynamics_parameters)
ContinuousForcing{NamedTuple{(:μ₀, :λ, :m),Tuple{Float64,Int64,Float64}}}
├── func: growing_and_grazing
├── parameters: (μ₀ = 1.1574074074074073e-5, λ = 5, m = 1.1574074074074074e-6)
└── field dependencies: (:P,)

The model

The name "P" for phytoplankton is specified in the constructor for IncompressibleModel. We additionally specify a fifth-order advection scheme, third-order Runge-Kutta time-stepping, isotropic viscosity and diffusivities, and Coriolis forces appropriate for planktonic convection at mid-latitudes on Earth.

using Oceananigans.Advection

model = IncompressibleModel(
                   grid = grid,
              advection = UpwindBiasedFifthOrder(),
            timestepper = :RungeKutta3,
                closure = IsotropicDiffusivity(ν=1e-4, κ=1e-4),
               coriolis = FPlane(f=1e-4),
                tracers = (:b, :P), # P for Plankton
               buoyancy = BuoyancyTracer(),
                forcing = (P=plankton_dynamics,),
    boundary_conditions = (b=buoyancy_bcs,)
)
IncompressibleModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=64, Ny=1, Nz=64)
├── tracers: (:b, :P)
├── closure: IsotropicDiffusivity{Float64,NamedTuple{(:b, :P),Tuple{Float64,Float64}}}
├── buoyancy: BuoyancyTracer
└── coriolis: FPlane{Float64}

Initial condition

We set the initial phytoplankton at $P = 1 \, \rm{μM}$. For buoyancy, we use a stratification that's mixed near the surface and linearly stratified below, superposed with surface-concentrated random noise.

mixed_layer_depth = 32 # m

stratification(z) = z < -mixed_layer_depth ? N² * z : - N² * mixed_layer_depth

noise(z) = 1e-4 * N² * grid.Lz * randn() * exp(z / 4)

initial_buoyancy(x, y, z) = stratification(z) + noise(z)

set!(model, b=initial_buoyancy, P=1)

Adaptive time-stepping, logging, output and simulation setup

We use a TimeStepWizard that limits the time-step to 2 minutes, and adapts the time-step such that CFL (Courant-Freidrichs-Lewy) number hovers around 1.0,

wizard = TimeStepWizard(cfl=1.0, Δt=2minutes, max_change=1.1, max_Δt=2minutes)
TimeStepWizard{Float64}(1.0, Inf, 1.1, 0.5, 120.0, 0.0, 120.0)

We also write a function that prints the progress of the simulation

using Printf

progress(sim) = @printf("Iteration: %d, time: %s, Δt: %s\n",
                        sim.model.clock.iteration,
                        prettytime(sim.model.clock.time),
                        prettytime(sim.Δt.Δt))

simulation = Simulation(model, Δt=wizard, stop_time=24hour,
                        iteration_interval=20, progress=progress)
Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0 
├── Next time step (TimeStepWizard{Float64}): 2 minutes 
├── Iteration interval: 20
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0 seconds, wall time limit: Inf
├── Stop time: 1 day, stop iteration: Inf
├── Diagnostics: OrderedCollections.OrderedDict with 1 entry:
│   └── nan_checker => NaNChecker
└── Output writers: OrderedCollections.OrderedDict with no entries

We add a basic JLD2OutputWriter that writes velocities and both the two-dimensional and horizontally-averaged plankton concentration,

using Oceananigans.OutputWriters, Oceananigans.Fields

averaged_plankton = AveragedField(model.tracers.P, dims=(1, 2))

outputs = (w = model.velocities.w,
           plankton = model.tracers.P,
           averaged_plankton = averaged_plankton)

simulation.output_writers[:simple_output] =
    JLD2OutputWriter(model, outputs,
                     schedule = TimeInterval(20minutes),
                     prefix = "convecting_plankton",
                     force = true)
JLD2OutputWriter scheduled on TimeInterval(20 minutes):
├── filepath: ./convecting_plankton.jld2
├── 3 outputs: (:w, :plankton, :averaged_plankton)
├── field slicer: FieldSlicer(:, :, :, with_halos=false)
├── array type: Array{Float32}
├── including: [:grid, :coriolis, :buoyancy, :closure]
└── max filesize: Inf YiB
Using multiple output writers

Because each output writer is associated with a single output schedule, it often makes sense to use different output writers for different types of output. For example, reduced fields like AveragedField usually consume less disk space than two- or three-dimensional fields, and can thus be output more frequently without blowing up your hard drive. An arbitrary number of output writers may be added to simulation.output_writers.

The simulation is set up. Let there be plankton:

run!(simulation)
Iteration: 20, time: 40 minutes, Δt: 2 minutes
Iteration: 40, time: 1.333 hours, Δt: 2 minutes
Iteration: 60, time: 2 hours, Δt: 2 minutes
Iteration: 80, time: 2.372 hours, Δt: 1.155 minutes
Iteration: 100, time: 2.766 hours, Δt: 1.193 minutes
Iteration: 120, time: 3.133 hours, Δt: 1.138 minutes
Iteration: 140, time: 3.478 hours, Δt: 1.082 minutes
Iteration: 160, time: 3.821 hours, Δt: 1.032 minutes
Iteration: 180, time: 4.151 hours, Δt: 1.010 minutes
Iteration: 200, time: 4.457 hours, Δt: 55.702 seconds
Iteration: 220, time: 4.761 hours, Δt: 56.389 seconds
Iteration: 240, time: 5.081 hours, Δt: 58.239 seconds
Iteration: 260, time: 5.401 hours, Δt: 1.008 minutes
Iteration: 280, time: 5.759 hours, Δt: 1.109 minutes
Iteration: 300, time: 6.045 hours, Δt: 53.557 seconds
Iteration: 320, time: 6.349 hours, Δt: 54.947 seconds
Iteration: 340, time: 6.667 hours, Δt: 59.537 seconds
Iteration: 360, time: 6.973 hours, Δt: 55.100 seconds
Iteration: 380, time: 7.231 hours, Δt: 48.832 seconds
Iteration: 400, time: 7.496 hours, Δt: 48.856 seconds
Iteration: 420, time: 7.766 hours, Δt: 50.927 seconds
Iteration: 440, time: 8.062 hours, Δt: 56.020 seconds
Iteration: 460, time: 8.365 hours, Δt: 56.252 seconds
Iteration: 480, time: 8.683 hours, Δt: 58.215 seconds
Iteration: 500, time: 8.989 hours, Δt: 55.094 seconds
Iteration: 520, time: 9.309 hours, Δt: 58.572 seconds
Iteration: 540, time: 9.645 hours, Δt: 1.039 minutes
Iteration: 560, time: 9.971 hours, Δt: 1.016 minutes
Iteration: 580, time: 10.298 hours, Δt: 59.518 seconds
Iteration: 600, time: 10.617 hours, Δt: 1.001 minutes
Iteration: 620, time: 10.952 hours, Δt: 1.005 minutes
Iteration: 640, time: 11.310 hours, Δt: 1.096 minutes
Iteration: 660, time: 11.687 hours, Δt: 1.205 minutes
Iteration: 680, time: 12.080 hours, Δt: 1.206 minutes
Iteration: 700, time: 12.507 hours, Δt: 1.301 minutes
Iteration: 720, time: 12.977 hours, Δt: 1.431 minutes
Iteration: 740, time: 13.401 hours, Δt: 1.350 minutes
Iteration: 760, time: 13.848 hours, Δt: 1.359 minutes
Iteration: 780, time: 14.287 hours, Δt: 1.324 minutes
Iteration: 800, time: 14.632 hours, Δt: 1.053 minutes
Iteration: 820, time: 15 hours, Δt: 1.147 minutes
Iteration: 840, time: 15.417 hours, Δt: 1.262 minutes
Iteration: 860, time: 15.804 hours, Δt: 1.175 minutes
Iteration: 880, time: 16.207 hours, Δt: 1.240 minutes
Iteration: 900, time: 16.587 hours, Δt: 1.172 minutes
Iteration: 920, time: 16.987 hours, Δt: 1.202 minutes
Iteration: 940, time: 17.399 hours, Δt: 1.322 minutes
Iteration: 960, time: 17.856 hours, Δt: 1.419 minutes
Iteration: 980, time: 18.359 hours, Δt: 1.561 minutes
Iteration: 1000, time: 18.882 hours, Δt: 1.616 minutes
Iteration: 1020, time: 19.359 hours, Δt: 1.531 minutes
Iteration: 1040, time: 19.872 hours, Δt: 1.540 minutes
Iteration: 1060, time: 20.333 hours, Δt: 1.520 minutes
Iteration: 1080, time: 20.814 hours, Δt: 1.469 minutes
Iteration: 1100, time: 21.294 hours, Δt: 1.471 minutes
Iteration: 1120, time: 21.763 hours, Δt: 1.447 minutes
Iteration: 1140, time: 22.242 hours, Δt: 1.450 minutes
Iteration: 1160, time: 22.717 hours, Δt: 1.513 minutes
Iteration: 1180, time: 23.203 hours, Δt: 1.525 minutes
Iteration: 1200, time: 23.667 hours, Δt: 1.531 minutes
Iteration: 1213, time: 1 day, Δt: 1.635 minutes
[ Info: Simulation is stopping. Model time 1 day has hit or exceeded simulation stop time 1 day.

Notice how the time-step is reduced at early times, when turbulence is strong, and increases again towards the end of the simulation when turbulence fades.

Visualizing the solution

We'd like to a make a plankton movie. First we load the output file and build a time-series of the buoyancy flux,

using JLD2

file = jldopen(simulation.output_writers[:simple_output].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

times = [file["timeseries/t/$iter"] for iter in iterations]

buoyancy_flux_time_series = [buoyancy_flux(0, 0, t, buoyancy_flux_parameters) for t in times]

and then we construct the $x, z$ grid,

using Oceananigans.Grids: nodes

xw, yw, zw = nodes(model.velocities.w)
xp, yp, zp = nodes(model.tracers.P)

Finally, we animate plankton mixing and blooming,

using Plots

@info "Making a movie about plankton..."

w_lim = 0   # the maximum(abs(w)) across the whole timeseries

for (i, iteration) in enumerate(iterations)
    w = file["timeseries/w/$iteration"][:, 1, :]

    global w_lim = maximum([w_lim, maximum(abs.(w))])
end

anim = @animate for (i, iteration) in enumerate(iterations)

    @info "Plotting frame $i from iteration $iteration..."

    t = file["timeseries/t/$iteration"]
    w = file["timeseries/w/$iteration"][:, 1, :]
    P = file["timeseries/plankton/$iteration"][:, 1, :]
    averaged_P = file["timeseries/averaged_plankton/$iteration"][1, 1, :]

    P_min = minimum(P) - 1e-9
    P_max = maximum(P) + 1e-9
    P_lims = (0.95, 1.1)

    w_levels = range(-w_lim, stop=w_lim, length=20)

    P_levels = collect(range(P_lims[1], stop=P_lims[2], length=20))
    P_lims[1] > P_min && pushfirst!(P_levels, P_min)
    P_lims[2] < P_max && push!(P_levels, P_max)

    kwargs = (xlabel="x (m)", ylabel="y (m)", aspectratio=1, linewidth=0, colorbar=true,
              xlims=(0, model.grid.Lx), ylims=(-model.grid.Lz, 0))

    w_contours = contourf(xw, zw, w';
                          color = :balance,
                          levels = w_levels,
                          clims = (-w_lim, w_lim),
                          kwargs...)

    P_contours = contourf(xp, zp, clamp.(P, P_lims[1], P_lims[2])';
                          color = :matter,
                          levels = P_levels,
                          clims = P_lims,
                          kwargs...)

    P_profile = plot(averaged_P, zp,
                     linewidth = 2,
                     label = nothing,
                     xlims = (0.9, 1.3),
                     ylabel = "z (m)",
                     xlabel = "Plankton concentration (μM)")

    flux_plot = plot(times ./ hour, buoyancy_flux_time_series,
                     linewidth = 1,
                     label = "Buoyancy flux time series",
                     color = :black,
                     alpha = 0.4,
                     legend = :topright,
                     xlabel = "Time (hours)",
                     ylabel = "Buoyancy flux (m² s⁻³)",
                     ylims = (0.0, 1.1 * buoyancy_flux_parameters.initial_buoyancy_flux))

    plot!(flux_plot, times[1:i] ./ hour, buoyancy_flux_time_series[1:i],
          color = :steelblue,
          linewidth = 6,
          label = nothing)

    scatter!(flux_plot, times[i:i] / hour, buoyancy_flux_time_series[i:i],
             markershape = :circle,
             color = :steelblue,
             markerstrokewidth = 0,
             markersize = 15,
             label = "Current buoyancy flux")

    layout = Plots.grid(2, 2, widths=(0.7, 0.3))

    w_title = @sprintf("Vertical velocity (m s⁻¹) at %s", prettytime(t))
    P_title = @sprintf("Plankton concentration (μM) at %s", prettytime(t))

    plot(w_contours, flux_plot, P_contours, P_profile,
         title=[w_title "" P_title ""],
         layout=layout, size=(1000.5, 1000.5))
end

This page was generated using Literate.jl.