Flushing minerals
Chapter 13.3 of Bethke (2007) describes a "flush" process, in which a specified solution is added to the system, and at the same time an equal volume of the equilibrium solution, consisting of solvent water, primary aqueous species, dissolved minerals and dissolved gases, is removed. In the code, this is achieved by altering , , and appropriately.
Section 30.2 of Bethke (2007) describes an example of this involving alkali flooding of a petroleum reservoir. It is assumed that quartz dissolves and precipitates with rate where
[cm] is set through a specific surface area of cm/g(quartz);
mol.cm.smol.cm.day.
is the activity of H.
The initial configuration is given by:
1 molal Na;
0.2 molal Ca;
1 molal Cl (charge-balance species);
C;
pH .
The initial minerals are:
9.88249mol (cm) of free Calcite (in place of HCO in the basis);
3.6525mol (cm) of free Dolomite-ord (in place of Mg);
1.27923mol (cm) of free Muscovite (in place of K);
1.20579mol (cm) of free Kaolinite (in place of Al);
226.99243mol (cm) of free Quartz (in place of SiO(aq)).
The simulation contains only the following minerals: Analcime, Calcite, Dawsonite, Dolomite-ord, Gibbsite, Kaolinite, Muscovite, Paragonite and Phlogopite
Over the course of 20 days, the following is flushed through the system:
10kg of HO;
5moles of Na;
5moles of OH.
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<<<{"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 H+ Cl- Na+ Ca++ HCO3- Mg++ K+ Al+++ SiO2(aq)"
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"}>>> = "Calcite Dolomite-ord Muscovite Kaolinite Quartz"
[]
[]
[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
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-"
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"}>>> = "Calcite Dolomite-ord Muscovite Kaolinite 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"}>>> = "HCO3- Mg++ K+ Al+++ SiO2(aq)"
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+ Ca++ Calcite Dolomite-ord Muscovite Kaolinite Quartz"
constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = " 1.0 -5 1.0 1.0 0.2 9.88249 3.652471265 1.2792268 1.2057878 226.99243"
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 bulk_composition free_mineral free_mineral free_mineral free_mineral 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 dimensionless moles moles moles moles moles moles moles moles"
temperature<<<{"description": "The temperature (degC) of the aqueous solution"}>>> = 70.0
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"}>>> = ''
abs_tol<<<{"description": "If the residual of the algebraic system (measured in mol) is lower than this value, the Newton process (that finds the aqueous configuration) is deemed to have converged"}>>> = 1E-14
[]
[GlobalParams<<<{"href": "../../../syntax/GlobalParams/index.html"}>>>]
point = '0 0 0'
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
[pH]
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."}>>> = "pH"
[]
[molality_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."}>>> = 'molal_SiO2(aq)'
[]
[]
[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = 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_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<<<{"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-"
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"}>>> = "Calcite Dolomite-ord Muscovite Kaolinite"
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"}>>> = "HCO3- Mg++ K+ Al+++"
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+ Ca++ Calcite Dolomite-ord Muscovite Kaolinite SiO2(aq)"
constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = " 1.0 -5 1.0 1.0 0.2 9.88249 3.652471265 1.2792268 1.2057878 0.000301950628974"
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 bulk_composition free_mineral free_mineral free_mineral free_mineral 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 moles moles moles moles moles 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"}>>> = 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"}>>> = 226.992243
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"}>>> = moles
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
close_system_at_time<<<{"description": "Time at which to 'close' the system, that is, change a kg_solvent_water constraint to moles_bulk_water, and all free_molality and free_moles_mineral_species to moles_bulk_species"}>>> = 0.0
remove_fixed_activity_name<<<{"description": "The name of the species that should have their activity or fugacity constraint removed at time given in remove_fixed_activity_time. There should be an equal number of these names as times given in remove_fixed_activity_time. Each of these must be in the basis and have an activity or fugacity constraint"}>>> = "H+"
remove_fixed_activity_time<<<{"description": "The times at which the species in remove_fixed_activity_name should have their activity or fugacity constraint removed."}>>> = 0.0
mode<<<{"description": "This may vary temporally. If mode=1 then 'dump' mode is used, which means all non-kinetic mineral masses are removed from the system before the equilibrium solution is sought (ie, removal occurs at the beginning of the time step). If mode=2 then 'flow-through' mode is used, which means all mineral masses are removed from the system after it the equilbrium solution has been found (ie, at the end of a time step). If mode=3 then 'flush' mode is used, then before the equilibrium solution is sought (ie, at the start of a time step) water+species is removed from the system at the same rate as pure water + non-mineral solutes are entering the system (specified in source_species_rates). If mode=4 then 'heat-exchanger' mode is used, which means the entire current aqueous solution is removed, then the source_species are added, then the temperature is set to 'cold_temperature', the system is solved and any precipitated minerals are removed, then the temperature is set to 'temperature', the system re-solved and any precipitated minerals are removed. If mode is any other number, no special mode is active (the system simply responds to the source_species_rates, controlled_activity_value, etc)."}>>> = 3 # flush through the NaOH solution specified below:
source_species_names<<<{"description": "The name of the species that are added at rates given in source_species_rates. There must be an equal number of these as source_species_rates."}>>> = "H2O Na+ OH-"
source_species_rates<<<{"description": "Rates, in mols/time_unit, of addition of the species with names given in source_species_names. A negative value corresponds to removing a species: be careful that you don't cause negative mass problems!"}>>> = "27.75 0.25 0.25" # 1kg water/2days = 27.75moles/day. 0.5mol Na+/2days = 0.25mol/day
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
execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output for this test
abs_tol<<<{"description": "If the residual of the algebraic system (measured in mol) is lower than this value, the Newton process (that finds the aqueous configuration) is deemed to have converged"}>>> = 1E-12
[]The Executioner defines the time-stepping and provides physical meaning to "time"
[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
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 noticeable difference to the result, but users desiring an exact correspondence 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
- Craig M. Bethke.
Geochemical and Biogeochemical Reaction Modeling.
Cambridge University Press, 2 edition, 2007.
doi:10.1017/CBO9780511619670.[BibTeX]