Chemical model of Amazon river water

This example closely follows Section 6.2 of Bethke (2007).

A chemical analysis of the major element composition of water in the Amazon river Table 1. The "free" O $var element = document.getElementById("moose-equation-43505ccd-084c-48d9-9486-db4bbf54b13d");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ (aq) concentration of 5.8 $var element = document.getElementById("moose-equation-43d3a6e5-19c3-4770-8767-90c41419403a");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mg.kg $var element = document.getElementById("moose-equation-0152e449-3ac0-4e5a-bf72-b6926e29c2be");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ means its bulk composition will be more than this. In addition, the pH is 6.5. Assume that the extent of the system is 1 $var element = document.getElementById("moose-equation-884088d4-59ba-4724-9d46-e0060ade19e9");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ kg of solvent water along with the solutes dissolved in it.

Table 1: Major element composition of Amazon river water

Species	Concentration (mg.kg $var element = document.getElementById("moose-equation-d018494d-abd3-4f27-944e-52a44aa3c91f");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ )
HCO $var element = document.getElementById("moose-equation-f4733670-0f4b-48ef-a320-f45293ac22cb");katex.render("_{3}^{-}", element, {displayMode:false,throwOnError:false});$	19
SiO $var element = document.getElementById("moose-equation-39eacfbc-b13a-4f93-a387-0cdcdc419674");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ (aq)	7
O $var element = document.getElementById("moose-equation-c2739615-371b-4c20-9513-4358dc7f95dd");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ (aq)	5.8 free
Ca $var element = document.getElementById("moose-equation-26d192be-737a-4e7a-86a5-a7aec6e7be57");katex.render("^{2+}", element, {displayMode:false,throwOnError:false});$	4.3
SO $var element = document.getElementById("moose-equation-8c5e7819-428b-4a65-a1b1-f6ced1dca880");katex.render("_{4}^{2-}", element, {displayMode:false,throwOnError:false});$	3
Cl $var element = document.getElementById("moose-equation-7440edfe-2df7-4303-bf42-8231d56ccd60");katex.render("^{-}", element, {displayMode:false,throwOnError:false});$	4.9
Na $var element = document.getElementById("moose-equation-e0755765-88e9-464c-9d71-175cfac90afa");katex.render("^{+}", element, {displayMode:false,throwOnError:false});$	1.8
Mg $var element = document.getElementById("moose-equation-37f396c4-cc36-453f-8060-0fe04b6f249e");katex.render("^{2+}", element, {displayMode:false,throwOnError:false});$	1.1
Al $var element = document.getElementById("moose-equation-80e15c5e-3cd7-4915-835c-fe1e258ccd63");katex.render("^{3+}", element, {displayMode:false,throwOnError:false});$	0.07
Fe $var element = document.getElementById("moose-equation-0f2ee1a4-956d-4060-99ce-299aa4d5dcf7");katex.render("^{2+}", element, {displayMode:false,throwOnError:false});$	0.06

