How to determine the maximum stable time step for biogeochemistry dynamics?

[ali-ramadhan]

I was just playing around a LOBSTER box model (basically the first example in the docs) and noticed that I was getting some super negative values. Turns out it's probably numerical instability due to the time step because reducing the time step eventually leads to convergence (although sometimes blowup at intermediate time step sizes).

I'm curious if there's a way to determine the maximum stable time step for biogeochemistry dynamics (e.g. LOBSTER). This way we could take it into account when computing time step sizes for adaptive time-stepping.

If not, maybe the solution is an equivalent of the NaNChecker that stops a simulation if a tracer becomes too negative?

cc @suki-atdepthmrv @almacarolina

---

MWE:

using CairoMakie
using Oceananigans
using Oceananigans.Units
using Oceananigans.Fields: FunctionField

using OceanBioME

const year = years = 365day

grid = BoxModelGrid()
clock = Clock(time = zero(grid))

PAR⁰(t) = 100 * (1 - cos((t + 150days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2

const z = -10
PAR_func(t) = PAR⁰(t) * exp(0.2z)

PAR = FunctionField{Center, Center, Center}(PAR_func, grid; clock)

biogeochemistry = LOBSTER(;
    grid,
    light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation(PAR),
    carbonates = true,
    oxygen = true,
    variable_redfield = true,
    scale_negatives = false
)

model = BoxModel(; biogeochemistry, clock)

set!(model,
    NO₃ = 50,
    DIC = 2300,
    Alk = 2400,
    O₂ = 250,
    NH₄ = 0.1,
    DOC = 100,
    DON = 10,
    sPOC = 10,
    bPOC = 1,
    sPON = 1,
    bPON = 1,
    P = 5,
    Z = 2
)

simulation = Simulation(model; Δt = 10seconds, stop_time = 30days)

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, model.fields;
        filename = "box.jld2",
        schedule = TimeInterval(1hour),
        overwrite_existing = true
    )

prog(sim) = @info "$(prettytime(time(sim))) in $(prettytime(simulation.run_wall_time))"

simulation.callbacks[:progress] = Callback(prog, TimeInterval(1day))

@info "Running the model..."
run!(simulation)

times = FieldTimeSeries("box.jld2", "P").times

time_series = NamedTuple{keys(model.fields)}(FieldTimeSeries("box.jld2", "$field")[1, 1, 1, :] for field in keys(model.fields))

fig = Figure(size = (1200, 1200), fontsize = 24)

axs = []
for (name, tracer) in pairs(time_series)
    idx = (length(axs))
    push!(axs, Axis(fig[floor(Int, idx/2), Int(idx%2)], ylabel = "$name", xlabel = "Time (days)"))
    lines!(axs[end], times / day, tracer, linewidth = 3)
end

save("box_model.png", fig)

With a time step of 15 mins (super negative NO3, Alk, and DIC):

With a time step of 5 mins (crash and NaNs after ~0.5 days):

With a time step of 1 min (seems fine):

With a time step of 10 seconds (seems converged):

---

Environment:

OceanBioME.jl main branch with

Julia Version 1.10.7
Commit 4976d05258e (2024-11-26 15:57 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 48 × AMD Ryzen Threadripper 7960X 24-Cores
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 16 default, 0 interactive, 8 GC (on 48 virtual cores)
Environment:
  LD_PRELOAD = /usr/NX/lib/libnxegl.so

[jagoosw]

Hi Ali,

Sorry, I completely missed this before.

As far as I'm aware it's not generally possible to determine a stable timestep and it would be very involved to compute it.

This is why I came up with the ScaleNegativeTracers modifier, but I think in a box model it makes sense to error on negatives. Maybe it would be good to have a callback for that in the source code.