System analysis: Residence times
\(\leftarrow\) BACK TO TEACHING SECTION OF MY HOMEPAGE
1 Definition and relevance
The term residence time refers to the time that a molecule or particle spends within a system (e.g. a water body) subject to in- and outflow. It is often denoted by the greek symbol \(\tau\). In the context of lakes and seas, it is usually measured in years. For technical systems like bioreactors or wastewater treatment plants, the unit of \(\tau\) is usually hours.
Why is knowledge of residence times important?
In a lake, the residence time determines the possible proliferation of planktonic organisms, i.e. phytoplankton and zooplankton.
In a wastewater treatment plant, \(\tau\) determines the time available for bacteria to degrade organic matter.
More generally, in biotechnical systems, manipulating \(\tau\) is usually a means to tune system performance, e.g. to maximize the yield of a product.
Residence times can be specified for arbitrary systems with in- and outflows. Also, \(\tau\) cann not only be defined for water but also for liquids, gas, or for solids (e.g. sludge in wastewater treatment facilities).
In the sections below, we will put the sole focus on the residence time in continuous flow stirred tank reactors (CFST).
2 Residence times in the CFSTR
2.1 Theoretical vs. actual residence time
Recall the characteristic differential equation for the concentration \(c\) of a substance in a CFSTR (Equation 1).
\[ \frac{d}{dt} c = \frac{q_{in}}{v} \cdot (c_{in} - c) \tag{1}\]
Looking at the term \(\frac{q_{in}}{v}\), you realize that its inverse, \(\frac{v}{q_{in}}\), has the unit of a time. In fact, this is the theoretical residence time of the CFSTR and, consequently, Equation 1 could alternatively expressed as Equation 2.
\[ \frac{d}{dt} c = \frac{1}{\tau} \cdot (c_{in} - c) \tag{2}\]
The fact that \(\frac{v}{q_{in}} = \tau\) is also in agreement with common sense: If, for example, a reactor with a constant volume of 10 m³ receives an inflow of 2 m³/h, the entire initial contents should be fully replaced within 5 hours of time from a budget point of view. But that view is a bit too simple …
By definition, particles and molecules in a CFSTR are subject to complete mixing. This also means that the actual position of a particular molecule or particle is a matter of coincidence. This consideration suggests that the estimate \(\frac{v}{q_{in}}\) can hardly apply to all particles or molecules. Hence, \(\tau\) appears to be kind of an average estimate only and actual residence times must follow some probability distribution.
The graphics below (Figure 1) illustrates an arbitrarily chosen probability distribution just to remind you of the basic statistical concept of distributions. Recall that both the left and the right hand side plot represent different but equivalent views on the same data.
But how does the distribution of actual residence times look like in case of a CFSTR? Does it follow a normal distribution as in Figure 1 or not?
A little literature survey would provide a quick answer. However, I will demonstrate first how some knowledge in programming can be used to find a valid answer in an empirical manner. If you understand how to implement algorithms like this, you would be able to study residence time distributions in arbitrary systems for which no theoretical model is available or applicable.
2.2 An empirical approach using agent-based simulation
The script below is a basic example of an agent-based model. In that case, the agents are individual particles that enter a continuous flow stirred tank reactor (e.g. a polymictic reservoir) via the inflow. In steady state, we know that both the reactor volume and the concentration of particles are constant and the rate of inflow equals the outflow rate.
We start with a definition of the system and its boundary conditions. Thereafter, we compute the distribution of particle ages over a series of distinct time steps. Specifically, we iterate over time steps using a for loop. Individual particles with their respective ages are stored in the vector age. In each iteration, we append new particles of zero age to this vector. Likewise, in each time step, we drop the appropriate amount of elements from the vector to represent outflow from the system. The elements to be dropped are chosen randomly so as to mimic coincidence inherent to a fully mixed system.
rm(list=ls()) # initial cleanup
set.seed(99) # optional; makes random numbers repeatable
# description of the system
v <- 100 # volume of the system (L^3)
q <- 5 # flow rate (L^3 / time)
c.in <- 20 # particle concentration in inflow (e.g. particles / L^3)
# settings for the simulation
dt <- 1 # length of an individual time step (time)
nt <- 200 # number of time steps to consider
# check settings
stopifnot(q * dt <= 0.1 * v) # time step should be < residence time
stopifnot(q * c.in >= 10) # too low input may result in bias
# initialize simulation
time <- 0
age <- numeric(0) # system has no particles initially
x <- NULL # empty variable to store results
# run simulation by iterating over time
for (time in seq(from=1, by=dt, length.out=nt)) {
age <- age + dt # update particle ages
lost <- runif(length(age)) < (dt * q / v) # select outgoing particles
age <- age[!lost] # drop outgoing particles
age <- c(age, rep(0, q * c.in * dt)) # add incoming particles
x <- rbind(x, # update result table
cbind(time=time, age=age)
)
}At this point, we have collected the result in a lengthy matrix x with two columns (“time” and “age”). So let’s start visualizing the outcome. First of all, we need to check whether the system reasonably approached a steady state indicated by a constant number of particles in the system (Figure 2). Because of random effects, we shouldn’t expect the number to be exactly constant.
# check how the number of particles in the system stabilizes over time
z <- aggregate(list(particles=1:nrow(x)), list(time=x[,"time"]), length)
with(z, plot(time, particles, type="l", ylab="Particles in system"))
We now look at the distribution of particle ages at distinct time steps. The hist() function is convenient for that purpose. Obviously, the distribution of particle ages is far from being normal (Gaussian). The results rather suggest a right-skewed distribution with a long tail (Figure 3).
# check the distribution of particle ages at particular points in time
par(mfrow=c(2,2))
for (t in c(2, 10, 100, 200)) {
hist(x[x[,"time"] == t, "age"], xlab="Particle age", main=paste("time:",t))
}
par(mfrow=c(1,1))
For a time-continuous view on the distribution of residence times, we can visualize the quantiles of the particle age (Figure 4). Again, we can see that the distribution is right-skewed as suggested by the histograms in Figure 3.
# statistics of particle ages over the entire simulation period
stats <- function(x) {
setNames(
c(mean(x), quantile(x, probs=c(0.1, 0.5, 0.9))),
c("mean", "10% percentile","median","90% percentile")
)
}
z <- array2DF(by(x[,"age"], list(time=x[,"time"]), stats))
z[,"time"] <- as.numeric(z[,"time"])
ycols <- colnames(z)[colnames(z) != "time"]
plot(range(z[,"time"]), range(z[,ycols]), type="n",
xlab="time", ylab="Particle age")
for (i in 1:length(ycols)) {
lines(z[,"time"], z[,ycols[i]], col=i)
}
legend("topleft", bty="n", lty=1, col=1:length(ycols), legend=ycols)
abline(h=v/q, lty=3)
legend(x=dt*nt/2, y=v/q, bty="n", yjust=0, lty=3, legend="Theoretical estimate (v/q)")
2.3 Formal solution for residence times in a CFSTR
The empirical results of the agent-based simulation from Section 2.2 suggest an extremely right-skewed distribution. Most of the particles pass the reactor quickly within a time span less than \(\tau\). At the same time, a few particles remain in the system for very long. Overall, the shape of the curve reminds of an exponential function with a negative argument.
In fact, a quick search brings up this paper by Nauman (2008) providing an analytical expression for the residence time distribution in a CFSTR (Equation 3). In that formula, \(F\) denotes cumulated distribution function of \(\tau\) and \(T\) represents the mean residence time as estimated by the ratio of volume and flow (cf. Equation 1 and Equation 2).
\[ F(\tau) = 1 - e^{-\tau/T} \tag{3}\]
The piece of R code below illustrates the prediction by Equation 3 together with the outcome of the agent-based model (using the age distribution of the final time step). The empirical finding is obviously very close to expectation (Figure 5).
# note: variables 'q', 'v' and 'x' defined in the agent-based model
F <- function(tau, T) {
1 - exp(-tau/T)
}
tau <- 1:100
T <- v/q
plot(tau, F(tau=tau, T=T), ylim=0:1, type="l", ylab="F(tau)")
p <- c(1, 5, 10, 25, 50, 75, 90, 95, 99) / 100
points(quantile(x[x[,"time"] == max(x[,"time"]), "age"], probs=p),
p, pch=20)
legend("bottomright", bty="n", lty=c(1,NA), pch=c(NA,21),
legend=c("Expression by Nauman 2008",
"Outcome of agent-based simulation")
)
2.4 Aspects to consider in applications to real systems
One should keep in mind that complete mixing is always an ideal case which we will hardly see in real-world natural systems, e.g. lakes. Thus, even in well mixed systems we should expect slight deviations between the actual and the theoretical distribution of residence times.
Also, even if mixing was actually complete, none of the particles entering a reactor will appear at the outflow instantly, because the transport velocity is limited. Thus, there will always be a delay not covered by Equation 3 which may be particularly pronounced in reactors of large dimensions (e.g. lakes).