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.

grid = RectilinearGrid(size = (32, 32, 8), extent = (1kilometer, 1kilometer, 100meters))

32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU 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"##296".T), NamedTuple{(:M, :f, :N, :H, :g, :α), NTuple{6, 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)

LOBSTER{Float64} with carbonates ✅, oxygen ❌, variable Redfield ratio ❌and (:sPOM, :bPOM) sinking 
 Light attenuation: Two-band light attenuation model (Float64)
 Sediment: Nothing
 Particles: Nothing
 Modifiers: Nothing

To-do: change to a buoyancy parameterisation so we don't have to fake the temperature and salinity.

DIC_bcs = FieldBoundaryConditions(top = GasExchange(; gas = :CO₂))

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: ContinuousBoundaryFunction (::GasExchange{Val{:CO₂}, NamedTuple{(:A, :B, :C, :D), NTuple{4, Float64}}, NamedTuple{(:A₁, :A₂, :A₃, :B₁, :B₂, :B₃), NTuple{6, Float64}}, Float64, Float64, Float64, OceanBioME.Boundaries.pCO₂{NamedTuple{(:C, :invT, :logCT, :ClogT, :T², :ST², :ST, :S), NTuple{8, Float64}}, NamedTuple{(:C, :S, :S², :invT, :logT), NTuple{5, Float64}}, NamedTuple{(:C, :S, :S², :invT, :logT), NTuple{5, Float64}}, NamedTuple{(:C, :invT, :invTsqrtS, :invTS, :invTS¹⁵, :invTS², :sqrtS, :S, :logT, :logTsqrtS, :logTS, :TsqrtS), NTuple{12, Float64}}, NamedTuple{(:C, :invT, :logT, :sqrtSinvT, :sqrtS, :sqrtSlogT, :S), NTuple{7, Float64}}, Float64}}) at (Nothing, Nothing, Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

Model instantiation

model = NonhydrostaticModel(; grid,
                              biogeochemistry,
                              boundary_conditions = (DIC = DIC_bcs, ),
                              advection = WENO(grid),
                              timestepper = :RungeKutta3,
                              coriolis,
                              tracers = (:T, :S),
                              buoyancy,
                              background_fields = (T = T_field, v = V_field),
                              closure = vertical_diffusivity)

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── 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{CPU, 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
├── Callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
├── 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 Oceananigans.Architectures.CPU with 3×3×3 halo:
├── u: 32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on Oceananigans.Architectures.CPU
├── v: 32×32×8 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on Oceananigans.Architectures.CPU
└── w: 32×32×9 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on Oceananigans.Grids.RectilinearGrid on Oceananigans.Architectures.CPU

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 Oceananigans.Architectures.CPU
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 38×38×14 OffsetArray(::Array{Float64, 3}, -2:35, -2:35, -2:11) with eltype Float64 with indices -2:35×-2:35×-2:11
    └── max=0.0, min=0.0, mean=0.0

Periodically save the velocities and vorticity to a file.

simulation.output_writers[:fields] = JLD2OutputWriter(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.27e-01
[ Info:     ... simulation initialization complete (6.143 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (10.318 seconds).
i:     20, sim time:    6 hours, wall time: 18.779 seconds, Δt: 25.164 minutes, CFL: 4.74e-01
i:     40, sim time:   15 hours, wall time: 21.312 seconds, Δt: 30 minutes, CFL: 5.78e-01
i:     60, sim time: 1.042 days, wall time: 23.833 seconds, Δt: 30 minutes, CFL: 6.15e-01
i:     80, sim time: 1.458 days, wall time: 26.365 seconds, Δt: 30 minutes, CFL: 6.35e-01
i:    100, sim time: 1.875 days, wall time: 28.895 seconds, Δt: 30 minutes, CFL: 6.68e-01
i:    120, sim time: 2.292 days, wall time: 31.442 seconds, Δt: 29.425 minutes, CFL: 7.50e-01
i:    140, sim time: 2.621 days, wall time: 33.995 seconds, Δt: 26.070 minutes, CFL: 7.50e-01
i:    160, sim time: 2.933 days, wall time: 36.535 seconds, Δt: 22.899 minutes, CFL: 7.50e-01
i:    180, sim time: 3.210 days, wall time: 39.059 seconds, Δt: 20.464 minutes, CFL: 7.50e-01
i:    200, sim time: 3.486 days, wall time: 41.576 seconds, Δt: 18.964 minutes, CFL: 7.50e-01
i:    220, sim time: 3.715 days, wall time: 44.105 seconds, Δt: 16.861 minutes, CFL: 7.50e-01
i:    240, sim time: 3.928 days, wall time: 46.632 seconds, Δt: 15.505 minutes, CFL: 7.50e-01
i:    260, sim time: 4.121 days, wall time: 49.170 seconds, Δt: 13.340 minutes, CFL: 7.50e-01
i:    280, sim time: 4.295 days, wall time: 51.713 seconds, Δt: 12.893 minutes, CFL: 7.50e-01
i:    300, sim time: 4.461 days, wall time: 54.245 seconds, Δt: 12.870 minutes, CFL: 7.50e-01
i:    320, sim time: 4.631 days, wall time: 56.763 seconds, Δt: 13.847 minutes, CFL: 7.50e-01
i:    340, sim time: 4.814 days, wall time: 59.292 seconds, Δt: 12.996 minutes, CFL: 7.50e-01
i:    360, sim time: 4.975 days, wall time: 1.030 minutes, Δt: 11.849 minutes, CFL: 7.50e-01
i:    380, sim time: 5.127 days, wall time: 1.074 minutes, Δt: 10.534 minutes, CFL: 7.50e-01
i:    400, sim time: 5.265 days, wall time: 1.116 minutes, Δt: 10.396 minutes, CFL: 7.50e-01
i:    420, sim time: 5.407 days, wall time: 1.159 minutes, Δt: 10.576 minutes, CFL: 7.50e-01
i:    440, sim time: 5.542 days, wall time: 1.202 minutes, Δt: 10.064 minutes, CFL: 7.50e-01
i:    460, sim time: 5.674 days, wall time: 1.244 minutes, Δt: 9.861 minutes, CFL: 7.50e-01
i:    480, sim time: 5.805 days, wall time: 1.286 minutes, Δt: 9.628 minutes, CFL: 7.50e-01
i:    500, sim time: 5.930 days, wall time: 1.329 minutes, Δt: 9.718 minutes, CFL: 7.50e-01
i:    520, sim time: 6.062 days, wall time: 1.372 minutes, Δt: 9.949 minutes, CFL: 7.50e-01
i:    540, sim time: 6.188 days, wall time: 1.415 minutes, Δt: 10.059 minutes, CFL: 7.50e-01
i:    560, sim time: 6.330 days, wall time: 1.457 minutes, Δt: 10.468 minutes, CFL: 7.50e-01
i:    580, sim time: 6.469 days, wall time: 1.499 minutes, Δt: 10.910 minutes, CFL: 7.50e-01
i:    600, sim time: 6.605 days, wall time: 1.541 minutes, Δt: 10.324 minutes, CFL: 7.50e-01
i:    620, sim time: 6.746 days, wall time: 1.583 minutes, Δt: 10.431 minutes, CFL: 7.50e-01
i:    640, sim time: 6.885 days, wall time: 1.625 minutes, Δt: 10.639 minutes, CFL: 7.50e-01
i:    660, sim time: 7.015 days, wall time: 1.668 minutes, Δt: 10.736 minutes, CFL: 7.50e-01
i:    680, sim time: 7.158 days, wall time: 1.710 minutes, Δt: 10.893 minutes, CFL: 7.50e-01
i:    700, sim time: 7.295 days, wall time: 1.752 minutes, Δt: 10.534 minutes, CFL: 7.50e-01
i:    720, sim time: 7.424 days, wall time: 1.795 minutes, Δt: 10.389 minutes, CFL: 7.50e-01
i:    740, sim time: 7.563 days, wall time: 1.837 minutes, Δt: 9.627 minutes, CFL: 7.50e-01
i:    760, sim time: 7.685 days, wall time: 1.879 minutes, Δt: 8.875 minutes, CFL: 7.50e-01
i:    780, sim time: 7.805 days, wall time: 1.921 minutes, Δt: 8.765 minutes, CFL: 7.50e-01
i:    800, sim time: 7.923 days, wall time: 1.963 minutes, Δt: 8.729 minutes, CFL: 7.50e-01
i:    820, sim time: 8.044 days, wall time: 2.005 minutes, Δt: 9.355 minutes, CFL: 7.50e-01
i:    840, sim time: 8.167 days, wall time: 2.047 minutes, Δt: 10.167 minutes, CFL: 7.50e-01
i:    860, sim time: 8.301 days, wall time: 2.089 minutes, Δt: 10.343 minutes, CFL: 7.50e-01
i:    880, sim time: 8.440 days, wall time: 2.131 minutes, Δt: 10.806 minutes, CFL: 7.50e-01
i:    900, sim time: 8.581 days, wall time: 2.173 minutes, Δt: 10.526 minutes, CFL: 7.50e-01
i:    920, sim time: 8.710 days, wall time: 2.215 minutes, Δt: 10.145 minutes, CFL: 7.50e-01
i:    940, sim time: 8.847 days, wall time: 2.257 minutes, Δt: 10.212 minutes, CFL: 7.50e-01
i:    960, sim time: 8.982 days, wall time: 2.299 minutes, Δt: 10.450 minutes, CFL: 7.50e-01
i:    980, sim time: 9.122 days, wall time: 2.341 minutes, Δt: 11.325 minutes, CFL: 7.50e-01
i:   1000, sim time: 9.266 days, wall time: 2.383 minutes, Δt: 11.137 minutes, CFL: 7.50e-01
i:   1020, sim time: 9.416 days, wall time: 2.425 minutes, Δt: 10.677 minutes, CFL: 7.50e-01
i:   1040, sim time: 9.544 days, wall time: 2.467 minutes, Δt: 10.561 minutes, CFL: 7.50e-01
i:   1060, sim time: 9.674 days, wall time: 2.509 minutes, Δt: 10.361 minutes, CFL: 7.50e-01
i:   1080, sim time: 9.813 days, wall time: 2.551 minutes, Δt: 10.004 minutes, CFL: 7.50e-01
i:   1100, sim time: 9.952 days, wall time: 2.593 minutes, Δt: 10.257 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 2.608 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(resolution = (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

This page was generated using Literate.jl.