Install MOOSE link

New Users link

Examples and Tutorials link

Application Usage link

Physics and Syntax

Application Development link

Framework Development link

MOOSEDocs link

Infrastructure

Questions

Information and Tools

INL Applications and Remote Access link

Flow-through reactions that remove minerals

Chapter 13.3 of Bethke (2007) describes a "flow-through" process, whereby a mineral is removed from the system at each stage in a reaction process (such as progressively adding chemicals, changing the temperature, changing the pH, etc). In the code, this is achieved by the following process, after each stage's equilibrium configuration is computed:

  • subtracting from ;

  • then setting to a tiny number (but not zero, otherwise it might be swapped out of the basis);

  • setting the mole numbers of any surface components to zero, , as well as the molalities of unoccupied surface sites, and adsorbed species, .

Section 24.3 of Bethke (2007) describes an example of this involving the evaporation of seawater. Note that Bethke (2007) uses the HMW activity model, but geochemistry uses the Debye-Huckel approach, so the MOOSE results are not identical to Bethke (2007). Nevertheless, it is hoped that the following description helps users set up similar models.

The initial composition of the seawater is shown in Table 1. In addition:

  • CO(g) is used in the basis instead of H;

  • the fugacity of CO(g) is fixed to ;

  • after equilibration, all minerals are dumped from the system.

Table 1: Major element composition of seawater

SpeciesMolesApprox conc (mg.kg)
Cl0.565619350
Na0.485010760
SO0.029242710
Mg0.055011290
Ca0.01063411
K0.010576399
HCO0.002412142

Two cases are run. Each involves gradually removing 55mol of HO from the solution. The cases are:

  1. As the water is removed, minerals are allowed to precipitate, and then potentially dissolve, to maintain equilibrium.

  2. As the water is removed, minerals are allowed to precipitate, but when they precipitate they are immediately removed from the system following the "flow-through" method.

Only 6g of solvent water remains that the end of each simulation.

MOOSE input file: flow through

The GeochemicalModelDefinition defines the basis species, the equilibrium minerals and gases:

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Cl- Ca++ Mg++ Na+ K+ SO4-- HCO3-"
    equilibrium_minerals = "Dolomite Epsomite Gypsum Halite Magnesite Mirabilite Sylvite"
    equilibrium_gases = "CO2(g)"
    piecewise_linear_interpolation = true # for precise agreement with GWB
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/seawater_evaporation_flow_through.i)

The TimeDependentReactionSolver defines:

  • the swap;

  • the bulk composition of the aqueous species as well as the CO(g) fugacity via its activity;

  • that the system is closed at (the default) which, in this case, means that no more solvent water may be added to the system to maintain the 1.0kg of solvent water;

  • that the fugacity constraint remains throughout the simulation (there is are remove_activity lines);

  • that HO is removed from the system at a rate of 1.0mol/s.

  • the mode is defined via the AuxVariable called "mode"

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  swap_out_of_basis = "H+"
  swap_into_basis = "  CO2(g)"
  charge_balance_species = "Cl-" # this means the bulk moles of Cl- will not be exactly as set below
  constraint_species = "H2O              CO2(g)        Cl-              Na+              SO4--            Mg++             Ca++             K+               HCO3-"
  constraint_value = "  1.0              -3.5          0.5656           0.4850           0.02924          0.05501          0.01063          0.010576055      0.002412"
  constraint_meaning = "kg_solvent_water log10fugacity bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition"
  constraint_unit = "   kg               dimensionless moles            moles            moles            moles            moles            moles            moles"
  source_species_names = "H2O"
  source_species_rates = "-1.0" # 1kg H2O = 55.51 moles, each time step removes 1 mole
  mode = mode
  ramp_max_ionic_strength_initial = 0 # not needed in this simple example
  stoichiometric_ionic_str_using_Cl_only = true # for precise agreement with GWB
  execute_console_output_on = '' # only CSV output for this example
[]
(modules/geochemistry/test/tests/time_dependent_reactions/seawater_evaporation_flow_through.i)

The TimeStepper and Executioner define the meaning of time. In particular, the simulation is run for 55s, so that 55mol of water is removed::

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0 50 55'
    y = '5 5 1'
  []
[]
[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 55
[]

[AuxVariables]
  [mode]
  []
  [dolomite_mol]
  []
  [halite_mol]
  []
  [gypsum_mol]
  []
  [mirabilite_mol]
  []
[]

[AuxKernels]
  [mode_auxk]
    type = FunctionAux
    variable = mode
    function = 'if(t<=1.0, 1.0, 2.0)' # initial "dump" then "flow_through"
    execute_on = 'timestep_begin'
  []
  [dolomite_mol_auxk]
    type = GeochemistryQuantityAux
    reactor = reactor
    variable = dolomite_mol
    species = Dolomite
    quantity = moles_dumped
  []
  [gypsum_mol_auxk]
    type = GeochemistryQuantityAux
    reactor = reactor
    variable = gypsum_mol
    species = Gypsum
    quantity = moles_dumped
  []
  [halite_mol]
    type = GeochemistryQuantityAux
    reactor = reactor
    variable = halite_mol
    species = Halite
    quantity = moles_dumped
  []
  [mirabilite_mol]
    type = GeochemistryQuantityAux
    reactor = reactor
    variable = mirabilite_mol
    species = Mirabilite
    quantity = moles_dumped
  []
[]

[GlobalParams]
  point = '0 0 0'
[]
[Postprocessors]
  [solvent_kg]
    type = PointValue
    variable = 'kg_solvent_H2O'
  []
  [dolomite_mol]
    type = PointValue
    variable = dolomite_mol
  []
  [gypsum_mol]
    type = PointValue
    variable = 'gypsum_mol'
  []
  [halite_mol]
    type = PointValue
    variable = 'halite_mol'
  []
  [mirabilite_mol]
    type = PointValue
    variable = 'mirabilite_mol'
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/seawater_evaporation_flow_through.i)

