Solubility of gypsum and associated activities

This example closely follows Section 8.3 of Bethke (2007).

The solubility of gypsum (CaSO $var element = document.getElementById("moose-equation-e2e4a6f4-3e40-4cb5-a3a3-289a0b8f0cf4");katex.render("_{4}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .2H $var element = document.getElementById("moose-equation-a31beffd-668e-4826-97a3-822cc38e5bfd");katex.render("_{2}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ O) as a function of NaCl concentration is explored. A simulation is initialized with some undissolved Gypsum and small quantities of Na $var element = document.getElementById("moose-equation-2ecaa28f-db70-4e7a-9489-89f27119e048");katex.render("^{+}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and Cl $var element = document.getElementById("moose-equation-f3477090-6b89-4c5a-85e0-7f32875dbc71");katex.render("^{-}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , and then NaCl is progressively added. The following assumptions are made:

Gypsum is used in the basis instead of the aqueous species Ca $var element = document.getElementById("moose-equation-14fa0dd8-30a2-479b-acca-0f907d395d3a");katex.render("^{2+}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .
The initial amount of free gypsum is 0.5814 $var element = document.getElementById("moose-equation-1da01ad7-e62a-42b2-860e-3f66c5322784");katex.render("\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ mol ( $var element = document.getElementById("moose-equation-6c55ff4a-2c87-4ae9-ac77-e4519fee93d8");katex.render("\\approx 100\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ g).
Charge balance is performed on SO $var element = document.getElementById("moose-equation-6f2d3cbf-7b08-47bf-aa13-6ce5d71988c4");katex.render("_{4}^{2-}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .

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<<<{"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 Cl- Na+ SO4-- Ca++"
    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"}>>> = "Gypsum"
    piecewise_linear_interpolation<<<{"description": "If true then use a piecewise-linear interpolation of logK and Debye-Huckel parameters, regardless of the interpolation type specified in the database file.  This can be useful for comparing with results using other geochemistry software"}>>> = true # for comparison with GWB
  []
[]

A TimeDependentReactionSolver is used to simulate the adding of NaCl:

Ca $var element = document.getElementById("moose-equation-01d8f6fc-f8fb-473d-9ec8-e2d7c17a2dbe");katex.render("^{2+}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is swapped out of the basis in favor of Gypsum.
Na $var element = document.getElementById("moose-equation-a6744e44-ff49-4c65-97f3-364dfc3a290a");katex.render("^{+}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and Cl $var element = document.getElementById("moose-equation-04c16598-d272-412f-8779-66915976dd70");katex.render("^{+}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ are provided with a small initial molality. The molality for SO $var element = document.getElementById("moose-equation-016cf31c-ffa7-40ed-8bc4-116a84d459a8");katex.render("_{4}^{2-}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ will actually be controlled by charge neutrality.
0.5814 free moles of Gypsum are introduced into the initial system
During the initial setup, the geochemistry module will find the equilibrium configuration corresponding to these initial conditions. By default, it will then close the system (the close_system_at_time input defaults to 0.0). This stops any further Gypsum from entering the system.
The rate of addition of NaCl is defined to be 1.0 $var element = document.getElementById("moose-equation-aeb833dc-526d-4374-a52d-5f991440338c");katex.render("\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ mol/s
The other flags enable an accurate comparison with the Geochemists Workbench software.

[TimeDependentReactionSolver<<<{"href": "../../../syntax/TimeDependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  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"}>>> = "Ca++"
  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"}>>> = "Gypsum"
  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."}>>> = "SO4--"
  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              Cl-                Na+                SO4--            Gypsum"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              1E-10              1E-10              1E-6             0.5814"
  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 free_concentration free_concentration bulk_composition 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               molal              molal              moles            moles"
  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."}>>> = 'NaCl'
  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!"}>>> = '1.0'
  add_aux_pH<<<{"description": "Add AuxVariable, called pH, that records the pH of solvent water"}>>> = false # there is no H+ in the problem
  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 # not needed in this simple problem
  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
  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
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # 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<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  [TimeStepper<<<{"href": "../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = 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<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [dissolved_gypsum_moles]
  []
[]

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [dissolved_gypsum_moles]
    type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../../source/auxkernels/ParsedAux.html"}>>>
    coupled_variables<<<{"description": "Vector of coupled variable names"}>>> = 'bulk_moles_Gypsum free_mg_Gypsum'
    expression<<<{"description": "Parsed function expression to compute"}>>> = 'bulk_moles_Gypsum - free_mg_Gypsum / 1000 / 172.168 '
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = dissolved_gypsum_moles
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
  []
[]

Finally, Postprocessors allow the desired information to be written to the CSV output file

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [cl_molal]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'molal_Cl-'
  []
  [dissolved_gypsum_mol]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 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 1 $var element = document.getElementById("moose-equation-74d5d11d-8415-4486-8c6e-7b03c4551561");katex.render("\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ mol/s over 3 $var element = document.getElementById("moose-equation-4d526675-8a35-4180-b0d1-d8e6fc5fbf8a");katex.render("\\,", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ s used by the MOOSE input file. To extract the results from Geochemists Workbench its inbuilt graphing capability may be utilized, or the Cl $var element = document.getElementById("moose-equation-c04ce604-a3c1-41c9-b7b5-05a92bac76af");katex.render("^{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ molality vs the moles of Ca $var element = document.getElementById("moose-equation-b3b3c7a1-069f-4dea-b20c-390e04a32596");katex.render("^{2+}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ 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.

@book{bethke_2007,
    author = "Bethke, Craig M.",
    place = "Cambridge",
    edition = "2",
    title = "Geochemical and Biogeochemical Reaction Modeling",
    DOI = "10.1017/CBO9780511619670",
    publisher = "Cambridge University Press",
    year = "2007"
}

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.i)

[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
[]

[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
  []
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0    0.1'
    y = '0.01 0.1'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs]
  csv = true
[]

[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'
  []
[]

[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
  []
[]

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.i)

[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
[]

[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
  []
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0    0.1'
    y = '0.01 0.1'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs]
  csv = true
[]

[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'
  []
[]

[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
  []
[]

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.i)

[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
[]

[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
  []
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0    0.1'
    y = '0.01 0.1'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs]
  csv = true
[]

[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'
  []
[]

[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
  []
[]

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.i)

[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
[]

[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
  []
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0    0.1'
    y = '0.01 0.1'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs]
  csv = true
[]

[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'
  []
[]

[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
  []
[]

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.i)

[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
[]

[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
  []
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0    0.1'
    y = '0.01 0.1'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs]
  csv = true
[]

[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'
  []
[]

[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
  []
[]

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.i)

[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
[]

[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
  []
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0    0.1'
    y = '0.01 0.1'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs]
  csv = true
[]

[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'
  []
[]

[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
  []
[]

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.i)

[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
[]

[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
  []
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0    0.1'
    y = '0.01 0.1'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 3
[]

[Outputs]
  csv = true
[]

[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'
  []
[]

[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
  []
[]

(modules/geochemistry/test/tests/solubilities_and_activities/gypsum_solubility.rea)

# 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

The MOOSE input file
GWB input file
Results
References