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, \(\tau\) is typically given in hours.
Why is knowledge of residence times important?
In a lake, the residence time determines the possible proliferation of suspended organisms, i.e. phytoplankton and zooplankton.
In a wastewater treatment plant, \(\tau\) determines the time available for bacteria to degrade organic matter and increase their own biomass.
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\) can 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 (CFSTR).
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 will find 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 an average estimate only while 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 has been published (and even AI would not have the answer).
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 while 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 clean-up
# assumptions on the reactor and its boundary conditions
v <- 100 # storage volume (m³)
q <- 5 # constant flow rate (m³/h)
c.in <- 10 # inflow concentration (particles / m³)
# simulation settings (time unit to be consistent with the flow rate)
dt <- 1 # length of a simulation time step (h)
tmax <- 240 # let simulation end at this time (h)
# check inputs: we want only a minor fraction of the reactor's volume
# to be replaced in each time step; in case of failure, decrease "dt"
stopifnot(q * dt <= v / 10)
# initialization
t <- 0 # let the time start at zero
age <- c() # vector of particle ages is initially empty
res <- NULL # a matrix to store results; this is initially empty as well
# simulation
while (t < tmax) { # enter time loop
t <- t + dt # advance the simulation time
if (length(age) > 0) { # increment age of particles
age <- age + dt # in the system, if any
}
age <- c(age, rep(0, q * c.in * dt)) # add entering particles
stay <- runif(n=length(age)) > (q * dt / v) # which particles stay ?
age <- age[stay] # keep those, drop leaving ones
if (length(age) > 0) { # save results for this time step
res <- rbind(res, cbind(time=t, age=age)) # (append rows to matrix "res")
}
} # end of time loopAt this point, we have collected the result in a matrix res with two columns (“time” and “age”) but many rows. So let’s visualize the outcome using graphics. First of all, we need to check whether the system reasonably approached a steady state within the chosen simulation time (tmax). We can use two possible indicators (Figure 2):
the number of particles in the system must be constant with some random variation around the expected value (equal to concentration in inflow).
the mean particle age should have converged to a noisy constant level without notable trend.
layout(matrix(1:2, ncol=2))
# Dynamics of particle concentration
x <- aggregate(list(count=1:nrow(res)), list(time=res[,"time"]), length)
plot(x[,"time"], x[,"count"] / v, xlab="Hours", ylab="Particles / m³")
# Dynamics of mean particle age
x <- aggregate(age ~ time, data=res, mean) # using aggregate with formula
plot(x[,"time"], x[,"age"], xlab="Time (h)", ylab="Mean particle age (h)")
layout(1)
Instead of plotting the mean particle age over time (Figure 2), we could also plot quantiles of the particle age (Figure 3). Such a plot gives as an indication that the age distribution is not normal but skewed. There are two obvious indications:
the median (50% quantile) and the mean value differ,
lower and higher quantiles of corresponding probability levels (e.g. 25% and 75%; same color) are not symmetric about the median.
# empty plot with appropriate axes
plot(age ~ time, data=res, type="n", xlab="Time (h)", ylab="Particle age (h)")
# add selected quantiles
p <- c(grey80=1, grey50=25, black=50, grey50=75, grey80=99)
for (i in 1:length(p))
lines(age ~ time, data=aggregate(age ~ time, res, quantile, probs=p[i]/100),
col=names(p[i]))
# add mean with different line style
lines(age ~ time, data=aggregate(age ~ time, res, mean), lty=2)
legend("topleft", bty="n", title="Particle model", lty=c(rep(1, length(p)), 2),
col=c(names(p), "black"), legend=c(paste0(p,"% quantile"), "mean"))
# add a line for the theoretical residence time (volume divided by flow)
abline(h=v/q, lty=4, col="red")
legend("topright", bty="n", lty=4, col="red", legend="Quotient v/q")
In Figure 3, I also added a straight horizontal line to mark the theoretical residence time computed as the quotient of volume and flow rate (recall Equation 1 and Equation 2). Obviously, this estimate is perfectly representative of the mean particle age. But a particular particle can stay much shorter or longer in the system.
We now look at the distribution of particle ages at the very end of the simulation period. The hist() function is convenient for that purpose. Obviously, the distribution of particle ages is truly far from being normal (Gaussian). The results rather suggest a heavily right-skewed distribution with a long tail (Figure 4).
rows <- res[,"time"] == max(res[,"time"])
hist(res[rows, "age"], xlab="Particle age (h)", main="") 
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. The majority of particles passes 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 'res' defined in the agent-based model above
F <- function(time, v_over_q) {
1 - exp(-time / v_over_q)
}
some_ages <- 1:100
plot(some_ages, F(time=some_ages, v_over_q=v/q), ylim=0:1, type="l",
xlab="Particle age (h)", ylab="Probability (AGE < age)")
p <- c(1, 5, 10, 25, 50, 75, 90, 95, 99) / 100
points(quantile(res[res[,"time"] == max(res[,"time"]), "age"], probs=p),
p, pch=20)
legend("bottomright", bty="n", lty=c(1,NA), pch=c(NA,20),
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).