The concentration of Cl $var element = document.getElementById("moose-equation-898b7e56-8b58-4e30-99b4-4765caf1e56c");katex.render("^{-}", element, {displayMode:false,throwOnError:false});$ is listed as 1.9 $var element = document.getElementById("moose-equation-d03b2399-9236-4e09-9d28-db3b836376f6");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mg.kg $var element = document.getElementById("moose-equation-82f6efff-c2a5-4689-bd68-45629f819416");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ in Bethke (2007), but its molality is computed to be 1.38E-4 $var element = document.getElementById("moose-equation-74a22fc2-18c3-4577-bc78-913ae212e909");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mol/kg(solvent water). Hence its bulk composition must have been at least 1.38E-4 $var element = document.getElementById("moose-equation-a9d0e0a4-bdf1-42c9-81b2-818a7b61f371");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mol since all the stoichiometric coefficients of Cl- in all the secondary species are positive. This implies the bulk composition of Cl $var element = document.getElementById("moose-equation-4d62a712-929f-4ba2-90fe-471ae1cba166");katex.render("^{-}", element, {displayMode:false,throwOnError:false});$ must be more than 1.9 $var element = document.getElementById("moose-equation-5b7baab3-331c-47c1-8238-fbf19225553c");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mg/kg(soln). There is probably a typo in the book. Hence the bulk concentration is set to 4.9 $var element = document.getElementById("moose-equation-35d89f78-89f4-40f3-9f30-c7770c66ea9c");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mg/kg(soln), leading to 1.383E-4 $var element = document.getElementById("moose-equation-70f44d0f-a16e-4bb7-9db0-84ece365bed6");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ moles. In any case, charge balance is performed on Cl $var element = document.getElementById("moose-equation-9bc85ef8-8795-41ff-93c5-2a23466235f9");katex.render("^{-}", element, {displayMode:false,throwOnError:false});$ .

MOOSE input file: no minerals

The MOOSE input file contains the usual GeochemicalModelDefinition that specifies the database file to use, and in this case just the basis species and equilibrium minerals (that are prevented from precipitating, but for which the saturation indices are of interest). The flag piecewise_linear_interpolation = true in order to compare with the Geochemists Workbench result.

