Chemical Reactions Module

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 $j$, $\Psi_j$, and has the form

$$\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, $D$ is the hydrodynamic dispersion tensor, $N_m$ is the number of minerals formed by kinetic reactions, $\nu_{jm}$ is the stoichiometric coefficient for the $j^{\mathrm{th}}$ species in the $m^{\mathrm{th}}$ kinetic reaction, $I_m$ is the reaction rate for the $m^{\mathrm{th}}$ mineral, and $Q_j$ is a source (sink) of the $j^{\mathrm{th}}$ species, and $\mathbf{q}$ is the Darcy flux

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

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

$$\Psi_j = C_j + \sum_{i=1}^{N_s} \nu_{ji} C_i,$$ where $C_j$ is the concentration of the primary species, $C_i$ is the concentration of the $i^{\mathrm{th}}$ secondary species, $N_s$ is the total number of secondary species, and $\nu_{ji}$ are the stoichiometric coefficients.

Implementation details

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

[Materials]
  [./porous]
    type = GenericConstantMaterial
    prop_names = 'diffusivity conductivity porosity'
    prop_values = '1e-4 1e-4 0.2'
  [../]
[]
Copy

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:

$$\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:

$$\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:

$$Ca^{2+} + HCO_3^- \rightleftharpoons H^+ + CaCO_3(s)$$

with equilibrium constant $K_{eq} = 10^{1.8487}$, specific reactive surface area $a_m = 0.461$ m$^2$/L, kinetic rate constant $k_m = 6.456542 \times 10^{-2}$ mol/m$^2$ and activation energy $E_a = 15,000$ J/mol.

[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
  [../]
[]
Copy

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

• Chemical Reactions App
• AqueousEquilibriumRxnAuxConcentration of secondary equilibrium species
• EquilibriumConstantAuxEquilibrium constant for a given equilibrium species (in form log10(Keq))
• KineticDisPreConcAuxConcentration of secondary kinetic species
• KineticDisPreRateAuxKinetic rate of secondary kinetic species
• PHAuxpH of solution
• TotalConcentrationAuxTotal concentration of primary species (including stoichiometric contribution to secondary equilibrium species)

• Chemical Reactions App
• ChemicalOutFlowBCChemical flux boundary condition

• 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.)

• Chemical Reactions App
• CoupledBEEquilibriumSubDerivative of equilibrium species concentration wrt time
• CoupledBEKineticDerivative of kinetic species concentration wrt time
• CoupledConvectionReactionSubConvection of equilibrium species
• CoupledDiffusionReactionSubDiffusion of equilibrium species
• DarcyFluxPressureDarcy flux: - cond * (Grad P - rho * g) where cond is the hydraulic conductivity, P is fluid pressure, rho is fluid density and g is gravity
• DesorptionFromMatrixMass flow rate from the matrix to the porespace. Add this to TimeDerivative kernel to get complete DE for the fluid adsorbed in the matrix
• DesorptionToPorespaceMass flow rate to the porespace from the matrix. Add this to the other kernels for the porepressure variable to form the complete DE
• PrimaryConvectionConvection of primary species
• PrimaryDiffusionDiffusion of primary species
• PrimaryTimeDerivativeDerivative of primary species concentration wrt time

• 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

• Chemical Reactions App
• TotalMineralVolumeFractionTotal volume fraction of coupled mineral species

• Chemical Reactions App
• AqueousEquilibriumReactions
• SolidKineticReactions

• Chemical Reactions App
• AddCoupledEqSpeciesActionAdds coupled equilibrium Kernels and AuxKernels for primary species
• AddPrimarySpeciesActionAdds Variables for all primary species
• AddSecondarySpeciesActionAdds AuxVariables for all secondary species

• Chemical Reactions App
• AddCoupledSolidKinSpeciesActionAdds solid kinetic Kernels and AuxKernels for primary species
• AddPrimarySpeciesActionAdds Variables for all primary species
• AddSecondarySpeciesActionAdds AuxVariables for all secondary species

• 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