# Ideal age

Note

All the AIBECS tutorials and how-to guides are available as Jupyter notebooks. You can execute them online with binder or just view them with nbviewer by clicking on the badges above!

The tracer equation for the ideal age is

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

where the sink term on the right clamps the age to $0$ at the surface (where $\boldsymbol{z} < 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 = OCIM2.load()
(,
[1     ,      1]  =  0.000197784
[2     ,      1]  =  2.34279e-9
[10384 ,      1]  =  -1.95995e-7
[10442 ,      1]  =  -0.000191612
[10443 ,      1]  =  4.80961e-9
[20825 ,      1]  =  -1.83059e-9
[20883 ,      1]  =  5.00768e-9
[1     ,      2]  =  -5.02516e-8
[2     ,      2]  =  0.000187531
⋮
[200160, 200159]  =  -2.19656e-8
[197886, 200160]  =  1.08199e-10
[199766, 200160]  =  6.70981e-9
[199777, 200160]  =  -1.26352e-9
[199778, 200160]  =  -3.39279e-9
[199779, 200160]  =  7.59316e-9
[199790, 200160]  =  -7.41018e-9
[200156, 200160]  =  -3.44106e-8
[200159, 200160]  =  -2.00303e-8
[200160, 200160]  =  5.27945e-8)
Note

If it's your first time, Julia will ask you to download the OCIM2, in which case you should accept (i.e., type y and "return"). Once downloaded, AIBECS will remember where it downloaded the file and it will only load it from your laptop.

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

We define the function T(p) as

T(p) = TOCIM2
T (generic function with 1 method)

(It turns out the circulation T(p) does not effectively depend on p but that's how we must define it anyway, i.e., as a function of p.)

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)
200160-element Array{Float64,1}:
18.0675569520817
18.0675569520817
18.0675569520817
18.0675569520817
18.0675569520817
18.0675569520817
18.0675569520817
18.0675569520817
18.0675569520817
18.0675569520817
⋮
5433.2531421838175
5433.2531421838175
5433.2531421838175
5433.2531421838175
5433.2531421838175
5433.2531421838175
5433.2531421838175
5433.2531421838175
5433.2531421838175

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₀ = 20.0 (m) .

p = IdealAgeParameters(1.0, 20.0)

│ Row │ Symbol │ Value   │
│     │ Symbol │ Float64 │
├─────┼────────┼─────────┤
│ 1   │ τ      │ 1.0     │
│ 2   │ z₀     │ 20.0    │

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

F, ∇ₓF = state_function_and_Jacobian(T, G)
(AIBECS.var"#F#31"{typeof(Main.ex-1_ideal_age.T),typeof(Main.ex-1_ideal_age.G)}(Main.ex-1_ideal_age.T, Main.ex-1_ideal_age.G), AIBECS.var"#∇ₓF#33"{typeof(Main.ex-1_ideal_age.T),AIBECS.var"#∇ₓG#32"{typeof(Main.ex-1_ideal_age.G)}}(Main.ex-1_ideal_age.T, AIBECS.var"#∇ₓG#32"{typeof(Main.ex-1_ideal_age.G)}(Main.ex-1_ideal_age.G)))

(∇ₓF is the Jacobian of the state function $\nabla_{\boldsymbol{x}}\boldsymbol{F}$, calculated automatically using dual numbers.)

Now that F(x,p), and p are defined, we are going to solve for the steady-state. But first, we must create a SteadyStateProblem object that contains F, ∇ₓF, p, and an initial guess x_init for the age. (SteadyStateProblem is specialized from DiffEqBase for AIBECS models.)

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
200160-element Array{Float64,1}:
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, ∇ₓF, x_init, p)
SteadyStateProblem with uType Array{Float64,1}
u0: [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, 0.0]

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

age = solve(prob, CTKAlg())
u: 200160-element Array{Float64,1}:
804.9793529618016
2704.2854677669357
559.4624759776626
277.6810176373408
571.2245067452857
406.14915296475465
416.88383377410986
645.378543317121
411.6472087934807
242.71573068291767
⋮
1.2204098653557419e10
1.1820841832109285e10
1.3783960333917171e10
1.308318053973743e10
1.1999572353055996e10
1.122761721268902e10
1.3402055608211155e10
1.3322655154055176e10
1.2121012215573376e10

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"
200160-element Array{Quantity{Float64,𝐓,Unitful.FreeUnits{(yr,),𝐓,nothing}},1}:
2.5508256425133772e-5 yr
8.5693635376801e-5 yr
1.7728296067434234e-5 yr
8.799180471180977e-6 yr
1.8101012331269982e-5 yr
1.2870090024740621e-5 yr
1.321025153288304e-5 yr
2.0450811953923014e-5 yr
1.3044312900647725e-5 yr
7.691197387726494e-6 yr
⋮
386.72454982499994 yr
374.57987401162586 yr
436.7873454862591 yr
414.5809738299943 yr
380.24350245443236 yr
355.7817201779926 yr
424.6855150014942 yr
422.16946643772576 yr
384.0916994820067 yr

And we take a horizontal slice at about 2000m.

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

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