Adding fluids with different temperatures

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

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.

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

The TimeDependentReactionSolver needs some discussion:

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  swap_into_basis = "H2S(aq)"
  swap_out_of_basis = "O2(aq)"
  charge_balance_species = "Cl-"
  constraint_species = "H2O              H+            Cl-              Na+              Mg++             SO4--            Ca++             K+               HCO3-            Ba++            SiO2(aq)          Sr++             Zn++             Cu+              Al+++            Fe++             Mn++             H2S(aq)"
  constraint_value = "  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 = "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 = "   kg               dimensionless moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles            moles"
  close_system_at_time = -0.01
  remove_fixed_activity_name = 'H+'
  remove_fixed_activity_time = -0.01
  initial_temperature = 273
  temperature = 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 = "H2O Al+++ Ba++ Ca++ Cl- Cu+ Fe++ H+ HCO3- K+ Mg++ Mn++ Na+ O2(aq) SO4-- SiO2(aq) Sr++ Zn++"
  source_species_rates = "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 = mode
  execute_console_output_on = '' # only CSV output needed for this example
  stoichiometric_ionic_str_using_Cl_only = true # for comparison with GWB
[]

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]
  type = Transient
  start_time = -0.01 # to allow initial dump to occur
  [TimeStepper]
    type = FunctionDT
    function = timestepper
  []
  end_time = 10
[]

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

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

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

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

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