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.034 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.707 seconds).
i: 20, sim time: 6 hours, wall time: 22.919 seconds, Δt: 25.164 minutes, CFL: 4.70e-01
i: 40, sim time: 15 hours, wall time: 24.061 seconds, Δt: 30 minutes, CFL: 5.88e-01
i: 60, sim time: 1.042 days, wall time: 25.114 seconds, Δt: 30 minutes, CFL: 6.06e-01
i: 80, sim time: 1.458 days, wall time: 26.170 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 27.221 seconds, Δt: 30 minutes, CFL: 6.52e-01
i: 120, sim time: 2.292 days, wall time: 28.274 seconds, Δt: 30 minutes, CFL: 7.16e-01
i: 140, sim time: 2.687 days, wall time: 29.370 seconds, Δt: 28.626 minutes, CFL: 7.50e-01
i: 160, sim time: 3.018 days, wall time: 30.427 seconds, Δt: 24.961 minutes, CFL: 7.50e-01
i: 180, sim time: 3.323 days, wall time: 31.336 seconds, Δt: 19.356 minutes, CFL: 7.50e-01
i: 200, sim time: 3.559 days, wall time: 32.384 seconds, Δt: 16.277 minutes, CFL: 7.50e-01
i: 220, sim time: 3.750 days, wall time: 33.427 seconds, Δt: 14.399 minutes, CFL: 7.50e-01
i: 240, sim time: 3.936 days, wall time: 34.571 seconds, Δt: 13.104 minutes, CFL: 7.50e-01
i: 260, sim time: 4.101 days, wall time: 35.621 seconds, Δt: 13.037 minutes, CFL: 7.50e-01
i: 280, sim time: 4.268 days, wall time: 36.666 seconds, Δt: 12.875 minutes, CFL: 7.50e-01
i: 300, sim time: 4.434 days, wall time: 37.716 seconds, Δt: 12.278 minutes, CFL: 7.50e-01
i: 320, sim time: 4.600 days, wall time: 38.765 seconds, Δt: 12.097 minutes, CFL: 7.50e-01
i: 340, sim time: 4.767 days, wall time: 39.808 seconds, Δt: 11.919 minutes, CFL: 7.50e-01
i: 360, sim time: 4.933 days, wall time: 40.856 seconds, Δt: 11.734 minutes, CFL: 7.50e-01
i: 380, sim time: 5.083 days, wall time: 41.800 seconds, Δt: 11.608 minutes, CFL: 7.50e-01
i: 400, sim time: 5.241 days, wall time: 42.849 seconds, Δt: 12.021 minutes, CFL: 7.50e-01
i: 420, sim time: 5.404 days, wall time: 43.896 seconds, Δt: 13.007 minutes, CFL: 7.50e-01
i: 440, sim time: 5.569 days, wall time: 44.940 seconds, Δt: 12.059 minutes, CFL: 7.50e-01
i: 460, sim time: 5.722 days, wall time: 45.983 seconds, Δt: 11.125 minutes, CFL: 7.50e-01
i: 480, sim time: 5.863 days, wall time: 47.030 seconds, Δt: 10.535 minutes, CFL: 7.50e-01
i: 500, sim time: 6.000 days, wall time: 48.062 seconds, Δt: 9.654 minutes, CFL: 7.50e-01
i: 520, sim time: 6.120 days, wall time: 49.102 seconds, Δt: 8.777 minutes, CFL: 7.50e-01
i: 540, sim time: 6.237 days, wall time: 50.142 seconds, Δt: 8.368 minutes, CFL: 7.50e-01
i: 560, sim time: 6.345 days, wall time: 51.279 seconds, Δt: 8.016 minutes, CFL: 7.50e-01
i: 580, sim time: 6.449 days, wall time: 52.220 seconds, Δt: 7.549 minutes, CFL: 7.50e-01
i: 600, sim time: 6.550 days, wall time: 53.255 seconds, Δt: 6.965 minutes, CFL: 7.50e-01
i: 620, sim time: 6.644 days, wall time: 54.287 seconds, Δt: 6.712 minutes, CFL: 7.50e-01
i: 640, sim time: 6.735 days, wall time: 55.321 seconds, Δt: 6.509 minutes, CFL: 7.50e-01
i: 660, sim time: 6.822 days, wall time: 56.353 seconds, Δt: 6.519 minutes, CFL: 7.50e-01
i: 680, sim time: 6.913 days, wall time: 57.389 seconds, Δt: 7.131 minutes, CFL: 7.50e-01
i: 700, sim time: 7.010 days, wall time: 58.527 seconds, Δt: 7.463 minutes, CFL: 7.50e-01
i: 720, sim time: 7.115 days, wall time: 59.483 seconds, Δt: 7.581 minutes, CFL: 7.50e-01
i: 740, sim time: 7.221 days, wall time: 1.009 minutes, Δt: 7.896 minutes, CFL: 7.50e-01
i: 760, sim time: 7.329 days, wall time: 1.026 minutes, Δt: 8.367 minutes, CFL: 7.50e-01
i: 780, sim time: 7.440 days, wall time: 1.044 minutes, Δt: 8.879 minutes, CFL: 7.50e-01
i: 800, sim time: 7.563 days, wall time: 1.061 minutes, Δt: 9.111 minutes, CFL: 7.50e-01
i: 820, sim time: 7.679 days, wall time: 1.080 minutes, Δt: 8.773 minutes, CFL: 7.50e-01
i: 840, sim time: 7.797 days, wall time: 1.096 minutes, Δt: 8.530 minutes, CFL: 7.50e-01
i: 860, sim time: 7.913 days, wall time: 1.113 minutes, Δt: 9.154 minutes, CFL: 7.50e-01
i: 880, sim time: 8.030 days, wall time: 1.131 minutes, Δt: 8.279 minutes, CFL: 7.50e-01
i: 900, sim time: 8.138 days, wall time: 1.148 minutes, Δt: 7.768 minutes, CFL: 7.50e-01
i: 920, sim time: 8.240 days, wall time: 1.165 minutes, Δt: 7.486 minutes, CFL: 7.50e-01
i: 940, sim time: 8.344 days, wall time: 1.184 minutes, Δt: 7.766 minutes, CFL: 7.50e-01
i: 960, sim time: 8.450 days, wall time: 1.200 minutes, Δt: 8.009 minutes, CFL: 7.50e-01
i: 980, sim time: 8.565 days, wall time: 1.217 minutes, Δt: 8.932 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.686 days, wall time: 1.236 minutes, Δt: 9.549 minutes, CFL: 7.50e-01
i: 1020, sim time: 8.817 days, wall time: 1.252 minutes, Δt: 9.640 minutes, CFL: 7.50e-01
i: 1040, sim time: 8.944 days, wall time: 1.270 minutes, Δt: 9.679 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.071 days, wall time: 1.287 minutes, Δt: 9.128 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.190 days, wall time: 1.304 minutes, Δt: 8.391 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.302 days, wall time: 1.322 minutes, Δt: 8.244 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.415 days, wall time: 1.339 minutes, Δt: 8.567 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.530 days, wall time: 1.356 minutes, Δt: 8.817 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.651 days, wall time: 1.374 minutes, Δt: 8.921 minutes, CFL: 7.50e-01
i: 1180, sim time: 9.769 days, wall time: 1.392 minutes, Δt: 9.070 minutes, CFL: 7.50e-01
i: 1200, sim time: 9.890 days, wall time: 1.408 minutes, Δt: 9.073 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.425 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.