\(\leftarrow\) BACK TO TEACHING SECTION OF MY HOMEPAGE
Lakes and seas typically receive inflows from the catchment and they are thus subject to inputs of matter. This includes, e.g., nutrients like nitrogen and phosphorus which cause eutrophication or harmful contaminants like pesticides. Managers and/or scientists are then confronted with questions like:
How fast does the in-lake concentration respond to altered input loads?
How long does it take until the concentration of a pollutant falls below the sub-critical level again?
The term continuous flow stirred tank reactor, short CFSTR, denotes water bodies which are subject to both complete mixing as well as permanent in- and outflow (Figure 1).
Both natural and technical systems can be regarded as CFSTR: polymictic lakes, tanks used in conventional wastewater treatment, and bioreactors as employed in industry and research (Figure 2).
Besides actual systems, the CFSTR concept is of general importance in process-based modeling. Namely in spatially discretized models, the CFSTR concept applies to each single control volume.
To answer the questions raised in the motivation section, we need an equation that allows the dynamics of the concentration to be computed. As usual, such an equation is derived from mass balances. In fact, as long as we consider a conservative constituent (i.e. one that does not undergo decay, transformation, or growth), all we need to consider is the mass balance of the constituent (Equation 1) and the mass balance of water (Equation 2). Using the symbols introduced in Figure 1, we can write:
\[ \frac{d}{dt} (c \cdot v) = q_{in} \cdot c_{in} - q_{ex} \cdot c \tag{1}\]
\[ \frac{d}{dt} v = q_{in} - q_{ex} \tag{2}\]
Note that the product \(c \cdot v\) in Equation 1 represents the total mass in the system and the products on the right hand side are mass fluxes (a.k.a. loads).
That’s it. When we consider the transport of reactive substances and suspended organisms later, these two differential equations will still form the basic physical framework.
While Equation 1 correctly expresses the mass balance, we actually seek an expression for the change in concentration. We can obtain it by applying the product rule of differential calculus \((u v)' = u v' + u' v\) as follows:
\[ \frac{d}{dt} (c \cdot v) = c \cdot \frac{d}{dt} v + v \cdot \frac{d}{dt} c \]
Considering that we can replace the term \(\frac{d}{dt}v\) by the right hand side of Equation 2, a bit of rearrangement yields the characteristic expression below (Equation 3):
\[ \frac{d}{dt} c = \frac{q_{in}}{v} \cdot (c_{in} - c) \tag{3}\]
What does Equation 3 tell us? Essentially, two things:
If \(c_{in} > c\), the concentration will increase and, if \(c_{in} < c\), it will decline in agreement with common sense.
How fast the concentration \(c\) responds to an external perturbation (i.e. a situation where \(c_{in} <> c\)) depends on the term \(q_{in}/v\).
What does \(q_{in}/v\) actually represent? Looking at the unit, it turns out that it is the inverse of the theoretical water residence time (often denoted by the greek letter tau). We can conclude that, for a given inflow rate, the concentration changes faster, the smaller the lake volume is.
In any case, by setting the left hand side of Equation 3 to zero, we can see that the steady-state solution is \(c = c_{in}\). That is, for any constant value of \(c_{in}\), the concentration in the reactor will sooner or later converge to that value.
To use Equation 3 for the prediction of future concentrations in the CFSTR, we must get rid of the differential: hence, we need integration. Fortunately, this equation is simple enough to allow for an analytical solution with moderate skills in math.
We start with a little rearrangement
\[ \frac{d}{dt} c = \frac{q_{in}}{v} \cdot (c_{in} - c) = -\frac{q_{in}}{v} \cdot c + \frac{q_{in}}{v} \cdot c_{in} \]
and then introduce two new constants, \(a\) and \(b\), merely to simplify writing and improve clarity. With the two new constants \(a = -q_{in}/v\) and \(b = q_{in}/v \cdot c_{in}\), the differential equation to be integrated becomes
\[ \frac{d}{dt} c = a \cdot c + b \tag{4}\]
with the right hand side being a simple linear equation. We will now use the substitution method of integration. To do so, we substitute the entire right hand side by a new variable \(z\) related to the original variable by \(z = a \cdot c + b\):
\[ \frac{d}{dt} c = z \tag{5}\]
In addition, we must consider the relation between differentials \(dz\) and \(dc\). According to the basic rules of differentiation, the following holds (written very explicitly):
\[ \frac{dz}{dc} = \frac{d}{dc}z = \frac{d}{dc} (a \cdot c + b) = a \]
It follows from the last expression that \(dc = 1/a \cdot dz\). We can now use this to replace the \(dc\) in Equation 5 to finally obtain a differential equation containing \(z\) as the single variable.
\[ \frac{dz}{dt} = a \cdot z \]
Separation of variables and putting the integral operators leads to
\[ \int{\frac{1}{z}} \cdot dz = \int{a} \cdot dt \]
Looking into a table of fundamental integrals, we find that \(\int{(1/z)} \cdot dt = ln(z)\) and \(\int{a} \cdot dt = a \cdot t\). Remembering that indefinite integrals involve an unknown integration constant (which I will call \(x\)), we arrive at
\[ ln(z) = a \cdot t + x \]
This expression can be rearranged as follows where, in the last step, the unknown integration constant was replaced by another constant with a new symbol (\(x_2 = e^{x_1}\)) for convenience.
\[ \begin{align} z &= e^{a \cdot t + x} \\ &= e^{a \cdot t} \cdot e^{x} \\ &= e^{a \cdot t} \cdot x_2 \end{align} \]
We can now revert the variable substitution that took us from Equation 4 to Equation 5 and obtain
\[ a \cdot c + b = e^{a \cdot t + x} \]
which can be rearranged to
\[ \begin{align} c &= \frac{1}{a} \cdot \left( e^{a \cdot t} \cdot x_2 - b \right) \\ &= e^{a \cdot t} \cdot x_3 - \frac{b}{a} \end{align} \tag{6}\]
where the unknown integration constant was again redefined for convenience (\(x_3 = x_2 / a\)). At this step, we can determine the value of that unknown constant by exploiting information on the initial condition. The initial condition is the state of the system at zero time. If, we take Equation 6, set \(t=0\) and denote the respective initial concentration \(c_0\), we find for the unknown constant:
\[ x_3 = c_0 + \frac{b}{a} \tag{7}\]
Introducing Equation 7 into Equation 6 essentially took us to the solution of the definite integral:
\[ c = e^{a \cdot t} \cdot \left(c_0 + \frac{b}{a} \right) - \frac{b}{a} \]
The final step of the procedure is to replace the auxiliary constants \(a\) and \(b\), first introduced in Equation 4, by the respective original expressions. Doing so, we end up with the solution
\[ c = \left( c_0 - c_{in} \right) \cdot e^{\left( -\frac{q_{in}}{v} \cdot t \right)} + c_{in} \tag{8}\]
Whenever a new mathematical expression is obtained by a somewhat longer sequence of calculations, it needs to be verified. At least, one should perform all possible checks for general plausibility. In our case, it helps to consider the predicted value of the concentration in the reactor, \(c\), for a set of specific assumptions where the true outcome is known from common sense. The list below demonstrates typical considerations:
If \(t = 0\) we predict \(c = c_0\): OK, consistent with our initial condition.
If \(t \rightarrow \infty\) we predict \(c \rightarrow c_{in}\): OK, after a very long time, all the water initially present in the reactor was replaced by incoming water.
If \(v \rightarrow \infty\) we predict \(c \rightarrow c_0\): OK, in a very large water body, the concentration will hardly change whatever the input is.
If \(v = 0\) we predict \(c \rightarrow c_{in}\): OK, the contents of a negligbly small water body is replaced be the incoming water instantly.
After such very basic checks, it is always a good idea to create plots that help understanding the obtained equation better. For example, we can plot the dynamics of the concentration in an initially clean reactor in response to a sudden, permanent pollution of the inflow. The values of volume (v), flow (q.in), and incoming concentration (c.in) use in the code below are plain assumptions.
ct <- function(t, v, c.0, q.in, c.in) { # our analytical solution
(c.0 - c.in) * exp(-q.in/v*t) + c.in
}
times <- 0:30 # desired time points
prd <- ct(t=times, c.0=0, v=250, q.in=50, c.in=1) # predicted values
plot(times, prd, type="l", xlab="time", ylab="c(t)") # creates the plot
Once an analytical solution has been verified, it can be used in serious applications. However, you will find that there is a tendency to generally solve integration problems numerically rather than analytically. There are probably two reasons. First, many differential equations encountered in practice are too complex to be mathematically tractable because expressions do not match known fundamental integrals. Second, the use of numerical solvers is nowadays very convenient, probably more convenient than deriving lengthy expression on paper. But still, analytical solutions are relevant for two reasons:
They are exact solutions, so they can be used to check numerical integration algorithms for accuracy.
They are way cheaper in terms of computational resources as compared to numerical solutions, hence they are preferred if speed or memory matters.
A simple physical setup for simulating a CFSTR is illustrated in Figure 4. Since a low but truly constant inflow rate is somewhat difficult to achieve, water was transferred into the reactor in a semi-continuous manner using a syringe with 10 mL capacity. Further, plain salt (NaCl) was used as the model substance, because electric conductivity could then be used as a proxy for concentration. As long as the added amount of salt remains far below solubility, the relation between NaCl concentration and electric conductivity is reasonably linear.
Using the setup from Figure 4, one can run experiments to cover two contrasting scenarios.
The water in the blue beaker representing the CFSTR could be initially filled with tap water. The clean tap water is then gradually replaced with salty water from the large green reservoir shown on the right.
The water in the CFSTR (blue beaker) couly initially contain salt which is continuously diluted by clean tap water from the reservoir (green beaker).
The R code sections below hold data from actual experiments covering both of the two scenarios. The first section of R code builds a data frame holding the observed values of electric conductivity (column conObs) in the blue beaker from Figure 4 over time. Instead of times, the table rather stores the number of injections carried out with the syringe (column injections). Note that, in this piece of R code, the data.frame constructor was used instead of read.table merely for a more compact notation.
dynam <- data.frame(
experiment=c(rep("initiallyPolluted", 26), rep("initiallyClean", 27)),
injections=c(0:25, 0:26),
conObs=c(7910,7500,7100,6770,6440,6080,5780,5510,5250,4970,4720,4490,4270,
4060,3880,3710,3550,3340,3180,3050,2910,2760,2640,2520,2420,2300,514,1139,
1656,2210,2690,3140,3580,3990,4370,4720,5060,5380,5670,5950,6220,6460,6700,
6920,7130,7320,7510,7690,7860,8010,8170,8310,8430)
)The second section of R code builds a data frame holding constant metadata needed to evaluate the observed dynamics. Note that the abbreviation “con” is used in the names of parameters to denote conductivity (as a proxy for concentration).
const <- read.table(header=T, sep="", text='
experiment parameter value unit comment
initiallyPolluted reactorVolume 175 mL full_capacity_of_blue_beaker
initiallyPolluted conInitial 7910 µS/cm initial_con._in_blue_beaker
initiallyPolluted conInflow 514 µS/cm con._in_green_reservoir
initiallyPolluted injectionVolume 10 mL capacity_of_single_syringe
initiallyPolluted injectionInterval 15 s seconds_between_injections
initiallyClean reactorVolume 175 mL ...
initiallyClean conInitial 514 µS/cm ...
initiallyClean conInflow 10700 µS/cm ...
initiallyClean injectionVolume 10 mL ...
initiallyClean injectionInterval 15 s ...
')The R code below visualizes the observed dynamics for the two experiments. Along with this, the code also shows the corresponding expected dynamics as predicted by Equation 8. Apparently, the prediction is in close agreement with observations.
# load a package to facilitate pivoting of tables by 'dcast' and 'meld';
# this particular packages comes with a minimum of dependencies
library("data.table")
# implements the analytical solution for the concentration in the CFSTR
cfstr <- function(
t, # time (we will assign a vector of times as actual argument)
c.0, c.in, # c.0: initial concentration, c.in: concentration of inflow
q.in, v # q.in: inflow rate, v: volume of reactor
) {
(c.0 - c.in) * exp(-q.in/v * t) + c.in # solution found by integration
}
# transform the static data from long to wide format; using package 'reshape2'
const <- data.table::dcast(experiment ~ parameter, value.var="value",
data=as.data.table(const))
# merge dynamic and constant data into a single table for convenience
x <- merge(dynam, const, by="experiment")
rm(dynam, const) # not needed any longer
# translate the number of injections into times
x <- cbind(x, seconds=with(x, injections*injectionInterval))
# add expected concentrations according to the CFSTR model
x <- cbind(x, conMod=with(x, cfstr(t=seconds,
c.0=conInitial, c.in=conInflow,
q.in=injectionVolume/injectionInterval, v=reactorVolume)))
# compare observations against model predictions for individual experiments
# step 1: a function to plot a subset of the table representing an experiment
fn <- function(x) {
with(x, {
plot(range(seconds), range(c(conInitial, conInflow)), # empty plot with
type="n", xlab="seconds", # proper scales
ylab="concentration (as conductivity)")
points(seconds, conObs) # add observations
lines(seconds, conMod) # add predictions
lines(seconds, conInitial, lty=2) # as a reference
lines(seconds, conInflow, lty=3) # as a reference
mtext(side=3, unique(experiment), cex=par("cex")) # label on top
if (unique(experiment) == "initiallyClean") {
legend(x=20, y=1e4, bty="n", pch=c(1,NA,NA,NA), lty=c(NA,1,2,3),
legend=c("observed","predicted","initial","in inflow")
)
}
})
NULL # no return value; only to be called for the side effect of plotting
}
# step 2: actual plotting
par(mfcol=c(1,2)) # side-by-side layout for two experiments
x <- by(x, x[,"experiment"], fn) # subset table 'x' by experiment and call 'fn'
par(mfcol=c(1,1)) # reset layout
Note that the above code uses some advanced features of the R language. If you feel that the brief explanations below are not sufficient, have a look into R’s help pages.
with: Used to break up lists into its components (i.e. columns if applied to data frames). The components can then be used like normal variables in subsequent statements. Using with is never strictly necessary but it makes the code shorter and more readable.
dcast: A function imported from the data.table package. It transforms a table from long format to wide format. Just put a print(const) before and after the line containing the call to dcast and you will instantly understand what it does.
merge: Used to join two tables based on a common column. Here, we use it to associate the static meta data with every single observation. Although not very efficient in terms of memory, this helps to make the subsequent code straightforward.
par: The key to tuning many features of plots. Here, we use it to arrange multiple plots into a single figure.
by: Meant to split a data frame (i.e. a table) into horizontal slices (i.e. subsets of rows) based on unique entries in a particular column(s). Each of the created subsets is then pushed into a user-selected function for further processing. Here, we use by to create an individual plot for each experiment.