| Application | Example of state variable | Example |
|---|---|---|
| Readioactivity | Amount of isotope Cs137 | k=-0.023 year-1 |
| Microbial degradation | Leaf-litter biomass | link to publication |
| Decay of toxic substrance | Pesticide concentration | link to publication |
| Settling in mixed water column | Suspended detritus concentration | \(k=-u/z\) with settling velocity \(u\) and depth \(z\) |
| Growth of microorganisms | Density of \(E. coli\) culture | \(k=2.3\) hour-1 for optimum temperature and substrate supply |
System analysis: First-order processes
\(\leftarrow\) BACK TO TEACHING SECTION OF MY HOMEPAGE
1 Definition and examples
The term order is traditionally employed in chemical engineering. It is used to describe the dependency of the reaction rate on the concentration of reactants (basics can be found on this wikipedia page, for example). We will focus primarily on first-order processes characterized by Equation 1
\[ \frac{d}{dt} c = k \cdot c \tag{1}\]
where \(c\) denotes the concentration of a substance and \(k\) is constant characterizing the reaction velocity. In the original context, the differential equation states that, in a chemical reaction, the rate of consumption of a reactant (expressed by \(-\frac{d}{dt} c\); note the negative sign) is proportional to the (remaining) amount of that reactant (expressed by concentration \(c\)). Obviously, the proportionality constant \(k\) has the unit of “per time” (e.g. hour-1).
It turns out that the concept has many applications in ecology, some being related to biochemical reactions, others being unrelated to chemistry at all. Some prominent examples are listed in Table 1. Note that \(k\) can take positive or negative values, depending on whether the state variable increases or undergoes a reduction.
2 Analytical solution
The governing differential equation (Equation 1) can be integrated instantly by separation of variables. The following sequence of operations takes us to the indefinite integral
\[ \begin{aligned} \frac{d}{dt} c &= k \cdot c \\ \frac{1}{c} \cdot dc &= k \cdot dt \\ \int \frac{1}{c} \cdot dc &= \int k \cdot dt \\ ln(c) &= k \cdot t + x \\ c &= e^{k \cdot t} \cdot e^x \\ c &= e^{k \cdot t} \cdot x_2 \\ \end{aligned} \]
where the integration constant is denoted \(x\). Since \(x\) is unknown, the same applies to \(e^x\) and we can thus replace the latter by a new constant \(x_2\) for clarity. To eliminate the latter, we specify an initial condition. That is, we define the concentration at time zero to be some known value \(c_0\). Setting \(t=0\) in the latest expression from above, we find \(c_0 = x_2\). Consequently, the solution of Equation 1 is the well-known expression shown as Equation 2.
\[ c(t) = c_0 \cdot e^\left( k \cdot t \right) \tag{2}\]
3 Proof for the case of particle settling
In Table 1, the settling of particles in a turbulent environment was listed as an example of a first order loss process. According to Equation 2, the concentration of particles should converge to zero in an exponential manner. Since this is not necessarily intuitive, the little agent-based model below provides an empirical proof. In each time step, particle are randomly distributed over depth to simulated the effect of turbulent mixing. Once a particles touches the bottom in response to settling, it is no longer counted. Note that, in this example, the value of the rate constant \(k\) is equal to the quotient of the settling velocity (taken negative) and the depth.
rm(list=ls())
set.seed(999)
depth <- 1 # assumed water depth (m)
velo <- 0.2 # assumed settling velocity (m/day)
ninit <- 1e3 # initial amount of particles in system
dt <- 1 # duration of a time step (day)
ntimes <- 32 # number of time steps considered
n <- ninit # initialize
result <- c(time=0, n=n)
for (i in 1:ntimes) { # enter time loop
z <- runif(n=n) * depth # push particles to a new random depths (turbulence)
z <- z + velo * dt # let particles sink
n <- sum(z < depth) # count particles being still in suspension
result <- rbind(result, # record new state in result matrix
c(time=i*dt, n=n))
} # done with simulation
par(mfcol=c(1,2))
plot(n~time, data=result, xlab="Time", ylab="Susp. particles") # empirical
lines(result[,"time"], ninit*exp(result[,"time"] * -velo/depth)) # theory
legend("topright", bty="n", pch=c(1,NA), lty=c(NA, 1),
legend=c("Particle simulation","Analytical solution"))
plot(n~time, data=result, xlab="Time", ylab="Susp. particles", # log scale
log="y")
par(mfcol=c(1,1))
4 Converting rate constants to half-life or doubling time
Sometimes, it is more convenient to express the rate of a first-order process in units of a time (rather than time-1). This leads to the concept of the half-life time or doubling time (\(T\)). The latter is defined as the time required to achieve a change in the value of the state variable by factor 1/2 or 2, respectively.
The relation between \(k\) and \(T\) can simply be found by assuming \(c_0 = 1\) in Equation 2 and assuming either \(c(T) = 2\) or \(c(T) = 1/2\), depending on whether a growth or loss process is studied. Doing so leads to Equation 3.
\[ \begin{aligned} T = & \frac{ln(1/2)}{k} & \text{for } k < 0 \text{, i.e. decay} \\ T = & \frac{ln(2)}{k} & \text{for } k > 0 \text{, i.e. growth} \end{aligned} \tag{3}\]
5 First-order processes in the CFSTR
5.1 Basic equation and solution strategies
We will now combine the basic equation for the transport of a conservative substance in a stirred tank reactor with the equation for a first-order process (Equation 4).
\[ \begin{aligned} \frac{d}{dt}c & = \frac{q_{in}}{v} \cdot (c_{in} - c) + k \cdot c \\ & = \frac{1}{\tau} \cdot (c_{in} - c) + k \cdot c \end{aligned} \tag{4}\]
This equation enables us to predict the concentration of a constituent \(c\) in a perfectly mixed water body in response to:
the external input (expressed by the load \(c_{in} \cdot q_{in}\)),
the physical dimensions of the reactor (affecting the residence time \(\tau = v/q_{in}\)),
the turnover of the constituent according to first-order kinetics (expressed by the rate constant \(k\)).
As always, Equation 4 requires integration if one wants to predict the dynamics of the concentration \(c\) in response to altered boundary conditions. Good to know that this equation can still be solved analytically by the substitution method. However, in the subsequent sections we will rather focus on a steady-state solution of Equation 4. In later chapters of the course, we will then integrate equations like Equation 4 using numerical algorithms.
5.2 Derivation of the steady-state solution
The steady-state concentration in a CFSTR subject to a first-order reaction process is obtained by setting the left hand side of Equation 4 to zero. A few steps of simple rearrangement yield:
\[ c = \frac{c_{in}}{1 - k \cdot \tau} \tag{5}\]
5.2.1 Solutions for negative rate constants (decay)
We consider the common case of a negative constant \(k\) representing, e.g., microbial degradation of and organic constituent or elimination due to settling. Let’s understand how steady-state concentrations depend on the value of \(\tau\) or \(k\), respectively, by plotting Equation 5. As illustrated by Figure 2, the longer the residence time, the more pronounced is the “cleaning” effect. Similarly, larger absolute values of \(k\) result in a more efficient removal of the constituent.
rm(list=ls())
c.std <- function(c.in, tau, k) { # the steady-state solution
c.in / (1 - k * tau)
}
c.in <- 1 # assumed inflow concentration
k <- c(0, -0.01, -0.1, -1) # distinct values of 'k'
tau <- 0:50 # vector of residence times
plot(range(tau), c(0, c.in), type="n", # empty plot
xlab="Residence time",
ylab="Steady state concentration")
for (i in 1:length(k)) # add prediction for each 'k'
lines(tau, c.std(c.in, tau, k[i]), col=i) # using integer color codes
legend("right", bty="n", lty=1,
col=1:length(k), legend=paste("k =",k))
5.2.2 Solutions for positive rate constants (growth)
If we consider plankton organisms like algae or bacteria as a constituent, their growth can be represented by a positive value of the rate constant \(k\). Let’s visualize the steady-state solutions for positive \(k\) using a few lines of R code (Figure 3). In this case, I kept also the residence time at fixed value.
c.in <- 1 # assumed inflow conc.
tau <- 2 # assumed residence time
k <- seq(0, 2, 0.01) # set of values for 'k'
plot(k, c.std(c.in, tau, k), type="l",
xlab="k", ylab="Steady state concentration")
Figure 3 reveals that Equation 5 obviously has a discontinuity. Looking a bit closer at the plot and the structure of Equation 5, it becomes clear that the critical value is \(k=\frac{1}{tau}\). At this point, the denominator of Equation 5 becomes zero resulting in an infinite result. Moreover, for any value of \(k\) greater than this threshold, the predicted steady states are implausible, since negative concentrations just don’t make sense. On the other hand, why should large positive values of \(k\) be suspicious in general? Couldn’t such values occur for fast-growing plankton organisms?
What the negative concentrations highlighted in red in Figure 3 tell us is the following: A steady state does not exist if \(k > 1/tau\). In other words, for any \(k > 1/tau\), the concentration in the reactor would rise continuously after an inoculation event. It does not converge to some finite steady-state value.
While negative concentrations are clearly impossible, close-to-infinite concentrations are obviously problematic either. Common sense tells us that, in reality, there must be a mechanism that keeps population density within feasible bounds. The key is an inevitable limitation of resources. If we consider the growth of plankton algae, potentially limiting resources would be, e.g., dissolved nutrients like P and N or the available light. Once the nutrients are depleted or light is reduced too much by self-shading, the organisms cannot grow any further.
We will deal with the topic of resource limitation in a separate session. Hence, the remainder of this section is only for the impatient who cannot wait to see how the concept of limitation resolves the issue of implausible steady states illustrated by Figure 3.
Let’s adopt a simple extension of Equation 4 which implements limitation in a very generic manner without introducing the resource as another state variable (Equation 6).
\[ \frac{d}{dt}c = \frac{1}{\tau} \cdot (c_{in} - c) + k \cdot c \cdot \left( 1 - \frac{c}{c_{max}}\right) \tag{6}\]
In Equation 6, the parameter \(c_{max}\) specifies a maximum population density also known as carrying capacity. As \(c \rightarrow c_{max}\), the last term approaches zero and the concentration cannot increase any further. Again, we can set the left hand side to zero to find the steady-state solution of Equation 6. A sequence of rearrangements leads to a quadratic expression (Equation 7).
\[ 0 = c^2 + \left( \frac{c_{max}}{k \cdot \tau} \right) \cdot c - \left( \frac{c_{max} \cdot c_{in}}{k \cdot \tau} \right) \tag{7}\]
Since Equation 7 is written in the so-called reduced form (coefficient of the \(c^2\) term being 1), we can apply the well-known quadratic formula to find the two candidate solutions. (Equation 8).
\[ x = -\frac{p}{2} \pm \sqrt{\left( \frac{p}{2} \right)^2 - q} \tag{8}\]
You probably remember from school times that \(p\) and \(q\) refer to the coefficients of the reduced quadratic equation like so: \(0 = x^2 + px + q\).
The code below visualizes the two candidate solutions with only one of them being valid (blue graph in Figure 4).
rm(list=ls())
c.in <- 1 # assumed inflow conc.
tau <- 2 # assumed residence time
k <- seq(0.01, 10, by=0.01) # set of positive 'k'
c.max <- 5 # assumed carrying capacity
p <- c.max / k / tau - c.max # coefficients of the ...
q <- -c.max * c.in / k / tau # ... reduced quadratic equation
c.1 <- -p/2 + sqrt((p/2)^2 - q) # the two candidate solutions
c.2 <- -p/2 - sqrt((p/2)^2 - q)
plot(range(k), range(-1, c.max), # basic plot
type="n", log="x",
xlab="k (log scaled)",
ylab="Steady state concentration")
lines(k, c.1, col="blue")
lines(k, c.2, col="red")
abline(h=c(0, c.in, c.max), lty=1:3) # some lines for orientation
legend(x=0.01, y=4, bty="n", title="Assumed settings",
legend=c(paste("tau =",tau), paste("c.in =",c.in), paste("c.max =",c.max))
)
legend(x=1, y=2, bty="n", lty=c(1,1,2,3),
col=c("blue","red",rep("black", 2)),
legend=c("Valid solution","Invalid solution","c = c.in","c = c.max")
)
5.3 Application to algae growth in a shallow lake
In June 2023, samples were taken from the inflow and outlet of the Radeburg I reservoir in Saxony, Germany. This very shallow reservoir can reasonably be regarded as polymictic. Sampling at the in- and outlet is facilitated through two bridges labeled S1 and S2 in Figure 5.
The following data were obtained in the field or looked up in online systems run by the state water authorities (Table 2).
| Item | Value | Unit | Data source |
|---|---|---|---|
| Chlorophyll-a at S1 | 4.0 | µg/L | Measured with fluorescence probe |
| Chlorophyll-a at S2 | 40.0 | µg/L | Measured with fluorescence probe |
| Water temperature S2 | 19.3 | °C | Measured |
| Inflow rate | 1.2 | m³/s | Upstream gage (see link in caption) |
| Storage volume | 360000.0 | m³ | Based on water level (see link in caption) |
According to the Chlorophyll-a data (Table 2), the concentration of suspended phytoplankton increased by factor 10 during the passage of the reservoir. Thus, phytoplankton obviously grew. If we assume that the system was in steady state, we can use the concept of the first-order reactive CFSTR, to learn a bit more about the situation. For that purpose, we start with Equation 4 but expand the first-order term for a moment into two separate terms (Equation 9).
\[ \frac{d}{dt}c = \frac{1}{\tau} \cdot (c_{in} - c) + g \cdot c - e \cdot c \tag{9}\]
This expanded notation helps us to distinguish growth, represented by the positive constant \(g\), from any losses of phytoplankton expressed by the negative constant \(e\). Let’s look at these two constants a bit closer:
\(g\) represents the actual growth rate of phytoplankton. Hence, it represents the efficiency of proliferation under given conditions which likely include a limitation by nutrients and possibly light.
\(e\) summarizes all processes resulting in an elimination of algae cells from the water column. Grazing by zooplankton and settling are likely the two most prominent loss processes.
Unfortunately, without more detailed and sophisticated measurements, we cannot disentangle \(g\) and \(e\). Hence, we are bound to looking at the net growth by collapsing growth and loss into a single constant again, which we will denote \(k\) for simplicity (Equation 10).
\[ \begin{aligned} \frac{d}{dt}c & = \frac{1}{\tau} \cdot (c_{in} - c) + (g-e) \cdot c \\ & = \frac{1}{\tau} \cdot (c_{in} - c) + k \cdot c \end{aligned} \tag{10}\]
Let’s formulate the steady state solution of Equation 10 by setting the left hand side to zero and solve for the constant \(k\) (Equation 11):
\[ k = \frac{1 - \frac{c_{in}}{c}}{\tau} \tag{11}\]
Inserting the Chlorophyll-a, volume, and stream flow data from Table 2 and accounting for the number of seconds per day (86400), we find the value of \(k\) in the unit of day-1. For comparison with the water residence time \(\tau\), we can express \(k\) as a doubling time (see Section 4).
tau <- 0.36e6 / (1.2*86400)
k <- (1 - 4/40) / tau
print(paste(
"tau:", round(tau, 2),
", k:", k,
", doubling time:", round(log(2)/k, 2)
))[1] "tau: 3.47 , k: 0.2592 , doubling time: 2.67"What can be learnt from these numbers? First of all, the numbers demonstrate that a phytoplankton cell needs to divide at least every 2.7 days if the algae density is to remain constant in spite of permanent losses (expressed by \(\tau\)).
However, this doubling time would only apply if the algae would not undergo any losses due to settling and grazing. Since this is clearly unrealistic, the actual doubling time is certainly shorter than 2.7 days.
Another interesting question in this context would be: How close is the observed value of \(k\) to the maximum intrinsic growth rate of plankton algae? A quick search brings up the work of König-Rinke, 2008 which provides a collection of intrinsic growth rates (Table 3).
| Algae group | Max. growth rate (day-1) | Optimum temperature (°C) | Susceptibility to grazing |
|---|---|---|---|
| Small green algae and diatoms | 2.40 | 25 | high |
| Large diatoms | 1.20 | 20 | low |
| Cyanobacteria (not N2-fixing) | 1.10 | 25 | very low |
| Cyanobacteria (N2-fixing) | 0.84 | 25 | very low |
According to Table 3, the observed net growth rate of about 0.26 day-1 at close to 20 °C (Table 2) is far below the maximum intrinsic growth rates of algae. We can conclude that …
either nutrients or light conditions are severely limiting the algae growth rates,
or algae undergo massive elimination due to predation, settling, and possibly diseases,
or both limitation and elimination occur simultaneously.
A look at detailed output of the fluourescence probe suggests that the increase in total chlorophyll was primarily due to the proliferation of green algae and diatoms (Figure 6).
5.4 Application to phosphorus concentrations in lakes
In the majority of freshwater systems, primary production is limited by the availability of phosphorus. Excessive phosphorus input to rivers and lakes is the main cause of anthropogenic eutrophication. At the same time, the control of phosphorus emissions is the key to water quality management. But much much phosphorus can be discharged into a lake until water quality deteriorates notably?
This question has been tackled by many researchers and the topic was in the center of limnological research since the 1970’s. Prominent comprehensive publications of that time were
Vollenweider, R.A. and Kerekes, J. (1982) Eutrophication of Waters: Monitoring, Assessment and Control. Report of the Cooperative Program on Eutrophication. Organization for the Economic Development and Cooperation, Paris. (found here)
Sas, H. (1990) Lake restoration by reduction of nutrient loading: Expectations, experiences, extrapolations. SIL Proceedings, 1922-2010, 24:1, 247-251. (found here; this article is a summary of a related book)
and a comprehensive summary can be found in the review of Brett & Benjamin, 2008. The main question researchers attempted to answer was:
What is the relation between the in-lake P concentration and the P input if all boundary conditions (including external P loads) are reasonably constant?
Obviously, this question addresses a steady-state condition and calls for a mass balance approach. To understand the main components of a lake phosphorus balance, consider Figure 7. Importantly, a lake is generally (but not necessarily at all times) a sink for phosphorus. This is due to the fact that organic and inorganic particles permanently undergo settling but the contained P is only partly released back into the water. The other part which is bound in slowly degradable organic matter or minerals is effectively removed from the system as it is trapped in deeper sediment layers during diagenesis.
How can we compile all components from Figure 7 into a model?
The best approach would probably be a model with multiple state variables to represent the dissolved and particulate, organic and inorganic fractions of phosphorus in the water body as well as in different sediment layers. While this is perfectly feasible, the model would contain many unknown parameters and require a numerical solution because of simultaneous equations. We’ll keep it a bit simpler for now.
Specifically, we consider total phosphorus (TP) in the water column as the only state variable and describe the loss of P in response to settling as a first order process (Equation 12). Note that we can use a term \(u/z\) instead of a plain rate constant \(k\) to conceptually attribute the loss to settling. In that expression, \(u\) denotes the effective settling velocity (m/year) and \(z\) is the mean lake depth (m).
\[ \begin{aligned} \frac{d}{dt}TP_{lake} &= \frac{1}{\tau} \cdot (TP_{in} - TP_{lake}) - k \cdot TP_{lake} \\ &= \frac{1}{\tau} \cdot (TP_{in} - TP_{lake}) - \frac{u}{z} \cdot TP_{lake} \end{aligned} \tag{12}\]
By setting the left hand side of Equation 12 to zero, we obtain the steady state solution (Equation 13):
\[ TP_{lake} = \frac{TP_{in}}{1 + \frac{u}{z} \cdot \tau} \tag{13}\]
Note the different sign in the denominator of Equation 13 compared to Equation 5. This is because the process of settling was explicitly equipped with a negative sign in Equation 12 to let \(u\) be positive.
Let’s check whether Equation 13 reasonably explains TP concentrations in a set of lakes in North-East Germany. With the R code below, we import the data set and compute a few derived variables such as the residence time, the mean lake depth, and the apparent input concentration from given variables.
rm(list=ls())
lakes <- read.table(sep="", header=T, text='
# Water & phosphorus balance data for lakes in NE Germany.
# Values of variable TP_lake represent annual averages.
# For original data sources, see Table 5 in Grüneberg et al. (2023)
# Phosphor-Retentionsmodelle für pH-neutrale Tagebauseen (DBU, Endbericht).
# https://www.dbu.de/OPAC/ab/DBU-Abschlussbericht-AZ-34919_02-Hauptbericht.pdf
# units m2 m3 m3/a g/a g/m3
lake area volume inflow TP_load TP_lake
Ploetze. 77300 277000 7.91E+04 2.84E+04 8.50E-02
G.Glien. 670000 4530000 2.80E+05 5.69E+04 2.27E-02
Scharm. 12100000 108000000 1.04E+07 2.06E+06 5.30E-02
Langer 1300000 2840000 3.25E+07 3.90E+06 6.77E-02
Mueggel. 7660000 36500000 1.39E+08 9.62E+06 6.87E-02
Neuendf. 2960000 6000000 3.00E+08 1.93E+07 5.64E-02
Stechl. 4250000 96900000 1.14E+06 1.76E+05 1.40E-02
Nehmitz. 1710000 10900000 1.11E+06 9.70E+04 1.90E-02
')
# append residence time and inflow concentration as derived variables
lakes <- cbind(lakes,
tau = with(lakes, volume / inflow),
TP_in = with(lakes, TP_load / inflow),
z = with(lakes, volume / area) # mean depth
)
plot(TP_lake ~ TP_in, data=lakes, xlab="TP in inflow (g/m3)",
ylab="TP in lake (g/m3)")
When inspecting the newly created column z, you will realize that the data set includes both shallow and deep water bodies (Figure 9) some of which will certainly undergo seasonal stratification. However, we are looking at the lakes from a long-term perspective and consider all of them as mixed systems, independent of mixing frequency.
barplot(-lakes[,"z"], names.arg=lakes[,"lake"],
ylab="Mean depth (m)", las=2)
Looking at Equation 13, we realize that all variables but the apparent settling velocity \(u\) are specified through data frame lakes. The value of \(u\) thus needs to be estimated from data. This process is generally known as parameter fitting and falls into the class of (non-linear) optimization problems. Because we estimate a single parameter only, it is called a one-dimensional problem.
To get familiar with the concept of fitting, we take the intuitive way and simply try how well the model performs for different guesses of \(u\). As a measure of model performance, we use the root-mean-square error (RMSE). The larger the RMSE, the worse is the fit.
# Function to predict TP_lake for given values, including the unknown u (m/a)
TP_lake_est <- function(u, z, cin, tau) { cin / (1 + u/z * tau) }
# Table with guessed values of the apparent settling velocity u for which we
# want to inspect the performance; the other column takes is for the results
test <- data.frame(
u = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50),
rmse = NA
)
# Generic function to compute the RMSE
RMSE <- function(obs, mod) { sqrt(mean((obs - mod)^2)) }
# Measure and plot the model performance for each guessed u
for (i in 1:nrow(test)) {
mod <- with(lakes, TP_lake_est(u=test[i,"u"], z=z, cin=TP_in, tau=tau))
test[i,"rmse"] <- RMSE(obs=lakes[,"TP_lake"], mod=mod)
}
plot(rmse ~ u, data=test, type="b", xlab="Tested u", ylab="RMSE (low is best)")
A function like the RMSE, which measures the performance in relation to the (unknown) parameters, is generally called an objective function. The goal of sophisticated optimization algorithms is to find a (set of) parameters such that the objective function is minimized. In R you will find different built-in algorithms for non-linear optimization like
optimize: handles one-dimensional problems only (like this one)optim: for estimating multiple parameters simultaneouslynls: likeoptimbut more convenient for simple models
The code below uses the latter option to pinpoint the optimum value of parameter \(u\) for the given example.
# find value of 'u' that provides the best fit
fit <- nls(formula=TP_lake ~ TP_in / (1 + u/z * tau), data=lakes, start=c(u=0))
u <- coef(fit)["u"]
# append column with modeled concentration
lakes$TP_lake_mod <- with(lakes, TP_in / (1 + u/z * tau))
# plot fitted model
plot(TP_lake_mod ~ TP_lake, data=lakes, pch=20, asp=1,
xlab="TP observed (mg/L)", ylab="TP modeled (mg/L)")
abline(0,1)
# add a measure of the goodness-of-fit (R squared)
R2 <- function (obs, sim) {round( 1 - var(sim-obs) / var(obs) ,2) }
legend("topleft", bty="n",
legend=paste0("R² = ",with(lakes, R2(obs=TP_lake, sim=TP_lake_mod))),
)
According to Figure 11, the CFSTR model with a first order loss term explains the variability among observed TP concentrations pretty well. A stronger mismatch is observed for only one out of the eight lakes. Hence, the long-term level of total phosphorus in a lake is apparently indeed co-controlled by (1) the inflow concentration, (2) the residence time, and (3) depth.
For the sake of comparison, we also test a model where the square root of the water residence time \(\sqrt{\tau}\) is employed as a predictor. This modification was proposed in one of Vollenweider’s earlier works and it is listed as hypothesis 4 in the Brett & Benjamin, 2008 paper.
It is important to note that the introduction of the square root or any other exponent around \(\tau\) conflicts with physical principles (as indicated by the units) and makes this a (semi)-empirical model.
# alternative model with parameter 'k' and the square root of tau
fit <- nls(formula=TP_lake ~ TP_in / (1 + k * sqrt(tau)), data=lakes, start=c(k=0))
k <- coef(fit)["k"]
# append another column with modeled concentration
lakes$TP_lake_mod2 <- with(lakes, TP_in / (1 + k * sqrt(tau)))
# compare predictions by the two models
with(lakes, {
plot(range(TP_lake), range(c(TP_lake_mod, TP_lake_mod2)),
type="n", asp=1, xlab="TP observed (mg/L)", ylab="TP modeled (mg/L)"
)
points(TP_lake, TP_lake_mod, pch=20)
points(TP_lake, TP_lake_mod2, pch=4)
legend("topleft", bty="n", pch=c(20,4),
legend=c(
paste0("First order model, R² = ",R2(obs=TP_lake, sim=TP_lake_mod)),
paste0("Vollenweider's m., R² = ",R2(obs=TP_lake, sim=TP_lake_mod2))
)
)
})
abline(0,1)
According to Figure 12, the model formulation using \(\sqrt{\tau}\) comes with a slightly higher value of R squared. However, the gain in explanatory power appears to be rather marginal in this case and comes at the price of the missing mechanistic basis.