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(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
├── 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 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 (31.225 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (9.059 seconds).
i: 20, sim time: 6 hours, wall time: 0 seconds, Δt: 25.164 minutes, CFL: 4.74e-01
i: 40, sim time: 15 hours, wall time: 0 seconds, Δt: 30 minutes, CFL: 5.77e-01
i: 60, sim time: 1.042 days, wall time: 0 seconds, Δt: 30 minutes, CFL: 6.12e-01
i: 80, sim time: 1.458 days, wall time: 0 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 0 seconds, Δt: 30 minutes, CFL: 6.69e-01
i: 120, sim time: 2.292 days, wall time: 0 seconds, Δt: 29.577 minutes, CFL: 7.50e-01
i: 140, sim time: 2.621 days, wall time: 0 seconds, Δt: 26.294 minutes, CFL: 7.50e-01
i: 160, sim time: 2.933 days, wall time: 0 seconds, Δt: 23.119 minutes, CFL: 7.50e-01
i: 180, sim time: 3.211 days, wall time: 0 seconds, Δt: 20.878 minutes, CFL: 7.50e-01
i: 200, sim time: 3.487 days, wall time: 0 seconds, Δt: 19.433 minutes, CFL: 7.50e-01
i: 220, sim time: 3.729 days, wall time: 0 seconds, Δt: 17.105 minutes, CFL: 7.50e-01
i: 240, sim time: 3.939 days, wall time: 0 seconds, Δt: 15.493 minutes, CFL: 7.50e-01
i: 260, sim time: 4.131 days, wall time: 0 seconds, Δt: 13.540 minutes, CFL: 7.50e-01
i: 280, sim time: 4.314 days, wall time: 0 seconds, Δt: 13.106 minutes, CFL: 7.50e-01
i: 300, sim time: 4.478 days, wall time: 0 seconds, Δt: 12.763 minutes, CFL: 7.50e-01
i: 320, sim time: 4.649 days, wall time: 0 seconds, Δt: 13.968 minutes, CFL: 7.50e-01
i: 340, sim time: 4.832 days, wall time: 0 seconds, Δt: 12.660 minutes, CFL: 7.50e-01
i: 360, sim time: 4.991 days, wall time: 0 seconds, Δt: 12.120 minutes, CFL: 7.50e-01
i: 380, sim time: 5.138 days, wall time: 0 seconds, Δt: 10.965 minutes, CFL: 7.50e-01
i: 400, sim time: 5.288 days, wall time: 0 seconds, Δt: 11.202 minutes, CFL: 7.50e-01
i: 420, sim time: 5.440 days, wall time: 0 seconds, Δt: 11.052 minutes, CFL: 7.50e-01
i: 440, sim time: 5.583 days, wall time: 0 seconds, Δt: 10.787 minutes, CFL: 7.50e-01
i: 460, sim time: 5.726 days, wall time: 0 seconds, Δt: 10.563 minutes, CFL: 7.50e-01
i: 480, sim time: 5.862 days, wall time: 0 seconds, Δt: 10.516 minutes, CFL: 7.50e-01
i: 500, sim time: 6 days, wall time: 0 seconds, Δt: 10.328 minutes, CFL: 7.50e-01
i: 520, sim time: 6.140 days, wall time: 0 seconds, Δt: 10.214 minutes, CFL: 7.50e-01
i: 540, sim time: 6.278 days, wall time: 0 seconds, Δt: 10.213 minutes, CFL: 7.50e-01
i: 560, sim time: 6.417 days, wall time: 0 seconds, Δt: 10.509 minutes, CFL: 7.50e-01
i: 580, sim time: 6.553 days, wall time: 0 seconds, Δt: 10.780 minutes, CFL: 7.50e-01
i: 600, sim time: 6.688 days, wall time: 0 seconds, Δt: 10.419 minutes, CFL: 7.50e-01
i: 620, sim time: 6.830 days, wall time: 0 seconds, Δt: 10.627 minutes, CFL: 7.50e-01
i: 640, sim time: 6.968 days, wall time: 0 seconds, Δt: 10.646 minutes, CFL: 7.50e-01
i: 660, sim time: 7.105 days, wall time: 0 seconds, Δt: 10.692 minutes, CFL: 7.50e-01
i: 680, sim time: 7.242 days, wall time: 0 seconds, Δt: 10.766 minutes, CFL: 7.50e-01
i: 700, sim time: 7.378 days, wall time: 0 seconds, Δt: 10.861 minutes, CFL: 7.50e-01
i: 720, sim time: 7.523 days, wall time: 0 seconds, Δt: 10.851 minutes, CFL: 7.50e-01
i: 740, sim time: 7.653 days, wall time: 0 seconds, Δt: 9.594 minutes, CFL: 7.50e-01
i: 760, sim time: 7.781 days, wall time: 0 seconds, Δt: 8.594 minutes, CFL: 7.50e-01
i: 780, sim time: 7.900 days, wall time: 0 seconds, Δt: 8.651 minutes, CFL: 7.50e-01
i: 800, sim time: 8.018 days, wall time: 0 seconds, Δt: 8.673 minutes, CFL: 7.50e-01
i: 820, sim time: 8.140 days, wall time: 0 seconds, Δt: 9.218 minutes, CFL: 7.50e-01
i: 840, sim time: 8.264 days, wall time: 0 seconds, Δt: 10.123 minutes, CFL: 7.50e-01
i: 860, sim time: 8.409 days, wall time: 0 seconds, Δt: 11.070 minutes, CFL: 7.50e-01
i: 880, sim time: 8.560 days, wall time: 0 seconds, Δt: 10.373 minutes, CFL: 7.50e-01
i: 900, sim time: 8.687 days, wall time: 0 seconds, Δt: 9.740 minutes, CFL: 7.50e-01
i: 920, sim time: 8.819 days, wall time: 0 seconds, Δt: 10.165 minutes, CFL: 7.50e-01
i: 940, sim time: 8.944 days, wall time: 0 seconds, Δt: 9.777 minutes, CFL: 7.50e-01
i: 960, sim time: 9.073 days, wall time: 0 seconds, Δt: 9.635 minutes, CFL: 7.50e-01
i: 980, sim time: 9.200 days, wall time: 0 seconds, Δt: 9.783 minutes, CFL: 7.50e-01
i: 1000, sim time: 9.333 days, wall time: 0 seconds, Δt: 10.662 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.475 days, wall time: 0 seconds, Δt: 10.249 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.611 days, wall time: 0 seconds, Δt: 10.047 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.750 days, wall time: 0 seconds, Δt: 10.026 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.875 days, wall time: 0 seconds, Δt: 9.852 minutes, CFL: 7.50e-01
i: 1100, sim time: 10.000 days, wall time: 0 seconds, Δt: 9.835 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 0 seconds.
[ 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.