Introduction

The chemistry algorithm implemented in PFLOTRAN employs a multicomponent formulation of aqueous species, gases and minerals. The traansport equations are formulated in terms of the total concentration of a chosen set of primary or basis species. A feature of the code is that primary species which form the independent variables can be chosen arbitrarily (so long as they form an independent set of aqueous species), and need not reflect the form of the reactions as wwritten in the thermodynamic database. For example, the sets of primary species \(\{\rm K^+, Al^{3+},SiO_2,H^+\}\) and \(\{\rm K^+, AlOH_4^-,SiO_2,OH^-\}\) can be used interchangeably giving identical results.

Chemical reactions included in the code are aqueous complexing for both local equilibrium and kintic formulations (however the rate law is restricted to elementary reactions), kinetic reaction of minerals using the usual transition state rate law, gaseous reactions and various formulations of sorption. The Debye-Hueckel activity coefficient algorithm is implemented and can be invoked during various stages of stepping forward in time.

The reactive transport mode may be used either as a standalone application or sequentially coupled to a flow mode which can include both mass and heat flow. This latter capability provides feedback between flow, heat and transport allowing chemical reactions to alter material properties such as porosity, permeability, and tortuosity thereby altering the flow field. To include temperature effects on chemistry the database must provide log Ks relevant to the system at hand.

Governing Equations

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

