Ideal age

The tracer equation for the ideal age is

\[\left(\partial_t + \mathbf{T}\right) \boldsymbol{a} = 1 - \frac{\boldsymbol{a}}{τ} \, (\boldsymbol{z} \le z_0),\]

where the sink term on the right clamps the age to $0$ at the surface (where $\boldsymbol{z} \le z_0$). The smaller the timescale $\tau$, the quicker $\boldsymbol{a}$ is restored to $0$ at the surface.

AIBECS can interpret tracer equations as long as you arrange them under the generic form:

\[\big(\partial_t + \mathbf{T}(\boldsymbol{p}) \big) \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}),\]

where $\mathbf{T}(\boldsymbol{p})$ is the transport, $\boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p})$ is the net local sources and sinks, and $\boldsymbol{p}$ is the vector of model parameters. We will then use the AIBECS to simulate the ideal age by finding the steady-state of the system, i.e., the solution of

\[\partial_t \boldsymbol{x} = \boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}) - \mathbf{T}(\boldsymbol{p}) \, \boldsymbol{x} = 0.\]

In this tutorial, we will simulate the ideal age by

1. defining functions for T(p) and G(x,p),
2. defining the parameters p,
3. generating the state function F(x,p) and solving the associated steady-state problem,
4. and finally making a plot of our simulated ideal age.

We start by telling Julia that we want to use the AIBECS package and the OCIM2 circulation (the Ocean Circulation Inverse Model^[1]).

using AIBECS
grd, TOCIM2_48L = OCIM2_48L.load()

(, sparse([1, 2, 10385, 10443, 10444, 20827, 1, 2, 3, 62  …  405384, 405385, 405386, 404655, 405279, 405283, 405284, 405383, 405385, 405386], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2  …  405385, 405385, 405385, 405386, 405386, 405386, 405386, 405386, 405386, 405386], [0.0005175223035904686, -2.2822124779127418e-8, -3.0194836257899404e-7, -0.0005095966035957547, 1.0484064034372429e-8, 1.6632929441401007e-8, -5.357525815913356e-8, 0.0005298878550923177, 3.5862540303001864e-9, -7.257834319066547e-8  …  -1.630240831994384e-8, 4.1471866886309256e-8, -1.3121456023816401e-8, 2.585754349113357e-11, 1.5323197420723728e-9, 2.8868819287010998e-9, -7.724482133397593e-9, -1.3957072091171105e-8, -1.1289022605819024e-8, 2.8541808520184807e-8], 405386, 405386))

grd is an OceanGrid object containing information about the 3D grid of the OCIM2 circulation and TOCIM2_48L is the transport matrix representing advection and diffusion.

The local sources and sinks for the age take the form

function G(x,p)
    @unpack τ, z₀ = p
    return @. 1 - x / τ * (z ≤ z₀)
end

G (generic function with 1 method)

as per the tracer equation. The @unpack line unpacks the parameters τ and z₀. The return line returns the net sources and sinks. (The @. "macro" tells Julia that the operations apply to every element.)

We can define the vector z of depths with depthvec.

z = depthvec(grd)

405386-element Vector{Float64}:
    4.934540873153541
    4.934540873153541
    4.934540873153541
    4.934540873153541
    4.934540873153541
    4.934540873153541
    4.934540873153541
    4.934540873153541
    4.934540873153541
    4.934540873153541
    ⋮
 5581.917201508207
 5581.917201508207
 5581.917201508207
 5581.917201508207
 5581.917201508207
 5581.917201508207
 5581.917201508207
 5581.917201508207
 5581.917201508207

Now we must construct a type for p the parameters. This type must contain our parameters τ and z₀.

struct IdealAgeParameters{U} <: AbstractParameters{U}
    τ::U
    z₀::U
end

The type is now ready for us to generate an instance of the parameter p. Let's use τ = 1.0 (s) and z₀ the minimum depth of the model.

p = IdealAgeParameters(1.0, 30.0)

 Row │ Symbol  Value
     │ Symbol  Float64
─────┼─────────────────
   1 │ τ           1.0
   2 │ z₀         30.0

We now use the AIBECS to generate the state function $\boldsymbol{F}$ (and its Jacobian) via

F = AIBECSFunction(TOCIM2_48L, G)

(::SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, AIBECS.var"#f#58"{Tuple{AIBECS.var"#49#50"{SparseMatrixCSC{Float64, Int64}}}, Vector{Int64}, AIBECS.var"#G#56"{Tuple{typeof(Main.G)}, AIBECS.var"#tracers#54"{Int64, Int64}}, AIBECS.var"#tracer#55"{Int64, Int64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, AIBECS.var"#jac#63"{AIBECS.var"#T#60"{Tuple{AIBECS.var"#49#50"{SparseMatrixCSC{Float64, Int64}}}, Int64, Vector{Int64}}, AIBECS.var"#∇ₓG#59"{Tuple{typeof(Main.G)}, Int64, Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}) (generic function with 1 method)

Now that F and p are defined, we are going to solve for the steady-state. But first, we must create a SteadyStateProblem object that contains F, p, and an initial guess x_init for the age. (SteadyStateProblem comes from SciMLBase.)

Let's make a vector of 0's for our initial guess.

nb = sum(iswet(grd))  # number of wet boxes
x_init = zeros(nb)    # Start with age = 0 everywhere

405386-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

Now we can create our SteadyStateProblem instance

prob = SteadyStateProblem(F, x_init, p)

SteadyStateProblem with uType Vector{Float64}. In-place: false
u0: 405386-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

And finally, we can solve this problem, using the AIBECS CTKAlg() algorithm,

age = solve(prob, CTKAlg())

retcode: Success
u: 405386-element Vector{Float64}:
 1.0032939873142694
 1.0050746742544991
 1.0002618470038183
 1.0004491845445038
 1.0003315941268929
 1.0000482918824332
 1.0002351572130561
 1.0003995888227937
 1.0006776625690539
 1.0022588888620927
 ⋮
 3.529870757202301e10
 1.2433539018396751e10
 1.2416113239157198e10
 1.236564076423028e10
 1.2339193695370579e10
 1.2272466964580017e10
 1.230063368641268e10
 1.228925485656132e10
 1.2292310143856483e10

This should take a few seconds.

To conclude this tutorial, let's have a look at the age using AIBECS' plotting recipes and Plots.jl.

using Plots

We first convert the age in years (because the default SI unit we used, i.e., seconds, is a bit small relative to global ocean timescales).

age_in_yrs = age * u"s" .|> u"yr"

405386-element Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(yr,), 𝐓, nothing}}}:
 3.1792467973301814e-8 yr
  3.184889453743311e-8 yr
  3.169638524487978e-8 yr
  3.170232161331989e-8 yr
  3.169859539784055e-8 yr
 3.1689618091440196e-8 yr
  3.169553949644637e-8 yr
  3.170075001973514e-8 yr
  3.170956164502541e-8 yr
  3.175966768265307e-8 yr
                        ⋮
    1118.5485452639937 yr
    393.99507625411155 yr
     393.4428866313407 yr
     391.8435104136651 yr
    391.00545337321523 yr
     388.8910108683809 yr
     389.7835604232476 yr
     389.4229870636968 yr
    389.51980327580304 yr

And we take a horizontal slice at about 2000m.

plothorizontalslice(age_in_yrs, grd, depth=2000u"m", color=:magma)

Or look at the horiontal mean

plothorizontalmean(age_in_yrs, grd)

That's it for this tutorial... Good job!

This page was generated using Literate.jl.