An acidic solution
A particularly simple geochemistry
example involves computing the molality, etc, of species in an acidic solution with fixed pH.
The MOOSE input file: step 1
The database and basis species must be specified for all geochemistry
simulations. In this case, the full moose_geochemdb.json
database is used, and the basis species are simply HO, H and Cl. A piecewise linear interpolation for the equilibrium constants and the Debye-Huckel activity coefficients is used so that the result agrees exactly with the analogous GWB simulation:
[UserObjects]
[definition]
type = GeochemicalModelDefinition
database_file = "../../../database/moose_geochemdb.json"
basis_species = "H2O H+ Cl-"
piecewise_linear_interpolation = true # to reproduce the GWB result
[]
[]
(modules/geochemistry/test/tests/equilibrium_models/HCl.i)The MOOSE input file: step 2
The next piece of the input file involves specifying the initial conditions and the simulation type. This is a time-independent solve (just the equilibrium configuration is sought). The charge-balance species is chosen to be Cl, there is 1kg of solvent water and the pH is fixed to 2 (via log10activity = -2
):
[TimeIndependentReactionSolver]
model_definition = definition
geochemistry_reactor_name = reactor
charge_balance_species = "Cl-"
constraint_species = "H2O H+ Cl-"
constraint_value = " 1.0 -2 1E-2"
constraint_meaning = "kg_solvent_water log10activity bulk_composition"
constraint_unit = " kg dimensionless moles"
ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
execute_console_output_on = initial
abs_tol = 1E-15
[]
(modules/geochemistry/test/tests/equilibrium_models/HCl.i)The MOOSE input file: optional step 3
Postprocessors may be used to produce some CSV output. In this example, all the information is already obtained in the console output, but by way of example, the following postprocessors output the activity of water, the pH, the solvent-water mass, the Cl molality, etc:
[Postprocessors]
[pH]
type = PointValue
variable = 'pH'
[]
[solvent_mass]
type = PointValue
variable = 'kg_solvent_H2O'
[]
[molal_Cl-]
type = PointValue
variable = 'molal_Cl-'
[]
[mg_per_kg_HCl]
type = PointValue
variable = 'mg_per_kg_HCl'
[]
[activity_OH-]
type = PointValue
variable = 'activity_OH-'
[]
[bulk_H+]
type = PointValue
variable = 'bulk_moles_H+'
[]
[temperature]
type = PointValue
variable = 'solution_temperature'
[]
[]
(modules/geochemistry/test/tests/equilibrium_models/HCl.i)In addition:
[GlobalParams]
point = '0 0 0'
[]
(modules/geochemistry/test/tests/equilibrium_models/HCl.i)[Outputs]
csv = true
[]
(modules/geochemistry/test/tests/equilibrium_models/HCl.i)The equivalent GWB input file
The results may be compared with those generated by Geochemists Workbench. The GWB input file is
# React script that is analogous to HCl.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O = 1 free kg
Cl- = .01 molal
balance on Cl-
pH = 2
printout species = long
epsilon = 1e-15
(modules/geochemistry/test/tests/equilibrium_models/HCl.rea)Results
The geochemistry
module produces the same result as GWB (up to precision):
Summary:
Total number of iterations required = 9
Error in calculation = 6.888e-17mol
Charge of solution = 0mol (charge-balance species = Cl-)
Mass of solvent water = 1kg
Mass of aqueous solution = 1kg
pH = 2
Ionic strength = 0.01097mol/kg(solvent water)
Stoichiometric ionic strength = 0.01097mol/kg(solvent water)
Activity of water = 0.9996
Temperature = 25
Basis Species:
H+; bulk_moles = 0.01097mol; bulk_conc = 11.05mg/kg(solution); molality = 0.01097mol/kg(solvent water); free_conc = 11.06mg/kg(solvent water); act_coeff = 0.9114; log10(a) = -2
Cl-; bulk_moles = 0.01097mol; bulk_conc = 388.8mg/kg(solution); molality = 0.01097mol/kg(solvent water); free_conc = 389mg/kg(solvent water); act_coeff = 0.8956; log10(a) = -2.008
Equilibrium Species:
HCl; molality = 7.805e-11mol/kg(solvent water); free_conc = 2.846e-06mg/kg(solvent water); act_coeff = 1; log10(a) = -10.11; HCl = 1*H+ + 1*Cl-; log10K = 6.1
OH-; molality = 1.149e-12mol/kg(solvent water); free_conc = 1.954e-08mg/kg(solvent water); act_coeff = 0.8971; log10(a) = -11.99; OH- = 1*H2O - 1*H+; log10K = 13.99
Without an Action
Although most users will want to use the Action system, it is possible to build all geochemistry
models without using Actions for slightly more fine-grained control. The equivalent input file is:
# This is an example of an input file that does not utilize an action. Its functionality is the same as HCl.i
# This solves for molalities in a system just containing HCl
[GlobalParams]
point = '0 0 0'
[]
[Mesh]
type = GeneratedMesh
dim = 1
nx= 1
[]
[Variables]
[u]
[]
[]
[Kernels]
[u]
type = Diffusion
variable = u
[]
[]
[AuxVariables]
[solution_temperature]
[]
[kg_solvent_H2O]
[]
[activity_H2O]
[]
[bulk_moles_H2O]
[]
[pH]
[]
[molal_H+]
[]
[molal_Cl-]
[]
[molal_HCl]
[]
[molal_OH-]
[]
[mg_per_kg_H+]
[]
[mg_per_kg_Cl-]
[]
[mg_per_kg_HCl]
[]
[mg_per_kg_OH-]
[]
[activity_H+]
[]
[activity_Cl-]
[]
[activity_HCl]
[]
[activity_OH-]
[]
[bulk_moles_H+]
[]
[bulk_moles_Cl-]
[]
[bulk_moles_HCl]
[]
[bulk_moles_OH-]
[]
[]
[AuxKernels]
[solution_temperature]
type = GeochemistryQuantityAux
species = 'H+'
reactor = reactor
variable = solution_temperature
quantity = temperature
[]
[kg_solvent_H2O]
type = GeochemistryQuantityAux
species = 'H2O'
reactor = reactor
variable = kg_solvent_H2O
quantity = molal
[]
[activity_H2O]
type = GeochemistryQuantityAux
species = 'H2O'
reactor = reactor
variable = activity_H2O
quantity = activity
[]
[bulk_moles_H2O]
type = GeochemistryQuantityAux
species = 'H2O'
reactor = reactor
variable = bulk_moles_H2O
quantity = bulk_moles
[]
[pH]
type = GeochemistryQuantityAux
species = 'H+'
reactor = reactor
variable = pH
quantity = neglog10a
[]
[molal_H+]
type = GeochemistryQuantityAux
species = 'H+'
reactor = reactor
variable = 'molal_H+'
quantity = molal
[]
[molal_Cl-]
type = GeochemistryQuantityAux
species = 'Cl-'
reactor = reactor
variable = 'molal_Cl-'
quantity = molal
[]
[molal_HCl]
type = GeochemistryQuantityAux
species = 'HCl'
reactor = reactor
variable = 'molal_HCl'
quantity = molal
[]
[molal_OH-]
type = GeochemistryQuantityAux
species = 'OH-'
reactor = reactor
variable = 'molal_OH-'
quantity = molal
[]
[mg_per_kg_H+]
type = GeochemistryQuantityAux
species = 'H+'
reactor = reactor
variable = 'mg_per_kg_H+'
quantity = mg_per_kg
[]
[mg_per_kg_Cl-]
type = GeochemistryQuantityAux
species = 'Cl-'
reactor = reactor
variable = 'mg_per_kg_Cl-'
quantity = mg_per_kg
[]
[mg_per_kg_HCl]
type = GeochemistryQuantityAux
species = 'HCl'
reactor = reactor
variable = 'mg_per_kg_HCl'
quantity = mg_per_kg
[]
[mg_per_kg_OH-]
type = GeochemistryQuantityAux
species = 'OH-'
reactor = reactor
variable = 'mg_per_kg_OH-'
quantity = mg_per_kg
[]
[activity_H+]
type = GeochemistryQuantityAux
species = 'H+'
reactor = reactor
variable = 'activity_H+'
quantity = activity
[]
[activity_Cl-]
type = GeochemistryQuantityAux
species = 'Cl-'
reactor = reactor
variable = 'activity_Cl-'
quantity = activity
[]
[activity_HCl]
type = GeochemistryQuantityAux
species = 'HCl'
reactor = reactor
variable = 'activity_HCl'
quantity = activity
[]
[activity_OH-]
type = GeochemistryQuantityAux
species = 'OH-'
reactor = reactor
variable = 'activity_OH-'
quantity = activity
[]
[bulk_moles_H+]
type = GeochemistryQuantityAux
species = 'H+'
reactor = reactor
variable = 'bulk_moles_H+'
quantity = bulk_moles
[]
[bulk_moles_Cl-]
type = GeochemistryQuantityAux
species = 'Cl-'
reactor = reactor
variable = 'bulk_moles_Cl-'
quantity = bulk_moles
[]
[bulk_moles_HCl]
type = GeochemistryQuantityAux
species = 'HCl'
reactor = reactor
variable = 'bulk_moles_HCl'
quantity = bulk_moles
[]
[bulk_moles_OH-]
type = GeochemistryQuantityAux
species = 'OH-'
reactor = reactor
variable = 'bulk_moles_OH-'
quantity = bulk_moles
[]
[]
[Postprocessors]
[pH]
type = PointValue
variable = 'pH'
[]
[solvent_mass]
type = PointValue
variable = 'kg_solvent_H2O'
[]
[molal_Cl-]
type = PointValue
variable = 'molal_Cl-'
[]
[mg_per_kg_HCl]
type = PointValue
variable = 'mg_per_kg_HCl'
[]
[activity_OH-]
type = PointValue
variable = 'activity_OH-'
[]
[bulk_H+]
type = PointValue
variable = 'bulk_moles_H+'
[]
[temperature]
type = PointValue
variable = 'solution_temperature'
[]
[]
[Executioner]
type = Steady
[]
[UserObjects]
[definition]
type = GeochemicalModelDefinition
database_file = "../../../database/moose_geochemdb.json"
basis_species = "H2O H+ Cl-"
piecewise_linear_interpolation = true # to reproduce the GWB result
[]
[reactor]
type = GeochemistryTimeDependentReactor
model_definition = definition
charge_balance_species = "Cl-"
constraint_species = "H2O H+ Cl-"
constraint_value = " 1.0 -2 1E-2"
constraint_meaning = "kg_solvent_water log10activity bulk_composition"
constraint_unit = " kg dimensionless moles"
ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
abs_tol = 1E-15
[]
[nnn]
type = NearestNodeNumberUO
[]
[]
[Outputs]
csv = true
[console_output]
type = GeochemistryConsoleOutput
geochemistry_reactor = reactor
nearest_node_number_UO = nnn
solver_info = true
execute_on = initial
[]
[]
(modules/geochemistry/test/tests/equilibrium_models/HCl_no_action.i)