AIBECS interface (functions and types)

Model construction

AIBECS.AIBECSFunction Function

AIBECSFunction(T, G, nb; ∂G = nothing)
AIBECSFunction(Ts, Gs, nb, [nt, Tidx]; ∂Gs = nothing)
AIBECSFunction(Ts, Gs, nb, ::Type{P}; ∂Gs = nothing) where {P <: AbstractParameters}

Build a SciMLBase.ODEFunction for the AIBECS state equation

together with its Jacobian   , ready to plug into SteadyStateProblem or any SciML solver.

T is the linear transport operator: a sparse matrix (single tracer, constant in p) or a function p -> T(p) (depends on parameters). G is the local sources-and-sinks function G(x, p). For multi-tracer models pass tuples Ts and Gs; the optional Tidx maps each of the nt tracers to one of the operators in Ts (defaults to one operator per tracer). When a parameter type P <: AbstractParameters is supplied, the returned function also accepts a flat Vector parameter (converted internally via λ2p).

By default, ∇ₓG is built with ForwardDiff under the assumption that each ∂Gᵢ/∂xⱼ block is diagonal. To skip AD entirely, pass an analytical ∂Gs (or ∂G in the single-tracer call): a tuple-of-tuples that mirrors Gs, where ∂Gs[i][j](x₁, …, xₙₜ, p) returns the length-nb diagonal vector of . Every nt × nt block must be supplied (use zero(xⱼ) for blocks known to be zero). See check_∂Gs to validate analytical derivatives against ForwardDiff.

Notation convention: ∇ₓG, ∇ₓF, etc. denote Jacobian matrices; ∂G, ∂Gs[i][j] denote the diagonal vectors of the corresponding partial derivatives — both representations describe the same pointwise dependence, just packed differently.

T(p) = transportoperator(grd, OCIM2.load()[2])
G(x, p) = -x ./ p.τ
∂G(x, p) = -inv.(p.τ) .* one.(x)
fun = AIBECSFunction(T, G, count(iswet(grd)); ∂G = ∂G)
prob = SteadyStateProblem(fun, x₀, p)

AIBECS.check_∂Gs Function

check_∂Gs(Gs, ∂Gs, x, p, nb) -> NamedTuple

Compare user-supplied analytical derivatives ∂Gs against the ForwardDiff reference at state x and parameters p, for a model with nb boxes. Returns a NamedTuple (; maxabs, ad_jac, analytical_jac) with the overall max-abs difference and both Jacobian matrices. Use this while developing analytical derivatives; it is not called by the solver.

AIBECS.split_state_function_and_Jacobian Function

F, L, NL, ∇ₓF, ∇ₓL, ∇ₓNL, T = split_state_function_and_Jacobian(Ts, Ls, NLs, nb)

Returns the state function F and its jacobian, ∇ₓF, as well as a collection of split operators. This is experimental. Use at your own risk!

AIBECS.LinearOperators Type

LinearOperators(ops...)

Container for a tuple of linear operators that acts as their sum.

Lets AIBECS represent a transport operator built from several pieces (e.g. a sparse advection matrix plus a matrix-free diffusion operator) without materialising the sum eagerly. Supports *, \, + with sparse matrices, and in-place mul!. factorize, *, and \ materialise the sum on demand.

A = LinearOperators(T1, T2)   # behaves like T1 + T2
A * x                          # equivalent to (T1 + T2) * x
A \ b                         # equivalent to (T1 + T2) \ b

Parameters and metadata

AIBECS.AbstractParameters Type

AbstractParameters{T} <: AbstractVector{T}

An abstract type for AIBECS model parameters.

Parameters in AIBECS use the following convenience packages:

• Parameters

• FieldMetadata

• FieldDefaults

• Flatten

• Unitful

• DataFrames

• Distributions

These aim to allow for some nice features, which include

• nice syntax for unpacking parameters in functions via the @unpack macro (fron UnPack.jl)

• additional metadata on parameters

• easy conversion to and from vectors

• use of units and automatic conversions if necessary

• pretty table-format displays

• loading and saving to and from CSV files

• prior estimates for bayesian inference and optimization

See the list of examples to get an idea of how to generate parameters for your model.

Examples

Generate a simple parameter type via

julia> struct SimpleParams{T} <: AbstractParameters{T}
           α::T
           β::T
           γ::T
       end
SimpleParams

To create an instance of the SimpleParams(Float64) type, you can do

julia> p = SimpleParams(1.0, 2.0, 3.0)
SimpleParams{Float64}
│ Row │ Symbol │ Value   │
│     │ Symbol │ Float64 │
├─────┼────────┼─────────┤
│ 1   │ α      │ 1.0     │
│ 2   │ β      │ 2.0     │
│ 3   │ γ      │ 3.0     │

One of the core features from Parameters is unpacking in functions, e.g.,

julia> function simplef(p)
           @unpack α, γ = p
           return α + γ
       end
