The governing mass conservation equations for the geochemical transport mode for a multiphase system written in terms of a set of independent aqueous primary or basis species with the form

(1)\[\frac{{{\partial}}}{{{\partial}}t}\big(\varphi \sum_{{\alpha}}s_{{\alpha}}\Psi_j^{{\alpha}}\big) + \nabla\cdot\sum_{{\alpha}}{\boldsymbol{\Omega}}_j^{{\alpha}}= Q_j - \sum_m\nu_{jm} I_m -\frac{{{\partial}}S_j}{{{\partial}}t},\]


(2)\[\frac{{{\partial}}\varphi_m}{{{\partial}}t} = \overline{V}_m I_m,\]

for minerals with molar volume \(\overline{V}_m\), mineral reaction rate \(I_m\) and mineral volume fraction \(\varphi_m\) referenced to an REV. Sums over \({{\alpha}}\) in Eqn. (1) are over all fluid phases in the system. The quantity \(\Psi_j^{{\alpha}}\) denotes the total concentration of the \(j\)th primary species \({{\mathcal A}}_j^{\rm pri}\) in the \({{\alpha}}\)th fluid phase defined by

(3)\[\Psi_j^{{\alpha}}= \delta_{l{{\alpha}}}^{} C_j^l + \sum_{i=1}^{N_{\rm sec}}\nu_{ji}^{{{\alpha}}} C_i^{{\alpha}},\]

In this equation the subscript \(l\) represents the aqueous electrolyte phase from which the primary species are chosen. The secondary species concentrations \(C_i^{{\alpha}}\) are obtained from mass action equations corresponding to equilibrium conditions of the reactions

(4)\[\sum_j\nu_{ji}^{{\alpha}}{{\mathcal A}}_j^l {~\rightleftharpoons~}{{\mathcal A}}_i^{{\alpha}},\]

yielding the mass action equations

(5)\[C_i^{{\alpha}}= \frac{K_i^{{\alpha}}}{\gamma_i^{{\alpha}}} \prod_j \Big(\gamma_j^l C_j^l\Big)^{\nu_{ji}^{{\alpha}}},\]

with equilibrium constant \(K_i^{{\alpha}}\), and activity coefficients \(\gamma_k^{{\alpha}}\). For the molality of the \(k\)th aqueous species, the Debye-Hückel activity coefficient algorithm is given by

(6)\[\log\,\gamma_k = -\frac{z_k^2 A \sqrt{I}}{1+B \stackrel{\circ}{a}_k \sqrt{I}}+\dot b I,\]

and the Davies algorithm by the expression

(7)\[\log\,\gamma_k = -\frac{z_k^2}{2}\left[\frac{\sqrt{I}}{1+ \sqrt{I}}-0.3 I\right].\]

with valence \(z_k\), Debye-Hückel parameters \(A\), \(B\), and ionic radius \(\stackrel{\circ}{a}_k\), and ionic strength \(I\) defined as

(8)\[I = \frac{1}{2}\sum_{j=1}^{N_c} m_j^2 + \frac{1}{2}\sum_{i=1}^{N_{\rm sec}} m_i^2,\]

for molality \(m_j\) and \(m_i\) of primary and secondary species, respectively (note: \(C_i^l = \rho_l y_w^l m_i \simeq \rho_l m_i\), \(\rho_l\) = fluid density, \(y_w^l\) = mass fraction of \(\mathrm{H_2O}\)). For high-ionic strength solutions (approximately above 0.1 M) the Pitzer model should be used. Currently, however, only the Debye-Hückel algorithm is implemented in PFLOTRAN.

The total flux \({\boldsymbol{\Omega}}_j^{{\alpha}}\) for species-independent diffusion is given by

(9)\[{\boldsymbol{\Omega}}_j^{\alpha}= \big({\boldsymbol{q}}_{\alpha}- \varphi s_{\alpha}{\boldsymbol{D}}_{\alpha} \cdot {\boldsymbol{\nabla}}\big)\Psi_j^{\alpha}.\]

The diffusion/dispersion tensor \({\boldsymbol{D}}_{\alpha}\) may be different for different phases, e.g. an aqueous electrolyte solution or gas phase, but is assumed to be species independent. Dispersivity currently must be described through a diagonal dispersion tensor.

The Darcy velocity \({\boldsymbol{q}}_{{\alpha}}\) for phase \({{\alpha}}\) is given by

(10)\[{\boldsymbol{q}}_a = -\frac{kk_{{\alpha}}}{\mu_{{\alpha}}} {\boldsymbol{\nabla}}\big(p_{{\alpha}}-\rho_{{\alpha}}g z\big),\]

with bulk permeability of the porous medium \(k\) and relative permeability \(k_{{\alpha}}\), fluid viscosity \(\mu_{{\alpha}}\), pressure \(p_{{\alpha}}\), density \(\rho_{{\alpha}}\), and acceleration of gravity \(g\). The diffusivity/dispersivity tensor \({\boldsymbol{D}}_{{\alpha}}\) is the sum of contributions from molecular diffusion and dispersion which for an isotropic medium has the form

(11)\[{\boldsymbol{D}}_{{\alpha}}= \tau D_m {\boldsymbol{I}}+ a_T v{\boldsymbol{I}}+ \big(a_L-a_T\big)\frac{{\boldsymbol{v}}{\boldsymbol{v}}}{v},\]

with longitudinal and transverse dispersivity coefficients \(a_L\), \(a_T\), respectively, \(\tau\) refers to tortuosity, and \(D_m\) to the molecular diffusion coefficient. Currently, only a diagonal dispersion tensor with principal axes aligned with the grid for longitudinal and transverse dispersion is implemented in PFLOTRAN.

The porosity may be calculated from the mineral volume fractions according to the relation

(12)\[\varphi = 1 - \sum_m \varphi_m.\]

The temperature dependence of the diffusion coefficient is defined through the relation

(13)\[D_m(T) = D_m^\circ\exp\left[\frac{A_D}{R}\left(\frac{1}{T_0}-\frac{1}{T}\right)\right],\]

with diffusion activation energy \(A_D\) in kJ/mol. The quantity \(D_m^\circ\) denotes the diffusion coefficient at the reference temperature \(T_0\) taken as 25\(^\circ\)C and the quantity \(R\) denotes the gas constant (\(8.317\times 10^{-3}\) kJ/mol/K). The temperature \(T\) is in Kelvin.

The quantity \(Q_j\) denotes a source/sink term

(14)\[Q_j = \sum_n\frac{q_M}{\rho}\Psi_j \delta({\boldsymbol{r}}-{\boldsymbol{r}}_{n}),\]

where \(q_M\) denotes a mass rate in units of kg/s, \(\rho\) denotes the fluid density in kg/m\(^3\), and \({\boldsymbol{r}}_{n}\) refers to the location of the \(n\)th source/sink. The quantity \(S_j\) represents the sorbed concentration of the \(j\)th primary species considered in more detail in the next section.

