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, 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.5Set 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, 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{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 entriesAdapt 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`sNamedTuple 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 CUDAGPUand 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=1.613e-23Periodically 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.139 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.771 seconds).
i: 20, sim time: 6 hours, wall time: 18.130 seconds, Δt: 25.164 minutes, CFL: 4.72e-01
i: 40, sim time: 15 hours, wall time: 19.284 seconds, Δt: 30 minutes, CFL: 5.87e-01
i: 60, sim time: 1.042 days, wall time: 20.341 seconds, Δt: 30 minutes, CFL: 6.08e-01
i: 80, sim time: 1.458 days, wall time: 21.396 seconds, Δt: 30 minutes, CFL: 6.40e-01
i: 100, sim time: 1.875 days, wall time: 22.460 seconds, Δt: 30 minutes, CFL: 6.70e-01
i: 120, sim time: 2.292 days, wall time: 23.525 seconds, Δt: 30 minutes, CFL: 7.32e-01
i: 140, sim time: 2.667 days, wall time: 24.483 seconds, Δt: 27.945 minutes, CFL: 7.50e-01
i: 160, sim time: 3 days, wall time: 25.544 seconds, Δt: 23.481 minutes, CFL: 7.50e-01
i: 180, sim time: 3.278 days, wall time: 26.694 seconds, Δt: 19.529 minutes, CFL: 7.50e-01
i: 200, sim time: 3.500 days, wall time: 27.649 seconds, Δt: 16.271 minutes, CFL: 7.50e-01
i: 220, sim time: 3.707 days, wall time: 28.705 seconds, Δt: 14.239 minutes, CFL: 7.50e-01
i: 240, sim time: 3.888 days, wall time: 29.755 seconds, Δt: 13.035 minutes, CFL: 7.50e-01
i: 260, sim time: 4.051 days, wall time: 30.803 seconds, Δt: 11.847 minutes, CFL: 7.50e-01
i: 280, sim time: 4.207 days, wall time: 31.854 seconds, Δt: 11.462 minutes, CFL: 7.50e-01
i: 300, sim time: 4.357 days, wall time: 32.952 seconds, Δt: 11.360 minutes, CFL: 7.50e-01
i: 320, sim time: 4.500 days, wall time: 33.949 seconds, Δt: 10.550 minutes, CFL: 7.50e-01
i: 340, sim time: 4.639 days, wall time: 34.997 seconds, Δt: 9.921 minutes, CFL: 7.50e-01
i: 360, sim time: 4.764 days, wall time: 36.140 seconds, Δt: 9.868 minutes, CFL: 7.50e-01
i: 380, sim time: 4.893 days, wall time: 37.091 seconds, Δt: 9.389 minutes, CFL: 7.50e-01
i: 400, sim time: 5.021 days, wall time: 38.193 seconds, Δt: 9.917 minutes, CFL: 7.50e-01
i: 420, sim time: 5.153 days, wall time: 39.194 seconds, Δt: 10.001 minutes, CFL: 7.50e-01
i: 440, sim time: 5.277 days, wall time: 40.246 seconds, Δt: 9.545 minutes, CFL: 7.50e-01
i: 460, sim time: 5.405 days, wall time: 41.281 seconds, Δt: 9.457 minutes, CFL: 7.50e-01
i: 480, sim time: 5.533 days, wall time: 42.325 seconds, Δt: 9.315 minutes, CFL: 7.50e-01
i: 500, sim time: 5.658 days, wall time: 43.363 seconds, Δt: 8.760 minutes, CFL: 7.50e-01
i: 520, sim time: 5.773 days, wall time: 44.415 seconds, Δt: 8.245 minutes, CFL: 7.50e-01
i: 540, sim time: 5.883 days, wall time: 45.449 seconds, Δt: 7.867 minutes, CFL: 7.50e-01
i: 560, sim time: 5.988 days, wall time: 46.492 seconds, Δt: 7.994 minutes, CFL: 7.50e-01
i: 580, sim time: 6.095 days, wall time: 47.639 seconds, Δt: 8.611 minutes, CFL: 7.50e-01
i: 600, sim time: 6.219 days, wall time: 48.584 seconds, Δt: 9.650 minutes, CFL: 7.50e-01
i: 620, sim time: 6.347 days, wall time: 49.729 seconds, Δt: 9.596 minutes, CFL: 7.50e-01
i: 640, sim time: 6.479 days, wall time: 50.676 seconds, Δt: 10.265 minutes, CFL: 7.50e-01
i: 660, sim time: 6.618 days, wall time: 51.728 seconds, Δt: 10.027 minutes, CFL: 7.50e-01
i: 680, sim time: 6.757 days, wall time: 52.912 seconds, Δt: 9.816 minutes, CFL: 7.50e-01
i: 700, sim time: 6.891 days, wall time: 53.816 seconds, Δt: 10.624 minutes, CFL: 7.50e-01
i: 720, sim time: 7.030 days, wall time: 54.865 seconds, Δt: 10.787 minutes, CFL: 7.50e-01
i: 740, sim time: 7.167 days, wall time: 55.902 seconds, Δt: 10.884 minutes, CFL: 7.50e-01
i: 760, sim time: 7.308 days, wall time: 56.951 seconds, Δt: 10.184 minutes, CFL: 7.50e-01
i: 780, sim time: 7.445 days, wall time: 58.000 seconds, Δt: 9.865 minutes, CFL: 7.50e-01
i: 800, sim time: 7.576 days, wall time: 59.032 seconds, Δt: 9.874 minutes, CFL: 7.50e-01
i: 820, sim time: 7.702 days, wall time: 1.001 minutes, Δt: 10.146 minutes, CFL: 7.50e-01
i: 840, sim time: 7.833 days, wall time: 1.019 minutes, Δt: 9.922 minutes, CFL: 7.50e-01
i: 860, sim time: 7.963 days, wall time: 1.036 minutes, Δt: 9.174 minutes, CFL: 7.50e-01
i: 880, sim time: 8.083 days, wall time: 1.054 minutes, Δt: 9.225 minutes, CFL: 7.50e-01
i: 900, sim time: 8.211 days, wall time: 1.071 minutes, Δt: 8.908 minutes, CFL: 7.50e-01
i: 920, sim time: 8.332 days, wall time: 1.088 minutes, Δt: 9.463 minutes, CFL: 7.50e-01
i: 940, sim time: 8.459 days, wall time: 1.106 minutes, Δt: 10.472 minutes, CFL: 7.50e-01
i: 960, sim time: 8.598 days, wall time: 1.125 minutes, Δt: 10.599 minutes, CFL: 7.50e-01
i: 980, sim time: 8.742 days, wall time: 1.140 minutes, Δt: 10.916 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.887 days, wall time: 1.158 minutes, Δt: 10.598 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.021 days, wall time: 1.176 minutes, Δt: 10.180 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.163 days, wall time: 1.192 minutes, Δt: 10.228 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.298 days, wall time: 1.210 minutes, Δt: 9.592 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.423 days, wall time: 1.230 minutes, Δt: 9.525 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.555 days, wall time: 1.245 minutes, Δt: 10.185 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.689 days, wall time: 1.263 minutes, Δt: 10.594 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.830 days, wall time: 1.279 minutes, Δt: 10.259 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.966 days, wall time: 1.297 minutes, Δt: 10.284 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.301 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
endThis page was generated using Literate.jl.