Chemical Reactions Module

The chemical reactions module provides a set of tools for the calculation of multicomponent aqueous reactive transport in porous media, originally developed as the MOOSE application RAT (Guo et al., 2013).

Theory

The first part of defining the chemistry of a problem is to choose a set of independent primary species from which every other species (including minerals) can be expressed in terms of. Other chemical species that can be expressed as combinations of primary species are termed secondary species.

Following Lichtner (1996), the mass conservation equation is formulated in terms of the total concentration of a primary species jj, Ψj\Psi_j, and has the form

t(ϕΨj)+(qϕD)Ψj+m=1NmνjmImQj=0,j=1,2,,Nc,\frac{\partial}{\partial t} \left(\phi \Psi_j \right) + \nabla \cdot \left(\mathbf{q} - \phi D \nabla \right)\Psi_j + \sum_{m=1}^{N_m} \nu_{jm} I_m - Q_j= 0, \quad j = 1, 2, \ldots, N_c, where ϕ\phi is porosity, DD is the hydrodynamic dispersion tensor, NmN_m is the number of minerals formed by kinetic reactions, νjm\nu_{jm} is the stoichiometric coefficient for the jthj^{\mathrm{th}} species in the mthm^{\mathrm{th}} kinetic reaction, ImI_m is the reaction rate for the mthm^{\mathrm{th}} mineral, and QjQ_j is a source (sink) of the jthj^{\mathrm{th}} species, and q\mathbf{q} is the Darcy flux

q=Kμ(Pρg),\mathbf{q} = - \frac{K}{\mu} \left(\nabla P - \rho \mathbf{g}\right), where KK is permeability, μ\mu is fluid viscosity, PP is pressure, ρ\rho is fluid density, and g\mathbf{g} is gravity.

The total concentration of the jthj^{\mathrm{th}} primary species, Ψj\Psi_j, is defined as

Ψj=Cj+i=1NsνjiCi,\Psi_j = C_j + \sum_{i=1}^{N_s} \nu_{ji} C_i, where CjC_j is the concentration of the primary species, CiC_i is the concentration of the ithi^{\mathrm{th}} secondary species, NsN_s is the total number of secondary species, and νji\nu_{ji} are the stoichiometric coefficients.

Aqueous equilibrium reactions

The concentration of the ithi^{\mathrm{th}} secondary species, CiC_i, is calculated from mass action equations corresponding to equilibrium reactions

jνjiAjAi,\sum_j \nu_{ji} \mathcal{A}_j \rightleftharpoons \mathcal{A}_i, where Aj\mathcal{A}_j refers to the jthj^{\mathrm{th}} species. This yields

Ci=Kiγij=1Nc(γjCj)νji,C_i = \frac{K_i}{\gamma_i} \prod_{j=1}^{N_c} \left(\gamma_j C_j\right)^{\nu_{ji}}, where KiK_i is equilibrium constant, γi\gamma_i is the activity coefficient, and NcN_c is the number of primary species.

Solid kinetic reactions

Mineral precipitation/dissolution is possible via kinetic reactions of the form

jνjiAjMm,\sum_j \nu_{ji} \mathcal{A}_j \rightleftharpoons \mathcal{M}_m, where Mm\mathcal{M}_m refers to the mthm^{\mathrm{th}} mineral species.

The reaction rate ImI_m is based on transition state theory, a simple form of which gives

Im=±kmam1Ωmθη,I_m = \pm k_m a_m \left|1 - \Omega_m^{\theta}\right|^{\eta}, where ImI_m is positive for dissolution and negative for precipitation, kmk_m is the rate constant, ama_m is the specific reactive surface area, Ωm\Omega_m is termed the mineral saturation ratio, expressed as

Ωm=1Kmj(γjCj)νjm,\Omega_m = \frac{1}{K_m} \prod_{j}(\gamma_j C_j)^{\nu_{jm}}, where KmK_m is the equilibrium constant for mineral mm.

The rate constant kmk_m is typically reported at a reference temperature (commonly 25^{\circ}C). Using an Arrhenius relation, the temperature dependence of kmk_m is given as

km(T)=km,0exp[EaR(1T01T)],k_m(T) = k_{m,0} \exp\left[\frac{E_a}{R} \left(\frac{1}{T_0} - \frac{1}{T}\right)\right], where km,0k_{m,0} is the rate constant at reference temperature T0T_0, EaE_a is the activation energy, RR is the gas constant.

