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)

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 20 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")

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

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"

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.