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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.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{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.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{GPU, 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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.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=6.41337e-5, min=-5.24045e-5, mean=-8.27181e-24Periodically 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.25e-01
[ Info: ... simulation initialization complete (17.236 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.767 seconds).
i: 20, sim time: 6 hours, wall time: 23.762 seconds, Δt: 25.164 minutes, CFL: 4.70e-01
i: 40, sim time: 15 hours, wall time: 25.668 seconds, Δt: 30 minutes, CFL: 5.88e-01
i: 60, sim time: 1.042 days, wall time: 27.419 seconds, Δt: 30 minutes, CFL: 6.06e-01
i: 80, sim time: 1.458 days, wall time: 29.176 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 30.933 seconds, Δt: 30 minutes, CFL: 6.52e-01
i: 120, sim time: 2.292 days, wall time: 32.683 seconds, Δt: 30 minutes, CFL: 7.16e-01
i: 140, sim time: 2.687 days, wall time: 34.521 seconds, Δt: 28.626 minutes, CFL: 7.50e-01
i: 160, sim time: 3.018 days, wall time: 36.276 seconds, Δt: 24.961 minutes, CFL: 7.50e-01
i: 180, sim time: 3.323 days, wall time: 37.797 seconds, Δt: 19.356 minutes, CFL: 7.50e-01
i: 200, sim time: 3.559 days, wall time: 39.551 seconds, Δt: 16.277 minutes, CFL: 7.50e-01
i: 220, sim time: 3.750 days, wall time: 41.299 seconds, Δt: 14.399 minutes, CFL: 7.50e-01
i: 240, sim time: 3.936 days, wall time: 43.208 seconds, Δt: 13.104 minutes, CFL: 7.50e-01
i: 260, sim time: 4.101 days, wall time: 44.957 seconds, Δt: 13.037 minutes, CFL: 7.50e-01
i: 280, sim time: 4.268 days, wall time: 46.712 seconds, Δt: 12.875 minutes, CFL: 7.50e-01
i: 300, sim time: 4.434 days, wall time: 48.470 seconds, Δt: 12.278 minutes, CFL: 7.50e-01
i: 320, sim time: 4.600 days, wall time: 50.220 seconds, Δt: 12.097 minutes, CFL: 7.50e-01
i: 340, sim time: 4.767 days, wall time: 51.978 seconds, Δt: 11.919 minutes, CFL: 7.50e-01
i: 360, sim time: 4.933 days, wall time: 53.728 seconds, Δt: 11.734 minutes, CFL: 7.50e-01
i: 380, sim time: 5.083 days, wall time: 55.321 seconds, Δt: 11.608 minutes, CFL: 7.50e-01
i: 400, sim time: 5.241 days, wall time: 57.079 seconds, Δt: 12.021 minutes, CFL: 7.50e-01
i: 420, sim time: 5.404 days, wall time: 58.837 seconds, Δt: 13.007 minutes, CFL: 7.50e-01
i: 440, sim time: 5.569 days, wall time: 1.010 minutes, Δt: 12.059 minutes, CFL: 7.50e-01
i: 460, sim time: 5.722 days, wall time: 1.039 minutes, Δt: 11.125 minutes, CFL: 7.50e-01
i: 480, sim time: 5.863 days, wall time: 1.068 minutes, Δt: 10.535 minutes, CFL: 7.50e-01
i: 500, sim time: 6.000 days, wall time: 1.097 minutes, Δt: 9.654 minutes, CFL: 7.50e-01
i: 520, sim time: 6.120 days, wall time: 1.126 minutes, Δt: 8.777 minutes, CFL: 7.50e-01
i: 540, sim time: 6.237 days, wall time: 1.155 minutes, Δt: 8.368 minutes, CFL: 7.50e-01
i: 560, sim time: 6.345 days, wall time: 1.187 minutes, Δt: 8.016 minutes, CFL: 7.50e-01
i: 580, sim time: 6.449 days, wall time: 1.214 minutes, Δt: 7.549 minutes, CFL: 7.50e-01
i: 600, sim time: 6.550 days, wall time: 1.243 minutes, Δt: 6.965 minutes, CFL: 7.50e-01
i: 620, sim time: 6.644 days, wall time: 1.272 minutes, Δt: 6.712 minutes, CFL: 7.50e-01
i: 640, sim time: 6.735 days, wall time: 1.301 minutes, Δt: 6.509 minutes, CFL: 7.50e-01
i: 660, sim time: 6.822 days, wall time: 1.330 minutes, Δt: 6.519 minutes, CFL: 7.50e-01
i: 680, sim time: 6.913 days, wall time: 1.359 minutes, Δt: 7.131 minutes, CFL: 7.50e-01
i: 700, sim time: 7.010 days, wall time: 1.391 minutes, Δt: 7.463 minutes, CFL: 7.50e-01
i: 720, sim time: 7.115 days, wall time: 1.417 minutes, Δt: 7.581 minutes, CFL: 7.50e-01
i: 740, sim time: 7.221 days, wall time: 1.447 minutes, Δt: 7.896 minutes, CFL: 7.50e-01
i: 760, sim time: 7.329 days, wall time: 1.476 minutes, Δt: 8.367 minutes, CFL: 7.50e-01
i: 780, sim time: 7.440 days, wall time: 1.505 minutes, Δt: 8.879 minutes, CFL: 7.50e-01
i: 800, sim time: 7.563 days, wall time: 1.535 minutes, Δt: 9.111 minutes, CFL: 7.50e-01
i: 820, sim time: 7.679 days, wall time: 1.567 minutes, Δt: 8.773 minutes, CFL: 7.50e-01
i: 840, sim time: 7.797 days, wall time: 1.593 minutes, Δt: 8.530 minutes, CFL: 7.50e-01
i: 860, sim time: 7.913 days, wall time: 1.623 minutes, Δt: 9.154 minutes, CFL: 7.50e-01
i: 880, sim time: 8.030 days, wall time: 1.652 minutes, Δt: 8.279 minutes, CFL: 7.50e-01
i: 900, sim time: 8.138 days, wall time: 1.681 minutes, Δt: 7.768 minutes, CFL: 7.50e-01
i: 920, sim time: 8.240 days, wall time: 1.711 minutes, Δt: 7.486 minutes, CFL: 7.50e-01
i: 940, sim time: 8.344 days, wall time: 1.742 minutes, Δt: 7.766 minutes, CFL: 7.50e-01
i: 960, sim time: 8.450 days, wall time: 1.769 minutes, Δt: 8.009 minutes, CFL: 7.50e-01
i: 980, sim time: 8.565 days, wall time: 1.799 minutes, Δt: 8.932 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.686 days, wall time: 1.829 minutes, Δt: 9.549 minutes, CFL: 7.50e-01
i: 1020, sim time: 8.817 days, wall time: 1.857 minutes, Δt: 9.640 minutes, CFL: 7.50e-01
i: 1040, sim time: 8.944 days, wall time: 1.887 minutes, Δt: 9.679 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.071 days, wall time: 1.916 minutes, Δt: 9.128 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.190 days, wall time: 1.945 minutes, Δt: 8.391 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.302 days, wall time: 1.975 minutes, Δt: 8.244 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.415 days, wall time: 2.004 minutes, Δt: 8.567 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.530 days, wall time: 2.033 minutes, Δt: 8.817 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.651 days, wall time: 2.063 minutes, Δt: 8.921 minutes, CFL: 7.50e-01
i: 1180, sim time: 9.769 days, wall time: 2.093 minutes, Δt: 9.070 minutes, CFL: 7.50e-01
i: 1200, sim time: 9.890 days, wall time: 2.122 minutes, Δt: 9.073 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 2.149 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.