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
and
for minerals with molar volume \(\overline{V}_m\), reaction rate \(I_m\) and volume fraction \(\varphi_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
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
yielding the mass action equations
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)
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)
The \(\dot b\) term is from Helgeson (1969) given by
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
with 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}\)).
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
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
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
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
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
The temperature dependence of the diffusion coefficient is defined through the relation
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
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
The activity of water is calculated from the approximate relation
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
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
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
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
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
with rate-limiter \(f_{m}^{\rm lim}\). In the limit \(K_m Q_m\rightarrow\infty\), the rate becomes
Defining the affinity factor
or
the rate may be expressed alternatively as
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
Change in permeability involves a phenomenological relation with porosity
with
and
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 > 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}\).
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
and
The mineral volume fraction can be written in terms of the grain size as
where the grain density given by
is assumed to be constant. It follows that solving for \(\ell_m\) gives
and thus squaring yields
Therefore the mineral surface area \(A_m\) is given by
A similar expression can be written for the initial surface area
using the same grain density \(\eta_m\) by assumption. Taking their ratio then gives the desired result
which is independent of the grain density. It should be noted, however, that this result only applies to primary minerals because of the restriction \(\phi_m^0 > 0\). For secondary minerals, or a primary mineral which has completely dissolved at a grid cell, Eqn. (40) 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 \(\varphi\) and the volume fraction \(\varphi_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
with \(\sum_m\hat\varphi_m=1\). The two formulations are related by the porosity as given by
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
The inverse relation is given by
and similarly for the mass fraction, where
and the solid molar density \(\eta_s\) is given by
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
with the mean molecular weight \(W_s\) of the solid phase equal to
Mass and mole fractions are related by the expression
Variable Surface Area¶
An semi-analytical solution can be derived for the mineral volume fraction mass balance equation
for stationary state conditions. The reaction rate \(I_m\) is assumed to have the typical form based on transition state theory
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
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
where \(\zeta_m\) is defined as the ratio
where \(\phi_m^0\) refers to the initial mineral volume fraction at \(t=0\) and \(\alpha_m\) is given by
The equation for \(\zeta_m\) can be solved analytically with the initial condition \(\zeta(0)=1\) to give
This solution breaks down if \(n=1\), in which case one can solve for \(\zeta_m\) directly to give the exponential relation
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
with saturation \(s_l\). Combining time derivative terms the transport equations become
This equation can be rewritten as
where the local retardation factor \(R_j\) is defined in terms of the distribution coefficient \(K_j^D\) as
For the case when \(R_j\) = constant, the transport equation can be written in the form
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
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/\rho_w\) represents molality, where \(\rho_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
by writing
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 \(\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\)
The local retardation coefficient \(R_j\) can be expressed in the alternative forms
Three distinct models are available for the sorption isotherm \(S_j\) in PFLOTRAN:
linear \(K_D\) model:
(70)¶\[S_j = \varphi s_l K_j^D C_j = \hat K_j^D m_j,\]with distribution coefficient \(\hat K_j^D\).
Langmuir isotherm:
(71)¶\[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\).
Freundlich isotherm:
(72)¶\[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 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
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
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
with corresponding selectivity coefficient \(k_j^{{\alpha}}\). The half-reaction selectivity coefficients are related to the \(K_{ij}^{{\alpha}}\) by
or
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 (73). 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. (77) yields
An alternative form of reactions (73) often found in the literature is
obtained by dividing reaction (73) through by the product \(z_i z_j\). The mass action equations corresponding to reactions (79) have the form
The selectivity coefficients corresponding to the two forms are related by the expression
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
The ion exchange reaction rates may be eliminated from the aqueous transport equations to yield
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
assuming material properties are not altered by mineral precipitation/dissolution reactions. It follows that
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
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
The site concentration \(\omega_{{\alpha}}\) is related to the sorbed concentrations \(S_k^{{\alpha}}\) by the expression
Selectivity Coefficient Relations¶
The selectivity coefficients satisfy the relations
and from the identity
the following relation is obtained
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
and
This latter relation follows from adding the two reactions
to give
with \({{\mathcal A}}_k^{z_k+}\) as reference cation.
In terms of the selectivity coefficients \(K_{ij}^{{\alpha}}\) it follows that
or
In terms of the coefficients \(k_i^{\alpha}\) and \(\overline k_i^{{\alpha}}\) corresponding to reference cation \({\mathcal A}_k\) the transformation becomes
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
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
Taking the reference cation as \({\mathcal A}_j\) then \(k_i^{{\alpha}}\) is given by
As an example consider the ion-exchange reactions with Ca\(^{2+}\) as reference cation
with selectivity coefficients \(K_{\rm NaCa}\) and \(K_{\rm MgCa}\). Alternatively, using Na\(^+\) as reference cation gives
with selectivity coefficients \(K_{\rm CaNa}\) and \(K_{\rm MgNa}\). The selectivity coefficients are related by the equations
Gaines-Thomas Exchange¶
The Gaines-Thomas convention (Gaines and Thomas, 1953), is based on the equi-valent fractions \(X_k^{{\alpha}}\) defined by
with
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
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
for the reference cation mole fraction \(X_j\). From the mass action equation Eqn. (74) it follows that
Defining the function
its derivative is given by
The reference mole fraction is then obtained by Newton-Raphson iteration
The sorbed concentration for the \(j\)th cation appearing in the accumulation term is given by
with the derivatives for \(j\ne l\)
and for \(j=l\)
Surface Complexation¶
Surface complexation reactions are assumed to have the form
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
and the empty site concentration by
where the ion activity product \(Q_i\) is defined by
The site concentration \(\omega_{{\alpha}}\) satisfies the relation
and is constant. The equilibrium sorbed concentration \(S_{j{{\alpha}}}^{\rm eq}\) is defined as
Multirate Sorption¶
In the multirate model the rates of sorption reactions are described through a kinetic relation given by
for surface complexes, and
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
The fraction of sites \(f_{{\alpha}}\) belonging to site \({{\alpha}}\) is determined from the relation
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
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
Then the transport equations become
The total sorbed concentrations are obtained from the equations
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
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
At equilibrium \(\Gamma_r=0\) and the equilibrium mass action equation is retrieved
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\))
where \(\Psi_j\) and \(\vec\Omega_j\) are the total concentration and flux, respectively, defined as
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
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
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
The transport equations for primary species, mobile and immobile colloids, read
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
In the kinetic case either form of the primary species transport equations given by Eqn. (139) or (143) can be used provided it is coupled with the appropriate kinetic equations Eqns. (140) – (142). 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
with rate constant \(k_k\) and where \(S_{km}^{\rm eq}\) is a known function of the solute concentrations. In this case, Eqn. (140) 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
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
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
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
an equally acceptable reaction is the scaled reaction
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
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
The saturation index \({\rm SI}_m\) transforms according to
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
This requires that the reaction rate \(I_m\) transform as
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
where the last result is obtained by scaling Temkin’s constant according to
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
with the corresponding equilibrium constant fit to the Maier-Kelly expansion Eqn. (146). 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 (157).
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!
