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 , , and has the form
where is porosity, is the hydrodynamic dispersion tensor, is the number of minerals formed by kinetic reactions, is the stoichiometric coefficient for the species in the kinetic reaction, is the reaction rate for the mineral, and is a source (sink) of the species, and is the Darcy flux
where is permeability, is fluid viscosity, is pressure, is fluid density, and is gravity.
The total concentration of the primary species, , is defined as
where is the concentration of the primary species, is the concentration of the secondary species, is the total number of secondary species, and are the stoichiometric coefficients.
Aqueous equilibrium reactions
The concentration of the secondary species, , is calculated from mass action equations corresponding to equilibrium reactions
where refers to the species. This yields
where is equilibrium constant, is the activity coefficient, and is the number of primary species.
Solid kinetic reactions
Mineral precipitation/dissolution is possible via kinetic reactions of the form
where refers to the mineral species.
The reaction rate is based on transition state theory, a simple form of which gives
where is positive for dissolution and negative for precipitation, is the rate constant, is the specific reactive surface area, is termed the mineral saturation ratio, expressed as
where is the equilibrium constant for mineral .
The rate constant is typically reported at a reference temperature (commonly 25C). Using an Arrhenius relation, the temperature dependence of is given as
where is the rate constant at reference temperature , is the activation energy, is the gas constant.
The exponents and 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, , of the mineral species is
This can be expressed in terms of the mineral volume fraction
where is the molar volume of the mineral species,
The total porosity, , can be defined in terms of the volume fraction of all mineral species
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
:
PrimaryTimeDerivative
Rate of change of primary species concentrationPrimaryConvection
Convective flux of primary speciesPrimaryDiffusion
Diffusion of primary species
The transport of primary species present in a secondary species is included using:
CoupledBEEquilibriumSub
Rate of change of primary species concentration in an equilibrium secondary speciesCoupledConvectionReactionSub
Convection of primary species concentration in an equilibrium secondary speciesCoupledDiffusionReactionSub
Diffusion of primary species concentration in an equilibrium secondary species
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
DarcyFluxPressure
Darcy flux
AuxKernels
The following AuxKernels
are used to calculate secondary species and mineral concentrations, as well as total primary species concentration and solution pH.
AqueousEquilibriumRxnAux
The concentration of secondary species for the equilibrium reactionKineticDisPreRateAux
The kinetic rate of the kinetic reactionKineticDisPreConcAux
The concentration of mineral speciesTotalConcentrationAux
The total concentration of a given primary speciesPHAux
The pH of the solutionEquilibriumConstantAux
Temperature-dependent equilibrium constant
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:
[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 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:
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:
Kinetic reaction:
with equilibrium constant , specific reactive surface area m/L, kinetic rate constant mol/m and activation energy 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
[../]
[]
(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
- 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
- ChemicalCompositionActionSets 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
- ThermochimicaElementDataProvides access to Thermochimica-calculated data at elements.
- ThermochimicaNodalDataProvides access to Thermochimica-calculated data at nodes.
References
- 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]
- Peter C Lichtner.
Continuum formulation of multicomponent-multiphase reactive transport.
Reviews in Mineralogy and Geochemistry, 34:1–81, 1996.[BibTeX]