Dumping minerals and then adding chemicals
Chapter 13.3 of Bethke (2007) describes a "dump" process, whereby a mineral is added to the initial solution, equilibrium sought, then any undissolved mineral is removed. In the code, this is achieved by:
subtracting from ;
then setting ;
setting the mole numbers of any surface components to zero, , as well as the molalities of unoccupied surface sites, and adsorbed species, ;
finally, swapping an appropriate aqueous species into the basis in place of .
After this process, other chemicals can be added to this system, and/or the temperature varied, and/or the pH varied, etc.
Section 15.2 of Bethke (2007) provides an example of this "dump" process. Assume:
the mineral calcite is used in the basis instead of HCO, and that the free mole number of the calcite is 0.2708 (corresponding to a free volume of calcite is cm);
the pH is initially 8;
the concentration of Na is 100mmolal;
the concentration of Ca is 10mmolal;
charge balance is enforced on Cl.
MOOSE input file
The GeochemicalModelDefinition defines the basis species, the equilibrium mineral and gas of interest in this problem, as well as using the piecewise_linear_interpolation flag for accurate comparisons with the Geochemists Workbench software:
[UserObjects]
[definition]
type = GeochemicalModelDefinition
database_file = "../../../database/moose_geochemdb.json"
basis_species = "H2O H+ Na+ Cl- Ca++ HCO3-"
equilibrium_minerals = "Calcite"
equilibrium_gases = "CO2(g)"
piecewise_linear_interpolation = true # for comparison with GWB
[]
[]
The TimeDependentReactionSolver defines:
the swaps;
the initial bulk mole number of the aqueous species in the current basis, the free mole number of Calcite and the pH (via the H
activity);that the system is closed at (the default), so that after this time no more Calcite can be added to the system to maintain its free mole number;
that at the activity constraint is removed on H so that the pH can vary;
that HCl is added at a rate of 0.001mol/s;
using
mode = 1, that Calcite is dumped from the system after each time-step (calcite never precipitates in this example, so the dumping actually only happens at the first time-step).
[TimeDependentReactionSolver]
model_definition = definition
geochemistry_reactor_name = reactor
swap_out_of_basis = "HCO3-"
swap_into_basis = "Calcite"
charge_balance_species = "Cl-"
constraint_species = "H2O Calcite Ca++ Na+ Cl- H+"
constraint_value = " 1.0 10 0.01 0.1 0.11 -8"
constraint_meaning = "kg_solvent_water free_mineral bulk_composition bulk_composition bulk_composition log10activity"
constraint_unit = " kg cm3 moles moles moles dimensionless"
ramp_max_ionic_strength_initial = 10
remove_fixed_activity_name = 'H+'
remove_fixed_activity_time = 0
source_species_names = 'HCl'
source_species_rates = 1E-3
mode = 1 # in this case, Calcite never re-precipitates, so never need to turn the dump option off
stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
execute_console_output_on = '' # only CSV output for this test
[]
The Executioner provides meaning to time, in particular that 100s is used, so a total of 0.1mol of HCl is added:
[Executioner]
type = Transient
dt = 10
end_time = 100
[]
The results presented below were run with a time-step of 1 to improve the appearance of the graphs.
A set of Postprocessors record the desired information using AuxVariables automatically added by the TimeDependentReactionSolver:
[Postprocessors]
[cm3_Calcite]
type = PointValue
point = '0 0 0'
variable = 'free_cm3_Calcite'
[]
[pH]
type = PointValue
point = '0 0 0'
variable = 'pH'
[]
[molal_CO2aq]
type = PointValue
point = '0 0 0'
variable = 'molal_CO2(aq)'
[]
[molal_CaCl+]
type = PointValue
point = '0 0 0'
variable = 'molal_CaCl+'
[]
[molal_HCO3-]
type = PointValue
point = '0 0 0'
variable = 'molal_HCO3-'
[]
[molal_Ca++]
type = PointValue
point = '0 0 0'
variable = 'molal_Ca++'
[]
[fugacity_CO2]
type = PointValue
point = '0 0 0'
variable = 'activity_CO2(g)'
[]
[]
GWB input file
The equivalent Geochemists Workbench file is
# React script that is equivalent to calcite_dumping.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O = 1 free kg
H+ = 8 pH
Na+ = 0.1 mol
Ca++ = 0.01024844 mol
Cl- = 0.12 mol
balance on Cl-
swap Calcite for HCO3-
Calcite = 0.2708 free mol
suppress all
unsuppress Calcite
printout species = long
epsilon = 1e-13
dump
react 0.1 mol of HCl
go
Note the slightly increased Ca bulk composition: this is explained here.
Results
Two cases are run:
The system is brought to equilibrium and then the calcite is "dumped" using
mode = 1in the input file. After this, 100mmol of HCl is added to the system.The system is brought to equilibrium and but the calcite is not "dumped" using
mode = 0in the input file. After this, 100mmol of HCl is added to the system.
Bethke (2007) presents the results in Figures 15.6 and 15.7. Both GWB and the geochemistry module produce the same results, as shown below.
Figure 1: CO2(g) fugacity as HCl is added to the fluid. Compare with Bethke's Figure 15.6
Figure 2: pH as HCl is added to the fluid. Compare with Bethke's Figure 15.6
Figure 3: Species concentrations as HCl is added to the fluid. Compare with Bethke's Figure 15.7
References
- Craig M. Bethke.
Geochemical and Biogeochemical Reaction Modeling.
Cambridge University Press, 2 edition, 2007.
doi:10.1017/CBO9780511619670.[BibTeX]