Differences between MOOSE and the Geochemist's Workbench

The Geochemists Workbench (GWB) software is a popular "gold-standard" geochemistry solver, and to help new users of the geochemistry module, many tests and examples are provided with GWB equivalents.

On the other hand, the C++ code behind the geochemistry solver is different than GWB and Bethke (2007) because it is part of the MOOSE framework, and this leads to a number of differences that advanced users should be aware of. In almost all cases, these differences make little difference to the final results (certainly within the accuracy of the database and experimental observations) and are only of importance if exact equivalence is required (such as for benchmarking).

Temperature dependence

Special temperature

Inspecting the tests and examples will reveal that the flag piecewise_linear_interpolation is frequently set true in the GeochemicalModelDefinition. This is because that at the special temperatures 0, 25, 60, 100 150, 200, 250 and 300C, GWB does not perform a fourth-order least-squares fit to the equilibrium-constant data. Instead, it uses the data exactly as tabulated in the database file.

High temperatures

In all benchmarks run, the geochemistry module results agree exactly with GWB when the temperature is less than roughly 150C. However, the results sometimes differ above this value. The small discrepancies are probably due to slightly different interpolations of equilibrium constants, although it is difficult to be completely sure without being able to inspect the GWB code.

Stoichiometric ionic strength

The GWB software calculates the stoichiometric ionic strength (for computing the water activity via the Debye-Huckel model) using the Cl molality only. Hence, most of the tests and examples use the flag stoichiometric_ionic_str_using_Cl_only = true.

Setting bulk constraints

When the user sets moles_bulk_species, this means the mole number in the current basis (after the swaps). For instance, the following means there is:

  • 1 mole of Na. There is also additional Na inside the other basis species Mirabilite. Some of the 1mol will be free, while the remaining part will be bound into secondary species such as NaCl.

  • 1 mole of Cl. Some of this will be free, while the other part will be bound into secondary species such as NaCl.

  • 1 free mole of Mirabilite precipitate floating around in the solution, which contains 2 moles of Na, 1 mole of SO and 10 moles of HO

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O Cl- Na+ SO4--"
    equilibrium_minerals = "Mirabilite"
    piecewise_linear_interpolation = true # for comparison with GWB
  []
[]

[TimeIndependentReactionSolver]
  model_definition = definition
  swap_out_of_basis = "SO4--"
  swap_into_basis = "  Mirabilite"
  charge_balance_species = "Cl-"
  constraint_species = "H2O    Na+ Cl- Mirabilite"
  constraint_value = "  1.0    1.0 1.0 1.0"
  constraint_meaning = "kg_solvent_water bulk_composition bulk_composition free_mineral"
  constraint_unit = "kg moles moles moles"
  ramp_max_ionic_strength_initial = 0 # not needed in this simple problem
  stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
  add_aux_pH = false
  mol_cutoff = 1E-5
  abs_tol = 1E-15
[]
(modules/geochemistry/test/tests/equilibrium_models/bulk_constraints.i)

The result is shown in Table 1. Note that Mirabilite has a bulk mole number of 1.895. One mol of this Mirabilite is the free quantity, and 0.895mol is dissolved into solution and has added to the molality of the other species (Na, NaSO, SO and NaCl).

Table 1: Species molality for the bulk_constraints simulation

SpeciesMolality
Na2.266
Cl0.976
NaSO0.4991
SO0.3957
NaCl0.02402
Mirabilite1.895 bulk mol

A GWB input file that looks equivalent is

# React script that appears equivalent to bulk_constraints.i but is not
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O          = 1 free kg
Na+          = 1 mol
Cl-          = 1 mol
balance on Cl-
swap Mirabilite for SO4--
Mirabilite   = 1 free mol
suppress all
unsuppress Mirabilite
printout  species = long
epsilon = 1e-13
go
(modules/geochemistry/test/tests/equilibrium_models/bulk_constraints.rea)

In fact, this does not converge in GWB (unless the charge-balance species is moved to Na). The file that is actually equivalent sets


Na+ = 2.79 mol

The additional 1.79mol is twice the 0.895mol of dissolved Mirabilite (each mole of Mirabilite contains two moles of Na).

The difference is that GWB expects the bulk composition defined in the bases before the swaps, while geochemistry expects it after the swaps.

In almost all instances this is not an issue because there are no swaps, or because the differences are small (the above example is an extreme case). For instance, none of the results of the tests and examples are greatly impacted by these differences. If the differences prove to be annoying to users then the GeochemicalSolver could be modified with a flag that allows the GWB-style system to be solved, but this change is more than a mere cosmetic change.

Setting activity

When the user sets a constraint on a species' activity then MOOSE respects that constraint until it is turned off. For instance, in the following input file the pH is set to 6 and remains so throughout the simulation, even as the Calcite is precipitating.

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O Cl- Ca++ H+ HCO3-"
    equilibrium_minerals = "Calcite"
    piecewise_linear_interpolation = true # for comparison with GWB
  []
[]

[TimeIndependentReactionSolver]
  model_definition = definition
  charge_balance_species = "Cl-"
  constraint_species = "H2O              H+       Ca++               Cl-                HCO3-"
  constraint_value = "  1.0              1E-6     1                  1                  1"
  constraint_meaning = "kg_solvent_water activity bulk_composition bulk_composition bulk_composition"
  constraint_unit = "   kg             dimensionless moles           moles              moles"
  ramp_max_ionic_strength_initial = 0 # not needed in this simple problem
  stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
  mol_cutoff = 1E-5
  abs_tol = 1E-13
[]
(modules/geochemistry/test/tests/equilibrium_models/ph_constraint.i)

The final result is 97.92g of free calcite in the solution at pH 6.

This is different from the seemingly-equivalent GWB input file:

# React script that appears equivalent to ph_constraint.i but is not
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O          = 1 free kg
H+           = 6 pH
Ca++         = 1 mol
Cl-          = 1 mol
balance on Cl-
HCO3-        = 1 mol
suppress all
unsuppress Calcite
printout  species = long
epsilon = 1e-13
go
(modules/geochemistry/test/tests/equilibrium_models/ph_constraint.rea)

Note that there is no fixing of the pH, so the Geochemists Workbench software:

  1. Equilibrates the solution without precipitating minerals at pH 6

  2. Closes the system so that the bulk mole composition of H is fixed

  3. Allows the Calcite precipitates to form, which alters the pH

The final result is 25.98g of free calcite precipitate at pH 4.689.

If a user wants to replicate the GWB procedure, two models must be run in succession: the first has no minerals and the user records the bulk composition of each species; the second has minerals and the bulk composition fixed. The geochemistry module can be enhanced if this proves to be an annoying feature.

  • A close_before_precipitating flag can be added to the TimeIndependentReactionSolver.

  • A prevent_precipitation_on input can be added to the TimeDependentReactionSolver. It is an AuxVariable: when it is 1.0 (the default) then precipitation of minerals specified in prevent_precipitation is prevented, otherwise precipitates can form (ie, the prevent_precipitation list is ignored).

  • Both of these involve simple manipulations of the _prevent_precipitation data of the GeochemicalSolver.

References

  1. Craig M. Bethke. Geochemical and Biogeochemical Reaction Modeling. Cambridge University Press, 2 edition, 2007. doi:10.1017/CBO9780511619670.[BibTeX]