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 for each reaction. Considering the stoichiometry of Reactions 1--3, 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[OH^-]}{dt} & = & R_3 & \qquad (1d)\[5mm] \displaystyle \frac{d[H^+]}{dt} & = & R_1 + R_2 + R_3 & \qquad (1e)\[5mm] \displaystyle \frac{d[H_2O]}{dt} & = & -R_3 & \qquad (1f)\[5mm] \end{array} $$

where the individual net rates ($R_1$, $R_2$ and $R_3$) are defined as the difference between the forward and backward reaction rates:

$$ \begin{array}{rcll} R_1 & = & k_{1f} \cdot [H_2CO_3] - k_{1b} \cdot [HCO_3^-] \cdot [H^+] & \hspace{2cm} (2a) \[5mm] R_2 & = & k_{2f} \cdot [HCO_3^{-}] - k_{2b} \cdot [CO_3^{2-}] \cdot [H^+] & \hspace{2cm} (2b) \[5mm] R_3 & = & k_{3f} \cdot [H_2O] - k_{3b} \cdot [OH^{-}] \cdot [H^+] & \hspace{2cm} (2c) \end{array} $$

Note that the rate of the slow process, $R_C$, only affects $H_2CO_3$.

Task 2: Equilibrium conditions

In the absence of the slow process (i.e., $R_C=0$), the carbonate species and $OH^-$ and $H^+$ will reach an equilibrium as a result of the reversible reactions 1--3. In this 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 equation 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) $$ Finally, the equilibrium condition implies $R_3=0$, which combined with Eq. 2c yields $$ K_w' = \frac{ [OH^-] \cdot [H^+] }{[H_2O]}. \qquad\hspace{2cm} (3c) $$ Here, we defined the dissociation constants according to $K_1 = k_{1f}/k_{1b}$, $K_2 = k_{2f}/k_{2b}$, and $K_w' = k_{3f}/k_{3b}$. If relationships 3a--3c are valid, then the time derivatives in equations 1b, 1d and 1e are zero as well. This means that these latter three equations do not provide any extra information for the equilibrium concentrations.

Note: In literature, the ionic product of water, $K_w = [OH^-]\cdot [H^+] \approx 10^{-14}~mol^2~L^{-2}$ is more commonly reported in the context of water self-ionization. This value should, however, not be confused with the equilibrium constant $K_w'$. The correct value of $K_w'$ is calculated by dividing $K_w$ with the concentration of water, $[H_2O]\approx 55.3~mol~L^{-1}$, as shown in Eq. 3c (see table in Task 2).

Task 3: Lump-sum species

By summing equations 1a, 1b and 1c, and additionally considering the proposed definition of the lump-sum variable $[DIC]$, we obtain $$ \frac{d[DIC]}{dt} = R_C. \qquad\hspace{2cm} (4a) $$ Similarly, by summing equation 1b with equation 1d and equation 1c multiplied with a factor of 2, then subtracting equation 1e, and additionally considering the proposed definition of the lump-sum variable $[ALK]$, we obtain $$ \frac{d[ALK]}{dt} = 0. \qquad\hspace{2cm} (4b) $$ Finally, by summing equation 1d with equation 1f and additionally considering the proposed definition of the lump-sum variable $[H_xO]$, we obtain $$ \frac{d[H_xO]}{dt} = 0. \qquad\hspace{2cm} (4c) $$

Key take-home message

These equations show that under the assumption of a local chemical equilibrium, $DIC$ is added at a rate of the slow process ($R_C$), while $ALK$ and $H_xO$ are conserved (not changing over time). Notably, $DIC$, $ALK$ and $H_xO$ are not affected by the fast reversible reactions!

Task 4: From lump-sum species to original species

By solving for $[H_2CO_3]$ and $[CO_3^{2-}]$ in equations 3a and 3b, and by substituting the result to the definition of $[DIC]$, we obtain

$$ [DIC] = [HCO_3^-] + [H^+]\cdot \frac{[HCO_3^-]}{K_{1}} + K_{2}\cdot \frac{ [HCO_3^-]}{[H^+]} $$

When we express each term on the right-hand side using the common denominator ($K_1\cdot [H^+]$), we obtain

$$ [DIC] = \frac{K_1\cdot[H^+]\cdot[HCO_3^-]}{K_1\cdot[H^+]} + \frac{ [H^+]^2\cdot[HCO_3^-]}{K_1 \cdot[H^+]} + \frac{ K_1 \cdot K_{2} \cdot[HCO_3^-]}{K_1 \cdot [H^+]} $$

This equation can easily be solved for $[HCO_3^-]$:

