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.113 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.668 seconds).
i: 20, sim time: 6 hours, wall time: 22.952 seconds, Δt: 25.164 minutes, CFL: 4.70e-01
i: 40, sim time: 15 hours, wall time: 24.111 seconds, Δt: 30 minutes, CFL: 5.88e-01
i: 60, sim time: 1.042 days, wall time: 25.172 seconds, Δt: 30 minutes, CFL: 6.06e-01
i: 80, sim time: 1.458 days, wall time: 26.232 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 27.281 seconds, Δt: 30 minutes, CFL: 6.52e-01
i: 120, sim time: 2.292 days, wall time: 28.343 seconds, Δt: 30 minutes, CFL: 7.16e-01
i: 140, sim time: 2.687 days, wall time: 29.442 seconds, Δt: 28.626 minutes, CFL: 7.50e-01
i: 160, sim time: 3.018 days, wall time: 30.500 seconds, Δt: 24.961 minutes, CFL: 7.50e-01
i: 180, sim time: 3.323 days, wall time: 31.414 seconds, Δt: 19.356 minutes, CFL: 7.50e-01
i: 200, sim time: 3.559 days, wall time: 32.465 seconds, Δt: 16.277 minutes, CFL: 7.50e-01
i: 220, sim time: 3.750 days, wall time: 33.512 seconds, Δt: 14.399 minutes, CFL: 7.50e-01
i: 240, sim time: 3.936 days, wall time: 34.659 seconds, Δt: 13.104 minutes, CFL: 7.50e-01
i: 260, sim time: 4.101 days, wall time: 35.705 seconds, Δt: 13.037 minutes, CFL: 7.50e-01
i: 280, sim time: 4.268 days, wall time: 36.758 seconds, Δt: 12.875 minutes, CFL: 7.50e-01
i: 300, sim time: 4.434 days, wall time: 37.801 seconds, Δt: 12.278 minutes, CFL: 7.50e-01
i: 320, sim time: 4.600 days, wall time: 38.847 seconds, Δt: 12.097 minutes, CFL: 7.50e-01
i: 340, sim time: 4.767 days, wall time: 39.894 seconds, Δt: 11.919 minutes, CFL: 7.50e-01
i: 360, sim time: 4.933 days, wall time: 40.942 seconds, Δt: 11.734 minutes, CFL: 7.50e-01
i: 380, sim time: 5.083 days, wall time: 41.887 seconds, Δt: 11.608 minutes, CFL: 7.50e-01
i: 400, sim time: 5.241 days, wall time: 42.932 seconds, Δt: 12.021 minutes, CFL: 7.50e-01
i: 420, sim time: 5.404 days, wall time: 43.984 seconds, Δt: 13.007 minutes, CFL: 7.50e-01
i: 440, sim time: 5.569 days, wall time: 45.028 seconds, Δt: 12.059 minutes, CFL: 7.50e-01
i: 460, sim time: 5.722 days, wall time: 46.074 seconds, Δt: 11.125 minutes, CFL: 7.50e-01
i: 480, sim time: 5.863 days, wall time: 47.122 seconds, Δt: 10.535 minutes, CFL: 7.50e-01
i: 500, sim time: 6.000 days, wall time: 48.159 seconds, Δt: 9.654 minutes, CFL: 7.50e-01
i: 520, sim time: 6.120 days, wall time: 49.196 seconds, Δt: 8.777 minutes, CFL: 7.50e-01
i: 540, sim time: 6.237 days, wall time: 50.241 seconds, Δt: 8.368 minutes, CFL: 7.50e-01
i: 560, sim time: 6.345 days, wall time: 51.376 seconds, Δt: 8.016 minutes, CFL: 7.50e-01
i: 580, sim time: 6.449 days, wall time: 52.319 seconds, Δt: 7.549 minutes, CFL: 7.50e-01
i: 600, sim time: 6.550 days, wall time: 53.356 seconds, Δt: 6.965 minutes, CFL: 7.50e-01
i: 620, sim time: 6.644 days, wall time: 54.393 seconds, Δt: 6.712 minutes, CFL: 7.50e-01
i: 640, sim time: 6.735 days, wall time: 55.431 seconds, Δt: 6.509 minutes, CFL: 7.50e-01
i: 660, sim time: 6.822 days, wall time: 56.472 seconds, Δt: 6.519 minutes, CFL: 7.50e-01
i: 680, sim time: 6.913 days, wall time: 57.515 seconds, Δt: 7.131 minutes, CFL: 7.50e-01
i: 700, sim time: 7.010 days, wall time: 58.654 seconds, Δt: 7.463 minutes, CFL: 7.50e-01
i: 720, sim time: 7.115 days, wall time: 59.600 seconds, Δt: 7.581 minutes, CFL: 7.50e-01
i: 740, sim time: 7.221 days, wall time: 1.011 minutes, Δt: 7.896 minutes, CFL: 7.50e-01
i: 760, sim time: 7.329 days, wall time: 1.028 minutes, Δt: 8.367 minutes, CFL: 7.50e-01
i: 780, sim time: 7.440 days, wall time: 1.046 minutes, Δt: 8.879 minutes, CFL: 7.50e-01
i: 800, sim time: 7.563 days, wall time: 1.063 minutes, Δt: 9.111 minutes, CFL: 7.50e-01
i: 820, sim time: 7.679 days, wall time: 1.082 minutes, Δt: 8.773 minutes, CFL: 7.50e-01
i: 840, sim time: 7.797 days, wall time: 1.098 minutes, Δt: 8.530 minutes, CFL: 7.50e-01
i: 860, sim time: 7.913 days, wall time: 1.115 minutes, Δt: 9.154 minutes, CFL: 7.50e-01
i: 880, sim time: 8.030 days, wall time: 1.133 minutes, Δt: 8.279 minutes, CFL: 7.50e-01
i: 900, sim time: 8.138 days, wall time: 1.150 minutes, Δt: 7.768 minutes, CFL: 7.50e-01
i: 920, sim time: 8.240 days, wall time: 1.167 minutes, Δt: 7.486 minutes, CFL: 7.50e-01
i: 940, sim time: 8.344 days, wall time: 1.186 minutes, Δt: 7.766 minutes, CFL: 7.50e-01
i: 960, sim time: 8.450 days, wall time: 1.202 minutes, Δt: 8.009 minutes, CFL: 7.50e-01
i: 980, sim time: 8.565 days, wall time: 1.220 minutes, Δt: 8.932 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.686 days, wall time: 1.238 minutes, Δt: 9.549 minutes, CFL: 7.50e-01
i: 1020, sim time: 8.817 days, wall time: 1.254 minutes, Δt: 9.640 minutes, CFL: 7.50e-01
i: 1040, sim time: 8.944 days, wall time: 1.272 minutes, Δt: 9.679 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.071 days, wall time: 1.289 minutes, Δt: 9.128 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.190 days, wall time: 1.307 minutes, Δt: 8.391 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.302 days, wall time: 1.324 minutes, Δt: 8.244 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.415 days, wall time: 1.341 minutes, Δt: 8.567 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.530 days, wall time: 1.359 minutes, Δt: 8.817 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.651 days, wall time: 1.376 minutes, Δt: 8.921 minutes, CFL: 7.50e-01
i: 1180, sim time: 9.769 days, wall time: 1.394 minutes, Δt: 9.070 minutes, CFL: 7.50e-01
i: 1200, sim time: 9.890 days, wall time: 1.411 minutes, Δt: 9.073 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.427 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.