So, lets omit mixing and taphomic decay (for simplicity) but retain the burial term. As before, let's assume the independent variables (e.g. $a$, $R$) are constant (do not vary with depth, $x$, or time, $t$). So we have, from the General Equation,

$$\frac {\partial C}{\partial t} = aR - w\frac {\partial C}{\partial x}$$

Notice that this has both $x$ and $t$ dimensions and therefore test concentations vary with both time and depth into the sediment. This means we do not simply have one value for the concentration at each point in time. Rather we have a profile of dead test concentration with depth at each point in time. Also, notice that, like before, the rate at which concentration changes through time ($\frac {\partial C}{\partial t}$) is described by the production rate minus some other term. The other term is the rate of change in test concentration due to burial. But are tests always lost due to burial? Well, no. That depends on whether the $\frac {\partial C}{\partial x}$ term is positive or negative, in other words, if dead test concentrations are increasing or decreasing with depth.

Well, the simplest way to solve this partial differential equation is to do it numerically via discretization and difference equations. So let's discretize it by approximating the derivatives,

$$\frac {C_{x, t+1} - C_{x,t}}{\Delta t} = aR - w \bigg( \frac {C_{x,t} - C_{x-1,t}}{\Delta x} \bigg)$$

Now solving for $C_{x,t+1}$ gives us an expression which relates the concentration at one point in time to the concentraton at the previous point in time. We can use this to propagate the model through time and generate a solution.

$$C_{x,t+1} = C_{x,t} + \Delta t \bigg[aR - w \bigg( \frac {C_{x,t} - C_{x-1,t}}{\Delta x} \bigg) \bigg]$$

In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [3]:
x_max = 10.0     # total depth into sediment to model
nx    = 100      # number of depth steps
dx    = x_max/nx # size of depth step

t_max = 100      # total time to model
dt    = dx*0.1   # size of time step
nt    = t_max/dt # number of time steps

C = np.zeros((nt, nx))

In [4]:
a = 25.0 # standing crop size
R = 2.0  # reproduction rate
w = 0.5  # sedimentation rate

for ts in np.arange(1,nt-1):
    C[ts,0]  = dt*a*R # dC/dx = 0 at x = 0 
    C[ts,1:] = C[ts-1,1:] + dt * (a*R - w*(C[ts-1,1:] - C[ts-1,0:-1])/dx)

In [5]:
fig = plt.figure()

p = fig.add_subplot(111, xlim=(0, nx), ylim=(np.min(C)*1.1, np.max(C)*1.1))
plt.plot(C[200,:], lw=3)
plt.plot(C[500,:], lw=3)
plt.plot(C[1000,:], lw=3)
plt.plot(C[1500,:], lw=3)
plt.plot(C[5000,:], lw=3)
plt.grid()
plt.xlabel('depth, x')
plt.ylabel('dead test concentration, C')

plt.show()



In [0]:


In [0]:


In [0]:


In [0]:


In [ ]: