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{Int64, OceanBioME.Models.GasExchangeModel.ScaledGasTransferVelocity.SchmidtScaledTransferVelocity{OceanBioME.Models.GasExchangeModel.PolynomialParameterisation{2, Tuple{Int64, Int64, 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}, Int64})
└── 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 (3.930 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (8.519 seconds).
i: 20, sim time: 6 hours, wall time: 13.665 seconds, Δt: 25.164 minutes, CFL: 4.74e-01
i: 40, sim time: 15 hours, wall time: 14.988 seconds, Δt: 30 minutes, CFL: 5.77e-01
i: 60, sim time: 1.042 days, wall time: 16.307 seconds, Δt: 30 minutes, CFL: 6.12e-01
i: 80, sim time: 1.458 days, wall time: 17.623 seconds, Δt: 30 minutes, CFL: 6.33e-01
i: 100, sim time: 1.875 days, wall time: 18.925 seconds, Δt: 30 minutes, CFL: 6.69e-01
i: 120, sim time: 2.292 days, wall time: 20.225 seconds, Δt: 29.577 minutes, CFL: 7.50e-01
i: 140, sim time: 2.621 days, wall time: 21.521 seconds, Δt: 26.295 minutes, CFL: 7.50e-01
i: 160, sim time: 2.933 days, wall time: 22.812 seconds, Δt: 23.121 minutes, CFL: 7.50e-01
i: 180, sim time: 3.196 days, wall time: 24.091 seconds, Δt: 20.903 minutes, CFL: 7.50e-01
i: 200, sim time: 3.472 days, wall time: 25.373 seconds, Δt: 19.466 minutes, CFL: 7.50e-01
i: 220, sim time: 3.692 days, wall time: 26.658 seconds, Δt: 17.454 minutes, CFL: 7.50e-01
i: 240, sim time: 3.912 days, wall time: 27.937 seconds, Δt: 15.801 minutes, CFL: 7.50e-01
i: 260, sim time: 4.093 days, wall time: 29.215 seconds, Δt: 13.710 minutes, CFL: 7.50e-01
i: 280, sim time: 4.268 days, wall time: 30.491 seconds, Δt: 13.195 minutes, CFL: 7.50e-01
i: 300, sim time: 4.434 days, wall time: 31.776 seconds, Δt: 12.676 minutes, CFL: 7.50e-01
i: 320, sim time: 4.602 days, wall time: 33.060 seconds, Δt: 13.634 minutes, CFL: 7.50e-01
i: 340, sim time: 4.788 days, wall time: 34.348 seconds, Δt: 12.943 minutes, CFL: 7.50e-01
i: 360, sim time: 4.950 days, wall time: 35.633 seconds, Δt: 11.716 minutes, CFL: 7.50e-01
i: 380, sim time: 5.099 days, wall time: 36.920 seconds, Δt: 11.171 minutes, CFL: 7.50e-01
i: 400, sim time: 5.250 days, wall time: 38.201 seconds, Δt: 10.998 minutes, CFL: 7.50e-01
i: 420, sim time: 5.403 days, wall time: 39.480 seconds, Δt: 11.142 minutes, CFL: 7.50e-01
i: 440, sim time: 5.546 days, wall time: 40.755 seconds, Δt: 10.918 minutes, CFL: 7.50e-01
i: 460, sim time: 5.689 days, wall time: 42.033 seconds, Δt: 10.593 minutes, CFL: 7.50e-01
i: 480, sim time: 5.831 days, wall time: 43.310 seconds, Δt: 10.472 minutes, CFL: 7.50e-01
i: 500, sim time: 5.967 days, wall time: 44.592 seconds, Δt: 10.324 minutes, CFL: 7.50e-01
i: 520, sim time: 6.098 days, wall time: 45.874 seconds, Δt: 10.261 minutes, CFL: 7.50e-01
i: 540, sim time: 6.237 days, wall time: 47.154 seconds, Δt: 10.151 minutes, CFL: 7.50e-01
i: 560, sim time: 6.377 days, wall time: 48.440 seconds, Δt: 10.515 minutes, CFL: 7.50e-01
i: 580, sim time: 6.515 days, wall time: 49.718 seconds, Δt: 10.836 minutes, CFL: 7.50e-01
i: 600, sim time: 6.657 days, wall time: 50.995 seconds, Δt: 10.386 minutes, CFL: 7.50e-01
i: 620, sim time: 6.794 days, wall time: 52.273 seconds, Δt: 10.534 minutes, CFL: 7.50e-01
i: 640, sim time: 6.931 days, wall time: 53.542 seconds, Δt: 10.626 minutes, CFL: 7.50e-01
i: 660, sim time: 7.074 days, wall time: 54.813 seconds, Δt: 10.611 minutes, CFL: 7.50e-01
i: 680, sim time: 7.211 days, wall time: 56.085 seconds, Δt: 10.726 minutes, CFL: 7.50e-01
i: 700, sim time: 7.348 days, wall time: 57.365 seconds, Δt: 10.816 minutes, CFL: 7.50e-01
i: 720, sim time: 7.492 days, wall time: 58.638 seconds, Δt: 10.977 minutes, CFL: 7.50e-01
i: 740, sim time: 7.633 days, wall time: 59.914 seconds, Δt: 9.755 minutes, CFL: 7.50e-01
i: 760, sim time: 7.756 days, wall time: 1.020 minutes, Δt: 8.771 minutes, CFL: 7.50e-01
i: 780, sim time: 7.876 days, wall time: 1.041 minutes, Δt: 8.620 minutes, CFL: 7.50e-01
i: 800, sim time: 7.995 days, wall time: 1.063 minutes, Δt: 8.727 minutes, CFL: 7.50e-01
i: 820, sim time: 8.115 days, wall time: 1.084 minutes, Δt: 9.071 minutes, CFL: 7.50e-01
i: 840, sim time: 8.240 days, wall time: 1.105 minutes, Δt: 9.919 minutes, CFL: 7.50e-01
i: 860, sim time: 8.379 days, wall time: 1.127 minutes, Δt: 10.960 minutes, CFL: 7.50e-01
i: 880, sim time: 8.531 days, wall time: 1.148 minutes, Δt: 10.663 minutes, CFL: 7.50e-01
i: 900, sim time: 8.667 days, wall time: 1.170 minutes, Δt: 9.910 minutes, CFL: 7.50e-01
i: 920, sim time: 8.798 days, wall time: 1.191 minutes, Δt: 10.140 minutes, CFL: 7.50e-01
i: 940, sim time: 8.931 days, wall time: 1.212 minutes, Δt: 9.854 minutes, CFL: 7.50e-01
i: 960, sim time: 9.060 days, wall time: 1.234 minutes, Δt: 9.585 minutes, CFL: 7.50e-01
i: 980, sim time: 9.187 days, wall time: 1.255 minutes, Δt: 9.740 minutes, CFL: 7.50e-01
i: 1000, sim time: 9.321 days, wall time: 1.276 minutes, Δt: 10.596 minutes, CFL: 7.50e-01
i: 1020, sim time: 9.460 days, wall time: 1.298 minutes, Δt: 10.290 minutes, CFL: 7.50e-01
i: 1040, sim time: 9.597 days, wall time: 1.319 minutes, Δt: 10.070 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.736 days, wall time: 1.341 minutes, Δt: 9.977 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.861 days, wall time: 1.362 minutes, Δt: 9.870 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.992 days, wall time: 1.383 minutes, Δt: 9.819 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 1.385 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.