simplef (generic function with 1 method)

julia> simplef(p) # 1.0 + 3.0
4.0

More complex examples are permitted by adding metadata (thanks to FieldMetadata.jl). You can add units

julia> @units struct UnitParams{T} <: AbstractParameters{T}
           α::T | u"km"
           β::T | u"hr"
           γ::T | u"m/s"
       end ;

julia> p = UnitParams(1.0, 2.0, 3.0)
UnitParams{Float64}
│ Row │ Symbol │ Value   │ Unit     │
│     │ Symbol │ Float64 │ Unitful… │
├─────┼────────┼─────────┼──────────┤
│ 1   │ α      │ 1.0     │ km       │
│ 2   │ β      │ 2.0     │ hr       │
│ 3   │ γ      │ 3.0     │ m s^-1   │

Note that when adding units to your parameters, they will be converted to SI when unpacked, as in, e.g.,

julia> function speed(p)
           @unpack α, β, γ = p
           return α / β + γ
       end
speed (generic function with 1 method)

julia> speed(p) # (1.0 km / 2.0 hr + 3 m/s) in m/s
3.138888888888889

Another example for optimizable/flattenable parameters

julia> @initial_value @units @flattenable struct OptParams{T} <: AbstractParameters{T}
           α::T | 3.6 | u"km"  | true
           β::T | 1.0 | u"hr"  | false
           γ::T | 1.0 | u"m/s" | true
       end ;

julia> p = OptParams(initial_value(OptParams)...)
OptParams{Float64}
│ Row │ Symbol │ Value   │ Initial value │ Unit     │ Optimizable │
│     │ Symbol │ Float64 │ Float64       │ Unitful… │ Bool        │
├─────┼────────┼─────────┼───────────────┼──────────┼─────────────┤
│ 1   │ α      │ 3.6     │ 3.6           │ km       │ 1           │
│ 2   │ β      │ 1.0     │ 1.0           │ hr       │ 0           │
│ 3   │ γ      │ 1.0     │ 1.0           │ m s^-1   │ 1           │

Thanks to the FieldMetaData interface, you can chain the following preloaded metadata:

• initial_value

• units (from Unitful.jl)

• prior (from Distributions.jl)

• description (String)

• bounds (2-element Tuple)

• logscaled (Bool)

• flattenable (to convert to vectors of optimizable parameters only)

• reference (String)

Here is an example of parameter with all the possible metadata available in AIBECS:

julia> @initial_value @units @prior @description @bounds @logscaled @flattenable @reference struct FullParams{T} <: AbstractParameters{T}
           α::T | 1.0 | u"km"  | Normal(0,1)    | "The distance"   | (-Inf, Inf) | false | false | "Jean et al., 2042"
           β::T | 2.0 | u"hr"  | LogNormal(0,1) | "The time"       | (   0, Inf) | true  | true  | "Claude et al. 1983"
           γ::T | 3.0 | u"mol" | Normal(1,2)    | "The # of moles" | (  -1,   1) | false | true  | "Dusse et al. 2000"
       end ;

julia> FullParams(4.0, 5.0, 6.0)
FullParams{Float64}
│ Row │ Symbol │ Value   │ Initial value │ Unit     │ Prior                            │ Description    │ Bounds      │ Logscaled │ Optimizable │ Reference          │
│     │ Symbol │ Float64 │ Float64       │ Unitful… │ Distribu…                        │ String         │ Tuple…      │ Bool      │ Bool        │ String             │
├─────┼────────┼─────────┼───────────────┼──────────┼──────────────────────────────────┼────────────────┼─────────────┼───────────┼─────────────┼────────────────────┤
│ 1   │ α      │ 4.0     │ 1.0           │ km       │ Normal{Float64}(μ=0.0, σ=1.0)    │ The distance   │ (-Inf, Inf) │ 0         │ 0           │ Jean et al., 2042  │
│ 2   │ β      │ 5.0     │ 2.0           │ hr       │ LogNormal{Float64}(μ=0.0, σ=1.0) │ The time       │ (0, Inf)    │ 1         │ 1           │ Claude et al. 1983 │
│ 3   │ γ      │ 6.0     │ 3.0           │ mol      │ Normal{Float64}(μ=1.0, σ=2.0)    │ The # of moles │ (-1, 1)     │ 0         │ 1           │ Dusse et al. 2000  │

Note that there is no check that the metadata you give is consistent. These metadata will hopefully be useful for advanced usage of AIBECS, e.g., using prior information and/or bounds for optimization.

AIBECS.table Function

table(p)

Returns a DataFrame representation of p. Useful for printing and saving as a text/LaTeX table. Requires using DataFrames so the AIBECSDataFramesExt extension is activated.

AIBECS.latex Function

latex(p)

Returns a LaTeX-formatted table of the parameters.

Requires using DataFrames.

AIBECS.mismatch Function

