AIBECS interface (functions and types)

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.

AgedJacobianFactors

Type containing the Jacobian Factors and the age of the Jacobian. This allows for the Shamanskii method to not update the Jacobian at each iterate.

DIVO(grd)

Build the DIVO operator such that

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

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.

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.)

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; 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)

flattenable_symbols(p)

Returns the flattenable symbols in p.

The flattenable symbols are those symbols that are kepth when using p as a vector, e.g., when doing vec(p). (Useful when passing parameters to an optimization routine expeting a vector of optimizable parameters.)

Can also be used directly on the type of p (because the flattenable symbols of p::T are contained in the type T).

flattenable_values(p::T) where {T <: AbstractParameters}

Returns a vector of the flattenable values of p.

Note that vec(p) is different from values(p).

latex(p)

Returns a LaTeX-formatted table of the parameters.

localderivative(G, x, p)
localderivative(Gᵢ, xs, i, p)
localderivative(Gᵢ, dx, xs, i, p)

Returns the "local" derivative of G (or Gᵢ), i.e., equivalent to the vector

∇ₓG(x,p) * ones(size(x))

but using ForwardDiff's Jacobian instead.

minimap(grd; central_longitude=200°)

Plots a surface map of grd with no ticks, labels, or colorbar.

The goal of this function is to provide a way to plot both a minimap and a cruise track when plotting a transect:

plottransect(dummy, grd; ct=sort(ct))
plot!(inset_subplots=bbox(0.73,0.73,0.15,0.15), subplot=2)
minimap!(grd; subplot=2)
plotcruisetrack!(sort(ct), subplot=2)

minimap(grd; central_longitude=200°)

Plots a surface map of grd with no ticks, labels, or colorbar.

The goal of this function is to provide a way to plot both a minimap and a cruise track when plotting a transect:

plottransect(dummy, grd; ct=sort(ct))
plot!(inset_subplots=bbox(0.73,0.73,0.15,0.15), subplot=2)
minimap!(grd; subplot=2)
plotcruisetrack!(sort(ct), subplot=2)

minimap(grd; central_longitude=200°)

Plots a surface map of grd with no ticks, labels, or colorbar.

The goal of this function is to provide a way to plot both a minimap and a cruise track when plotting a transect:

plottransect(dummy, grd; ct=sort(ct))
plot!(inset_subplots=bbox(0.73,0.73,0.15,0.15), subplot=2)
minimap!(grd; subplot=2)
plotcruisetrack!(sort(ct), subplot=2)

mismatch(x, xobs, σ²xobs, v)

Volume-weighted mismatch of modelled tracer x against observed mean, xobs, given observed variance, σ²xobs, and volumes v.

onlykeep(x::MetadataArray, idx)

Returns x[idx] but also applies the index to the metadata that is originally of the same length as x.

p2λ(p::AbstractParameters)

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

plotdepthprofile(x, grd; lonlat)

Plots the profile of tracer x interpolated at lonlat=(x,y) coordinates.

plothorizontalaverage(x, grd; mask=1)

Plots a horizontal average of tracer x.

plothorizontalslice(x, grd; depth)

Plots a heatmap of the horizontal slice of tracer x at depth depth.

plotmeridionalmean(x, grd; mask=1)

Plots a meridional average of tracer x.

plotmeridionalslice(x, grd; lon)

Plots a meridional slice of tracer x at longitude lon.

PlotParameter(p::AbstractParameters, s)

Plots the PDF of parameter p with symbol s

PlotParameter(p::AbstractParameters, s)

Plots the PDF of parameter p with symbol s

PlotParameter(p::AbstractParameters, s)

Plots the PDF of parameter p with symbol s

PlotParameters(p::AbstractParameters)

Plots the PDF of all the flattenable parameters in p.

PlotParameters(p::AbstractParameters)

Plots the PDF of all the flattenable parameters in p.

PlotParameters(p::AbstractParameters)

Plots the PDF of all the flattenable parameters in p.

PlotStencil(st)

Plots the stencil st.

PlotStencil(st)

Plots the stencil st.

PlotStencil(st)

Plots the stencil st.

plottransect(x, grd; ct=ct)

Plots the transect of tracer x along transect ct.

plotverticalaverage(x, grd, mask=1)

Plots the vertical average of tracer x.

plotzonalmean(x, grd; mask=1)

Plots a zonal average of tracer x.

plotzonalslice(x, grd; lat)

Plots a zonal slice of tracer x at latitude lat.

plotzonalintegral(x, grd; mask=1)

Plots a zonal integral of tracer x.

horizontalintegral(x, grd; mask=1)

Plots a horizontal integral of tracer x.

plotmeridionalintegral(x, grd; mask=1)

Plots a meridional integral of tracer x.

plotverticalintegral(x, grd, mask=1)

Plots the vertical integral of tracer x.

RatioAtStation(x, y, grd, station, depthlims=(0,Inf))

