rxn¶

Chemical reactions between components are invoked with this control statement. It is used in conjunction with control statement trac. For facilitating the construction of the rxn input, a header describing the input is required before each group whether data is entered for that group or not (see examples) unless otherwise noted.

The header is an arbitrary text string that must be contained on a single line. Note that for components that do not react with other components, rxn is unnecessary. Specifically, conservative tracers or tracers that follow the equilibrium sorption isotherms can be modeled with just the trac macro.

Note the parameters NCPNT, NIMM, NVAP (used as indices for input) are determined by the code using information input for the trac macro. NCPNT is equal to the number of liquid components, NIMM is equal to the number of immobile components, and NVAP is equal to the number of vapor components specified in the trac macro.

• Group 1 - NCPLX, NUMRXN
• Group 2 - NGROUPS
  • GROUP (ICPNT), ICPNT = 1, NCPNT (repeated NGROUPS times, once for each group).
• Group 3 - IDCPNT, CPNTNAM, IFXCONC, CPNTPRT, CPNTGS (repeated NCPNT times)
• Group 4 - IDCPLX, CPLXNAM, CPLXPRT (repeated NCPLX times)
• Group 5 - IDIMM, IMMNAM, IMMPRT (repeated NIMM times;)
• Group 6 - IDVAP, VAPNAM, VAPPRT (repeated NVAP times)
• Group 7 - ISKIP
• Group 8 - RSDMAX
  • HEADING

Note that this is an additional heading line that precedes the LOGKEQ heading.

• Group 9 - LOGKEQ
• Group 10 - CKEQ, HEQ (NCPLX times)

or

• Group 10 - KEYWORD, NUM_TEMPS
  • EQTEMP(I), I = 1, NUM_TEMPS
  • LKEQ(I), I = 1, NUM_TEMPS

For group 10, the keyword ‘lookup’ can be used in place of CKEQ and HEQ. LOOKUP allows a lookup table to be used to describe the equilibrium constant (K) as a function of temperature. After the keyword, the user must specify the number of values of temperature and K that will be used to describe K as a function of temperature. On the next lines the temperatures, and then the K values are entered. FEHM performs a piecewise linear interpolation between the values given.

• Group 11 - STOIC(ICPNT), ICPNT = 1, NCPNT (repeated NCPLX times, once for each aqueous complex)

Input for groups 9, 10 and 11 is omitted if NCPLX, the number of aqueous complexes, is zero.

The remaining groups are entered as a unit for each kinetic reaction. If there are no kinetic reactions specified, none of the following groups (or their headers) are input. The input for groups 14 and on, depend on the kinetic reaction type specified in group 12. The input for each kinetic reaction type is described below.

• Group 12 - IDRXN
• Group 13 - JA, JB, JC (JA, JB, JC - defined on page 22)
• IDRXN = 1: Linear kinetic reaction
• Group 14 - IAQUEOUS, IIMMOBILE
• Group 15 - KD

or

• Group 15 - KEYWORD, NUM_TEMPS
  • EQTEMP(I), I = 1, NUM_TEMPS
  • TCOEFF(I), I = 1, NUM_TEMPS

• Group 16 - RATE
• IDRXN = 2: Langmuir kinetic reaction
• Group 14 - IAQUEOUS, IIMMOBILE
• Group 15 - DISTCOEEF

or

• Group 15 - KEYWORD, NUM_TEMPS
  • EQTEMP(I), I = 1, NUM_TEMPS
  • TCOEFF(I), I = 1, NUM_TEMPS

• Group 16 - RATE
• Group 17 - MAXCONC
• IDRXN = 3: General reaction
• Group 14 - NIMMOBILE, NAQUEOUS, NVAPOR
• Group 15 - KFOR, KREV
• Group 16 - IIMMOBILE (I = 1, NIMMOBILE)
• Group 17 - IMSTOIC (I = 1, NIMMOBILE)

Omit groups 16 and 17 (including headers) if NIMMOBILE is zero.

• Group 18 - IAQUEOUS (I = 1, NAQUEOUS)
• Group 19 - AQSTOIC (I = 1, NAQUEOUS)

