Cooling a solution in contact with feldspars

MOOSE input file: model definition

The MOOSE input file defines the basis species and equilibrium minerals via the GeochemicalModelDefinition. Note the appearance of the remove_all_extrapolated_secondary_species flag, which is necessary to ensure convergence of this problem. The offending secondary species is Al13O4(OH)24(7+) that has defined only at low temperatures in the original database. The default process to use for such species is to extrapolate the equilibrium constants to higher temperatures. In this case at C which prevents geochemistry from converging, so the species is removed (it is also removed by the Geochemists Workbench software).

[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+ Cl- Al+++ K+ SiO2(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"}>>> = "Albite Maximum Muscovite Quartz"
    remove_all_extrapolated_secondary_species<<<{"description": "After reading the database file, immediately remove all secondary species that have extrapolated equilibrium constants.  Sometimes these extrapolations are completely crazy and those secondary species greatly impact the results"}>>> = true # this removes Al13O4(OH)24(7+) that has extreme logK values
  []
[]

MOOSE input file: model time-dependence

The TimeDependentReactionSolver defines:

• the swaps mentioned above

• the charge-balance species

• the constraints on the species, which in this case hold just at the initial time because the system is closed after that time (no more minerals are added after the initialization, but they could be added during the initialization to ensure the correct number of free moles)

• the initial temperature

• the temperature as a function of time

• other flags that allow accurate comparison with the Geochemists Workbench software

[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"}>>> = "Al+++ K+ H+ SiO2(aq)"
  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"}>>> = "Albite Maximum Muscovite Quartz"
  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              Muscovite    Na+              Cl-              Albite       Maximum      Quartz"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              5            1.14093          1.14093          20           10           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 bulk_composition bulk_composition free_mineral 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               cm3          moles            moles            cm3          cm3          cm3"
  initial_temperature<<<{"description": "The initial aqueous solution is equilibrated at this system before adding reactants, changing temperature, etc."}>>> = 300
  temperature<<<{"description": "Temperature.  This has two different meanings if mode!=4.  (1) If no species are being added to the solution (no source_species_rates are positive) then this is the temperature of the aqueous solution.  (2) If species are being added, this is the temperature of the species being added.  In case (2), the final aqueous-solution temperature is computed assuming the species are added, temperature is equilibrated and then, if species are also being removed, they are removed.  If you wish to add species and simultaneously alter the temperature, you will have to use a sequence of heat-add-heat-add, etc steps.  In the case that mode=4, temperature is the final temperature of the aqueous solution"}>>> = temperature
  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-14
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output for this example
[]

The time-dependent temperature is defined through an AuxVariable:

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [temperature]
  []
[]
[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [temperature]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = temperature
    function<<<{"description": "The function to use as the value"}>>> = '300 - t'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin' # so that it is correct when we solve the system
  []
[]

along with:

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  start_time = -10 # so that the output at 300degC is easily captured
  dt = 10
  end_time = 275
[]

MOOSE input file: recording the desired quantities

The desired mineral volumes and solution temperature are recorded into Postprocessors using the AuxVariables automatically included by the TimeDependentReactionSolver:

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [solution_temperature]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'temperature'
  []
  [cm3_Max_Micro]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Maximum'
  []
  [cm3_Albite]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Albite'
  []
  [cm3_Muscovite]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Muscovite'
  []
  [cm3_Quartz]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Quartz'
  []
[]

GWB input file

The equivalent Geochemists Workbench input file is

# React script that is equivalent to cooling.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature initial = 300 C, final = 25 C
H2O          = 1 free kg
swap Albite for Al+++
Albite       = .19986 free mol
swap "Maximum Microcline" for K+
"Maximum Microcline" = .09196 free mol
swap Muscovite for H+
Muscovite    = .03553 free mol
swap Quartz for SiO2(aq)
Quartz       = .08815 free mol
Na+          = 1 molal
Cl-          = 1 molal
balance on Cl-

suppress  "Albite low"
printout  species = long
epsilon = 1e-14

Results

Figure 1: Precipitated volumes as a function of temperature

References