Molality \(m_i\) and molarity \(C_i\) are related by the density of water \(\rho_w\) according to

(15)\[C_i = \rho_w m_i.\]

The activity of water is calculated from the approximate relation

(16)\[a_{\rm H_2O}^{} = 1 - 0.017 \sum_i m_i.\]

Mineral Precipitation and Dissolution

The reaction rate \(I_m\) is based on transition state theory taken as positive for precipitation and negative for dissolution, with the form

(17)\[I_m = -a_m\left(\sum_l k_{ml}(T) {{{\mathcal P}}}_{ml}\right) \Big|1-\big(K_m Q_m\big)^{1/\sigma_m}\Big|^{\beta_m} {\rm sign}(1-K_mQ_m),\]

where the sum over \(l\) represents contributions from parallel reaction mechanisms such as pH dependence etc., and where \(K_m\) denotes the equilibrium constant, \(\sigma_m\) refers to Temkin’s constant which is defined as the average stoichiometric coefficient of the overall reaction (Lichtner, 1996b; see also Section [thermo:database]), \(\beta_m\) denotes the affinity power, \(a_m\) refers to the specific mineral surface area, and the ion activity product \(Q_m\) is defined as

(18)\[Q_m = \prod_j \big(\gamma_j m_j\big)^{\nu_{jm}},\]

with molality \(m_j\) of the \(j\)th primary species. The rate constant \(k_{ml}\) is a function of temperature given by the Arrhenius relation

(19)\[k_{ml} (T) = k_{ml}^0 \exp\left[\frac{E_{ml}}{R}\Big(\frac{1}{T_0}-\frac{1}{T}\Big)\right],\]

where \(k_{ml}^0\) refers to the rate constant at the reference temperature \(T_0\) taken as 25\(^\circ\)C, with \(T\) in units of Kelvin, \(E_{ml}\) denotes the activation energy (kJ/mol), and the quantity \({{{\mathcal P}}}_{ml}\) denotes the prefactor for the \(l\)th parallel reaction with the form

(20)\[{{{\mathcal P}}}_{ml} = \prod_i\dfrac{\big(\gamma_i m_i\big)^{{{\alpha}}_{il}^m}}{1+K_{ml}\big(\gamma_i m_i\big)^{{{\beta}}_{il}^m} },\]

where the product index \(i\) generally runs over both primary and secondary species, the quantities \(\alpha_{il}^m\) and \(\beta_{il}^m\) refer to prefactor coefficients, and \(K_{ml}\) is an attenuation factor. The quantity \(R\) denotes the gas constant (\(8.317 \times 10^{-3}\) kJ/mol/K).

Rate Limiter

A rate-limited form of the mineral kinetic rate law can be devised according to the expression

(21)\[\widehat I_m = -a_m^{} \sum_l \Bigg(\dfrac{{\mathcal P}_{ml}^{} k_{ml}^{} }{1+\dfrac{k_{ml}^{}}{k_{ml}^{\rm lim}} \left(K_m Q_m\right)^{1/\sigma_m}} \Bigg) \Big[ 1-\left(K_m Q_m\right)^{1/\sigma_m} \Big],\]

with rate-limiter \(r_{\rm lim}\). In the limit \(K_mQ_m\rightarrow\infty\), the rate becomes

(22)\[\lim_{K_m Q_m\rightarrow\infty}\widehat I_m = a_m^{}\sum_l k_{ml}^{\rm lim} {\mathcal P}_{ml}^{}.\]

Defining the affinity factor

(23)\[\Omega_m = 1-\left(K_m Q_m\right)^{1/\sigma_m},\]


(24)\[K_mQ_m = \Big(1-\Omega_m\Big)^{\sigma_m},\]

the rate may be expressed alternatively as

(25)\[\widehat I_m = -a_m^{} \sum_l \Bigg( \dfrac{{\mathcal P}_{ml}^{} k_{ml}^{} }{1+\dfrac{k_{ml}^{}}{k_{ml}^{\rm lim}} \big(1-\Omega_m\big)} \Bigg) \Omega_m.\]

Changes in Material Properties

Porosity, permeability, tortuosity and mineral surface area may be updated optionally due to mineral precipitation and dissolution reactions according to the relations

(26)\[\varphi = 1-\sum_m\varphi_m,\]
(27)\[k = k_0 f(\varphi,\,\varphi_0,\,\varphi_c,\,a),\]


(28)\[f = \left(\frac{\varphi-\varphi_c}{\varphi_0-\varphi_c}\right)^a,\]
(29)\[= f_{\rm min} \ \ \ \text{if} \ \ \ \varphi \leq \varphi_c,\]
(30)\[\tau = \tau_0 \left(\frac{\varphi}{\varphi_0}\right)^b,\]