mismatch(p [, k])
mismatch(::Type{T}, λ)
mismatch(d::Distribution, λ)

Parameter-space and λ-space mismatch (negative log prior).

The methods that touch a Distribution or Bijector are added by the AIBECSDistributionsExt and AIBECSBijectorsExt extensions (activated by using Distributions, Bijectors).

AIBECS.f_and_∇ₓf Function

f_and_∇ₓf(ωs, μx, σ²x, v, ωp, ::Type{P})
f_and_∇ₓf(ωs, ωp, grd, obs, ::Type{P}; kwargs...)
f_and_∇ₓf(ωs, ωp, grd, modify, obs, ::Type{P})

Build a (f, ∇ₓf) pair where f(x, λorp) is a volume-weighted mismatch objective combining tracer-vs-observation misfit with a parameter prior mismatch, and ∇ₓf is its analytical state Jacobian.

ωs weights each tracer in the sum, ωp weights the parameter-prior term, and P <: AbstractParameters is the parameter type. The first form takes volume-mean / variance / volume vectors directly; the second takes a grid plus observation tables; the third additionally lets the user transform tracers before the misfit is evaluated. Used as the model entry point for Newton-style optimisation (see the parameter-optimisation how-to).

AIBECS.p2λ Function

p2λ(p::AbstractParameters)

Converts p to a real-valued vector for optimization (referred to as the λ space in AIBECS).

Requires using Bijectors.

AIBECS.λ2p Function

λ2p(T::Type{AbstractParameters}, λ::Vector)

Converts real-valued vector λ back to a parameters object p.

Requires using Bijectors.

Tracers and state vector

AIBECS.state_to_tracers Function

state_to_tracers(x, nb, nt)
state_to_tracers(x, grd)

Split the flat state vector x into a tuple of nt tracer views of length nb.

The returned objects are views into x, so mutating them mutates x. The 2-argument form infers nb from grd and nt from length(x) / nb.

xs = state_to_tracers(x, grd)
DIP, DOP = xs                    # for a 2-tracer model

AIBECS.state_to_tracer Function

state_to_tracer(x, nb, nt, i)

View into the ith tracer of the flat state vector x (length nb*nt). Equivalent to state_to_tracers(x, nb, nt)[i].

AIBECS.tracers_to_state Function

tracers_to_state(xs)

Concatenate a tuple/array of tracer vectors xs into a flat state vector, the inverse of state_to_tracers.

AIBECS.tracer_indices Function

tracer_indices(nb, nt, i)

Range of indices in the flat state vector that correspond to the ith tracer of an nt-tracer, nb-box model.

AIBECS.unpack_tracers Function

unpack_tracers(x, grd)
unpack_tracers(x, nb, nt)

Alias for state_to_tracers with a more user-friendly name.

Circulations

AIBECS.TwoBoxModel Module

TwoBoxModel

Pedagogical 2-box (surface + deep) circulation and grid from Sarmiento and Gruber (2006), useful as a minimal sanity-check model. Use TwoBoxModel.load to obtain (grd, T).

┌────────────────────────────────┐
│                        Surface │
├────────────────────────────────┤
│                           Deep │
└────────────────────────────────┘

AIBECS.TwoBoxModel.load Method

load

Returns grd and T (in that order).

AIBECS.Primeau_2x2x2 Module

Primeau_2x2x2

Pedagogical 2×2×2 "shoebox" circulation by François Primeau, with 5 wet boxes and 3 dry boxes (see his Intro2TransportOperators notebook and boxmodel.png). Use Primeau_2x2x2.load to obtain (grd, T).

AIBECS.Primeau_2x2x2.load Method

load

Returns grd and T (in that order).

AIBECS.Archer_etal_2000 Module

Archer_etal_2000

3-box circulation of Archer et al. (2000), implemented as a regular 6-box grid so that AIBECS can build particulate sinking-flux operators on it. The boxes group as: high-latitude surface (1+3), low-latitude surface (2), deep (4+5+6). Use Archer_etal_2000.load to obtain (grd, T).

┌───────┬────────────────────────┐
│ 1(H)  │         2(L)           │
├───────┼────────────────────────┤
│ 3(H)  │         4(D)           │
├───────┼────────────────────────┤
│ 5(D)  │         6(D)           │
└───────┴────────────────────────┘

AIBECS.Archer_etal_2000.load Method

load

Returns grd and T (in that order).

AIBECS.Haine_and_Hall_2025 Module

Haine_and_Hall_2025

9-box pedagogical circulation of Haine and Hall (2002), as described in Figure 6 of Haine, Griffies, Gebbie, and Jiang (2025), A review of Green's function methods for tracer timescales and pathways in ocean models. This implementation translates the reference Julia version in OceanGreensFunctionMethods.jl (file src/pedagogical_tracer_box_model.jl) into AIBECS's T_advection / T_diffusion machinery. Use Haine_and_Hall_2025.load to obtain (grd, T).

