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.505 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (5.355 seconds).
i:     20, sim time:    6 hours, wall time: 19.235 seconds, Δt: 25.164 minutes, CFL: 4.71e-01
i:     40, sim time:   15 hours, wall time: 20.639 seconds, Δt: 30 minutes, CFL: 5.84e-01
i:     60, sim time: 1.042 days, wall time: 21.925 seconds, Δt: 30 minutes, CFL: 6.05e-01
i:     80, sim time: 1.458 days, wall time: 23.212 seconds, Δt: 30 minutes, CFL: 6.46e-01
i:    100, sim time: 1.875 days, wall time: 24.502 seconds, Δt: 30 minutes, CFL: 6.63e-01
i:    120, sim time: 2.292 days, wall time: 25.792 seconds, Δt: 30 minutes, CFL: 7.45e-01
i:    140, sim time: 2.667 days, wall time: 26.962 seconds, Δt: 29.225 minutes, CFL: 7.50e-01
i:    160, sim time:     3 days, wall time: 28.249 seconds, Δt: 25.239 minutes, CFL: 7.50e-01
i:    180, sim time: 3.312 days, wall time: 29.531 seconds, Δt: 21.746 minutes, CFL: 7.50e-01
i:    200, sim time: 3.580 days, wall time: 30.820 seconds, Δt: 17.778 minutes, CFL: 7.50e-01
i:    220, sim time: 3.791 days, wall time: 32.104 seconds, Δt: 14.330 minutes, CFL: 7.50e-01
i:    240, sim time: 3.970 days, wall time: 33.391 seconds, Δt: 12.587 minutes, CFL: 7.50e-01
i:    260, sim time: 4.134 days, wall time: 34.678 seconds, Δt: 12.084 minutes, CFL: 7.50e-01
i:    280, sim time: 4.301 days, wall time: 35.962 seconds, Δt: 12.624 minutes, CFL: 7.50e-01
i:    300, sim time: 4.471 days, wall time: 37.242 seconds, Δt: 12.926 minutes, CFL: 7.50e-01
i:    320, sim time: 4.636 days, wall time: 38.526 seconds, Δt: 12.282 minutes, CFL: 7.50e-01
i:    340, sim time: 4.790 days, wall time: 39.803 seconds, Δt: 11.240 minutes, CFL: 7.50e-01
i:    360, sim time: 4.939 days, wall time: 41.140 seconds, Δt: 10.940 minutes, CFL: 7.50e-01
i:    380, sim time: 5.083 days, wall time: 42.354 seconds, Δt: 10.901 minutes, CFL: 7.50e-01
i:    400, sim time: 5.227 days, wall time: 43.630 seconds, Δt: 10.882 minutes, CFL: 7.50e-01
i:    420, sim time: 5.363 days, wall time: 44.912 seconds, Δt: 10.295 minutes, CFL: 7.50e-01
i:    440, sim time: 5.498 days, wall time: 46.181 seconds, Δt: 9.293 minutes, CFL: 7.50e-01
i:    460, sim time: 5.613 days, wall time: 47.463 seconds, Δt: 7.964 minutes, CFL: 7.50e-01
i:    480, sim time: 5.720 days, wall time: 48.731 seconds, Δt: 7.400 minutes, CFL: 7.50e-01
i:    500, sim time: 5.823 days, wall time: 50.002 seconds, Δt: 7.348 minutes, CFL: 7.50e-01
i:    520, sim time: 5.917 days, wall time: 51.274 seconds, Δt: 7.127 minutes, CFL: 7.50e-01
i:    540, sim time: 6.015 days, wall time: 52.603 seconds, Δt: 7.234 minutes, CFL: 7.50e-01
i:    560, sim time: 6.116 days, wall time: 53.836 seconds, Δt: 8.254 minutes, CFL: 7.50e-01
i:    580, sim time: 6.238 days, wall time: 55.114 seconds, Δt: 9.325 minutes, CFL: 7.50e-01
i:    600, sim time: 6.367 days, wall time: 56.403 seconds, Δt: 9.702 minutes, CFL: 7.50e-01
i:    620, sim time: 6.494 days, wall time: 57.681 seconds, Δt: 9.308 minutes, CFL: 7.50e-01
i:    640, sim time: 6.624 days, wall time: 58.966 seconds, Δt: 9.526 minutes, CFL: 7.50e-01
i:    660, sim time: 6.750 days, wall time: 1.004 minutes, Δt: 9.211 minutes, CFL: 7.50e-01
i:    680, sim time: 6.879 days, wall time: 1.025 minutes, Δt: 9.466 minutes, CFL: 7.50e-01
i:    700, sim time: 7.007 days, wall time: 1.049 minutes, Δt: 9.524 minutes, CFL: 7.50e-01
i:    720, sim time: 7.140 days, wall time: 1.068 minutes, Δt: 10.609 minutes, CFL: 7.50e-01
i:    740, sim time: 7.277 days, wall time: 1.089 minutes, Δt: 9.414 minutes, CFL: 7.50e-01
i:    760, sim time: 7.408 days, wall time: 1.111 minutes, Δt: 9.364 minutes, CFL: 7.50e-01
i:    780, sim time: 7.525 days, wall time: 1.132 minutes, Δt: 9.126 minutes, CFL: 7.50e-01
i:    800, sim time: 7.648 days, wall time: 1.153 minutes, Δt: 9.257 minutes, CFL: 7.50e-01
i:    820, sim time: 7.769 days, wall time: 1.176 minutes, Δt: 8.825 minutes, CFL: 7.50e-01
i:    840, sim time: 7.887 days, wall time: 1.196 minutes, Δt: 8.660 minutes, CFL: 7.50e-01
i:    860, sim time: 8.006 days, wall time: 1.220 minutes, Δt: 9.337 minutes, CFL: 7.50e-01
i:    880, sim time: 8.140 days, wall time: 1.239 minutes, Δt: 10.339 minutes, CFL: 7.50e-01
i:    900, sim time: 8.279 days, wall time: 1.260 minutes, Δt: 10.321 minutes, CFL: 7.50e-01
i:    920, sim time: 8.417 days, wall time: 1.281 minutes, Δt: 10.204 minutes, CFL: 7.50e-01
i:    940, sim time: 8.557 days, wall time: 1.303 minutes, Δt: 10.549 minutes, CFL: 7.50e-01
i:    960, sim time: 8.697 days, wall time: 1.324 minutes, Δt: 10.967 minutes, CFL: 7.50e-01
i:    980, sim time: 8.849 days, wall time: 1.347 minutes, Δt: 11.350 minutes, CFL: 7.50e-01
i:   1000, sim time:     9 days, wall time: 1.367 minutes, Δt: 11.231 minutes, CFL: 7.50e-01
i:   1020, sim time: 9.156 days, wall time: 1.388 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.452 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.521 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.