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.

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

We start by selecting the circulation for Radiocarbon. .) (And this time, we are using the OCCA matrix by Forget ^[1].)

using AIBECS
grd, T_OCCA = OCCA.load()

(, sparse([1, 2, 9551, 9606, 1, 2, 3, 57, 9552, 9607  …  79928, 84632, 84659, 84660, 84661, 79891, 79929, 84633, 84660, 84661], [1, 1, 1, 1, 2, 2, 2, 2, 2, 2  …  84660, 84660, 84660, 84660, 84660, 84661, 84661, 84661, 84661, 84661], [2.2827161967151829e-7, 1.8030621712982388e-10, -2.202759122121374e-7, -1.669794081839543e-8, -8.32821319857397e-8, 3.645587487556812e-7, -2.4345873581518986e-8, 2.346272528650749e-8, -3.294980961099172e-7, 6.798120022902578e-9  …  -2.337773104642457e-8, -1.0853897038905248e-8, -2.4801404742669764e-8, 7.658935248883193e-8, -2.49410804832536e-8, -2.2788300189327118e-9, -3.386139290468414e-8, -2.0539768803799346e-8, -1.8193948479846445e-8, 5.86467452504204e-8], 84661, 84661))

The local sources and sinks are simply given by

function G(R,p)
    @unpack λ, h, Ratm, τ = p
    return @. λ / h * (Ratm - R) * (z ≤ h) - R / τ
end

G (generic function with 1 method)

We can define z via

z = depthvec(grd)

84661-element Vector{Float64}:
   25.0
   25.0
   25.0
   25.0
   25.0
   25.0
   25.0
   25.0
   25.0
   25.0
    ⋮
 5062.25
 5062.25
 5062.25
 5062.25
 5062.25
 5062.25
 5062.25
 5062.25
 5062.25

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"
    Ratm::U | u"M"
end

units (generic function with 46 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),
                          Ratm = 42.0u"nM")

 Row │ Symbol  Value      Unit
     │ Symbol  Float64    FreeUnit…
─────┼──────────────────────────────
   1 │ λ          5.0     m yr⁻¹
   2 │ h         50.0     m
   3 │ τ       8266.64    yr
   4 │ Ratm       4.2e-8  M

We build the state function F and the corresponding steady-state problem (and solve it) via

F = AIBECSFunction(T_OCCA, G)
x = zeros(length(z)) # an initial guess
prob = SteadyStateProblem(F, x, p)
R = solve(prob, CTKAlg()).u

84661-element Vector{Float64}:
 3.7813682936931906e-5
 3.775789970720447e-5
 3.660570790298678e-5
 3.707009029891997e-5
 3.715744797475306e-5
 3.71583992727843e-5
 3.707936439247658e-5
 3.7149283799005414e-5
 3.7232859237789926e-5
 3.735496095490033e-5
 ⋮
 3.638493414383584e-5
 3.6381532923815344e-5
 3.637308898920196e-5
 3.638148729939853e-5
 3.641671381563023e-5
 3.6455997988650354e-5
 3.6506957827145925e-5
 3.6522182554320104e-5
 3.6535884095673e-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 τ, Ratm = p
C14age = @. log(Ratm / R) * τ * u"s" |> u"yr"

84661-element Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(yr,), 𝐓, nothing}}}:
  867.9858872976015 yr
  880.1899484002798 yr
 1136.3776304295714 yr
 1032.1660922537071 yr
  1012.708218981491 yr
 1012.4965806538096 yr
  1030.098224932373 yr
 1014.5247521962259 yr
   995.948020099204 yr
  968.8826854822956 yr
                     ⋮
  1186.385775950517 yr
 1187.1585679085028 yr
 1189.0774285370644 yr
 1187.1689347334582 yr
 1179.1685985315798 yr
 1170.2558452939318 yr
 1158.7084282578335 yr
 1155.2616567592718 yr
 1152.1609520256272 yr

and plot it at 700 m using the horizontalslice Plots recipe

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

look at a zonal average using the zonalaverage plot recipe

plotzonalaverage(C14age, grd; color=:viridis)

or look at a meridional slice through the Atlantic at 30°W using the meridionalslice plot recipe

plotmeridionalslice(C14age, grd, lon=-30, color=:viridis)

This page was generated using Literate.jl.