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.
grid = RectilinearGrid(size = (32, 32, 8), extent = (1kilometer, 1kilometer, 100meters))32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU 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"##227".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,
carbonates = true,
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{Float64, OceanBioME.Models.GasExchangeModel.ScaledGasTransferVelocity.SchmidtScaledTransferVelocity{OceanBioME.Models.GasExchangeModel.PolynomialParameterisation{2, Tuple{Float64, Float64, Float64}}, OceanBioME.Models.GasExchangeModel.PolynomialParameterisation{4, NTuple{5, Float64}}}, OceanBioME.Models.GasExchangeModel.CarbonDioxideConcentration{CarbonChemistry{OceanBioME.Models.CarbonChemistryModel.K0{Float64}, @NamedTuple{K1::OceanBioME.Models.CarbonChemistryModel.K1{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, K2::OceanBioME.Models.CarbonChemistryModel.K2{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}}, OceanBioME.Models.CarbonChemistryModel.KB{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.KS{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.KF{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, OceanBioME.Models.CarbonChemistryModel.KS{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, @NamedTuple{KP1::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, KP2::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, KP3::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}}, OceanBioME.Models.CarbonChemistryModel.KSi{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64}, OceanBioME.Models.CarbonChemistryModel.KW{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, OceanBioME.Models.CarbonChemistryModel.KSP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, typeof(OceanBioME.Models.teos10_polynomial_approximation)}, OceanBioME.Models.GasExchangeModel.PolynomialVirialCoefficientForCarbonDioxide{Float64}, OceanBioME.Models.GasExchangeModel.CrossVirialCoefficientForCarbonDioxide{Float64}, Int64, Nothing}, Float64})
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Model instantiation
model = NonhydrostaticModel(; grid,
biogeochemistry,
boundary_conditions = (DIC = DIC_bcs, ),
advection = WENO(grid),
timestepper = :RungeKutta3,
coriolis,
tracers = (:T, :S),
buoyancy,
background_fields = (T = T_field, v = V_field),
closure = vertical_diffusivity)NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: WENO reconstruction 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{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 15.625 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 10 days
├── Stop iteration : Inf
├── Wall time limit: Inf
├── 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 entries
└── Diagnostics: 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 Oceananigans.Architectures.CPU with 3×3×3 halo:
├── u: 32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on RectilinearGrid on Oceananigans.Architectures.CPU
├── v: 32×32×8 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on RectilinearGrid on Oceananigans.Architectures.CPU
└── w: 32×32×9 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on RectilinearGrid on Oceananigans.Architectures.CPUand also calculate the vertical vorticity.
ζ = Field(∂x(v) - ∂y(u))32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on RectilinearGrid on Oceananigans.Architectures.CPU
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 38×38×14 OffsetArray(::Array{Float64, 3}, -2:35, -2:35, -2:11) with eltype Float64 with indices -2:35×-2:35×-2:11
└── max=0.0, min=0.0, mean=0.0Periodically save the velocities and vorticity to a file.
simulation.output_writers[:fields] = JLD2OutputWriter(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 (33.215 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (9.389 seconds).
i: 20, sim time: 6 hours, wall time: 43.929 seconds, Δt: 25.164 minutes, CFL: 4.74e-01
i: 40, sim time: 15 hours, wall time: 45.368 seconds, Δt: 30 minutes, CFL: 5.77e-01
i: 60, sim time: 1.042 days, wall time: 46.716 seconds, Δt: 30 minutes, CFL: 6.12e-01
i: 80, sim time: 1.458 days, wall time: 48.113 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 49.531 seconds, Δt: 30 minutes, CFL: 6.69e-01
i: 120, sim time: 2.292 days, wall time: 50.855 seconds, Δt: 29.577 minutes, CFL: 7.50e-01
i: 140, sim time: 2.602 days, wall time: 52.175 seconds, Δt: 26.318 minutes, CFL: 7.50e-01
i: 160, sim time: 2.917 days, wall time: 53.496 seconds, Δt: 23.040 minutes, CFL: 7.50e-01
i: 180, sim time: 3.196 days, wall time: 54.855 seconds, Δt: 20.903 minutes, CFL: 7.50e-01
i: 200, sim time: 3.458 days, wall time: 56.211 seconds, Δt: 19.496 minutes, CFL: 7.50e-01
i: 220, sim time: 3.692 days, wall time: 57.536 seconds, Δt: 17.453 minutes, CFL: 7.50e-01
i: 240, sim time: 3.912 days, wall time: 58.836 seconds, Δt: 15.803 minutes, CFL: 7.50e-01
i: 260, sim time: 4.103 days, wall time: 1.002 minutes, Δt: 13.631 minutes, CFL: 7.50e-01
i: 280, sim time: 4.277 days, wall time: 1.024 minutes, Δt: 12.951 minutes, CFL: 7.50e-01
i: 300, sim time: 4.444 days, wall time: 1.046 minutes, Δt: 12.758 minutes, CFL: 7.50e-01
i: 320, sim time: 4.611 days, wall time: 1.068 minutes, Δt: 13.512 minutes, CFL: 7.50e-01
i: 340, sim time: 4.796 days, wall time: 1.090 minutes, Δt: 13.000 minutes, CFL: 7.50e-01
i: 360, sim time: 4.958 days, wall time: 1.112 minutes, Δt: 11.773 minutes, CFL: 7.50e-01
i: 380, sim time: 5.107 days, wall time: 1.134 minutes, Δt: 11.086 minutes, CFL: 7.50e-01
i: 400, sim time: 5.258 days, wall time: 1.156 minutes, Δt: 11.009 minutes, CFL: 7.50e-01
i: 420, sim time: 5.411 days, wall time: 1.178 minutes, Δt: 11.241 minutes, CFL: 7.50e-01
i: 440, sim time: 5.561 days, wall time: 1.200 minutes, Δt: 10.797 minutes, CFL: 7.50e-01
i: 460, sim time: 5.696 days, wall time: 1.223 minutes, Δt: 10.701 minutes, CFL: 7.50e-01
i: 480, sim time: 5.833 days, wall time: 1.246 minutes, Δt: 10.466 minutes, CFL: 7.50e-01
i: 500, sim time: 5.974 days, wall time: 1.270 minutes, Δt: 10.336 minutes, CFL: 7.50e-01
i: 520, sim time: 6.112 days, wall time: 1.293 minutes, Δt: 10.249 minutes, CFL: 7.50e-01
i: 540, sim time: 6.250 days, wall time: 1.316 minutes, Δt: 10.178 minutes, CFL: 7.50e-01
i: 560, sim time: 6.391 days, wall time: 1.339 minutes, Δt: 10.543 minutes, CFL: 7.50e-01
i: 580, sim time: 6.522 days, wall time: 1.362 minutes, Δt: 10.866 minutes, CFL: 7.50e-01
i: 600, sim time: 6.664 days, wall time: 1.384 minutes, Δt: 10.386 minutes, CFL: 7.50e-01
i: 620, sim time: 6.801 days, wall time: 1.407 minutes, Δt: 10.549 minutes, CFL: 7.50e-01
i: 640, sim time: 6.931 days, wall time: 1.428 minutes, Δt: 10.627 minutes, CFL: 7.50e-01
i: 660, sim time: 7.074 days, wall time: 1.450 minutes, Δt: 10.613 minutes, CFL: 7.50e-01
i: 680, sim time: 7.211 days, wall time: 1.472 minutes, Δt: 10.728 minutes, CFL: 7.50e-01
i: 700, sim time: 7.341 days, wall time: 1.493 minutes, Δt: 10.785 minutes, CFL: 7.50e-01
i: 720, sim time: 7.485 days, wall time: 1.515 minutes, Δt: 10.991 minutes, CFL: 7.50e-01
i: 740, sim time: 7.626 days, wall time: 1.537 minutes, Δt: 9.820 minutes, CFL: 7.50e-01
i: 760, sim time: 7.756 days, wall time: 1.559 minutes, Δt: 8.774 minutes, CFL: 7.50e-01
i: 780, sim time: 7.876 days, wall time: 1.582 minutes, Δt: 8.622 minutes, CFL: 7.50e-01
i: 800, sim time: 7.995 days, wall time: 1.605 minutes, Δt: 8.735 minutes, CFL: 7.50e-01
i: 820, sim time: 8.115 days, wall time: 1.628 minutes, Δt: 9.073 minutes, CFL: 7.50e-01
i: 840, sim time: 8.240 days, wall time: 1.650 minutes, Δt: 9.916 minutes, CFL: 7.50e-01
i: 860, sim time: 8.379 days, wall time: 1.672 minutes, Δt: 10.964 minutes, CFL: 7.50e-01
i: 880, sim time: 8.531 days, wall time: 1.694 minutes, Δt: 10.674 minutes, CFL: 7.50e-01
i: 900, sim time: 8.667 days, wall time: 1.716 minutes, Δt: 9.912 minutes, CFL: 7.50e-01
i: 920, sim time: 8.798 days, wall time: 1.738 minutes, Δt: 10.138 minutes, CFL: 7.50e-01
i: 940, sim time: 8.931 days, wall time: 1.759 minutes, Δt: 9.865 minutes, CFL: 7.50e-01
i: 960, sim time: 9.060 days, wall time: 1.781 minutes, Δt: 9.584 minutes, CFL: 7.50e-01
i: 980, sim time: 9.187 days, wall time: 1.803 minutes, Δt: 9.731 minutes, CFL: 7.50e-01
i: 1000, sim time: 9.321 days, wall time: 1.824 minutes, Δt: 10.590 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.460 days, wall time: 1.846 minutes, Δt: 10.296 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.597 days, wall time: 1.867 minutes, Δt: 10.072 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.736 days, wall time: 1.889 minutes, Δt: 9.978 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.861 days, wall time: 1.910 minutes, Δt: 9.877 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.992 days, wall time: 1.932 minutes, Δt: 9.819 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.934 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.