Cooling a solution in contact with feldspars
This page follows section 14.1 of Bethke (2007) and involves computing the equilibrium system as a solution in contact with feldspars is cooled.
The initial composition of the water is shown in Table 1. In addition, the "albite low" mineral is supressed (its precipitation is disallowed). The system is initialised at 300C and slowly cooled to 25C.
Table 1: Initial composition
| Species | Concentration |
|---|---|
| Na | 1 molal (bulk) |
| Cl | 1 molal (bulk) |
| Albite (in place of Al) | 20 cm (free) |
| Maximum microcline (in place of K) | 10 cm (free) |
| Muscovite (in place of H) | 5 cm (free) |
| Quartz (in place of SiO(aq) | 2 cm (free) |
Bethke (2007) predicts that microcline slowly precipitates, albite slowly dissolves, while muscovite and quartz remain largely unchanged (see Bethke (2007) Figure 14.1)
MOOSE input file: model definition
The MOOSE input file defines the basis species and equilibrium minerals via the GeochemicalModelDefinition. Note the appearance of the remove_all_extrapolated_secondary_species flag, which is necessary to ensure convergence of this problem. The offending secondary species is Al13O4(OH)24(7+) that has defined only at low temperatures in the original database. The default process to use for such species is to extrapolate the equilibrium constants to higher temperatures. In this case at C which prevents geochemistry from converging, so the species is removed (it is also removed by the Geochemists Workbench software).
[UserObjects]
[definition]
type = GeochemicalModelDefinition
database_file = "../../../database/moose_geochemdb.json"
basis_species = "H2O H+ Na+ Cl- Al+++ K+ SiO2(aq)"
equilibrium_minerals = "Albite Maximum Muscovite Quartz"
remove_all_extrapolated_secondary_species = true # this removes Al13O4(OH)24(7+) that has extreme logK values
[]
[]
MOOSE input file: model time-dependence
The TimeDependentReactionSolver defines:
the swaps mentioned above
the charge-balance species
the constraints on the species, which in this case hold just at the initial time because the system is closed after that time (no more minerals are added after the initialization, but they could be added during the initialization to ensure the correct number of free moles)
the initial temperature
the temperature as a function of time
other flags that allow accurate comparison with the Geochemists Workbench software
[TimeDependentReactionSolver]
model_definition = definition
geochemistry_reactor_name = reactor
swap_out_of_basis = "Al+++ K+ H+ SiO2(aq)"
swap_into_basis = "Albite Maximum Muscovite Quartz"
charge_balance_species = "Cl-"
constraint_species = "H2O Muscovite Na+ Cl- Albite Maximum Quartz"
constraint_value = " 1.0 5 1.14093 1.14093 20 10 2"
constraint_meaning = "kg_solvent_water free_mineral bulk_composition bulk_composition free_mineral free_mineral free_mineral"
constraint_unit = "kg cm3 moles moles cm3 cm3 cm3"
initial_temperature = 300
temperature = temperature
ramp_max_ionic_strength_initial = 0 # not needed in this simple problem
stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
abs_tol = 1E-14
execute_console_output_on = '' # only CSV output for this example
[]
The total bulk mole number of Na+ and Cl- is set at 1.14093 instead of 1.0 as set above. This is to compare with the GWB result (see below).
The time-dependent temperature is defined through an AuxVariable:
[AuxVariables]
[temperature]
[]
[]
[AuxKernels]
[temperature]
type = FunctionAux
variable = temperature
function = '300 - t'
execute_on = 'timestep_begin' # so that it is correct when we solve the system
[]
[]
along with:
[Executioner]
type = Transient
start_time = -10 # so that the output at 300degC is easily captured
dt = 10
end_time = 275
[]
MOOSE input file: recording the desired quantities
The desired mineral volumes and solution temperature are recorded into Postprocessors using the AuxVariables automatically included by the TimeDependentReactionSolver:
[Postprocessors]
[solution_temperature]
type = PointValue
point = '0 0 0'
variable = 'temperature'
[]
[cm3_Max_Micro]
type = PointValue
point = '0 0 0'
variable = 'free_cm3_Maximum'
[]
[cm3_Albite]
type = PointValue
point = '0 0 0'
variable = 'free_cm3_Albite'
[]
[cm3_Muscovite]
type = PointValue
point = '0 0 0'
variable = 'free_cm3_Muscovite'
[]
[cm3_Quartz]
type = PointValue
point = '0 0 0'
variable = 'free_cm3_Quartz'
[]
[]
GWB input file
The equivalent Geochemists Workbench input file is
# React script that is equivalent to cooling.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature initial = 300 C, final = 25 C
H2O = 1 free kg
swap Albite for Al+++
Albite = .19986 free mol
swap "Maximum Microcline" for K+
"Maximum Microcline" = .09196 free mol
swap Muscovite for H+
Muscovite = .03553 free mol
swap Quartz for SiO2(aq)
Quartz = .08815 free mol
Na+ = 1 molal
Cl- = 1 molal
balance on Cl-
suppress "Albite low"
printout species = long
epsilon = 1e-14
Although the bulk composition for Na+ and Cl- are specified as 1 molal, GWB changes the mole number of Cl- to 1.14 to ensure charge balance. That is the reason why the bulk composition is set at 1.14 in the MOOSE input file. In this case the results are virtually unchanged.
Results
Figure 1: Precipitated volumes as a function of temperature
References
- Craig M. Bethke.
Geochemical and Biogeochemical Reaction Modeling.
Cambridge University Press, 2 edition, 2007.
doi:10.1017/CBO9780511619670.[BibTeX]