(31)\[a_m = a_m^0 \left(\frac{\varphi_m}{\varphi_m^0}\right)^n \left(\frac{1-\varphi}{1-\varphi_0}\right)^{n'},\]

where the super/subscript 0 denotes initial values, with a typical value for \(n\) of \(2/3\) reflecting the surface to volume ratio. Note that this relation only applies to primary minerals \((\varphi_m^0\ne 0)\). The quantity \(\varphi_c\) refers to a critical porosity below which the permeability is assumed to be constant with scale factor \(f_{\rm min}\).

In PFLOTRAN the solid is represented as an aggregate of minerals described quantitatively by specifying its porosity \(\varphi\) and the volume fraction \(\varphi_m\) of each primary mineral. It is not necessary that Eqn. (26) relating porosity and mineral volume fractions holds. Typically, however, the solid composition is specified by giving the mass fraction \(y_m\) of each of the primary minerals making up the solid phase. The volume fraction is related to mole \(x_m\) and mass \(y_m\) fractions by the expressions

(32)\[\begin{split}\varphi_m &= (1-\varphi) \frac{x_m \overline V_m}{\sum_{m'} x_{m'} \overline V_{m'}},\\ &= (1-\varphi) \frac{y_m^{} \rho_m^{-1}}{\sum_{m'} y_{m'}^{} \rho_{m'}^{-1}},\end{split}\]

with inverse relation

(33)\[x_m = \frac{\varphi_m}{\overline V_m \eta_s(1-\varphi)},\]

and similarly for the mass fraction, where

(34)\[\rho_m^{} = W_m^{} \overline V_m^{-1},\]

and the solid molar density \(\eta_s\) is given by

(35)\[\eta_s = \frac{1}{\sum_m x_m \overline V_m}.\]

In these relations \(W_m\) refers to the formula weight and \(\overline{V_m}\) the molar volume of the \(m\)th mineral. The solid molar density is related to the mass density \(\rho_s\) by

(36)\[\rho_s = W_s \eta_s,\]

with the mean molecular weight \(W_s\) of the solid phase equal to

(37)\[W_s = \sum_m x_m W_m = \frac{1}{\sum_m W_m^{-1} y_m^{}}.\]

Mass and mole fractions are related by the expression

(38)\[W_m x_m = W_s y_m.\]

Affinity Threshold

An affinity threshold \(f\) for precipitation may be introduced which only allows precipitation to occur if \(K_m Q_m > f > 1\).

Surface Armoring

Surface armoring occurs when one mineral precipitates on top of another mineral, blocking that mineral from reacting. Thus suppose mineral \({{\mathcal M}}_m\) is being replaced by the secondary mineral \({{\mathcal M}}_{m'}\). Blocking may be described phenomenologically by the surface area relation

(39)\[a_m(t) = a_m^0 \left(\frac{\varphi_m}{\varphi_m^0}\right)^n \left(\frac{1-\varphi}{1-\varphi_0}\right)^{n'} \left(\frac{\varphi_{m'}^c - \varphi_{m'}}{\varphi_{m'}^c}\right)^{n''},\]

for \(\varphi_{m'} < \varphi_{m'}^c\), and

(40)\[a_m = 0,\]

if \(\varphi_{m'}(t) \geq \varphi_{m'}^c\), where \(\varphi_{m'}^c\) represents the critical volume fraction necessary for complete blocking of the reaction of mineral \({{\mathcal M}}_m\).


Sorption reactions incorporated into PFLOTRAN consist of ion exchange and surface complexation reactions for both equilibrium and multirate formulations.

Ion Exchange

Ion exchange reactions may be represented either in terms of bulk- or mineral-specific rock properties. Changes in bulk sorption properties can be expected as a result of mineral reactions. However, only the mineral-based formulation enables these effects to be captured in the model. The bulk rock sorption site concentration \(\omega_{{\alpha}}\), in units of moles of sites per bulk sediment volume (mol/dm:math:^3), is related to the bulk cation exchange capacity \(Q_{{\alpha}}\) (mol/kg) by the expression

(41)\[\omega_{{\alpha}}= \frac{N_{\rm site}}{V} = \frac{N_{\rm site}}{M_s} \frac{M_s}{V_s} \frac{V_s}{V} = (1-\phi) \rho_s Q_{{\alpha}},\]

with solid density \(\rho_s\) and porosity \(\varphi\). The cation exchange capacity associated with the \(m\)th mineral is defined on a molar basis as

(42)\[\omega_m^{\rm CEC} = \frac{N_m}{V} = \frac{N_m}{M_m} \frac{M_m}{V_m} \frac{V_m}{V} = Q_m^{\rm CEC} \rho_m \phi_m.\]

The site concentration \(\omega_{{\alpha}}\) is related to the sorbed concentrations \(S_k^{{\alpha}}\) by the expression

(43)\[\omega_{{\alpha}}^{} = \sum_k z_k^{} S_k^{{\alpha}}.\]

In PFLOTRAN ion exchange reactions are expressed in the form with the reference cation denoted by \({{\mathcal A}}_j^{z_j+}\) on the right-hand side of the reaction

(44)\[z_j^{} {{\mathcal A}}_i^{z_i+} + z_i^{} (\chi_{{\alpha}})_{z_j} {{\mathcal A}}_j {~\rightleftharpoons~}z_i^{} {{\mathcal A}}_j^{z_j+} + z_j^{} (\chi_{{\alpha}})_{z_i} {{\mathcal A}}_i,\]

with valencies \(z_j\), \(z_i\) of cations \({\mathcal A}_j^{z_j+}\) and \({\mathcal A}_i^{z_i+}\), respectively, and exchange site \(\chi_{{\alpha}}^-\) on the solid surface. The cations \({{\mathcal A}}_i^{z_i+}, \,i\ne j\) represents all other cations besides the reference cation. The corresponding mass action equation is given by

(45)\[K_{ij}^{{\alpha}}= \left(\frac{\lambda_i^{{\alpha}}X_i^{{\alpha}}}{a_i}\right)^{z_j} \left(\frac{a_j}{\lambda_j^{{\alpha}}X_j^{{\alpha}}}\right)^{z_i},\]

with selectivity coefficient \(K_{ij}^{{\alpha}}\), solid phase activity coefficients \(\lambda_l^{{\alpha}}\) (taken as unity in what follows), and mole fraction \(X_l^{{\alpha}}\) of the \(l\)th cation of site \({{\alpha}}\). For \(N_c\) cations participating in exchange reactions, there are \(N_c-1\) independent reactions and thus \(N_c-1\) independent selectivity coefficients.

The exchange reactions may also be expressed as half reactions in the form

(46)\[z_j^{} \chi_{{\alpha}}^- + {{\mathcal A}}_j^{z_j+} {~\rightleftharpoons~}(\chi_{{\alpha}})_{z_j} {{\mathcal A}}_j^{},\]

with corresponding selectivity coefficient \(k_j^{{\alpha}}\). The half-reaction selectivity coefficients are related to the \(K_{ij}^{{\alpha}}\) by

(47)\[\log K_{ij}^{{\alpha}}= z_j^{} \log k_i^{{\alpha}}- z_i^{} \log k_j^{{\alpha}},\]


(48)\[K_{ij}^{{\alpha}}= \frac{(k_i^{{\alpha}})^{z_j}}{(k_j^{{\alpha}})^{z_i}}.\]

This relation is obtained by multiplying the half reaction for cation \({{\mathcal A}}_j\) by the valence \(z_i\) and subtracting from the half reaction for \({{\mathcal A}}_i\) multiplied by \(z_j\), resulting in cancelation of the empty site \(X^{{\alpha}}\), to obtain the complete exchange reaction (44). It should be noted that the coefficients \(k_l^{{\alpha}}\) are not unique since, although there are \(N_c\) coefficients in number, only \(N_c-1\) are independent and one may be chosen arbitrarily, usually taken as unity. Thus for \(k_j^{{\alpha}}=1\), Eqn. (48) yields

(49)\[k_i^{{\alpha}}= \big(K_{ij}^{{\alpha}}\big)^{1/z_j}.\]

An alternative form of reactions (44) often found in the literature is

(50)\[\frac{1}{z_i} \,{{\mathcal A}}_i + \frac{1}{z_j}\, (\chi_{{\alpha}})_{z_j} {{\mathcal A}}_j {~\rightleftharpoons~}\frac{1}{z_j} \,{{\mathcal A}}_j + \frac{1}{z_i}\, (\chi_{{\alpha}})_{z_i} {{\mathcal A}}_i,\]

obtained by dividing reaction (44) through by the product \(z_i z_j\). The mass action equations corresponding to reactions (50) have the form

(51)\[{\widetilde K}_{ij}^{{\alpha}}= \frac{({\widetilde k}_i^{{\alpha}})^{1/z_i}}{({\widetilde k}_j^{{\alpha}})^{1/z_j}} = \left(\frac{a_j}{X_j^{{\alpha}}}\right)^{1/z_j} \left(\frac{X_i^{{\alpha}}}{a_i}\right)^{1/z_i}.\]

The selectivity coefficients corresponding to the two forms are related by the expression

(52)\[{\widetilde K}_{ij}^{{\alpha}}= \left(K_{ij}^{{\alpha}}\right)^{1/(z_i z_j)},\]

and similarly for \(k_i^{{\alpha}}\), \(k_j^{{\alpha}}\). When comparing with other formulations it is important that the user determine which form of the ion exchange reactions are being used and make the appropriate transformations.

The selectivity coefficients satisfy the relations

(53)\[K_{ji}^{{\alpha}}= \big(K_{ij}^{{\alpha}}\big)^{-1},\]

and from the identity

(54)\[\left(\frac{X_i^{{\alpha}}}{a_i}\right)^{z_j}\left(\frac{a_j}{X_j^{{\alpha}}}\right)^{z_i} = \left[ \left(\frac{X_i^{{\alpha}}}{a_i}\right)^{z_l} \left(\frac{a_l}{X_l^{{\alpha}}}\right)^{z_i} \right]^{z_j/z_l} \left[ \left(\frac{X_l^{{\alpha}}}{a_l}\right)^{z_j}\left(\frac{a_j}{X_j^{{\alpha}}}\right)^{z_l} \right]^{z_i/z_l},\]

the following relation is obtained

(55)\[K_{ij}^{{\alpha}}= \big(K_{il}^{{\alpha}}\big)^{z_j/z_l}\big(K_{lj}^{{\alpha}}\big)^{z_i/z_l}.\]

To see how the selectivity coefficients change when changing the reference cation from \({{\mathcal A}}_j^{z_j+}\) to \({{\mathcal A}}_k^{z_k+}\) note that

(56)\[\widetilde K_{jk}^{{\alpha}}= \big(\widetilde K_{kj}^{{\alpha}}\big)^{-1},\]


(57)\[\widetilde K_{ik}^{{\alpha}}= \widetilde K_{ij}^{{\alpha}}\, \widetilde K_{jk}^{{\alpha}}.\]

This latter relation follows from adding the two reactions

(58)\[\begin{split}\frac{1}{z_i} \,{{\mathcal A}}_i + \frac{1}{z_j}\, (\chi_{{\alpha}})_{z_j} {{\mathcal A}}_j &{~\rightleftharpoons~}\frac{1}{z_j} \,{{\mathcal A}}_j + \frac{1}{z_i}\, (\chi_{{\alpha}})_{z_i} {{\mathcal A}}_i,\\ \frac{1}{z_j} \,{{\mathcal A}}_j + \frac{1}{z_k}\, (\chi_{{\alpha}})_{z_k} {{\mathcal A}}_k &{~\rightleftharpoons~}\frac{1}{z_k} \,{{\mathcal A}}_k + \frac{1}{z_j}\, (\chi_{{\alpha}})_{z_j} {{\mathcal A}}_j,\end{split}\]

to give

(59)\[\frac{1}{z_i} \,{{\mathcal A}}_i + \frac{1}{z_k}\, (\chi_{{\alpha}})_{z_k} {{\mathcal A}}_k {~\rightleftharpoons~}\frac{1}{z_k} \,{{\mathcal A}}_k + \frac{1}{z_i}\, (\chi_{{\alpha}})_{z_i} {{\mathcal A}}_i,\]

with \({{\mathcal A}}_k^{z_k+}\) as reference species.

In terms of the selectivity coefficients \(K_{ij}^{{\alpha}}\) it follows that

(60)\[\big(K_{ik}^{{\alpha}}\big)^{1/(z_i z_k)} = \big(K_{ij}^{{\alpha}}\big)^{1/(z_i z_j)} \big(K_{jk}^{{\alpha}}\big)^{1/(z_j z_k)},\]


(61)\[K_{ik}^{{\alpha}}= \big(K_{ij}^{{\alpha}}\big)^{z_k /z_j} \big(K_{jk}^{{\alpha}}\big)^{z_i/ z_j}.\]

In terms of the coefficients \(k_i^{{\alpha}}\) and \(\overline k_i^{{\alpha}}\) corresponding to reference cation \({{\mathcal A}}_k\) the transformation becomes

(62)\[\frac{\big(\overline k_i^{{\alpha}}\big)^{z_k}}{\big(\overline k_i^{{\alpha}}\big)^{z_i}} = \left(\frac{\big(k_i^{{\alpha}}\big)^{z_j}}{\big(k_i^{{\alpha}}\big)^{z_j}}\right)^{z_k/z_j} \left(\frac{\big(k_j^{{\alpha}}\big)^{z_k}}{\big(k_k^{{\alpha}}\big)^{z_j}}\right)^{z_i/z_j}.\]

In terms of the coefficients \(k_l^{{\alpha}}\) the sorbed concentration for the \(i\)th cation can be expressed as a function of the reference cation from the mass action equations according to

(63)\[X_i^{{\alpha}}= k_i^{{\alpha}}a_i^{} \left(\frac{X_j^{{\alpha}}}{k_j^{{\alpha}}a_j^{}}\right)^{z_i/z_j}.\]

For a given reference species \({{\mathcal A}}_{J_0}\) the coefficients \(K_{iJ_0}\) are uniquely determined. For some other choice of reference species, say \({{\mathcal A}}_{I_0}\), the coefficients \(K_{iI_0}\) are related to the original coefficients by the expressions

(64)\[\begin{split}\log K_{J_0I_0} &= -\log K_{I_0J_0},\\\end{split}\]

Taking the reference cation as \({{\mathcal A}}_j\) then \(k_i^{{\alpha}}\) is given by

(65)\[\begin{split}k_i^{{\alpha}}&= \big(K_{ij}^{{\alpha}}(k_j^{{\alpha}})^{z_i}\big)^{1/z_j},\\ &= (K_{ij}^{{\alpha}})^{1/z_j}, \ \ \ \ \ \ \ \ \ \ \ \ (k_j^{{\alpha}}=1),\\ &= K_{ij}^{{\alpha}}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (z_j=1).\end{split}\]

As an example consider the ion-exchange reactions with Ca\(^{2+}\) as reference species

(66)\[\begin{split}\rm 2 \, Na^+ + \chi_2 Ca &{~\rightleftharpoons~}\rm Ca^{2+} + 2 \, \chi Na,\\ \rm 2\,Mg^{2+} + 2\,\chi_2 Ca &{~\rightleftharpoons~}\rm 2\,Ca^{2+} + 2\,\chi_2 Mg,\end{split}\]

with selectivity coefficients \(K_{\rm NaCa}\) and \(K_{\rm MgCa}\). Alternatively, using Na\(^+\) as reference cation gives

(67)\[\begin{split}\rm Ca^{2+} + 2 \, \chi Na &{~\rightleftharpoons~}\rm 2 \, Na^+ + \chi_2 Ca,\\ \rm Mg^{2+} + 2 \, \chi Na &{~\rightleftharpoons~}\rm 2 \, Na^{+} + \chi_2 Mg,\end{split}\]

with selectivity coefficients \(K_{\rm CaNa}\) and \(K_{\rm MgNa}\). The selectivity coefficients are related by the equations

(68)\[\begin{split}\log K_{\rm CaNa} & = -\log K_{\rm NaCa},\\ \log K_{\rm MgNa} &= \frac{1}{2} \, \log K_{\rm MgCa} - \log K_{\rm NaCa}.\end{split}\]

Using the Gaines-Thomas convention, the equivalent fractions \(X_k^{{\alpha}}\) are defined by

(69)\[X_k^{{\alpha}}= \frac{z_k S_k^{{\alpha}}}{\displaystyle\sum_l z_l S_l^{{\alpha}}} = \frac{z_k}{\omega_{{\alpha}}}S_k^{{\alpha}},\]


(70)\[\sum_k X_k^{{\alpha}}= 1.\]

For equivalent exchange \((z_j=z_i=z)\), an explicit expression exists for the sorbed concentrations given by

(71)\[S_j^{{\alpha}}= \frac{\omega_{{\alpha}}}{z} \frac{k_j^{{\alpha}}\gamma_j m_j^{}}{\displaystyle\sum_l k_l^{{\alpha}}\gamma_l m_l^{}},\]

where \(m_k\) denotes the \(k\)th cation molality. This expression follows directly from the mass action equations for the sorbed cations and conservation of exchange sites.

In the more general case \((z_i\ne z_j)\) it is necessary to solve the nonlinear equation

(72)\[X_j^{{\alpha}}+ \sum_{i\ne j} X_i^{{\alpha}}= 1,\]

for the reference cation mole fraction \(X_j\). From the mass action equation Eqn. (45) it follows that

(73)\[X_i^{{\alpha}}= k_i^{{\alpha}}a_i^{} \left(\frac{X_j^{{\alpha}}}{k_j^{{\alpha}}a_j^{}}\right)^{z_i/z_j}.\]

Defining the function

(74)\[f(X_j^{{\alpha}}) = X_j^{{\alpha}}+ \sum_{i\ne j}X_i^{{\alpha}}(X_j^{{\alpha}})-1,\]

its derivative is given by

(75)\[\frac{df}{dX_j^{{\alpha}}} = 1 - \frac{1}{z_j^{} X_j^{{\alpha}}}\sum_{i\ne j} z_i^{} k_i^{{\alpha}}a_i^{} \left(\frac{X_j^{{\alpha}}}{k_j^{{\alpha}}a_j^{}}\right)^{z_i/z_j}.\]

The reference mole fraction is then obtained by Newton-Raphson iteration

(76)\[(X_j^{{\alpha}})^{k+1} = (X_j^{{\alpha}})^k -\dfrac{f[(X_j^{{\alpha}})^k]}{\dfrac{df[(X_j^{{\alpha}})^k]}{dX_j^{{\alpha}}}}.\]

The sorbed concentration for the \(j\)th cation appearing in the accumulation term is given by

(77)\[S_j^{{\alpha}}= \frac{\omega_{{\alpha}}}{z_j} X_j^{{\alpha}},\]

with the derivatives for \(j\ne l\)

(78)\[\begin{split}\dfrac{{{\partial}}S_j^{{\alpha}}}{{{\partial}}m_l} &= -\frac{\omega_{{\alpha}}}{m_l} \dfrac{X_j^{{\alpha}}X_l^{{\alpha}}}{\displaystyle\sum_l z_l X_l^{{\alpha}}},\\ &= -\frac{1}{m_l} \dfrac{z_jz_lS_j^{{\alpha}}S_l^{{\alpha}}}{\displaystyle\sum_l z_l^2 S_l^{{\alpha}}},\end{split}\]

and for \(j=l\)

(79)\[\begin{split}\dfrac{{{\partial}}S_j^{{\alpha}}}{{{\partial}}m_j} &= \frac{\omega_{{\alpha}}X_j^{{\alpha}}}{z_j m_j} \left(1-\dfrac{z_j X_j^{{\alpha}}}{\displaystyle\sum_{l} z_{l} X_{l}^{{\alpha}}}\right),\\ &= \frac{S_j^{{\alpha}}}{m_j} \left(1-\dfrac{z_j^2 S_j^{{\alpha}}}{\displaystyle\sum_{l} z_{l}^2 S_{l}^{{\alpha}}}\right).\end{split}\]

Surface Complexation

Surface complexation reactions are assumed to have the form

(80)\[\nu_{{\alpha}}>\chi_{{\alpha}}+ \sum_j\nu_{ji} {{\mathcal A}}_j {~\rightleftharpoons~}> {{\mathcal S}}_{i{{\alpha}}},\]

for the \(i\)th surface complex \(>{{\mathcal S}}_{i{{\alpha}}}\) on site \({{\alpha}}\) and empty site \(>\chi_{{\alpha}}\). As follows from the corresponding mass action equation the equilibrium sorption concentration \(S_{i{{\alpha}}}^{\rm eq}\) is given by

(81)\[S_{i{{\alpha}}}^{\rm eq}= \frac{\omega_{{\alpha}}K_i Q_i}{1+\sum_l K_lQ_l},\]

and the empty site concentration by

(82)\[S_{{\alpha}}^{\rm eq}= \frac{\omega_{{\alpha}}}{1+\sum_l K_lQ_l},\]

where the ion activity product \(Q_i\) is defined by

(83)\[Q_i= \prod_j\big(\gamma_jC_j\big)^{\nu_{ji}}.\]

The site concentration \(\omega_{{\alpha}}\) satisfies the relation

(84)\[\omega_{{\alpha}}= S_{{\alpha}}+ \sum_i S_{i{{\alpha}}},\]

and is constant. The equilibrium sorbed concentration \(S_{j{{\alpha}}}^{\rm eq}\) is defined as

(85)\[S_{j{{\alpha}}}^{\rm eq} = \sum_i \nu_{ji}^{} S_{i{{\alpha}}}^{\rm eq}= \frac{\omega_{{\alpha}}}{1+\sum_l K_lQ_l} \sum_i \nu_{ji}K_i Q_i.\]

Multirate Sorption

In the multirate model the rates of sorption reactions are described through a kinetic relation given by

(86)\[\frac{{{\partial}}S_{i{{\alpha}}}}{{{\partial}}t} = k_{{\alpha}}^{} \big(S_{i{{\alpha}}}^{\rm eq}-S_{i{{\alpha}}}\big),\]

for surface complexes, and

(87)\[\begin{split}\frac{{{\partial}}S_{{{\alpha}}}}{{{\partial}}t} &= -\sum_i k_{{\alpha}}^{} \big(S_{i{{\alpha}}}^{\rm eq}-S_{i{{\alpha}}}\big),\\ &= k_{{\alpha}}\big(S_{{\alpha}}^{\rm eq}-S_{{{\alpha}}}\big),\end{split}\]

for empty sites, where \(S_{{\alpha}}^{\rm eq}\) denotes the equilibrium sorbed concentration. For simplicity, in what follows it is assumed that \(\nu_{{\alpha}}=1\). With each site \({{\alpha}}\) is associated a rate constant \(k_{{\alpha}}\) and site concentration \(\omega_{{\alpha}}\). These quantities are defined through a given distribution of sites \(\wp({{\alpha}})\), such that

(88)\[\int_0^\infty \wp(k_{{\alpha}})dk_{{\alpha}}= 1.\]

The fraction of sites \(f_{{\alpha}}\) belonging to site \({{\alpha}}\) is determined from the relation

(89)\[f_{{\alpha}}= \int_{k_{{\alpha}}-\Delta k_{{\alpha}}/2}^{k_{{\alpha}}+\Delta k_{{\alpha}}/2} \wp(k_{{\alpha}})dk_{{\alpha}}\simeq \wp(k_{{\alpha}})\Delta k_{{\alpha}},\]

with the property that


Given that the total site concentration is \(\omega\), then the site concentration \(\omega_{{\alpha}}\) associated with site \({{\alpha}}\) is equal to

(91)\[\omega_{{\alpha}}= f_{{\alpha}}\omega.\]

An alternative form of these equations is obtained by introducing the total sorbed concentration for the \(j\)th primary species for each site defined as

(92)\[S_{j{{\alpha}}}= \sum_i \nu_{ji}S_{i{{\alpha}}}.\]

Then the transport equations become

(93)\[\frac{{{\partial}}}{{{\partial}}t}\left(\varphi \Psi_j + \sum_{{{\alpha}}}S_{j{{\alpha}}}\right) + {\boldsymbol{\nabla}}\cdot{\boldsymbol{\Omega}}_j = - \sum_m\nu_{jm}I_m.\]

The total sorbed concentrations are obtained from the equations

(94)\[\frac{{{\partial}}S_{j{{\alpha}}}}{{{\partial}}t} = k_{{\alpha}}^{} \big(S_{j{{\alpha}}}^{\rm eq}-S_{j{{\alpha}}}\big).\]

Sorption Isotherm <Under Revision>

The distribution coefficient \(\widetilde K_j^D\) [m:math:^3 kg\(^{-1}\)] is customarily defined as the ratio of sorbed to aqueous concentrations with the sorbed concentration referenced to the mass of solid as given by

(95)\[\begin{split}\widetilde K_j^D &= \frac{M_j^s/M_s}{M_j^{\rm aq}/V_l},\\ &= \frac{N_j^s/M_s}{N_j^{\rm aq}/V_l},\\ &= \frac{\widetilde S_j}{C_j} = \frac{1}{\rho_w}\frac{\widetilde S_j}{m_j},\end{split}\]

where \(M_j^s = W_j N_j^s\), \(M_j^{\rm aq}=W_j N_j^{\rm aq}\), refers to the mass and number of moles of sorbed and aqueous solute related by the formula weight \(W_j\) of the \(j\)th species, \(M_s\) refers to the mass of the solid, \(V_l\) denotes the aqueous volume, \(\widetilde S_j=N_j^s/M_s\) [mol kg\(^{-1}\)] represents the sorbed concentration referenced to the mass of solid, \(C_j=N_j^{\rm aq}/V_l\) denotes molarity, and \(m_j=C_j/\rho_w\) represents molality, where \(\rho_w\) is the density of pure water.

The distribution coefficient \(\widetilde K_j^D\) may be related to its dimensionless counterpart \(K_j^D\) [—] defined by

(96)\[K_j^D = \frac{N_i^s}{N_i^{\rm aq}} = \frac{N_i^s/V}{N_i^{\rm aq}/V}= \frac{1}{\varphi s_l}\frac{S_j}{C_j},\]

by writing

(97)\[\begin{split}K_j^D &= \frac{N_i^s}{M_s} \frac{M_s}{V_s} \frac{V_s}{V_p} \frac{V_p}{V_l} \frac{V_l}{N_i^{\rm aq}},\\ &= \rho_s \frac{1-\varphi}{\varphi s_l} \widetilde K_j^D = \frac{\rho_b}{\varphi s_l} \widetilde K_j^D,\end{split}\]

with grain density \(\rho_s=M_s/V_s\), bulk density \(\rho_b=(1-\varphi)\rho_s\), porosity \(\varphi=V_p/V\), and saturation \(s_l=V_l/V_p\).

An alternative definition of the distribution coefficient denoted by \(\widehat K_j^D\) [kg m\(^{-3}\)] is obtained by using molality to define the solute concentration and referencing the sorbed concentration to the bulk volume \(V\)

(98)\[\widehat K_j^D = \frac{N_j^s/V}{N_j^{\rm aq}/M_w} = \frac{S_j}{m_j}.\]

A sorption isotherm \(S_j\) may be specified for any primary species \({{\mathcal A}}_j\) resulting in the transport equation

(99)\[\frac{{{\partial}}}{{{\partial}}t} \varphi s_l C_j + {\boldsymbol{\nabla}}\cdot{\boldsymbol{F}}_j = -\frac{{{\partial}}S_j}{{{\partial}}t},\]

for a partially saturated medium. Substituting \(S_j=\varphi s_l K_j^D C_j\) from Eqn. (96) and introducing the retardation \({{{\mathcal R}}}_j\) gives

(100)\[\frac{{{\partial}}}{{{\partial}}t} R_j \varphi s_l C_j + {\boldsymbol{\nabla}}\cdot{\boldsymbol{F}}_j = 0,\]

with the retardation given by the alternative forms

(101)\[\begin{split}R_j &= 1 + K_j^D, \ \ \ \ \ \ (\text{dimensionless)},\\ &= 1+ \frac{\rho_b}{\varphi s_l} \widetilde K_j^D, \ \ \ \ \ \ (\text{conventional}),\\ &= 1+ \frac{1}{\varphi s_l \rho_w} \widehat K_j^D, \ \ \ \ \ \ (\text{molality-based}).\end{split}\]

Three distinct models are available for the sorption isotherm \(S_j\) in PFLOTRAN:

  • linear \(K_D\) model:

    (102)\[S_j = \varphi s_l K_j^D C_j = \widehat K_j^D m_j,\]

    with distribution coefficient \(\widehat K_j^D\),

  • Langmuir isotherm:

    (103)\[S_j= \frac{K_j^L b_j^L C_j/ \rho_w}{1+K_j^L C_j/ \rho_w} = \frac{K_j^L b_j^L m_j}{1+K_j^L m_j},\]

    with Langmuir coefficients \(K_j^L\) and \(b_j^L\), and

  • Freundlich isotherm:

    (104)\[S_j = K_j^F \left(\frac{C_j}{\rho_w}\right)^{(1/n_j^F)} = K_j^F \big(m_j\big)^{(1/n_j^F)},\]

    with coefficients \(K_j^F\) and \(n_j^F\).

Colloid-Facilitated Transport

Colloid-facilitated transport is implemented into PFLOTRAN based on surface complexation reactions. Competition between mobile and immobile colloids and stationary mineral surfaces is taken into account. Colloid filtration processes are not currently implemented into PFLOTRAN. A colloid is treated as a solid particle suspended in solution or attached to a mineral surface. Colloids may be generated through nucleation of minerals in solution, although this effect is not included currently in the code.

Three separate reactions may take place involving competition between mobile and immobile colloids and mineral surfaces

(105)\[\begin{split}\mathrm{>} X_k^{{\rm m}}+ \sum_j\nu_{jk}{{\mathcal A}}_j &{~\rightleftharpoons~} \mathrm{>} S_k^{{\rm m}},\\ \mathrm{>} X_k^{{\rm im}}+ \sum_j\nu_{jk}{{\mathcal A}}_j &{~\rightleftharpoons~} \mathrm{>} S_k^{{\rm im}},\\ \mathrm{>} X_k^s + \sum_j\nu_{jk}{{\mathcal A}}_j &{~\rightleftharpoons~} \mathrm{>} S_k^s,\end{split}\]

with corresponding reaction rates \(I_k^{{\rm m}}\), \(I_k^{{\rm im}}\), and \(I_k^s\), where the superscripts \(s\), \(m\), and \(im\) denote mineral surfaces, and mobile and immobile colloids, respectively. In addition, reaction with minerals \({{\mathcal M}}_s\) may occur according to the reaction

(106)\[\sum_j\nu_{js}{{\mathcal A}}_j {~\rightleftharpoons~}{{\mathcal M}}_s.\]

The transport equations for primary species, mobile and immobile colloids, read

(107)\[\frac{{{\partial}}}{{{\partial}}t} \varphi s_l \Psi_j^l + {\boldsymbol{\nabla}}\cdot{\boldsymbol{\Omega}}_j^l = -\sum_k\nu_{jk}\big(I_k^{{\rm m}}+ I_k^{{\rm im}}+ \sum_s I_k^s\big) - \sum_s \nu_{js} I_s,\]
(108)\[\frac{{{\partial}}}{{{\partial}}t} \varphi s_l S_k^{{\rm m}} + {\boldsymbol{\nabla}}\cdot{\boldsymbol{q}}_c S_k^{{\rm m}} = I_k^{{\rm m}},\]
(109)\[\frac{{{\partial}}}{{{\partial}}t} S_k^{{\rm im}} = I_k^{{\rm im}},\]
(110)\[\frac{{{\partial}}}{{{\partial}}t} S_k^s = I_k^s,\]

where \({\boldsymbol{q}}_c\) denotes the colloid Darcy velocity which may be greater than the fluid velocity \({\boldsymbol{q}}\). For conditions of local equilibrium the sorption reaction rates may be eliminated and replaced by algebraic sorption isotherms to yield

(111)\[\frac{{{\partial}}}{{{\partial}}t}\Big[ \varphi s_l \Psi_j^l + \sum_k \nu_{jk} \big(\varphi s_l S_k^{{\rm m}}+ S_k^{{\rm im}}+ \sum_s S_k^s\big) \Big] + {\boldsymbol{\nabla}}\cdot\Big({\boldsymbol{\Omega}}_j^l + {\boldsymbol{q}}_c \sum_k \nu_{jk} S_k^{{\rm m}}\Big) = - \sum_s \nu_{js} I_s.\]

In the kinetic case either form of the primary species transport equations given by Eqn. (107) or (111) can be used provided it is coupled with the appropriate kinetic equations Eqns. (108)(110). The mobile case leads to additional equations that must be solved simultaneously with the primary species equations. A typical expression for \(I_k^m\) might be

(112)\[I_k^m = k_k\big(S_k^m - S_{km}^{\rm eq}\big),\]

with rate constant \(k_k\) and where \(S_{km}^{\rm eq}\) is a known function of the solute concentrations. In this case, Eqn. (108) must be added to the primary species transport equations. Further reduction of the transport equations for the case where a flux term is present in the kinetic equation is not possible in general for complex flux terms.

Tracer Mean Age

PFLOTRAN implements the Eulerian formulation of solute age for a nonreactive tracer following Goode (1996). PFLOTRAN solves the advection-diffusion/dispersion equation for the mean age given by

(113)\[\frac{{{\partial}}}{{{\partial}}t} \varphi s AC + {\boldsymbol{\nabla}}\cdot\Big({\boldsymbol{q}}AC - \varphi s D {\boldsymbol{\nabla}}(AC)\Big) = \varphi s C,\]

where \(A\) denotes the mean age of the tracer with concentration \(C\). Other quantities appearing in the age equation are identical to the tracer transport equation for a partially saturated porous medium with saturation state \(s\). The age and tracer transport equations are solved simultaneously for the age-concentration \(\alpha = A C\) and tracer concentration \(C\). The age-concentration \({{\alpha}}\) satisfies the usual advection-diffusion-dispersion equation with a source term on the right-hand side.

The mean tracer age is calculated in PFLOTRAN by adding the species Tracer_Age together with Tracer to the list of primary species


including sorption through a constant \(K_d\) model if desired

      DISTRIBUTION_COEFFICIENT 500. ! kg water/m^3 bulk
      DISTRIBUTION_COEFFICIENT 500. ! kg water/m^3 bulk

and specifying these species in the initial and boundary CONSTRAINT condition as e.g.:

    Tracer     1.e-8        F
    Tracer_Age 1.e-16       F

Output is given in terms of \(\alpha\) and \(C\) from which the mean age \(A\) can be obtained as \(A= \alpha/C\).

Thermodynamic Database

PFLOTRAN reads thermodynamic data from a database that may be customized by the user. Reactions included in the database consist of aqueous complexation, mineral precipitation and dissolution, gaseous reactions, and surface complexation. Ion exchange reactions and their selectivity coefficients are entered directly from the input file. A standard database supplied with the code is referred to as hanford.dat and is found in the ./database directory in the PFLOTRAN Git repository. This database is an ascii text file that can be edited by any editor and is equivalent to the EQ3/6 database:
generated by GEMBOCHS.V2-Jewel.src.R5 03-dec-1996 14:19:25

The database provides equilibrium constants in the form of log \(K\) values at a specified set of temperatures listed in the top line of the database. A least squares fit is used to interpolate the log \(K\) values between the database temperatures using a Maier-Kelly expansion of the form

(114)\[\log K = c_{-1} \ln T + c_0 + c_1 T + \frac{c_2}{T} + \frac{c_3}{T^2},\]

with fit coefficients \(c_i\). The thermodynamic database stores all chemical reaction properties (equilibrium constant \(\log K_r\), reaction stoichiometry \(\nu_{ir}\), species valence \(z_i\), Debye parameter \(a_i\), mineral molar volume \(\overline V_m\), and formula weight \(w_i\)) used in PFLOTRAN. The database is divided into 5 blocks as listed in Table [tdatabase], consisting of database primary species, aqueous complex reactions, gaseous reactions, mineral reactions, and surface complexation reactions. Each block is terminated by a line beginning with ’null’. The quantity \(N_{\rm temp}\) refers to the number of temperatures at which log \(K\) values are stored in the database. In the hanford.dat database \(N_{\rm temp}=8\) with equilibrium constants stored at the temperatures: 0, 25, 60, 100, 150, 200, 250, and 300\(^\circ\)C. The pressure is assumed to lie along the saturation curve of pure water for temperatures above 25\(^\circ\)C and is equal to 1 bar at lower temperatures. Reactions in the database are assumed to be written in the canonical form

(115)\[{{\mathcal A}}_r {~\rightleftharpoons~}\sum_{i=1}^{\rm nspec} \nu_{ir}{{\mathcal A}}_i,\]

for species \({{\mathcal A}}_r\), where nspec refers to the number of aqueous or gaseous species \({{\mathcal A}}_i\) on the right-hand side of the reaction. Redox reactions in the standard database hanford.dat are usually written in terms of O\(_{2(g)}\). Complexation reactions involving redox sensitive species are written in such a manner as to preserve the redox state.

Primary Species: name, \(a_0\), \(z\), \(w\)
Secondary Species: name, nspec, (\(\nu\)(n), name(\(n\)), \(n\)=1, nspec), log\(K\)(1:\(N_{\rm temp}\)), \(a_0\), \(z\), \(w\)
Gaseous Species: name, \(\overline V\), nspec, (\(\nu\)(n), name(\(n\)), \(n\)=1, nspec), log\(K\)(1:\(N_{\rm temp}\)), \(w\)
Minerals: name, \(\overline V\), nspec, (\(\nu\)(n), name(\(n\)), \(n\)=1, nspec), log\(K\)(1:\(N_{\rm temp}\)), \(w\)
Surface Complexes: \(>\)name, nspec, \(\nu\), \(>\)site, (\(\nu\)(n), name(\(n\)), \(n\)=1, nspec-1), log\(K\)(1:\(N_{\rm temp}\)), \(z\), \(w\)

Note that chemical reactions are not unique. For example, given a particular mineral reaction

(116)\[\sum_j \nu_{jm} {{\mathcal A}}_j {~\rightleftharpoons~}{{\mathcal M}}_m,\]

and equally acceptable reaction is the scaled reaction

(117)\[\sum_j \lambda_m\nu_{jm} {\mathcal A}_j {~\rightleftharpoons~}\lambda_m {\mathcal M}_m,\]

with scale factor \(\lambda_m\) corresponding to a different choice of formula unit. A different scale factor may be used for each mineral. The scaled reaction corresponds to

(118)\[\sum_j \nu_{jm}' {\mathcal A}_j {~\rightleftharpoons~} {\mathcal M}_m',\]

with \({\mathcal M}_m' = \lambda_m{\mathcal M}_m\), \(\nu_{jm}' = \lambda_m\nu_{jm}\). In addition, the mineral molar volume \(\overline V_m\), formula weight \(W_m\), and equilibrium constant \(K_m\) are scaled according to

(119)\[\begin{split}\overline V_m' &= \lambda_m\overline V_m,\\ W_m' &= \lambda_m W_m,\\ \log K_m' &= \lambda_m \log K_m.\end{split}\]

The saturation index \({\rm SI}_m\) transforms according to

(120)\[{\rm SI}_m' = K_m' Q_m' = \big(K_m Q_m\big)^{\lambda_m} = ({\rm SI}_m)^{\lambda_m}.\]

Consequently, equilibrium is not affected as is to be expected. However, a more general form for the reaction rate is needed involving Temkin’s constant [see Eqn. (17)], in order to ensure that the identical solution to the reactive transport equations is obtained using the scaled reaction (Lichtner, 2016). Thus it is necessary that the following conditions hold

(121)\[\begin{split}{\overline V}_m' I_m' &= \overline V_m I_m,\\ \nu_{jm}' I_m' &= \nu_{jm} I_m.\end{split}\]

This requires that the reaction rate \(I_m\) transform as

(122)\[I_m' = \frac{1}{\lambda_m} I_m,\]

which guarantees that mineral volume fractions and solute concentrations are identical to that obtained from solving the reactive transport equations using the unscaled reaction.

From the above relations it is found that the reaction rate transforms according to

(123)\[\begin{split}I_m' &= -\frac{k_m A_m}{\lambda_m} \big(1-(K_m'Q_m')^{1/\sigma_m'}\big),\\ &= -\frac{k_m A_m}{\lambda_m} \big(1-(K_m Q_m)^{1/(\lambda_m \sigma_m)} \big),\\ &= \frac{1}{\lambda_m} I_m,\end{split}\]

where the last result is obtained by scaling Temkin’s constant according to

(124)\[\sigma_m' = \lambda_m\sigma_m.\]

It should be noted that the mineral concentration \((C_m' =({\overline V}_m^{-1})^{'} \phi_m = \lambda_m^{-1} C_m)\), differs in the two formulations; however, mass density \((\rho_m = W_m \overline V_m^{-1})\) is an invariant, unlike molar density \((\eta_m=\overline V_m^{-1})\). The scaling factor \(\lambda_m\) can be found under MINERAL_KINETICS with the option MINERAL_SCALE_FACTOR.

Eh, pe

Output for Eh and pe is calculated from the half-cell reaction

(125)\[\rm 2 \, H_2O - 4 \, H^+ - 4\,e^- {~\rightleftharpoons~}\rm O_2,\]

with the corresponding equilibrium constant fit to the Maier-Kelly expansion Eqn. (114). The fit coefficients are listed in Table below.

coefficient value
\(c_{-1}\) 6.745529048
\(c_0\) -48.295936593
\(c_1\) -0.000557816
\(c_2\) 27780.749538022
\(c_3\) 4027.337694858

Table: Fit coefficients for log \(K\) of reaction (125).