In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises
Answer
Task 1: Mass balance equations of the carbonate system
Assuming that the forward and backward reactions are elementary, we can use the first-order kinetics to describe the rate laws for each reaction. Considering the stoichiometries of reactions 1--4, we obtain the following mass balance equations:
$$
\begin{array}{rclr}
\displaystyle
\frac{d[H_2CO_3]}{dt} & = & -R_1 + \color{red}{R_C} & \qquad (1a)\[5mm]
\displaystyle
\frac{d[HCO_3^-]}{dt} & = & R_1 - R_2 & \qquad (1b) \[5mm]
\displaystyle
\frac{d[CO_3^{2-}]}{dt} & = & R_2 & \qquad (1c)\[5mm]
\displaystyle
\frac{d[NH_4^+]}{dt} & = & -R_3 & \qquad (1d) \[5mm]
\displaystyle
\frac{d[NH_3]}{dt} & = & R_3 + \color{red}{R_N} & \qquad (1e) \[5mm]
\displaystyle
\frac{d[H_2O]}{dt} & = & -R_4 & \qquad (1f) \[5mm]
\displaystyle
\frac{d[OH^-]}{dt} & = & R_4 & \qquad (1g) \[5mm]
\displaystyle
\frac{d[H^+]}{dt} & = & R_1 + R_2 + R_3 + R_4 & \qquad (1h) \[5mm]
\end{array}
$$
where the individual rates ($R_1$, $R_2$, $R_3$ and $R_4$) are defined as the difference between the forward and backward reaction rates:
$$
\begin{array}{rclr}
R_1 &=& k_{1f} \cdot [H_2CO_3] - k_{1b} \cdot [HCO_3^-] \cdot [H^+] &
\qquad\hspace{2cm} (2a) \[4mm]
R_2 &=& k_{2f} \cdot [HCO_3^{-}] - k_{2b} \cdot [CO_3^{2-}] \cdot [H^+] &
\qquad\hspace{2cm} (2b) \[4mm]
R_3 &=& k_{3f} \cdot [NH_4^{+}] - k_{3b} \cdot [NH_3] \cdot [H^+] &
\qquad\hspace{2cm} (2c) \[4mm]
R_4 &=& k_{4f} \cdot [H_2O] - k_{4b} \cdot [OH^-] \cdot [H^+] &
\qquad\hspace{2cm} (2d)
\end{array}
$$
Note that the slow processes, $R_C$ and $R_N$, only affect $H_2CO_3$ and $NH_3$.
Task 2: Equilibrium conditions
In an equilibrium, the time derivative of each species in the system is zero. Based on equations 1a and 2a, we obtain $R_1=0$, which yields the relationship
$$
K_{1} = \frac{ [HCO_3^-] \cdot [H^+] }{ [H_2CO_3] }. \qquad\hspace{2cm} (3a)
$$
Similarly, based on Eq. 1c and 2b, we obtain $R_2=0$, which yields the relationship
$$
K_{2} = \frac{ [CO_3^{2-}] \cdot [H^+] }{[HCO_3^-]}.
\qquad\hspace{2cm} (3b)
$$
Based on Eq. 1d and 2c, we obtain $R_3=0$, which yields the relationship
$$
K_{n} = \frac{ [NH_3] \cdot [H^+] }{[NH_4^+]}.
\qquad\hspace{2cm} (3c)
$$
Finally, based on Eq. 1f and 2d, we obtain $R_4=0$, which yields the relationship
$$
K_{w}' = \frac{ [OH^-] \cdot [H^+] }{[H_2O]}.
\qquad\hspace{2cm} (3d)
$$
Here, we defined the dissociation constants $K_1$, $K_2$, $K_n$, $K_w'$ according to $K_1 = k_{1f}/k_{1b}$, $K_2 = k_{2f}/k_{2b}$, $K_n = k_{3f}/k_{3b}$, and $K_w' = k_{4f}/k_{4b}$. If these relationships are valid, then the time derivatives in equations 1b, 1e, 1g and 1h are zero. This means that these four equations do not provide any extra information for the equilibrium concentrations.
Task 3: Lump-sum species
Here we combine equations 1a--1h and see what they imply for the suggested lump-sum species. By summing equations 1a, 1b and 1c, we obtain
$$
\frac{d[DIC]}{dt} = {\color{red}{R_C}} \qquad\hspace{2cm} (4a)
$$
Similarly, by summing equations 1d and 1e, we obtain
$$
\frac{d[NH_x]}{dt} = {\color{red}{R_N}} \qquad\hspace{2cm} (4b)
$$
Similarly, by summing equations 1f and 1g, we obtain
$$
\frac{d[H_xO]}{dt} = 0 \qquad\hspace{2cm} (4c)
$$
Finally, by summing equations according to (1b) + $2\times$(1c) + (1g) $-$ (1d) $-$ (1h), we obtain
$$
\frac{d[ALK]}{dt} = 0 \qquad\hspace{2cm} (4d)
$$
These equations show that DIC is added at a rate of the slow process, $R_C$, $NH_x$ is added at a rate of the slow process, $R_N$, while $H_xO$ and $ALK$ are conserved (i.e., not changing over time). Specifically, none of the lump-sum species $DIC$, $NH_x$, $H_xO$ and $ALK$ is affected by the fast reversible reactions.
Task 4: From lump-sum species to the original species
This task has been solved in parts I and II of this reader series. Here we only summarize the results:
$$
\begin{array}{rclr}
[HCO_3^-] &=& \displaystyle \frac{K_{1}[H^+]}{[H^+][H^+] + K_{1}[H^+] + K_{1} K_{2}}\cdot [DIC] &
\qquad\hspace{2cm} (5a) \[5mm]
[H_2CO_3] &=& \displaystyle \frac{[H^+][H^+]}{[H^+][H^+] + K_{1}[H^+] + K_{1} K_{2}}\cdot [DIC] &
\qquad\hspace{2cm} (5b) \[5mm]
[CO_3^{2-}] &=& \displaystyle \frac{K_1 K_2}{[H^+][H^+] + K_{1}[H^+] + K_{1} K_{2}}\cdot [DIC] &
\qquad\hspace{2cm} (5c) \[5mm]
[NH_3] &=& \displaystyle \frac{K_n}{K_n + [H^+]}\cdot [NH_x] & \qquad\hspace{2cm} (5d) \[5mm]
[NH_4^+] &=& \displaystyle \frac{[H^+]}{K_n + [H^+]} \cdot [NH_x] & \qquad\hspace{2cm} (5e) \[5mm]
[OH^-] &=& \displaystyle \frac{K_w'}{K_w' + [H^+]}\cdot [H_xO] & \qquad\hspace{2cm} (5f) \[5mm]
[H_2O] &=& \displaystyle \frac{[H^+]}{K_w' + [H^+]}\cdot [H_xO] & \qquad\hspace{2cm} (5g)
\end{array}
$$
Once the concentrations $[HCO_3^-]$, $[CO_3^{2-}]$, $[NH_4^+]$ and $[OH^-]$ are calculated from equations 5a, 5c, 5e and 5f, respectively, the alkalinity is calculated from the definition:
$$
[ALK] = [HCO_3^-] + 2\cdot [CO_3^{2-}] + [OH^-] - [NH_4^+] - [H^+]
\qquad\hspace{2cm} (5h)
$$
Because $[HCO_3^-]$ and $[CO_3^{2-}]$ depend on $[DIC]$ and $[H^+]$ (Eqs. 5a and 5c), $[NH_4^+]$ depends on $[NH_x]$ and $[H^+]$ (Eq. 5d), and $[OH^-]$ depends on $[H_xO]$ and $[H^+]$ (Eq. 5f), equation 5h provides an intimate relationship between $[ALK]$, $[DIC]$, $[NH_x]$, $[H_xO]$ and $[H^+]$ when the carbonate, ammonium and water species are in equilibrium. Using this relationship, one can calculate any of the five species provided that the other \emph{four} species are known. Because equation 5h is highly non-linear, these calculations are done numerically.
Task 5: Application to organic matter mineralization
R implementation
First, we expand the function that calculates alkalinity by including the concentration of $NH_4^+$ (Eq. $5e$).
solveALK <- function(K1=4.44e-07, K2=4.67e-11, Kn=5.68e-10, Kw.prime=1.79e-16,
HxO=55.3, DIC=0, NHx=0, pH=7){
H <- 10^(-pH)
HCO3 <- K1*H / (H^2 + K1*H + K1*K2) * DIC # eq 5a
CO3 <- K1*K2 / (H^2 + K1*H + K1*K2) * DIC # eq 5c
NH4 <- H / (H + Kn) * NHx # eq 5e
OH <- Kw.prime / (Kw.prime + H) * HxO # eq 5f
ALK <- HCO3 + 2*CO3 + OH - NH4 - H # eq 5h
return(ALK)
}
Next, we update the function that calculates $[H^+]$ from given concentrations of $ALK$, $DIC$, $NH_x$ and $H_xO$ (Eq. $5h$). This is done by finding a root of a function defined as the left-hand-side minus the right-hand-side of Eq. $5h$.
solveH <- function(K1=4.44e-07, K2=4.67e-11, Kn=5.68e-10, Kw.prime=1.79e-16,
HxO=55.3, DIC=0, NHx=0, ALK=0){
# function whose root has to be sought
rootFun <- function(H) {
ALK.est <- solveALK(K1=K1, K2=K2, Kn=Kn, Kw.prime=Kw.prime,
HxO=HxO, DIC=DIC, NHx=NHx, pH=-log10(H))
return(ALK.est - ALK)
}
# uniroot will find the root
h <- uniroot(f = rootFun, lower = 0, upper = 1, tol = 1e-30)
return(h$root)
}
It is always good to test the functions. We see that the results are consistent with expectations.
pH.new <- -log10(solveH(DIC=2e-3, NHx=1e-5, ALK=2e-3))
ALK <- solveALK(DIC=2e-3, NHx=1e-5, pH = pH.new)
pH.neutral <- -log10(solveH(DIC=0, NHx=0, ALK=0))
c(pH.new=pH.new, ALK=ALK, pH.neutral=pH.neutral)
Next, we implement the model function based on the local equilibrium assumption. Note that
-
The equilibrium constants (dissociation constants) are in $mol~L^{-1}$, so all species must be in this unit, too!
-
$OM$ is added to the model as a state variable removed by the slow process.
-
Because we will additionally want to use the function for modelling ammonia degassing (Task 6), the effect of this process is added, too. Its rate constant ($\lambda$) will be set to 0 by default.
OMMpHmodel <- function(t, state, parms) {
with (as.list(c(state, parms)), {
# variables needed as output and for regulation of process rates
H <- solveH(K1=K1, K2=K2, Kn=Kn, Kw.prime=Kw.prime,
HxO=HxO, DIC=DIC, NHx=NHx, ALK=ALK)
HCO3 <- K1*H / (H^2 + K1*H + K1*K2) * DIC # eq 5a
CO3 <- K1*K2 / (H^2 + K1*H + K1*K2) * DIC # eq 5c
NH4 <- H / (H + Kn) * NHx # eq 5e
OH <- Kw.prime / (Kw.prime + H) * HxO # eq 5f
NH3 <- NHx - NH4 # from definition of NHx
H2CO3 <- DIC - HCO3 - CO3 # from definition of DIC
H2O <- HxO - OH # from definition of HxO
# process rates
RC <- rC*OM # OM degradation
Rdegas <- lambda * NH3 # ammonia degassing (specifically removes NH3)
# mass balance equations (assuming local equilibrium)
dOM <- -RC
dDIC <- RC # Eq. 4a
dNHx <- RC*NC - Rdegas # Eq. 4b - ammonia degassing
dHxO <- 0 # Eq. 4c
dALK <- 0 # Eq. 4d
return(list(c(dOM, dDIC, dNHx, dHxO, dALK),
pH = -log10(H),
H = H, OH = OH, H2O = H2O,
H2CO3 = H2CO3, HCO3 = HCO3, CO3 = CO3,
NH4 = NH4, NH3 = NH3))
})
}
We run the model for 10 days, using the initial conditions given in the task.
require(deSolve)
pars <- c(rC = 0.1, # [1/day] OM degradation
lambda = 0, # [1/day] ammonia degasing
NC = 16/106, # molN/molC
K1 = 4.44e-07, # mol/L
K2 = 4.67e-11, # mol/L
Kn = 5.68e-10, # mol/L
Kw.prime = 1.79e-16) # mol/L
pH.ini <- 7.9 # H.ini based on pH
OM.ini <- 100*1e-6 # mol/L
DIC.ini <- 2000*1e-6 # mol/L
NHx.ini <- 10*1e-6 # mol/L
HxO.ini <- 55.3 # mol/L
ALK.ini <- solveALK(HxO = HxO.ini, DIC = DIC.ini, NHx = NHx.ini, pH = pH.ini)
state.ini <- c(OM=OM.ini, DIC=DIC.ini, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini)
times <- seq(from=0, to=10, length.out=100)
out <- ode(y=state.ini, times=times, func=OMMpHmodel, parms=pars)
plot(out, mfrow = c(4,4), las = 1)
We see that mineralization of about $60~\mu mol~L^{-1}$ of organic matter over the period of 10 days leads to a decrease in $pH$ from 7.9 to about 7.63. As expected, $ALK$ and $H_xO$ are conserved and changes in the concentration of $H_2O$ are negligible due to vastly larger concentrations of water compared to all the other species in the system.
Task 6: Does the formulation of slow reversible processes matter?
If the stoichiometry of OM mineralization is represented by the reaction
$$
\quad CH_2O(NH_3)_{NC} + O_2 \rightarrow HCO_3^- + H^- + {NC} \cdot NH_3,
$$
the mass balance equations 1c--g do not change while equations 1a, 1b, and 1h change in the following way:
$$
\begin{array}{rclr}
\displaystyle
\frac{d[H_2CO_3]}{dt} & = & R_1 & \qquad (1a')\[5mm]
\displaystyle
\frac{d[HCO_3^-]}{dt} & = & -R_1 + R_2 + \color{red}{R_C}& \qquad (1b') \[5mm]
\displaystyle
\frac{d[H^+]}{dt} & = & -R_1 - R_2 - R_3 - R_4 + \color{red}{R_C} & \qquad (1h') \[5mm]
\end{array}
$$
This change, however, has no effect on the differential equations for the lump-sum variables $DIC$, $NH_x$, $H_xO$ and $ALK$, as can easily be verified by following the same derivation steps as those leading to equations 4a--d. This can also be shown more easily: although $R_C$ adds $HCO_3^-$ instead of $H_2CO_3$, these two species are both included in the $DIC$ pool, thus summing them up yields the same result for the time-derivative of $DIC$. Similarly, although the process increases $ALK$ at a rate $R_C$ due to the addition of $HCO_3^-$, it also decreases $ALK$ at the same rate $R_C$ due to the addition of $H^+$, as follows from the negative sign in front of $[H^+]$ in the definition of $ALK$.
This result has important implications: the specific "chemical" formulation of the slow reversible process has no effect on the dynamics of the lump-sum variables, as long as the chemical reactions are correctly balanced in terms of the elements and charges and the local equilibrium conditions can be assumed.
Task 7: Application to ammonia degassing
R implementation
First, we reproduce here the ammonia degassing model developed in Part I:
solveB <- function(Keq, Atot, Btot){
rootFun <- function(B) B + B/(Keq+B)*Atot - Btot
return(uniroot(f = rootFun, lower = 0, upper = Btot, tol = 1e-20)$root)
}
AmmoniaDegassing <- function(t, state, parms) {
with (as.list(c(state, parms)), {
H <- solveB(Keq = Kn, Atot = NHx, Btot = Hx)
NH3 <- Kn/(Kn+H) * NHx
NH4 <- NHx - NH3
# mass balance equations for the lump-sum variables:
dNHx <- -lambda * NH3
dHx <- 0
return(list(c(dNHx, dHx),
H = H, NH3 = NH3, NH4 = NH4, # NH4 corresponds to [NH4+]
pH = -log10(H), # return also pH
Q = NH3*H/NH4 )) # return also Q
})
}
We define the initial conditions and run the model. In this step, we reuse the parameters defined for the OMMpHmodel model with the value of the degassing rate constant, $\lambda$, changed to $1~d^{-1}$.
# update model parameters
pars1 <- pars
pars1["lambda"] <- 1
# initial conditions
NHx.ini <- 1e-3 # mol/L
pH.ini <- 8
H.ini <- 10^(-pH.ini) # mol/L
Hx.ini <- H.ini + H.ini/(pars1[["Kn"]] + H.ini) * NHx.ini
state.ini1 <- c(NHx = NHx.ini, Hx = Hx.ini)
# run the model
out1 <- ode(y=state.ini1, times=times, func=AmmoniaDegassing, parms=pars1)
Next, we run the OMMpHmodel model with the same initial conditions. By setting the OM degradation rate constant to 0, and further setting the concentrations of $DIC$ and $H_xO$ to 0, we expect that the model becomes equivalent to the AmmoniaDegassing model.
# update model parameters
pars2 <- pars1
pars2["rC"] <- 0
# initial conditions
ALK.ini <- solveALK(HxO=0, DIC=0, NHx=NHx.ini, pH=pH.ini)
state.ini2 <- c(OM=OM.ini, DIC=0, NHx=NHx.ini, HxO=0, ALK=ALK.ini)
# run the model
out2 <- ode(y=state.ini2, times=times, func=OMMpHmodel, parms=pars2)
Thus, we expect to obtain the same results, at least for the state variables that are common for both models. This is indeed the case when we plot the results.
plot(out2, obs=out1, which=c("NHx", "NH4", "NH3", "pH"),
mfrow = c(1,4), lwd=2, lty=1, col=1,
obspar = list(lty=2,lwd=2,col=2,type="l"))
Now that we have verified that the model OMMpHmodel can be used for solving the ammonia degassing problem, too, we can illustrate the $pH$ buffering role of $DIC$ in natural waters. We will use the same initial conditions and the ammonia degassing rate constant as above and perform the calculations for $DIC$ concentrations of 0, 1, and $2~mmol~L^{-1}$.
# without DIC (but with ammonia and water!)
ALK.ini3 <- solveALK(DIC=0, HxO=HxO.ini, NHx=NHx.ini, pH=pH.ini)
state.ini3 <- c(OM=OM.ini, DIC=0, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini3)
out3 <- ode(y=state.ini3, times=times, func=OMMpHmodel, parms=pars2)
# model with DIC=1 mmol/L
ALK.ini4 <- solveALK(DIC=1e-3, HxO=HxO.ini, NHx=NHx.ini, pH=pH.ini)
state.ini4 <- c(OM=OM.ini, DIC=1e-3, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini4)
out4 <- ode(y=state.ini4, times=times, func=OMMpHmodel, parms=pars2)
# model with DIC=2 mmol/L
ALK.ini5 <- solveALK(DIC=2e-3, HxO=HxO.ini, NHx=NHx.ini, pH=pH.ini)
state.ini5 <- c(OM=OM.ini, DIC=2e-3, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini5)
out5 <- ode(y=state.ini5, times=times, func=OMMpHmodel, parms=pars2)
# plot results
plot(out3, out4, out5,
which=c("pH", "H", "OH", "HxO", "DIC", "ALK"),
mfrow = c(2,3), las = 1, col=1:3, lty=1:3, lwd=2)
legend("topright", legend = c("DIC=0", "DIC=1 mmol/L", "DIC=2 mmol/L"), bty="n",
col = 1:3, lwd = 2, lty = 1:3)
# plot results
plot(out3, out4, out5,
which=c("NHx", "NH4", "NH3", "H2CO3", "HCO3", "CO3"),
mfrow = c(2,3), las = 1, col=1:3, lty=1:3, lwd=2)
legend("topright", legend = c("DIC=0", "DIC=1 mmol/L", "DIC=2 mmol/L"), bty="n",
col = 1:3, lwd = 2, lty = 1:3)
The results show that when $DIC$ is present, the protons released due to ammonia degassing are scavenged by both $CO_3^{2-}$ and $HCO_3^-$ in the solution, leading to an increase in $H_2CO_3$. Consequently, the increase in the concentration of protons, and thus the decrease in $pH$, is smaller in the presence of $DIC$, illustrating the $pH$ buffering effect by the carbonate system. The magnitude of the effect does not, however, scale linearly with the $DIC$ concentration.
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.
Answer
Task 1: Mass balance equations of the carbonate system
Assuming that the forward and backward reactions are elementary, we can use the first-order kinetics to describe the rate laws for each reaction. Considering the stoichiometries of reactions 1--4, we obtain the following mass balance equations:
$$ \begin{array}{rclr} \displaystyle \frac{d[H_2CO_3]}{dt} & = & -R_1 + \color{red}{R_C} & \qquad (1a)\[5mm] \displaystyle \frac{d[HCO_3^-]}{dt} & = & R_1 - R_2 & \qquad (1b) \[5mm] \displaystyle \frac{d[CO_3^{2-}]}{dt} & = & R_2 & \qquad (1c)\[5mm] \displaystyle \frac{d[NH_4^+]}{dt} & = & -R_3 & \qquad (1d) \[5mm] \displaystyle \frac{d[NH_3]}{dt} & = & R_3 + \color{red}{R_N} & \qquad (1e) \[5mm] \displaystyle \frac{d[H_2O]}{dt} & = & -R_4 & \qquad (1f) \[5mm] \displaystyle \frac{d[OH^-]}{dt} & = & R_4 & \qquad (1g) \[5mm] \displaystyle \frac{d[H^+]}{dt} & = & R_1 + R_2 + R_3 + R_4 & \qquad (1h) \[5mm] \end{array} $$
where the individual rates ($R_1$, $R_2$, $R_3$ and $R_4$) are defined as the difference between the forward and backward reaction rates:
$$ \begin{array}{rclr} R_1 &=& k_{1f} \cdot [H_2CO_3] - k_{1b} \cdot [HCO_3^-] \cdot [H^+] & \qquad\hspace{2cm} (2a) \[4mm] R_2 &=& k_{2f} \cdot [HCO_3^{-}] - k_{2b} \cdot [CO_3^{2-}] \cdot [H^+] & \qquad\hspace{2cm} (2b) \[4mm] R_3 &=& k_{3f} \cdot [NH_4^{+}] - k_{3b} \cdot [NH_3] \cdot [H^+] & \qquad\hspace{2cm} (2c) \[4mm] R_4 &=& k_{4f} \cdot [H_2O] - k_{4b} \cdot [OH^-] \cdot [H^+] & \qquad\hspace{2cm} (2d) \end{array} $$
Note that the slow processes, $R_C$ and $R_N$, only affect $H_2CO_3$ and $NH_3$.
Task 2: Equilibrium conditions
In an equilibrium, the time derivative of each species in the system is zero. Based on equations 1a and 2a, we obtain $R_1=0$, which yields the relationship
$$ K_{1} = \frac{ [HCO_3^-] \cdot [H^+] }{ [H_2CO_3] }. \qquad\hspace{2cm} (3a) $$
Similarly, based on Eq. 1c and 2b, we obtain $R_2=0$, which yields the relationship
$$ K_{2} = \frac{ [CO_3^{2-}] \cdot [H^+] }{[HCO_3^-]}. \qquad\hspace{2cm} (3b) $$ Based on Eq. 1d and 2c, we obtain $R_3=0$, which yields the relationship
$$ K_{n} = \frac{ [NH_3] \cdot [H^+] }{[NH_4^+]}. \qquad\hspace{2cm} (3c) $$
Finally, based on Eq. 1f and 2d, we obtain $R_4=0$, which yields the relationship
$$ K_{w}' = \frac{ [OH^-] \cdot [H^+] }{[H_2O]}. \qquad\hspace{2cm} (3d) $$
Here, we defined the dissociation constants $K_1$, $K_2$, $K_n$, $K_w'$ according to $K_1 = k_{1f}/k_{1b}$, $K_2 = k_{2f}/k_{2b}$, $K_n = k_{3f}/k_{3b}$, and $K_w' = k_{4f}/k_{4b}$. If these relationships are valid, then the time derivatives in equations 1b, 1e, 1g and 1h are zero. This means that these four equations do not provide any extra information for the equilibrium concentrations.
Task 3: Lump-sum species
Here we combine equations 1a--1h and see what they imply for the suggested lump-sum species. By summing equations 1a, 1b and 1c, we obtain
$$ \frac{d[DIC]}{dt} = {\color{red}{R_C}} \qquad\hspace{2cm} (4a) $$
Similarly, by summing equations 1d and 1e, we obtain
$$ \frac{d[NH_x]}{dt} = {\color{red}{R_N}} \qquad\hspace{2cm} (4b) $$
Similarly, by summing equations 1f and 1g, we obtain
$$ \frac{d[H_xO]}{dt} = 0 \qquad\hspace{2cm} (4c) $$
Finally, by summing equations according to (1b) + $2\times$(1c) + (1g) $-$ (1d) $-$ (1h), we obtain $$ \frac{d[ALK]}{dt} = 0 \qquad\hspace{2cm} (4d) $$ These equations show that DIC is added at a rate of the slow process, $R_C$, $NH_x$ is added at a rate of the slow process, $R_N$, while $H_xO$ and $ALK$ are conserved (i.e., not changing over time). Specifically, none of the lump-sum species $DIC$, $NH_x$, $H_xO$ and $ALK$ is affected by the fast reversible reactions.
Task 4: From lump-sum species to the original species
This task has been solved in parts I and II of this reader series. Here we only summarize the results:
$$ \begin{array}{rclr} [HCO_3^-] &=& \displaystyle \frac{K_{1}[H^+]}{[H^+][H^+] + K_{1}[H^+] + K_{1} K_{2}}\cdot [DIC] & \qquad\hspace{2cm} (5a) \[5mm] [H_2CO_3] &=& \displaystyle \frac{[H^+][H^+]}{[H^+][H^+] + K_{1}[H^+] + K_{1} K_{2}}\cdot [DIC] & \qquad\hspace{2cm} (5b) \[5mm] [CO_3^{2-}] &=& \displaystyle \frac{K_1 K_2}{[H^+][H^+] + K_{1}[H^+] + K_{1} K_{2}}\cdot [DIC] & \qquad\hspace{2cm} (5c) \[5mm] [NH_3] &=& \displaystyle \frac{K_n}{K_n + [H^+]}\cdot [NH_x] & \qquad\hspace{2cm} (5d) \[5mm] [NH_4^+] &=& \displaystyle \frac{[H^+]}{K_n + [H^+]} \cdot [NH_x] & \qquad\hspace{2cm} (5e) \[5mm] [OH^-] &=& \displaystyle \frac{K_w'}{K_w' + [H^+]}\cdot [H_xO] & \qquad\hspace{2cm} (5f) \[5mm] [H_2O] &=& \displaystyle \frac{[H^+]}{K_w' + [H^+]}\cdot [H_xO] & \qquad\hspace{2cm} (5g) \end{array} $$
Once the concentrations $[HCO_3^-]$, $[CO_3^{2-}]$, $[NH_4^+]$ and $[OH^-]$ are calculated from equations 5a, 5c, 5e and 5f, respectively, the alkalinity is calculated from the definition: $$ [ALK] = [HCO_3^-] + 2\cdot [CO_3^{2-}] + [OH^-] - [NH_4^+] - [H^+] \qquad\hspace{2cm} (5h) $$ Because $[HCO_3^-]$ and $[CO_3^{2-}]$ depend on $[DIC]$ and $[H^+]$ (Eqs. 5a and 5c), $[NH_4^+]$ depends on $[NH_x]$ and $[H^+]$ (Eq. 5d), and $[OH^-]$ depends on $[H_xO]$ and $[H^+]$ (Eq. 5f), equation 5h provides an intimate relationship between $[ALK]$, $[DIC]$, $[NH_x]$, $[H_xO]$ and $[H^+]$ when the carbonate, ammonium and water species are in equilibrium. Using this relationship, one can calculate any of the five species provided that the other \emph{four} species are known. Because equation 5h is highly non-linear, these calculations are done numerically.
Task 5: Application to organic matter mineralization
R implementation
First, we expand the function that calculates alkalinity by including the concentration of $NH_4^+$ (Eq. $5e$).
solveALK <- function(K1=4.44e-07, K2=4.67e-11, Kn=5.68e-10, Kw.prime=1.79e-16, HxO=55.3, DIC=0, NHx=0, pH=7){ H <- 10^(-pH) HCO3 <- K1*H / (H^2 + K1*H + K1*K2) * DIC # eq 5a CO3 <- K1*K2 / (H^2 + K1*H + K1*K2) * DIC # eq 5c NH4 <- H / (H + Kn) * NHx # eq 5e OH <- Kw.prime / (Kw.prime + H) * HxO # eq 5f ALK <- HCO3 + 2*CO3 + OH - NH4 - H # eq 5h return(ALK) }
Next, we update the function that calculates $[H^+]$ from given concentrations of $ALK$, $DIC$, $NH_x$ and $H_xO$ (Eq. $5h$). This is done by finding a root of a function defined as the left-hand-side minus the right-hand-side of Eq. $5h$.
solveH <- function(K1=4.44e-07, K2=4.67e-11, Kn=5.68e-10, Kw.prime=1.79e-16, HxO=55.3, DIC=0, NHx=0, ALK=0){ # function whose root has to be sought rootFun <- function(H) { ALK.est <- solveALK(K1=K1, K2=K2, Kn=Kn, Kw.prime=Kw.prime, HxO=HxO, DIC=DIC, NHx=NHx, pH=-log10(H)) return(ALK.est - ALK) } # uniroot will find the root h <- uniroot(f = rootFun, lower = 0, upper = 1, tol = 1e-30) return(h$root) }
It is always good to test the functions. We see that the results are consistent with expectations.
pH.new <- -log10(solveH(DIC=2e-3, NHx=1e-5, ALK=2e-3)) ALK <- solveALK(DIC=2e-3, NHx=1e-5, pH = pH.new) pH.neutral <- -log10(solveH(DIC=0, NHx=0, ALK=0)) c(pH.new=pH.new, ALK=ALK, pH.neutral=pH.neutral)
Next, we implement the model function based on the local equilibrium assumption. Note that
-
The equilibrium constants (dissociation constants) are in $mol~L^{-1}$, so all species must be in this unit, too!
-
$OM$ is added to the model as a state variable removed by the slow process.
-
Because we will additionally want to use the function for modelling ammonia degassing (Task 6), the effect of this process is added, too. Its rate constant ($\lambda$) will be set to 0 by default.
OMMpHmodel <- function(t, state, parms) { with (as.list(c(state, parms)), { # variables needed as output and for regulation of process rates H <- solveH(K1=K1, K2=K2, Kn=Kn, Kw.prime=Kw.prime, HxO=HxO, DIC=DIC, NHx=NHx, ALK=ALK) HCO3 <- K1*H / (H^2 + K1*H + K1*K2) * DIC # eq 5a CO3 <- K1*K2 / (H^2 + K1*H + K1*K2) * DIC # eq 5c NH4 <- H / (H + Kn) * NHx # eq 5e OH <- Kw.prime / (Kw.prime + H) * HxO # eq 5f NH3 <- NHx - NH4 # from definition of NHx H2CO3 <- DIC - HCO3 - CO3 # from definition of DIC H2O <- HxO - OH # from definition of HxO # process rates RC <- rC*OM # OM degradation Rdegas <- lambda * NH3 # ammonia degassing (specifically removes NH3) # mass balance equations (assuming local equilibrium) dOM <- -RC dDIC <- RC # Eq. 4a dNHx <- RC*NC - Rdegas # Eq. 4b - ammonia degassing dHxO <- 0 # Eq. 4c dALK <- 0 # Eq. 4d return(list(c(dOM, dDIC, dNHx, dHxO, dALK), pH = -log10(H), H = H, OH = OH, H2O = H2O, H2CO3 = H2CO3, HCO3 = HCO3, CO3 = CO3, NH4 = NH4, NH3 = NH3)) }) }
We run the model for 10 days, using the initial conditions given in the task.
require(deSolve) pars <- c(rC = 0.1, # [1/day] OM degradation lambda = 0, # [1/day] ammonia degasing NC = 16/106, # molN/molC K1 = 4.44e-07, # mol/L K2 = 4.67e-11, # mol/L Kn = 5.68e-10, # mol/L Kw.prime = 1.79e-16) # mol/L pH.ini <- 7.9 # H.ini based on pH OM.ini <- 100*1e-6 # mol/L DIC.ini <- 2000*1e-6 # mol/L NHx.ini <- 10*1e-6 # mol/L HxO.ini <- 55.3 # mol/L ALK.ini <- solveALK(HxO = HxO.ini, DIC = DIC.ini, NHx = NHx.ini, pH = pH.ini) state.ini <- c(OM=OM.ini, DIC=DIC.ini, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini) times <- seq(from=0, to=10, length.out=100) out <- ode(y=state.ini, times=times, func=OMMpHmodel, parms=pars)
plot(out, mfrow = c(4,4), las = 1)
We see that mineralization of about $60~\mu mol~L^{-1}$ of organic matter over the period of 10 days leads to a decrease in $pH$ from 7.9 to about 7.63. As expected, $ALK$ and $H_xO$ are conserved and changes in the concentration of $H_2O$ are negligible due to vastly larger concentrations of water compared to all the other species in the system.
Task 6: Does the formulation of slow reversible processes matter?
If the stoichiometry of OM mineralization is represented by the reaction $$ \quad CH_2O(NH_3)_{NC} + O_2 \rightarrow HCO_3^- + H^- + {NC} \cdot NH_3, $$ the mass balance equations 1c--g do not change while equations 1a, 1b, and 1h change in the following way:
$$ \begin{array}{rclr} \displaystyle \frac{d[H_2CO_3]}{dt} & = & R_1 & \qquad (1a')\[5mm] \displaystyle \frac{d[HCO_3^-]}{dt} & = & -R_1 + R_2 + \color{red}{R_C}& \qquad (1b') \[5mm] \displaystyle \frac{d[H^+]}{dt} & = & -R_1 - R_2 - R_3 - R_4 + \color{red}{R_C} & \qquad (1h') \[5mm] \end{array} $$ This change, however, has no effect on the differential equations for the lump-sum variables $DIC$, $NH_x$, $H_xO$ and $ALK$, as can easily be verified by following the same derivation steps as those leading to equations 4a--d. This can also be shown more easily: although $R_C$ adds $HCO_3^-$ instead of $H_2CO_3$, these two species are both included in the $DIC$ pool, thus summing them up yields the same result for the time-derivative of $DIC$. Similarly, although the process increases $ALK$ at a rate $R_C$ due to the addition of $HCO_3^-$, it also decreases $ALK$ at the same rate $R_C$ due to the addition of $H^+$, as follows from the negative sign in front of $[H^+]$ in the definition of $ALK$.
This result has important implications: the specific "chemical" formulation of the slow reversible process has no effect on the dynamics of the lump-sum variables, as long as the chemical reactions are correctly balanced in terms of the elements and charges and the local equilibrium conditions can be assumed.
Task 7: Application to ammonia degassing
R implementation
First, we reproduce here the ammonia degassing model developed in Part I:
solveB <- function(Keq, Atot, Btot){ rootFun <- function(B) B + B/(Keq+B)*Atot - Btot return(uniroot(f = rootFun, lower = 0, upper = Btot, tol = 1e-20)$root) } AmmoniaDegassing <- function(t, state, parms) { with (as.list(c(state, parms)), { H <- solveB(Keq = Kn, Atot = NHx, Btot = Hx) NH3 <- Kn/(Kn+H) * NHx NH4 <- NHx - NH3 # mass balance equations for the lump-sum variables: dNHx <- -lambda * NH3 dHx <- 0 return(list(c(dNHx, dHx), H = H, NH3 = NH3, NH4 = NH4, # NH4 corresponds to [NH4+] pH = -log10(H), # return also pH Q = NH3*H/NH4 )) # return also Q }) }
We define the initial conditions and run the model. In this step, we reuse the parameters defined for the OMMpHmodel model with the value of the degassing rate constant, $\lambda$, changed to $1~d^{-1}$.
# update model parameters pars1 <- pars pars1["lambda"] <- 1 # initial conditions NHx.ini <- 1e-3 # mol/L pH.ini <- 8 H.ini <- 10^(-pH.ini) # mol/L Hx.ini <- H.ini + H.ini/(pars1[["Kn"]] + H.ini) * NHx.ini state.ini1 <- c(NHx = NHx.ini, Hx = Hx.ini) # run the model out1 <- ode(y=state.ini1, times=times, func=AmmoniaDegassing, parms=pars1)
Next, we run the OMMpHmodel model with the same initial conditions. By setting the OM degradation rate constant to 0, and further setting the concentrations of $DIC$ and $H_xO$ to 0, we expect that the model becomes equivalent to the AmmoniaDegassing model.
# update model parameters pars2 <- pars1 pars2["rC"] <- 0 # initial conditions ALK.ini <- solveALK(HxO=0, DIC=0, NHx=NHx.ini, pH=pH.ini) state.ini2 <- c(OM=OM.ini, DIC=0, NHx=NHx.ini, HxO=0, ALK=ALK.ini) # run the model out2 <- ode(y=state.ini2, times=times, func=OMMpHmodel, parms=pars2)
Thus, we expect to obtain the same results, at least for the state variables that are common for both models. This is indeed the case when we plot the results.
plot(out2, obs=out1, which=c("NHx", "NH4", "NH3", "pH"), mfrow = c(1,4), lwd=2, lty=1, col=1, obspar = list(lty=2,lwd=2,col=2,type="l"))
Now that we have verified that the model OMMpHmodel can be used for solving the ammonia degassing problem, too, we can illustrate the $pH$ buffering role of $DIC$ in natural waters. We will use the same initial conditions and the ammonia degassing rate constant as above and perform the calculations for $DIC$ concentrations of 0, 1, and $2~mmol~L^{-1}$.
# without DIC (but with ammonia and water!) ALK.ini3 <- solveALK(DIC=0, HxO=HxO.ini, NHx=NHx.ini, pH=pH.ini) state.ini3 <- c(OM=OM.ini, DIC=0, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini3) out3 <- ode(y=state.ini3, times=times, func=OMMpHmodel, parms=pars2) # model with DIC=1 mmol/L ALK.ini4 <- solveALK(DIC=1e-3, HxO=HxO.ini, NHx=NHx.ini, pH=pH.ini) state.ini4 <- c(OM=OM.ini, DIC=1e-3, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini4) out4 <- ode(y=state.ini4, times=times, func=OMMpHmodel, parms=pars2) # model with DIC=2 mmol/L ALK.ini5 <- solveALK(DIC=2e-3, HxO=HxO.ini, NHx=NHx.ini, pH=pH.ini) state.ini5 <- c(OM=OM.ini, DIC=2e-3, NHx=NHx.ini, HxO=HxO.ini, ALK=ALK.ini5) out5 <- ode(y=state.ini5, times=times, func=OMMpHmodel, parms=pars2)
# plot results plot(out3, out4, out5, which=c("pH", "H", "OH", "HxO", "DIC", "ALK"), mfrow = c(2,3), las = 1, col=1:3, lty=1:3, lwd=2) legend("topright", legend = c("DIC=0", "DIC=1 mmol/L", "DIC=2 mmol/L"), bty="n", col = 1:3, lwd = 2, lty = 1:3)
# plot results plot(out3, out4, out5, which=c("NHx", "NH4", "NH3", "H2CO3", "HCO3", "CO3"), mfrow = c(2,3), las = 1, col=1:3, lty=1:3, lwd=2) legend("topright", legend = c("DIC=0", "DIC=1 mmol/L", "DIC=2 mmol/L"), bty="n", col = 1:3, lwd = 2, lty = 1:3)
The results show that when $DIC$ is present, the protons released due to ammonia degassing are scavenged by both $CO_3^{2-}$ and $HCO_3^-$ in the solution, leading to an increase in $H_2CO_3$. Consequently, the increase in the concentration of protons, and thus the decrease in $pH$, is smaller in the presence of $DIC$, illustrating the $pH$ buffering effect by the carbonate system. The magnitude of the effect does not, however, scale linearly with the $DIC$ concentration.
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.