Roadmap and ways to contribute

AIBECS is maintained by a few people. A lot of the machinery is in place, but several directions remain wide open. This page lists the ones we think would have the biggest impact, partly to record where we'd like to take the package, and partly as an invitation: if any of these resonate, please open an issue or start a PR.

Parameter estimation beyond the current example

The parameter-optimisation how-to wires AIBECS' parameter struct, F-1's cached adjoint, a shared LinearSolve.jl cache, and the Optimization.jl frontend together by hand. It works, but it is one of several plausible designs. A few we'd like to see explored:

• Probabilistic programming via Turing.jl. AIBECS lets you declare priors with @prior (or by overloading prior directly), and the optimiser in the how-to minimises -log posterior up to a constant. But this is bespoke code and not a careful, expert-written Bayesian workflow. Wiring AIBECS parameters into a Turing @model would unlock full Bayesian inference (MCMC, variational families) over the steady-state solution, with the same priors driving both regimes. The F-1 gradient/Hessian would compose naturally with Turing's HMC / NUTS samplers. A small worked example fitting a single tracer would be a great first contribution.

• Lean further on SciML parameter-estimation idioms. The SciML stack (SciMLSensitivity, DiffEqParamEstim, Optimization.jl) has a well-trodden pipeline for fitting ODEProblem / SteadyStateProblem parameters from data. Right now AIBECS bolts F1Method.jl onto Optimization by hand. Restructuring the AIBECS objective surface so it slots into SciML's build_loss_objective-style helpers would let users mix AIBECS models with the rest of the SciML param-fitting ecosystem without bespoke glue.

Alternative sensitivity backends

The how-to uses F1Method because it caches the steady-state factors and exposes both gradient and Hessian. Two natural alternatives that would compose well with the AIBECS objective surface:

• SciMLSensitivity.SteadyStateAdjoint computes a per-call adjoint gradient with no cross-call caching and no Hessian — appropriate when the steady-state solve dominates and you only need the gradient.

• ImplicitDifferentiation.jl provides general implicit-function differentiation; no F-1 Hessian today, but the package is actively developed.

Either would be a meaningful upstream PR — slotting underneath the AIBECS objective / gradient surface and letting users pick the sensitivity strategy at solve time.

Simpler parameter representation

The current AbstractParameters machinery stacks several macros from FieldMetadata.jl (@units, @prior, @description, @bounds, @logscaled, @flattenable, @reference, plus an @initial_value we define ourselves) on top of a plain struct. It is expressive — every field carries units, a prior, a description, bounds, etc. — but the macro layering has real costs:

• Reading and writing a parameter value goes through the FieldMetadata/Flatten plumbing rather than ordinary struct access, which makes it hard to plug AIBECS parameters straight into packages that expect a NamedTuple, ComponentArrays.jl, or a Turing @model.

• The metadata is attached at the type level via macros, so building a parameter set programmatically (e.g. from a CSV, or generated by another package) is awkward.

• New contributors have to learn the macro stack before they can add a parameter.

A nice direction would be a flatter representation — one small Parameter value per field carrying value, unit, prior, bounds, logscaled, description, reference, gathered into a NamedTuple or ComponentArray. Standard introspection (propertynames, getproperty) then works without macros, the optimisable subset is a filter rather than a @flattenable annotation, and downstream packages can read/write parameters with no AIBECS-specific glue. The existing pretty-printing, CSV round-tripping, and unit conversions would need to be re-implemented on top, but the surface area shrinks considerably.

More circulations

AIBECS currently ships OCIM0/1/2/2-48L, OCCA, Archer et al. (2000), the two-box and Primeau 2×2×2 toy circulations, and Haine & Hall (2025). Two natural extensions:

• More published transport matrices. ECCO-derived and CMIP-derived matrices, ACCESS-OM2, MOM6, and any of the offline products coming out of OceanTransportMatrixBuilder.jl would each be a self-contained PR — see the existing OCIM*.jl / OCCA.jl files for the shape (download via DataDeps, expose T, grid, wet3D, and the associated diagnostics).

