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.304 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.371 seconds).
i: 20, sim time: 6 hours, wall time: 23.159 seconds, Δt: 25.164 minutes, CFL: 4.70e-01
i: 40, sim time: 15 hours, wall time: 24.312 seconds, Δt: 30 minutes, CFL: 5.88e-01
i: 60, sim time: 1.042 days, wall time: 25.365 seconds, Δt: 30 minutes, CFL: 6.06e-01
i: 80, sim time: 1.458 days, wall time: 26.416 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 27.477 seconds, Δt: 30 minutes, CFL: 6.52e-01
i: 120, sim time: 2.292 days, wall time: 28.531 seconds, Δt: 30 minutes, CFL: 7.16e-01
i: 140, sim time: 2.687 days, wall time: 29.629 seconds, Δt: 28.626 minutes, CFL: 7.50e-01
i: 160, sim time: 3.018 days, wall time: 30.683 seconds, Δt: 24.961 minutes, CFL: 7.50e-01
i: 180, sim time: 3.323 days, wall time: 31.596 seconds, Δt: 19.356 minutes, CFL: 7.50e-01
i: 200, sim time: 3.559 days, wall time: 32.646 seconds, Δt: 16.277 minutes, CFL: 7.50e-01
i: 220, sim time: 3.750 days, wall time: 33.695 seconds, Δt: 14.399 minutes, CFL: 7.50e-01
i: 240, sim time: 3.936 days, wall time: 34.843 seconds, Δt: 13.104 minutes, CFL: 7.50e-01
i: 260, sim time: 4.101 days, wall time: 35.891 seconds, Δt: 13.037 minutes, CFL: 7.50e-01
i: 280, sim time: 4.268 days, wall time: 36.934 seconds, Δt: 12.875 minutes, CFL: 7.50e-01
i: 300, sim time: 4.434 days, wall time: 37.987 seconds, Δt: 12.278 minutes, CFL: 7.50e-01
i: 320, sim time: 4.600 days, wall time: 39.034 seconds, Δt: 12.097 minutes, CFL: 7.50e-01
i: 340, sim time: 4.767 days, wall time: 40.078 seconds, Δt: 11.919 minutes, CFL: 7.50e-01
i: 360, sim time: 4.933 days, wall time: 41.130 seconds, Δt: 11.734 minutes, CFL: 7.50e-01
i: 380, sim time: 5.083 days, wall time: 42.077 seconds, Δt: 11.608 minutes, CFL: 7.50e-01
i: 400, sim time: 5.241 days, wall time: 43.189 seconds, Δt: 12.021 minutes, CFL: 7.50e-01
i: 420, sim time: 5.404 days, wall time: 44.231 seconds, Δt: 13.007 minutes, CFL: 7.50e-01
i: 440, sim time: 5.569 days, wall time: 45.283 seconds, Δt: 12.059 minutes, CFL: 7.50e-01
i: 460, sim time: 5.722 days, wall time: 46.332 seconds, Δt: 11.125 minutes, CFL: 7.50e-01
i: 480, sim time: 5.863 days, wall time: 47.380 seconds, Δt: 10.535 minutes, CFL: 7.50e-01
i: 500, sim time: 6.000 days, wall time: 48.414 seconds, Δt: 9.654 minutes, CFL: 7.50e-01
i: 520, sim time: 6.120 days, wall time: 49.459 seconds, Δt: 8.777 minutes, CFL: 7.50e-01
i: 540, sim time: 6.237 days, wall time: 50.495 seconds, Δt: 8.368 minutes, CFL: 7.50e-01
i: 560, sim time: 6.345 days, wall time: 51.630 seconds, Δt: 8.016 minutes, CFL: 7.50e-01
i: 580, sim time: 6.449 days, wall time: 52.569 seconds, Δt: 7.549 minutes, CFL: 7.50e-01
i: 600, sim time: 6.550 days, wall time: 53.604 seconds, Δt: 6.965 minutes, CFL: 7.50e-01
i: 620, sim time: 6.644 days, wall time: 54.639 seconds, Δt: 6.712 minutes, CFL: 7.50e-01
i: 640, sim time: 6.735 days, wall time: 55.674 seconds, Δt: 6.509 minutes, CFL: 7.50e-01
i: 660, sim time: 6.822 days, wall time: 56.706 seconds, Δt: 6.519 minutes, CFL: 7.50e-01
i: 680, sim time: 6.913 days, wall time: 57.745 seconds, Δt: 7.131 minutes, CFL: 7.50e-01
i: 700, sim time: 7.010 days, wall time: 58.880 seconds, Δt: 7.463 minutes, CFL: 7.50e-01
i: 720, sim time: 7.115 days, wall time: 59.823 seconds, Δt: 7.581 minutes, CFL: 7.50e-01
i: 740, sim time: 7.221 days, wall time: 1.014 minutes, Δt: 7.896 minutes, CFL: 7.50e-01
i: 760, sim time: 7.329 days, wall time: 1.032 minutes, Δt: 8.367 minutes, CFL: 7.50e-01
i: 780, sim time: 7.440 days, wall time: 1.049 minutes, Δt: 8.879 minutes, CFL: 7.50e-01
i: 800, sim time: 7.563 days, wall time: 1.066 minutes, Δt: 9.111 minutes, CFL: 7.50e-01
i: 820, sim time: 7.679 days, wall time: 1.085 minutes, Δt: 8.773 minutes, CFL: 7.50e-01
i: 840, sim time: 7.797 days, wall time: 1.101 minutes, Δt: 8.530 minutes, CFL: 7.50e-01
i: 860, sim time: 7.913 days, wall time: 1.118 minutes, Δt: 9.154 minutes, CFL: 7.50e-01
i: 880, sim time: 8.030 days, wall time: 1.136 minutes, Δt: 8.279 minutes, CFL: 7.50e-01
i: 900, sim time: 8.138 days, wall time: 1.153 minutes, Δt: 7.768 minutes, CFL: 7.50e-01
i: 920, sim time: 8.240 days, wall time: 1.170 minutes, Δt: 7.486 minutes, CFL: 7.50e-01
i: 940, sim time: 8.344 days, wall time: 1.189 minutes, Δt: 7.766 minutes, CFL: 7.50e-01
i: 960, sim time: 8.450 days, wall time: 1.205 minutes, Δt: 8.009 minutes, CFL: 7.50e-01
i: 980, sim time: 8.565 days, wall time: 1.222 minutes, Δt: 8.932 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.686 days, wall time: 1.241 minutes, Δt: 9.549 minutes, CFL: 7.50e-01
i: 1020, sim time: 8.817 days, wall time: 1.257 minutes, Δt: 9.640 minutes, CFL: 7.50e-01
i: 1040, sim time: 8.944 days, wall time: 1.275 minutes, Δt: 9.679 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.071 days, wall time: 1.292 minutes, Δt: 9.128 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.190 days, wall time: 1.309 minutes, Δt: 8.391 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.302 days, wall time: 1.327 minutes, Δt: 8.244 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.415 days, wall time: 1.344 minutes, Δt: 8.567 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.530 days, wall time: 1.361 minutes, Δt: 8.817 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.651 days, wall time: 1.379 minutes, Δt: 8.921 minutes, CFL: 7.50e-01
i: 1180, sim time: 9.769 days, wall time: 1.397 minutes, Δt: 9.070 minutes, CFL: 7.50e-01
i: 1200, sim time: 9.890 days, wall time: 1.414 minutes, Δt: 9.073 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.430 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.