Accuracy of the radiocarbon method (from a mathematical point of view)

Consider the following situation: We have an ice core and we know that the age distribution of air bubbles entrapped at depth $d$ is $p(a,d)$ with $\int_{0}^\infty p(a,d)\ \text{d}a = 1$ for all $d$. The distribution may be arbitrary, i.e. doesn't have to be normal.

We then have:

• The mean age of air bubbles at depth $d$ is

$$\overline{A}(d) = \int_{0}^\infty a\ p(a,d)\ \text{d}a$$

• The concentration of carbon $^{12}C$ at depth $d$ is

$$\overline{^{12}C}(d) = \int_{0}^\infty\ ^{12}C(-a) \ p(a,d) \ \text{d}a$$

where $^{12}C(t)=\,^{12}C(-a)$ is the concentration of $^{12}C$ in the atmosphere at time $t\leq 0$ with $t=0$ meaning "today".

• The concentration of $^{14}C$ at depth $d$ is

$$\overline{^{14}C}(d) = \int_{0}^\infty\ ^{14}C(-a) e^{-\lambda\ a} \ p(a,d) \ \text{d}a = \alpha \int_{0}^\infty\ ^{12}C(-a)\ e^{- \lambda\ a } \ p(a,d) \ \text{d}a$$

with $\alpha$ the ratio of $^{14}C/^{12}C$ in the atmosphere (which is assumed to be constant) and $\lambda$ the reciprocal of mean-life of $^{14}C$.

Applying the radiocarbon method means to calculate an age $\overline{^{14}A}(d)$ from the concentrations $\overline{^{14}C}(d)$ and $\overline{^{12}C}(d)$. This means to calculate

$$\overline{^{14}A}(d) = \frac{1}{\lambda}\ln\Big(\frac{\alpha\ \overline{^{12}C}(d)}{\overline{^{14}C}(d)}\Big)$$

Intuition tells that the age $\overline{^{14}A}(d)$ determined by the radiocarbon method must be approximately the mean age $\overline{A}(d)$ of air bubbles at depth $d$. But how to show this mathematically, and how to estimate the difference $\epsilon$ between the two depending on the shapes of $^{12}C(t)$ and $p(a,d)$ and on depth $d$?

So I want to show that

$$\overline{^{14}A}(d) = \overline{A}(d) + \epsilon(d)$$

with some (hopefully) small $\epsilon(d)$ or, respectively, find the $\epsilon(d)$ with

$$\frac{1}{\lambda}\Bigg(\ln\alpha + \ln\Big(\int_{0}^\infty f(a) \ p(a,d) \ \text{d}a\Big) - \ln\Big(\int_{0}^\infty f(a) \ p(a,d) \ e^{-\lambda\ a}\ \text{d}a\Big)\Bigg)\\ = \int_{0}^\infty a\ p(a,d)\ \text{d}a + \epsilon(d)$$

for good-natured functions $^{12}C(t)$ and $p(a,d)$. Let be $f(a) = ^{12}\kern-0.24emC(-a)$.

If I already made mistakes in the derivation of the problem, please excuse me and drop me a note. Otherwise I am keen to get a solution. How can – for example – it be shown that the difference $\epsilon$ is small when $p(a,d)$ is a narrow distribution around a mean, and that it is even smaller when it is a symmetric or even normal distribution?

Addendum:

The question is simpler when we consider the specific function $^{12}C(t) = k\ t$. Then we have

$$\frac{1}{\lambda}\Bigg(\ln\alpha + \ln\Big(\int_{0}^\infty a \ p(a,d) \ \text{d}a\Big) - \ln\Big(\int_{0}^\infty a \ p(a,d) \ e^{-\lambda\ a}\ \text{d}a\Big)\Bigg)\\ = \int_{0}^\infty a\ p(a,d)\ \text{d}a + \epsilon(d)$$

Can $\epsilon(d)$ be given exactly?

The error $\epsilon(d)$ can in fact be given exactly when $p(a,d)$ is a delta function, i.e. $p(a,d) = \delta(d - a\rho)$. Then we have

$$\frac{1}{\lambda}\Big(\ln\alpha + \ln f(d/\rho) - \ln\big(f(d/\rho)\ e^{-\lambda\ d/\rho}\big)\Big) = \ln\alpha/\lambda + d/\rho \\= d/\rho + \epsilon(d),$$