The exponents θ\theta and η\eta in the reaction rate equation are specific to each mineral reaction, and should be measured experimentally. For simplicity, they are set to unity in this module.

The rate of change in molar concentration, CmC_m, of the mthm^{\mathrm{th}} mineral species is

Cmt=Im.\frac{\partial C_m}{\partial t} = I_m.

This can be expressed in terms of the mineral volume fraction

ϕm=VmCm,\phi_m = \overline{V}_m C_m, where Vm\overline{V}_m is the molar volume of the mthm^{\mathrm{th}} mineral species, ϕmt=VmIm.\frac{\partial \phi_m}{\partial t} = \overline{V}_m I_m.

The total porosity, ϕt\phi_t, can be defined in terms of the volume fraction of all mineral species

ϕt=1mϕm.\phi_t = 1 - \sum_m \phi_m.

In this way, the change in porosity due to mineral precipitation or dissolution can be calculated, which can then be used to calculate changes in other material properties such as permeability.

Implementation details

The physics described above is implemented in a number of Kernels and AuxKernels

Kernels

The transport of each primary species is calculated using the following Kernels:

The transport of primary species present in a secondary species is included using:

The amount of primary species converted to an immobile mineral phase is given by

  • CoupledBEKinetic Rate of change of primary species concentration in mineral species due to dissolution or precipitation

The Darcy flux is calculated using

AuxKernels

The following AuxKernels are used to calculate secondary species and mineral concentrations, as well as total primary species concentration and solution pH.

Material properties

The Kernels above require several material properties to be defined using the following names: porosity, diffusivity and conductivity. These can be defined using one of the Materials available in the framework. For example, constant properties can be implemented using a GenericConstantMaterial with the following:

Required material properties

[Materials]
  [./porous]
    type = GenericConstantMaterial
    prop_names = 'diffusivity conductivity porosity'
    prop_values = '1e-4 1e-4 0.2'
  [../]
[]
(modules/chemical_reactions/test/tests/aqueous_equilibrium/1species.i)

More complicated formulations can be added by creating new Materials as required.

Boundary condition

A flux boundary condition, ChemicalOutFlowBC is provided to define u\nabla u on a boundary.

Postprocessors

The total volume fraction of a given mineral species can be calculated using a TotalMineralVolumeFraction postprocessor.

Reaction network parser

The chemical reactions module includes a reaction network parser in the Actions system that enables chemical reactions to be specified in a natural form in the input file. The parser then adds all required Variables, AuxVariables, Kernels and AuxKernels to represent the total geochemical model. To use the reaction network parser, a ReactionNetwork block can be added to the input file.

Equilibrium reactions can be entered in the ReactionNetwork block using an AqueousEquilibriumReactions sub-block, while kinetic reactions are entered in a SolidKineticReactions sub-block.

The input file syntax for equilibrium reactions has to following form:

ν11A1+ν21A2+=Aeq1log10(Keq),ν12A1+ν22A2+=Aeq2log10(Keq)\begin{aligned} \nu_{11} \mathcal{A}_1 + \nu_{21} \mathcal{A}_2 + \ldots =& \mathcal{A}_{eq1} \log_{10}(K_{eq}),\\ \nu_{12} \mathcal{A}_1 + \nu_{22} \mathcal{A}_2 + \ldots =& \mathcal{A}_{eq2} \log_{10}(K_{eq}) \end{aligned}

Individual equilibrium reactions are provided with the primary species on the left hand side, while the equilibrium species follows the = sign, followed by the log of the equilibrium constant. A comma is used to delimit reactions, so that multiple equilibrium reactions can be entered.

The syntax for solid kinetic reactions is similar, except that no equilibrium constant is entered in the reactions block.

To demonstrate the use of the reaction network parser, consider the geochemical model used in Guo et al. (2013), which features aqueous equilibrium reactions as well as kinetic mineral dissolution and precipitation.

Equilibrium reactions:

