One dimensional column forced by external data with carbonate chemistry

In this example we setup a simple 1D column with the LOBSTER biogeochemical model and observe its evolution. This example demonstrates:

  • How to setup OceanBioME's biogeochemical models
  • How to load external forcing data
  • How to run with optional tracer sets such as carbonate chemistry
  • How to setup a non-uniform grid for better near surface resolution
  • How to visualise results

For this example we use force by mixing layer depth and surface photosynthetically available radiation (PAR) data from the Mercator Ocean model and NASA VIIRS observations.

Install dependencies

First we check we have the dependencies installed

using Pkg
pkg"add OceanBioME, Oceananigans, NetCDF, Interpolations, DataDeps, CairoMakie"

Model setup

First load the required packages

using OceanBioME
using Oceananigans, Random, Printf, NetCDF, Interpolations, DataDeps
using Oceananigans.Units
using Oceananigans.Fields: FunctionField

import Oceananigans.TurbulenceClosures: maximum_numeric_diffusivity

maximum_numeric_diffusivity(κ::NamedTuple) = maximum(maximum.(values(κ)))
maximum_numeric_diffusivity(κ::FunctionField) = maximum(κ)

const year = years = 365days # just for these idealised cases

Load external forcing data

Loading the forcing data from our online copy

dd = DataDep(
    "example_data",
    "example data from subpolar re-analysis and observational products",
    "https://github.com/OceanBioME/OceanBioME_example_data/raw/main/subpolar.nc"
)
register(dd)
filename = datadep"example_data/subpolar.nc"
times = ncread(filename, "time")
temp = ncread(filename, "temp")
salinity = ncread(filename, "so")
mld = ncread(filename, "mld")
par = ncread(filename, "par")

temperature_itp = LinearInterpolation(times, temp)
salinity_itp = LinearInterpolation(times, salinity)
mld_itp = LinearInterpolation(times, mld)
PAR_itp = LinearInterpolation(times, par)

t_function(x, y, z, t) = temperature_itp(mod(t, 364days))
s_function(x, y, z, t) = salinity_itp(mod(t, 364days))
surface_PAR(x, y, t) = PAR_itp(mod(t, 364days))
κₜ(x, y, z, t) = 2e-2 * max(1 - (z + mld_itp(mod(t, 364days)) / 2)^2 / (mld_itp(mod(t, 364days)) / 2)^2, 0) + 1e-4

Grid and diffusivity field

Define the grid (in this case a non uniform grid for better resolution near the surface) and an extra Oceananigans field for the PAR to be stored in

Nz = 33
Lz = 600meters
refinement = 10
stretching = 5.754
h(k) = (k - 1) / Nz
ζ₀(k) = 1 + (h(k) - 1) / refinement
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

grid = RectilinearGrid(size = (1, 1, Nz), x = (0, 20meters), y = (0, 20meters), z = z_faces)

clock = Clock(; time = 0.0)

κ = FunctionField{Center, Center, Center}(κₜ, grid; clock)
FunctionField located at (Center, Center, Center)
├── func: κₜ (generic function with 1 method)
├── grid: 1×1×33 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 1×1×3 halo
├── clock: Clock{Float64, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days)
└── parameters: nothing

Biogeochemical and Oceananigans model

Here we instantiate the LOBSTER model with carbonate chemistry and a surface flux of DIC (CO₂)

biogeochemistry = LOBSTER(; grid,
                            surface_photosynthetically_active_radiation = surface_PAR,
                            carbonates = true,
                            scale_negatives = true)

CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition()

T = FunctionField{Center, Center, Center}(t_function, grid; clock)
S = FunctionField{Center, Center, Center}(s_function, grid; clock)

model = NonhydrostaticModel(; grid, clock,
                              closure = ScalarDiffusivity(ν = κ, κ = κ),
                              biogeochemistry,
                              boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),),
                              auxiliary_fields = (; T, S))

set!(model, P = 0.03, Z = 0.03, NO₃ = 11.0, NH₄ = 0.05, DIC = 2200.0, Alk = 2400.0)

Simulation

Next we setup the simulation along with some callbacks that:

  • Show the progress of the simulation
  • Store the output
  • Prevent the tracers from going negative from numerical error (see discussion of this in the positivity preservation implementation page)
  • Adapt the timestep length to reduce the run time
simulation = Simulation(model, Δt = 1minutes, stop_time = 100days)

progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n",
                                iteration(sim),
                                prettytime(sim),
                                prettytime(sim.Δt),
                                prettytime(sim.run_wall_time))

simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(500))

filename = "data_forced"
simulation.output_writers[:profiles] = JLD2OutputWriter(model,
                                                        model.tracers,
                                                        filename = "$filename.jld2",
                                                        schedule = TimeInterval(1day),
                                                        overwrite_existing = true)

wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2,
                        max_change = 1.5, min_change = 0.75,
                        cell_advection_timescale = column_advection_timescale)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
┌ Warning: Attempting to store typeof(Main.var"##234".κₜ).
│ JLD2 only stores functions by name.
│  This may not be useful for anonymous functions.
└ @ JLD2 ~/.julia/packages/JLD2/phaon/src/data/writing_datatypes.jl:447
┌ Warning: Attempting to store typeof(Main.var"##234".κₜ).
│ JLD2 only stores functions by name.
│  This may not be useful for anonymous functions.
└ @ JLD2 ~/.julia/packages/JLD2/phaon/src/data/writing_datatypes.jl:447

Run!

We are ready to run the simulation

run!(simulation)
[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 1 minute, wall time: 0 seconds
[ Info:     ... simulation initialization complete (29.439 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (7.356 seconds).
Iteration: 0500, time: 6.598 hours, Δt: 47.460 seconds, wall time: 0 seconds
Iteration: 1000, time: 13.187 hours, Δt: 47.420 seconds, wall time: 0 seconds
Iteration: 1500, time: 19.771 hours, Δt: 47.386 seconds, wall time: 0 seconds
Iteration: 2000, time: 1.098 days, Δt: 47.380 seconds, wall time: 0 seconds
Iteration: 2500, time: 1.372 days, Δt: 47.427 seconds, wall time: 0 seconds
Iteration: 3000, time: 1.647 days, Δt: 47.487 seconds, wall time: 0 seconds
Iteration: 3500, time: 1.921 days, Δt: 47.481 seconds, wall time: 0 seconds
Iteration: 4000, time: 2.196 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 4500, time: 2.470 days, Δt: 47.360 seconds, wall time: 0 seconds
Iteration: 5000, time: 2.744 days, Δt: 47.449 seconds, wall time: 0 seconds
Iteration: 5500, time: 3.018 days, Δt: 47.405 seconds, wall time: 0 seconds
Iteration: 6000, time: 3.292 days, Δt: 47.445 seconds, wall time: 0 seconds
Iteration: 6500, time: 3.566 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 7000, time: 3.841 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 7500, time: 4.115 days, Δt: 47.427 seconds, wall time: 0 seconds
Iteration: 8000, time: 4.389 days, Δt: 47.347 seconds, wall time: 0 seconds
Iteration: 8500, time: 4.663 days, Δt: 47.314 seconds, wall time: 0 seconds
Iteration: 9000, time: 4.937 days, Δt: 47.320 seconds, wall time: 0 seconds
Iteration: 9500, time: 5.210 days, Δt: 47.335 seconds, wall time: 0 seconds
Iteration: 10000, time: 5.484 days, Δt: 47.349 seconds, wall time: 0 seconds
Iteration: 10500, time: 5.758 days, Δt: 47.366 seconds, wall time: 0 seconds
Iteration: 11000, time: 6.032 days, Δt: 47.388 seconds, wall time: 0 seconds
Iteration: 11500, time: 6.307 days, Δt: 47.442 seconds, wall time: 0 seconds
Iteration: 12000, time: 6.581 days, Δt: 47.510 seconds, wall time: 0 seconds
Iteration: 12500, time: 6.856 days, Δt: 47.494 seconds, wall time: 0 seconds
Iteration: 13000, time: 7.131 days, Δt: 47.500 seconds, wall time: 0 seconds
Iteration: 13500, time: 7.406 days, Δt: 47.481 seconds, wall time: 0 seconds
Iteration: 14000, time: 7.680 days, Δt: 47.401 seconds, wall time: 0 seconds
Iteration: 14500, time: 7.954 days, Δt: 47.345 seconds, wall time: 0 seconds
Iteration: 15000, time: 8.228 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 15500, time: 8.502 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 16000, time: 8.776 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 16500, time: 9.049 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 17000, time: 9.323 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 17500, time: 9.597 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 18000, time: 9.871 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 18500, time: 10.145 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 19000, time: 10.419 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 19500, time: 10.693 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 20000, time: 10.966 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 20500, time: 11.240 days, Δt: 47.338 seconds, wall time: 0 seconds
Iteration: 21000, time: 11.514 days, Δt: 47.340 seconds, wall time: 0 seconds
Iteration: 21500, time: 11.788 days, Δt: 47.341 seconds, wall time: 0 seconds
Iteration: 22000, time: 12.061 days, Δt: 47.330 seconds, wall time: 0 seconds
Iteration: 22500, time: 12.335 days, Δt: 47.314 seconds, wall time: 0 seconds
Iteration: 23000, time: 12.609 days, Δt: 47.376 seconds, wall time: 0 seconds
Iteration: 23500, time: 12.884 days, Δt: 47.531 seconds, wall time: 0 seconds
Iteration: 24000, time: 13.158 days, Δt: 47.431 seconds, wall time: 0 seconds
Iteration: 24500, time: 13.432 days, Δt: 47.385 seconds, wall time: 0 seconds
Iteration: 25000, time: 13.707 days, Δt: 47.349 seconds, wall time: 0 seconds
Iteration: 25500, time: 13.981 days, Δt: 47.325 seconds, wall time: 0 seconds
Iteration: 26000, time: 14.254 days, Δt: 47.332 seconds, wall time: 0 seconds
Iteration: 26500, time: 14.528 days, Δt: 47.345 seconds, wall time: 0 seconds
Iteration: 27000, time: 14.802 days, Δt: 47.360 seconds, wall time: 0 seconds
Iteration: 27500, time: 15.076 days, Δt: 47.408 seconds, wall time: 0 seconds
Iteration: 28000, time: 15.351 days, Δt: 47.467 seconds, wall time: 0 seconds
Iteration: 28500, time: 15.625 days, Δt: 47.334 seconds, wall time: 0 seconds
Iteration: 29000, time: 15.899 days, Δt: 47.317 seconds, wall time: 0 seconds
Iteration: 29500, time: 16.173 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 30000, time: 16.447 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 30500, time: 16.720 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 31000, time: 16.994 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 31500, time: 17.268 days, Δt: 47.333 seconds, wall time: 0 seconds
Iteration: 32000, time: 17.542 days, Δt: 47.329 seconds, wall time: 0 seconds
Iteration: 32500, time: 17.816 days, Δt: 47.326 seconds, wall time: 0 seconds
Iteration: 33000, time: 18.089 days, Δt: 47.323 seconds, wall time: 0 seconds
Iteration: 33500, time: 18.363 days, Δt: 47.321 seconds, wall time: 0 seconds
Iteration: 34000, time: 18.637 days, Δt: 47.319 seconds, wall time: 0 seconds
Iteration: 34500, time: 18.911 days, Δt: 47.318 seconds, wall time: 0 seconds
Iteration: 35000, time: 19.185 days, Δt: 47.320 seconds, wall time: 0 seconds
Iteration: 35500, time: 19.458 days, Δt: 47.325 seconds, wall time: 0 seconds
Iteration: 36000, time: 19.732 days, Δt: 47.330 seconds, wall time: 0 seconds
Iteration: 36500, time: 20.006 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 37000, time: 20.280 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 37500, time: 20.554 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 38000, time: 20.828 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 38500, time: 21.101 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 39000, time: 21.375 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 39500, time: 21.649 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 40000, time: 21.923 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 40500, time: 22.197 days, Δt: 47.335 seconds, wall time: 0 seconds
Iteration: 41000, time: 22.471 days, Δt: 47.332 seconds, wall time: 0 seconds
Iteration: 41500, time: 22.745 days, Δt: 47.329 seconds, wall time: 0 seconds
Iteration: 42000, time: 23.018 days, Δt: 47.328 seconds, wall time: 0 seconds
Iteration: 42500, time: 23.292 days, Δt: 47.400 seconds, wall time: 0 seconds
Iteration: 43000, time: 23.567 days, Δt: 47.522 seconds, wall time: 0 seconds
Iteration: 43500, time: 23.842 days, Δt: 47.420 seconds, wall time: 0 seconds
Iteration: 44000, time: 24.116 days, Δt: 47.336 seconds, wall time: 0 seconds
Iteration: 44500, time: 24.389 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 45000, time: 24.663 days, Δt: 47.345 seconds, wall time: 0 seconds
Iteration: 45500, time: 24.938 days, Δt: 47.428 seconds, wall time: 0 seconds
Iteration: 46000, time: 25.211 days, Δt: 47.461 seconds, wall time: 0 seconds
Iteration: 46500, time: 25.486 days, Δt: 47.468 seconds, wall time: 0 seconds
Iteration: 47000, time: 25.761 days, Δt: 47.475 seconds, wall time: 0 seconds
Iteration: 47500, time: 26.035 days, Δt: 47.482 seconds, wall time: 0 seconds
Iteration: 48000, time: 26.310 days, Δt: 47.481 seconds, wall time: 0 seconds
Iteration: 48500, time: 26.585 days, Δt: 47.481 seconds, wall time: 0 seconds
Iteration: 49000, time: 26.859 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 49500, time: 27.134 days, Δt: 47.481 seconds, wall time: 0 seconds
Iteration: 50000, time: 27.409 days, Δt: 47.483 seconds, wall time: 0 seconds
Iteration: 50500, time: 27.684 days, Δt: 47.485 seconds, wall time: 0 seconds
Iteration: 51000, time: 27.958 days, Δt: 47.487 seconds, wall time: 0 seconds
Iteration: 51500, time: 28.233 days, Δt: 47.489 seconds, wall time: 0 seconds
Iteration: 52000, time: 28.508 days, Δt: 47.491 seconds, wall time: 0 seconds
Iteration: 52500, time: 28.783 days, Δt: 47.492 seconds, wall time: 0 seconds
Iteration: 53000, time: 29.057 days, Δt: 47.497 seconds, wall time: 0 seconds
Iteration: 53500, time: 29.332 days, Δt: 47.515 seconds, wall time: 0 seconds
Iteration: 54000, time: 29.607 days, Δt: 47.533 seconds, wall time: 0 seconds
Iteration: 54500, time: 29.882 days, Δt: 47.550 seconds, wall time: 0 seconds
Iteration: 55000, time: 30.157 days, Δt: 47.502 seconds, wall time: 0 seconds
Iteration: 55500, time: 30.432 days, Δt: 47.443 seconds, wall time: 0 seconds
Iteration: 56000, time: 30.706 days, Δt: 47.396 seconds, wall time: 0 seconds
Iteration: 56500, time: 30.981 days, Δt: 47.360 seconds, wall time: 0 seconds
Iteration: 57000, time: 31.254 days, Δt: 47.358 seconds, wall time: 0 seconds
Iteration: 57500, time: 31.528 days, Δt: 47.359 seconds, wall time: 0 seconds
Iteration: 58000, time: 31.802 days, Δt: 47.360 seconds, wall time: 0 seconds
Iteration: 58500, time: 32.076 days, Δt: 47.372 seconds, wall time: 0 seconds
Iteration: 59000, time: 32.350 days, Δt: 47.426 seconds, wall time: 0 seconds
Iteration: 59500, time: 32.625 days, Δt: 47.498 seconds, wall time: 0 seconds
Iteration: 60000, time: 32.900 days, Δt: 47.519 seconds, wall time: 0 seconds
Iteration: 60500, time: 33.174 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 61000, time: 33.449 days, Δt: 47.466 seconds, wall time: 0 seconds
Iteration: 61500, time: 33.724 days, Δt: 47.451 seconds, wall time: 0 seconds
Iteration: 62000, time: 33.998 days, Δt: 47.438 seconds, wall time: 0 seconds
Iteration: 62500, time: 34.272 days, Δt: 47.454 seconds, wall time: 0 seconds
Iteration: 63000, time: 34.547 days, Δt: 47.471 seconds, wall time: 0 seconds
Iteration: 63500, time: 34.822 days, Δt: 47.489 seconds, wall time: 0 seconds
Iteration: 64000, time: 35.096 days, Δt: 47.434 seconds, wall time: 0 seconds
Iteration: 64500, time: 35.370 days, Δt: 47.317 seconds, wall time: 0 seconds
Iteration: 65000, time: 35.644 days, Δt: 47.361 seconds, wall time: 0 seconds
Iteration: 65500, time: 35.919 days, Δt: 47.488 seconds, wall time: 0 seconds
Iteration: 66000, time: 36.193 days, Δt: 47.399 seconds, wall time: 0 seconds
Iteration: 66500, time: 36.467 days, Δt: 47.376 seconds, wall time: 0 seconds
Iteration: 67000, time: 36.741 days, Δt: 47.357 seconds, wall time: 0 seconds
Iteration: 67500, time: 37.015 days, Δt: 47.344 seconds, wall time: 0 seconds
Iteration: 68000, time: 37.290 days, Δt: 47.499 seconds, wall time: 0 seconds
Iteration: 68500, time: 37.564 days, Δt: 47.383 seconds, wall time: 0 seconds
Iteration: 69000, time: 37.838 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 69500, time: 38.112 days, Δt: 47.336 seconds, wall time: 0 seconds
Iteration: 70000, time: 38.386 days, Δt: 47.367 seconds, wall time: 0 seconds
Iteration: 70500, time: 38.660 days, Δt: 47.410 seconds, wall time: 0 seconds
Iteration: 71000, time: 38.935 days, Δt: 47.464 seconds, wall time: 0 seconds
Iteration: 71500, time: 39.209 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 72000, time: 39.484 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 72500, time: 39.758 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 73000, time: 40.033 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 73500, time: 40.308 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 74000, time: 40.583 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 74500, time: 40.857 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 75000, time: 41.132 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 75500, time: 41.407 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 76000, time: 41.681 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 76500, time: 41.956 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 77000, time: 42.231 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 77500, time: 42.506 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 78000, time: 42.780 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 78500, time: 43.055 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 79000, time: 43.330 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 79500, time: 43.604 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 80000, time: 43.879 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 80500, time: 44.154 days, Δt: 47.421 seconds, wall time: 0 seconds
Iteration: 81000, time: 44.428 days, Δt: 47.344 seconds, wall time: 0 seconds
Iteration: 81500, time: 44.702 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 82000, time: 44.976 days, Δt: 47.332 seconds, wall time: 0 seconds
Iteration: 82500, time: 45.249 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 83000, time: 45.523 days, Δt: 47.333 seconds, wall time: 0 seconds
Iteration: 83500, time: 45.797 days, Δt: 47.400 seconds, wall time: 0 seconds
Iteration: 84000, time: 46.071 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 84500, time: 46.346 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 85000, time: 46.621 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 85500, time: 46.896 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 86000, time: 47.170 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 86500, time: 47.445 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 87000, time: 47.720 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 87500, time: 47.995 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 88000, time: 48.269 days, Δt: 47.475 seconds, wall time: 0 seconds
Iteration: 88500, time: 48.544 days, Δt: 47.470 seconds, wall time: 0 seconds
Iteration: 89000, time: 48.819 days, Δt: 47.466 seconds, wall time: 0 seconds
Iteration: 89500, time: 49.093 days, Δt: 47.461 seconds, wall time: 0 seconds
Iteration: 90000, time: 49.367 days, Δt: 47.456 seconds, wall time: 0 seconds
Iteration: 90500, time: 49.642 days, Δt: 47.451 seconds, wall time: 0 seconds
Iteration: 91000, time: 49.917 days, Δt: 47.446 seconds, wall time: 0 seconds
Iteration: 91500, time: 50.191 days, Δt: 47.448 seconds, wall time: 0 seconds
Iteration: 92000, time: 50.466 days, Δt: 47.454 seconds, wall time: 0 seconds
Iteration: 92500, time: 50.740 days, Δt: 47.459 seconds, wall time: 0 seconds
Iteration: 93000, time: 51.015 days, Δt: 47.464 seconds, wall time: 0 seconds
Iteration: 93500, time: 51.290 days, Δt: 47.462 seconds, wall time: 0 seconds
Iteration: 94000, time: 51.564 days, Δt: 47.460 seconds, wall time: 0 seconds
Iteration: 94500, time: 51.839 days, Δt: 47.458 seconds, wall time: 0 seconds
Iteration: 95000, time: 52.113 days, Δt: 47.456 seconds, wall time: 0 seconds
Iteration: 95500, time: 52.388 days, Δt: 47.452 seconds, wall time: 0 seconds
Iteration: 96000, time: 52.662 days, Δt: 47.449 seconds, wall time: 0 seconds
Iteration: 96500, time: 52.937 days, Δt: 47.445 seconds, wall time: 0 seconds
Iteration: 97000, time: 53.211 days, Δt: 47.445 seconds, wall time: 0 seconds
Iteration: 97500, time: 53.486 days, Δt: 47.445 seconds, wall time: 0 seconds
Iteration: 98000, time: 53.761 days, Δt: 47.445 seconds, wall time: 0 seconds
Iteration: 98500, time: 54.035 days, Δt: 47.445 seconds, wall time: 0 seconds
Iteration: 99000, time: 54.309 days, Δt: 47.451 seconds, wall time: 0 seconds
Iteration: 99500, time: 54.584 days, Δt: 47.457 seconds, wall time: 0 seconds
Iteration: 100000, time: 54.858 days, Δt: 47.463 seconds, wall time: 0 seconds
Iteration: 100500, time: 55.133 days, Δt: 47.468 seconds, wall time: 0 seconds
Iteration: 101000, time: 55.408 days, Δt: 47.472 seconds, wall time: 0 seconds
Iteration: 101500, time: 55.682 days, Δt: 47.475 seconds, wall time: 0 seconds
Iteration: 102000, time: 55.957 days, Δt: 47.479 seconds, wall time: 0 seconds
Iteration: 102500, time: 56.232 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 103000, time: 56.507 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 103500, time: 56.781 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 104000, time: 57.056 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 104500, time: 57.331 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 105000, time: 57.606 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 105500, time: 57.880 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 106000, time: 58.155 days, Δt: 47.457 seconds, wall time: 0 seconds
Iteration: 106500, time: 58.429 days, Δt: 47.418 seconds, wall time: 0 seconds
Iteration: 107000, time: 58.704 days, Δt: 47.385 seconds, wall time: 0 seconds
Iteration: 107500, time: 58.978 days, Δt: 47.358 seconds, wall time: 0 seconds
Iteration: 108000, time: 59.252 days, Δt: 47.379 seconds, wall time: 0 seconds
Iteration: 108500, time: 59.526 days, Δt: 47.411 seconds, wall time: 0 seconds
Iteration: 109000, time: 59.800 days, Δt: 47.448 seconds, wall time: 0 seconds
Iteration: 109500, time: 60.075 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 110000, time: 60.350 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 110500, time: 60.624 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 111000, time: 60.899 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 111500, time: 61.174 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 112000, time: 61.448 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 112500, time: 61.723 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 113000, time: 61.998 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 113500, time: 62.273 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 114000, time: 62.547 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 114500, time: 62.822 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 115000, time: 63.097 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 115500, time: 63.371 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 116000, time: 63.646 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 116500, time: 63.921 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 117000, time: 64.196 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 117500, time: 64.470 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 118000, time: 64.745 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 118500, time: 65.020 days, Δt: 47.472 seconds, wall time: 0 seconds
Iteration: 119000, time: 65.294 days, Δt: 47.353 seconds, wall time: 0 seconds
Iteration: 119500, time: 65.568 days, Δt: 47.311 seconds, wall time: 0 seconds
Iteration: 120000, time: 65.842 days, Δt: 47.364 seconds, wall time: 0 seconds
Iteration: 120500, time: 66.116 days, Δt: 47.453 seconds, wall time: 0 seconds
Iteration: 121000, time: 66.390 days, Δt: 47.345 seconds, wall time: 0 seconds
Iteration: 121500, time: 66.664 days, Δt: 47.332 seconds, wall time: 0 seconds
Iteration: 122000, time: 66.939 days, Δt: 47.455 seconds, wall time: 0 seconds
Iteration: 122500, time: 67.213 days, Δt: 47.313 seconds, wall time: 0 seconds
Iteration: 123000, time: 67.487 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 123500, time: 67.761 days, Δt: 47.355 seconds, wall time: 0 seconds
Iteration: 124000, time: 68.035 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 124500, time: 68.310 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 125000, time: 68.585 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 125500, time: 68.859 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 126000, time: 69.134 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 126500, time: 69.409 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 127000, time: 69.684 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 127500, time: 69.958 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 128000, time: 70.233 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 128500, time: 70.508 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 129000, time: 70.783 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 129500, time: 71.057 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 130000, time: 71.332 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 130500, time: 71.607 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 131000, time: 71.881 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 131500, time: 72.156 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 132000, time: 72.431 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 132500, time: 72.706 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 133000, time: 72.980 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 133500, time: 73.255 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 134000, time: 73.530 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 134500, time: 73.805 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 135000, time: 74.079 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 135500, time: 74.354 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 136000, time: 74.629 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 136500, time: 74.903 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 137000, time: 75.178 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 137500, time: 75.453 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 138000, time: 75.728 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 138500, time: 76.002 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 139000, time: 76.277 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 139500, time: 76.552 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 140000, time: 76.827 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 140500, time: 77.101 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 141000, time: 77.376 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 141500, time: 77.651 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 142000, time: 77.925 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 142500, time: 78.200 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 143000, time: 78.475 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 143500, time: 78.750 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 144000, time: 79.024 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 144500, time: 79.299 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 145000, time: 79.574 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 145500, time: 79.848 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 146000, time: 80.123 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 146500, time: 80.398 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 147000, time: 80.673 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 147500, time: 80.947 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 148000, time: 81.222 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 148500, time: 81.497 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 149000, time: 81.772 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 149500, time: 82.046 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 150000, time: 82.321 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 150500, time: 82.596 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 151000, time: 82.870 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 151500, time: 83.145 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 152000, time: 83.420 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 152500, time: 83.695 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 153000, time: 83.969 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 153500, time: 84.244 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 154000, time: 84.519 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 154500, time: 84.794 days, Δt: 47.480 seconds, wall time: 0 seconds
Iteration: 155000, time: 85.068 days, Δt: 47.479 seconds, wall time: 0 seconds
Iteration: 155500, time: 85.343 days, Δt: 47.477 seconds, wall time: 0 seconds
Iteration: 156000, time: 85.618 days, Δt: 47.475 seconds, wall time: 0 seconds
Iteration: 156500, time: 85.892 days, Δt: 47.473 seconds, wall time: 0 seconds
Iteration: 157000, time: 86.167 days, Δt: 47.321 seconds, wall time: 0 seconds
Iteration: 157500, time: 86.441 days, Δt: 47.321 seconds, wall time: 0 seconds
Iteration: 158000, time: 86.715 days, Δt: 47.392 seconds, wall time: 0 seconds
Iteration: 158500, time: 86.989 days, Δt: 47.355 seconds, wall time: 0 seconds
Iteration: 159000, time: 87.263 days, Δt: 47.459 seconds, wall time: 0 seconds
Iteration: 159500, time: 87.538 days, Δt: 47.489 seconds, wall time: 0 seconds
Iteration: 160000, time: 87.812 days, Δt: 47.528 seconds, wall time: 0 seconds
Iteration: 160500, time: 88.085 days, Δt: 47.325 seconds, wall time: 0 seconds
Iteration: 161000, time: 88.359 days, Δt: 47.316 seconds, wall time: 0 seconds
Iteration: 161500, time: 88.633 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 162000, time: 88.907 days, Δt: 47.313 seconds, wall time: 0 seconds
Iteration: 162500, time: 89.180 days, Δt: 47.346 seconds, wall time: 0 seconds
Iteration: 163000, time: 89.455 days, Δt: 47.465 seconds, wall time: 0 seconds
Iteration: 163500, time: 89.729 days, Δt: 47.448 seconds, wall time: 0 seconds
Iteration: 164000, time: 90.003 days, Δt: 47.340 seconds, wall time: 0 seconds
Iteration: 164500, time: 90.277 days, Δt: 47.339 seconds, wall time: 0 seconds
Iteration: 165000, time: 90.551 days, Δt: 47.338 seconds, wall time: 0 seconds
Iteration: 165500, time: 90.825 days, Δt: 47.338 seconds, wall time: 0 seconds
Iteration: 166000, time: 91.098 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 166500, time: 91.372 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 167000, time: 91.646 days, Δt: 47.338 seconds, wall time: 0 seconds
Iteration: 167500, time: 91.920 days, Δt: 47.338 seconds, wall time: 0 seconds
Iteration: 168000, time: 92.193 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 168500, time: 92.467 days, Δt: 47.363 seconds, wall time: 0 seconds
Iteration: 169000, time: 92.742 days, Δt: 47.527 seconds, wall time: 0 seconds
Iteration: 169500, time: 93.016 days, Δt: 47.352 seconds, wall time: 0 seconds
Iteration: 170000, time: 93.290 days, Δt: 47.333 seconds, wall time: 0 seconds
Iteration: 170500, time: 93.564 days, Δt: 47.319 seconds, wall time: 0 seconds
Iteration: 171000, time: 93.838 days, Δt: 47.400 seconds, wall time: 0 seconds
Iteration: 171500, time: 94.112 days, Δt: 47.404 seconds, wall time: 0 seconds
Iteration: 172000, time: 94.386 days, Δt: 47.383 seconds, wall time: 0 seconds
Iteration: 172500, time: 94.660 days, Δt: 47.376 seconds, wall time: 0 seconds
Iteration: 173000, time: 94.934 days, Δt: 47.830 seconds, wall time: 0 seconds
Iteration: 173500, time: 95.210 days, Δt: 47.431 seconds, wall time: 0 seconds
Iteration: 174000, time: 95.484 days, Δt: 47.442 seconds, wall time: 0 seconds
Iteration: 174500, time: 95.758 days, Δt: 47.337 seconds, wall time: 0 seconds
Iteration: 175000, time: 96.032 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 175500, time: 96.306 days, Δt: 47.348 seconds, wall time: 0 seconds
Iteration: 176000, time: 96.580 days, Δt: 47.354 seconds, wall time: 0 seconds
Iteration: 176500, time: 96.854 days, Δt: 47.323 seconds, wall time: 0 seconds
Iteration: 177000, time: 97.129 days, Δt: 47.430 seconds, wall time: 0 seconds
Iteration: 177500, time: 97.403 days, Δt: 47.352 seconds, wall time: 0 seconds
Iteration: 178000, time: 97.677 days, Δt: 47.316 seconds, wall time: 0 seconds
Iteration: 178500, time: 97.950 days, Δt: 47.316 seconds, wall time: 0 seconds
Iteration: 179000, time: 98.223 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 179500, time: 98.497 days, Δt: 47.320 seconds, wall time: 0 seconds
Iteration: 180000, time: 98.771 days, Δt: 47.351 seconds, wall time: 0 seconds
Iteration: 180500, time: 99.045 days, Δt: 47.497 seconds, wall time: 0 seconds
Iteration: 181000, time: 99.319 days, Δt: 47.312 seconds, wall time: 0 seconds
Iteration: 181500, time: 99.593 days, Δt: 47.311 seconds, wall time: 0 seconds
Iteration: 182000, time: 99.868 days, Δt: 47.369 seconds, wall time: 0 seconds
[ Info: Simulation is stopping after running for 0 seconds.
[ Info: Simulation time 100.000 days equals or exceeds stop time 100 days.

Load output and plot

Now we can visualise the results with some post processing to diagnose the air-sea CO₂ flux

   P = FieldTimeSeries("$filename.jld2", "P")
 NO₃ = FieldTimeSeries("$filename.jld2", "NO₃")
   Z = FieldTimeSeries("$filename.jld2", "Z")
sPOM = FieldTimeSeries("$filename.jld2", "sPOM")
bPOM = FieldTimeSeries("$filename.jld2", "bPOM")
 DIC = FieldTimeSeries("$filename.jld2", "DIC")
 Alk = FieldTimeSeries("$filename.jld2", "Alk")

x, y, z = nodes(P)
times = P.times

We compute the air-sea CO₂ flux at the surface (corresponding to vertical index k = grid.Nz) and the carbon export by computing how much carbon sinks below some arbitrary depth; here we use depth that corresponds to k = grid.Nz - 20.

air_sea_CO₂_flux = zeros(length(times))
carbon_export = zeros(length(times))

using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity

for (n, t) in enumerate(times)
    clock.time = t

    air_sea_CO₂_flux[n] = CO₂_flux.condition.func(1, 1, grid, clock, (; DIC = DIC[n], Alk = Alk[n], T, S))

    carbon_export[n] = (sPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] +
                        bPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry)
end

Both air_sea_CO₂_flux and carbon_export are in units mmol CO₂ / (m² s).

using CairoMakie

fig = Figure(size = (1000, 1500), fontsize = 20)

axis_kwargs = (xlabel = "Time (days)", ylabel = "z (m)", limits = ((0, times[end] / days), (-150meters, 0)))

axP = Axis(fig[1, 1]; title = "Phytoplankton concentration (mmol N/m³)", axis_kwargs...)
hmP = heatmap!(times / days, z, interior(P, 1, 1, :, :)', colormap=:batlow)
Colorbar(fig[1, 2], hmP)

axNO₃ = Axis(fig[2, 1]; title = "Nitrate concentration (mmol N/m³)", axis_kwargs...)
hmNO₃ = heatmap!(times / days, z, interior(NO₃, 1, 1, :, :)', colormap=:batlow)
Colorbar(fig[2, 2], hmNO₃)

axZ = Axis(fig[3, 1]; title = "Zooplankton concentration (mmol N/m³)", axis_kwargs...)
hmZ = heatmap!(times / days, z, interior(Z, 1, 1, :, :)', colormap=:batlow)
Colorbar(fig[3, 2], hmZ)

axD = Axis(fig[4, 1]; title = "Detritus concentration (mmol N/m³)", axis_kwargs...)
hmD = heatmap!(times / days, z, interior(sPOM, 1, 1, :, :)' .+ interior(bPOM, 1, 1, :, :)', colormap=:batlow)
Colorbar(fig[4, 2], hmD)

CO₂_molar_mass = (12 + 2 * 16) * 1e-3 # kg / mol

axfDIC = Axis(fig[5, 1], xlabel = "Time (days)", ylabel = "Flux (kgCO₂/m²/year)",
                         title = "Air-sea CO₂ flux and Sinking", limits = ((0, times[end] / days), nothing))
lines!(axfDIC, times / days, cumsum(air_sea_CO₂_flux) /1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Air-sea flux")
lines!(axfDIC, times / days, cumsum(carbon_export) /1e3    * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export")
Legend(fig[5, 2], axfDIC, framevisible = false)

fig


This page was generated using Literate.jl.