Omit groups 18 and 19 (including headers) if NAQUEOUS is zero.

• Group 20 - IVAPOR (I = 1, NVAPOR)
• Group 21 - IVSTOIC (I = 1, NVAPOR)

Omit groups 20 and 21 (including headers) if NVAPOR is zero.

• IDRXN = 4: Dual Monod kinetics biodegradation reaction
• Group 14 - NAQUEOUS (must be >2 and < 5), NIMMOBILE
• Group 15 - SUBSTRATE, ELECACC, COMP3, COMP4, COMP5

• Group 16 - BIOMASS
• Group 17 - KS
• Group 18 - KA
• Group 19 - DECAY
• Group 20 - ELECACCSTOIC
• Group 21 - COMP3STOIC
• Group 22 - COMP4STOIC
• Group 23 - COMP5STOIC
• Group 24 - PHTHRESH
• Group 25 - QM
• Group 26 - YIELD
• Group 27 - XMINIT
• Group 28 - NBIOFRM

Omit Group 29 if NBIOFRM = 0

• Group 29 - ICBIO (I = 1, NBIOFRM)
• IDRXN=5: Radioactive Decay Reaction
• Group 14 - HALFLIFE
• Group 15 - RXNTYPE
• Group 16 - PARENT, DAUGHTER
• IDRXN=6: Kinetic Henry’s Law reaction
• Group 14 - IAQUEOUS, IVAPOR
• Group 15 - KH
• Group 16 - RATE
• IDRXN = 7: Kinetic precipitation/dissolution reaction for total component concentrations
• Group 14 - IIMMOBILE
• Group 15 - NAQUEOUS
• Group 16 - IAQUEOUS (I = 1, NAQUEOUS)
• Group 17 - IMSTOIC
• Group 18 - AQSTOIC (I = 1, NAQUEOUS)
• Group 19 - SOLUBILITY

or

• Group 19 - KEYWORD, NUM_TEMPS
  • EQTEMP(I), I = 1, NUM_TEMPS
  • TCOEFF(I), I = 1, NUM_TEMPS

• Group 20 - RATE
• Group 21 - SAREA
• IDRXN = 8: Kinetic precipitation/dissolution reaction for total component concentrations with rates based on free-ion concentrations

In general, the following examples illustrate only portions of the rxn macro and putting all of these example inputs together will not result in a working FEHM rxn macro. However, the dissolution example (the last example in this section) provides an example of a complete rxn macro and corresponds to the first example given for the trac macro. In addition, the Reactive Transport Example and the Validation Test Plan for the FEHM Application Version 2.30 (10086-VTP-2.30-00) include full example problems with input files which demonstrate the use of rxn. These input files can be used to see how rxn fits in with the other macros.

Specifically, the information in rxn must be consistent with the trac macro. For example, if a linear kinetic sorption reaction is invoked by rxn, a liquid component and solid component must be specified in the trac macro.

General Reaction Parameters. In the following example two aqueous complexes and four kinetic reactions are specified. Three liquid components, Co, Fe, and EDTA, two complexes, CoEDTA and FeEDTA, and three immobile species, Co, CoEDTA and FeEDTA, are identified. Aqueous components 1 and 3 are coupled during solution while 2 is solved independently. Note that group 6 data has been omitted since there are no vapor species in this example.

[example]

Equilibrium Reaction Parameters. Equilibrium speciation reactions modeled by FEHM can be written in the following general form:

\(\sum_{j=1}^{N_c} a_{ij}\hat{C_j} \Leftrightarrow \hat{X_i} \;\;\;\; i = 1,...,N_x\)

where \(\hat{C_j}\) is the chemical formula for the aqueous component j, and \(\hat{X_i}\) is the cchemical formula for the aqueous complex i, \(a_{ij}\) is a stoichiometric coefficient {STOIC} giving the number of moles of component j in complex i, and \(N_x\) is the number of aaqueous components.

Here is a simple example of an equilibrium speciation reaction:

\(aA + bB \Leftrightarrow A_aB_b\)

