Kinetically-controlled dissolution of albite into an acidic solution

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 Albite will be controlled by a kinetic rate law.

  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Cl- Na+ SiO2(aq) Al+++"
    kinetic_minerals = "Albite"
    kinetic_rate_descriptions = "rate_albite"

The rate law for Albite is defined by a GeochemistryKineticRate UserObject (note the promoting_species):

  [rate_albite]
    type = GeochemistryKineticRate
    kinetic_species_name = Albite
    intrinsic_rate_constant = 5.4432E-8 # 6.3E-13mol/s/cm^2 = 5.4432E-8mol/day/cm^2
    multiply_by_mass = true
    area_quantity = 1000
    promoting_species_names = "H+"
    promoting_indices = "1.0"
  []

The TimeDependentReactionSolver defines the following.

• The initial concentration of the species (see note below)

• The initial mole number for Albite

• That the system is closed at time zero (by default) so the free_molality constraints becomes inactive (no SiO(aq) or Al are added or removed from the system by an external agent after this time)

• The pH, via the activity constraint on H. This constraint is not removed, so this effectively means HCl is continually added or removed from the system to maintain the pH (remember Cl is the charge-balance species)

• That the kinetic rates are updated during the Newton solve that finds the equilibrium configuration at each time step. This helps with stability since it is an implicit approach. Geochemists Workbench appears to only compute the kinetic rates at the beginning of the time step (an explicit approach).

[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-              Na+              SiO2(aq)           Al+++"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -1.5          0.1              0.1              1E-6               1E-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 free_concentration 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               dimensionless moles            moles            molal              molal"
  initial_temperature<<<{"description": "The initial aqueous solution is equilibrated at this system before adding reactants, changing temperature, etc."}>>> = 70.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"}>>> = 70.0
  kinetic_species_name<<<{"description": "Names of the kinetic species given initial values in kinetic_species_initial_value"}>>> = Albite
  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"}>>> = 250
  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"}>>> = g
  evaluate_kinetic_rates_always<<<{"description": "If true, then, evaluate the kinetic rates at every Newton step during the solve using the current values of molality, activity, etc (ie, implement an implicit solve).  If false, then evaluate the kinetic rates using the values of molality, activity, etc, at the start of the current time step (ie, implement an explicit solve)"}>>> = true # implicit time-marching used for stability
  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
[]

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

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  dt = 5
  end_time = 30 # measured in days
[]

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

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [mole_change_albite]
  []
[]
[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [mole_change_albite]
    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_Albite
    expression<<<{"description": "Parsed function expression to compute"}>>> = 'moles_Albite - 0.953387'
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = mole_change_albite
  []
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [mole_change_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."}>>> = "mole_change_albite"
  []
[]
[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 kinetic_albite.i
time begin = 0 days, end = 30 days
T = 70
pH = 1.5
0.1 molal Cl-
0.1 molal Na+
1 umolal SiO2(aq)
1 umolal Al+++

react 250 g "Albite low"
kinetic "Albite low" rate_con = 6.3e-13, apower(H+) = 1 \
   surface = 1000
fix pH
go

Note that the bulk composition for Na+ is 0.1, as mentioned above. When running, GWB increases the bulk composition of Cl- in order to enforce charge neutrality (which looks strange because the initial composition is uncharged) but that does not impact the results presented below.

Results

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

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

References