that means $\epsilon(d) = \epsilon = \ln\alpha / \lambda$. For small $\alpha$ and large $\lambda$ (as in the case of $^{14}C$), this error may be neglected. But I am interested in the general case.

There are two issues here. The first is that the carbon method gives a certain sort of average with density $\asymp C(-a)p(a)$ (I omit $d$ because it is the same everywhere), not $p(a)$, so even if the logarithm of the ratio were a perfect mean age estimator, it would not give the mean for the true distribution, but the mean for the weighted distribution. The second issue is the logarithm and integrals itself. Your only real hope is that $p$ is essentially concentrated on an interval over which neither $C$, not $e^{-\lambda t}$ vary noticeably. Then you are in good shape.

Let's say that $p$ is normal with mean $0$ and variance $1$ (I leave the scaling to you), $\log C$ is $\mu$-Lipshitz with small $\mu$ and $\lambda$ is also small. Then the worst case scenarios are when you take the extreme cases $C(a)=e^{\pm\mu a}$ resulting in the deviations between $\frac\lambda 2-\mu$ and $\frac \lambda 2+\mu$.

Unfortunately, the picture here shows that this Lipschitzness of $\log C(a)$ is a too optimistic assumption and one can only rely on the statement that the concentration was changing between some level and $1.5$ times that level (going down was not too sharp, but going up was rather dramatic each time). In that case the worst case scenarios are when $C(a)$ jumps from $1$ to $1.5$ at some point. This doesn't result in any closed form formula, but, just running the honest computation, we see that the error now is within $0.3$ units for $\lambda$ up to $0.3$ and stays at $0.15$ even when $\lambda$ approaches $0$, which by itself isn't that bad (any amount up to one standard deviation in the true age distribution, which is normalized to $1$ here shouldn't bother us too much and that will be attained only at $\lambda\approx 1.8$).

So, the key questions are

1. Is normality of the true age distribution a reasonable assumption?

2. What is the (estimated) standard deviation in the age distribution?

3. What is the half-life of ${}^{14}C$? (yeah, I can look that one up but you should know it by heart).

Once we know that, we can, probably, make some decent conclusions.

Edit: Let me just elaborate on what can be concluded about the error without any assumption whatsoever. As I said in the comments, you can then underestimate the true mean age by an arbitrarily large amount. However, assuming that the atmospheric carbon concentration $C$ was bounded between some level and $K$ times that level over the accumulation period, we can bound the possible overestimate. I will denote by $p$ the density of the deviation from the true mean, so $\int xp=0$ and normalize $\lambda$ to $1$ (so our unit of time is about 8267 years). Then we want to figure out how negative the quantity $\log\frac{\int Cpe^{x}}{\int Cp}$ be.

Note that if $C$ were constant, we would have the ratio of integrals at least $1$ by Jensen, so no overestimate would be possible at all. Now if we can find a probability distribution with $\int Cpe^x-u\int Cp=-v$ and $\int px=0$ for some positive $u<1$ and $v$, then we can satisfy the same equations by a linear combination of 2 Dirac measures, which will be non-degenerate because $v>0$, $u<1$ (we may lose the normalization $\int p=1$ but it doesn't matter). Thus, it is enough to find the minimum of the ratio $$ \frac{Kbe^{-a}+ae^{b}}{Kb+a}, a,b>0\,. $$ Using AM-GM, we see that it is at least $$ K^{\frac b{a+b}}\frac{a+b}{Kb+a}\,. $$ Normalizing to $a+b=1$ here, we get $ K^b/[(K-1)b+1] $ to minimize over $[0,1]$. This leads to the equation $$ \log K[(K-1)b+1]-(K-1)=0\,, $$ i.e., $$ b=\frac 1{\log K}-\frac 1{K-1} $$ and finally to the logarithm of the minimum equal to $$ \log\frac{\log K}{K-1}-\frac{\log K}{K-1}+1 $$ units of time. For $K=1.5$, that is mere $-0.02$ units (about $170$ years). When $K=4$ (million years back in time), then it is $-0.234$ units (under $2000$ years). Even beyond that the belief is that $K\le 10$, which gives $-0.62$ units (about $5000$ years). Finally, even if we assume $K=100$ (quite an unbelievable atmosphere), then we still only get $-2.11$ units (under $18000$ years).

The moral is that an essential overestimate of the mean with the carbon method is pretty much impossible, so if you measured $X$ years, you can quite safely say that the ice is at least $X$ years old, but, unfortunately, this is where your knowledge ends. No upper bound is feasible without making assumptions.