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}, Float64, Nothing}, Float64})
└── 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{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{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{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 => 4
│ ├── stop_iteration_exceeded => -
│ ├── wall_time_limit_exceeded => e
│ └── nan_checker => }
├── 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=5.35045e-5, min=-4.94356e-5, mean=-2.33679e-23Periodically 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.27e-01
[ Info: ... simulation initialization complete (25.949 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (12.404 seconds).
i: 20, sim time: 6 hours, wall time: 39.436 seconds, Δt: 25.164 minutes, CFL: 4.74e-01
i: 40, sim time: 15 hours, wall time: 40.259 seconds, Δt: 30 minutes, CFL: 5.77e-01
i: 60, sim time: 1.042 days, wall time: 41.087 seconds, Δt: 30 minutes, CFL: 6.12e-01
i: 80, sim time: 1.458 days, wall time: 41.915 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 42.730 seconds, Δt: 30 minutes, CFL: 6.69e-01
i: 120, sim time: 2.292 days, wall time: 43.544 seconds, Δt: 29.577 minutes, CFL: 7.50e-01
i: 140, sim time: 2.621 days, wall time: 44.356 seconds, Δt: 26.294 minutes, CFL: 7.50e-01
i: 160, sim time: 2.933 days, wall time: 45.162 seconds, Δt: 23.119 minutes, CFL: 7.50e-01
i: 180, sim time: 3.196 days, wall time: 45.958 seconds, Δt: 20.903 minutes, CFL: 7.50e-01
i: 200, sim time: 3.459 days, wall time: 46.763 seconds, Δt: 19.497 minutes, CFL: 7.50e-01
i: 220, sim time: 3.692 days, wall time: 47.588 seconds, Δt: 17.453 minutes, CFL: 7.50e-01
i: 240, sim time: 3.912 days, wall time: 48.393 seconds, Δt: 15.802 minutes, CFL: 7.50e-01
i: 260, sim time: 4.103 days, wall time: 49.203 seconds, Δt: 13.633 minutes, CFL: 7.50e-01
i: 280, sim time: 4.277 days, wall time: 50.011 seconds, Δt: 12.955 minutes, CFL: 7.50e-01
i: 300, sim time: 4.444 days, wall time: 50.810 seconds, Δt: 12.758 minutes, CFL: 7.50e-01
i: 320, sim time: 4.611 days, wall time: 51.608 seconds, Δt: 13.516 minutes, CFL: 7.50e-01
i: 340, sim time: 4.796 days, wall time: 52.410 seconds, Δt: 13.001 minutes, CFL: 7.50e-01
i: 360, sim time: 4.958 days, wall time: 53.209 seconds, Δt: 11.774 minutes, CFL: 7.50e-01
i: 380, sim time: 5.107 days, wall time: 54.009 seconds, Δt: 11.085 minutes, CFL: 7.50e-01
i: 400, sim time: 5.258 days, wall time: 54.821 seconds, Δt: 11.009 minutes, CFL: 7.50e-01
i: 420, sim time: 5.411 days, wall time: 55.614 seconds, Δt: 11.242 minutes, CFL: 7.50e-01
i: 440, sim time: 5.553 days, wall time: 56.404 seconds, Δt: 10.848 minutes, CFL: 7.50e-01
i: 460, sim time: 5.689 days, wall time: 57.211 seconds, Δt: 10.594 minutes, CFL: 7.50e-01
i: 480, sim time: 5.831 days, wall time: 58.014 seconds, Δt: 10.474 minutes, CFL: 7.50e-01
i: 500, sim time: 5.967 days, wall time: 58.819 seconds, Δt: 10.327 minutes, CFL: 7.50e-01
i: 520, sim time: 6.105 days, wall time: 59.627 seconds, Δt: 10.246 minutes, CFL: 7.50e-01
i: 540, sim time: 6.244 days, wall time: 1.007 minutes, Δt: 10.159 minutes, CFL: 7.50e-01
i: 560, sim time: 6.384 days, wall time: 1.020 minutes, Δt: 10.530 minutes, CFL: 7.50e-01
i: 580, sim time: 6.522 days, wall time: 1.034 minutes, Δt: 10.866 minutes, CFL: 7.50e-01
i: 600, sim time: 6.664 days, wall time: 1.047 minutes, Δt: 10.383 minutes, CFL: 7.50e-01
i: 620, sim time: 6.794 days, wall time: 1.061 minutes, Δt: 10.531 minutes, CFL: 7.50e-01
i: 640, sim time: 6.924 days, wall time: 1.074 minutes, Δt: 10.623 minutes, CFL: 7.50e-01
i: 660, sim time: 7.067 days, wall time: 1.087 minutes, Δt: 10.612 minutes, CFL: 7.50e-01
i: 680, sim time: 7.204 days, wall time: 1.101 minutes, Δt: 10.777 minutes, CFL: 7.50e-01
i: 700, sim time: 7.341 days, wall time: 1.114 minutes, Δt: 10.789 minutes, CFL: 7.50e-01
i: 720, sim time: 7.485 days, wall time: 1.127 minutes, Δt: 11.004 minutes, CFL: 7.50e-01
i: 740, sim time: 7.626 days, wall time: 1.140 minutes, Δt: 9.836 minutes, CFL: 7.50e-01
i: 760, sim time: 7.756 days, wall time: 1.154 minutes, Δt: 8.779 minutes, CFL: 7.50e-01
i: 780, sim time: 7.876 days, wall time: 1.167 minutes, Δt: 8.625 minutes, CFL: 7.50e-01
i: 800, sim time: 7.995 days, wall time: 1.180 minutes, Δt: 8.748 minutes, CFL: 7.50e-01
i: 820, sim time: 8.115 days, wall time: 1.194 minutes, Δt: 9.076 minutes, CFL: 7.50e-01
i: 840, sim time: 8.240 days, wall time: 1.207 minutes, Δt: 9.916 minutes, CFL: 7.50e-01
i: 860, sim time: 8.379 days, wall time: 1.220 minutes, Δt: 10.973 minutes, CFL: 7.50e-01
i: 880, sim time: 8.531 days, wall time: 1.234 minutes, Δt: 10.690 minutes, CFL: 7.50e-01
i: 900, sim time: 8.667 days, wall time: 1.247 minutes, Δt: 9.917 minutes, CFL: 7.50e-01
i: 920, sim time: 8.798 days, wall time: 1.261 minutes, Δt: 10.134 minutes, CFL: 7.50e-01
i: 940, sim time: 8.931 days, wall time: 1.274 minutes, Δt: 9.879 minutes, CFL: 7.50e-01
i: 960, sim time: 9.060 days, wall time: 1.287 minutes, Δt: 9.580 minutes, CFL: 7.50e-01
i: 980, sim time: 9.187 days, wall time: 1.301 minutes, Δt: 9.717 minutes, CFL: 7.50e-01
i: 1000, sim time: 9.321 days, wall time: 1.314 minutes, Δt: 10.580 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.460 days, wall time: 1.328 minutes, Δt: 10.306 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.597 days, wall time: 1.341 minutes, Δt: 10.076 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.736 days, wall time: 1.354 minutes, Δt: 9.980 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.861 days, wall time: 1.367 minutes, Δt: 9.886 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.992 days, wall time: 1.381 minutes, Δt: 9.821 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.382 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.