$$ \boxed{~[HCO_3^-]= \frac{K_{1}\cdot [H^+]}{[H^+]\cdot [H^+] + K_{1}\cdot [H^+] + K_{1}\cdot K_{2}}\cdot [DIC]~} \qquad\hspace{2cm} (5a) $$

Combining this result with equation 3a, we obtain

$$ \boxed{~[H_2CO_3] = \frac{[H^+]\cdot [H^+]}{[H^+]\cdot [H^+] + K_{1}\cdot [H^+] + K_{1}\cdot K_{2}}\cdot [DIC]~} \qquad\hspace{2cm} (5b) $$

Similarly, combining the result with equation 3b, we obtain

$$ \boxed{~[CO_3^{2-}] = \frac{K_1\cdot K_2}{[H^+]\cdot [H^+] + K_{1}\cdot [H^+] + K_{1}\cdot K_{2}}\cdot [DIC]~} \qquad\hspace{2cm} (5c) $$

Regarding the concentrations of $[OH^-]$ and $[H_2O]$, they can be calculated from the relationship 3c. First, by substituting $[H_2O]$ based on the definition of the lump-sum $H_xO$, we obtain $$ K_w' = \frac{[OH^-]\cdot [H^+]}{[H_xO] - [OH^-]}, $$ which upon rearrangement yields the expression for $[OH^-]$ as a function of the lump-sum $[H_x0]$ and $[H^+]$: $$ \boxed{~[OH^-] = \frac{K_w'}{K_w' + [H^+]}\cdot [H_xO]~} \qquad\hspace{2cm} (5d) $$ Finally, by combining this last result with the definition of $H_xO$, we obtain $$ [H_2O] = [H_xO] - [OH^-] = \left( 1-\frac{K_w'}{K_w' + [H^+]} \right) \cdot [H_xO], $$ which yields the expression for $[H_2O]$ as a function of the lump-sum $[H_x0]$ and $[H^+]$: $$ \boxed{~[H_2O] = \frac{[H^+]}{K_w' + [H^+]}\cdot [H_xO].~} \qquad\hspace{2cm} (5e) $$

Once the concentrations $[HCO_3^-]$, $[CO_3^{2-}]$ and $[OH^-]$ are calculated from equations 5a--5d based on the values of $[DIC]$, $[H_xO]$ and $[H^+]$, the alkalinity is calculated from the definition:

$$ [ALK] = [HCO_3^-] + 2\cdot [CO_3^{2-}] + [OH^-] - [H^+] \qquad\hspace{2cm} (5f) $$

Key take-home message

Because $[HCO_3^-]$ and $[CO_3^{2-}]$ depend on $[DIC]$ and $[H^+]$ (Eqs. $5a$ and $5c$), and $[OH^-]$ depends on $[H_xO]$ and $[H^+]$ (Eq. $5d$), equation 5f provides an intimate relationship between $[ALK]$, $[DIC]$, $[H_xO]$ and $[H^+]$ when the carbonate and water species are in equilibrium. Using this relationship, one can calculate

1. $[ALK]$ if the values of $[DIC]$, $[H_xO]$ and $[H^+]$ are known, or

2. $[DIC]$ if the values of $[H^+]$, $[H_xO]$ and $[ALK]$ are known, or

3. $[H^+]$ if the values of $[DIC]$, $[ALK]$ and $[H_xO]$ are known.

Note that equation $5f$ is highly non-linear; thus these calculations need to be done numerically. This contrasts to part I of this reader series, where it was possible to also find an analytical solution. Although this is no longer possible when the number of fast acid-base reactions (e.g., Reactions 1--3) increases, we can proceed thanks to numerical solvers (see next section).

Another important point to note here is the statement 3: $[H^+]$ can be calculated if the values of $[DIC]$, $[ALK]$ and $[H_xO]$ are known. This statement is the basis for the general approach in $pH$ modeling: one needs to first identify the appropriate lump-sum species (e.g., $DIC$, $ALK$ and $H_xO$) and then formulate the effect of each slow process on these lump-sum species. For example, in this particular exercise, the effects are: (i) the addition or removal of $DIC$ (Eq. 4a), (ii) conservation of $ALK$ (Eq. 4b), and (iii) conservation of $H_xO$ (Eq. 4c). Then, by solving the differential equations for the lump-sum species numerically, one can quantify the corresponding evolution of $pH$ based on the statement 3, i.e., by numerically solving for $[H^+]$ from Eq. 5f. This approach will be illustrated in the following section.

Task 5: Application to atmospheric CO$_2$ exchange

