Chemical model of Red Sea brine

MOOSE input file: no precipitation or dissolution

In order to estimate the brine's oxidation state, Bethke (2007) recommends using the mineral Sphalerite instead of primary species O2(aq), and the mineral Barite instead of primary aqueous species Ba++. Bethke (2007) uses a small initial value of g of free amounts of each of these minerals to avoid dissolution when considering precipitation/dissolution later in the calculation.

The MOOSE input file contains the usual GeochemicalModelDefinition that specifies the database file to use, and in this case the basis species and equilibrium minerals. Various minerals are included in order to make the comparison to the precipitation+dissolution case (below) clearer. The TimeIndependentReactionSolver will prevent all mineral precipitation will be forbidden, except the Sphalerite and Barite which are constrained to have a fixed number of free moles, so the additional minerals do not impact the system whatsoever, except their saturation indices will be computed after the solution is found. 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+ Na+ K+ Mg++ Ca++ Cl- SO4-- HCO3- Cu+ F- Fe++ Pb++ Zn++ O2(aq) Ba++"
    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"}>>> = "Sphalerite Barite Fluorite Chalcocite Bornite Chalcopyrite Pyrite Galena Covellite"
    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 swaps are defined so that the minerals Sphalerite and Barite can be provided with a small number of free moles.

• The pH is fixed using the activity of H.

• There is 1kg of solvent water.

• The bulk mole number of the aqueous species is also fixed appropriately. The numbers are different than the concentration in mg.kg 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
  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) Ba++"
  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"}>>> = "Sphalerite Barite"
  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)"}>>> = "Sphalerite Barite Fluorite Chalcocite Bornite Chalcopyrite Pyrite Galena Covellite"
  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+            Na+              K+               Mg++             Ca++             Cl-              SO4--            HCO3-            Cu+              F-               Fe++            Pb++              Zn++             Sphalerite   Barite"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -5.6          5.42             0.0643           0.0423           0.173            5.89             0.0118           0.00309          5.50E-06         0.000354         0.00195          4.09E-06         0.000111         1E-11        0.5E-11"
  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 bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition free_mineral free_mineral"
  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 moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles        moles"
  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 example
  temperature<<<{"description": "The temperature (degC) of the aqueous solution"}>>> = 60
  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-12
[]

Geochemists Workbench input file: no precipitation or dissolution

The analogous Geochemists Workbench input file for this case is

# React script that is equivalent to red_sea_no_precip.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 60 C
H2O          = 1 free kg
Cl-          = 5.89 mol
balance on Cl-
Na+          = 5.42 mol
Ca++         = 0.173 mol
K+           = 0.0643 mol
SO4--        = 0.0118 mol
Mg++         = 0.0423 mol
HCO3-        = 0.00309 mol
Fe++         = 0.00195 mol
Zn++         = 0.000111 mol
F-           = 0.000354 mol
swap Barite for Ba++
Barite       = 0.5e-11 free mol
Pb++         = 4.09e-06 mol
Cu+          = 5.50e-6 mol
pH           = 5.6
swap Sphalerite for O2(aq)
Sphalerite   = 1e-11 free mol
suppress ALL
unsuppress  Barite Sphalerite
printout  species = long
epsilon = 1e-12

Results: no precipitation or dissolution

The geochemistry results match those from Geochemists Workbench exactly.

Mineral precipitation and dissolution

The results are slightly different when minerals are allowed to precipitate and the Sphalerite and Barite are allowed to dissolve.

MOOSE input file: allowing precipitation or dissolution

The MOOSE input file is almost identical to the one above

• the prevent_precipitation option is removed

• a bulk number of moles are provided to Sphalerite and Barite so that they can dissolve. If a free number of moles had been specified then geochemistry assumes that moles of the mineral are added or removed by an external agent in order to keep the free number as specified. Hence, specifying free mole numbers prevents any dissolution. The bulk number of moles is set to the amount predicted by red_sea_no_precip.i

[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"}>>> = "O2(aq) Ba++"
  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"}>>> = "Sphalerite Barite"
  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+            Na+              K+               Mg++             Ca++             Cl-              SO4--            HCO3-            Cu+              F-               Fe++            Pb++              Zn++             Sphalerite       Barite"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -5.6          5.42             0.0643           0.0423           0.173            5.89             0.0118           0.00309          5.50E-06         0.000354         0.00195          4.09E-06         0.000111         5.87E-8          9.772E-6"
  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 bulk_composition bulk_composition bulk_composition 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 moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles"
  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
  temperature<<<{"description": "The temperature (degC) of the aqueous solution"}>>> = 60
  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-12
[]

Geochemists Workbench input file: allowing precipitation or dissolution

The analogous Geochemists Workbench input file for this case is

# React script that is equivalent to red_sea_precip.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 60 C
H2O          = 1 free kg
Cl-          = 5.89 mol
balance on Cl-
Na+          = 5.42 mol
Ca++         = 0.173 mol
K+           = 0.0643 mol
SO4--        = 0.0118 mol
Mg++         = 0.0423 mol
HCO3-        = 0.00309 mol
Fe++         = 0.00195 mol
Zn++         = 0.000111 mol
F-           = 0.000354 mol
swap Barite for Ba++
Barite       = 0.5e-11 free mol
Pb++         = 4.09e-06 mol
Cu+          = 5.50e-6 mol
pH           = 5.6
swap Sphalerite for O2(aq)
Sphalerite   = 1e-11 free mol
printout  species = long
epsilon = 1e-14

Results: allowing precipitation and dissolution

Allowing the minerals to precipitate, both codes and Bethke (2007) predicts that the Sphalerite dissolves and that 3 minerals precipitate in small quantities, as shown in Table 3

References