Box layout: 1 longitude × 3 latitudes × 3 depths. Linear (column-major) indexing of the wet3D = trues(3, 1, 3) array gives:

╔════════════════╤════════════════╤════════════════╗
║                │                │                ║
║  1: H/Therm.   │  2: M/Therm.   │  3: L/Therm.   ║
║                │                │                ║
╟────────────────┼────────────────┼────────────────╢
║                │                │                ║
║  4: H/Deep     │  5: M/Deep     │  6: L/Deep     ║
║                │                │                ║
╟────────────────┼────────────────┼────────────────╢
║                │                │                ║
║  7: H/Abyss.   │  8: M/Abyss.   │  9: L/Abyss.   ║
║                │                │                ║
╚════════════════╧════════════════╧════════════════╝

• H = high latitudes

• M = mid latitudes

• L = low latitudes

• Therm. = thermocline

• Abyss. = abyssal

• Only boxes 1 and 2 are surface boxes (box 3 is considered to not exchange with the atmosphere).

Transport (all rates in Sv = 1e6 m³/s):

• abyssal overturning Ψ_a = 20 Sv on the closed loop 3 → 2 → 1 → 4 → 7 → 8 → 9 → 6 → 3

• intermediate overturning Ψ_i = 10 Sv on the closed loop 9 → 8 → 7 → 4 → 5 → 6 → 9

• vertical diffusive exchange ν = 5 Sv between thermocline–deep and deep–abyssal at every latitude.

Geometry: each box has a uniform volume of 300 Sv yr (≈ 9.5e15 m³), matching the upstream pedagogical choice. We pick realistic depths (0–3700 m, three equal layers) and equal-area latitude bands over [-80°, 0°] in the southern hemisphere, then solve for the longitude span Δlon ≈ 33° that yields the target volume.

AIBECS.Haine_and_Hall_2025.load Method

load

Returns grd and T (in that order).

AIBECS.OCIM0 Module

The OCIM0 module is used to load the OCIM v0 matrix and grid for use in AIBECS.

Tip

To load the OCIM v0 matrix and grid, do

julia> grd, T = OCIM0.load()

See DeVries and Primeau (2011) and Primeau et al. (2013) for more details.

Note

The files, that are downloaded from a public and persistent URL on Zenodo, were created with the code available at https://github.com/briochemc/OceanCirculations.

AIBECS.OCIM0.load Method

grd, T = load()

Returns the grid and the transport matrix.

Requires using JLD2 so that the AIBECSJLD2Ext extension is activated.

See DeVries and Primeau (2011) and Primeau et al. (2013) for more details.

AIBECS.OCIM1 Module

The OCIM1 module is used to load the OCIM v1 matrix and grid for use in AIBECS.

Tip

To load the OCIM v1 matrix and grid, do

julia> grd, T = OCIM1.load()

See DeVries and Primeau (2011) and Primeau et al. (2013) for more details.

Note

The files, that are downloaded from a public and persistent URL on Zenodo, were created with the code available at https://github.com/briochemc/OceanCirculations.

AIBECS.OCIM1.load Method

grd, T = load(; version=VERSIONS[1])

Returns the grid and the transport matrix.

Requires using JLD2 so that the AIBECSJLD2Ext extension is activated.

See DeVries and Primeau (2011) and DeVries (2014) for more details.

AIBECS.OCIM2 Module

The OCIM2 module is used to load the OCIM2 matrices and grid for use in AIBECS.

Tip

To load the default OCIM2 matrix and grid, do

julia> grd, T = OCIM2.load()

But you can also load the other matrices by specifying which version you want, e.g.,

julia> grd, T = OCIM2.load(version="OCIM2_KiHIGH_noHe")

See DeVries and Holzer (2019) for more details

Note

The files, that are downloaded from a public and persistent URL on Zenodo, were created with the code available at https://github.com/briochemc/OceanCirculations.

AIBECS.OCIM2.load Method

load

Returns the grid, the transport matrix, and the He fluxes.

Tip

To load the default OCIM2 matrix and grid, do

julia> grd, T = OCIM2.load()

But you can also load the other matrices by specifying which version you want, e.g.,

julia> grd, T = OCIM2.load(version="KiHIGH_noHe")

Add the HeFluxes=true keyword argument if you want the OCIM-produced He fields with ³He and ⁴He as 3rd and 4th arguments. (3rd and 4th output are returned as nothing for "noHe" versions).

See DeVries and Holzer (2019) for more details.

To see all available versions, do

julia> OCIM2.VERSIONS

AIBECS.OCIM2_48L Module

The OCIM2_48L module is used to load the OCIM2-48L matrices and grid for use in AIBECS.

Tip

To load the default OCIM2_48L matrix and grid, do

julia> grd, T = OCIM2_48L.load()

But you can also load the other matrices by specifying which version you want, e.g.,

julia> grd, T = OCIM2_48L.load(version="OCIM2_48L_KiHIGH_noHe")

