Adding fluids with different temperatures

Section 22.2 of Bethke (2007) describes the mixing of different fluids, which Bethke terms as "pickup".

Step 1

The simulation begins with forming an equilibrium solution of seawater. The major element composition is shown in Table 1. In addition:

  • C;

  • pH ;

  • the following minerals are prevented from precipitating: Quartz, Tridymite, Cristobalite, Chalcedony and Hematite.

The system is equilibrated, and the final composition without precipitated minerals is saved

Table 1: Major element composition of seawater

SpeciesConcentration (mmolal)
Cl559
Na480
Mg54.5
SO29.5
Ca10.5
K10.1
HCO2.4
SiO(aq)0.17
Sr0.09
Ba
Zn
Al
Cu
Fe
Mn
O(aq)0.123 (free)

Step 2

The hot hydrothermal fluid has composition shown in Table 2. In addition:

  • the temperature is 273C;

  • the pH is 4;

  • HS(aq) is used in the basis instead of O(aq).

The system is brought to equilibrium, then the minerals are allowed to precipitate, and then all precipitated minerals are dumped.

Table 2: Major element composition of hydrothermal fluid

SpeciesConcentration (mmolal)
Cl600
Na529
Mg
SO
Ca21.6
K26.7
HCO2.0
Ba
SiO(aq)20.2
Sr
Zn
Cu
Al
Fe
Mn
HS(aq)6.81

Step 3

The seawater at 4C is slowly mixed into the geothermal fluid at 273C until the final ratio is 10(seawater):1(geothermal). A constant heat capacity is assumed.

MOOSE input file

MOOSE input files dealing with the equilibration (seawater_mixing_step1.i and seawater_mixing_step2.i) may be found in the test suite. The focus here is on the mixing. The GeochemicalModelDefinition defines the basis species, the relevant minerals and the gas (so that its fugacity is computed by the simulation):

