Quartz deposition in a fracture

Section 26.2 of Bethke (2007) describes quartz deposition in a hydrothermal fracture. The setup is similar to quartz dissolution but with a different form of the reaction rate. The reaction is QuartzSiO2(aq) ,\mathrm{Quartz} \rightarrow \mathrm{SiO}_{2}\mathrm{(aq)} \ , with rate r=AseE/(RT)(1QK) .r = A_{s} e^{-E/(RT)}\left( 1 - \frac{Q}{K} \right) \ . It is assumed that:

  • there is 1\,kg of water;

  • the initial temperature is 300^{\circ}C and this steadily reduces to 25^{\circ}C over the course of 1 year;

  • the mineral quartz is used instead of SiO2_{2}(aq) in the basis initially while quartz is an equilibrium mineral (before it is kinetically-controlled);

  • initially 400\,g (6.657313\,mol) of quartz is added to the water;

  • the specific surface area is As=2.35×105A_{s} = 2.35\times 10^{-5}\,cm2^{2}/g(quartz);

  • the activation energy is E=72800E=72800\,J.mol1^{-1};

  • all other silica-containing minerals are prevented from precipitating.

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 300^{\circ}C, 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]
  model_definition = definition
  geochemistry_reactor_name = reactor
  swap_into_basis = "Quartz"
  swap_out_of_basis = "SiO2(aq)"
  charge_balance_species = "Cl-"
  constraint_species = "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 = "  1.0              1E-10            1E-10            0.4"
  constraint_meaning = "kg_solvent_water bulk_composition bulk_composition free_mineral"
  constraint_unit = "   kg               moles            moles            kg"
  temperature = 300.0
  add_aux_pH = false # there is no H+ in this system
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
[]

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O SiO2(aq) Na+ Cl-"
    equilibrium_minerals = "Quartz"
  []
[]
(moose/modules/geochemistry/test/tests/kinetics/quartz_equilibrium_at300degC.i)

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

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
[]
(moose/modules/geochemistry/test/tests/kinetics/quartz_deposition.i)

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]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = "Cl-"
  constraint_species = "H2O              Na+              Cl-              SiO2(aq)"
  constraint_value = "  1.0              1E-10            1E-10            0.009722905"
  constraint_meaning = "kg_solvent_water bulk_composition bulk_composition free_concentration"
  constraint_unit = "   kg               moles            moles            molal"
  initial_temperature = 300.0
  temperature = temp_controller
  kinetic_species_name = Quartz
  kinetic_species_initial_value = 400
  kinetic_species_unit = g
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  add_aux_pH = false # there is no H+ in this system
  evaluate_kinetic_rates_always = true # implicit time-marching used for stability
  execute_console_output_on = '' # only CSV output used in this example
[]
(moose/modules/geochemistry/test/tests/kinetics/quartz_deposition.i)

The temperature controller is:

[temp_controller_auxk]
  type = FunctionAux
  function = '300 - 275 * t'
  variable = temp_controller
  execute_on = 'timestep_begin'
[]
(moose/modules/geochemistry/test/tests/kinetics/quartz_deposition.i)

with time defined through:

[Executioner]
  type = Transient
  dt = 0.02
  end_time = 1 # measured in years
[]
(moose/modules/geochemistry/test/tests/kinetics/quartz_deposition.i)

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

[AuxVariables]
  [temp_controller]
  []
  [diss_rate]
  []
[]
[AuxKernels]
  [temp_controller_auxk]
    type = FunctionAux
    function = '300 - 275 * t'
    variable = temp_controller
    execute_on = 'timestep_begin'
  []
  [diss_rate]
    type = ParsedAux
    coupled_variables = mol_change_Quartz
    expression = '-mol_change_Quartz / 0.02' # 0.02 = timestep size
    variable = diss_rate
  []
[]
[Postprocessors]
  [mg_per_kg_sio2]
    type = PointValue
    variable = "mg_per_kg_SiO2(aq)"
  []
  [rate_mole_per_year]
    type = PointValue
    variable = diss_rate
  []
  [temperature]
    type = PointValue
    variable = "solution_temperature"
  []
[]
[Outputs]
  csv = true
[]
(moose/modules/geochemistry/test/tests/kinetics/quartz_deposition.i)

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
(moose/modules/geochemistry/test/tests/kinetics/quartz_deposition.rea)

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.001\,yr.

References

  1. Craig M. Bethke. Geochemical and Biogeochemical Reaction Modeling. Cambridge University Press, 2 edition, 2007. doi:10.1017/CBO9780511619670.[BibTeX]