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.458 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.691 seconds).
i: 20, sim time: 6 hours, wall time: 18.857 seconds, Δt: 25.164 minutes, CFL: 4.72e-01
i: 40, sim time: 15 hours, wall time: 20.256 seconds, Δt: 30 minutes, CFL: 5.87e-01
i: 60, sim time: 1.042 days, wall time: 21.589 seconds, Δt: 30 minutes, CFL: 6.08e-01
i: 80, sim time: 1.458 days, wall time: 22.907 seconds, Δt: 30 minutes, CFL: 6.40e-01
i: 100, sim time: 1.875 days, wall time: 24.218 seconds, Δt: 30 minutes, CFL: 6.70e-01
i: 120, sim time: 2.292 days, wall time: 25.552 seconds, Δt: 30 minutes, CFL: 7.32e-01
i: 140, sim time: 2.667 days, wall time: 26.736 seconds, Δt: 27.945 minutes, CFL: 7.50e-01
i: 160, sim time: 3 days, wall time: 28.047 seconds, Δt: 23.481 minutes, CFL: 7.50e-01
i: 180, sim time: 3.278 days, wall time: 29.491 seconds, Δt: 19.529 minutes, CFL: 7.50e-01
i: 200, sim time: 3.500 days, wall time: 30.661 seconds, Δt: 16.271 minutes, CFL: 7.50e-01
i: 220, sim time: 3.707 days, wall time: 31.967 seconds, Δt: 14.239 minutes, CFL: 7.50e-01
i: 240, sim time: 3.888 days, wall time: 33.265 seconds, Δt: 13.035 minutes, CFL: 7.50e-01
i: 260, sim time: 4.051 days, wall time: 34.556 seconds, Δt: 11.847 minutes, CFL: 7.50e-01
i: 280, sim time: 4.207 days, wall time: 35.852 seconds, Δt: 11.462 minutes, CFL: 7.50e-01
i: 300, sim time: 4.357 days, wall time: 37.198 seconds, Δt: 11.360 minutes, CFL: 7.50e-01
i: 320, sim time: 4.500 days, wall time: 38.429 seconds, Δt: 10.550 minutes, CFL: 7.50e-01
i: 340, sim time: 4.639 days, wall time: 39.731 seconds, Δt: 9.921 minutes, CFL: 7.50e-01
i: 360, sim time: 4.764 days, wall time: 41.145 seconds, Δt: 9.868 minutes, CFL: 7.50e-01
i: 380, sim time: 4.893 days, wall time: 42.319 seconds, Δt: 9.389 minutes, CFL: 7.50e-01
i: 400, sim time: 5.021 days, wall time: 43.671 seconds, Δt: 9.917 minutes, CFL: 7.50e-01
i: 420, sim time: 5.153 days, wall time: 44.923 seconds, Δt: 10.001 minutes, CFL: 7.50e-01
i: 440, sim time: 5.277 days, wall time: 46.225 seconds, Δt: 9.545 minutes, CFL: 7.50e-01
i: 460, sim time: 5.405 days, wall time: 47.501 seconds, Δt: 9.457 minutes, CFL: 7.50e-01
i: 480, sim time: 5.533 days, wall time: 48.784 seconds, Δt: 9.315 minutes, CFL: 7.50e-01
i: 500, sim time: 5.658 days, wall time: 50.058 seconds, Δt: 8.760 minutes, CFL: 7.50e-01
i: 520, sim time: 5.773 days, wall time: 51.352 seconds, Δt: 8.245 minutes, CFL: 7.50e-01
i: 540, sim time: 5.883 days, wall time: 52.640 seconds, Δt: 7.867 minutes, CFL: 7.50e-01
i: 560, sim time: 5.988 days, wall time: 53.920 seconds, Δt: 7.994 minutes, CFL: 7.50e-01
i: 580, sim time: 6.095 days, wall time: 55.324 seconds, Δt: 8.611 minutes, CFL: 7.50e-01
i: 600, sim time: 6.219 days, wall time: 56.496 seconds, Δt: 9.650 minutes, CFL: 7.50e-01
i: 620, sim time: 6.347 days, wall time: 57.915 seconds, Δt: 9.596 minutes, CFL: 7.50e-01
i: 640, sim time: 6.479 days, wall time: 59.092 seconds, Δt: 10.265 minutes, CFL: 7.50e-01
i: 660, sim time: 6.618 days, wall time: 1.007 minutes, Δt: 10.027 minutes, CFL: 7.50e-01
i: 680, sim time: 6.757 days, wall time: 1.031 minutes, Δt: 9.816 minutes, CFL: 7.50e-01
i: 700, sim time: 6.891 days, wall time: 1.050 minutes, Δt: 10.624 minutes, CFL: 7.50e-01
i: 720, sim time: 7.030 days, wall time: 1.072 minutes, Δt: 10.787 minutes, CFL: 7.50e-01
i: 740, sim time: 7.167 days, wall time: 1.093 minutes, Δt: 10.884 minutes, CFL: 7.50e-01
i: 760, sim time: 7.308 days, wall time: 1.115 minutes, Δt: 10.184 minutes, CFL: 7.50e-01
i: 780, sim time: 7.445 days, wall time: 1.136 minutes, Δt: 9.865 minutes, CFL: 7.50e-01
i: 800, sim time: 7.576 days, wall time: 1.158 minutes, Δt: 9.874 minutes, CFL: 7.50e-01
i: 820, sim time: 7.702 days, wall time: 1.179 minutes, Δt: 10.146 minutes, CFL: 7.50e-01
i: 840, sim time: 7.833 days, wall time: 1.201 minutes, Δt: 9.922 minutes, CFL: 7.50e-01
i: 860, sim time: 7.963 days, wall time: 1.222 minutes, Δt: 9.174 minutes, CFL: 7.50e-01
i: 880, sim time: 8.083 days, wall time: 1.244 minutes, Δt: 9.225 minutes, CFL: 7.50e-01
i: 900, sim time: 8.211 days, wall time: 1.265 minutes, Δt: 8.908 minutes, CFL: 7.50e-01
i: 920, sim time: 8.332 days, wall time: 1.287 minutes, Δt: 9.463 minutes, CFL: 7.50e-01
i: 940, sim time: 8.459 days, wall time: 1.308 minutes, Δt: 10.472 minutes, CFL: 7.50e-01
i: 960, sim time: 8.598 days, wall time: 1.332 minutes, Δt: 10.599 minutes, CFL: 7.50e-01
i: 980, sim time: 8.742 days, wall time: 1.352 minutes, Δt: 10.916 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.887 days, wall time: 1.374 minutes, Δt: 10.598 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.021 days, wall time: 1.396 minutes, Δt: 10.180 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.163 days, wall time: 1.417 minutes, Δt: 10.228 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.298 days, wall time: 1.439 minutes, Δt: 9.592 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.423 days, wall time: 1.464 minutes, Δt: 9.525 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.555 days, wall time: 1.482 minutes, Δt: 10.185 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.689 days, wall time: 1.505 minutes, Δt: 10.594 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.830 days, wall time: 1.525 minutes, Δt: 10.259 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.966 days, wall time: 1.547 minutes, Δt: 10.284 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.552 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.