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, CairoMakie
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, Periodic, Periodic, 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"##277".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,
                            carbonate_system = CarbonateSystem(),
                            detritus = TwoParticleAndDissolved(grid; 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 GasExchange
└── 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{CUDAGPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CUDAGPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(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{CUDAGPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 15.625 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / 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 => 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

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, Periodic, Periodic, 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, Periodic, Periodic, 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=5.56792e-5, min=-5.24045e-5, mean=-2.31611e-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.27e-01
[ Info:     ... simulation initialization complete (12.686 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (5.183 seconds).
i:     20, sim time:    6 hours, wall time: 19.251 seconds, Δt: 25.164 minutes, CFL: 4.71e-01
i:     40, sim time:   15 hours, wall time: 20.666 seconds, Δt: 30 minutes, CFL: 5.84e-01
i:     60, sim time: 1.042 days, wall time: 21.967 seconds, Δt: 30 minutes, CFL: 6.05e-01
i:     80, sim time: 1.458 days, wall time: 23.265 seconds, Δt: 30 minutes, CFL: 6.46e-01
i:    100, sim time: 1.875 days, wall time: 24.564 seconds, Δt: 30 minutes, CFL: 6.63e-01
i:    120, sim time: 2.292 days, wall time: 25.867 seconds, Δt: 30 minutes, CFL: 7.45e-01
i:    140, sim time: 2.667 days, wall time: 27.050 seconds, Δt: 29.225 minutes, CFL: 7.50e-01
i:    160, sim time:     3 days, wall time: 28.346 seconds, Δt: 25.239 minutes, CFL: 7.50e-01
i:    180, sim time: 3.312 days, wall time: 29.660 seconds, Δt: 21.746 minutes, CFL: 7.50e-01
i:    200, sim time: 3.580 days, wall time: 30.957 seconds, Δt: 17.778 minutes, CFL: 7.50e-01
i:    220, sim time: 3.791 days, wall time: 32.264 seconds, Δt: 14.330 minutes, CFL: 7.50e-01
i:    240, sim time: 3.970 days, wall time: 33.558 seconds, Δt: 12.587 minutes, CFL: 7.50e-01
i:    260, sim time: 4.134 days, wall time: 34.857 seconds, Δt: 12.084 minutes, CFL: 7.50e-01
i:    280, sim time: 4.301 days, wall time: 36.150 seconds, Δt: 12.624 minutes, CFL: 7.50e-01
i:    300, sim time: 4.471 days, wall time: 37.450 seconds, Δt: 12.926 minutes, CFL: 7.50e-01
i:    320, sim time: 4.636 days, wall time: 38.744 seconds, Δt: 12.282 minutes, CFL: 7.50e-01
i:    340, sim time: 4.790 days, wall time: 40.032 seconds, Δt: 11.240 minutes, CFL: 7.50e-01
i:    360, sim time: 4.939 days, wall time: 41.364 seconds, Δt: 10.940 minutes, CFL: 7.50e-01
i:    380, sim time: 5.083 days, wall time: 42.585 seconds, Δt: 10.901 minutes, CFL: 7.50e-01
i:    400, sim time: 5.227 days, wall time: 43.861 seconds, Δt: 10.882 minutes, CFL: 7.50e-01
i:    420, sim time: 5.363 days, wall time: 45.136 seconds, Δt: 10.295 minutes, CFL: 7.50e-01
i:    440, sim time: 5.498 days, wall time: 46.409 seconds, Δt: 9.293 minutes, CFL: 7.50e-01
i:    460, sim time: 5.613 days, wall time: 47.686 seconds, Δt: 7.964 minutes, CFL: 7.50e-01
i:    480, sim time: 5.720 days, wall time: 48.962 seconds, Δt: 7.400 minutes, CFL: 7.50e-01
i:    500, sim time: 5.823 days, wall time: 50.239 seconds, Δt: 7.348 minutes, CFL: 7.50e-01
i:    520, sim time: 5.917 days, wall time: 51.518 seconds, Δt: 7.127 minutes, CFL: 7.50e-01
i:    540, sim time: 6.015 days, wall time: 52.853 seconds, Δt: 7.234 minutes, CFL: 7.50e-01
i:    560, sim time: 6.116 days, wall time: 54.067 seconds, Δt: 8.254 minutes, CFL: 7.50e-01
i:    580, sim time: 6.238 days, wall time: 55.339 seconds, Δt: 9.325 minutes, CFL: 7.50e-01
i:    600, sim time: 6.367 days, wall time: 56.623 seconds, Δt: 9.702 minutes, CFL: 7.50e-01
i:    620, sim time: 6.494 days, wall time: 57.889 seconds, Δt: 9.308 minutes, CFL: 7.50e-01
i:    640, sim time: 6.624 days, wall time: 59.168 seconds, Δt: 9.526 minutes, CFL: 7.50e-01
i:    660, sim time: 6.750 days, wall time: 1.007 minutes, Δt: 9.211 minutes, CFL: 7.50e-01
i:    680, sim time: 6.879 days, wall time: 1.029 minutes, Δt: 9.466 minutes, CFL: 7.50e-01
i:    700, sim time: 7.007 days, wall time: 1.053 minutes, Δt: 9.524 minutes, CFL: 7.50e-01
i:    720, sim time: 7.140 days, wall time: 1.071 minutes, Δt: 10.609 minutes, CFL: 7.50e-01
i:    740, sim time: 7.277 days, wall time: 1.092 minutes, Δt: 9.414 minutes, CFL: 7.50e-01
i:    760, sim time: 7.408 days, wall time: 1.113 minutes, Δt: 9.364 minutes, CFL: 7.50e-01
i:    780, sim time: 7.525 days, wall time: 1.134 minutes, Δt: 9.126 minutes, CFL: 7.50e-01
i:    800, sim time: 7.648 days, wall time: 1.156 minutes, Δt: 9.257 minutes, CFL: 7.50e-01
i:    820, sim time: 7.769 days, wall time: 1.178 minutes, Δt: 8.825 minutes, CFL: 7.50e-01
i:    840, sim time: 7.887 days, wall time: 1.198 minutes, Δt: 8.660 minutes, CFL: 7.50e-01
i:    860, sim time: 8.006 days, wall time: 1.222 minutes, Δt: 9.337 minutes, CFL: 7.50e-01
i:    880, sim time: 8.140 days, wall time: 1.240 minutes, Δt: 10.339 minutes, CFL: 7.50e-01
i:    900, sim time: 8.279 days, wall time: 1.261 minutes, Δt: 10.321 minutes, CFL: 7.50e-01
i:    920, sim time: 8.417 days, wall time: 1.283 minutes, Δt: 10.204 minutes, CFL: 7.50e-01
i:    940, sim time: 8.557 days, wall time: 1.304 minutes, Δt: 10.549 minutes, CFL: 7.50e-01
i:    960, sim time: 8.697 days, wall time: 1.325 minutes, Δt: 10.967 minutes, CFL: 7.50e-01
i:    980, sim time: 8.849 days, wall time: 1.348 minutes, Δt: 11.350 minutes, CFL: 7.50e-01
i:   1000, sim time:     9 days, wall time: 1.368 minutes, Δt: 11.231 minutes, CFL: 7.50e-01
i:   1020, sim time: 9.156 days, wall time: 1.389 minutes, Δt: 11.662 minutes, CFL: 7.50e-01
i:   1040, sim time: 9.304 days, wall time: 1.410 minutes, Δt: 10.689 minutes, CFL: 7.50e-01
i:   1060, sim time: 9.445 days, wall time: 1.431 minutes, Δt: 10.222 minutes, CFL: 7.50e-01
i:   1080, sim time: 9.583 days, wall time: 1.453 minutes, Δt: 10.077 minutes, CFL: 7.50e-01
i:   1100, sim time: 9.713 days, wall time: 1.474 minutes, Δt: 9.551 minutes, CFL: 7.50e-01
i:   1120, sim time: 9.840 days, wall time: 1.498 minutes, Δt: 9.631 minutes, CFL: 7.50e-01
i:   1140, sim time: 9.972 days, wall time: 1.516 minutes, Δt: 10.063 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.522 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.

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

This page was generated using Literate.jl.