Solubility of gypsum and associated activities
This example closely follows Section 8.3 of Bethke (2007).
The solubility of gypsum (CaSO.2HO) as a function of NaCl concentration is explored. A simulation is initialized with some undissolved Gypsum and small quantities of Na and Cl, and then NaCl is progressively added. The following assumptions are made:
Gypsum is used in the basis instead of the aqueous species Ca.
The initial amount of free gypsum is 0.5814mol (g).
Charge balance is performed on SO.
Hence, the basis is (H20, Na+, Cl-, SO4–, gypsum). The simulation computes how much gypsum is dissolved as a function of the concentration of NaCl.
The results are shown in Figure 8.6 of Bethke (2007).
The MOOSE input file
The input file begins by defining the basis aqueous species and the equilibrium minerals (only Gypsum) for the model. The piecewise_linear_interpolation flag is used for accurate comparison with the Geochemists Workbench software.
[UserObjects]
[definition]
type = GeochemicalModelDefinition
database_file = "../../../database/moose_geochemdb.json"
basis_species = "H2O Cl- Na+ SO4-- Ca++"
equilibrium_minerals = "Gypsum"
piecewise_linear_interpolation = true # for comparison with GWB
[]
[]
A TimeDependentReactionSolver is used to simulate the adding of NaCl:
Ca is swapped out of the basis in favor of Gypsum.
Na and Cl are provided with a small initial molality. The molality for SO will actually be controlled by charge neutrality.
0.5814 free moles of Gypsum are introduced into the initial system
During the initial setup, the
geochemistrymodule will find the equilibrium configuration corresponding to these initial conditions. By default, it will then close the system (theclose_system_at_timeinput defaults to 0.0). This stops any further Gypsum from entering the system.The rate of addition of NaCl is defined to be 1.0mol/s
The other flags enable an accurate comparison with the Geochemists Workbench software.
[TimeDependentReactionSolver]
model_definition = definition
swap_out_of_basis = "Ca++"
swap_into_basis = "Gypsum"
charge_balance_species = "SO4--"
constraint_species = "H2O Cl- Na+ SO4-- Gypsum"
constraint_value = " 1.0 1E-10 1E-10 1E-6 0.5814"
constraint_meaning = "kg_solvent_water free_concentration free_concentration bulk_composition free_mineral"
constraint_unit = " kg molal molal moles moles"
source_species_names = 'NaCl'
source_species_rates = '1.0'
add_aux_pH = false # there is no H+ in the problem
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-12
execute_console_output_on = '' # only CSV output in this example
[]
The time-stepping and output are defined in the usual MOOSE way. In this case, a FunctionDT timestepper is used to capture the nonlinear behaviour at the beginning of the simulation:
[Executioner]
type = Transient
[TimeStepper]
type = FunctionDT
function = timestepper
[]
end_time = 3
[]
[Outputs]
csv = true
[]
An AuxVariable captures the amount of dissolved gypsum by using the bulk_moles_Gypsum and free_mg_Gypsum AuxVariables that are automatically added by the TimeDependentReactionSolver is used to simulate the adding of NaCl:
[AuxVariables]
[dissolved_gypsum_moles]
[]
[]
[AuxKernels]
[dissolved_gypsum_moles]
type = ParsedAux
coupled_variables = 'bulk_moles_Gypsum free_mg_Gypsum'
expression = 'bulk_moles_Gypsum - free_mg_Gypsum / 1000 / 172.168 '
variable = dissolved_gypsum_moles
execute_on = 'timestep_end'
[]
[]
Finally, Postprocessors allow the desired information to be written to the CSV output file
[Postprocessors]
[cl_molal]
type = PointValue
point = '0 0 0'
variable = 'molal_Cl-'
[]
[dissolved_gypsum_mol]
type = PointValue
point = '0 0 0'
variable = dissolved_gypsum_moles
[]
[]
GWB input file
An equivalent Geochemists Workbench input file is
# React script that is equivalent to gypsum_solubility.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O = 1 free kg
Cl- = 1e-10 free mol
Na+ = 1e-10 free mol
SO4-- = 1e-6 mol
balance on SO4--
swap Gypsum for Ca++
Gypsum = .5814 free mol
react 3 mol of NaCl
suppress ALL
unsuppress Gypsum
dxprint = .01
A total of 3 moles of NaCl is added to the system, which is equivalent to the 1mol/s over 3s used by the MOOSE input file. To extract the results from Geochemists Workbench its inbuilt graphing capability may be utilized, or the Cl molality vs the moles of Ca in the solution may be extracted from plaintext output file.
Results
The results are shown in Figure 1.
Figure 1: Gypsum solubility as a function of chlorine molality.
References
- Craig M. Bethke.
Geochemical and Biogeochemical Reaction Modeling.
Cambridge University Press, 2 edition, 2007.
doi:10.1017/CBO9780511619670.[BibTeX]