Chemical models of Morro de Ferro groundwater

MOOSE input file

The MOOSE input file contains the GeochemicalModelDefinition and TimeIndependentReactionSolver. The bulk mole number of the aqueous species is also fixed appropriately in the latter. The numbers are different than the concentration in mg.l given in the above table, and may be worked out using the TDS. The other flags options ensure nice convergence and 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
  temperature<<<{"description": "The temperature (degC) of the aqueous solution"}>>> = 22
  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-" # this means the bulk moles of Cl- will not be exactly as set below
  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)             Cl-               HCO3-             Ca++              Mg++              Na+               K+                Fe++              Mn++              Zn++              SO4--"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -6.05         0.13438E-3         3.041E-5          0.0295E-3         0.005938E-3       0.01448E-3        0.0018704E-3      0.005115E-3       0.01307E-3        0.005042E-3       0.001897E-3       0.01562E-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  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             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"}>>> = 1E-2
  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."}>>> = 10
  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-5
  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+ Cl- O2(aq) HCO3- Ca++ Mg++ Na+ K+ Fe++ Mn++ Zn++ SO4--"
  []
[]

Geochemists Workbench input file

The equivalent GWB input file is

# React script that is equivalent to the morro.i MOOSE input file
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 22 C
H2O          = 1 free kg
Cl-          = 3.041E-5 mol
balance on Cl-
H+           = 6.05 pH
O2(aq)       = 0.13438E-3 free molal
HCO3-        = 0.0295E-3 mol
Ca++         = 0.005938E-3 mol
Mg++         = 0.01448E-3 mol
Na+          = 0.0018704E-3 mol
K+           = 0.005115E-3 mol
Fe++         = 0.01307E-3 mol
Mn++         = 0.005042E-3 mol
Zn++         = 0.001897E-3 mol
SO4--        = 0.01562E-4 mol
printout  species = long
suppress all
epsilon = 1e-13
go

MOOSE input file

The MOOSE input file is very similar to the redox-equilibrium case. The differences are:

• Fe is included in the basis;

• the sum of the bulk composition for Fe and Fe equals the bulk composition of Fe in the equilibrium case

[TimeIndependentReactionSolver<<<{"href": "../../../syntax/TimeIndependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  temperature<<<{"description": "The temperature (degC) of the aqueous solution"}>>> = 22
  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-" # this means the bulk moles of Cl- will not be exactly as set below
  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)             Cl-              HCO3-            Ca++             Mg++             Na+              K+               Fe++            Fe+++             Mn++             Zn++             SO4--"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -6.05         0.13438E-3         3.041E-5         0.0295E-3        0.005938E-3      0.01448E-3       0.0018704E-3     0.005115E-3      0.012534E-3     0.0005372E-3      0.005042E-3      0.001897E-3      0.01562E-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 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            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"}>>> = 1E-2
  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."}>>> = 10
  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-5
  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+ Cl- O2(aq) HCO3- Ca++ Mg++ Na+ K+ Fe++ Fe+++ Mn++ Zn++ SO4--"
  []
[]

Geochemists Workbench input file

The equivalent GWB input file is

# React script that is equivalent to the morro_disequilibrium.i MOOSE input file
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 22 C
decouple Fe+++
H2O          = 1 free kg
Cl-          = 3.041E-5 mol
balance on Cl-
H+           = 6.05 pH
O2(aq)       = 0.13438E-3 free molal
HCO3-        = 0.0295E-3 mol
Ca++         = 0.005938E-3 mol
Mg++         = 0.01448E-3 mol
Na+          = 0.0018704E-3 mol
K+           = 0.005115E-3 mol
Fe++         = 0.012534E-3 mol
Fe+++        = 5.372E-7 mol
Mn++         = 0.005042E-3 mol
Zn++         = 0.001897E-3 mol
SO4--        = 0.01562E-4 mol
printout  species = long
suppress all
epsilon = 1e-13
go

Error and charge-neutrality error

The geochemistry simulations report an error of less than mol, and that the charge of the solution has magnitude less than mol.

Solution mass

The solution mass is 1.000kg.

Ionic strength and water activity

The ionic strength is 1.01E-4mol/kg(solvent water) for the equilibrium case and 1.29E-4mol/kg(solvent water) for the disequilibrium case. The water activity is 1.000 in both cases

pH and pe

The pH is 6.05 in both cases (as specified by the constraint), the pe is 14.7 in both cases.

Nernst Eh values

Bethke (2007) computes the Nernst Eh values for the redox half-reactions as: Both codes produce this result. The latter one is only relevant in the redox-disequilibrium case

Aqueous species distribution

The species distribution predicted by Bethke (2007) is shown in the right-hand column of Table 2. This result is obtained by ignoring any potential precipitation of minerals. Both codes predict the same result, up to precision.