See DeVries and Holzer (2019) for more details

Note

The files, that are downloaded from a public and persistent URL on Zenodo, were created with the code available at https://github.com/briochemc/OceanCirculations.

AIBECS.OCIM2_48L.load Method

grd, T = load()

Returns the grid and the transport matrix of the OCIM2_48L model.

See DeVries and Holzer (2019) and Holzer et al. (2021) for more details.

Requires using MAT, NCDatasets so that the AIBECSOCIM2_48LExt extension is activated.

AIBECS.OCCA Module

The OCCA module is used to load the OCCA matrix and grid for use in AIBECS.

Tip

To load the OCCA matrix and grid, do

julia> grd, T = OCCA.load()

See Forget (2010) for more details

Note

The files, that are downloaded from a public and persistent URL on Zenodo, were created with the code available at https://github.com/briochemc/OceanCirculations.

Mass conservation

The raw OCCA matrix is not mass-conserving, which can cause issues. To avoid these issues, AIBECS's load() call will enforce mass conservation by default. This is done by subtracting Diagonal((Tᵀ v) ./ v) from T. Pass conservemass = false to recover the raw matrix; in that case tracer budgets will drift unless the geological restoring is modeled explicitly.

AIBECS.OCCA.load Method

grd, T = load()

Returns the grid and the transport matrix.

Requires using JLD2 so that the AIBECSJLD2Ext extension is activated.

See Forget (2010) for more details.

AIBECS.CirculationGeneration Module

CirculationGeneration

Tools for assembling simple advection/diffusion transport operators from sequences of inter-box flow rates, with T_advection and T_diffusion as the two building blocks. Used by the pedagogical circulation modules (TwoBoxModel, Primeau_2x2x2, Archer_etal_2000, Haine_and_Hall_2025).

AIBECS.CirculationGeneration.T_advection Method

T_advection(ϕ, sequence, v3D, nb)

Returns the sparse matrix transport operator, T.

T is such that it gives the flux divergence due to the volumetric flow rate ϕ (in m³ s⁻¹) going through all the boxes in sequence.

AIBECS.CirculationGeneration.T_advection Method

T_advection(ϕ, orig, dest, v3D, nb)

Returns the sparse matrix transport operator, T.

T is such that it gives the flux divergence due to the volumetric flow rate ϕ (in m³ s⁻¹) from box orig to box dest.

AIBECS.CirculationGeneration.T_diffusion Method

T_diffusion(ν, i, j, v3D, nb)

Returns the sparse matrix transport operator, T.

T is such that it gives the flux divergence due to the volumetric mixing rate ν (in m³ s⁻¹) between boxes i and j.

Sinking particles

AIBECS.transportoperator Function

transportoperator(grd, w)

Returns the transportoperator for the given settling velocity w.

The settling velocity can be provided as either a scalar (e.g., w = 100.0 in units of meters per second) or as a function of depth (e.g., w(z) = 2z + 1).

Examples

Create the particle flux divergence with settling velocity of 100m/s

julia> T = transportoperator(grd, 100.0)

Or with settling velocity function w(z) = 2z + 1

julia> T = transportoperator(grd, z -> 2z + 1)

By default, the seafloor flux is set to zero, so that all the particles that reach it are remineralized there. You can let particles go through by setting fsedremin=0.0, via, e.g.,

julia> T = transportoperator(grd, z -> 2z + 1; fsedremin=0.0)

For finer control and advanced use, see the particle-flux divergence operator function, PFDO.

AIBECS.PFDO Function

PFDO(grd; w_top)

Builds the particle-flux-divergence operator PFDO for a given particle sinking speed (w_top).

Schematic of a grid cell:

     top  ┌─────────────────────────────────┐   ┬
          │      ↓ w_top          ↓ Φ_top   │   │
          │ (settling velovity)    (flux)   │   │
          │                                 │   │
          │            x                    │   δz
          │      (particle conc.)           │   │
          │                                 │   │
          │                                 │   │
  bottom  └─────────────────────────────────┘   ┴

• δz is the height of grid cell [m]

• w_top is the particle sinking speed at the top of the grid cell [m s⁻¹]

• Φ_top is the flux at the top of the grid cell [mol m⁻² s⁻¹]

• x is the particle concentration of the cell [mol m⁻³]

The PFDO is defined by

PFDO * x = δΦ/δz ≈ dΦ/dz,

i.e., so that applied to x it approximates the flux divergence of x, dΦ/δz. It is calculated as the matrix product of DIVO and FATO:

PFDO = DIVO * FATO.

where the divergence operator DIVO, is defined by (z increasing downward)

DIVO * ϕ_top = 1/δz * (ϕ_top[below] - ϕ_top) = δϕ/δz ≈ dΦ/δz.

and FATO, the flux-at-top operator, is defined by

FATO * x = w_top * x[above] ≈ ϕ_top.

Example