• A generic loader. Most additions today copy an existing OCIM2.jl-style file. A small load_circulation(path) that reads a standard layout (transport operator + grid + masks) would cut the boilerplate and make user-contributed circulations a one-file affair.

Regional circulation models

All circulations shipped today are global. Ingesting natively regional transport matrices — e.g. from regional MOM6/ROMS runs with explicit open boundary conditions — would let AIBECS run biogeochemistry on a subdomain at much finer resolution than any global matrix can offer. Schoonover et al. (2023) is a good reference for the kind of regional (and high-resolution) offline transport this would target. The main scaffolding work is generalising boundary-condition handling so a prescribed open-boundary inflow (possibly time-varying) becomes a first-class concept alongside the existing surface and sediment boundaries.

High-resolution circulation models

Whether AIBECS itself should handle eddy-resolving (¼° or finer) global circulations is an open question — at those resolutions direct factorisation is out of reach and GPU-resident iterative solvers become mandatory (see also the iterative-solvers section). The ACCESS-OM2 × Oceananigans project demonstrates feasibility: GPU + iterative solves of single tracers at ¼° in a periodic setting, with a target of 0.1° × 75 levels — roughly 1200× the size of OCIM2. Even if the full eddy-resolving regime ends up living outside AIBECS, the choices made there (matrix storage, solver/preconditioner pairs, GPU offload) should feed back into AIBECS once iterative solvers land.

Coarsening via lump-and-spray

OceanTransportMatrixBuilder.jl provides lump_and_spray(wet3D, vol, T; di, dj, dk) which returns LUMP, SPRAY, and a coarsened volume vector such that LUMP * T * SPRAY is a volume-conserving coarsened transport operator (and LUMP * x / SPRAY * x_c move tracer fields between the two grids). Wiring this into AIBECS so any shipped circulation can be coarsened on demand would unlock fast iteration on parameter sweeps and sensitivity studies before refining on the native grid.

Non-ocean components

Several biogeochemical problems need state variables that aren't 3D ocean tracers:

• Slab atmospheres. A single well-mixed atmospheric reservoir for CO₂ (or O₂) is required to close the carbon cycle — pCO₂(atm) evolves in time and exchanges with surface DIC via air–sea flux. See PCO2 model for an example reference where pCO2atm is carried as a 0D state variable alongside the 3D tracers.

• Surface-mean intermediates. Quantities like a surface-mean DIC or ALK (like the DICsrf/ALKsrf variables in the PCO2 model) are useful to apply precipitation/evaporation virtual fluxes consistently across the surface without recomputing the mean inside every right-hand-side call.

• Other "auxiliary" state. Sediment pools, ice reservoirs, simple land boxes — anything that couples to the 3D field through a small number of boundary fluxes.

The common requirement is letting AIBECS's AIBECSFunction carry state vectors that aren't shaped like a tracer-on-the-grid, with the Jacobian assembly extended to handle the resulting block-sparse couplings (3D ↔ 2D surface, 3D ↔ 0D atmosphere). Writing the automatic scaffolding for building the Jacobian rows/columns for these different-sized blocks is the bulk of the work.

Iterative linear solvers and preconditioners

Steady-state solves currently rely on direct sparse solvers on the CPU LinearSolve.jl (UMFPACK by default). That is the right choice for OCIM-sized problems, but it does not scale well, as both memory and factorisation time grow quickly with increasing resolution, and there is no GPU path. Nothing in the call chain — SteadyStateProblem → AIBECS state function → linear solve — is inherently CPU-only, so a GPU-aware factorisation like CUDSS.jl should in principle drop in, but this path has not been exercised. Krylov-subspace solvers (GMRES, BiCGStab) from Krylov.jl or IterativeSolvers.jl would address both, but only with effective preconditioning: ILU(0), algebraic multigrid (via AlgebraicMultigrid.jl), or physics-based block preconditioners that exploit the advection–diffusion structure. A useful first step is to benchmark a few solver/preconditioner pairs against the existing LU baseline on OCIM2 and the 48-layer variant; from there, exposing the choice through the solve(...; alg=...) interface follows the standard LinearSolve pattern.

