Flushing minerals

MOOSE: stage 1

This is run in MOOSE using a 2-stage approach. The first stage finds the free molality of SiO2(aq) in the equilibrium configuration before the system is flushed. It is necessary to find this because since Quartz is a kinetic mineral it cannot be part of the basis in the second stage, so SiO2(aq) takes its place in the basis in that stage. The first stage is a time-independent simulation:

# Alkali flushing of a reservoir (an example of flushing)
# This input file finds the molality of SiO2(aq), which is needed because when Quartz is specified as kinetic it cannot appear in the basis
[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Cl- Na+ Ca++ HCO3- Mg++ K+ Al+++ SiO2(aq)"
    equilibrium_minerals = "Calcite Dolomite-ord Muscovite Kaolinite Quartz"
  []
[]

[TimeIndependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = "Cl-"
  swap_into_basis = "Calcite Dolomite-ord Muscovite Kaolinite Quartz"
  swap_out_of_basis = "HCO3- Mg++ K+ Al+++ SiO2(aq)"
  constraint_species = "H2O              H+            Cl-              Na+              Ca++             Calcite      Dolomite-ord Muscovite    Kaolinite    Quartz"
  constraint_value = "  1.0              -5            1.0              1.0              0.2              9.88249      3.652471265  1.2792268    1.2057878    226.99243"
  constraint_meaning = "kg_solvent_water log10activity bulk_composition bulk_composition bulk_composition free_mineral free_mineral free_mineral free_mineral free_mineral"
  constraint_unit = "   kg               dimensionless moles            moles            moles            moles        moles        moles        moles        moles"
  temperature = 70.0
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
  execute_console_output_on = ''
  abs_tol = 1E-14
[]

[GlobalParams]
  point = '0 0 0'
[]
[Postprocessors]
  [pH]
    type = PointValue
    variable = "pH"
  []
  [molality_sio2]
    type = PointValue
    variable = 'molal_SiO2(aq)'
  []
[]

[Outputs]
  csv = true
[]

MOOSE: stage 2

This stage adds the NaOH solution and uses a kinetically-controlled quartz. The GeochemicalModelDefinition defines the basis species and equilibrium minerals. It also specifies that Quartz is a kinetic mineral:

  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Cl- Na+ Ca++ HCO3- Mg++ K+ Al+++ SiO2(aq)"
    equilibrium_minerals = "Analcime Calcite Dawsonite Dolomite-ord Gibbsite Kaolinite Muscovite Paragonite Phlogopite"
    kinetic_minerals = "Quartz"
    kinetic_rate_descriptions = "rate_quartz"

The kinetic rate for Quartz is specified by a GeochemistryKineticRate UserObject

  [rate_quartz]
    type = GeochemistryKineticRate
    kinetic_species_name = Quartz
    intrinsic_rate_constant = 1.5552E-13 # 1.8E-18mol/s/cm^2 = 1.5552E-13mol/day/cm^2
    multiply_by_mass = true
    area_quantity = 1000
    promoting_species_names = "H+"
    promoting_species_indices = "-0.5"
  []

The TimeDependentReactionSolver specifies:

• the swaps;

• the initial pH (via the H activity), the bulk composition of the aqueous species (note the difference for Ca++ compared with GWB) and the free mole number for the minerals (not Quartz because it is a kinetic species) and free molality of SiO2(aq) from stage 1;

• the temperature;

• the initial number of moles for Quartz

• that the kinetic rate should be updated during the Newton process that finds the equilibrium configuration at each time-step (which is unnecessary and computationally inefficient for this example, but is included for illustration)

• that the system is closed at , which is the default, which means that no more SiO2(aq) or HO can be added/removed from the system to maintain the free mole numbers of these species;

• that the fixed pH is removed at ;

• that "flush" mode is active;

• the reactant rate: HO is added at rate 27.72mol/day, and NaOH is added at rate 0.25mol/day.

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = "Cl-"
  swap_into_basis = "Calcite Dolomite-ord Muscovite Kaolinite"
  swap_out_of_basis = "HCO3- Mg++ K+ Al+++"
  constraint_species = "H2O              H+            Cl-              Na+              Ca++             Calcite      Dolomite-ord Muscovite    Kaolinite    SiO2(aq)"
  constraint_value = "  1.0              -5            1.0              1.0              0.2              9.88249      3.652471265  1.2792268    1.2057878    0.000301950628974"
  constraint_meaning = "kg_solvent_water log10activity bulk_composition bulk_composition bulk_composition free_mineral free_mineral free_mineral free_mineral free_concentration"
  constraint_unit = "   kg               dimensionless moles            moles            moles            moles        moles        moles        moles        molal"
  initial_temperature = 70.0
  temperature = 70.0
  kinetic_species_name = Quartz
  kinetic_species_initial_value = 226.992243
  kinetic_species_unit = moles
  evaluate_kinetic_rates_always = true # implicit time-marching used for stability
  close_system_at_time = 0.0
  remove_fixed_activity_name = "H+"
  remove_fixed_activity_time = 0.0
  mode = 3 # flush through the NaOH solution specified below:
  source_species_names = "H2O   Na+  OH-"
  source_species_rates = "27.75 0.25 0.25" # 1kg water/2days = 27.75moles/day.  0.5mol Na+/2days = 0.25mol/day
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  stoichiometric_ionic_str_using_Cl_only = true
  execute_console_output_on = '' # only CSV output for this test
  abs_tol = 1E-13
[]

The Executioner defines the time-stepping and provides physical meaning to "time"

[Executioner]
  type = Transient
  dt = 1
  end_time = 20 # measured in days
[]

The graphs below were generated using dt = 0.2 but the file in the test suite uses a larger time-step for efficiency.

A set of Postprocessors record the desired information using the AuxVariables automatically added by the TimeDependentReactionSolver. For example:

  [pH]
    type = PointValue
    variable = "pH"
  []

GWB input file

The equivalent Geochemists Workbench input file is

# React script that is equivalent to flushing.i
T = 70
H+ = 5 pH
Na+ = 1 molal
Ca++ = 0.288445 molal
Cl- = 1 molal
balance on Cl-
swap Dolomite-ord for Mg++
swap Muscovite for K+
swap Kaolinite for Al+++
swap Quartz for SiO2(aq)
swap Calcite for HCO3-
Calcite = 9.88249 free mol
Dolomite-ord = 3.652471265 free mol
Muscovite = 1.2792268 free mol
Kaolinite = 1.2057878 free mol
Quartz = 226.99243 free mol

react 1 kg H2O
react 0.5 moles Na+
react 0.5 moles OH-

suppress all
unsuppress Analcime Calcite Dawsonite Dolomite-ord Gibbsite Kaolinite Muscovite Paragonite Phlogopite Quartz Saponite-Ca
kinetic Quartz rate_con = 1.8e-18, apower(H+) = -0.5 surface = 1000

reactants times 10
time end = 20 days
flush
go

Note the different initial molality specified for Ca. In this case it makes no noticable difference to the result, but users desiring an exact correspondance between the two software packages should understand the differences.

Results

Figures 30.3 and 30.4 in Bethke (2007) shows the results. The GWB and geochemistry module results are shown in the figures below.

Figure 1: pH during the alkali flooding experiment. Compare with Bethke's Figure 30.3

Figure 2: Volume change of the precipitating minerals during the alkali flooding experiment. Compare with Bethke's Figure 30.4

Figure 3: Volume change of the dissolving minerals during the alkali flooding experiment. Compare with Bethke's Figure 30.4

Figure 4: Volume change of the kinetically-controlled quartz mineral during the alkali flooding experiment. Compare with Bethke's Figure 30.4

References