Modelling foraminiferal test accumulation in Python by Andrew Berkeley

A simple model of test production and loss

Here, we'll put test production and taphononic decay together for the first time. To keep things simple, we'll drop the mixing and burial terms and consider only a surface sediment (or any fixed sediment layer). This means we don't have to consider variations in any of the variables with depth and we'll also assume they don't vary with time. Although this is not very realistic, it gives us a good starting point for understanding how test production and loss interact with one another. With these assumptions, the general equation becomes,

$$\frac {d C}{d t} = aR - \lambda C \hspace{1cm} (1) $$

This is a fairly simple differential equation which states that the rate of change in dead test concentration is equal to the rate of test production ($aR$) minus the rate of test loss ($\lambda C$). Right away we can make some simple observations. The rate of change in test concentrations will be positive if $aR > \lambda C$ and therefore concentrations will increase through time. Conversely, the rate of test accumulation will be negative if $aR < \lambda C$ and test concentrations will be decreasing. And if $aR = \lambda C$, then the rate of change in test concentrations be zero, i.e. test concentrations will stay the same.

But let's try to lok at this in more detail. The solution to equation (1) is,

$$C(t) = a \frac {R}{\lambda} \Big[1 - \exp(-\lambda t) \Big] \hspace{1cm} (2)$$

which describes how tests accumulate as function of time within a single layer (ignoring burial, mixing, etc.). Let use some Python code to plot the time evolution of equation (2) to see what happens to test accumulation.

First we need to import some useful libraries. numpy allows us to easily work with arrays of numbers, while matplotlib gives us plotting functionality

We'll evaluate equation (2) for a series of values for time. We'll use the numpy .arange() function to produce a list of numbers from 0 to 200 to represent each successive time, $t$.

Now we need to set values for the independent variables, population size ($a$), the reproduction rate ($R$) and the taphonomic decay rate ($\lambda$).

And now we can simply evaluate the equation (2) for each time in our t variable. Because t is a numpy array, we can pass it in to basic mathematical operations (e.g. multiplication, *) and the operation is evaluated for every value in the array - that is, the output is an array of results.

So now we have a new array, C, which contains the modelled dead test concentrations at each time according to the equation above. Let's plot it out to see what we notice.