[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- Na+ Mg++ SO4-- Ca++ K+ HCO3- Ba++ SiO2(aq) Sr++ Zn++ Cu+ Al+++ Fe++ Mn++ O2(aq)"
    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"}>>> = "Anhydrite Pyrite Talc Amrph^silica Barite Dolomite-ord Muscovite Nontronit-Na Pyrolusite Strontianite"
    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"}>>> = "O2(g)"
  []
[]
(modules/geochemistry/test/tests/time_dependent_reactions/mixing.i)

The TimeDependentReactionSolver needs some discussion:

[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_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"}>>> = "H2S(aq)"
  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"}>>> = "O2(aq)"
  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-"
  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              H+            Cl-              Na+              Mg++             SO4--            Ca++             K+               HCO3-            Ba++            SiO2(aq)          Sr++             Zn++             Cu+              Al+++            Fe++             Mn++             H2S(aq)"
  constraint_value<<<{"description": "Numerical value of the containts on constraint_species"}>>> = "  1.0              6.309573E-5   600E-3           529E-3           0.01E-6          0.01E-6          21.6E-3          26.7E-3          2.0E-3           15E-6           20.2E-3           100.5E-6         41E-6            0.02E-6          4.1E-6           903E-6           1039E-6          6.81E-3"
  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 activity      bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition bulk_composition 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            moles            moles            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.01
  remove_fixed_activity_name<<<{"description": "The name of the species that should have their activity or fugacity constraint removed at time given in remove_fixed_activity_time.  There should be an equal number of these names as times given in remove_fixed_activity_time.  Each of these must be in the basis and have an activity or fugacity constraint"}>>> = 'H+'
  remove_fixed_activity_time<<<{"description": "The times at which the species in remove_fixed_activity_name should have their activity or fugacity constraint removed."}>>> = -0.01
  initial_temperature<<<{"description": "The initial aqueous solution is equilibrated at this system before adding reactants, changing temperature, etc."}>>> = 273
  temperature<<<{"description": "Temperature.  This has two different meanings if mode!=4.  (1) If no species are being added to the solution (no source_species_rates are positive) then this is the temperature of the aqueous solution.  (2) If species are being added, this is the temperature of the species being added.  In case (2), the final aqueous-solution temperature is computed assuming the species are added, temperature is equilibrated and then, if species are also being removed, they are removed.  If you wish to add species and simultaneously alter the temperature, you will have to use a sequence of heat-add-heat-add, etc steps.  In the case that mode=4, temperature is the final temperature of the aqueous solution"}>>> = T
  # The following source species and rates are taken from the Geochemists Workbench (see output from mixing.rea)
  # An alternative is to run the seawater_mixing MOOSE input files and extract the source species and rates
  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 Al+++ Ba++ Ca++ Cl- Cu+ Fe++ H+ HCO3- K+ Mg++ Mn++ Na+ O2(aq) SO4-- SiO2(aq) Sr++ Zn++"
  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!"}>>> = "H2O_rate Al+++_rate Ba++_rate Ca++_rate Cl-_rate Cu+_rate Fe++_rate H+_rate HCO3-_rate K+_rate Mg++_rate Mn++_rate Na+_rate O2aq_rate SO4--_rate SiO2aq_rate Sr++_rate Zn++_rate"
  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
  execute_console_output_on<<<{"description": "When to execute the geochemistry console output"}>>> = '' # only CSV output needed for this 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 comparison with GWB
[]
(modules/geochemistry/test/tests/time_dependent_reactions/mixing.i)

From the top to the bottom of this block:

  • the swaps are defined

  • the composition of the fluid is defined via the bulk mole composition of the species as well as the activity of H which sets the pH

  • the system is closed at . This means that all "free" type constraints are converted into "bulk" type constraints. The only such constraint in this case is on HO, so after no more HO can be added to the system to maintain the kg_solvent_water constraint.

  • the activity constraint on H is also removed at . This means that no more H can be added to the system to maintain the pH constraint

  • the initial temperature is 273C, which is the temperature at which the system is initially equilibrated at

  • the subsequent temperature is controlled by the T variable (see below)

  • various source species are provided with rates (see below)

  • the mode is set equal to the mode variable (see below)

The Executioner defines the meaning of time:

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  start_time = -0.01 # to allow initial dump to occur
  [TimeStepper<<<{"href": "../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = FunctionDT
    function = timestepper
  []
  end_time = 10
[]
(modules/geochemistry/test/tests/time_dependent_reactions/mixing.i)

The temperature is controlled using:

  [T_auxk]
    type = FunctionAux
    variable = T
    function = 'if(t<=0, 273, 4)' # during initialisation and dumping, T=273, while during adding T=temperature of reactants
    execute_on = timestep_begin
  []
  [H2O_rate_auxk]
    type = FunctionAux
    variable = H2O_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 55.510000000000005)'
  []
  [Al+++_rate]
    type = FunctionAux
    variable = Al+++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 3.643e-10)'
  []
  [Ba++_rate]
    type = FunctionAux
    variable = Ba++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 8.831e-08)'
  []
  [Ca++_rate]
    type = FunctionAux
    variable = Ca++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.0104)'
  []
  [Cl-_rate]
    type = FunctionAux
    variable = Cl-_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.559)'
  []
  [Cu+_rate]
    type = FunctionAux
    variable = Cu+_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 7.000000000000001e-09)'
  []
  [Fe++_rate]
    type = FunctionAux
    variable = Fe++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 4.746e-15)'
  []
  [H+_rate]
    type = FunctionAux
    variable = H+_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.0002005)'
  []
  [HCO3-_rate]
    type = FunctionAux
    variable = HCO3-_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.002153)'
  []
  [K+_rate]
    type = FunctionAux
    variable = K+_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.010100000000000001)'
  []
  [Mg++_rate]
    type = FunctionAux
    variable = Mg++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.054400000000000004)'
  []
  [Mn++_rate]
    type = FunctionAux
    variable = Mn++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 6.79e-14)'
  []
  [Na+_rate]
    type = FunctionAux
    variable = Na+_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.48019999999999996)'
  []
  [O2aq_rate]
    type = FunctionAux
    variable = O2aq_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.000123)'
  []
  [SO4--_rate]
    type = FunctionAux
    variable = SO4--_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.0295)'
  []
  [SiO2aq_rate]
    type = FunctionAux
    variable = SiO2aq_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 0.00017)'
  []
  [Sr++_rate]
    type = FunctionAux
    variable = Sr++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 3.8350000000000004e-05)'
  []
  [Zn++_rate]
    type = FunctionAux
    variable = Zn++_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 1e-08)'
  []
[]

[Postprocessors]
  [temperature]
    type = PointValue
    point = '0 0 0'
    variable = "solution_temperature"
  []
  [fugactity_O2]
    type = PointValue
    point = '0 0 0'
    variable = "activity_O2(g)"
  []
  [molal_SO4--]
    type = PointValue
    point = '0 0 0'
    variable = "molal_SO4--"
  []
  [molal_NaSO4]
    type = PointValue
    point = '0 0 0'
    variable = "molal_NaSO4-"
  []
  [molal_H2Saq]
    type = PointValue
    point = '0 0 0'
    variable = "molal_H2S(aq)"
  []
  [molal_HSO4-]
    type = PointValue
    point = '0 0 0'
    variable = "molal_HSO4-"
  []
  [cm3_Anhydrite]
    type = PointValue
    point = '0 0 0'
    variable = "free_cm3_Anhydrite"
  []
  [cm3_Pyrite]
    type = PointValue
    point = '0 0 0'
    variable = "free_cm3_Pyrite"
  []
  [cm3_Talc]
    type = PointValue
    point = '0 0 0'
    variable = "free_cm3_Talc"
  []
  [cm3_AmSil]
    type = PointValue
    point = '0 0 0'
    variable = "free_cm3_Amrph^silica"
  []