Time stepping

AIBECS is built around steady-state solves. A first-class time-stepping interface — integrating the same AIBECSFunction forward in time via the SciML ODE stack — would unlock transient runs, spin-up diagnostics, seasonal forcings, and natural coupling with non-steady components (slab atmospheres, sediments). Open design questions: how to expose the transport split (linear T + nonlinear sources) cleanly to IMEX schemes, how to share Jacobian assembly with the existing steady-state solver, and how to use a precomputed steady state as the initial condition. This entry is meant to open a discussion, not prescribe a design.

Pseudo-transient solvers for stiff problems

Newton's method on F(x) = 0 works well when the initial guess is good and the Jacobian is well-conditioned, but stiff source/sink terms, strongly nonlinear closures, or multiple co-existing equilibria can make it stall or pick the wrong branch. A robust fallback is pseudo-transient continuation: integrate dx/dt = F(x) toward steady state with an implicit scheme whose time step starts small and grows as the residual drops, approaching a Newton step in the limit. SIAMFANLEquations.jl already implements it (and Newton-Krylov cousins) and is a natural backend; SciML's DynamicSS family is another option. Either could slot in as an alternative algorithm under solve(SteadyStateProblem(...)).

Periodic (seasonally cyclic) solvers

Steady-state solves assume time-invariant forcing. A complementary regime — much cheaper than fully time-stepping — is the periodic steady state: given a seasonally-cyclic transport (e.g. 12 monthly matrices) and seasonally-cyclic forcings, solve directly for the periodic orbit. Two references worth building on:

• CYCLOCIM (Huang, Primeau & DeVries, 2021) — a 4-D variational assimilation system for the climatological seasonal cycle of the ocean circulation, producing 12 monthly transport matrices.

• Pasquier et al. (2025) — uses a periodic-solve framework to compute deep-ocean transit-time distributions and sequestration efficiency offline from climate-model archives.

A solve_periodic interface alongside the existing steady-state path (and any future time-stepping path) would open seasonally-resolved BGC studies that currently force users into long transient integrations.

Circulation diagnostic tutorials

AIBECS already has everything needed to diagnose a circulation; three short tutorials would make this much more visible:

• Water-mass fractions. Solve for a passive tracer with a Dirichlet/restoring boundary condition (= 1 in a chosen surface region, = 0 elsewhere); the steady interior values give the fraction of each gridcell last in contact with the surface in that region.

• Adjoint tracers and mean time to reemergence. Build the true adjoint transport V⁻¹ Tᵀ V (the adjoint with respect to the volume-weighted inner product ⟨x, y⟩ = xᵀ V y, the discrete equivalent of ∫ x y dV), with a unit source everywhere and a fast surface relaxation as sink; the steady solution is the mean time to reemergence. The forward analogue with the same source/sink gives the ideal age.

• Ventilation per unit area. From the age a and reemergence-time r fields, compute calV↑ = vol · a · 1_surf / area and calV↓ = vol · r · 1_surf / area. Pasquier et al. (2024) uses calV↓ to diagnose deoxygenation pathways.

Visualisation

• A GeoMakie-based how-to (real projections, coastlines, country borders) on top of the existing Makie how-to — the current guide deliberately stays on Cartesian Axis to keep the dependency footprint small.

• Makie recipes to mirror the existing Plots.jl recipes (plothorizontalslice, plotzonalmean, plotdepthprofile, …) so Makie users don't have to compute slices / means by hand.

General

• More how-tos, especially around coupling AIBECS to external forcings, working with custom grids, or porting a published model.

• More tutorials: anything that walks through a real published setup end-to-end will help new users a lot.

• Performance: profiling reports, faster sparse Jacobian construction, GPU experiments — all welcome (see also the solver/preconditioner section above).

If you have an idea that doesn't fit any of the buckets above, that's also fine; please open an issue first and we'll figure out the right shape together.