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()

(, 
⠻⣦⡀⠀⠘⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠲⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠀⠙⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠈⠻⣦⡀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⢤⡀⠈⠻⣦⡀⠈⠳⣄⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠈⠻⣦⡀⠈⠳⣄⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠈⠻⣦⡀⠈⠳⣄⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠈⠻⣦⡀⠘⢦⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⣀⠈⠻⣦⡈⢧
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠓⠦⣌⠻⣦)

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.781368293693143e-5
 3.7757899707203986e-5
 3.6605707902986216e-5
 3.707009029891931e-5
 3.715744797475248e-5
 3.7158399272783746e-5
 3.707936439247613e-5
 3.7149283799004885e-5
 3.723285923778942e-5
 3.7354960954899674e-5
 ⋮
 3.638493414383458e-5
 3.638153292381407e-5
 3.637308898920067e-5
 3.638148729939723e-5
 3.6416713815628925e-5
 3.645599798864905e-5
 3.650695782714463e-5
 3.652218255431881e-5
 3.65358840956717e-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.9858872977056 yr
  880.1899484003853 yr
 1136.3776304296994 yr
 1032.1660922538545 yr
 1012.7082189816208 yr
 1012.4965806539328 yr
 1030.0982249324732 yr
 1014.5247521963447 yr
  995.9480200993163 yr
   968.882685482441 yr
                     ⋮
  1186.385775950803 yr
 1187.1585679087918 yr
 1189.0774285373586 yr
  1187.168934733754 yr
  1179.168598531876 yr
 1170.2558452942267 yr
 1158.7084282581272 yr
 1155.2616567595655 yr
 1152.1609520259226 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.