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<<<{"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 H+ Cl- Ca++ Mg++ Na+ K+ SO4-- HCO3-"
    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"}>>> = "Dolomite Epsomite Gypsum Halite Magnesite Mirabilite Sylvite"
    equilibrium_gases<<<{"description": "A list of gases that are in equilibrium with the aqueous solution and can have their fugacities fixed, at least for some time and spatial location.  All members of this list must be in the 'gas' section of the database file"}>>> = "CO2(g)"
    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 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<<<{"href": "../../../syntax/TimeDependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  geochemistry_reactor_name<<<{"description": "The name that will be given to the GeochemistryReactor UserObject built by this action"}>>> = reactor
  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"}>>> = "H+"
  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"}>>> = "  CO2(g)"
  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."}>>> = "Cl-" # this means the bulk moles of Cl- will not be exactly as set below
  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              CO2(g)        Cl-              Na+              SO4--            Mg++             Ca++             K+               HCO3-"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -3.5          0.5656           0.4850           0.02924          0.05501          0.01063          0.010576055      0.002412"
  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 log10fugacity bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition"
  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               dimensionless moles            moles            moles            moles            moles            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."}>>> = "H2O"
  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" # 1kg H2O = 55.51 moles, each time step removes 1 mole
  mode<<<{"description": "This may vary temporally.  If mode=1 then 'dump' mode is used, which means all non-kinetic mineral masses are removed from the system before the equilibrium solution is sought (ie, removal occurs at the beginning of the time step).  If mode=2 then 'flow-through' mode is used, which means all mineral masses are removed from the system after it the equilbrium solution has been found (ie, at the end of a time step).  If mode=3 then 'flush' mode is used, then before the equilibrium solution is sought (ie, at the start of a time step) water+species is removed from the system at the same rate as pure water + non-mineral solutes are entering the system (specified in source_species_rates).  If mode=4 then 'heat-exchanger' mode is used, which means the entire current aqueous solution is removed, then the source_species are added, then the temperature is set to 'cold_temperature', the system is solved and any precipitated minerals are removed, then the temperature is set to 'temperature', the system re-solved and any precipitated minerals are removed.  If mode is any other number, no special mode is active (the system simply responds to the source_species_rates, controlled_activity_value, etc)."}>>> = mode
  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 example
  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 precise agreement with GWB
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # 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<<<{"href": "../../../syntax/Functions/index.html"}>>>]
  [timestepper]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../../source/functions/PiecewiseLinear.html"}>>>
    x<<<{"description": "The abscissa values"}>>> = '0 50 55'
    y<<<{"description": "The ordinate values"}>>> = '5 5 1'
  []
[]
[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  [TimeStepper<<<{"href": "../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = timestepper
  []
  end_time = 55
[]

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [mode]
  []
  [dolomite_mol]
  []
  [halite_mol]
  []
  [gypsum_mol]
  []
  [mirabilite_mol]
  []
[]

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [mode_auxk]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = mode
    function<<<{"description": "The function to use as the value"}>>> = 'if(t<=1.0, 1.0, 2.0)' # initial "dump" then "flow_through"
    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_begin'
  []
  [dolomite_mol_auxk]
    type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
    reactor<<<{"description": "The name of the Geochemistry*Reactor user object."}>>> = reactor
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = dolomite_mol
    species<<<{"description": "Species name"}>>> = Dolomite
    quantity<<<{"description": "Type of quantity to output.  These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis).  These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water).  These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active.  These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V).  These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep).  If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = moles_dumped
  []
  [gypsum_mol_auxk]
    type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
    reactor<<<{"description": "The name of the Geochemistry*Reactor user object."}>>> = reactor
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = gypsum_mol
    species<<<{"description": "Species name"}>>> = Gypsum
    quantity<<<{"description": "Type of quantity to output.  These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis).  These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water).  These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active.  These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V).  These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep).  If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = moles_dumped
  []
  [halite_mol]
    type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
    reactor<<<{"description": "The name of the Geochemistry*Reactor user object."}>>> = reactor
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = halite_mol
    species<<<{"description": "Species name"}>>> = Halite
    quantity<<<{"description": "Type of quantity to output.  These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis).  These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water).  These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active.  These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V).  These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep).  If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = moles_dumped
  []
  [mirabilite_mol]
    type = GeochemistryQuantityAux<<<{"description": "Extracts information from the Reactor and records it in the AuxVariable", "href": "../../../source/auxkernels/GeochemistryQuantityAux.html"}>>>
    reactor<<<{"description": "The name of the Geochemistry*Reactor user object."}>>> = reactor
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = mirabilite_mol
    species<<<{"description": "Species name"}>>> = Mirabilite
    quantity<<<{"description": "Type of quantity to output.  These are available for non-kinetic species: activity (which equals fugacity for gases); bulk moles (this will be zero if the species is not in the basis); neglog10a = -log10(activity); transported_moles_in_original_basis (will throw an error if species is not in original basis).  These are available only for non-kinetic non-minerals: molal = mol(species)/kg(solvent_water); mg_per_kg = mg(species)/kg(solvent_water).  These are available only for minerals: free_mg = free mg; free_cm3 = free cubic-centimeters; moles_dumped = moles dumped when special dump and flow-through modes are active.  These are available for minerals that host sorbing sites: surface_charge (C/m^2); surface_potential (V).  These are available for kinetic species: kinetic_moles; kinetic_additions (-dt * rate = mole increment in kinetic species for this timestep).  If quantity=temperature, then the 'species' is ignored and the AuxVariable will record the aqueous solution temperature in degC"}>>> = moles_dumped
  []
[]

[GlobalParams<<<{"href": "../../../syntax/GlobalParams/index.html"}>>>]
  point = '0 0 0'
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [solvent_kg]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'kg_solvent_H2O'
  []
  [dolomite_mol]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = dolomite_mol
  []
  [gypsum_mol]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'gypsum_mol'
  []
  [halite_mol]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'halite_mol'
  []
  [mirabilite_mol]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = '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<<<{"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 H+ Cl- Ca++ Mg++ Na+ K+ SO4-- HCO3-"
    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"}>>> = "Dolomite Epsomite Gypsum Halite Magnesite Mirabilite Sylvite"
    equilibrium_gases<<<{"description": "A list of gases that are in equilibrium with the aqueous solution and can have their fugacities fixed, at least for some time and spatial location.  All members of this list must be in the 'gas' section of the database file"}>>> = "CO2(g)"
    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 precise agreement with GWB
  []
[]

[TimeDependentReactionSolver<<<{"href": "../../../syntax/TimeDependentReactionSolver/index.html"}>>>]
  model_definition<<<{"description": "The name of the GeochemicalModelDefinition user object (you must create this UserObject yourself)"}>>> = definition
  geochemistry_reactor_name<<<{"description": "The name that will be given to the GeochemistryReactor UserObject built by this action"}>>> = reactor
  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"}>>> = "H+"
  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"}>>> = "  CO2(g)"
  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."}>>> = "Cl-" # this means the bulk moles of Cl- will not be exactly as set below
  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              CO2(g)        Cl-              Na+              SO4--            Mg++             Ca++             K+               HCO3-"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              -3.5          0.5656           0.4850           0.02924          0.05501          0.01063          0.010576055      0.002412"
  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 log10fugacity bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition"
  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               dimensionless moles            moles            moles            moles            moles            moles            moles"
  close_system_at_time<<<{"description": "Time at which to 'close' the system, that is, change a kg_solvent_water constraint to moles_bulk_water, and all free_molality and free_moles_mineral_species to moles_bulk_species"}>>> = 0
  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."}>>> = "H2O"
  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" # 1kg H2O = 55.51 moles, each time step removes 1 mole
  mode<<<{"description": "This may vary temporally.  If mode=1 then 'dump' mode is used, which means all non-kinetic mineral masses are removed from the system before the equilibrium solution is sought (ie, removal occurs at the beginning of the time step).  If mode=2 then 'flow-through' mode is used, which means all mineral masses are removed from the system after it the equilbrium solution has been found (ie, at the end of a time step).  If mode=3 then 'flush' mode is used, then before the equilibrium solution is sought (ie, at the start of a time step) water+species is removed from the system at the same rate as pure water + non-mineral solutes are entering the system (specified in source_species_rates).  If mode=4 then 'heat-exchanger' mode is used, which means the entire current aqueous solution is removed, then the source_species are added, then the temperature is set to 'cold_temperature', the system is solved and any precipitated minerals are removed, then the temperature is set to 'temperature', the system re-solved and any precipitated minerals are removed.  If mode is any other number, no special mode is active (the system simply responds to the source_species_rates, controlled_activity_value, etc)."}>>> = mode
  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 example
  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 precise agreement with GWB
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output for this example
[]

[Functions<<<{"href": "../../../syntax/Functions/index.html"}>>>]
  [timestepper]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../../../source/functions/PiecewiseLinear.html"}>>>
    x<<<{"description": "The abscissa values"}>>> = '0 50 55'
    y<<<{"description": "The ordinate values"}>>> = '5 5 1'
  []
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  [TimeStepper<<<{"href": "../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = timestepper
  []
  end_time = 55
[]

[AuxVariables<<<{"href": "../../../syntax/AuxVariables/index.html"}>>>]
  [mode]
  []
[]

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [mode]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = mode
    function<<<{"description": "The function to use as the value"}>>> = 'if(t<=1.0, 1.0, 0.0)' # initial "dump" then "normal"
    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_begin'
  []
[]

[GlobalParams<<<{"href": "../../../syntax/GlobalParams/index.html"}>>>]
  point = '0 0 0'
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [solvent_kg]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'kg_solvent_H2O'
  []
  [dolomite]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Dolomite'
  []
  [gypsum]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Gypsum'
  []
  [halite]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Halite'
  []
  [mirabilite]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = 'free_cm3_Mirabilite'
  []
[]

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