Biogeochemistry in submesoscale eddies in the Eady model

In this example we setup a 3D model with a constant background buoyancy gradient with associated thermal wind (the Eady model) with the LOBSTER biogeochemical model. This example demonstrates how to use biogeochemistry in a more complicated physical model. The parameters in this example correspond roughly to those used by Taylor (2016) and result to the generation of a single submesoscale eddy.

Install dependencies

First we ensure we have the required dependencies installed

using Pkg
pkg "add OceanBioME, Oceananigans, CairoMakie"

Model setup

We load the required packages. Although not required, we also set the random seed to ensure reproducibility of the results.

using OceanBioME, Oceananigans, Printf
using Oceananigans.Units

using Random
Random.seed!(11)
Random.TaskLocalRNG()

Construct a grid with uniform grid spacing on GPU (just remove GPU() to run on CPU)

grid = RectilinearGrid(GPU(); size = (32, 32, 8), extent = (1kilometer, 1kilometer, 100meters))
32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1000.0) regularly spaced with Δx=31.25
├── Periodic y ∈ [0.0, 1000.0) regularly spaced with Δy=31.25
└── Bounded  z ∈ [-100.0, 0.0] regularly spaced with Δz=12.5

Set the Coriolis and buoyancy models.

coriolis = FPlane(f = 1e-4) # [s⁻¹]
buoyancy = SeawaterBuoyancy()
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078)

Specify parameters that are used to construct the background state.

background_state_parameters = (M = 1e-4,       # s⁻¹, geostrophic shear
                               f = coriolis.f, # s⁻¹, Coriolis parameter
                               N = 1e-4,       # s⁻¹, buoyancy frequency
                               H = grid.Lz,
                               g = buoyancy.gravitational_acceleration,
                               α = buoyancy.equation_of_state.thermal_expansion)
(M = 0.0001, f = 0.0001, N = 0.0001, H = 100.0, g = 9.80665, α = 0.000167)

We assume a background buoyancy $B$ with a constant stratification and also a constant lateral gradient (in the zonal direction). The background velocity components $U$ and $V$ are prescribed so that the thermal wind relationship is satisfied, that is, $f \partial_z U = - \partial_y B$ and $f \partial_z V = \partial_x B$.

T(x, y, z, t, p) = (p.M^2 * x + p.N^2 * (z + p.H)) / (p.g * p.α)
V(x, y, z, t, p) = p.M^2 / p.f * (z + p.H)

V_field = BackgroundField(V, parameters = background_state_parameters)
T_field = BackgroundField(T, parameters = background_state_parameters)
BackgroundField{typeof(Main.var"##225".T), @NamedTuple{M::Float64, f::Float64, N::Float64, H::Float64, g::Float64, α::Float64}}
├── func: T (generic function with 1 method)
└── parameters: (M = 0.0001, f = 0.0001, N = 0.0001, H = 100.0, g = 9.80665, α = 0.000167)

Specify some horizontal and vertical viscosity and diffusivity.

νᵥ = κᵥ = 1e-4 # [m² s⁻¹]
vertical_diffusivity = VerticalScalarDiffusivity(ν = νᵥ, κ = κᵥ)
VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=0.0001)

Setup the biogeochemical model with optional carbonate chemistry turned on.

biogeochemistry = LOBSTER(; grid,
                            carbonates = true,
                            open_bottom = true)

