Radiocarbon

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
grd, T_OCIM2 = OCIM2.load()
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[1],
                          τ = 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)
50 100 150 200 250 300 350 -60 -30 0 30 60 Longitude (°) Latitude (°) 750 1000 1250 1500 1750 yr

look at a zonal average using the zonalaverage plot recipe

zonalaverage(C14age, grd; color=:viridis)
-60 -30 0 30 60 0 1000 2000 3000 4000 5000 Latitude (°) Depth (m)