where a and b are the stoichiometric coefficients of components A and B, respectively. At equilibrium, the concentrations of A, B, and \(A_aB_b\) must satisfy the law of mass action for this reaction:

\(K = \frac{[A_aB_b]}{[A]^a[B]^b}\)

where K is the equilibrium constant {CKEQ} for the reaction and [X] is the molar concentration of species X. The total aqueous concentrations of components A and B are given by:

\(U_A = [A] + a[A_aB_b]\)

\(U_b = [B] + b[A_aB_b]\)

Therefore, the total aqueous concentrations, U, are functions of the uncomplexed component concentrations and the complex concentrations. Note that the kinetic reactions discussed in the next section are functions of the uncomplexed component concentrations and complex concentrations with the exception of the radioactive decay and precipitation/dissolution reaction in which the total aqueous concentration is used in the reaction rate law.

The following input describes the equilibrium speciation relations between the total aqueous components and aqueous complexes. Recall that the names of the components and aqueous complexes are given in Groups 3-6. The skip node option, ISKIP, is turned off and the tolerance for the speciation solver, RSDMAX, is set to 1e-9 which is the recommended value. The stoichiometry, STOIC, is specified so that the components Co and EDTA make up the first complex CoEDTA, and the components Fe and EDTA make up the second complex FeEDTA. The equilibrium constants, CKEQ, for CoEDTA and FeEDTA are 1.e18 and 6.31e27, respectively. The enthalpy change is 0.

Below is an example of using the ‘lookup’ keyword in place of CKEQ and HEQ. This example describes the temperature dependence of the equilibrium constant of complex HCO3-. Note that LOGKEQ = 1 in this example, so log K is specified in the input rather than K.

General Kinetic Parameters. In FEHM, eight kinetic reaction models are supported. Additional kinetic subroutines can be added without significant code development. The following is an example of input for the general kinetic parameters. A linear kinetic reaction, IDRXN = 1, is specified as occurring at each node in the problem (JA = 1, JB and JC = 0).

IDRXN = 1: Linear Kinetic Reaction. The retardation of contaminants due to adsorption/desorption can be modeled with a linear kinetic sorption/desorption expression. The rate of adsorption/desorption of component j is given by:

\(R_j = -k_m \left( c_j - \frac{m_j}{K_D} \right)\)

where \(c_j\) denotes the uncomplexed aqueous concentration of j {IAQUEOUS}, \(m_j\) denotes the adsorbed concentration of species j {IIMMOBILE}, \(k_m\) is the mass transfer coefficient {RATE}, and \(K_D\) is the distribution coefficient {KD}.

As \(k_m \rightarrow \infty\), this expression reduces to the linear equilibrium isotherm. The following example illustrates input for a linear kinetic reaction.

In this kinetic reaction, aqueous component 1 adsorbs to form the immobile component 3. The \(K_D\) for the reaction is 5.07 and the mass transfer coefficient is 1.0.

IDRXN = 2: Langmuir Kinetic Reaction. The Langmuir kinetic reaction rate law is given by:

\(R_j = -k_m \frac{\rho}{\theta} \left( K_L c_j \left( m_j^{MAX} - m_j \right) - m_j \right)\)

where \(k_m\) is the rate constant for desorption {RATE}, \(\rho\) is the bulk rock density {DENR}, \(\theta\) is the porosity {POR}. \(K_L\) is the distribution coefficient {DISTCOEFF}, and \(m_j^{MAX}\) is the maximum concentration that can adsorb onto the j solid {MAXCONC}. Ask \(k_m \rightarrow \infty\), this expression reduces to the Langmuir equilibrium isotherm. Example input for a Langmuir kinetic reaction follows. In this kinetic reaction, aqueous complex 101 sorbs to form immobile component 1. The distribution coefficient for the reaction is 2.e5 and the mass transfer coefficient is 0.05. The maximum sorbed concentration is 1.69e-5.

IDRXN = 3: General Kinetic Reaction. Many reactions fall under the category of the general kinetic reaction. The reaction is described by a forward rate constant {KFOR}, a reverse rate constant {KREV}, and a set of stoichiometric coefficients. The form of the general reversible reaction is given by:

\(\sum_{j=1}^{N_c+N_x+N_m+N_v} \eta_j \hat{Z}_j \Leftrightarrow \sum_{k=1}^{N_c+N_x+N_m+N_v} \eta_k' \hat{z}_k\)

where \(N_c\) is the number of aqueous components, \(N_x\) is the number of aqueous complexes { NAQUEOUS = Nc + Nx }, \(N_m\) is the number of immobile components {NIMMOBILE}, \(N_v\) is the number of vapor components {NVAPOR}, \(\eta_j\) are reactant stoichiometric coefficients, \(\eta_k'\) are product stoichiometric coefficients, and \(\hat{z}_i\) is the chemical formula for species i, which may be an uncomplexed aqueous component, aqueous complex, immobile component or vapor component.

The rate law for a general reversible reaction is given by the following expression:

\(R(z_i) = (\eta_i - \eta_i') \left[ k_f \prod_{j=1}^{N_c+N_x+N_m+N_v} z_j^{\eta_j} - k_r \prod_{k=1}^{N_c+N_x+N_m+N_v} z_k^{\eta_k'} \right]\)

where \(z_i\) is the concentration of species i.

The following is example input for a general kinetic reaction. Solid component 1 reacts to form solid components 2 and 3. The forward reaction rate is 1.26e-2 and the reverse reaction rate is 0. Therefore, this is an irreversible reaction. Note also that only solid components are reacting so groups 18-21 have been omitted.

IDRXN = 4: Dual Monod Biodegradation Reaction. Biodegradation is an irreversible process in which bacteria oxidize an organic substrate to produce energy and biomass. In addition to biomass, the biodegradation process requires the presence of an electron acceptor (e.g. oxygen, nitrate, etc.) and nutrients (e.g. nitrogen and phosphorous). An example of a biodegradation reaction is given by the following reaction:

\(\mathrm{Substrate} + \mathrm{Electron Acceptor} + H^+ + \mathrm{Nutrients} \rightarrow \mathrm{cells} + CO_2 + H_2O + NH_3\)

FEHM models the rate of biodegradation of a substrate with a multiplicative Monod model, which is given by:

\(R_s = -q_m m_b \frac{[S]}{K_s + [S]K_A} \frac{[A]}{K_A + [A]}\)

where [S] is the aqueous concentration of substrate (a.k.a the electron donor) {SUBSTRATE}, [A] is the aqueous concentration of the electron acceptor {ELECACC}, and \(m_b\) is the concentration of the immobile biomass {BIOMASS}.

The parameter qm is the maximum specific rate of substrate utilization {QM}, which represents the maximum amount of substrate that can be consumed per unit mass of bacteria per unit time. The parameters \(K_S\) {KS} and \(K_A\) {KA} are the Monod half-maximum-rate concentrations for the electron donor and electron acceptor, respectively. The rate of microbial growth is given by the synthesis rate (which is proportional to the rate of substrate degradation) minus a first-order decay rate.

\(R_{cells} = -YR_s - b(m_b - m_{b,init})\)

where Y is the microbial yield coefficient {YIELD} and b is the first-order microbial decay coefficient {DECAY}. In the above equation, the assumption is made that the background conditions are sufficient to sustain a microbial population of a given size; therefore, the biomass concentration is not allowed to fall below its initial background concentration (\(m_{b,init}\)) {XMINIT}. In the following example of input for a dual monod biodegradation reaction, aqueous component 3 is the substrate and aqueous component 4 is the electron acceptor. Note that there are only two entries in Group 15 and Groups 21-23 are omitted since NAQUEOUS = 2. In addition Group 29 data is left out since NBIOFRM = 0.

IDRXN = 5: Radioactive Decay Reaction. Radioactive decay is a simple first order decay process given by:

\(A \rightarrow B\)

where A is the parent {PARENT} and B is the daughter product {DAUGHTER}. The half life of the reaction is defined as the time it takes for the concentration of A to decrease by a factor of 2. In the following example of input for a radioactive decay reaction, aqueous component 2, the parent, reacts to form aqueous component 3, the daughter product. The half life for the reaction is 432.0 years.

IDRXN = 7: Kinetic Precipitation/Dissolution Reaction. A general reaction describing the precipitation/dissolution of a mineral p {IIMMOBILE} can be written in the following form:

\(m_p \Leftrightarrow \mu_{p1}c_1 + \mu_{p1}c_2 + ... + \mu_{p1}c_{N_c}\)

where the \(c_j\) are the aqueous concentrations {IAQUEOUS} and the \(\mu_{pj}\) are stoichiometric coefficients {AQSTOIC}. The equilibrium constant for this reaction is known as the solubility product. Since the activity of a pure solid is equal to one, the reaction quotient \(Q_p\) is defined as follows:

\(Q_p = \prod_{j=1}^{N_c} c_j^{\mu_{pj}}\)

At equilibrium, \(Q_p\) is equal to the solubility product. The surface-controlled rate of precipitation/dissolution of a mineral is given by:

\(R(m_p) = sign\left(\log\frac{Q_p}{K_{sp}} \right) A_p k_p \left| \frac{Q_p}{K_{sp}} -1 \right|\)

where Ap is reactive surface area of the mineral {AREA}, \(k_p\) is the precipitation rate constant {RATE}, and \(K_{sp}\) is the solubility product {SOLUBILITY}. Currently, this precipitation/dissolution subroutine only allows for the total aqueous concentration of a component to dissolve. The dissolution of uncomplexed aqueous concentration and complex concentrations is not currently supported.

The following is example input for a kinetic precipitation/dissolution reaction (calcite dissolution). This example corresponds to the example used for the trac macro (See In the following example of trac, calcite dissolution is simulated. The input groups are given to the right of the table to facilitate review of the example. The initial solute concentration is set to 0 (but is later overwritten by group 14 input), the implicitness factor is 1 resulting in a 1st order solution, the equation tolerance is 1.e-7, and the upstream weighting is set to 0.5. The solute transport solution is turned on as the heat and mass solution is turned off at day 1. The heat and mass solution resumes on day 1000. Two solutes are simulated in this example. Solute 1 is a nonsorbing (conservative) liquid species (a1 = a2 = 0., b = 1.) with a molecular diffusion coefficient of 1.e-9 m2/s, and dispersivity of 0.0067 m in the X-direction. This transport model applies over the entire problem domain. The initial solute concentration for solute 1 is 6.26e-5 mol/kg-water. Solute 2 is a solid species with an initial solute concentration of 2.e-5 mol/kg-solid. There is no solute source for either solute. The corresponding data for the rxn macro, that would complete this example, is given on page 159.). A single kinetic reaction is specified with no aqueous complexes. The liquid and immobile components are identified. Aqueous component 1 participates in the precipitation/dissolution reaction with immobile component 1. No data is input for groups 4, 6, 9, 10, and 11. The solubility product for the reaction is 6.26e-5 and the reaction rate constant is 100.

IDRXN = 8: Kinetic Precipitation/Dissolution Reaction (rates based on free-ion concentrations). The reaction modeled is analogous to that for IDRXN =7, except that rates are based on the uncomplexed (free-ion) concentration of the species. The total concentration is equal to the free ion concentration + all of the complex concentrations. For a more detailed discussion of the differences between total aqueous and free-ion concentration, see the “Models and Methods Summary” of the FEHM Application [Zyvoloski et al. 1999, [[FEHM-MMS.htm#14829|]], [[FEHM-MMS.htm#14829|]]].

For example, for a species such as Cobalt from the multisolute problem (see [[FEHM-UM.9.5.htm#76956|See Reactive Transport Example]]), the free ion concentration is simply the concentration of Cobalt in it’s uncomplexed state. The total Cobalt would be Free Ion Cobalt + all Cobalt Complexes (e.g. CoEDTA from the multisolute verification problem). Using IDRXN = 7 allows the total Cobalt to dissolve, while IDRXN= 8 allows only the free ion Cobalt to dissolve.