Plots a meridional transect of tracer x along cruise track ct.

The keyword argument zlims=(ztop, zbottom) can be provided if you only want to only plot for depths z ∈ (ztop, zbottom). (z is positive downwards in this case)

RatioAtStation(x, y, grd, station, depthlims=(0,Inf))

Plots a meridional transect of tracer x along cruise track ct.

The keyword argument zlims=(ztop, zbottom) can be provided if you only want to only plot for depths z ∈ (ztop, zbottom). (z is positive downwards in this case)

RatioAtStation(x, y, grd, station, depthlims=(0,Inf))

Plots a meridional transect of tracer x along cruise track ct.

The keyword argument zlims=(ztop, zbottom) can be provided if you only want to only plot for depths z ∈ (ztop, zbottom). (z is positive downwards in this case)

reconstruct(T, v)

Reconstructs the parameter of type T from flattenable vector v.

Note that T must be the generic type (or type wrapper if you will). I.e., T must not include the parametric scalar type, as in, use T == Params instead of T == Params{Float64}.

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.

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!

surfacemap(x, grd)

Plots a surface heatmap of tracer x.

symbols(p)

Returns the symbols in p.

Can also be used directly on the type of p (because the symbols of p::T are contained in the type T).

table(p)

Returns a DataFrame (a table) of p. Useful for printing and saving into an actual text/latex table.

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.

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

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

Note that the instance of your parameters type T is required here because it contains information on non-optimizable parameters and priors of optimizable parameters

∇mismatch(x, xobs, σ²xobs, v)

Adjoint of the gradient of mismatch(x, xobs, σ²xobs, v).

*(x::MetadataArray, q::Quantity)

Preserve metadata (and history of unit conversions) when multiplying `x` by a quantity.

+(p::T, v::Vector) where {T <: AbstractParameters}

Adds the flattened vector v to p.

Warning: This method for + is implemented only for differentiation using dual and hyperdual numbers. If you want to change the values of p, you should do so explicitly rather than use this + method.

\(Jf::JacobianFactors, y)

Dispatches backslash to work with all JacobianFactors subtypes.

getindex(p::T, i) where {T <: AbstractParameters}

Returns the i-th element of vec(p).

This is not efficient and only used for testing the derivatives with ForwardDiff.

length(p::AbstractParameter)

Returns the length of the flattened/optimizable vector of p.

May be different from the number of parameters. Can also be used directly on the type of p.

size(p::AbstractParameter)

Returns the size of the flattened/optimizable vector of p.

May be different from the number of parameters. Can also be used directly on the type of p.

values(p::T) where {T <: AbstractParameters}

Returns a vector of all the values of p.

Note that values(p) is different from vec(p).

vec(p::T) where {T <: AbstractParameters}

Returns a SI-unit-converted vector of flattenable values of p.

Note that vec(p) ≠ flattenable_values(p) if p has units.

bijector(T::AbstractParameters [, k::Symbol])

Returns the function for the change of variables of the parameters.

The function is a bijection from the supports/domains of the priors to ℝⁿ, from the Bijectors.jl package. You can specify the parameter symbol to get the bijector of that parameter.

solve(prob::SciMLBase.AbstractSteadyStateProblem,
      alg::CTKAlg;
      nrm=norm,
      τstop=1e12*365*24*60*60,
      preprint="",
      maxItNewton=50)

Solves prob using an AIBECS-customized version of C.T.Kelley Shamanskii-method algorithm.

unpack(p <: AbstractParameters, s)

Unpacks the parameter s from p.

Note this is specialized and will convert the parameter value to SI units.

upreferred(x::MetadataArray)

Converts x.parent to SI unit but keeps x.metadata for safekeeping.

ustrip(x::MetadataArray)

Strips unit from x.parent but stores it in x.metadata for safekeeping.

load

Returns grd and T (in that order).

load

Returns grd and T (in that order).

load

Returns grd and T (in that order).

load

Returns the grid and the transport matrix (in that order).

Usage

julia> grd, T = OCIM0.load()

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

load

Returns the grid and the transport matrix (in that order).

Usage

julia> grd, T = OCIM1.load()

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

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

julia> grd, T = OCIM2.load()

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

load

Returns the grid, the transport matrix, and the He fluxes (in that order).

julia> grd, T = OCIM2.load()

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

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

julia> grd, T = OCIM2_48L.load()

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

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.

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

julia> grd, T = OCCA.load()

load

Returns the grid and the transport matrix.

julia> grd, T = OCCA.load()

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.

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

load()

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

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

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

Z, lats, lons = load()

Returns the fine resolution topography from ETOPO.

regrid(gws, grd)

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

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

fractiontopo(grd)

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

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.

Z, lats, lons = ETOPO.load()

Returns the fine resolution topography from ETOPO.

It is the same as

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

If you want the ice elevation, use

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

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.

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.

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.