Implementation in R

First, we define a function that calculates alkalinity from the known concentrations of $DIC$ and $H_xO$ and $pH$ (equation 5f).

solveALK <- function(K1 = 4.44e-07, K2 = 4.67e-11, Kw.prime = 1.79e-16, 
                     HxO = 55.3, DIC, pH){

  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
  OH     <- Kw.prime / (Kw.prime + H) * HxO       # eq 5d
  ALK <- HCO3 + 2*CO3 + OH - H                    # eq 5f
  return(ALK)
}

Next, we define a function that calculates $[H^+]$ from the known concentrations of $DIC$, $ALK$ and $H_xO$ (equation 5f). This is done by finding the root of a function defined as the left-hand-side minus the right-hand-side of equation 5f, using a numerical solver implemented in uniroot.

solveH <- function(K1 = 4.44e-07, K2 = 4.67e-11, Kw.prime = 1.79e-16,
                   HxO = 55.3, DIC, ALK){

  # function whose root has to be sought
  rootFun <- function(H) {

    HCO3    <- K1*H  / (H^2 + K1*H + K1*K2) * DIC    # eq 5a
    CO3     <- K1*K2 / (H^2 + K1*H + K1*K2) * DIC    # eq 5c
    OH      <- Kw.prime / (Kw.prime + H) * HxO       # eq 5d
    ALK.est <- HCO3 + 2*CO3 + OH - H                 # eq 5f
    return(ALK.est - ALK)
  }

  # uniroot will find the root    
  h <- uniroot(f = rootFun, lower = 0, upper = 1e-3, tol = 1e-30)
  return(h$root)
}

It is always good to test the functions. For pure water with zero $[DIC]$ and $[ALK]$, we expect to get $pH$ of about 7 (i.e., neutral $pH$). This follows from Eqs. 3c and 5f, because at these conditions, we have $[OH^-]=[H^+]$ and $[OH^-]\cdot [H^+] = K_w'\cdot [H_2O]$, which yields $[H^+]=\sqrt{K_w'\cdot [H_xO]}\approx 10^{-7}$ and so $pH=-log_{10}([H^+])=7$. The result is as expected:

(pH.neutral  <- -log10(solveH(DIC = 0, ALK = 0)))

In contrast, for $[DIC] = 2~mmol~L^{-1}$ and $[ALK] = 2~mmol~L^{-1}$, we expect to get $pH$ of around $8$. Furthermore, when we substitute that $pH$ value to calculate back $[ALK]$, we expect to get the original value. Indeed, the results are as expected:

(pH.new  <- -log10(solveH(DIC = 2e-3, ALK = 2e-3)))
(ALK     <-  solveALK(DIC = 2e-3, pH = pH.new))

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!

CO2dissol <- function(t, state, parms) {
  with (as.list(c(state, parms)), {

    # mass balance equations
    dDIC <- RC
    dALK <- 0
    dHxO <- 0

    # calculate the relevant species for the output
    H    <- solveH(K1 = K1, K2 = K2, Kw.prime = Kw.prime, 
                   HxO = HxO, DIC = DIC, 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
    OH   <- Kw.prime / (Kw.prime + H) * HxO      # eq 5d

    return(list(c(dDIC, dALK, dHxO), 
                pH = -log10(H), 
                H = H, OH = OH, H2O = HxO - OH,
                H2CO3 = DIC - HCO3 - CO3,
                HCO3 = HCO3,
                CO3  = CO3))
  })
}

We run the model for 10 days, using the initial conditions given in the task.

require(deSolve)

pars    <- c(RC = 0.1e-3,          # [mol/L/day]
             K1 = 4.44e-07,        # [mol/L]
             K2 = 4.67e-11,        # [mol/L]
             Kw.prime = 1.79e-16)  # [mol/L]

pH.ini  <- 8
DIC.ini <- 2e-3
HxO.ini <- 55.3
ALK.ini <- solveALK(HxO = HxO.ini, DIC = DIC.ini, pH = pH.ini)

yini    <- c(DIC = DIC.ini, ALK = ALK.ini, HxO = HxO.ini)
times   <- seq(from=0, to=10, length.out=100)

out     <- ode(y=yini, times=times, func=CO2dissol, parms=pars)

plot(out, mfrow = c(2,5), las = 1)

We see that after 10 days, the addition of $CO_2$ at a rate or r pars["RC"] $mol~L^{-1}~d^{-1}$ leads to water acidification, i.e., a decrease in $pH$ from 8 to about 6.6.

dynamic-R/RTM documentation built on Feb. 28, 2025, 1:23 p.m.