AIBECS functions

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

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.

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.

julia> grd, T = OCCA.load()

load

Returns grd and T (in that order).

load

Returns grd and T (in that order).

load

Returns grd and T (in that order).

plothorizontalslice(x, grd; depth)

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

surfacemap(x, grd)

Plots a surface heatmap of tracer x.

plotverticalintegral(x, grd, mask=1)

Plots the vertical integral of tracer x.

plotverticalaverage(x, grd, mask=1)

Plots the vertical average of tracer x.

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)

plotmeridionalslice(x, grd; lon)

Plots a meridional slice of tracer x at longitude lon.

plotzonalmean(x, grd; mask=1)

Plots a zonal average of tracer x.

plotzonalintegral(x, grd; mask=1)

Plots a zonal integral of tracer x.

plotzonalslice(x, grd; lat)

Plots a zonal slice of tracer x at latitude lat.

plotmeridionalmean(x, grd; mask=1)

Plots a meridional average of tracer x.

plotmeridionalintegral(x, grd; mask=1)

Plots a meridional integral of tracer x.

horizontalintegral(x, grd; mask=1)

Plots a horizontal integral of tracer x.

plothorizontalaverage(x, grd; mask=1)

Plots a horizontal average of tracer x.

plotdepthprofile(x, grd; lonlat)

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

plottransect(x, grd; ct=ct)

Plots the transect of tracer x along transect ct.

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)

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.

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.

latex(p)

Returns a LaTeX-formatted table of the parameters.

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.

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!

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

Solves prob using the modified C.T.Kelley Shamanskii algorithm.

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.

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

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.