(1)\[\frac{\partial}{\partial t}\big(\porosity \sum_{\alpha} \saturation_{\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},\]

and

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

for minerals with molar volume \(\overline{V}_m\), reaction rate \(I_m\) and volume fraction \(\porosity_m\) referenced to an REV. The term involving \(S_j\) describes sorptive processes considered in more detail below. 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 index \(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 extended Debye-Hückel activity coefficient for an aqueous electrolyte solution with ionic strength \(I\) is given by (Debye and Hückel, 1923)

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

with valence \(z_k\), ionic radius \(\stackrel{\circ}{a}_k\) in angstroms, and where the Debye-Hückel parameters \(A\), \(B\) are defined by (Helgeson and Kirkham, 1974)

(7)\[\begin{split}A &= \frac{N_A^2 e^3\sqrt{2\pi}}{\ln 10 \sqrt{1000}\big(\epsilon(T,p)RT\big)^{3/2}},\\ B &= N_A e\sqrt{\frac{8\pi}{1000 \, \epsilon(T,p) RT}} \times 10^{-8}.\end{split}\]

The \(\dot b\) term is from Helgeson (1969) given by

(8)\[\dot b = 15698.4\, T^{-1} + 41.8088 \,\ln(T) - 0.0367626 \,T - 974169.0\, T^{-2} - 268.902,\]

The quantity \(\epsilon(T,p)\) is the dielectric constant of pure water which can be found in e.g. Johnson and Norton (1991). Ionic strength \(I\) is defined as

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

with molality \(m_j\) and \(m_i\) of primary and secondary species, respectively (note: \(C_i^l = \density_l y_w^l m_i \simeq \density_l m_i\), \(\density_l\) = fluid density, \(y_w^l\) = mass fraction of \(\mathrm{H_2O}\)).

Values in CGS units used for the various constants appearing in the expressions for A and B are based on the most recent values (2020) for Avogrado’s number (\(N_A = 6.0221409 \times 10^{23}\) 1/mole), charge (\(e = 4.80320425 \times 10^{-10}\) esu), Boltzmann’s constant (\(k_B=1.38064852\times 10^{-16}\) erg/K), gas constant (\(R=8.31446261815324 \times 10^7\) erg/K/mole = \(N_A k_B\)) and \(\pi=3.14159265359\). Density of pure water is based on the IFC97 EoS. Debye-Huckel coefficients are calculated at selected temperatures along the saturation curve of pure water and linearly interpolated at intermediate temperatures.

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.

Other forms for activity coefficients exist although not currently implemented. A simplified form is given by the Davies equation

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

taking \(A = 1/2\) and \(\stackrel{\circ}{a}_k B = 1\), and \(\dot b = 0.15\) in the extended Debye-Hückel equation.

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

(11)\[{\boldsymbol{\Omega}}_j^{\alpha}= \big({\boldsymbol{q}}_{\alpha}- \porosity \saturation_{\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

(12)\[{\boldsymbol{q}}_a = -\frac{kk_{{\alpha}}}{\mu_{{\alpha}}} {\boldsymbol{\nabla}}\big(p_{{\alpha}}-\density_{{\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 \(\density_{{\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

(13)\[{\boldsymbol{D}}_{{\alpha}}= \tortuosity 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, \(\tortuosity\) 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

(14)\[\porosity = 1 - \sum_m \porosity_m.\]

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

(15)\[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

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

where \(q_M\) denotes a mass rate in units of kg/s, \(\density\) 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 \(\density_w\) according to

(17)\[C_i = \density_w m_i.\]

The activity of water is calculated from the approximate relation

(18)\[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

(19)\[I_m = -A_m\Big(\sum_l k_{ml}(T) {{{\mathcal P}}}_{ml}\Big) \Big|1-\big(K_m Q_m\big)^{1/\sigma_m}\Big|^{\beta_m} {\rm sign}(1-K_m Q_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

(20)\[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

(21)\[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 298.15\(^\circ\)K, with \(T\) in units of Kelvin, \(E_{ml}\) denotes the activation energy (J/mol), and the quantity \({{{\mathcal P}}}_{ml}\) denotes the prefactor for the \(l\)th parallel reaction with the form

(22)\[{{{\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

In the case of precipitation the mineral reaction rate can grow to unreasonable values. In such casesd it may be necessary to limit the rate so that it approaches a constant value as \(K_m Q_m \rightarrow\infty\). A rate-limited form of the mineral kinetic rate law can be devised according to the expression

(23)\[I_m^{\rm RL} = -A_m^{} \Big( \sum_l k_{ml}^{} {\mathcal P}_{ml}^{} \Big) \Bigg|\frac{1-\big(K_m Q_m\big)^{1/\sigma_m}}{1+\dfrac{1}{f_{m}^{\rm lim}}\big(K_m Q_m\big)^{1/\sigma_m}} \Bigg|^{\beta_m} {\rm sign}(1-K_m Q_m),\]

with rate-limiter \(f_{m}^{\rm lim}\). In the limit \(K_m Q_m\rightarrow\infty\), the rate becomes

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

Defining the affinity factor

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

or

(26)\[K_m Q_m = \Big(1-\Omega_m\Big)^{\sigma_m},\]

the rate may be expressed alternatively as

(27)\[I_m^{\rm RL} = -A_m^{} \Big(\sum_l k_{ml}^{} {\mathcal P}_{ml}^{} \Big) \left|\frac{\Omega_m}{1+\frac{1}{f_m^{\rm lim}} \big(1-\Omega_m\big)}\right|^{\beta_m} {\rm sign}(1-K_m Q_m).\]

Changes in Material Properties

Permeability, tortuosity and mineral surface area may be updated optionally due to mineral precipitation and dissolution reactions through the change in porosity

(28)\[\porosity = 1-\sum_m\porosity_m.\]

In this case, the initial porosity is calculated as

(29)\[\porosity_0 = 1-\sum_m\porosity_m^0,\]

where the super/subscript 0 denotes initial values.

Change in permeability involves a phenomenological relation with porosity

(30)\[k = k_0 f(\porosity,\,\porosity_0,\,\porosity_c,\,n),\]

with

(31)\[f = \left(\frac{\porosity-\porosity_c}{\porosity_0-\porosity_c}\right)^n,\]
(32)\[= f_{\rm min} \ \ \ \text{if} \ \ \ \porosity \leq \porosity_c,\]
(33)\[\tortuosity = \tortuosity_0 \left(\frac{\porosity}{\porosity_0}\right)^b,\]

and

(34)\[A_m = A_m^0 \left(\frac{\porosity_m}{\porosity_m^0}\right)^n \left(\frac{1-\porosity}{1-\porosity_0}\right)^{n'},\]

with a typical value for \(n\) of \(2/3\) reflecting the surface to volume ratio. Note that this relation only applies to primary minerals \((\porosity_m^0 > 0)\). The quantity \(\porosity_c\) refers to a critical porosity below which the permeability is assumed to be constant with scale factor \(f_{\rm min}\).

The two-thirds power arises from the assumption that the number of reacting mineral grains contained in a REV remains constant. To see this consider cubical grains with the length of a side denoted by \(\ell_m\) (note that spheres could also be used without changing the result). Then the volume and surface area of an individual grain are given by

(35)\[v_m = \ell_m^3,\]

and

(36)\[a_m = 6 \ell_m^2.\]

The mineral volume fraction can be written in terms of the grain size as

(37)\[\porosity_m = \frac{V_m}{V} = \frac{N_m v_m}{V} = \eta_m \ell_m^3,\]

where the grain density given by

(38)\[\eta_m = \frac{N_m}{V}\]

is assumed to be constant. It follows that solving for \(\ell_m\) gives

(39)\[\ell_m = \left(\frac{\porosity_m}{\eta_m}\right)^{1/3},\]

and thus squaring yields

(40)\[\ell_m^2 = \left(\frac{\porosity_m}{\eta_m}\right)^{2/3}.\]

Therefore the mineral surface area \(A_m\) is given by

(41)\[A_m = \eta_m a_m = 6 \eta_m \ell_m^{2} = 6 \eta_m \left(\frac{\porosity_m}{\eta_m}\right)^{2/3} = 6 \eta_m^{1/3} \porosity_m^{2/3}.\]

A similar expression can be written for the initial surface area

(42)\[A_m^0 = 6 \eta_m \left(\frac{\porosity_m^0}{\eta_m}\right)^{2/3},\]

using the same grain density \(\eta_m\) by assumption. Taking their ratio then gives the desired result

(43)\[A_m = A_m^0 \left(\frac{\porosity_m}{\porosity_m^0}\right)^{2/3},\]

which is independent of the grain density. It should be noted, however, that this result only applies to primary minerals because of the restriction \(\porosity_m^0 > 0\). For secondary minerals, or a primary mineral which has completely dissolved at a grid cell, Eqn. (41) must be used (This formulation is currently not implemented in PFLOTRAN).

In PFLOTRAN the solid is represented as an aggregate of minerals described quantitatively by specifying its porosity \(\porosity\) and the volume fraction \(\porosity_m = V_m/V\) of each primary mineral referenced to the bulk volume \(V\) of the porous medium. It is not necessary that Eqn. (28) relating porosity and mineral volume fractions holds, and often the porosity is kept constant during the simulation.

An alternative formulation of the mineral volume fraction is to specify it relative to the total mineral volume rather than the bulk volume

(44)\[\hat\porosity_m = \frac{V_m}{V_s} = \frac{V_m}{\sum_{m'} V_{m'}},\]

with \(\sum_m\hat\porosity_m=1\). The two formulations are related by the porosity as given by

(45)\[\porosity_m = (1-\porosity) \hat\porosity_m.\]

The solid composition may also be specified by giving the mass or mole fractions \(y_m, x_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

(46)\[\begin{split}\porosity_m &= (1-\porosity) \frac{x_m \overline V_m}{\sum_{m'} x_{m'} \overline V_{m'}},\\ &= (1-\porosity) \frac{y_m^{} \density_m^{-1}}{\sum_{m'} y_{m'}^{} \density_{m'}^{-1}}.\end{split}\]

The inverse relation is given by

(47)\[x_m = \frac{\porosity_m}{\overline V_m \eta_s(1-\porosity)},\]

and similarly for the mass fraction, where

(48)\[\density_m^{} = W_m^{} \overline V_m^{-1},\]

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

(49)\[\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 \(\density_s\) by

(50)\[\density_s = W_s \eta_s,\]

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

(51)\[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

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

Variable Surface Area

An semi-analytical solution can be derived for the mineral volume fraction mass balance equation

(53)\[\frac{\partial\porosity_m}{\partial t} = \overline V_m I_m\]

for stationary state conditions. The reaction rate \(I_m\) is assumed to have the typical form based on transition state theory

(54)\[I_m = - k_m A_m \Omega_m\]

with affinity factor \(\Omega_m = 1-K_m Q_m\) assumed to be constant. The mineral surface area \(A_m\) is assumed to be a power law function of the mineral volume fraction

(55)\[A_m = A_m^0 \left(\frac{\porosity_m}{\porosity_m^0}\right)^n,\]

with constant \(n\). The affinity factor \(\Omega_m\) is constant, for example, for a stationary state or at the inlet boundary. The mineral mass balance equation can then be written in the form

(56)\[\frac{\partial\zeta_m}{\partial t} = -\alpha_m \zeta_m^n,\]

where \(\zeta_m\) is defined as the ratio

(57)\[\zeta_m = \frac{\porosity_m}{\porosity_m^0},\]

where \(\porosity_m^0\) refers to the initial mineral volume fraction at \(t=0\) and \(\alpha_m\) is given by

(58)\[\alpha_m = \frac{\overline V_m k_m A_m^0 \Omega_m}{\porosity_m^0}.\]

The equation for \(\zeta_m\) can be solved analytically with the initial condition \(\zeta(0)=1\) to give

(59)\[\zeta_m(t) = \left(1-(1-n) \alpha_m t \right)^{1/(1-n)}, \ \ \ (n\ne 1).\]

This solution breaks down if \(n=1\), in which case one can solve for \(\zeta_m\) directly to give the exponential relation

(60)\[\zeta_m(t) = {\rm e}^{-\alpha_m t}, \ \ \ (n=1).\]

Affinity Threshold

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

Sorption

Sorption reactions incorporated into PFLOTRAN consist of specifying a sorption isotherm, ion exchange reactions, and equilibrium and multirate formulations of surface complexation reactions. Each of these is dealt with in more detail below.

Sorption Isotherm

The sorption isotherm relates the sorbed concentration at the solid surface to the aqueous concentration in contact with the solid at constant temperature. It is a function of the free ion primary species concentrations \(S_j(c_1,\,\ldots, \,c_{N_c})\) (not total conentrations). It is a phenomenological formulation as opposed to a mechanisitc one and is typically not associated with an explicit chemical reaction. Finally, note that a sorption isotherm may represent equilibrium or kinetic processes depending on the data used to fit the isotherm.

The sorption isotherm appears as a source/sink term in the transport equations as given by

(61)\[\frac{\partial}{\partial t} \porosity \saturation_l \Psi_j + \vec\nabla\cdot\vec\Omega_j = -\frac{\partial S_j}{\partial t},\]

with saturation \(s_l\). Combining time derivative terms the transport equations become

(62)\[\frac{\partial}{\partial t} \big(\porosity \saturation_l\Psi_j + S_j \big) + \vec\nabla\cdot\vec\Omega_j = 0,\]

This equation can be rewritten as

(63)\[\frac{\partial}{\partial t} \Big[\porosity \saturation_l\Psi_j R_j \Big] + \vec\nabla\cdot\vec\Omega_j = 0,\]

where the local retardation factor \(R_j\) is defined in terms of the distribution coefficient \(K_j^D\) as

(64)\[\begin{split}R_j &= 1 + K_j^D,\\ K_j^D &= \frac{S_j}{\porosity \saturation_l\Psi_j}.\end{split}\]

For the case when \(R_j\) = constant, the transport equation can be written in the form

(65)\[\frac{\partial}{\partial t} \Big[\porosity \saturation_l\Psi_j\Big] + \vec\nabla\cdot\frac{1}{R_j}\vec\Omega_j = 0,\]

resulting in retarded advective and diffusive/dispersive transport. Note that the retardation varies inversely with the total concentration, not the free ion concentration, and thus aqueous complexing reactions lead to a reduction in the retardation. As a consequence strong complexing can reduce significantly the retardation coefficient compared to the value obtained using the free ion concentration.

The distribution coefficient \(\tilde K_j^D\) [m\(^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

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

where \(S_j^M = 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, \(\tilde 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/\density_w\) represents molality, where \(\density_w\) is the density of pure water.

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

(67)\[K_j^D = \frac{N_j^s}{N_j^{\rm aq}} = \frac{N_j^s/V}{N_j^{\rm aq}/V}= \frac{1}{\porosity \saturation_l}\frac{S_j}{C_j},\]

by writing

(68)\[\begin{split}K_j^D &= \frac{N_j^s}{M_s} \frac{M_s}{V_s} \frac{V_s}{V_p} \frac{V_p}{V_l} \frac{V_l}{N_j^{\rm aq}},\\ &= \density_s \frac{1-\porosity}{\porosity \saturation_l} \tilde K_j^D = \frac{\density_b}{\porosity \saturation_l} \tilde K_j^D,\end{split}\]

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

An alternative definition of the distribution coefficient denoted by \(\hat 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\)

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

The local retardation coefficient \(R_j\) can be expressed in the alternative forms

(70)\[\begin{split}R_j &= 1 + K_j^D, \ \ \ \ \ \ (\text{dimensionless)},\\ &= 1+ \frac{\density_b}{\porosity \saturation_l} \tilde K_j^D, \ \ \ \ \ \ (\text{conventional}),\\ &= 1+ \frac{1}{\porosity \saturation_l \density_w} \hat 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:

    (71)\[S_j = \porosity \saturation_l K_j^D C_j = \hat K_j^D m_j,\]

    with distribution coefficient \(\hat K_j^D\).

  • Langmuir isotherm:

    (72)\[S_j= \frac{K_j^L b_j^L C_j/ \density_w}{1+K_j^L C_j/ \density_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\).

  • Freundlich isotherm:

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

    with coefficient \(K_j^F\) and inverse power \(n_j^F\).

Ion Exchange

In PFLOTRAN ion exchange reactions are written in terms of a reference cation denoted by \({\mathcal A}_j^{z_j+}\) which appears on the right-hand side of the reaction

(74)\[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\) for cations \({\mathcal A}_j^{z_j+}\) and \({\mathcal A}_i^{z_i+}\), respectively, and exchange site \(\chi_{{\alpha}}^-\) of type \(\alpha\) on the solid surface. The cations \({{\mathcal A}}_i^{z_i+}, \,i\ne j\) represent all other cations besides the reference cation. The corresponding mass action equation is given by

(75)\[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 on 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

(76)\[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

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

or

(78)\[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^{z_j+}\) by the valence \(z_i\) and subtracting from the half reaction for \({\mathcal A}_i^{z_i+}\) multiplied by \(z_j\), resulting in cancelation of the empty site \(\chi_{\alpha}^-\), to obtain the complete exchange reaction (74). 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. (78) yields

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

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

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

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

(81)\[{\tilde K}_{ij}^{\alpha}= \frac{({\tilde k}_i^{{\alpha}})^{1/z_i}}{({\tilde 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

(82)\[{\tilde 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 governing equations incorporating homogeneous aqueous complexing reactions combined with ion exchange reactions with reaction rates \(\Gamma_{ji}\) and with reference cation \({\mathcal A}_j\) have the form

(83)\[\begin{split}\frac{\partial}{\partial t } \porosity \Psi_j + \vec\nabla\cdot\vec\Omega_j &= \sum_{i\ne j} z_i \Gamma_{ji},\\ \frac{\partial}{\partial t } \porosity \Psi_i + \vec\nabla\cdot\vec\Omega_i &= -z_j \Gamma_{ji},\\ \frac{\partial S_j}{\partial t} &= -\sum_{i\ne j} z_i \Gamma_{ji},\\ \frac{\partial S_i}{\partial t} &= z_j \Gamma_{ji}.\end{split}\]

The ion exchange reaction rates may be eliminated from the aqueous transport equations to yield

(84)\[\begin{split}\frac{\partial}{\partial t } \porosity \Psi_j + \vec\nabla\cdot\vec\Omega_j &= -\frac{\partial S_j}{\partial t},\\ \frac{\partial}{\partial t } \porosity \Psi_i + \vec\nabla\cdot\vec\Omega_i &= -\frac{\partial S_i}{\partial t}.\end{split}\]

Assuming conditions of local equilibrium the ion exchange reaction rates may be eliminated and replaced by isotherms.

It can be easily demonstrated that the governing equations conserve the exchange site density \(\omega\) given by

(85)\[\omega = z_j S_j + \sum_{i\ne j} z_i S_i,\]

assuming material properties are not altered by mineral precipitation/dissolution reactions. It follows that

(86)\[\begin{split}\frac{\partial\omega(\vec r, \, t)}{\partial t} &= z_j \sum_{i\ne j} z_i \Gamma_{ji} - z_j \sum_{i \ne j} z_i \Gamma_{ji},\\ &=0.\end{split}\]

Since charge is conserved by the ion exchange reactions, the transport equations coupled to ion exchange must also conserve charge and as a result no additional constraints are needed.

Exchange Capacity

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 rock volume (mol/dm\(^3\)), is related to the bulk cation exchange capacity \(Q_{\alpha}\) (mol/kg) by the expression

(87)\[\omega_{{\alpha}}= \frac{N_{\rm site}}{V} = \frac{N_{\rm site}}{M_s} \frac{M_s}{V_s} \frac{V_s}{V} = (1-\porosity) \density_s Q_{{\alpha}},\]

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

(88)\[\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} \density_m \porosity_m.\]

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

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

Selectivity Coefficient Relations

The selectivity coefficients satisfy the relations

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

and from the identity

(91)\[\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

(92)\[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

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

and

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

This latter relation follows from adding the two reactions

(95)\[\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

(96)\[\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 cation.

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

(97)\[\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)},\]

or

(98)\[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

(99)\[\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

(100)\[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 cation \({\mathcal A}_{J_0}\) the coefficients \(K_{iJ_0}\) are uniquely determined. For some other choice of reference cation, say \({\mathcal A}_{I_0}\), the coefficients \(K_{iI_0}\) are related to the original coefficients by the expressions

(101)\[\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

(102)\[\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 cation

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

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

(104)\[\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

(105)\[\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}\]

Gaines-Thomas Exchange

The Gaines-Thomas convention (Gaines and Thomas, 1953), is based on the equi-valent fractions \(X_k^{{\alpha}}\) defined by

(106)\[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}},\]

with

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

The index \(\alpha\) refers to distinct exchange sites.

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

(108)\[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

(109)\[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. (75) it follows that

(110)\[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

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

its derivative is given by

(112)\[\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

(113)\[(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

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

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

(115)\[\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\)

(116)\[\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

(117)\[\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

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

and the empty site concentration by

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

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

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

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

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

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

(122)\[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

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

for surface complexes, and

(124)\[\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

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

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

(126)\[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

(127)\[\sum_{{\alpha}}f_{{\alpha}}=1.\]

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

(128)\[\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

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

Then the transport equations become

(130)\[\frac{{{\partial}}}{{{\partial}}t}\left(\porosity \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

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

Aqueous Complexing Reaction Kinetics

PFLOTRAN allows the user to input kinetic reactions for homogeneous aqueous complexing reactions through the GENERAL_REACTION keyword. The reactions are treated as being elementary reactions with reaction rate expressions derived from the law of mass action. Use the sandbox for more general kinetic rate laws not limited to elementary reactions based on the law of mass action.

To develop the governing equations for this system, reactions are written for intrinsically fast and slow reactions corresponding to local equilibrium and kinetic rates of reaction according to

(132)\[\begin{split}\sum_j \nu_{ji}^{leq} {\mathcal A}_j &\rightleftharpoons {\mathcal A}_i, \ \ \ (\text{fast}),\\ \emptyset &\rightleftharpoons \sum_j \nu_{jr}^{kin} {\mathcal A}_j, \ \ \ (\text{slow}).\end{split}\]

The sums are over a set of independent primary species. In the expression for kinetic reactions all species are brought to the right-hand side with reactants having negative stoichiometric coefficients and products positive coefficients. The reaction rates corresponding to fast reactions are eliminated from the transport equations and replaced by algebraic mass action relations.

The kinetic rate expression is assumed to have the form of the difference between forward and backward reaction rates proportional to the product of concentrations of reactants and products, respectively, raised to the power of their stochiometric coefficients

(133)\[\Gamma_r = k_r^+ \prod_{\nu_{jr}^{kin}<0} (a_j)^{-\nu_{jr}^{kin}} - k_r^- \prod_{\nu_{jr}^{kin}>0} (a_j)^{\nu_{jr}^{kin}}.\]

At equilibrium \(\Gamma_r=0\) and the equilibrium mass action equation is retrieved

(134)\[K_r = \frac{k_r^+}{k_r^-} = \prod_j a_j^{\nu_{jr}^{kin}},\]

with the equilibrium constant \(K_r\) equal to the ratio of the forward to backward rate constants.

With the above reactions the transport equations for primary species have the form (including precipitation/disollution reactions with rates \(\Gamma_m\))

(135)\[\frac{\partial}{\partial t} \porosity \Psi_j + \vec\nabla\cdot\vec\Omega_j = \sum_r \nu_{jr}^{kin} \Gamma_r -\sum_m \nu_{jm} \Gamma_m,\]

where \(\Psi_j\) and \(\vec\Omega_j\) are the total concentration and flux, respectively, defined as

(136)\[\begin{split}\Psi_j = c_j + \sum_i \nu_{ji}^{leq} c_i,\\ \vec\Omega_j = \vec F_j + \sum_i \nu_{ji}^{leq} \vec F_i,\end{split}\]

where \(\vec F_k\) is the individual species flux consisting of contributions from advection, diffusion and dispersion, and the secondary species concentrations \(c_i\) are given by the mass action law

(137)\[c_i = \frac{K_i}{\gamma_i} \prod_j \big(\gamma_j c_j\big)^{\nu_{ji}^{leq}},\]

relating secondary species concentrations to primary species. Thus in this formulation the reaction rates for intrinsically fast reactions are replaced by mass action equations thereby reducing the number of partial differential equations that are necessary to solve.

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

(138)\[\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

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

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

(140)\[\frac{{{\partial}}}{{{\partial}}t} \porosity \saturation_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,\]
(141)\[\frac{{{\partial}}}{{{\partial}}t} \porosity \saturation_l S_k^{{\rm m}} + {\boldsymbol{\nabla}}\cdot{\boldsymbol{q}}_c S_k^{{\rm m}} = I_k^{{\rm m}},\]
(142)\[\frac{{{\partial}}}{{{\partial}}t} S_k^{{\rm im}} = I_k^{{\rm im}},\]
(143)\[\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

(144)\[\frac{{{\partial}}}{{{\partial}}t}\Big[ \porosity \saturation_l \Psi_j^l + \sum_k \nu_{jk} \big(\porosity \saturation_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. (140) or (144) can be used provided it is coupled with the appropriate kinetic equations Eqns. (141)(143). 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

(145)\[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. (141) 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

(146)\[\frac{{{\partial}}}{{{\partial}}t} \porosity \saturation AC + {\boldsymbol{\nabla}}\cdot\Big({\boldsymbol{q}}AC - \porosity \saturation D {\boldsymbol{\nabla}}(AC)\Big) = \porosity \saturation 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

PRIMARY_SPECIES
  Tracer
  Tracer_Age
/

Sorption may be included through a constant \(K_d\) model if desired.

SORPTION
  ISOTHERM_REACTIONS
    Tracer
      TYPE LINEAR
      DISTRIBUTION_COEFFICIENT 500. ! kg water/m^3 bulk
    /
    Tracer_Age
      TYPE LINEAR
      DISTRIBUTION_COEFFICIENT 500. ! kg water/m^3 bulk
    /
  /
/

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

CONSTRAINT initial
  CONCENTRATIONS
    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

Database Structure

PFLOTRAN reads thermodynamic data from a database file 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:

data0.com.V8.R6
CII: GEMBOCHS.V2-EQ8-data0.com.V8.R6
THERMODYNAMIC 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

(147)\[\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 form

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

as a dissasociation reaction 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, \(v\), nspec, (\(\nu\)(n), name(\(n\)), \(n\)=1, nspec), log\(K\)(1: \(N_{\rm temp}\)), \(w\)

Minerals:

name, \(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\)

The quantities: name, \(>\)name, \(a_0\), \(z\), \(w\), \(\nu\), \(\log K\), and \(v\) refer, respectively, to the aqueous or gas species, mineral or surface complex, Debye-Hueckel radius parameter, charge, formula weight [g/mol], stoichiometric coefficient, logarithm of the equilibrium constant to base 10, and molar volume [cm\(^3\)/mol].

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

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

an equally acceptable reaction is the scaled reaction

(150)\[\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

(151)\[\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

(152)\[\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

(153)\[{\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. (19)], 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

(154)\[\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

(155)\[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

(156)\[\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

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

It should be noted that the mineral concentration \((C_m' =({\overline V}_m^{-1})^{'} \porosity_m = \lambda_m^{-1} C_m)\), differs in the two formulations; however, mass density \((\density_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

(158)\[\rm 2 \, H_2O - 4 \, H^+ - 4\,e^- \rightleftharpoons \rm O_2,\]

with the corresponding equilibrium constant fit to the Maier-Kelly expansion Eqn. (147). 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 (158).

Python Script to Select Primary and Secondary Species from Thermodynamic Database

A python script is available to help the user extract secondary species, gases and minerals from the thermodynamic database for a given set of primary species. Surface complexation reactions are not included. The python script can be found in ./tools/contrib/sec_species/rxn.py in the PFLOTRAN Git repository. The current implementation is based on the hanford.dat database. Input files are aq_sec.dat, gases.dat and minerals.dat. In addition, for each of these files there is a corresponding file containing a list of species to be skipped: aq_skip.dat, gas_skip.dat and min.dat. Before running the script it is advisable to copy the entire directory sec_species to the local hard drive to avoid conflicts when updating the PFLOTRAN repository. To run the script simply type in a terminal window:

python rxn.py

The user has to edit the rxn.py file to set the list of primary species. For example,

pri=[’Fe++’,’Fe+++’,’H+’,’H2O’]

Note that the species H2O must be include in the list of primary species. Output appears on the screen and also in the file chem.out, a listing of which appears below. The number of primary and secondary species, gases and minerals is printed out at the end of the chem.out file.

chem.out

PRIMARY_SPECIES
Fe++
Fe+++
H+
H2O
/
SECONDARY_SPECIES
O2(aq)
H2(aq)
Fe(OH)2(aq)
Fe(OH)2+
Fe(OH)3(aq)
Fe(OH)3-
Fe(OH)4-
Fe(OH)4--
Fe2(OH)2++++
Fe3(OH)4(5+)
FeOH+
FeOH++
HO2-
OH-
/
GASES
H2(g)
H2O(g)
O2(g)
/
MINERALS
Fe
Fe(OH)2
Fe(OH)3
FeO
Ferrihydrite
Goethite
Hematite
Magnetite
Wustite
/
================================================
npri =  4  nsec =  14  ngas =  3  nmin =  9

Finished!