julia> PFDO(grd, w_top=1.0) # 1.0 m/s (SI units assumed)

PFDO(w, DIVO, Iabove)

Returns DIVO * FATO(w, Iabove)

This function is useful to avoid reconstructing DIVO and Iabove every time. It should allow for faster runs.

PFDO(grd, δz, w_top, w_bot, frac_seafloor, cumfrac_seafloor, fsedremin, Iabove)

Returns the particle-flux-divergence operator for a given sinking speed as a function of depth.

This is a slightly different construction where I take in top and bottom settling velocities, and where the bottom one can be modified to further allow a fraction of particles to sink through (buried into) the sea floor.

Below is a detailed explanation of how this function computes the particle-flux divergence. Take these 3 boxes on top of each other:

┌──────────────────┐
│ water            │
│←----c----→       │
├───────────┲━━━━━━┥ ←- Φ_top
│           ┃⣿⣿⣿⣿⣿⣿│
│           ┃⣿⣿⣿⣿⣿⣿│
│←-d-→ ←-a-→┃⣿⣿⣿⣿⣿⣿│
├─────┲━━━━━┹──────┤ ←- Φ_bot
│     ┃←----b-----→│
│     ┃⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿│
│     ┃⣿⣿⣿ land ⣿⣿⣿│
│     ┃⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿│
└━━━━━┹────────────┘

A part of each of these boxes is water, while the rest is land. For the middle box, we will use the fractional area a = frac_seafloor of seafloor, and the cumulated fractional area b = cumfrac_seafloor of seafloor, to determine the flux of particles at the top ϕ_top and bottom ϕ_bot.

From the top box, only the particles in the water can enter the middle box. So the flux at the top of the middle box, ϕ_top, is proportional to c = 1 - Iabove * b. At the bottom of the middle box, we have to take care of the case of particles hitting the sediments and potentially buried there, or kept in the water. For the part going through the area d, the flux is proportional to d = 1 - b. For the part going through a (hitting the sediments), the flux is proportional to a * (1 - f) where f is the fraction of particles forcibly kept in the water. (So if f = 0, the particles appear to move out of the middle box, but are not carried into the bottom one, and thus are removed from the system / buried.)

AIBECS.DIVO Function

DIVO(grd)

Build the DIVO operator such that

DIVO * ϕ_top = 1/δz * (ϕ_top - ϕ_top[below]) ≈ dΦ/δz.

AIBECS.FATO Function

FATO(w_top, Iabove)

Build the FATO operator for a particle sinking speed w_top

(w_top is the sinking speed at the top of each grid cell.)

The FATO operator is defined by

FATO * x = w_top * x(above) ≈ ϕ_top.

Diagnostics

AIBECS.stencil Function

stencil(T, grd)
stencil(grd, T)

Unique relative neighbour offsets (CartesianIndexes) used by the transport operator T on grid grd.

Useful for inspecting which directions a circulation matrix actually couples boxes along.

grd, T = OCCA.load()
stencil(T, grd)   # e.g. [CI(0,0,0), CI(0,1,0), CI(0,-1,0), ...]

AIBECS.directional_transport Function

directional_transport(T, grd, dir)

Sub-operator of T that contains only the transport along the relative direction dir (a CartesianIndex from stencil).

Useful for diagnosing how much of a circulation acts along, say, the vertical versus the zonal axis.

AIBECS.directional_transports Function

directional_transports(T, grd)

Dictionary mapping human-readable direction labels ("East", "AboveNorth", …) to the corresponding directional_transport sub-operators of T.

grd, T = OCCA.load()
parts = directional_transports(T, grd)
parts["East"]                 # zonal eastward component

AIBECS.smooth_operator Function

smooth_operator(grd, T; σs=(1.0, 1.0, 0.25))

return a matrix of the same size and sparsity as T that smoothes data using a Gaussian Kernel for values, but conserving mass.

This matrix can also likely be used as a covariance matrix for observations in a Bayesian framework.

Plot recipes

AIBECS.plotdepthprofile Function

plotdepthprofile(x, grd; lon, lat, kwargs...)

Plot the depth profile of tracer x at the given lon and lat. Requires using Plots.

AIBECS.plotdepthprofile! Function

plotdepthprofile!(plt, x, grd; lon, lat, kwargs...)

In-place version of plotdepthprofile. Requires using Plots.

AIBECS.plothorizontalmean Function

plothorizontalmean(x, grd; kwargs...)

Plot the area-weighted horizontal mean of tracer x as a depth profile. Aliased as plothorizontalaverage. Requires using Plots.

AIBECS.plothorizontalmean! Function

plothorizontalmean!(plt, x, grd; kwargs...)

In-place version of plothorizontalmean. Requires using Plots.

AIBECS.plothorizontalslice Function

plothorizontalslice(x, grd; depth, kwargs...)

Plot a horizontal (lon-lat) slice of tracer x at the given depth. Requires using Plots.

