In this tutorial, we will simulate the radiocarbon age using the AIBECS by

1. defining the transport T(p) and the sources and sinks 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 radiocarbon age.
Note

Although this tutorial is self-contained, it involves non-linearitiess and is slightly more complicated than the first tutorial for simulating the ideal age. (So do not hesitate to start with the ideal-age tutorial if you wish.)

The tracer equation for radiocarbon is

$\big(\partial_t + \mathbf{T} \big) \boldsymbol{R} = \frac{\lambda}{h} (\overline{\boldsymbol{R}}_\mathsf{atm} - \boldsymbol{R}) (\boldsymbol{z} ≤ h) - \boldsymbol{R} / \tau.$

where the first term on the right of the equal sign represents the air–sea gas exchange with a piston velocity $λ$ over a depth $h$ and the second term represents the radioactive decay of radiocarbon with timescale $\tau$.

Note

We need not specify the value of the atmospheric radiocarbon concentration because it is not important for determining the age of a water parcel — only the relative concentration $\boldsymbol{R}/\overline{\boldsymbol{R}}_\mathsf{atm}$ matters.

We start by selecting the circulation for Radiocarbon

using AIBECS
T(p) = T_OCIM2
T (generic function with 1 method)

The local sources and sinks are simply given by

function G(R,p)
@unpack λ, h, R̅atm, τ = p
return @. λ / h * (R̅atm - R) * (z ≤ h) - R / τ
end
G (generic function with 1 method)

We can define z via

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

In this tutorial we will specify some units for the parameters. Such features must be imported to be used

import AIBECS: @units, units

We define the parameters using the dedicated API from the AIBECS, including keyword arguments and units this time

@units struct RadiocarbonParameters{U} <: AbstractParameters{U}
λ::U    | u"m/yr"
h::U    | u"m"
τ::U    | u"yr"
R̅atm::U | u"M"
end
units (generic function with 39 methods)

For the air–sea gas exchange, we use a constant piston velocity $\lambda$ of 50m / 10years. And for the radioactive decay we use a timescale $\tau$ of 5730/log(2) years.

p = RadiocarbonParameters(λ = 50u"m"/10u"yr",
h = grd.δdepth,
τ = 5730u"yr"/log(2),
R̅atm = 42.0u"nM")

│ Row │ Symbol │ Value   │ Unit     │
│     │ Symbol │ Float64 │ Unitful… │
├─────┼────────┼─────────┼──────────┤
│ 1   │ λ      │ 5.0     │ m yr^-1  │
│ 2   │ h      │ 36.1351 │ m        │
│ 3   │ τ      │ 8266.64 │ yr       │
│ 4   │ R̅atm   │ 4.2e-8  │ M        │
Note

The parameters are converted to SI units when unpacked. When you specify units for your parameters, you must either supply their values in that unit.

We generate the state function and its Jacobian, generate the corresponding steady-state problem, and solve it, via

F, ∇ₓF = state_function_and_Jacobian(T, G)
x = zeros(length(z)) # an initial guess
prob = SteadyStateProblem(F, ∇ₓF, x, p)
R = solve(prob, CTKAlg()).u
200160-element Array{Float64,1}:
3.5560535238235804e-5
3.556277056874239e-5
3.550107499504647e-5
3.5523387657289385e-5
3.5824950925183595e-5
3.592509725690728e-5
3.581202090505318e-5
3.589200096626715e-5
3.6058151597808294e-5
3.628098578395497e-5
⋮
3.651437090109276e-5
3.660112493330399e-5
3.617455107082283e-5
3.6320739193845454e-5
3.655481631782142e-5
3.67343646405315e-5
3.625512385606775e-5
3.626845589255476e-5
3.6524638793488156e-5

This should take a few seconds on a laptop. Once the radiocarbon concentration is computed, we can convert it into the corresponding age in years, via

@unpack τ, R̅atm = p
C14age = @. log(R̅atm / R) * τ * u"s" |> u"yr"
200160-element Array{Quantity{Float64,𝐓,Unitful.FreeUnits{(yr,),𝐓,nothing}},1}:
1375.8434152993696 yr
1375.32379157963 yr
1389.6775147947856 yr
1384.4835068981781 yr
1314.6029561590803 yr
1291.5263314298802 yr
1317.5871095495904 yr
1299.1455539193744 yr
1260.9660566249188 yr
1210.0366022347087 yr
⋮
1157.0299807672932 yr
1137.4126642078122 yr
1234.3234518831982 yr
1200.9837246614934 yr
1147.8784395735056 yr
1107.3741042683553 yr
1215.9313579502964 yr
1212.8920379889166 yr
1154.7057160369573 yr

and plot it at 700 m using the horizontalslice Plots recipe

using Plots
horizontalslice(C14age, grd, depth=700u"m", color=:viridis)

look at a zonal average using the zonalaverage plot recipe

zonalaverage(C14age, grd; color=:viridis)