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.155 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (5.060 seconds).
i: 20, sim time: 6 hours, wall time: 18.399 seconds, Δt: 25.164 minutes, CFL: 4.72e-01
i: 40, sim time: 15 hours, wall time: 19.549 seconds, Δt: 30 minutes, CFL: 5.87e-01
i: 60, sim time: 1.042 days, wall time: 20.610 seconds, Δt: 30 minutes, CFL: 6.08e-01
i: 80, sim time: 1.458 days, wall time: 21.662 seconds, Δt: 30 minutes, CFL: 6.40e-01
i: 100, sim time: 1.875 days, wall time: 22.720 seconds, Δt: 30 minutes, CFL: 6.70e-01
i: 120, sim time: 2.292 days, wall time: 23.779 seconds, Δt: 30 minutes, CFL: 7.32e-01
i: 140, sim time: 2.667 days, wall time: 24.739 seconds, Δt: 27.945 minutes, CFL: 7.50e-01
i: 160, sim time: 3 days, wall time: 25.798 seconds, Δt: 23.481 minutes, CFL: 7.50e-01
i: 180, sim time: 3.278 days, wall time: 26.945 seconds, Δt: 19.529 minutes, CFL: 7.50e-01
i: 200, sim time: 3.500 days, wall time: 27.906 seconds, Δt: 16.271 minutes, CFL: 7.50e-01
i: 220, sim time: 3.707 days, wall time: 28.970 seconds, Δt: 14.239 minutes, CFL: 7.50e-01
i: 240, sim time: 3.888 days, wall time: 30.018 seconds, Δt: 13.035 minutes, CFL: 7.50e-01
i: 260, sim time: 4.051 days, wall time: 31.065 seconds, Δt: 11.847 minutes, CFL: 7.50e-01
i: 280, sim time: 4.207 days, wall time: 32.113 seconds, Δt: 11.462 minutes, CFL: 7.50e-01
i: 300, sim time: 4.357 days, wall time: 33.212 seconds, Δt: 11.360 minutes, CFL: 7.50e-01
i: 320, sim time: 4.500 days, wall time: 34.213 seconds, Δt: 10.550 minutes, CFL: 7.50e-01
i: 340, sim time: 4.639 days, wall time: 35.263 seconds, Δt: 9.921 minutes, CFL: 7.50e-01
i: 360, sim time: 4.764 days, wall time: 36.405 seconds, Δt: 9.868 minutes, CFL: 7.50e-01
i: 380, sim time: 4.893 days, wall time: 37.355 seconds, Δt: 9.389 minutes, CFL: 7.50e-01
i: 400, sim time: 5.021 days, wall time: 38.451 seconds, Δt: 9.917 minutes, CFL: 7.50e-01
i: 420, sim time: 5.153 days, wall time: 39.444 seconds, Δt: 10.001 minutes, CFL: 7.50e-01
i: 440, sim time: 5.277 days, wall time: 40.494 seconds, Δt: 9.545 minutes, CFL: 7.50e-01
i: 460, sim time: 5.405 days, wall time: 41.532 seconds, Δt: 9.457 minutes, CFL: 7.50e-01
i: 480, sim time: 5.533 days, wall time: 42.578 seconds, Δt: 9.315 minutes, CFL: 7.50e-01
i: 500, sim time: 5.658 days, wall time: 43.619 seconds, Δt: 8.760 minutes, CFL: 7.50e-01
i: 520, sim time: 5.773 days, wall time: 44.668 seconds, Δt: 8.245 minutes, CFL: 7.50e-01
i: 540, sim time: 5.883 days, wall time: 45.712 seconds, Δt: 7.867 minutes, CFL: 7.50e-01
i: 560, sim time: 5.988 days, wall time: 46.762 seconds, Δt: 7.994 minutes, CFL: 7.50e-01
i: 580, sim time: 6.095 days, wall time: 47.914 seconds, Δt: 8.611 minutes, CFL: 7.50e-01
i: 600, sim time: 6.219 days, wall time: 48.858 seconds, Δt: 9.650 minutes, CFL: 7.50e-01
i: 620, sim time: 6.347 days, wall time: 50.006 seconds, Δt: 9.596 minutes, CFL: 7.50e-01
i: 640, sim time: 6.479 days, wall time: 50.953 seconds, Δt: 10.265 minutes, CFL: 7.50e-01
i: 660, sim time: 6.618 days, wall time: 52.001 seconds, Δt: 10.027 minutes, CFL: 7.50e-01
i: 680, sim time: 6.757 days, wall time: 53.186 seconds, Δt: 9.816 minutes, CFL: 7.50e-01
i: 700, sim time: 6.891 days, wall time: 54.089 seconds, Δt: 10.624 minutes, CFL: 7.50e-01
i: 720, sim time: 7.030 days, wall time: 55.140 seconds, Δt: 10.787 minutes, CFL: 7.50e-01
i: 740, sim time: 7.167 days, wall time: 56.180 seconds, Δt: 10.884 minutes, CFL: 7.50e-01
i: 760, sim time: 7.308 days, wall time: 57.220 seconds, Δt: 10.184 minutes, CFL: 7.50e-01
i: 780, sim time: 7.445 days, wall time: 58.270 seconds, Δt: 9.865 minutes, CFL: 7.50e-01
i: 800, sim time: 7.576 days, wall time: 59.306 seconds, Δt: 9.874 minutes, CFL: 7.50e-01
i: 820, sim time: 7.702 days, wall time: 1.006 minutes, Δt: 10.146 minutes, CFL: 7.50e-01
i: 840, sim time: 7.833 days, wall time: 1.023 minutes, Δt: 9.922 minutes, CFL: 7.50e-01
i: 860, sim time: 7.963 days, wall time: 1.041 minutes, Δt: 9.174 minutes, CFL: 7.50e-01
i: 880, sim time: 8.083 days, wall time: 1.058 minutes, Δt: 9.225 minutes, CFL: 7.50e-01
i: 900, sim time: 8.211 days, wall time: 1.075 minutes, Δt: 8.908 minutes, CFL: 7.50e-01
i: 920, sim time: 8.332 days, wall time: 1.093 minutes, Δt: 9.463 minutes, CFL: 7.50e-01
i: 940, sim time: 8.459 days, wall time: 1.110 minutes, Δt: 10.472 minutes, CFL: 7.50e-01
i: 960, sim time: 8.598 days, wall time: 1.129 minutes, Δt: 10.599 minutes, CFL: 7.50e-01
i: 980, sim time: 8.742 days, wall time: 1.145 minutes, Δt: 10.916 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.887 days, wall time: 1.162 minutes, Δt: 10.598 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.021 days, wall time: 1.181 minutes, Δt: 10.180 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.163 days, wall time: 1.197 minutes, Δt: 10.228 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.298 days, wall time: 1.214 minutes, Δt: 9.592 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.423 days, wall time: 1.234 minutes, Δt: 9.525 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.555 days, wall time: 1.249 minutes, Δt: 10.185 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.689 days, wall time: 1.267 minutes, Δt: 10.594 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.830 days, wall time: 1.284 minutes, Δt: 10.259 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.966 days, wall time: 1.301 minutes, Δt: 10.284 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.305 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.