AIBECS.plothorizontalslice! Function

plothorizontalslice!(plt, x, grd; depth, kwargs...)

In-place version of plothorizontalslice. Requires using Plots.

AIBECS.plotmeridionalmean Function

plotmeridionalmean(x, grd; kwargs...)

Plot the meridionally averaged tracer x as a lon-depth map. Aliased as plotmeridionalaverage. Requires using Plots.

AIBECS.plotmeridionalslice Function

plotmeridionalslice(x, grd; lon, kwargs...)

Plot a meridional (lat-depth) slice of tracer x at longitude lon. Requires using Plots.

AIBECS.plotmeridionalslice! Function

plotmeridionalslice!(plt, x, grd; lon, kwargs...)

In-place version of plotmeridionalslice. Requires using Plots.

AIBECS.plotstencil Function

plotstencil(T, grd; kwargs...)

Visualize the stencil of transport operator T on grid grd. Requires using Plots.

AIBECS.plotstencil! Function

plotstencil!(plt, T, grd; kwargs...)

In-place version of plotstencil. Requires using Plots.

AIBECS.plottransect Function

plottransect(x, grd, ct; kwargs...)

Plot tracer x along the geographic transect ct. Requires using Plots.

AIBECS.plotverticalmean Function

plotverticalmean(x, grd; kwargs...)

Plot the depth-weighted mean of tracer x as a horizontal map. Aliased as plotverticalaverage. Requires using Plots.

AIBECS.plotzonalmean Function

plotzonalmean(x, grd; kwargs...)

Plot the zonally averaged tracer x as a lat-depth map. Aliased as plotzonalaverage. Requires using Plots.

AIBECS.plotzonalmean! Function

plotzonalmean!(plt, x, grd; kwargs...)

In-place version of plotzonalmean. Requires using Plots.

AIBECS.plotzonalslice Function

plotzonalslice(x, grd; lat, kwargs...)

Plot a zonal (lon-depth) slice of tracer x at latitude lat. Requires using Plots.

AIBECS.plot∫dx Function

plot∫dx(x, grd; kwargs...)

Plot the zonal integral of tracer x as a lat-depth map. Aliased as plotzonalintegral. Requires using Plots.

AIBECS.plot∫dxdy Function

plot∫dxdy(x, grd; kwargs...)

Plot the horizontally integrated tracer x as a depth profile. Aliased as plothorizontalintegral. Requires using Plots.

AIBECS.plot∫dxdy! Function

plot∫dxdy!(plt, x, grd; kwargs...)

In-place version of plot∫dxdy. Requires using Plots.

AIBECS.plot∫dy Function

plot∫dy(x, grd; kwargs...)

Plot the meridional integral of tracer x as a lon-depth map. Aliased as plotmeridionalintegral. Requires using Plots.

AIBECS.plot∫dz Function

plot∫dz(x, grd; kwargs...)

Plot the depth-integral of tracer x as a horizontal map. Aliased as plotverticalintegral. Requires using Plots.

AIBECS.plot∫dz! Function

plot∫dz!(plt, x, grd; kwargs...)

In-place version of plot∫dz. Requires using Plots.

AIBECS.ratioatstation Function

ratioatstation(x, y, grd, station; kwargs...)

Scatter plot of x against y for the column at station, useful for tracer-vs-tracer property plots. Requires using Plots.

AIBECS.ratioatstation! Function

ratioatstation!(plt, x, y, grd, station; kwargs...)

In-place version of ratioatstation. Requires using Plots.

AIBECS.surfacemap Function

surfacemap(x, grd; kwargs...)

Plot a horizontal map of tracer x at the surface. Requires using Plots.

AIBECS.surfacemap! Function

surfacemap!(plt, x, grd; kwargs...)

In-place version of surfacemap. Requires using Plots.

Datasets

AIBECS.AeolianSources Module

This AeolianSources module is to create aeolian sources for use in AIBECS-generated models.

AIBECS.AeolianSources.load Function

load()

Returns the 2D aeolian deposition fields from Chien et al. (2016).

You can specify a different dataset via, e.g.,

julia> AeolianSources.load("Kok")

At this stage, only two datasets are available:

• "Chien" (default) for different dust types (fires, biofuels, dust, ...)

• "Kok" for dust from different regions of origin

AIBECS.Rivers Module

Rivers