Note that a smaller time-step size (y = '1 0.1 0.001') was used to generate the figures below, but the test-suite file uses larger time-steps for efficiency.

The mode is defined by the AuxKernel

  [mode_auxk]
    type = FunctionAux
    variable = mode
    function = 'if(t<=1.0, 1.0, 2.0)' # initial "dump" then "flow_through"
    execute_on = 'timestep_begin'
  []
(modules/geochemistry/test/tests/time_dependent_reactions/seawater_evaporation_flow_through.i)

A set of AuxVariables and AuxKernels such as

  [dolomite_mol_auxk]
    type = GeochemistryQuantityAux
    reactor = reactor
    variable = dolomite_mol
    species = Dolomite
    quantity = moles_dumped
  []
(modules/geochemistry/test/tests/time_dependent_reactions/seawater_evaporation_flow_through.i)

record the desired quantities.

MOOSE input file: no flow through

The no flow-through case is similar (just varying in the mode, the AuxKernels and the Postprocessors):

#Progressively remove H2O until virtually none remains
[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Cl- Ca++ Mg++ Na+ K+ SO4-- HCO3-"
    equilibrium_minerals = "Dolomite Epsomite Gypsum Halite Magnesite Mirabilite Sylvite"
    equilibrium_gases = "CO2(g)"
    piecewise_linear_interpolation = true # for precise agreement with GWB
  []
[]

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  swap_out_of_basis = "H+"
  swap_into_basis = "  CO2(g)"
  charge_balance_species = "Cl-" # this means the bulk moles of Cl- will not be exactly as set below
  constraint_species = "H2O              CO2(g)        Cl-              Na+              SO4--            Mg++             Ca++             K+               HCO3-"
  constraint_value = "  1.0              -3.5          0.5656           0.4850           0.02924          0.05501          0.01063          0.010576055      0.002412"
  constraint_meaning = "kg_solvent_water log10fugacity bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition"
  constraint_unit = "   kg               dimensionless moles            moles            moles            moles            moles            moles            moles"
  close_system_at_time = 0
  source_species_names = "H2O"
  source_species_rates = "-1.0" # 1kg H2O = 55.51 moles, each time step removes 1 mole
  mode = mode
  ramp_max_ionic_strength_initial = 0 # not needed in this simple example
  stoichiometric_ionic_str_using_Cl_only = true # for precise agreement with GWB
  execute_console_output_on = '' # only CSV output for this example
[]

[Functions]
  [timestepper]
    type = PiecewiseLinear
    x = '0 50 55'
    y = '5 5 1'
  []
[]

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

[AuxVariables]
  [mode]
  []
[]

[AuxKernels]
  [mode]
    type = FunctionAux
    variable = mode
    function = 'if(t<=1.0, 1.0, 0.0)' # initial "dump" then "normal"
    execute_on = 'timestep_begin'
  []
[]

[GlobalParams]
  point = '0 0 0'
[]
[Postprocessors]
  [solvent_kg]
    type = PointValue
    variable = 'kg_solvent_H2O'
  []
  [dolomite]
    type = PointValue
    variable = 'free_cm3_Dolomite'
  []
  [gypsum]
    type = PointValue
    variable = 'free_cm3_Gypsum'
  []
  [halite]
    type = PointValue
    variable = 'free_cm3_Halite'
  []
  [mirabilite]
    type = PointValue
    variable = 'free_cm3_Mirabilite'
  []
[]

[Outputs]
  csv = true
[]
(modules/geochemistry/test/tests/time_dependent_reactions/seawater_evaporation_no_flow_through.i)

GWB input file

The equivalent Geochemists Workbench input file is

# React script that is equivalent to seawater_evaporation_flow_through.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
temperature = 25 C
H2O          = 1 free kg
Cl-          = 0.5656 mol
balance on Cl-
Na+          = 0.4850 mol
SO4--        = 0.02924 mol
Mg++         = 0.05501 mol
Ca++         = 0.01063 mol
K+           = 0.010576055 mol
HCO3-        = 0.002412
swap CO2(g) for H+
log f CO2(g) = -3.5
fix fugacity of CO2(g)
printout  species = long
epsilon = 1e-13
react -55 mol of H2O
delxi = 0.001
dump
go
flow-through
go
(modules/geochemistry/test/tests/time_dependent_reactions/seawater_evaporation.rea)

Results

Using the HMW activity model means that the results of case 1 are very different from case 2, as shown in Figures 24.7, 24.8 and 24.9 of Bethke (2007). Using the Debye-Huckel model means the results of both cases are quite similar: the impact of the flow-through is evident towards the end of the simulation on the mirabilite mineral in particular. The Geochemists Workbench and the geochemistry module results are shown below.

Figure 1: Minerals precipitated as seawater evaporates and flow-through is used.

Figure 2: Minerals precipitated and dissolved as seawater evaporates (no flow-through is used).

References

  1. Craig M. Bethke. Geochemical and Biogeochemical Reaction Modeling. Cambridge University Press, 2 edition, 2007. doi:10.1017/CBO9780511619670.[BibTeX]