Dissolution of pyrite, with and without fixed oxygen fugacity

MOOSE input file

The TimeDependentReactionSolver involves specifying:

• the swaps

• the initial free mole number of Hematite, the pH via the H activity, the bulk mole composition of the aqueous species, and the O2(g) fugacity

• the system is closed at (the default)

• that at the activity constraints are removed from H (so the pH can vary during the simulation because an external agent is no longer fixing it) and O2(g) (meaning the external agent controlling the O2(g) fugacity is no longer adding/removing O2(g) from the aqueous solution)

• that 1mg/second of pyrite is added to the system

[TimeDependentReactionSolver<<<{"href": "../../../syntax/TimeDependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  geochemistry_reactor_name<<<{"description": "The name that will be given to the GeochemistryReactor UserObject built by this action"}>>> = reactor
  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"}>>> = "O2(aq) 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"}>>> = "O2(g) 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              Hematite     H+            Ca++             Mg++             Na+              HCO3-            SO4--            Cl-              O2(g)"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              1            -6.5          4                1                2                18               3                5                0.2"
  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 free_mineral log10activity bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition fugacity"
  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               mg           dimensionless mg               mg               mg               mg               mg               mg               dimensionless"
  remove_fixed_activity_name<<<{"description": "The name of the species that should have their activity or fugacity constraint removed at time given in remove_fixed_activity_time.  There should be an equal number of these names as times given in remove_fixed_activity_time.  Each of these must be in the basis and have an activity or fugacity constraint"}>>> = "H+ O2(g)"
  remove_fixed_activity_time<<<{"description": "The times at which the species in remove_fixed_activity_name should have their activity or fugacity constraint removed."}>>> = '0  0'
  source_species_names<<<{"description": "The name of the species that are added at rates given in source_species_rates.  There must be an equal number of these as source_species_rates."}>>> = "Pyrite"
  source_species_rates<<<{"description": "Rates, in mols/time_unit, of addition of the species with names given in source_species_names.  A negative value corresponds to removing a species: be careful that you don't cause negative mass problems!"}>>> = 8.336E-6 # = 1mg(pyrite)/second, 1mg(pyrite) = 8.34E-6
  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
  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-13
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output is required
[]

GWB input file

The equivalent Geochemists Workbench input file is

# React script that is equivalent to dissolution_pyrite_1.i
Ca++ = 9.98e-5 mol
Mg++ = 4.114e-5 mol
Na+ = 8.7e-5 mol
HCO3- = 2.95e-4 mol
SO4-- = 3.123e-5 mol
Cl- = 1.4103e-5 mol
balance on Cl-
H+ = 6.5 pH
swap O2(g) for O2(aq)
f O2(g) = 0.2
swap Hematite for Fe++
Hematite = 6.26e-6 free mol
epsilon = 1e-13
react 8.336e-5 mol Pyrite
go

Results

The following graphs demonstrate agreement between the geochemistry module and the GWB software, which may also be compared with Bethke (2007) Figure 14.2 and 14.3.

Figure 1: Precipitated volumes as pyrite is reacted in a system that is closed to oxygen

Figure 2: Solution pH as pyrite is reacted in a system that is closed to oxygen

Figure 3: Species concentrations as pyrite is reacted in a system that is closed to oxygen

MOOSE input file

The TimeDependentReactionSolver is almost identical to Case 1. The differences are:

• that at the activity constraints are removed from H, but not O2(g)

• that 1000mg/second of pyrite is added to the system

[TimeDependentReactionSolver<<<{"href": "../../../syntax/TimeDependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  geochemistry_reactor_name<<<{"description": "The name that will be given to the GeochemistryReactor UserObject built by this action"}>>> = reactor
  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"}>>> = "O2(aq) 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"}>>> = "O2(g) 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              Hematite     H+            Ca++             Mg++             Na+              HCO3-            SO4--            Cl-              O2(g)"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              1            -6.5          4                1                2                18               3                5                0.2"
  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 free_mineral log10activity bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition fugacity"
  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               mg           dimensionless mg               mg               mg               mg               mg               mg               dimensionless"
  remove_fixed_activity_name<<<{"description": "The name of the species that should have their activity or fugacity constraint removed at time given in remove_fixed_activity_time.  There should be an equal number of these names as times given in remove_fixed_activity_time.  Each of these must be in the basis and have an activity or fugacity constraint"}>>> = "H+"
  remove_fixed_activity_time<<<{"description": "The times at which the species in remove_fixed_activity_name should have their activity or fugacity constraint removed."}>>> = '0 '
  source_species_names<<<{"description": "The name of the species that are added at rates given in source_species_rates.  There must be an equal number of these as source_species_rates."}>>> = "Pyrite"
  source_species_rates<<<{"description": "Rates, in mols/time_unit, of addition of the species with names given in source_species_names.  A negative value corresponds to removing a species: be careful that you don't cause negative mass problems!"}>>> = 8.336E-4 # = 0.1g(pyrite)/second over 10 seconds, 1g(pyrite) = 8.34E-3
  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
  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-13
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output is required
[]

Results

The following graphs demonstrate agreement between the geochemistry module and the GWB software, which may also be compared with Bethke (2007) Figure 14.4.

Figure 4: Precipitated volumes as pyrite is reacted in a system that is open to oxygen

Figure 5: Solution pH as pyrite is reacted in a system that is open to oxygen

The two cases are clearly dramatically different. In Case 1, about 8mg of pyrite dissolves, precipitating hematite, before hematite suddenly dissolves and no further pyrite dissolution occurs. In Case 2, pyrite dissolution and hematite precipitation continues indefinitely in a rather acidic solution.