Kinetically-controlled dissolution of quartz

MOOSE input file

The MOOSE input file defines the model using the GeochemicalModelDefinition. This defines the basis species as well as defining that the dynamics of the mineral Quartz will be controlled by a kinetic rate law. Note that this model contains H and Cl, which is a bit different than Bethke's set up (these species are provided with very small bulk composition so they don't impact the result). The reason for this is that geochemistry requires a charge-balance species to be defined.

  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O SiO2(aq) H+ Cl-"
    kinetic_minerals = "Quartz"
    kinetic_rate_descriptions = "rate_quartz"
    piecewise_linear_interpolation = true # for comparison with GWB

The rate law for Quartz is defined by a GeochemistryKineticRate UserObject:

  [rate_quartz]
    type = GeochemistryKineticRate
    kinetic_species_name = Quartz
    intrinsic_rate_constant = 1.728E-10 # 2.0E-15mol/s/cm^2 = 1.728E-10mol/day/cm^2
    multiply_by_mass = true
    area_quantity = 1000
  []

The TimeDependentReactionSolver defines the initial concentration of the species, including the initial mole number for Quartz. The system is closed at time zero (by default) so the free_molality constraint of SiO(aq) becomes inactive (no SiO(aq) is added or removed from the system by an external agent after this time):

[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
  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+               Cl-              SiO2(aq)"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              1E-10            1E-10            1E-9"
  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 bulk_composition bulk_composition free_concentration"
  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               moles            moles            molal"
  initial_temperature<<<{"description": "The initial aqueous solution is equilibrated at this system before adding reactants, changing temperature, etc."}>>> = 100.0
  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"}>>> = 100.0
  kinetic_species_name<<<{"description": "Names of the kinetic species given initial values in kinetic_species_initial_value"}>>> = Quartz
  kinetic_species_initial_value<<<{"description": "Initial number of moles, mass or volume (depending on kinetic_species_unit) for each of the species named in kinetic_species_name"}>>> = 5
  kinetic_species_unit<<<{"description": "Units of the numerical values given in kinetic_species_initial_value.  Moles: mole number.  kg: kilograms.  g: grams.  mg: milligrams.  ug: micrograms.  cm3: cubic centimeters"}>>> = kg
  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 # max_ionic_strength in such a simple problem does not need ramping
  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
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output for this example
[]

and the Executioner defines the time-stepping (time is measured in days in this input file)

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  [TimeStepper<<<{"href": "../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = timestepper
  []
  end_time = 5.0
[]

An AuxVariable, AuxKernel, Postprocessor and Output allow the mole number of the Quartz mineral to be recorded into a CSV file using the moles_Quartz AuxVariable added automatically by the TimeDependentReactionSolver

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [diss]
  []
[]
[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [diss]
    type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../../source/auxkernels/ParsedAux.html"}>>>
    coupled_variables<<<{"description": "Vector of coupled variable names"}>>> = moles_Quartz
    expression<<<{"description": "Parsed function expression to compute"}>>> = '83.216414271 - moles_Quartz'
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = diss
  []
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [dissolved_moles]
    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."}>>> = diss
  []
[]
[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

GWB input file

The equivalent Geochemists Workbench input file is

# React script that is equivalent to quartz_dissolution.i
time begin = 0 days, end = 5 days
T = 100
SiO2(aq) = 1 umolal
react 5000 g Quartz
kinetic Quartz rate_con = 2.0e-15 surface = 1000
go

Results

The results shown below can be compared with Bethke (2007) Figure 16.1.

Figure 1: Change in mole number of kinetically-controlled quartz. Compare with Bethke's Figure 16.1

References