In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises
Problem
Consider a system with chemical species $A$, $B$ and $AB$ (think, for example, of $A=NH_3$, $B=H^+$ and $AB=NH_4^+$) that undergo the following reactions:
- species $AB$ can dissociate into $A$ and $B$ at a rate controlled by the rate constant $k_f$: $$\mathrm{forward\ reaction}!!:\qquad AB \xrightarrow{k_f} A + B $$
- species $A$ and $B$ can react to form species $AB$ at a rate controlled by the rate constant $k_b$: $$\mathrm{backward\ reaction}!!:\qquad AB \xleftarrow{k_b} A + B $$
- species $A$ is net removed from the system at a rate controlled by the rate constant $\lambda$: $$\mathrm{net\ removal\ process}!!:\qquad A\xrightarrow{\lambda} \dots$$
In class you have seen that if the reversible process
$$AB \mathrel{\mathop{\rightleftarrows}^{{k_f}}_{{k_b}}} A + B$$
is vastly faster than the net removal process, then the evolution of the concentrations of $A$, $B$ and $AB$ can be modeled by assuming that on the time-scale of the slow removal process, the species $A$, $B$ and $AB$ are in an equilibrium (the so-called local equilibrium assumption).
This document will showcase how such local equilibrium chemical models can be implemented in R.
The full model
First, we model the dynamics of the chemical species $A$, $B$ and $AB$ without making the local equilibrium assumption. Thus, the state variables of the model are concentrations of the three species: $[A]$, $[B]$ and $[AB]$.
Mass balances
As shown in class, the mass balance equations for the concentrations of $A$, $B$ and $AB$ are
\begin{equation}
\begin{array}{rcl}
\displaystyle
\frac{d[A]}{dt} & = & R_f - R_b - \lambda \cdot [A], \[5mm]
\displaystyle
\frac{d[B]}{dt} & = & R_f -R_b, \[5mm]
\displaystyle
\frac{d[AB]}{dt} & = & - R_f + R_b,
\end{array}
\end{equation}
where the rates for the forward and backward reactions are
\begin{equation}\label{rates}
\begin{array}{rl}
\mathrm{forward\ reaction}!!: & R_f = k_f \cdot [AB], \[4mm]
\mathrm{backward\ reaction}!!: & R_b = k_b \cdot [A] \cdot [B].
\end{array}
\end{equation}
If the rate of the slow removal process is zero ($\lambda=0$), the system will reach an equilibrium. By definition, equilibrium corresponds to the situation where state variables do not change over time, i.e., their time derivatives are equal to zero. This occurs when the rates of the forward and backward reactions are equal, i.e., $R_f=R_b$. Using the rate expressions above, this equality yields the following relationship between the equilibrium concentrations of the involved species and the rate constants:
\begin{equation}
K_{eq}\mathop{\equiv}^\mathrm{def} \frac{k_f}{k_b} = \frac{[A]{eq}\cdot[B]{eq}}{[AB]{eq}}.
\end{equation}
In this equation, $K{eq}$ denotes the equilibrium constant, which is defined as the ratio of the rate constants for the forward and backward reactions. Equation 2 states that when the system reaches an equilibrium, the composition of the system is such that the reaction quotient, defined according to
\begin{equation}
Q= \frac{[A]\cdot[B]}{[AB]},
\end{equation}
is equal to the equilibrium constant $K_{eq}$. Thus, the reaction quotient $Q$ is a measure of how far from an equilibrium the system is.
Model implementation
To solve this model numerically, we will assume the following parameter values:
| parameter | Value | Unit |
|:-------------:|:--------------:|:--------------:|
| $k_f$ | 1000 | $s^{-1}$ |
| $k_b$ | $2 \times 10^6$| $(mol~m^{-3})^{-1}~s^{-1}$ |
| $\lambda$ | 0.1 | $s^{-1}$ |
These parameters imply the equilibrium constant of $K_{eq}=k_f/k_b = 0.5\times 10^{-3}~mol~m^{-3}$.
Additionally, we will assume that the total concentration of species $A$ in the system, denoted as $[A_{tot}]$ and equal to the sum $[A] + [AB]$, is $0.070~mol~m^{-3}$, the total concentration of species $B$, denoted as $[B_{tot}]$ and equal to the sum $[B] + [AB]$, is $0.060~mol~m^{-3}$, and the initial concentration of species (state variable) $B$ is $0.001~mol~m^{-3}$. Based on this information, the initial concentrations of the species (state variables) $AB$ and $A$ can easily be calculated according to $[AB]{ini} = [B{tot}]{ini} - [B]{ini}$ and $[A]{ini} = [A{tot}]{ini} - [AB]{ini}$.
Using this setup, we make a model that considers explicitly both the rapid equilibration process and the slower net removal process. That is, we model the dynamics of species $A$, $B$ and $AB$ based on the mass balance equations (1). We also calculate the total concentrations of elements A and B ($A_{tot}$ and $B_{tot}$), and the reaction quotient $Q=[A]\cdot[B]/[AB]$, as output variables. As noted above, the quotient is useful because it tells us how far from equilibrium the system is, where by "system" we mean the set of state variables $A$, $B$ and $AB$. We run the model for 1 millisecond and 10 seconds to illustrate what happens on the time scales of the fast equilibration and slow removal processes.
R-implementation
We need the R-package deSolve to solve this model:
library(deSolve)
The parameter values are in a vector called parms. The elements in this vector have a name; we will be able to use this name in the model function.
pars <- c (kf = 1000, # [/s]
kb = 2e6, # [/mol m3 /s]
lambda = 0.1) # [/s]
The initial conditions are calculated based on the values given above. The names of the state variables in the initial condition vector state.ini.full also determine the names that we can use in the model function.
Atot.ini <- 0.070 # [mol/m3] Initial concentration of total A = A + AB
Btot.ini <- 0.060 # [mol/m3] Initial concentration of total B = B + AB
B.ini <- 0.001 # [mol/m3] Initial concentration of B (STATE VARIABLE)
# initial concentrations of the other STATE VARIABLES
AB.ini <- Btot.ini - B.ini # [mol/m3]
A.ini <- Atot.ini - AB.ini # [mol/m3]
# initial conditions of all state variables
state.ini.full <- c(A = A.ini, B = B.ini, AB = AB.ini)
Now we define the model function based on the differential equations (1).
FullModel <- function(t, state, parms) {
with (as.list(c(state, parms)), {
# rate expressions
# AB <-> A+B
ratef <- kf * AB # forward AB -> A+B
rateb <- kb * A*B # backward A+B -> AB
rateA <- lambda*A # slow removal, specific to A
# mass balances
dA <- ratef - rateb - rateA
dB <- ratef - rateb
dAB <- - ratef + rateb
return(list(c(dA, dB, dAB), # the time derivatives
Atot = A+AB, Btot = B+AB, # total A and B
Q = A*B/AB)) # quotient
})
}
Model runs
To run the model, we define the time sequence for which we want the output to be generated (times) and solve the model using the ode function. Finally, we plot the dynamics of the state variables and of the three output variables.
First, we run the model for $1~ms$ to illustrate how fast the equilibrium is reached.
times <- seq(from=0, to=1e-3, length.out = 1000)
out.full <- ode(y=state.ini.full, times=times, func=FullModel, parms=pars)
plot(out.full, mfrow=c(2,3), ylab="mol/m3", xlab="time (s)")
We see that the system is initially in a disequilibrium, as indicated by the quotient $Q$ being different (lower) from the equilibrium constant $K_{eq}=0.5\times 10^{-3}~mol~m^{-3}$. However, we see that $Q$ reaches the value of $K_{eq}$ in about $0.2~ms$; this is the time scale on which the equilibrium between $A$, $B$ and $AB$ (due to the fast reversible reaction $AB\rightleftarrows A+B$) is reached. We also see that on this time scale the effect of the slow removal process is negligible (i.e., $[A]$ are practically does not change after reaching equilibrium).
Next, we run the model for $10~s$ to illustrate the effects on the time-scale of the slow removal process.
times <- seq(from=0, to=10, length.out=1000)
out.full <- ode(y=state.ini.full, times=times, func=FullModel, parms=pars)
plot(out.full, mfrow=c(2,3), ylab="mol/m3", xlab="time (s)")
As we can see from the quotient, on this time scale the species $A$, $B$ and $AB$ can be considered to be always in equilibrium.
Model based on the local equilibrium approximation
Now we model the system assuming a local equilibrium of the fast reversible reactions. As explained in class, this is done by considering different state variables, namely the total A ($A_{tot} = A+AB$) and total B ($B_{tot} = B+AB$), and by assuming that the species $A$, $B$ and $AB$ are always in an equilibrium, meaning that at any time the following relationship is valid (compare with equations 3 and 4):
\begin{equation}
\frac{[A]\cdot[B]}{[AB]} = K_{eq}.
\end{equation}
Mass balance equations and additional relationships
As explained in class, the local equilibrium condition implies the following mass balance equations for the new state variables:
\begin{equation}
\begin{array}{rcl}
\displaystyle
\frac{d[A_{tot}]}{dt} & = & -\lambda \cdot [A] \[5mm]
\displaystyle
\frac{d[B_{tot}]}{dt} & = & 0
\end{array}
\end{equation}
Additionally, it implies the following relationships between the new and the original state variables:
\begin{equation}
[B] + \frac{[B]}{K_{eq} + [B]} \cdot [A_{tot}] = [B_{tot}]
\end{equation}
\begin{equation}
[A] = \frac{K_{eq}}{K_{eq} + [B]} \cdot [A_{tot}]
\end{equation}
\begin{equation}
[AB] = \frac{[B]}{K_{eq} + [B]} \cdot [A_{tot}]
\end{equation}
\footnotesize
\textbf{Note:} You can use the following steps to derive the relationships 7--9. Using $K_{eq} = \frac{[A] \cdot [B]}{[AB]}$, we derive $[AB] = \frac{[A] \cdot [B]}{K_{eq}}$. Based on the definition of $A_{tot}$, we therefore have $[A_{tot}] = [A] + [AB] = [A] \cdot \left(1 +\frac{[B]}{K_{eq}}\right)$, which yields the expression $[A] = \frac{K_{eq}}{K_{eq} + [B]} \cdot [A_{tot}]$, which is Eq. 8. Substituting this expression back to $[AB] = \frac{[A] \cdot [B]}{K_{eq}}$, we obtain Eq. 9. Finally, substituting this result to the definition of $B_{tot}$ yields Eq. 7.
\normalsize
It is important to note the meaning of equations 6--9. Equation 7 shows how to calculate the concentration of $B$ for any given concentration of $A_{tot}$ and $B_{tot}$. Once this concentration is known, the concentrations of $A$ and $AB$ are calculated using equations 8 and 9, respectively. Then, using the concentration of $A$, the rate of change of $A_{tot}$ can be calculated from equation 6. These steps outline the approach of how to model the system under the local equilibrium assumption.
Solving for the equilibrium concentration of B
As explained above, finding the concentration of $B$ is the critical step in this approach. This is done by solving for $[B]$ in equation 7. This can be done in two ways: analytically or numerically.
An analytical approach
First, because the equation for finding $[B]$ is not that complicated, it is possible to find an analytical solution. Indeed, rearranging equation 7 yields a quadratic equation from which $[B]$ can readily be calculated. As shown in class, the solution is
\begin{equation}
[B] = \frac{1}{2}\left( -\beta + \sqrt{\beta^2 + 4 K_{eq}\cdot [B_{tot}]} \right),
\end{equation}
where
$$\beta = K_{eq} + [A_{tot}] - [B_{tot}].$$
A numerical approach
For more complex equilibria, such as those involving more than three chemical species (for example, carbonate equilibrium in water involves four species: $H_2CO_3$, $HCO_3^-$, $CO_3^{2-}$, and $H^+$) finding an analytical solution may be too difficult, or even impossible. In this case, numerical approximation can be used. For this particular case, we rearrange equation 7 according to
\begin{equation}
f([B]) = [B] + \frac{[B]}{K_{eq} + [B]} \cdot [A_{tot}] - [B_{tot}] = 0,
\end{equation}
which shows that $[B]$ can be found by finding the root (note: root of a function $f(x)$ is a value of $x$ where $f(x)=0$) of the function $f([B])$ defined on the left-hand side of this equation. In R, this is done using the uniroot function, as shown below.
R-implementation
Now we are ready to implement the model in R. We illustrate both approaches, one based on the analytical solution of $[B]$ and the other based on the numerical solution of $[B]$.
An analytical approach
Implementing this approach in R is easier because we have done the necessary work already by solving equation 7 for $[B]$. In addition to the time-derivatives of the state variables $A_{tot}$ and $B_{tot}$, we also return the values of $A$, $B$, $AB$, and $Q$ as the output.
EqModel_an <- function(t, state, parms) {
with (as.list(c(state, parms)), {
# equilibrium constant
Keq <- kf/kb
# calculate B from Atot and Btot (eq. 10)
beta <- Keq + Atot - Btot
B <- 0.5 * (-beta + sqrt(beta^2 + 4*Keq*Btot))
# calculate A from Atot and B (eq. 8)
A <- Keq /(Keq+B)*Atot
# mass balance equations for Atot and Btot (eq. 6)
dAtot <- -lambda * A
dBtot <- 0
return(list(c(dAtot, dBtot),
B = B, A = A, AB = Atot-A,
Q = A*B/(Atot-A) ))
})
}
A numerical approach
To implement this approach, we first define a function solveB for finding $[B]$ at equilibrium, i.e., the root of equation 11. Note that, within the scope of the function solveB, another function is defined (rootFun), which implements the function of $f([B])$ in equation 11. The root of this function is then found using uniroot. We need to input suitable lower and upper boundaries between which the root is sought. We choose 0 and $B_{tot}$, since $[B]$ cannot be negative or larger than $[B_{tot}]$.
solveB <- function(Keq, Atot, Btot){
# function whose root has to be sought (eq. 11)
rootFun <- function(B) {
return( B + B/(Keq+B)*Atot - Btot )
}
# uniroot will find the root; it returns a list with $root being the solution
r <- uniroot(f = rootFun, lower = 0, upper = Btot, tol = 1e-20)
return( r$root )
}
The remaining steps are similar as in the analytical approach. The only difference is that we use the above function solveB instead of calculating $B$ analytically.
EqModel_num <- function(t, state, parms) {
with (as.list(c(state, parms)), {
Keq <- kf/kb
# calculate B from Atot and Btot numerically, using solveB
B <- solveB(Keq = Keq, Atot = Atot, Btot = Btot)
# calculate A from Atot and B (eq. 8)
A <- Keq /(Keq+B)*Atot
# mass balance equations for Atot and Btot (eq. 6)
dAtot <- -lambda * A
dBtot <- 0
return(list(c(dAtot, dBtot),
B = B, A = A, AB = Atot-A,
Q = A*B/(Atot-A) ))
})
}
Finally, we run both models for 10 seconds.
state.ini.eq <- c(Atot = Atot.ini, Btot = Btot.ini)
out.eq_an <- ode(y=state.ini.eq, times=times, func=EqModel_an, parms=pars)
out.eq_num <- ode(y=state.ini.eq, times=times, func=EqModel_num, parms=pars)
Comparison of all models
We output the dynamics of the full model (black), alongside the equilibrium formulation (red and green). To optimize the plotting code, we introduce a function plot_AB that plots a specific variable ($A$, $B$, etc.) for all models in one graph. This function illustrates the usefulness of having named columns in matrices, as we do not need to worry about the different order of state variables in columns of the outputs out.full, out.eq_an and out.eq_num.
# plotting function
plot_AB <- function(varname, out1, out2, out3, legendloc="topright") {
plot (out1[,"time"], out1[,varname], type="l", main=varname,
ylab="mol/m3", xlab="time (s)", col=1, lty=1, lwd=1)
lines(out2[,"time"], out2[,varname], col=2, lty=2, lwd=2)
lines(out3[,"time"], out3[,varname], col=3, lty=3, lwd=4)
legend(legendloc, legend=c("Full", "Eq_an", "Eq_num"),
col=c(1,2,3), lty=c(1,2,3), lwd=c(1,2,4))
}
# plot variables Atot, Btot, Q, A, B and AB in separate graphs
par(mfrow = c(2,3))
plot_AB("Atot", out.full, out.eq_an, out.eq_num)
plot_AB("Btot", out.full, out.eq_an, out.eq_num)
plot_AB("Q", out.full, out.eq_an, out.eq_num, legendloc="bottomright")
plot_AB("A", out.full, out.eq_an, out.eq_num)
plot_AB("B", out.full, out.eq_an, out.eq_num, legendloc="bottomright")
plot_AB("AB", out.full, out.eq_an, out.eq_num)
The graphs show that in the full model, the initial concentrations of $A$, $B$ and $AB$ were not at equilibrium (the quotient $Q$ calculated from equation 4 is not equal to $K_{eq}$). However, the equilibrium concentration is rapidly reached, as indicated by the quotient $Q$ rapidly reaching $K_{eq}$. The full and approximate models yield indistinguishable results once the equilibrium has been reached.
Benchmarking
In class, it was mentioned that if we use the approximate model instead of the full model (because we are not interested in, or can ignore, the fast equilibration process), we can save considerable computational resources. To see the effects in this particular case, we compare the speed of solution for the three implementations: the full model, and the analytical and numerical versions of the local equilibrium model.
We run each model 100 times, so that random effets due to the computer being otherwise engaged are filtered out. The time output by R-function system.time is in seconds.
times <- seq(from=0, to=10, length.out=100)
cat("full model: \n")
print(system.time(
for (i in 1:100)
out.full <- ode(y=state.ini.full, times=times, func=FullModel, parms=pars)
)/100)
cat("analytical: \n")
print(system.time(
for (i in 1:100)
out.eq_an <- ode(y=state.ini.eq, times=times, func=EqModel_an, parms=pars)
)/100)
cat("numerical: \n")
print(system.time(
for (i in 1:100)
out.eq_num <- ode(y=state.ini.eq, times=times, func=EqModel_num, parms=pars)
)/100)
Two notes:
- The fact that the full model is so efficiently solved is because the default deSolve integrator (ode) is particularly well suited to solve so-called "stiff" problems (i.e., problems with largely different time scales). When using a traditional integration routine (Euler or Runge-Kutta), the solution would take several orders of magnitude longer.
- The apparant inefficiency of the numerically solved equilibrium model is due to the fact that uniroot is not the most efficient solver for this problem.
Application to ammonia degassing
Here, you will check your understanding of the concepts introduced above by modeling the effect of ammonia degassing on water $pH$. A possible environmental scenario could be a lake polluted with a spill of ammonia from an industrial source.
Task
Consider a water body with an initial $pH$ of 8 ($pH = -log([H^+])$, where $[H^+]$ is the concentration of protons in $mol~L^{-1}$!), and an initial concentration of total dissolved ammonia ($NH_x = NH_3 + NH_4^+$) of $1~mol~m^{-3}$. Write a model that will predict the change of water $pH$, along with the change of [$NH_x$], as a function of time due to ammonia ($NH_3$) degassing. Make the following assumptions in your model:
- The reversible (equilibration) process corresponds to the dissociation of ammonium: $$NH_4^+ \rightleftarrows NH_3 + H^+$$
- The qquilibrium constant for this process is approximately $K_{eq} = 5.38\times 10^{-10}\, mol \, L^{-1}$.
- Dissolved ammonia species $NH_3$ is degassing according to the first-order kinetics. Assume the rate constant of $1~d^{-1}$.
- The equilibration process is much faster than ammonia degassing. Thus, you can use the local equilibrium assumption to solve the model.
- For simplicity, assume that ammonia degassing is the only process that affects the water $pH$. In aquatic environments, water dissociation ($H_2O \rightleftarrows H^+ + OH^-$) also affects the $pH$, but you can ignore this process in this exercise. Additionally, there are other chemical species, such as dissolved CO$_2$, which strongly counter-act $pH$ changes induced by ammonia degassing. You can assume in this exercise that these species are not present in the water. These processes will be included in other exercises, where a more realistic (and also more complicated) model will be developed.
You can start with the R-markdown template file RTM_0D.Rmd to implement this model. You can obtain this file from Rstudio: File $\rightarrow$ new File $\rightarrow$ Rmarkdown $\rightarrow$ from template $\rightarrow$ RTM_0D. Save this file under a different name. Do not forget to change the heading of this file.
dynamic-R/RTM documentation built on Feb. 28, 2025, 1:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Problem
Consider a system with chemical species $A$, $B$ and $AB$ (think, for example, of $A=NH_3$, $B=H^+$ and $AB=NH_4^+$) that undergo the following reactions:
- species $AB$ can dissociate into $A$ and $B$ at a rate controlled by the rate constant $k_f$: $$\mathrm{forward\ reaction}!!:\qquad AB \xrightarrow{k_f} A + B $$
- species $A$ and $B$ can react to form species $AB$ at a rate controlled by the rate constant $k_b$: $$\mathrm{backward\ reaction}!!:\qquad AB \xleftarrow{k_b} A + B $$
- species $A$ is net removed from the system at a rate controlled by the rate constant $\lambda$: $$\mathrm{net\ removal\ process}!!:\qquad A\xrightarrow{\lambda} \dots$$
In class you have seen that if the reversible process $$AB \mathrel{\mathop{\rightleftarrows}^{{k_f}}_{{k_b}}} A + B$$
is vastly faster than the net removal process, then the evolution of the concentrations of $A$, $B$ and $AB$ can be modeled by assuming that on the time-scale of the slow removal process, the species $A$, $B$ and $AB$ are in an equilibrium (the so-called local equilibrium assumption).
This document will showcase how such local equilibrium chemical models can be implemented in R.
The full model
First, we model the dynamics of the chemical species $A$, $B$ and $AB$ without making the local equilibrium assumption. Thus, the state variables of the model are concentrations of the three species: $[A]$, $[B]$ and $[AB]$.
Mass balances
As shown in class, the mass balance equations for the concentrations of $A$, $B$ and $AB$ are \begin{equation} \begin{array}{rcl} \displaystyle \frac{d[A]}{dt} & = & R_f - R_b - \lambda \cdot [A], \[5mm] \displaystyle \frac{d[B]}{dt} & = & R_f -R_b, \[5mm] \displaystyle \frac{d[AB]}{dt} & = & - R_f + R_b, \end{array} \end{equation} where the rates for the forward and backward reactions are \begin{equation}\label{rates} \begin{array}{rl} \mathrm{forward\ reaction}!!: & R_f = k_f \cdot [AB], \[4mm] \mathrm{backward\ reaction}!!: & R_b = k_b \cdot [A] \cdot [B]. \end{array} \end{equation}
If the rate of the slow removal process is zero ($\lambda=0$), the system will reach an equilibrium. By definition, equilibrium corresponds to the situation where state variables do not change over time, i.e., their time derivatives are equal to zero. This occurs when the rates of the forward and backward reactions are equal, i.e., $R_f=R_b$. Using the rate expressions above, this equality yields the following relationship between the equilibrium concentrations of the involved species and the rate constants: \begin{equation} K_{eq}\mathop{\equiv}^\mathrm{def} \frac{k_f}{k_b} = \frac{[A]{eq}\cdot[B]{eq}}{[AB]{eq}}. \end{equation} In this equation, $K{eq}$ denotes the equilibrium constant, which is defined as the ratio of the rate constants for the forward and backward reactions. Equation 2 states that when the system reaches an equilibrium, the composition of the system is such that the reaction quotient, defined according to \begin{equation} Q= \frac{[A]\cdot[B]}{[AB]}, \end{equation} is equal to the equilibrium constant $K_{eq}$. Thus, the reaction quotient $Q$ is a measure of how far from an equilibrium the system is.
Model implementation
To solve this model numerically, we will assume the following parameter values:
| parameter | Value | Unit | |:-------------:|:--------------:|:--------------:| | $k_f$ | 1000 | $s^{-1}$ | | $k_b$ | $2 \times 10^6$| $(mol~m^{-3})^{-1}~s^{-1}$ | | $\lambda$ | 0.1 | $s^{-1}$ |
These parameters imply the equilibrium constant of $K_{eq}=k_f/k_b = 0.5\times 10^{-3}~mol~m^{-3}$.
Additionally, we will assume that the total concentration of species $A$ in the system, denoted as $[A_{tot}]$ and equal to the sum $[A] + [AB]$, is $0.070~mol~m^{-3}$, the total concentration of species $B$, denoted as $[B_{tot}]$ and equal to the sum $[B] + [AB]$, is $0.060~mol~m^{-3}$, and the initial concentration of species (state variable) $B$ is $0.001~mol~m^{-3}$. Based on this information, the initial concentrations of the species (state variables) $AB$ and $A$ can easily be calculated according to $[AB]{ini} = [B{tot}]{ini} - [B]{ini}$ and $[A]{ini} = [A{tot}]{ini} - [AB]{ini}$.
Using this setup, we make a model that considers explicitly both the rapid equilibration process and the slower net removal process. That is, we model the dynamics of species $A$, $B$ and $AB$ based on the mass balance equations (1). We also calculate the total concentrations of elements A and B ($A_{tot}$ and $B_{tot}$), and the reaction quotient $Q=[A]\cdot[B]/[AB]$, as output variables. As noted above, the quotient is useful because it tells us how far from equilibrium the system is, where by "system" we mean the set of state variables $A$, $B$ and $AB$. We run the model for 1 millisecond and 10 seconds to illustrate what happens on the time scales of the fast equilibration and slow removal processes.
R-implementation
We need the R-package deSolve to solve this model:
library(deSolve)
The parameter values are in a vector called parms. The elements in this vector have a name; we will be able to use this name in the model function.
pars <- c (kf = 1000, # [/s] kb = 2e6, # [/mol m3 /s] lambda = 0.1) # [/s]
The initial conditions are calculated based on the values given above. The names of the state variables in the initial condition vector state.ini.full also determine the names that we can use in the model function.
Atot.ini <- 0.070 # [mol/m3] Initial concentration of total A = A + AB Btot.ini <- 0.060 # [mol/m3] Initial concentration of total B = B + AB B.ini <- 0.001 # [mol/m3] Initial concentration of B (STATE VARIABLE) # initial concentrations of the other STATE VARIABLES AB.ini <- Btot.ini - B.ini # [mol/m3] A.ini <- Atot.ini - AB.ini # [mol/m3] # initial conditions of all state variables state.ini.full <- c(A = A.ini, B = B.ini, AB = AB.ini)
Now we define the model function based on the differential equations (1).
FullModel <- function(t, state, parms) { with (as.list(c(state, parms)), { # rate expressions # AB <-> A+B ratef <- kf * AB # forward AB -> A+B rateb <- kb * A*B # backward A+B -> AB rateA <- lambda*A # slow removal, specific to A # mass balances dA <- ratef - rateb - rateA dB <- ratef - rateb dAB <- - ratef + rateb return(list(c(dA, dB, dAB), # the time derivatives Atot = A+AB, Btot = B+AB, # total A and B Q = A*B/AB)) # quotient }) }
Model runs
To run the model, we define the time sequence for which we want the output to be generated (times) and solve the model using the ode function. Finally, we plot the dynamics of the state variables and of the three output variables.
First, we run the model for $1~ms$ to illustrate how fast the equilibrium is reached.
times <- seq(from=0, to=1e-3, length.out = 1000) out.full <- ode(y=state.ini.full, times=times, func=FullModel, parms=pars) plot(out.full, mfrow=c(2,3), ylab="mol/m3", xlab="time (s)")
We see that the system is initially in a disequilibrium, as indicated by the quotient $Q$ being different (lower) from the equilibrium constant $K_{eq}=0.5\times 10^{-3}~mol~m^{-3}$. However, we see that $Q$ reaches the value of $K_{eq}$ in about $0.2~ms$; this is the time scale on which the equilibrium between $A$, $B$ and $AB$ (due to the fast reversible reaction $AB\rightleftarrows A+B$) is reached. We also see that on this time scale the effect of the slow removal process is negligible (i.e., $[A]$ are practically does not change after reaching equilibrium).
Next, we run the model for $10~s$ to illustrate the effects on the time-scale of the slow removal process.
times <- seq(from=0, to=10, length.out=1000) out.full <- ode(y=state.ini.full, times=times, func=FullModel, parms=pars) plot(out.full, mfrow=c(2,3), ylab="mol/m3", xlab="time (s)")
As we can see from the quotient, on this time scale the species $A$, $B$ and $AB$ can be considered to be always in equilibrium.
Model based on the local equilibrium approximation
Now we model the system assuming a local equilibrium of the fast reversible reactions. As explained in class, this is done by considering different state variables, namely the total A ($A_{tot} = A+AB$) and total B ($B_{tot} = B+AB$), and by assuming that the species $A$, $B$ and $AB$ are always in an equilibrium, meaning that at any time the following relationship is valid (compare with equations 3 and 4): \begin{equation} \frac{[A]\cdot[B]}{[AB]} = K_{eq}. \end{equation}
Mass balance equations and additional relationships
As explained in class, the local equilibrium condition implies the following mass balance equations for the new state variables: \begin{equation} \begin{array}{rcl} \displaystyle \frac{d[A_{tot}]}{dt} & = & -\lambda \cdot [A] \[5mm] \displaystyle \frac{d[B_{tot}]}{dt} & = & 0 \end{array} \end{equation} Additionally, it implies the following relationships between the new and the original state variables:
\begin{equation} [B] + \frac{[B]}{K_{eq} + [B]} \cdot [A_{tot}] = [B_{tot}] \end{equation}
\begin{equation} [A] = \frac{K_{eq}}{K_{eq} + [B]} \cdot [A_{tot}] \end{equation}
\begin{equation} [AB] = \frac{[B]}{K_{eq} + [B]} \cdot [A_{tot}] \end{equation}
\footnotesize \textbf{Note:} You can use the following steps to derive the relationships 7--9. Using $K_{eq} = \frac{[A] \cdot [B]}{[AB]}$, we derive $[AB] = \frac{[A] \cdot [B]}{K_{eq}}$. Based on the definition of $A_{tot}$, we therefore have $[A_{tot}] = [A] + [AB] = [A] \cdot \left(1 +\frac{[B]}{K_{eq}}\right)$, which yields the expression $[A] = \frac{K_{eq}}{K_{eq} + [B]} \cdot [A_{tot}]$, which is Eq. 8. Substituting this expression back to $[AB] = \frac{[A] \cdot [B]}{K_{eq}}$, we obtain Eq. 9. Finally, substituting this result to the definition of $B_{tot}$ yields Eq. 7. \normalsize
It is important to note the meaning of equations 6--9. Equation 7 shows how to calculate the concentration of $B$ for any given concentration of $A_{tot}$ and $B_{tot}$. Once this concentration is known, the concentrations of $A$ and $AB$ are calculated using equations 8 and 9, respectively. Then, using the concentration of $A$, the rate of change of $A_{tot}$ can be calculated from equation 6. These steps outline the approach of how to model the system under the local equilibrium assumption.
Solving for the equilibrium concentration of B
As explained above, finding the concentration of $B$ is the critical step in this approach. This is done by solving for $[B]$ in equation 7. This can be done in two ways: analytically or numerically.
An analytical approach
First, because the equation for finding $[B]$ is not that complicated, it is possible to find an analytical solution. Indeed, rearranging equation 7 yields a quadratic equation from which $[B]$ can readily be calculated. As shown in class, the solution is
\begin{equation} [B] = \frac{1}{2}\left( -\beta + \sqrt{\beta^2 + 4 K_{eq}\cdot [B_{tot}]} \right), \end{equation} where $$\beta = K_{eq} + [A_{tot}] - [B_{tot}].$$
A numerical approach
For more complex equilibria, such as those involving more than three chemical species (for example, carbonate equilibrium in water involves four species: $H_2CO_3$, $HCO_3^-$, $CO_3^{2-}$, and $H^+$) finding an analytical solution may be too difficult, or even impossible. In this case, numerical approximation can be used. For this particular case, we rearrange equation 7 according to
\begin{equation}
f([B]) = [B] + \frac{[B]}{K_{eq} + [B]} \cdot [A_{tot}] - [B_{tot}] = 0,
\end{equation}
which shows that $[B]$ can be found by finding the root (note: root of a function $f(x)$ is a value of $x$ where $f(x)=0$) of the function $f([B])$ defined on the left-hand side of this equation. In R, this is done using the uniroot function, as shown below.
R-implementation
Now we are ready to implement the model in R. We illustrate both approaches, one based on the analytical solution of $[B]$ and the other based on the numerical solution of $[B]$.
An analytical approach
Implementing this approach in R is easier because we have done the necessary work already by solving equation 7 for $[B]$. In addition to the time-derivatives of the state variables $A_{tot}$ and $B_{tot}$, we also return the values of $A$, $B$, $AB$, and $Q$ as the output.
EqModel_an <- function(t, state, parms) { with (as.list(c(state, parms)), { # equilibrium constant Keq <- kf/kb # calculate B from Atot and Btot (eq. 10) beta <- Keq + Atot - Btot B <- 0.5 * (-beta + sqrt(beta^2 + 4*Keq*Btot)) # calculate A from Atot and B (eq. 8) A <- Keq /(Keq+B)*Atot # mass balance equations for Atot and Btot (eq. 6) dAtot <- -lambda * A dBtot <- 0 return(list(c(dAtot, dBtot), B = B, A = A, AB = Atot-A, Q = A*B/(Atot-A) )) }) }
A numerical approach
To implement this approach, we first define a function solveB for finding $[B]$ at equilibrium, i.e., the root of equation 11. Note that, within the scope of the function solveB, another function is defined (rootFun), which implements the function of $f([B])$ in equation 11. The root of this function is then found using uniroot. We need to input suitable lower and upper boundaries between which the root is sought. We choose 0 and $B_{tot}$, since $[B]$ cannot be negative or larger than $[B_{tot}]$.
solveB <- function(Keq, Atot, Btot){ # function whose root has to be sought (eq. 11) rootFun <- function(B) { return( B + B/(Keq+B)*Atot - Btot ) } # uniroot will find the root; it returns a list with $root being the solution r <- uniroot(f = rootFun, lower = 0, upper = Btot, tol = 1e-20) return( r$root ) }
The remaining steps are similar as in the analytical approach. The only difference is that we use the above function solveB instead of calculating $B$ analytically.
EqModel_num <- function(t, state, parms) { with (as.list(c(state, parms)), { Keq <- kf/kb # calculate B from Atot and Btot numerically, using solveB B <- solveB(Keq = Keq, Atot = Atot, Btot = Btot) # calculate A from Atot and B (eq. 8) A <- Keq /(Keq+B)*Atot # mass balance equations for Atot and Btot (eq. 6) dAtot <- -lambda * A dBtot <- 0 return(list(c(dAtot, dBtot), B = B, A = A, AB = Atot-A, Q = A*B/(Atot-A) )) }) }
Finally, we run both models for 10 seconds.
state.ini.eq <- c(Atot = Atot.ini, Btot = Btot.ini) out.eq_an <- ode(y=state.ini.eq, times=times, func=EqModel_an, parms=pars) out.eq_num <- ode(y=state.ini.eq, times=times, func=EqModel_num, parms=pars)
Comparison of all models
We output the dynamics of the full model (black), alongside the equilibrium formulation (red and green). To optimize the plotting code, we introduce a function plot_AB that plots a specific variable ($A$, $B$, etc.) for all models in one graph. This function illustrates the usefulness of having named columns in matrices, as we do not need to worry about the different order of state variables in columns of the outputs out.full, out.eq_an and out.eq_num.
# plotting function plot_AB <- function(varname, out1, out2, out3, legendloc="topright") { plot (out1[,"time"], out1[,varname], type="l", main=varname, ylab="mol/m3", xlab="time (s)", col=1, lty=1, lwd=1) lines(out2[,"time"], out2[,varname], col=2, lty=2, lwd=2) lines(out3[,"time"], out3[,varname], col=3, lty=3, lwd=4) legend(legendloc, legend=c("Full", "Eq_an", "Eq_num"), col=c(1,2,3), lty=c(1,2,3), lwd=c(1,2,4)) } # plot variables Atot, Btot, Q, A, B and AB in separate graphs par(mfrow = c(2,3)) plot_AB("Atot", out.full, out.eq_an, out.eq_num) plot_AB("Btot", out.full, out.eq_an, out.eq_num) plot_AB("Q", out.full, out.eq_an, out.eq_num, legendloc="bottomright") plot_AB("A", out.full, out.eq_an, out.eq_num) plot_AB("B", out.full, out.eq_an, out.eq_num, legendloc="bottomright") plot_AB("AB", out.full, out.eq_an, out.eq_num)
The graphs show that in the full model, the initial concentrations of $A$, $B$ and $AB$ were not at equilibrium (the quotient $Q$ calculated from equation 4 is not equal to $K_{eq}$). However, the equilibrium concentration is rapidly reached, as indicated by the quotient $Q$ rapidly reaching $K_{eq}$. The full and approximate models yield indistinguishable results once the equilibrium has been reached.
Benchmarking
In class, it was mentioned that if we use the approximate model instead of the full model (because we are not interested in, or can ignore, the fast equilibration process), we can save considerable computational resources. To see the effects in this particular case, we compare the speed of solution for the three implementations: the full model, and the analytical and numerical versions of the local equilibrium model.
We run each model 100 times, so that random effets due to the computer being otherwise engaged are filtered out. The time output by R-function system.time is in seconds.
times <- seq(from=0, to=10, length.out=100) cat("full model: \n") print(system.time( for (i in 1:100) out.full <- ode(y=state.ini.full, times=times, func=FullModel, parms=pars) )/100) cat("analytical: \n") print(system.time( for (i in 1:100) out.eq_an <- ode(y=state.ini.eq, times=times, func=EqModel_an, parms=pars) )/100) cat("numerical: \n") print(system.time( for (i in 1:100) out.eq_num <- ode(y=state.ini.eq, times=times, func=EqModel_num, parms=pars) )/100)
Two notes:
- The fact that the full model is so efficiently solved is because the default deSolve integrator (ode) is particularly well suited to solve so-called "stiff" problems (i.e., problems with largely different time scales). When using a traditional integration routine (Euler or Runge-Kutta), the solution would take several orders of magnitude longer.
- The apparant inefficiency of the numerically solved equilibrium model is due to the fact that uniroot is not the most efficient solver for this problem.
Application to ammonia degassing
Here, you will check your understanding of the concepts introduced above by modeling the effect of ammonia degassing on water $pH$. A possible environmental scenario could be a lake polluted with a spill of ammonia from an industrial source.
Task
Consider a water body with an initial $pH$ of 8 ($pH = -log([H^+])$, where $[H^+]$ is the concentration of protons in $mol~L^{-1}$!), and an initial concentration of total dissolved ammonia ($NH_x = NH_3 + NH_4^+$) of $1~mol~m^{-3}$. Write a model that will predict the change of water $pH$, along with the change of [$NH_x$], as a function of time due to ammonia ($NH_3$) degassing. Make the following assumptions in your model:
- The reversible (equilibration) process corresponds to the dissociation of ammonium: $$NH_4^+ \rightleftarrows NH_3 + H^+$$
- The qquilibrium constant for this process is approximately $K_{eq} = 5.38\times 10^{-10}\, mol \, L^{-1}$.
- Dissolved ammonia species $NH_3$ is degassing according to the first-order kinetics. Assume the rate constant of $1~d^{-1}$.
- The equilibration process is much faster than ammonia degassing. Thus, you can use the local equilibrium assumption to solve the model.
- For simplicity, assume that ammonia degassing is the only process that affects the water $pH$. In aquatic environments, water dissociation ($H_2O \rightleftarrows H^+ + OH^-$) also affects the $pH$, but you can ignore this process in this exercise. Additionally, there are other chemical species, such as dissolved CO$_2$, which strongly counter-act $pH$ changes induced by ammonia degassing. You can assume in this exercise that these species are not present in the water. These processes will be included in other exercises, where a more realistic (and also more complicated) model will be developed.
You can start with the R-markdown template file RTM_0D.Rmd to implement this model. You can obtain this file from Rstudio: File $\rightarrow$ new File $\rightarrow$ Rmarkdown $\rightarrow$ from template $\rightarrow$ RTM_0D. Save this file under a different name. Do not forget to change the heading of this file.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.