[]

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

[Executioner]
  type = Transient
  start_time = -0.01 # to allow initial dump to occur
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 10
[]

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Cl- Na+ Mg++ SO4-- Ca++ K+ HCO3- Ba++ SiO2(aq) Sr++ Zn++ Cu+ Al+++ Fe++ Mn++ O2(aq)"
    equilibrium_minerals = "Anhydrite Pyrite Talc Amrph^silica Barite Dolomite-ord Muscovite Nontronit-Na Pyrolusite Strontianite"
    equilibrium_gases = "O2(g)"
  []
[]

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

This means that during initialization the temperature is 273C, and whenever the source species' rates are nonzero (which is for in this case) then the reactant temperature is 4C.

The mode is set to "dump" for , and otherwise no special mode, by:

  [mode_auxk]
    type = FunctionAux
    variable = mode
    function = 'if(t<=0, 1, 0)' # dump at start of first timestep
    execute_on = timestep_begin
  []
(modules/geochemistry/test/tests/time_dependent_reactions/mixing.i)

All the source-species rates follow the same pattern:

  [H2O_rate_auxk]
    type = FunctionAux
    variable = H2O_rate
    execute_on = timestep_begin
    function = 'if(t<=0, 0, 55.510000000000005)'
  []
(modules/geochemistry/test/tests/time_dependent_reactions/mixing.i)

The numerical values of the rates were derived from the Geochemists Workbench input file (below) so that the results match as closely as possible to the GWB results, but they could equally be derived from seawater_mixing_step1.i and seawater_mixing_step2.i.

GWB input file

The equivalent Geochemists Workbench input file is

# React script, the second half of which is equivalent to mixing.i
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
suppress all
unsuppress Anhydrite Pyrite Talc Amrph^silica Barite Dolomite-ord Muscovite Nontronit-Na Pyrolusite Strontianite
temperature = 4 C
H2O          = 1 free kg
H+           = 8.1 pH
Cl-          = 559 mmolal
balance on Cl-
Na+          = 480.2 mmolal
Mg++         = 54.5 mmolal
SO4--        = 29.5 mmolal
Ca++         = 10.5 mmolal
K+           = 10.1 mmolal
HCO3-        = 2.4 mmolal
SiO2(aq)     = 0.17 mmolal
Sr++         = 0.09 mmolal

Ba++         = 0.2 umolal
Zn++         = 0.01 umolal
Al+++        = 0.005 umolal
Cu+          = 0.007 umolal
Fe++         = 0.001 umolal
Mn++         = 0.001 umolal
O2(aq)       = 123.0 free umolal
printout  species = long
epsilon = 1e-14
go

pickup reactant = fluid

# Below here is equivalent to mixing.i  The material above here provides the rates to mixing.i

reactant times 10

H+       = 4.2 pH
swap H2S(aq) for O2(aq)

Cl-      = 600.0 mmolal
Na+      = 529.0 mmolal
K+       = 26.7 mmolal
Ca++     = 21.6 mmolal
SiO2(aq) = 20.2 mmolal
H2S(aq)  = 6.81 mmolal
HCO3-    = 2.0 mmolal

Mn++  = 1039.0 umolal
Fe++  = 903.0 umolal
Sr++  = 100.5 umolal
Zn++  = 41.0 umolal
Ba++  = 15.0 umolal
Al+++ = 4.1 umolal
Cu+   = 0.02 umolal
Mg++  = 0.01 umolal
SO4-- = 0.01 umolal

dump
T initial 273, reactants = 4
delxi = 0.001
go
(modules/geochemistry/test/tests/time_dependent_reactions/mixing.rea)

Results

Figures 22.3 and 22.4 of Bethke (2007) show the results, which may be compared with the results below. Note, there are two reasons why the two software packages produce slightly different results.

  1. Most importantly, the interpolations of equilibrium constants at high temperature appears to be different in the two software packages. This may be verified by running the other simulations in the tests and examples at higher temperatures and observing that the results begin to differ at higher temperatures, despite being identical at lower temperatures. In the current situation, this mostly impacts the system at C and causes aspects of the models to differ.

  2. Less importantly, during the initial equilibration process, Geochemists Workbench first equilibrates the solution, then removes the pH constraint, then allows precipitates to form, then dumps the precipitates. Conversely, the geochemistry module first equilibrates the solution allowing precipitates to form, then dumps the precipitates and removes the pH constraint. The difference between these two approaches is that during the precipitation geochemistry is adding/removing HCl to the aqueous solution to maintain the user-specified pH constraint. This means the initial configuration before adding the seawater is slightly different, and explains the different behaviour of the Talc mineral during the early stages of the simulation (which may be verified by creating another model with fixed bulk mole number for H).

Figure 1: Aqueous solution temperature.

Figure 2: Precipitates formed.

Figure 3: Oxygen fugacity.

Figure 4: Species concentration.

References

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