[UserObjects<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
  [definition]
    type = GeochemicalModelDefinition<<<{"description": "User object that parses a geochemical database file, and only retains information relevant to the current geochemical model", "href": "../../../source/userobjects/GeochemicalModelDefinition.html"}>>>
    database_file<<<{"description": "The name of the geochemical database file"}>>> = "../../../database/moose_geochemdb.json"
    basis_species<<<{"description": "A list of basis components relevant to the aqueous-equilibrium problem. H2O must appear first in this list.  These components must be chosen from the 'basis species' in the database, the sorbing sites (if any) and the decoupled redox states that are in disequilibrium (if any)."}>>> = "H2O H+ SiO2(aq) Al+++ Fe++ Ca++ Mg++ Na+ HCO3- SO4-- Cl- O2(aq)"
    equilibrium_minerals<<<{"description": "A list of minerals that are in equilibrium with the aqueous solution.  All members of this list must be in the 'minerals' section of the database file"}>>> = "Nontronit-Ca Nontronit-Mg Nontronit-Na Hematite Kaolinite Beidellit-Ca Beidellit-H Beidellit-Mg Beidellit-Na Pyrophyllite Gibbsite Paragonite Quartz"
    piecewise_linear_interpolation<<<{"description": "If true then use a piecewise-linear interpolation of logK and Debye-Huckel parameters, regardless of the interpolation type specified in the database file.  This can be useful for comparing with results using other geochemistry software"}>>> = true # for comparison with GWB
  []
[]

To instruct MOOSE to find the equilibrium configuration, a TimeIndependentReactionSolver is used:

The pH is fixed by setting the activity of H $var element = document.getElementById("moose-equation-5e5d58e3-2241-4e0b-aa87-a3b77335f796");katex.render("^{+}", element, {displayMode:false,throwOnError:false});$ .
The bulk mole number of the aqueous species is also fixed appropriately. The numbers are different than the concentration in mg.kg $var element = document.getElementById("moose-equation-69228fa7-dc25-4337-a31f-40889748a2c4");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ given in the above table, and may be worked out using the TDS.
The prevent_precipitation input prevents any minerals from precipitating when finding the equilibrium configuration, even if their saturation indices are positive.
The other flags enable an accurate comparison with the Geochemists Workbench software.

[TimeIndependentReactionSolver<<<{"href": "../../../syntax/TimeIndependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  charge_balance_species<<<{"description": "Charge balance will be enforced on this basis species.  This means that its bulk mole number may be changed from the initial value you provide in order to ensure charge neutrality.  After the initial swaps have been performed, this must be in the basis, and it must be provided with a bulk_composition constraint_meaning."}>>> = "Cl-"
  constraint_species<<<{"description": "Names of the species that have their values fixed to constraint_value with meaning constraint_meaning.  All basis species (after swap_into_basis and swap_out_of_basis) must be provided with exactly one constraint.  These constraints are used to compute the configuration during the initial problem setup, and in time-dependent simulations they may be modified as time progresses."}>>> = "H2O              H+            O2(aq)             SiO2(aq)         Al+++            Fe++             Ca++             Mg++             Na+              HCO3-            SO4--            Cl-      "
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -6.5          1.813E-4           0.0001165        2.5945E-6        1.0744E-6        0.0001073        4.526E-5         7.83E-5          0.0003114        3.1233E-5        1.383E-4"
  constraint_meaning<<<{"description": "Meanings of the numerical values given in constraint_value.  kg_solvent_water: can only be applied to H2O and units must be kg.  bulk_composition: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species but not in kinetic species, and units must be moles or mass (kg, g, etc).  bulk_composition_with_kinetic: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species and in kinetic species, and units must be moles or mass (kg, g, etc).  free_concentration: can be applied to all basis species that are not gas and not H2O and not mineral, and represents the total amount of the basis species existing freely (not as secondary species) within the solution, and units must be molal or mass_per_kg_solvent.  free_mineral: can be applied to all mineral basis species, and represents the total amount of the mineral existing freely (precipitated) within the solution, and units must be moles, mass or cm3.  activity and log10activity: can be applied to basis species that are not gas and not mineral and not sorbing sites, and represents the activity of the basis species (recall pH = -log10activity), and units must be dimensionless.  fugacity and log10fugacity: can be applied to gases, and units must be dimensionless"}>>> = "kg_solvent_water log10activity free_concentration bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition"
  constraint_unit<<<{"description": "Units of the numerical values given in constraint_value.  Dimensionless: should only be used for activity or fugacity constraints.  Moles: mole number.  Molal: moles per kg solvent water.  kg: kilograms.  g: grams.  mg: milligrams.  ug: micrograms.  kg_per_kg_solvent: kilograms per kg solvent water.  g_per_kg_solvent: grams per kg solvent water.  mg_per_kg_solvent: milligrams per kg solvent water.  ug_per_kg_solvent: micrograms per kg solvent water.  cm3: cubic centimeters"}>>> = "   kg               dimensionless molal              moles            moles            moles            moles            moles            moles            moles            moles            moles"
  max_initial_residual<<<{"description": "Attempt to alter the initial-guess molalities so that the initial residual for the Newton process (that finds the aqueous configuration) is less than this number of moles"}>>> = 0.1 # this small value is needed so that the charge-balance species is not switched to Ca++ (which, by the way, does not make a huge difference to the result)
  prevent_precipitation<<<{"description": "Mineral species in this list will be prevented from precipitating, irrespective of their saturation index (unless they are initially in the basis)"}>>> = "Nontronit-Ca Nontronit-Mg Nontronit-Na Hematite Kaolinite Beidellit-Ca Beidellit-H Beidellit-Mg Beidellit-Na Pyrophyllite Gibbsite Paragonite Quartz"
  stoichiometric_ionic_str_using_Cl_only<<<{"description": "If set to true, the stoichiometric ionic strength will be set equal to Cl- molality (or max_ionic_strength if the Cl- molality is too big).  This flag overrides ionic_str_using_basis_molality_only"}>>> = true # for comparison with GWB
  mol_cutoff<<<{"description": "Information regarding species with molalities less than this amount will not be outputted"}>>> = 1E-7
  abs_tol<<<{"description": "If the residual of the algebraic system (measured in mol) is lower than this value, the Newton process (that finds the aqueous configuration) is deemed to have converged"}>>> = 1E-15
[]

GWB input file: no minerals

The analogous GWB input file is

# React script that is the GWB equivalent to amazon.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O          = 1 free kg
HCO3-        = .0003114 mol
SiO2(aq)     = .0001165 mol
O2(aq)       = .0001813 free molal
Ca++         = .0001073 mol
SO4--        = 3.1233e-5 mol
Cl-          = .0001383 mol
balance on Cl-
Na+          = 7.83e-5 mol
Mg++         = 4.526e-5 mol
Al+++        = 2.5945e-6 mol
Fe++         = 1.0744e-6 mol
pH           = 6.5
suppress all
printout  species = long
epsilon = 1e-12

Results: no minerals

Error and charge-neutrality error

The geochemistry simulation reports an error of 3.163e-16 $var element = document.getElementById("moose-equation-7539856e-d56a-43d4-96b1-6a17e3c0dd8a");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mol, and that the charge of the solution is 0 $var element = document.getElementById("moose-equation-c881cb28-4442-426e-a017-338ae858481b");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mol.

Solution mass

The solution mass is 1.000 $var element = document.getElementById("moose-equation-68b47873-e985-432e-9ad7-a0bf150f9f30");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ kg.

Ionic strength and water activity

The ionic strength is 0.0005644 $var element = document.getElementById("moose-equation-afacef0b-8084-4bd5-aad1-ed24d5540fec");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ mol/kg(solvent water), and the water activity is 1.000.

pH, pe and Eh

The pH is 6.5, the pe is 14.07, and Eh = 0.832 $var element = document.getElementById("moose-equation-fa67e6a2-5369-499f-b5a9-b9c111e07010");katex.render("\\,", element, {displayMode:false,throwOnError:false});$ V.

Species distribution

Ignoring all mineral supersaturation, Bethke (2007) states that an equilibrium calculation of the molalities of the most abundant species results in Table 2.

Identical results are produced by the geochemistry module and the Geochemists Workbench software.

Table 2: Calculated molalities, activity coefficients and activities of the most abundant species in Amazon river water

Species	Molality (mol.kg $var element = document.getElementById("moose-equation-373ec350-0897-40e5-8e96-b0c878c24c0e");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ )	Activity coeff	log $var element = document.getElementById("moose-equation-63571e13-9e52-4c2c-9d68-e2f000c6b571");katex.render("_{10}", element, {displayMode:false,throwOnError:false});$ a
HCO $var element = document.getElementById("moose-equation-bc2adc96-24d5-4568-877e-9073f6a833fe");katex.render("_{3}^{-}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-382a75e1-b69b-489a-aaf8-b3f58e60707c");katex.render("0.181\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	0.974	-3.75
O $var element = document.getElementById("moose-equation-e0d0d819-de5d-4ced-949f-44a1831ebf29");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ (aq)	$var element = document.getElementById("moose-equation-fadafe2d-982b-47bd-b0c5-9144dd68258d");katex.render("0.181\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	1.000	-3.74
Cl $var element = document.getElementById("moose-equation-760b0075-f237-47d0-9e30-f1ee0307372a");katex.render("^{-}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-026a2aa1-e5ea-470e-812c-d9a8a6ce5881");katex.render("0.138\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	0.973	-3.87
CO $var element = document.getElementById("moose-equation-0da99489-7565-489a-afde-35519c389dd5");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ (aq)	$var element = document.getElementById("moose-equation-ec98699d-1593-465c-a372-c7058f9251a2");katex.render("0.130\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	1.000	-3.89
SiO $var element = document.getElementById("moose-equation-e8ce12f8-a248-4c5e-92d4-82afdd3202a2");katex.render("_{2}", element, {displayMode:false,throwOnError:false});$ (aq)	$var element = document.getElementById("moose-equation-5f8ace94-302b-4fa2-b89e-1986f6472454");katex.render("0.116\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	1.000	-3.93
Ca $var element = document.getElementById("moose-equation-16da1961-a743-4cd1-9028-777e7251ed2e");katex.render("^{2+}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-20425682-6d1a-42e1-850a-fb67bf78446c");katex.render("0.106\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	0.899	-4.02
Na $var element = document.getElementById("moose-equation-a9237c1a-b6b1-4db3-85ef-29530a32657a");katex.render("^{+}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-0677744d-30a2-4266-bd5a-817b4e077d4b");katex.render("0.783\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$	0.973	-4.12
Mg $var element = document.getElementById("moose-equation-d46df615-86ec-4a5b-85a7-e370c2ec74ae");katex.render("^{2+}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-7b65d3d3-64db-4328-913e-2b0121fbb4e7");katex.render("0.450\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$	0.901	-4.39
SO $var element = document.getElementById("moose-equation-43d2c37f-870f-49de-a8df-e0bb80779261");katex.render("_{4}^{2-}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-e603df50-ac26-4f33-9534-20e4eff50c8c");katex.render("0.305\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$	0.898	-4.56

Minerals

Both the geochemistry module and the GWB software predict identical saturation indices for minerals. They predict the saturation indices of the equilibrium solution in Table 2 are greater than 0 for the minerals: nontronite, hematite, kaolinite, beidellite, pyrophyllite, gibbsite, paragonite and quartz (see Figure 6.3 of Bethke (2007)).

Fixing mineral volumes

An alternate model assumes the water is in equilibrium with kaolinite and hematite, and that these control the aluminum and iron concentrations. That is, instead of Al $var element = document.getElementById("moose-equation-49431081-2746-4395-b53f-0c7208f382a2");katex.render("^{3+}", element, {displayMode:false,throwOnError:false});$ , kaolinite is used in the basis, and instead of Fe $var element = document.getElementById("moose-equation-9abcd07e-a37c-4656-b954-28f609360c35");katex.render("^{2+}", element, {displayMode:false,throwOnError:false});$ , hematite is used in the basis. The free volumes of each of these is rather unimportant, but Bethke (2007) suggests using the values in Table 3.

Table 3: Alternate model of Amazon river water assumes equilibrium with kaoline and hematite with specified amount

Mineral	Free amount (cm $var element = document.getElementById("moose-equation-12b43e1e-8b6f-4a54-ae51-d596ef49c67d");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ )
Kaolinite	1
Hematite	1

MOOSE input file: with minerals

[TimeIndependentReactionSolver<<<{"href": "../../../syntax/TimeIndependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  swap_out_of_basis<<<{"description": "Species that should be removed from the model_definition's basis and be replaced with the swap_into_basis species"}>>> = "Al+++ Fe++"
  swap_into_basis<<<{"description": "Species that should be removed from the model_definition's equilibrium species list and added to the basis.  There must be the same number of species in swap_out_of_basis and swap_into_basis.  These swaps are performed before any other computations during the initial problem setup. If this list contains more than one species, the swapping is performed one-by-one, starting with the first pair (swap_out_of_basis[0] and swap_into_basis[0]), then the next pair, etc"}>>> = "Kaolinite Hematite"
  charge_balance_species<<<{"description": "Charge balance will be enforced on this basis species.  This means that its bulk mole number may be changed from the initial value you provide in order to ensure charge neutrality.  After the initial swaps have been performed, this must be in the basis, and it must be provided with a bulk_composition constraint_meaning."}>>> = "Cl-"
  constraint_species<<<{"description": "Names of the species that have their values fixed to constraint_value with meaning constraint_meaning.  All basis species (after swap_into_basis and swap_out_of_basis) must be provided with exactly one constraint.  These constraints are used to compute the configuration during the initial problem setup, and in time-dependent simulations they may be modified as time progresses."}>>> = "H2O              H+            O2(aq)             SiO2(aq)         Kaolinite    Hematite     Ca++             Mg++             Na+              HCO3-            SO4--            Cl-      "
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -6.5          1.813E-4            0.0001165       1E-2         0.033        0.0001073        4.526E-5         7.83E-5          0.0003114        3.1233E-5        1.383E-4"
  constraint_meaning<<<{"description": "Meanings of the numerical values given in constraint_value.  kg_solvent_water: can only be applied to H2O and units must be kg.  bulk_composition: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species but not in kinetic species, and units must be moles or mass (kg, g, etc).  bulk_composition_with_kinetic: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species and in kinetic species, and units must be moles or mass (kg, g, etc).  free_concentration: can be applied to all basis species that are not gas and not H2O and not mineral, and represents the total amount of the basis species existing freely (not as secondary species) within the solution, and units must be molal or mass_per_kg_solvent.  free_mineral: can be applied to all mineral basis species, and represents the total amount of the mineral existing freely (precipitated) within the solution, and units must be moles, mass or cm3.  activity and log10activity: can be applied to basis species that are not gas and not mineral and not sorbing sites, and represents the activity of the basis species (recall pH = -log10activity), and units must be dimensionless.  fugacity and log10fugacity: can be applied to gases, and units must be dimensionless"}>>> = "kg_solvent_water log10activity free_concentration bulk_composition free_mineral free_mineral bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition"
  constraint_unit<<<{"description": "Units of the numerical values given in constraint_value.  Dimensionless: should only be used for activity or fugacity constraints.  Moles: mole number.  Molal: moles per kg solvent water.  kg: kilograms.  g: grams.  mg: milligrams.  ug: micrograms.  kg_per_kg_solvent: kilograms per kg solvent water.  g_per_kg_solvent: grams per kg solvent water.  mg_per_kg_solvent: milligrams per kg solvent water.  ug_per_kg_solvent: micrograms per kg solvent water.  cm3: cubic centimeters"}>>> = "   kg               dimensionless molal              moles            moles        moles        moles            moles            moles            moles            moles            moles"
  prevent_precipitation<<<{"description": "Mineral species in this list will be prevented from precipitating, irrespective of their saturation index (unless they are initially in the basis)"}>>> = "Nontronit-Ca Nontronit-Mg Nontronit-Na Hematite Kaolinite Beidellit-Ca Beidellit-H Beidellit-Mg Beidellit-Na Pyrophyllite Gibbsite Paragonite Quartz"
  ramp_max_ionic_strength_initial<<<{"description": "The number of iterations over which to progressively increase the maximum ionic strength (from zero to max_ionic_strength) during the initial equilibration.  Increasing this can help in convergence of the Newton process, at the cost of spending more time finding the aqueous configuration."}>>> = 0 # not needed in this simple problem
  stoichiometric_ionic_str_using_Cl_only<<<{"description": "If set to true, the stoichiometric ionic strength will be set equal to Cl- molality (or max_ionic_strength if the Cl- molality is too big).  This flag overrides ionic_str_using_basis_molality_only"}>>> = true # for comparison with GWB
  mol_cutoff<<<{"description": "Information regarding species with molalities less than this amount will not be outputted"}>>> = 1E-7
  abs_tol<<<{"description": "If the residual of the algebraic system (measured in mol) is lower than this value, the Newton process (that finds the aqueous configuration) is deemed to have converged"}>>> = 1E-15
[]

[UserObjects<<<{"href": "../../../syntax/UserObjects/index.html"}>>>]
  [definition]
    type = GeochemicalModelDefinition<<<{"description": "User object that parses a geochemical database file, and only retains information relevant to the current geochemical model", "href": "../../../source/userobjects/GeochemicalModelDefinition.html"}>>>
    database_file<<<{"description": "The name of the geochemical database file"}>>> = "../../../database/moose_geochemdb.json"
    basis_species<<<{"description": "A list of basis components relevant to the aqueous-equilibrium problem. H2O must appear first in this list.  These components must be chosen from the 'basis species' in the database, the sorbing sites (if any) and the decoupled redox states that are in disequilibrium (if any)."}>>> = "H2O H+ SiO2(aq) Al+++ Fe++ Ca++ Mg++ Na+ HCO3- SO4-- Cl- O2(aq)"
    equilibrium_minerals<<<{"description": "A list of minerals that are in equilibrium with the aqueous solution.  All members of this list must be in the 'minerals' section of the database file"}>>> = "Nontronit-Ca Nontronit-Mg Nontronit-Na Hematite Kaolinite Beidellit-Ca Beidellit-H Beidellit-Mg Beidellit-Na Pyrophyllite Gibbsite Paragonite Quartz"
    piecewise_linear_interpolation<<<{"description": "If true then use a piecewise-linear interpolation of logK and Debye-Huckel parameters, regardless of the interpolation type specified in the database file.  This can be useful for comparing with results using other geochemistry software"}>>> = true # for comparison with GWB
  []
[]

GWB input file: with minerals

# React script that is the GWB equivalent to amazon_with_minerals.rea
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O          = 1 free kg
HCO3-        = .0003114 mol
SiO2(aq)     = .0001165 mol
O2(aq)       = .0001813 free molal
Ca++         = .0001073 mol
SO4--        = 3.1233e-5 mol
Cl-          = .0001383 mol
balance on Cl-
Na+          = 7.83e-5 mol
Mg++         = 4.526e-5 mol
swap Kaolinite for Al+++
Kaolinite    = .01 free mol
swap Hematite for Fe++
Hematite     = .033 free mol
pH           = 6.5
suppress ALL
unsuppress  Hematite Kaolinite
printout  species = long
epsilon = 1e-15

Results: with minerals

According to this model, only nontronite clay is supersaturated by any significant amount. The geochemistry module produces identical results to the GWB software.

References

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

@book{bethke_2007,
    author = "Bethke, Craig M.",
    place = "Cambridge",
    edition = "2",
    title = "Geochemical and Biogeochemical Reaction Modeling",
    DOI = "10.1017/CBO9780511619670",
    publisher = "Cambridge University Press",
    year = "2007"
}

(moose/modules/geochemistry/test/tests/equilibrium_models/amazon.i)

[TimeIndependentReactionSolver]
  model_definition = definition
  charge_balance_species = "Cl-"
  constraint_species = "H2O              H+            O2(aq)             SiO2(aq)         Al+++            Fe++             Ca++             Mg++             Na+              HCO3-            SO4--            Cl-      "
  constraint_value = "  1.0              -6.5          1.813E-4           0.0001165        2.5945E-6        1.0744E-6        0.0001073        4.526E-5         7.83E-5          0.0003114        3.1233E-5        1.383E-4"
  constraint_meaning = "kg_solvent_water log10activity free_concentration bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition"
  constraint_unit = "   kg               dimensionless molal              moles            moles            moles            moles            moles            moles            moles            moles            moles"
  max_initial_residual = 0.1 # this small value is needed so that the charge-balance species is not switched to Ca++ (which, by the way, does not make a huge difference to the result)
  prevent_precipitation = "Nontronit-Ca Nontronit-Mg Nontronit-Na Hematite Kaolinite Beidellit-Ca Beidellit-H Beidellit-Mg Beidellit-Na Pyrophyllite Gibbsite Paragonite Quartz"
  stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
  mol_cutoff = 1E-7
  abs_tol = 1E-15
[]

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ SiO2(aq) Al+++ Fe++ Ca++ Mg++ Na+ HCO3- SO4-- Cl- O2(aq)"
    equilibrium_minerals = "Nontronit-Ca Nontronit-Mg Nontronit-Na Hematite Kaolinite Beidellit-Ca Beidellit-H Beidellit-Mg Beidellit-Na Pyrophyllite Gibbsite Paragonite Quartz"
    piecewise_linear_interpolation = true # for comparison with GWB
  []
[]