Quartz deposition in a fracture

MOOSE input file: stage 1

The MOOSE simulation is in 2 stages. The first determines the molality of SiO2(aq) that is in equilibrium with the quartz. This is necessary because the problem description assumes that the water has had enough time to equilibrate with the quartz mineral at 300C, and in this stage the quartz mineral is not governed by a kinetic rate.

The system described in the input file also includes very small amounts of Na and Cl. These do not impact the results but are necessary because the geochemistry module requires a charge-balance species to be defined.

# Finds the equilibrium free-concentration of SiO2(aq) in contact with Quartz at 300degC
[TimeIndependentReactionSolver<<<{"href": "../../../syntax/TimeIndependentReactionSolver/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_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"}>>> = "Quartz"
  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"}>>> = "SiO2(aq)"
  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 Na+ Cl- Quartz"
# the amount of free quartz is irrelevant to the equilibrium system, provided it is large enough, but here 400g is used to make the connection with quartz_deposition.i
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              1E-10            1E-10            0.4"
  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_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               moles            moles            kg"
  temperature<<<{"description": "The temperature (degC) of the aqueous solution"}>>> = 300.0
  add_aux_pH<<<{"description": "Add AuxVariable, called pH, that records the pH of solvent water"}>>> = false # there is no H+ in this system
  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
[]

[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 SiO2(aq) Na+ Cl-"
    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"}>>> = "Quartz"
  []
[]

The output of this simulation is that the molality of SiO2(aq) is 0.009723mol.

MOOSE input file: stage 2

The second stage uses this molality and performs the time-dependent simulation, as the temperature is reduced. The GeochemistryKineticRate is defined:

  [rate_quartz]
    type = GeochemistryKineticRate
    kinetic_species_name = Quartz
    intrinsic_rate_constant = 7.4112E2 # 2.35E-5mol/s/cm^2 = 7.411E2mol/yr/cm^2
    multiply_by_mass = true
    area_quantity = 1
    activation_energy = 72800.0
  []

The TimeDependentReactionSolver defines the free molality of SiO2(aq) at the initial time, the initial mole number of quartz and that the temperature is controlled using the temp_controller AuxVariable:

[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              Na+              Cl-              SiO2(aq)"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              1E-10            1E-10            0.009722905"
  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."}>>> = 300.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"}>>> = temp_controller
  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"}>>> = 400
  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
  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
  add_aux_pH<<<{"description": "Add AuxVariable, called pH, that records the pH of solvent water"}>>> = false # there is no H+ in this system
  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
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output used in this example
[]

The temperature controller is:

  [temp_controller_auxk]
    type = FunctionAux
    function = '300 - 275 * t'
    variable = temp_controller
    execute_on = 'timestep_begin'
  []

with time defined through:

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  dt = 0.02
  end_time = 1 # measured in years
[]

The figures below were generated using a time-step of 0.001yr. A set of AuxVariables and Postprocessors define the desired output using the mg_per_kg_SiO2(aq) variable automatically included by the TimeDependentReactionSolver:

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [temp_controller]
  []
  [diss_rate]
  []
[]
[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [temp_controller_auxk]
    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"}>>>
    function<<<{"description": "The function to use as the value"}>>> = '300 - 275 * t'
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = temp_controller
    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'
  []
  [diss_rate]
    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"}>>> = mol_change_Quartz
    expression<<<{"description": "Parsed function expression to compute"}>>> = '-mol_change_Quartz / 0.02' # 0.02 = timestep size
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = diss_rate
  []
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [mg_per_kg_sio2]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = "mg_per_kg_SiO2(aq)"
  []
  [rate_mole_per_year]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = diss_rate
  []
  [temperature]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = "solution_temperature"
  []
[]
[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 file is

# React script that is equivalent to quartz_deposition.i
time end = 1 year
T initial = 300, final = 25

swap Quartz for SiO2(aq)
400 free grams Quartz
kinetic Quartz surface = 1
kinetic Quartz pre-exp = 2.35e-5, act_eng = 72800

suppress Tridymite Chalcedony Cristobalite Amrph^silica
delxi = 0.001
go

Results

Bethke (2007) presents results in Figures 26.3 and 26.4, which look like:

Figure 1: Change in free mass of SiO2(aq) as fluid flows through a fracture, changing temperature as it does so. Compare with Bethke's Figure 26.3

Figure 2: Quartz reaction rate as fluid flows through a fracture, changing temperature as it does so. Compare with Bethke's Figure 26.4

The accuracy of the geochemistry simulation depends on the time-step size. The above figures were generated using a step size of 0.001yr.

References