DIC_bcs = FieldBoundaryConditions(top = CarbonDioxideGasExchangeBoundaryCondition())
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: DiscreteBoundaryFunction with (::OceanBioME.Models.GasExchangeModel.GasExchange{Float64, OceanBioME.Models.GasExchangeModel.ScaledGasTransferVelocity.SchmidtScaledTransferVelocity{OceanBioME.Models.GasExchangeModel.PolynomialParameterisation{2, Tuple{Float64, Float64, Float64}}, OceanBioME.Models.GasExchangeModel.PolynomialParameterisation{4, NTuple{5, Float64}}}, OceanBioME.Models.GasExchangeModel.CarbonDioxideConcentration{OceanBioME.Models.CarbonChemistryModel.CarbonChemistry{OceanBioME.Models.CarbonChemistryModel.K0{Float64}, @NamedTuple{K1::OceanBioME.Models.CarbonChemistryModel.K1{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, K2::OceanBioME.Models.CarbonChemistryModel.K2{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}}, OceanBioME.Models.CarbonChemistryModel.KB{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.KS{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.KF{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, OceanBioME.Models.CarbonChemistryModel.KS{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, @NamedTuple{KP1::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, KP2::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, KP3::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}}, OceanBioME.Models.CarbonChemistryModel.KSi{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64}, OceanBioME.Models.CarbonChemistryModel.KW{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, OceanBioME.Models.CarbonChemistryModel.KSP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, typeof(OceanBioME.Models.teos10_polynomial_approximation)}, OceanBioME.Models.GasExchangeModel.PolynomialVirialCoefficientForCarbonDioxide{Float64}, OceanBioME.Models.GasExchangeModel.CrossVirialCoefficientForCarbonDioxide{Float64}, Float64, Nothing}, Float64})
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

Model instantiation

model = NonhydrostaticModel(; grid,
                              biogeochemistry,
                              boundary_conditions = (DIC = DIC_bcs, ),
                              advection = WENO(),
                              timestepper = :RungeKutta3,
                              coriolis,
                              tracers = (:T, :S),
                              buoyancy,
                              background_fields = (T = T_field, v = V_field),
                              closure = vertical_diffusivity)
NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: WENO{3, Float64, Float32}(order=5)
├── tracers: (T, S, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk)
├── closure: VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=(T=0.0001, S=0.0001, NO₃=0.0001, NH₄=0.0001, P=0.0001, Z=0.0001, sPOM=0.0001, bPOM=0.0001, DOM=0.0001, DIC=0.0001, Alk=0.0001))
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)

Initial conditions

Start with a bit of random noise added to the background thermal wind and an arbitary biogeochemical state.

Ξ(z) = randn() * z / grid.Lz * (z / grid.Lz + 1)

Ũ = 1e-3
uᵢ(x, y, z) = Ũ * Ξ(z)
vᵢ(x, y, z) = Ũ * Ξ(z)

set!(model, u=uᵢ, v=vᵢ, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2200.0, Alk = 2409.0, S = 35, T = 20)

Setup the simulation

Choose an appropriate initial timestep for this resolution and set up the simulation

Δx = minimum_xspacing(grid, Center(), Center(), Center())
Δy = minimum_yspacing(grid, Center(), Center(), Center())
Δz = minimum_zspacing(grid, Center(), Center(), Center())

Δt₀ = 0.75 * min(Δx, Δy, Δz) / V(0, 0, 0, 0, background_state_parameters)

simulation = Simulation(model, Δt = Δt₀, stop_time = 10days)
Simulation of NonhydrostaticModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 15.625 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 10 days
├── Stop iteration: Inf
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => 4
│   ├── stop_iteration_exceeded => -
│   ├── wall_time_limit_exceeded => e
│   └── nan_checker => }
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entries

Adapt the time step while keeping the CFL number fixed.

wizard = TimeStepWizard(cfl = 0.75, diffusive_cfl = 0.75, max_Δt = 30minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(5))

Create a progress message.

progress(sim) = @printf("i: % 6d, sim time: % 10s, wall time: % 10s, Δt: % 10s, CFL: %.2e\n",
                        sim.model.clock.iteration,
                        prettytime(sim.model.clock.time),
                        prettytime(sim.run_wall_time),
                        prettytime(sim.Δt),
                        AdvectiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(20))
Callback of progress on IterationInterval(20)

Here, we add some diagnostics to calculate and output.

u, v, w = model.velocities # unpack velocity `Field`s
NamedTuple with 3 Fields on 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo:
├── u: 32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── v: 32×32×8 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CUDAGPU
└── w: 32×32×9 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on Oceananigans.Grids.RectilinearGrid on CUDAGPU

and also calculate the vertical vorticity.

ζ = Field(∂x(v) - ∂y(u))
32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 38×38×14 OffsetArray(::CUDA.CuArray{Float64, 3, CUDA.DeviceMemory}, -2:35, -2:35, -2:11) with eltype Float64 with indices -2:35×-2:35×-2:11
    └── max=6.41337e-5, min=-5.24045e-5, mean=-1.28213e-23

Periodically save the velocities and vorticity to a file.

simulation.output_writers[:fields] = JLD2Writer(model, merge(model.tracers, (; u, v, w, ζ));
                                                schedule = TimeInterval(2hours),
                                                filename = "eady_turbulence_bgc",
                                                overwrite_existing = true)

Run the simulation

run!(simulation)
[ Info: Initializing simulation...
i:      0, sim time:  0 seconds, wall time:  0 seconds, Δt: 17.188 minutes, CFL: 3.25e-01
[ Info:     ... simulation initialization complete (3.052 minutes)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (17.017 seconds).
i:     20, sim time:    6 hours, wall time: 3.347 minutes, Δt: 25.164 minutes, CFL: 4.70e-01
i:     40, sim time:   15 hours, wall time: 3.351 minutes, Δt: 30 minutes, CFL: 5.87e-01
i:     60, sim time: 1.042 days, wall time: 3.355 minutes, Δt: 30 minutes, CFL: 6.06e-01
i:     80, sim time: 1.458 days, wall time: 3.358 minutes, Δt: 30 minutes, CFL: 6.34e-01
i:    100, sim time: 1.875 days, wall time: 3.362 minutes, Δt: 30 minutes, CFL: 6.52e-01
i:    120, sim time: 2.292 days, wall time: 3.366 minutes, Δt: 30 minutes, CFL: 7.16e-01
i:    140, sim time: 2.687 days, wall time: 3.369 minutes, Δt: 28.625 minutes, CFL: 7.50e-01
i:    160, sim time:     3 days, wall time: 3.373 minutes, Δt: 25.231 minutes, CFL: 7.50e-01
i:    180, sim time: 3.308 days, wall time: 3.377 minutes, Δt: 19.652 minutes, CFL: 7.50e-01
i:    200, sim time: 3.548 days, wall time: 3.380 minutes, Δt: 16.176 minutes, CFL: 7.50e-01
i:    220, sim time: 3.749 days, wall time: 3.384 minutes, Δt: 14.433 minutes, CFL: 7.50e-01
i:    240, sim time: 3.926 days, wall time: 3.387 minutes, Δt: 13.164 minutes, CFL: 7.50e-01
i:    260, sim time: 4.092 days, wall time: 3.391 minutes, Δt: 13.173 minutes, CFL: 7.50e-01
i:    280, sim time: 4.250 days, wall time: 3.394 minutes, Δt: 12.970 minutes, CFL: 7.50e-01
i:    300, sim time: 4.413 days, wall time: 3.398 minutes, Δt: 12.543 minutes, CFL: 7.50e-01
i:    320, sim time: 4.577 days, wall time: 3.401 minutes, Δt: 12.198 minutes, CFL: 7.50e-01
i:    340, sim time: 4.743 days, wall time: 3.405 minutes, Δt: 11.900 minutes, CFL: 7.50e-01
i:    360, sim time: 4.909 days, wall time: 3.408 minutes, Δt: 11.816 minutes, CFL: 7.50e-01
i:    380, sim time: 5.066 days, wall time: 3.412 minutes, Δt: 11.594 minutes, CFL: 7.50e-01
i:    400, sim time: 5.216 days, wall time: 3.415 minutes, Δt: 11.887 minutes, CFL: 7.50e-01
i:    420, sim time: 5.378 days, wall time: 3.418 minutes, Δt: 12.793 minutes, CFL: 7.50e-01
i:    440, sim time: 5.544 days, wall time: 3.422 minutes, Δt: 12.223 minutes, CFL: 7.50e-01
i:    460, sim time: 5.698 days, wall time: 3.425 minutes, Δt: 11.115 minutes, CFL: 7.50e-01
i:    480, sim time: 5.841 days, wall time: 3.429 minutes, Δt: 10.360 minutes, CFL: 7.50e-01
i:    500, sim time: 5.980 days, wall time: 3.431 minutes, Δt: 9.878 minutes, CFL: 7.50e-01
i:    520, sim time: 6.109 days, wall time: 3.434 minutes, Δt: 8.825 minutes, CFL: 7.50e-01
i:    540, sim time: 6.226 days, wall time: 3.436 minutes, Δt: 8.569 minutes, CFL: 7.50e-01
i:    560, sim time: 6.333 days, wall time: 3.439 minutes, Δt: 8.053 minutes, CFL: 7.50e-01
i:    580, sim time: 6.438 days, wall time: 3.441 minutes, Δt: 7.533 minutes, CFL: 7.50e-01
i:    600, sim time: 6.535 days, wall time: 3.444 minutes, Δt: 7.065 minutes, CFL: 7.50e-01
i:    620, sim time: 6.630 days, wall time: 3.446 minutes, Δt: 6.587 minutes, CFL: 7.50e-01
i:    640, sim time: 6.722 days, wall time: 3.449 minutes, Δt: 6.692 minutes, CFL: 7.50e-01
i:    660, sim time: 6.815 days, wall time: 3.452 minutes, Δt: 6.533 minutes, CFL: 7.50e-01
i:    680, sim time: 6.902 days, wall time: 3.454 minutes, Δt: 6.739 minutes, CFL: 7.50e-01
i:    700, sim time: 6.995 days, wall time: 3.457 minutes, Δt: 7.441 minutes, CFL: 7.50e-01
i:    720, sim time: 7.100 days, wall time: 3.459 minutes, Δt: 8.095 minutes, CFL: 7.50e-01
i:    740, sim time: 7.204 days, wall time: 3.462 minutes, Δt: 7.773 minutes, CFL: 7.50e-01
i:    760, sim time: 7.311 days, wall time: 3.465 minutes, Δt: 7.987 minutes, CFL: 7.50e-01
i:    780, sim time: 7.422 days, wall time: 3.468 minutes, Δt: 8.171 minutes, CFL: 7.50e-01
i:    800, sim time: 7.537 days, wall time: 3.470 minutes, Δt: 9.040 minutes, CFL: 7.50e-01
i:    820, sim time: 7.659 days, wall time: 3.473 minutes, Δt: 8.828 minutes, CFL: 7.50e-01
i:    840, sim time: 7.775 days, wall time: 3.477 minutes, Δt: 8.821 minutes, CFL: 7.50e-01
i:    860, sim time: 7.896 days, wall time: 3.480 minutes, Δt: 9.231 minutes, CFL: 7.50e-01
i:    880, sim time: 8.019 days, wall time: 3.483 minutes, Δt: 8.945 minutes, CFL: 7.50e-01
i:    900, sim time: 8.136 days, wall time: 3.486 minutes, Δt: 8.205 minutes, CFL: 7.50e-01
i:    920, sim time: 8.244 days, wall time: 3.489 minutes, Δt: 7.660 minutes, CFL: 7.50e-01
i:    940, sim time: 8.344 days, wall time: 3.493 minutes, Δt: 7.684 minutes, CFL: 7.50e-01
i:    960, sim time: 8.449 days, wall time: 3.496 minutes, Δt: 7.679 minutes, CFL: 7.50e-01
i:    980, sim time: 8.554 days, wall time: 3.499 minutes, Δt: 7.947 minutes, CFL: 7.50e-01
i:   1000, sim time: 8.660 days, wall time: 3.502 minutes, Δt: 8.957 minutes, CFL: 7.50e-01
i:   1020, sim time: 8.776 days, wall time: 3.505 minutes, Δt: 9.564 minutes, CFL: 7.50e-01
i:   1040, sim time: 8.909 days, wall time: 3.508 minutes, Δt: 9.797 minutes, CFL: 7.50e-01
i:   1060, sim time: 9.035 days, wall time: 3.511 minutes, Δt: 9.899 minutes, CFL: 7.50e-01
i:   1080, sim time: 9.162 days, wall time: 3.514 minutes, Δt: 9.016 minutes, CFL: 7.50e-01
i:   1100, sim time: 9.280 days, wall time: 3.518 minutes, Δt: 8.442 minutes, CFL: 7.50e-01
i:   1120, sim time: 9.391 days, wall time: 3.521 minutes, Δt: 8.132 minutes, CFL: 7.50e-01
i:   1140, sim time: 9.500 days, wall time: 3.524 minutes, Δt: 8.251 minutes, CFL: 7.50e-01
i:   1160, sim time: 9.613 days, wall time: 3.527 minutes, Δt: 8.571 minutes, CFL: 7.50e-01
i:   1180, sim time: 9.733 days, wall time: 3.530 minutes, Δt: 8.967 minutes, CFL: 7.50e-01
i:   1200, sim time: 9.853 days, wall time: 3.534 minutes, Δt: 9.340 minutes, CFL: 7.50e-01
i:   1220, sim time: 9.982 days, wall time: 3.536 minutes, Δt: 9.456 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 3.537 minutes.
[ Info: Simulation time 10 days equals or exceeds stop time 10 days.

Now load the saved output,

  ζ = FieldTimeSeries("eady_turbulence_bgc.jld2", "ζ")
  P = FieldTimeSeries("eady_turbulence_bgc.jld2", "P")
NO₃ = FieldTimeSeries("eady_turbulence_bgc.jld2", "NO₃")
NH₄ = FieldTimeSeries("eady_turbulence_bgc.jld2", "NH₄")
DIC = FieldTimeSeries("eady_turbulence_bgc.jld2", "DIC")

times = ζ.times

xζ, yζ, zζ = nodes(ζ)
xc, yc, zc = nodes(P)

and plot.

using CairoMakie

n = Observable(1)

  ζₙ = @lift interior(  ζ[$n], :, :, grid.Nz)
  Nₙ = @lift interior(NO₃[$n], :, :, grid.Nz) .+ interior(NH₄[$n], :, :, grid.Nz)
  Pₙ = @lift interior(  P[$n], :, :, grid.Nz)
DICₙ = @lift interior(DIC[$n], :, :, grid.Nz)

fig = Figure(size = (1600, 1600), fontsize = 20)

lims = [(minimum(T), maximum(T)) for T in (  ζ[:, :, grid.Nz, :],
                                           NO₃[:, :, grid.Nz, :] .+ NH₄[:, :, grid.Nz, :],
                                             P[:, :, grid.Nz, :],
                                           DIC[:, :, grid.Nz, :])]

axis_kwargs = (xlabel = "x (m)", ylabel = "y (m)", aspect = DataAspect())

ax1 = Axis(fig[1, 1]; title = "Vertical vorticity (1 / s)", axis_kwargs...)
hm1 = heatmap!(ax1, xζ, yζ, ζₙ, colormap = :balance, colorrange = lims[1])
Colorbar(fig[1, 2], hm1)

ax2 = Axis(fig[1, 3]; title = "Nutrient (NO₃ + NH₄) concentration (mmol N / m³)", axis_kwargs...)
hm2 = heatmap!(ax2, xc, yc, Nₙ, colormap = Reverse(:bamako), colorrange = lims[2])
Colorbar(fig[1, 4], hm2)

ax3 = Axis(fig[2, 1]; title = "Phytoplankton concentration (mmol N / m³)", axis_kwargs...)
hm3 = heatmap!(ax3, xc, yc, Pₙ, colormap = Reverse(:batlow), colorrange = lims[3])
Colorbar(fig[2, 2], hm3)

ax4 = Axis(fig[2, 3]; title = "Dissolved inorganic carbon (mmol C / m³)", axis_kwargs...)
hm4 = heatmap!(ax4, xc, yc, DICₙ, colormap = Reverse(:devon), colorrange = lims[4])
Colorbar(fig[2, 4], hm4)

title = @lift "t = $(prettytime(times[$n]))"
Label(fig[0, :], title, fontsize = 30)

record(fig, "eady.mp4", 1:length(times), framerate = 12) do i
    n[] = i
end
WARNING: using IntervalArithmetic.numtype in module TaylorSeries conflicts with an existing identifier.
┌ Warning: Error requiring `IntervalArithmetic` from `TaylorSeries`
│   exception =
│    LoadError: UndefVarError: `IntervalBox` not defined
│    Stacktrace:
│      [1] top-level scope
│        @ ~/.julia-oceananigans/packages/TaylorSeries/ZBRjU/src/intervals.jl:171
│      [2] include(mod::Module, _path::String)
│        @ Base ./Base.jl:495
│      [3] include(x::String)
│        @ TaylorSeries ~/.julia-oceananigans/packages/TaylorSeries/ZBRjU/src/TaylorSeries.jl:17
│      [4] top-level scope
│        @ ~/.julia-oceananigans/packages/Requires/1eCOK/src/Requires.jl:40
│      [5] eval
│        @ ./boot.jl:385 [inlined]
│      [6] eval
│        @ ~/.julia-oceananigans/packages/TaylorSeries/ZBRjU/src/TaylorSeries.jl:17 [inlined]
│      [7] (::TaylorSeries.var"#54#57")()
│        @ TaylorSeries ~/.julia-oceananigans/packages/Requires/1eCOK/src/require.jl:101
│      [8] macro expansion
│        @ timing.jl:395 [inlined]
│      [9] err(f::Any, listener::Module, modname::String, file::String, line::Any)
│        @ Requires ~/.julia-oceananigans/packages/Requires/1eCOK/src/require.jl:47
│     [10] (::TaylorSeries.var"#53#56")()
│        @ TaylorSeries ~/.julia-oceananigans/packages/Requires/1eCOK/src/require.jl:100
│     [11] withpath(f::Any, path::String)
│        @ Requires ~/.julia-oceananigans/packages/Requires/1eCOK/src/require.jl:37
│     [12] (::TaylorSeries.var"#52#55")()
│        @ TaylorSeries ~/.julia-oceananigans/packages/Requires/1eCOK/src/require.jl:99
│     [13] #invokelatest#2
│        @ ./essentials.jl:892 [inlined]
│     [14] invokelatest
│        @ ./essentials.jl:889 [inlined]
│     [15] foreach(f::typeof(invokelatest), itr::Vector{Function})
│        @ Base ./abstractarray.jl:3098
│     [16] loadpkg(pkg::Base.PkgId)
│        @ Requires ~/.julia-oceananigans/packages/Requires/1eCOK/src/require.jl:27
│     [17] #invokelatest#2
│        @ ./essentials.jl:892 [inlined]
│     [18] invokelatest
│        @ ./essentials.jl:889 [inlined]
│     [19] run_package_callbacks(modkey::Base.PkgId)
│        @ Base ./loading.jl:1229
│     [20] _tryrequire_from_serialized(modkey::Base.PkgId, path::String, ocachepath::String, sourcepath::String, depmods::Vector{Any})
│        @ Base ./loading.jl:1557
│     [21] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt128)
│        @ Base ./loading.jl:1644
│     [22] _require(pkg::Base.PkgId, env::String)
│        @ Base ./loading.jl:2008
│     [23] __require_prelocked(uuidkey::Base.PkgId, env::String)
│        @ Base ./loading.jl:1882
│     [24] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│     [25] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│     [26] _require_prelocked(uuidkey::Base.PkgId, env::String)
│        @ Base ./loading.jl:1873
│     [27] macro expansion
│        @ ./loading.jl:1860 [inlined]
│     [28] macro expansion
│        @ ./lock.jl:267 [inlined]
│     [29] __require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1823
│     [30] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│     [31] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│     [32] require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1816
│     [33] eval
│        @ ./boot.jl:385 [inlined]
│     [34] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
│        @ Base ./loading.jl:2146
│     [35] include_string
│        @ ./loading.jl:2156 [inlined]
│     [36] #47
│        @ ~/.julia-oceananigans/packages/Literate/TujCy/src/Literate.jl:934 [inlined]
│     [37] task_local_storage(body::Literate.var"#47#49"{String, Bool, Module, String}, key::Symbol, val::String)
│        @ Base ./task.jl:304
│     [38] #46
│        @ ~/.julia-oceananigans/packages/Literate/TujCy/src/Literate.jl:930 [inlined]
│     [39] (::IOCapture.var"#5#9"{Core.TypeofBottom, Literate.var"#46#48"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, Base.TTY, Base.TTY})()
│        @ IOCapture ~/.julia-oceananigans/packages/IOCapture/Y5rEA/src/IOCapture.jl:170
│     [40] with_logstate(f::Function, logstate::Any)
│        @ Base.CoreLogging ./logging.jl:517
│     [41] with_logger
│        @ ./logging.jl:630 [inlined]
│     [42] capture(f::Literate.var"#46#48"{String, Bool, Module, String}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer, io_context::Vector{Any})
│        @ IOCapture ~/.julia-oceananigans/packages/IOCapture/Y5rEA/src/IOCapture.jl:167
│     [43] execute_block(sb::Module, block::String; inputfile::String, fake_source::String, softscope::Bool, continue_on_error::Bool)
│        @ Literate ~/.julia-oceananigans/packages/Literate/TujCy/src/Literate.jl:928
│     [44] execute_block
│        @ ~/.julia-oceananigans/packages/Literate/TujCy/src/Literate.jl:912 [inlined]
│     [45] execute_markdown!(io::IOBuffer, sb::Module, block::String, outputdir::String; inputfile::String, fake_source::String, flavor::Literate.DocumenterFlavor, image_formats::Vector{Tuple{MIME, String}}, file_prefix::String, softscope::Bool, continue_on_error::Bool)
│        @ Literate ~/.julia-oceananigans/packages/Literate/TujCy/src/Literate.jl:669
│     [46] (::Literate.var"#30#32"{Dict{String, Any}, IOBuffer, Module, Literate.CodeChunk, Int64})()
│        @ Literate ~/.julia-oceananigans/packages/Literate/TujCy/src/Literate.jl:640
│     [47] cd(f::Literate.var"#30#32"{Dict{String, Any}, IOBuffer, Module, Literate.CodeChunk, Int64}, dir::String)
│        @ Base.Filesystem ./file.jl:112
│     [48] markdown(inputfile::String, outputdir::String; config::Dict{Any, Any}, kwargs::@Kwargs{flavor::Literate.DocumenterFlavor, repo_root_url::String, execute::Bool, postprocess::typeof(replace_silly_warning)})
│        @ Literate ~/.julia-oceananigans/packages/Literate/TujCy/src/Literate.jl:639
│     [49] top-level scope
│        @ ~/Oceananigans.jl-997/docs/make_examples.jl:12
│     [50] include(mod::Module, _path::String)
│        @ Base ./Base.jl:495
│     [51] exec_options(opts::Base.JLOptions)
│        @ Base ./client.jl:323
│    in expression starting at /var/lib/buildkite-agent/.julia-oceananigans/packages/TaylorSeries/ZBRjU/src/intervals.jl:171
└ @ Requires ~/.julia-oceananigans/packages/Requires/1eCOK/src/require.jl:51


This page was generated using Literate.jl.