H++HCO3CO2(aq)Keq=106.341HCO3H++CO32Keq=1010.325Ca2++HCO3H++CaCO3(aq)Keq=107.009Ca2++HCO3CaHCO3+Keq=100.653Ca2+H++CaOH+Keq=1012.85H+OHKeq=1013.991\begin{aligned} H^+ + HCO_3^- &\rightleftharpoons CO_2(aq) & K_{eq} &= 10^{6.341} \\ HCO_3^- &\rightleftharpoons H^+ + CO_3^{2-} & K_{eq} &= 10^{-10.325} \\ Ca^{2+} + HCO_3^- &\rightleftharpoons H^+ + CaCO_3(aq) & K_{eq} &= 10^{-7.009} \\ Ca^{2+} + HCO_3^- &\rightleftharpoons CaHCO_3^+ & K_{eq} &= 10^{-0.653} \\ Ca^{2+} &\rightleftharpoons H^+ + CaOH^+ & K_{eq} &= 10^{-12.85} \\ - H^+ &\rightleftharpoons OH^- & K_{eq} &= 10^{-13.991} \end{aligned}

Kinetic reaction:

Ca2++HCO3H++CaCO3(s)Ca^{2+} + HCO_3^- \rightleftharpoons H^+ + CaCO_3(s)

with equilibrium constant Keq=101.8487K_{eq} = 10^{1.8487}, specific reactive surface area am=0.461a_m = 0.461 m2^2/L, kinetic rate constant km=6.456542×102k_m = 6.456542 \times 10^{-2} mol/m2^2 and activation energy Ea=15,000E_a = 15,000 J/mol.

Example of AqueousEquilibriumReactions action.

[ReactionNetwork]
  [./AqueousEquilibriumReactions]
    primary_species = 'ca2+ hco3- h+'
    secondary_species = 'co2_aq co32- caco3_aq cahco3+ caoh+ oh-'
    pressure = pressure
    reactions = 'h+ + hco3- = co2_aq 6.341,
                 hco3- - h+ = co32- -10.325,
                 ca2+ + hco3- - h+ = caco3_aq -7.009,
                 ca2+ + hco3- = cahco3+ -0.653,
                 ca2+ - h+ = caoh+ -12.85,
                 - h+ = oh- -13.991'
  [../]
  [./SolidKineticReactions]
    primary_species = 'ca2+ hco3- h+'
    kin_reactions = 'ca2+ + hco3- - h+ = caco3_s'
    secondary_species = caco3_s
    log10_keq = 1.8487
    reference_temperature = 298.15
    system_temperature = 298.15
    gas_constant = 8.314
    specific_reactive_surface_area = 4.61e-4
    kinetic_rate_constant = 6.456542e-7
    activation_energy = 1.5e4
  [../]
[]
(modules/chemical_reactions/examples/calcium_bicarbonate/calcium_bicarbonate.i)

The reactive transport system above can be provided in the input file without using the reaction network parser. However, this adds more than 400 lines of input in this case, due to the large number of kernels that have to be provided!

Objects, Actions, and, Syntax

AuxKernelsAuxVariablesBCsChemicalComposition
  • Chemical Reactions App
  • CommonChemicalCompositionActionStore common ChemicalComposition action parameters
  • ChemicalCompositionActionTo use this object, you need to have the Thermochimica library installed. Refer to the documentation for guidance on how to enable it. (Original description: Sets up the thermodynamic model and variables for the thermochemistry solve using Thermochimica.)
KernelsMaterials
  • Chemical Reactions App
  • LangmuirMaterialMaterial type that holds info regarding Langmuir desorption from matrix to porespace and viceversa
  • MollifiedLangmuirMaterialMaterial type that holds info regarding MollifiedLangmuir desorption from matrix to porespace and viceversa
PostprocessorsReactionNetworkReactionNetwork/AqueousEquilibriumReactionsReactionNetwork/SolidKineticReactionsUserObjects
  • Chemical Reactions App
  • ThermochimicaElementDataTo use this object, you need to have the Thermochimica library installed. Refer to the documentation for guidance on how to enable it. (Original description: Provides access to Thermochimica-calculated data at elements.)
  • ThermochimicaNodalDataTo use this object, you need to have the Thermochimica library installed. Refer to the documentation for guidance on how to enable it. (Original description: Provides access to Thermochimica-calculated data at nodes.)

References

  1. Luanjing Guo, Hai Huang, Derek R Gaston, Cody J Permann, David Andrs, George D Redden, Chuan Lu, Don T Fox, and Yoshiko Fujita. A parallel, fully coupled, fully implicit solution to reactive transport in porous media using the preconditioned Jacobian-Free Newton-Krylov Method. Advances in Water Resources, 53:101–108, 2013.[BibTeX]
  2. Peter C Lichtner. Continuum formulation of multicomponent-multiphase reactive transport. Reviews in Mineralogy and Geochemistry, 34:1–81, 1996.[BibTeX]