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.378 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.986 seconds).
i: 20, sim time: 6 hours, wall time: 18.724 seconds, Δt: 25.164 minutes, CFL: 4.72e-01
i: 40, sim time: 15 hours, wall time: 20.132 seconds, Δt: 30 minutes, CFL: 5.87e-01
i: 60, sim time: 1.042 days, wall time: 21.416 seconds, Δt: 30 minutes, CFL: 6.08e-01
i: 80, sim time: 1.458 days, wall time: 22.701 seconds, Δt: 30 minutes, CFL: 6.40e-01
i: 100, sim time: 1.875 days, wall time: 23.992 seconds, Δt: 30 minutes, CFL: 6.70e-01
i: 120, sim time: 2.292 days, wall time: 25.290 seconds, Δt: 30 minutes, CFL: 7.32e-01
i: 140, sim time: 2.667 days, wall time: 26.462 seconds, Δt: 27.945 minutes, CFL: 7.50e-01
i: 160, sim time: 3 days, wall time: 27.749 seconds, Δt: 23.481 minutes, CFL: 7.50e-01
i: 180, sim time: 3.278 days, wall time: 29.151 seconds, Δt: 19.529 minutes, CFL: 7.50e-01
i: 200, sim time: 3.500 days, wall time: 30.321 seconds, Δt: 16.271 minutes, CFL: 7.50e-01
i: 220, sim time: 3.707 days, wall time: 31.612 seconds, Δt: 14.239 minutes, CFL: 7.50e-01
i: 240, sim time: 3.888 days, wall time: 32.892 seconds, Δt: 13.035 minutes, CFL: 7.50e-01
i: 260, sim time: 4.051 days, wall time: 34.176 seconds, Δt: 11.847 minutes, CFL: 7.50e-01
i: 280, sim time: 4.207 days, wall time: 35.459 seconds, Δt: 11.462 minutes, CFL: 7.50e-01
i: 300, sim time: 4.357 days, wall time: 36.800 seconds, Δt: 11.360 minutes, CFL: 7.50e-01
i: 320, sim time: 4.500 days, wall time: 38.019 seconds, Δt: 10.550 minutes, CFL: 7.50e-01
i: 340, sim time: 4.639 days, wall time: 39.297 seconds, Δt: 9.921 minutes, CFL: 7.50e-01
i: 360, sim time: 4.764 days, wall time: 40.688 seconds, Δt: 9.868 minutes, CFL: 7.50e-01
i: 380, sim time: 4.893 days, wall time: 41.855 seconds, Δt: 9.389 minutes, CFL: 7.50e-01
i: 400, sim time: 5.021 days, wall time: 43.203 seconds, Δt: 9.917 minutes, CFL: 7.50e-01
i: 420, sim time: 5.153 days, wall time: 44.429 seconds, Δt: 10.001 minutes, CFL: 7.50e-01
i: 440, sim time: 5.277 days, wall time: 45.719 seconds, Δt: 9.545 minutes, CFL: 7.50e-01
i: 460, sim time: 5.405 days, wall time: 47.001 seconds, Δt: 9.457 minutes, CFL: 7.50e-01
i: 480, sim time: 5.533 days, wall time: 48.288 seconds, Δt: 9.315 minutes, CFL: 7.50e-01
i: 500, sim time: 5.658 days, wall time: 49.561 seconds, Δt: 8.760 minutes, CFL: 7.50e-01
i: 520, sim time: 5.773 days, wall time: 50.849 seconds, Δt: 8.245 minutes, CFL: 7.50e-01
i: 540, sim time: 5.883 days, wall time: 52.134 seconds, Δt: 7.867 minutes, CFL: 7.50e-01
i: 560, sim time: 5.988 days, wall time: 53.414 seconds, Δt: 7.994 minutes, CFL: 7.50e-01
i: 580, sim time: 6.095 days, wall time: 54.816 seconds, Δt: 8.611 minutes, CFL: 7.50e-01
i: 600, sim time: 6.219 days, wall time: 55.982 seconds, Δt: 9.650 minutes, CFL: 7.50e-01
i: 620, sim time: 6.347 days, wall time: 57.376 seconds, Δt: 9.596 minutes, CFL: 7.50e-01
i: 640, sim time: 6.479 days, wall time: 58.544 seconds, Δt: 10.265 minutes, CFL: 7.50e-01
i: 660, sim time: 6.618 days, wall time: 59.829 seconds, Δt: 10.027 minutes, CFL: 7.50e-01
i: 680, sim time: 6.757 days, wall time: 1.021 minutes, Δt: 9.816 minutes, CFL: 7.50e-01
i: 700, sim time: 6.891 days, wall time: 1.040 minutes, Δt: 10.624 minutes, CFL: 7.50e-01
i: 720, sim time: 7.030 days, wall time: 1.061 minutes, Δt: 10.787 minutes, CFL: 7.50e-01
i: 740, sim time: 7.167 days, wall time: 1.083 minutes, Δt: 10.884 minutes, CFL: 7.50e-01
i: 760, sim time: 7.308 days, wall time: 1.104 minutes, Δt: 10.184 minutes, CFL: 7.50e-01
i: 780, sim time: 7.445 days, wall time: 1.126 minutes, Δt: 9.865 minutes, CFL: 7.50e-01
i: 800, sim time: 7.576 days, wall time: 1.147 minutes, Δt: 9.874 minutes, CFL: 7.50e-01
i: 820, sim time: 7.702 days, wall time: 1.168 minutes, Δt: 10.146 minutes, CFL: 7.50e-01
i: 840, sim time: 7.833 days, wall time: 1.190 minutes, Δt: 9.922 minutes, CFL: 7.50e-01
i: 860, sim time: 7.963 days, wall time: 1.211 minutes, Δt: 9.174 minutes, CFL: 7.50e-01
i: 880, sim time: 8.083 days, wall time: 1.232 minutes, Δt: 9.225 minutes, CFL: 7.50e-01
i: 900, sim time: 8.211 days, wall time: 1.254 minutes, Δt: 8.908 minutes, CFL: 7.50e-01
i: 920, sim time: 8.332 days, wall time: 1.275 minutes, Δt: 9.463 minutes, CFL: 7.50e-01
i: 940, sim time: 8.459 days, wall time: 1.296 minutes, Δt: 10.472 minutes, CFL: 7.50e-01
i: 960, sim time: 8.598 days, wall time: 1.320 minutes, Δt: 10.599 minutes, CFL: 7.50e-01
i: 980, sim time: 8.742 days, wall time: 1.339 minutes, Δt: 10.916 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.887 days, wall time: 1.361 minutes, Δt: 10.598 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.021 days, wall time: 1.383 minutes, Δt: 10.180 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.163 days, wall time: 1.404 minutes, Δt: 10.228 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.298 days, wall time: 1.425 minutes, Δt: 9.592 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.423 days, wall time: 1.450 minutes, Δt: 9.525 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.555 days, wall time: 1.468 minutes, Δt: 10.185 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.689 days, wall time: 1.491 minutes, Δt: 10.594 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.830 days, wall time: 1.511 minutes, Δt: 10.259 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.966 days, wall time: 1.533 minutes, Δt: 10.284 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.538 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.