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)