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.445 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.919 seconds).
i: 20, sim time: 6 hours, wall time: 20.130 seconds, Δt: 25.164 minutes, CFL: 4.72e-01
i: 40, sim time: 15 hours, wall time: 23.111 seconds, Δt: 30 minutes, CFL: 5.87e-01
i: 60, sim time: 1.042 days, wall time: 25.848 seconds, Δt: 30 minutes, CFL: 6.08e-01
i: 80, sim time: 1.458 days, wall time: 28.342 seconds, Δt: 30 minutes, CFL: 6.40e-01
i: 100, sim time: 1.875 days, wall time: 31.076 seconds, Δt: 30 minutes, CFL: 6.70e-01
i: 120, sim time: 2.292 days, wall time: 33.816 seconds, Δt: 30 minutes, CFL: 7.32e-01
i: 140, sim time: 2.667 days, wall time: 36.318 seconds, Δt: 27.945 minutes, CFL: 7.50e-01
i: 160, sim time: 3 days, wall time: 39.058 seconds, Δt: 23.481 minutes, CFL: 7.50e-01
i: 180, sim time: 3.278 days, wall time: 42.029 seconds, Δt: 19.529 minutes, CFL: 7.50e-01
i: 200, sim time: 3.500 days, wall time: 44.541 seconds, Δt: 16.271 minutes, CFL: 7.50e-01
i: 220, sim time: 3.707 days, wall time: 47.292 seconds, Δt: 14.239 minutes, CFL: 7.50e-01
i: 240, sim time: 3.888 days, wall time: 50.036 seconds, Δt: 13.035 minutes, CFL: 7.50e-01
i: 260, sim time: 4.051 days, wall time: 52.790 seconds, Δt: 11.847 minutes, CFL: 7.50e-01
i: 280, sim time: 4.207 days, wall time: 55.270 seconds, Δt: 11.462 minutes, CFL: 7.50e-01
i: 300, sim time: 4.357 days, wall time: 58.117 seconds, Δt: 11.360 minutes, CFL: 7.50e-01
i: 320, sim time: 4.500 days, wall time: 1.013 minutes, Δt: 10.550 minutes, CFL: 7.50e-01
i: 340, sim time: 4.639 days, wall time: 1.058 minutes, Δt: 9.921 minutes, CFL: 7.50e-01
i: 360, sim time: 4.764 days, wall time: 1.108 minutes, Δt: 9.868 minutes, CFL: 7.50e-01
i: 380, sim time: 4.893 days, wall time: 1.149 minutes, Δt: 9.389 minutes, CFL: 7.50e-01
i: 400, sim time: 5.021 days, wall time: 1.197 minutes, Δt: 9.917 minutes, CFL: 7.50e-01
i: 420, sim time: 5.153 days, wall time: 1.241 minutes, Δt: 10.001 minutes, CFL: 7.50e-01
i: 440, sim time: 5.277 days, wall time: 1.286 minutes, Δt: 9.545 minutes, CFL: 7.50e-01
i: 460, sim time: 5.405 days, wall time: 1.331 minutes, Δt: 9.457 minutes, CFL: 7.50e-01
i: 480, sim time: 5.533 days, wall time: 1.372 minutes, Δt: 9.315 minutes, CFL: 7.50e-01
i: 500, sim time: 5.658 days, wall time: 1.418 minutes, Δt: 8.760 minutes, CFL: 7.50e-01
i: 520, sim time: 5.773 days, wall time: 1.463 minutes, Δt: 8.245 minutes, CFL: 7.50e-01
i: 540, sim time: 5.883 days, wall time: 1.508 minutes, Δt: 7.867 minutes, CFL: 7.50e-01
i: 560, sim time: 5.988 days, wall time: 1.554 minutes, Δt: 7.994 minutes, CFL: 7.50e-01
i: 580, sim time: 6.095 days, wall time: 1.603 minutes, Δt: 8.611 minutes, CFL: 7.50e-01
i: 600, sim time: 6.219 days, wall time: 1.645 minutes, Δt: 9.650 minutes, CFL: 7.50e-01
i: 620, sim time: 6.347 days, wall time: 1.694 minutes, Δt: 9.596 minutes, CFL: 7.50e-01
i: 640, sim time: 6.479 days, wall time: 1.736 minutes, Δt: 10.265 minutes, CFL: 7.50e-01
i: 660, sim time: 6.618 days, wall time: 1.782 minutes, Δt: 10.027 minutes, CFL: 7.50e-01
i: 680, sim time: 6.757 days, wall time: 1.829 minutes, Δt: 9.816 minutes, CFL: 7.50e-01
i: 700, sim time: 6.891 days, wall time: 1.869 minutes, Δt: 10.624 minutes, CFL: 7.50e-01
i: 720, sim time: 7.030 days, wall time: 1.915 minutes, Δt: 10.787 minutes, CFL: 7.50e-01
i: 740, sim time: 7.167 days, wall time: 1.960 minutes, Δt: 10.884 minutes, CFL: 7.50e-01
i: 760, sim time: 7.308 days, wall time: 2.006 minutes, Δt: 10.184 minutes, CFL: 7.50e-01
i: 780, sim time: 7.445 days, wall time: 2.051 minutes, Δt: 9.865 minutes, CFL: 7.50e-01
i: 800, sim time: 7.576 days, wall time: 2.096 minutes, Δt: 9.874 minutes, CFL: 7.50e-01
i: 820, sim time: 7.702 days, wall time: 2.142 minutes, Δt: 10.146 minutes, CFL: 7.50e-01
i: 840, sim time: 7.833 days, wall time: 2.187 minutes, Δt: 9.922 minutes, CFL: 7.50e-01
i: 860, sim time: 7.963 days, wall time: 2.233 minutes, Δt: 9.174 minutes, CFL: 7.50e-01
i: 880, sim time: 8.083 days, wall time: 2.278 minutes, Δt: 9.225 minutes, CFL: 7.50e-01
i: 900, sim time: 8.211 days, wall time: 2.319 minutes, Δt: 8.908 minutes, CFL: 7.50e-01
i: 920, sim time: 8.332 days, wall time: 2.364 minutes, Δt: 9.463 minutes, CFL: 7.50e-01
i: 940, sim time: 8.459 days, wall time: 2.410 minutes, Δt: 10.472 minutes, CFL: 7.50e-01
i: 960, sim time: 8.598 days, wall time: 2.459 minutes, Δt: 10.599 minutes, CFL: 7.50e-01
i: 980, sim time: 8.742 days, wall time: 2.500 minutes, Δt: 10.916 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.887 days, wall time: 2.546 minutes, Δt: 10.598 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.021 days, wall time: 2.593 minutes, Δt: 10.180 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.163 days, wall time: 2.636 minutes, Δt: 10.228 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.298 days, wall time: 2.682 minutes, Δt: 9.592 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.423 days, wall time: 2.733 minutes, Δt: 9.525 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.555 days, wall time: 2.772 minutes, Δt: 10.185 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.689 days, wall time: 2.815 minutes, Δt: 10.594 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.830 days, wall time: 2.859 minutes, Δt: 10.259 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.966 days, wall time: 2.904 minutes, Δt: 10.284 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 2.916 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.