Annual-mean discharge locations and volumetric flow rates of the 200 largest rivers, extracted from the Dai and Trenberth global river-flow dataset (runoff-table-A.txt and runoff-table2-top50r.txt on UCAR's RDA, used with permission from Aiguo Dai). Use Rivers.load to obtain the list and regrid to project it onto an OceanGrid.

See Rivers.citation for citation strings.

AIBECS.Rivers.River Type

River(name, lon, lat, VFR)

A river entry: name, mouth lon/lat (degrees), and volumetric flow rate VFR (e.g. in km^3/yr or m^3/s via Unitful).

AIBECS.Rivers.citation Method

Rivers.citation()

Multi-line citation string for the Dai and Trenberth river-discharge dataset (plus its updates). BibTeX entries are stored in CITATION.bib under the keys Dai_2017, Dai_Trenberth_2002, Dai_etal_2009, and Dai_2016.

AIBECS.Rivers.load Function

Rivers.load(unit=m^3/s)

Return the 200 largest rivers as a Vector{River}, with volumetric flow rates converted to unit.

OceanGrids.regrid Method

regrid(R::Vector{River}, grd)

Bin a list of rivers R onto the surface boxes of grd, returning a vector of volumetric flow rates per wet box (zero for boxes with no river).

AIBECS.GroundWaters Module

GroundWaters

Coastal fresh-groundwater discharge dataset of Luijendijk et al. (2020), exposed as GroundWaterSource entries that can be regridded onto an OceanGrid. Loading the dataset requires the Shapefile and DataFrames packages so that the AIBECSShapefileExt extension is activated.

AIBECS.GroundWaters.GroundWaterSource Type

GroundWaterSource(lon, lat, VFR)

A point-source groundwater discharge: longitude, latitude, and volumetric flow rate VFR (Unitful quantity).

AIBECS.GroundWaters.load Method

gws = load()

Returns the coastal groundwater discharge dataset.

Requires using Shapefile, DataFrames so that the AIBECSShapefileExt extension is activated.

OceanGrids.regrid Method

regrid(gws, grd)

Regrids and bins the groundwater discharge values into grd surface boxes.

AIBECS.ETOPO Module

ETOPO

Access to the 1-arc-minute ETOPO1 global relief model (Amante and Eakins, 2009), either the bedrock (:bedrock) or ice-surface (:ice) variant. Provides ETOPO.load plus the bin-on-grid helpers fractiontopo and roughnesstopo for mapping subgrid topography onto an OceanGrid.

Loading the dataset requires Distances and NCDatasets so that the AIBECSETOPOExt extension is activated.

AIBECS.ETOPO.bintopo Method

bintopo(Z, lat, lon, grd)

Bins the depths Z into at locations (lat,lon) onto grd. To be used with the ETOPO dataset to dertermine the amount of subgrid topography in each box area

AIBECS.ETOPO.fractiontopo Function

fractiontopo(grd)

Returns the fraction of subgrid sediments in each wet box of grd.

Requires using Distances, NCDatasets so that the AIBECSETOPOExt extension is activated. Note that (i) points located above sea level are counted in the top layer, and (ii) points located in a dry box or below are counted in the deepest wet box.

AIBECS.ETOPO.load Method

Z, lats, lons = ETOPO.load()
Z, lats, lons = ETOPO.load(:bedrock)
Z, lats, lons = ETOPO.load(:ice)

Returns the fine resolution topography from ETOPO.

Requires using Distances, NCDatasets so that the AIBECSETOPOExt extension is activated.

AIBECS.ETOPO.roughnesstopo Function

roughnesstopo(grd)

Estimates topographic roughness from ETOPO.

Requires using Distances, NCDatasets so that the AIBECSETOPOExt extension is activated.

AIBECS.AO Module

AO

Helpers to download the AWESOME-OCIM (AO) MATLAB toolbox of John et al. (2020) as a zipped GitHub archive and unpack it locally via DataDeps. The unpacked tree contains MATLAB code, OCIM1 transport matrices, and auxiliary GEOTRACES/WOA datasets bundled by the AO authors.

Note: the upstream zip is served by GitHub directly (there is no official mirror), so downloads are subject to GitHub availability — see download_and_unpack.

AIBECS.AO.download_and_unpack Method

download_and_unpack

Downloads and unpacks the AO zip file from the MTEL website. Returns the local path of where the AO is installed for you.

Solver

AIBECS.CTKAlg Type

CTKAlg{L} <: SciMLBase.AbstractSteadyStateAlgorithm

Marker algorithm that selects AIBECS's customised Newton–Chord–Shamanskii solver when passed to SciMLBase.solve(::SteadyStateProblem, ::CTKAlg). The linsolve field is forwarded to LinearSolve.LinearProblem, so any SciMLLinearSolveAlgorithm (UMFPACK, KLU, MKL Pardiso, …) can be plugged in to swap the inner linear solver while keeping the Shamanskii nonlinear logic.

prob = SteadyStateProblem(fun, x₀, p)
sol  = solve(prob, CTKAlg(); maxItNewton = 50)
sol  = solve(prob, CTKAlg(linsolve = KLUFactorization()))

The underlying engine is AIBECS.NewtonChordShamanskii (with Armijo line search and lazy Jacobian refresh, after Kelley 2003).

The bound solve(::SteadyStateProblem, ::CTKAlg) method is documented under AIBECS.CTKAlg. For the upstream SteadyStateProblem constructor and generic solve entry point, see the SciMLBase documentation.