A Geothermal simulation involving 2D flow

Representative water composition

Water composition has been measured at Roosevelt Hot Springs, which is close to the FORGE site, and this may be representative of the type of water found in the FORGE fractured granite. The other type of water used in this presentation is a low TDS potable water, which is the injectate. The composition of both waters is shown in Table 1 (after Patel and Simmons (2020)). Note that there is no iron composition in either of the waters, so no iron-bearing minerals will be considered in this presentation.

Mineralogy

X-ray diffraction suggests the composition of the FORGE granite is as shown in Table 2.

Patel and Simmons (2020) interpret Plagioclase as a 9:1 mix of Albite and Anorthite, while Stuart Simmons suggests that Biotite could be Phlogopite. Further correspondence with Stuart Simmons has identified a set of representative minerals to be used in the simulations, and they are shown in Table 3

Water1 at equilibrium

A geochemistry input file that finds the equilibrium solution for Water1 is

# Equilibrium model "Water 1" from "Subtask 2C.4.7 Geochemical Modeling SSimmons-VPatil.pdf"
# At 60degC K-feldspar and Quartz both precipitate
[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = '../../../../geochemistry/database/moose_geochemdb.json'
    basis_species = 'H2O H+ Na+ K+ Ca++ Mg++ SiO2(aq) Al+++ Cl- SO4-- HCO3-'
    remove_all_extrapolated_secondary_species = true
    equilibrium_minerals = 'Albite Anhydrite Anorthite Calcite Chalcedony Clinochl-7A Illite K-feldspar Kaolinite Quartz Paragonite Phlogopite Zoisite Laumontite'
  []
[]

[TimeIndependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = 'Cl-'
  constraint_species = 'H2O H+      Na+  K+    Ca++    Mg++      SiO2(aq) Al+++    Cl-  SO4--  HCO3-'
  constraint_value = '  1.0 3.16E-8 0.12 0.016 0.68E-3 0.0008E-3 3.7E-3   0.004E-3 0.15 0.5E-3 1.4E-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'
  constraint_unit = 'kg dimensionless moles moles moles moles moles moles moles moles moles'
  temperature = 60
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
[]

[Postprocessors]
  [bulk_H+]
    type = PointValue
    point = '0 0 0'
    variable = 'bulk_moles_H+'
  []
[]

[Outputs]
  csv = true
[]

Note use of the flag

remove_all_extrapolated_secondary_species = true

This flag removes the odd secondary species whose equilibrium constants have only been measured for a small range of temperatures. By default, the geochemistry module extrapolates the equilibrium constants for such species, and frequently this means the equilibrium constants are unrealistically high or low (such as or ) at high temperatures. This leads to convergence problems, or the simulation converges to an unrealistic solution. Therefore, this flag is used for all input files in this presentation.

If minerals are allowed to precipitate from Water1, small amounts of K-feldspar and Quartz form.

Estimating the in-situ reservoir water

When performing reactive-transport simulations, it is necessary to start with a fluid-filled reservoir. For most purposes, the fluid composition is not critical, since the injectate will rapidly displace it. However, it is convenient to fill the initial reservoir with fluid that is at equilibrium with the reservoir rocks. This makes the interpretation of the reactive-transport simulations easier, for any changes in mineralogy can be attributed to the injectate, and not the in-situ reservoir water. Therefore, this section attempts to find a water composition that is at equilibrium with the reservoir.

In this section, all the minerals present in Table 3 are used, with the exception of Laumontite and Zoisite. These two are present only in trace quantities (if at all), but if they exist in the model, the rock mineralogy of Table 3 is unstable, so it is impossible to find an equilibrium. For instance, Anorthite will decompose into Laumontite. This suggests that the rock mineralogy is controlled by kinetics, which is only considered in later sections (where Laumontite and Zoisite are included). It also reveals that the modelling results could be significantly impacted by the choice of minerals.

The following method is used to estimate a water composition that is at equilibrium with the reservoir mineralogy.

1. Water1 is equilibrated at 60C, with pH fixed to 7.5, allowing any precipitates to form. Only Quartz and K-feldspar precipitate, as mentioned in the previous section.

2. The system is closed (at ), i.e. the pH is no longer fixed. The small amount of Quartz and K-feldspar precipitates are retained.

3. The temperature is raised to 220C (during ), allowing any precipitates to form or dissolve. Quartz dissolves entirely, K-feldspar precipitate remains, and Calcite and Phlogopite precipitate. The pH becomes 7.078. Note the use of remove_all_extrapolated_secondary_species = true in the GeochemicalModelDefinition. If the extrapolated secondary species are retained instead, the results are significantly different (and probably incorrect).

4. The following minerals are added (during ): Albite (16.8mol = 44% by weight), Anorthite (1.8mol = 5% by weight), K-feldspar (10.4mol = 29% by weight), Quartz (30.0mol = 18% by weight), Phlogopite (0.48mol = 2% by weight) and Illite (0.52mol = 2% by weight). The mol numbers are approximately what has been measured by X-ray diffraction as in Table 3.

# Minerals suggested by Stuart Simmons, but I do not include Laumontite and Zoisite as they are more stable than Anorthite so all Anorthite becomes one of these minerals which contradicts the XRD observations.  All minerals are considered in the kinetic models.
# Model of "Water 1" from "Subtask 2C.4.7 Geochemical Modeling SSimmons-VPatil.pdf" subjected to the following:
# (1) The system is equilibrated at 60deg, with pH fixed to 7.5, allowing any precipitates to form.  Note that the only minerals present in the system are those mentioned in "Subtask 2C.4.7 Geochemical Modeling SSimmons-VPatil.pdf".  If other minerals are present, the results change significantly.  Only Quartz and K-feldspar precipitate.
# (2) The system is closed (at time=0), ie the pH is no longer fixed.  The Quartz and K-feldspar precipitates are retained
# (3) The temperature is raised to 220degC (during 0<time<=1), allowing any precipitates to form or dissolve.  Quartz dissolves entirely, K-feldspar precipitate remains, and Calcite and Phlogopite precipitate.  The pH becomes 7.078.  Note the use of remove_all_extrapolated_secondary_species = true in the GeochemicalModelDefinition.  If the extrapolated secondary species are retained instead, the results are significantly different.
# (4) The following minerals are added (during 1<time<=2): Albite (16.8mol = 44% by weight), Anorthite (1.8mol = 5% by weight), K-feldspar (10.4mol = 29% by weight), Quartz (30.0mol = 18% by weight), Phlogopite (0.48mol = 2% by weight) and Illite (0.52mol = 2% by weight).  The mol numbers are approximately what has been measured by XRD, but it is not important to specify the exact composition of the rock (that will be done in the kinetic simulations): what is important here is that there is *some* precipitate.
# (5) The free moles precipitated are Albite 16.38, Anorthite 1.785, K-feldspar 10.68, Quartz 30.82, Phlogopite 0.52, Paragonite 0.44, Calcite 0.0004, Anhydrite 0.0004, Chalcedony 0, Illite 0, Kaolinite 0, Clinochl-7A 0.  Calcite is constrained by the initial HCO3- concentration and Anhydrite by the initial SO4-- concentration, and both have only been observed in trace quantities in agreement with this simulation
# (6) The free mole numbers of the basis species that are now in equilibrium with the minerals are extracted, which is the key output of this simulation.  Note that the original composition of "Water 1" is largely irrelevant.  As mentioned, the HCO3- and SO4-- concentrations constrain Calcite and Anhydrite.  Also, adding the minerals causes the pH to change to 6.16.
[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = '../../../../geochemistry/database/moose_geochemdb.json'
    basis_species = 'H2O H+ Na+ K+ Ca++ Mg++ SiO2(aq) Al+++ Cl- SO4-- HCO3-'
    equilibrium_minerals = 'Albite Anhydrite Anorthite Calcite Chalcedony Clinochl-7A Illite K-feldspar Kaolinite Quartz Paragonite Phlogopite'
    remove_all_extrapolated_secondary_species = true
  []
[]

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = 'Cl-'
  constraint_species = 'H2O H+      Na+  K+    Ca++    Mg++      SiO2(aq) Al+++    Cl-  SO4--  HCO3-'
  constraint_value = '  1.0 3.16E-8 0.12 0.016 0.68E-3 0.0008E-3 3.7E-3   0.004E-3 0.15 0.5E-3 1.4E-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'
  constraint_unit = 'kg dimensionless moles moles moles moles moles moles moles moles moles'
  initial_temperature = 60
  remove_fixed_activity_name = 'H+'
  remove_fixed_activity_time = 0
  temperature = 220
  source_species_names = 'Albite Anorthite K-feldspar Quartz Phlogopite Illite'
  source_species_rates = 'Albite_rate Anorthite_rate K-feldspar_rate Quartz_rate Phlogopite_rate Illite_rate'
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  mol_cutoff = 1E-100
  execute_console_output_on = 'timestep_end' # only CSV output
  solver_info = true
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = 'if(t<1, 1, if(t<1.01, 0.01, 1))'
  []
  end_time = 2
[]

[AuxVariables]
  [Albite_rate]
  []
  [Anorthite_rate]
  []
  [K-feldspar_rate]
  []
  [Quartz_rate]
  []
  [Phlogopite_rate]
  []
  [Illite_rate]
  []
  [transported_H2O]
  []
  [transported_H+]
  []
  [transported_Na+]
  []
  [transported_K+]
  []
  [transported_Ca++]
  []
  [transported_Mg++]
  []
  [transported_SiO2]
  []
  [transported_Al+++]
  []
  [transported_Cl-]
  []
  [transported_SO4--]
  []
  [transported_HCO3-]
  []
[]
[AuxKernels]
  [Albite_rate]
    type = FunctionAux
    variable = Albite_rate
    function = 'if(t>1, 16.8, 0)'
    execute_on = timestep_begin
  []
  [Anorthite_rate]
    type = FunctionAux
    variable = Anorthite_rate
    function = 'if(t>1, 1.8, 0)'
    execute_on = timestep_begin
  []
  [K-feldspar_rate]
    type = FunctionAux
    variable = K-feldspar_rate
    function = 'if(t>1, 10.4, 0)'
    execute_on = timestep_begin
  []
  [Quartz_rate]
    type = FunctionAux
    variable = Quartz_rate
    function = 'if(t>1, 30.0, 0)'
    execute_on = timestep_begin
  []
  [Phlogopite_rate]
    type = FunctionAux
    variable = Phlogopite_rate
    function = 'if(t>1, 0.48, 0)'
    execute_on = timestep_begin
  []
  [Illite_rate]
    type = FunctionAux
    variable = Illite_rate
    function = 'if(t>1, 0.52, 0)'
    execute_on = timestep_begin
  []
  [transported_H2O]
    type = GeochemistryQuantityAux
    species = 'H2O'
    variable = transported_H2O
    quantity = transported_moles_in_original_basis
  []
  [transported_H+]
    type = GeochemistryQuantityAux
    species = 'H+'
    variable = transported_H+
    quantity = transported_moles_in_original_basis
  []
  [transported_Na+]
    type = GeochemistryQuantityAux
    species = 'Na+'
    variable = transported_Na+
    quantity = transported_moles_in_original_basis
  []
  [transported_K+]
    type = GeochemistryQuantityAux
    species = 'K+'
    variable = transported_K+
    quantity = transported_moles_in_original_basis
  []
  [transported_Ca++]
    type = GeochemistryQuantityAux
    species = 'Ca++'
    variable = transported_Ca++
    quantity = transported_moles_in_original_basis
  []
  [transported_Mg++]
    type = GeochemistryQuantityAux
    species = 'Mg++'
    variable = transported_Mg++
    quantity = transported_moles_in_original_basis
  []
  [transported_SiO2]
    type = GeochemistryQuantityAux
    species = 'SiO2(aq)'
    variable = transported_SiO2
    quantity = transported_moles_in_original_basis
  []
  [transported_Al+++]
    type = GeochemistryQuantityAux
    species = 'Al+++'
    variable = transported_Al+++
    quantity = transported_moles_in_original_basis
  []
  [transported_Cl-]
    type = GeochemistryQuantityAux
    species = 'Cl-'
    variable = transported_Cl-
    quantity = transported_moles_in_original_basis
  []
  [transported_SO4--]
    type = GeochemistryQuantityAux
    species = 'SO4--'
    variable = transported_SO4--
    quantity = transported_moles_in_original_basis
  []
  [transported_HCO3-]
    type = GeochemistryQuantityAux
    species = 'HCO3-'
    variable = transported_HCO3-
    quantity = transported_moles_in_original_basis
  []
[]
[GlobalParams]
  point = '0 0 0'
  reactor = reactor
[]
[Postprocessors]
  [kg_solvent_water]
    type = PointValue
    variable = kg_solvent_H2O
  []
  [free_cm3_Albite]
    type = PointValue
    variable = free_cm3_Albite
  []
  [free_cm3_Anhydrite]
    type = PointValue
    variable = free_cm3_Anhydrite
  []
  [free_cm3_Anorthite]
    type = PointValue
    variable = free_cm3_Anorthite
  []
  [free_cm3_Calcite]
    type = PointValue
    variable = free_cm3_Calcite
  []
  [free_cm3_Chalcedony]
    type = PointValue
    variable = free_cm3_Chalcedony
  []
  [free_cm3_Clinochl-7A]
    type = PointValue
    variable = free_cm3_Clinochl-7A
  []
  [free_cm3_Illite]
    type = PointValue
    variable = free_cm3_Illite
  []
  [free_cm3_K-feldspar]
    type = PointValue
    variable = free_cm3_K-feldspar
  []
  [free_cm3_Kaolinite]
    type = PointValue
    variable = free_cm3_Kaolinite
  []
  [free_cm3_Quartz]
    type = PointValue
    variable = free_cm3_Quartz
  []
  [free_cm3_Paragonite]
    type = PointValue
    variable = free_cm3_Paragonite
  []
  [free_cm3_Phlogopite]
    type = PointValue
    variable = free_cm3_Phlogopite
  []
  [molal_H+]
    type = PointValue
    variable = molal_H+
  []
  [molal_Na+]
    type = PointValue
    variable = molal_Na+
  []
  [molal_K+]
    type = PointValue
    variable = molal_K+
  []
  [molal_Ca++]
    type = PointValue
    variable = molal_Ca++
  []
  [molal_Mg++]
    type = PointValue
    variable = molal_Mg++
  []
  [molal_SiO2]
    type = PointValue
    variable = molal_SiO2(aq)
  []
  [molal_Al+++]
    type = PointValue
    variable = molal_Al+++
  []
  [molal_SO4--]
    type = PointValue
    variable = molal_SO4--
  []
  [molal_HCO3-]
    type = PointValue
    variable = molal_HCO3-
  []
  [bulk_moles_Cl-]
    type = PointValue
    variable = bulk_moles_Cl-
  []
  [transported_H2O]
    type = PointValue
    variable = transported_H2O
  []
  [transported_H+]
    type = PointValue
    variable = transported_H+
  []
  [transported_Na+]
    type = PointValue
    variable = transported_Na+
  []
  [transported_K+]
    type = PointValue
    variable = transported_K+
  []
  [transported_Ca++]
    type = PointValue
    variable = transported_Ca++
  []
  [transported_Mg++]
    type = PointValue
    variable = transported_Mg++
  []
  [transported_SiO2]
    type = PointValue
    variable = transported_SiO2
  []
  [transported_Al+++]
    type = PointValue
    variable = transported_Al+++
  []
  [transported_Cl-]
    type = PointValue
    variable = transported_Cl-
  []
  [transported_SO4--]
    type = PointValue
    variable = transported_SO4--
  []
  [transported_HCO3-]
    type = PointValue
    variable = transported_HCO3-
  []
  [pH]
    type = PointValue
    variable = pH
  []
[]
[Outputs]
  csv = true
[]

The free moles precipitated are Albite 16.38, Anorthite 1.785, K-feldspar 10.68, Quartz 30.82, Phlogopite 0.52, Paragonite 0.44, Calcite 0.0004, Anhydrite 0.0004, while the rest of the minerals do not precipitate. Notice:

• The added Illite has dissolved in favor of Paragonite. This is important in interpreting reactive-transport results because if the reactive-transport's initial mineralogy is given by Table 3, Illite could dissolve regardless of the injectate.

• Calcite is constrained by the initial HCO3- concentration

• Anhydrite by the initial SO4– concentration

In-situ kinetics

The analysis of the previous section produced a water that is in equilibrium with the rock mineralogy, as well as a stable mineral assemblage. However:

• The analysis revealed that Table 3 may not be in equilibrium because Illite dissolves. Perhaps this does not matter greatly as Illite's mass fraction is only about 2%.

• The model did not include Laumontite and Zoisite since they are more stable than Anorthite, which appears in appreciable fractions in the X-ray analysis.

Presumably all these minerals are controlled by kinetic reactions, and using kinetics may result in quite different results. Therefore, this section performs an analysis of the in-situ reservoir kinetics to act as a baseline for the reactive-transport simulations.

In this presentation, kinetic rates are modelled using the formula (1) where

• [mol.s] is the reaction rate

• [g] is the mass of the mineral. This is the number of moles multiplied by the molar weight.

• [cm.g] is the specific surface area of the mineral

• [mol.s.m] is the reaction's intrinsic rate constant

• [J.mol] is the reaction's activation area

• J.K.mol is the gas constant

• [K] is the temperature

Palandri and Kharaka (2004) list parameters for a wide number of minerals. In many cases these are based on limited experimental data, so may be subject to large errors.

It is useful to analyse the rate of Illite and Anorthite conversion in order to act as a reference to the reactive-transport simulations

• The in-situ reservoir water from the previous section is used

• Minerals with ratios given in Table 3 are introduced

• Trace minerals are assumed to exist at a 0.1% volume fraction. It is necessary to have some mineral present in the simulation for the rate of precipitation and dissolution, Eq. (1), depends on the mineral mass: with zero mineral there can be no precipitation of that mineral!

• The quantities of the trace minerals Calcite and Anhydrite are constrained by the initial HCO3- and SO4– concentration, respectively

• For every litre of in-situ reservoir water, 99 litres of rock mineral is used, representing a porosity of 0.01

• After the minerals are introduced, the system is closed (no further water is allowed to enter the system) and its evolution is followed

The following input file models this situation

# Simulation to assess natural changes in the reservoir.  Recall that water_60_to_220degC has provided a stable mineral assemblage that is in agreement with XRD observations, and a water at equilibrium with that assemblage.  However, Stuart Simmons suggested including Laumontite and Zoisite into the simulation, and they were not included in water_60_to_220degC since they are more stable than Anorthite, so Anorthite completely dissolves when equilibrium is assumed.  Here, all minerals suggested by Stuart Simmons are added to the system and kinetics are used to determine the time scales of the mineral changes.  The initial water composition is the reservoir water from water_60_to_220degC.
# The initial mole numbers of the kinetic species are chosen to be such that:
# - the mass fractions are: Albite 0.44; Anorthite 0.05; K-feldspar 0.29; Quartz 0.18, Phlgoptite 0.02 and Illite 0.02 with trace amounts of the remaining minerals.  These are similar to that measured in bulk X-ray diffraction results of 10 samples from well 58-32, assuming that "plagioclase feldspar" consists of Albite and Anorthite in the ratio 9:1, and that Biotite is Phlogoptite.  The trace amounts of each mineral are necessary because of the way the kinetics works: precipitation rate is proportional to mineral-species mass, so without any mass, no precipitation is possible.  Precisely:
# - it is assumed that water_60_to_220degC consists of 1 litre of water (there is 1kg of solvent water) and that the porosity is 0.01, so the amount of rock should be 99000cm^3
# - the cm^3 of the trace minerals Calcite and Anhydrite is exactly given by water_60_to_220degC (0.016 and 0.018 respectively)
# - see initial_kinetic_moles.xlsx for the remaining mole numbers
# The results depend on the kinetic rates used and these are recognised to be poorly constrained by experiment
[UserObjects]
  [rate_Albite]
    type = GeochemistryKineticRate
    kinetic_species_name = Albite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 69.8E3
    one_over_T0 = 0.003354
  []
  [rate_Anhydrite]
    type = GeochemistryKineticRate
    kinetic_species_name = Anhydrite
    intrinsic_rate_constant = 1.0E-7
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 14.3E3
    one_over_T0 = 0.003354
  []
  [rate_Anorthite]
    type = GeochemistryKineticRate
    kinetic_species_name = Anorthite
    intrinsic_rate_constant = 1.0E-13
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 17.8E3
    one_over_T0 = 0.003354
  []
  [rate_Calcite]
    type = GeochemistryKineticRate
    kinetic_species_name = Calcite
    intrinsic_rate_constant = 1.0E-10
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 23.5E3
    one_over_T0 = 0.003354
  []
  [rate_Chalcedony]
    type = GeochemistryKineticRate
    kinetic_species_name = Chalcedony
    intrinsic_rate_constant = 1.0E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 90.1E3
    one_over_T0 = 0.003354
  []
  [rate_Clinochl-7A]
    type = GeochemistryKineticRate
    kinetic_species_name = Clinochl-7A
    intrinsic_rate_constant = 1.0E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 88.0E3
    one_over_T0 = 0.003354
  []
  [rate_Illite]
    type = GeochemistryKineticRate
    kinetic_species_name = Illite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 29E3
    one_over_T0 = 0.003354
  []
  [rate_K-feldspar]
    type = GeochemistryKineticRate
    kinetic_species_name = K-feldspar
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 38E3
    one_over_T0 = 0.003354
  []
  [rate_Kaolinite]
    type = GeochemistryKineticRate
    kinetic_species_name = Kaolinite
    intrinsic_rate_constant = 1E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22.2E3
    one_over_T0 = 0.003354
  []
  [rate_Quartz]
    type = GeochemistryKineticRate
    kinetic_species_name = Quartz
    intrinsic_rate_constant = 1E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 90.1E3
    one_over_T0 = 0.003354
  []
  [rate_Paragonite]
    type = GeochemistryKineticRate
    kinetic_species_name = Paragonite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22E3
    one_over_T0 = 0.003354
  []
  [rate_Phlogopite]
    type = GeochemistryKineticRate
    kinetic_species_name = Phlogopite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22E3
    one_over_T0 = 0.003354
  []
  [rate_Laumontite]
    type = GeochemistryKineticRate
    kinetic_species_name = Laumontite
    intrinsic_rate_constant = 1.0E-15
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 17.8E3
    one_over_T0 = 0.003354
  []
  [rate_Zoisite]
    type = GeochemistryKineticRate
    kinetic_species_name = Zoisite
    intrinsic_rate_constant = 1E-16
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 66.1E3
    one_over_T0 = 0.003354
  []
  [definition]
    type = GeochemicalModelDefinition
    database_file = '../../../../geochemistry/database/moose_geochemdb.json'
    basis_species = 'H2O H+ Na+ K+ Ca++ Mg++ SiO2(aq) Al+++ Cl- SO4-- HCO3-'
    remove_all_extrapolated_secondary_species = true
    kinetic_minerals = 'Albite Anhydrite Anorthite Calcite Chalcedony Clinochl-7A Illite K-feldspar Kaolinite Quartz Paragonite Phlogopite Zoisite Laumontite'
    kinetic_rate_descriptions = 'rate_Albite rate_Anhydrite rate_Anorthite rate_Calcite rate_Chalcedony rate_Clinochl-7A rate_Illite rate_K-feldspar rate_Kaolinite rate_Quartz rate_Paragonite rate_Phlogopite rate_Zoisite rate_Laumontite'
  []
[]

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = 'Cl-'
  constraint_species = 'H2O              H+                  Na+                K+                 Ca++               Mg++                SiO2(aq)           Al+++               Cl-              SO4--               HCO3-'
# Following numbers are from water_60_to_220degC_out.csv
  constraint_value = '  1.0006383866109  9.5165072498215e-07 0.100020379171     0.0059389061065    0.011570884507621  4.6626763057447e-06 0.0045110404925255 5.8096968688789e-17 0.13500708594394 6.6523540147676e-05 7.7361407898089e-05'
  constraint_meaning = 'kg_solvent_water free_concentration  free_concentration free_concentration free_concentration free_concentration  free_concentration free_concentration  bulk_composition free_concentration  free_concentration'
  constraint_unit = '   kg               molal               molal              molal              molal              molal               molal              molal               moles            molal               molal'
  initial_temperature = 220
  temperature = 220
  kinetic_species_name = '         Albite             Anorthite          K-feldspar         Quartz             Phlogopite         Paragonite         Calcite            Anhydrite          Chalcedony         Illite             Kaolinite          Clinochl-7A        Zoisite            Laumontite'
  kinetic_species_initial_value = '4.324073236492E+02 4.631370307325E+01 2.685015418378E+02 7.720095013956E+02 1.235192062541E+01 7.545461404965E-01 4.234651808835E-04 4.000485907930E-04 4.407616361072E+00 1.342524904876E+01 1.004823151125E+00 4.728132387707E-01 7.326007326007E-01 4.818116116598E-01'
  kinetic_species_unit = '         moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles'
  evaluate_kinetic_rates_always = true # otherwise will easily "run out" of dissolving species
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  execute_console_output_on = ''
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = 'max(1E2, 0.1 * t)'
  []
  end_time = 4E11
[]

[GlobalParams]
  point = '0 0 0'
[]

[Postprocessors]
  [cm3_Albite]
    type = PointValue
    variable = 'free_cm3_Albite'
  []
  [cm3_Anhydrite]
    type = PointValue
    variable = 'free_cm3_Anhydrite'
  []
  [cm3_Anorthite]
    type = PointValue
    variable = 'free_cm3_Anorthite'
  []
  [cm3_Calcite]
    type = PointValue
    variable = 'free_cm3_Calcite'
  []
  [cm3_Chalcedony]
    type = PointValue
    variable = 'free_cm3_Chalcedony'
  []
  [cm3_Clinochl-7A]
    type = PointValue
    variable = 'free_cm3_Clinochl-7A'
  []
  [cm3_Illite]
    type = PointValue
    variable = 'free_cm3_Illite'
  []
  [cm3_K-feldspar]
    type = PointValue
    variable = 'free_cm3_K-feldspar'
  []
  [cm3_Kaolinite]
    type = PointValue
    variable = 'free_cm3_Kaolinite'
  []
  [cm3_Quartz]
    type = PointValue
    variable = 'free_cm3_Quartz'
  []
  [cm3_Paragonite]
    type = PointValue
    variable = 'free_cm3_Paragonite'
  []
  [cm3_Phlogopite]
    type = PointValue
    variable = 'free_cm3_Phlogopite'
  []
  [cm3_Zoisite]
    type = PointValue
    variable = 'free_cm3_Zoisite'
  []
  [cm3_Laumontite]
    type = PointValue
    variable = 'free_cm3_Laumontite'
  []
  [cm3_mineral]
    type = LinearCombinationPostprocessor
    pp_names = 'cm3_Albite cm3_Anhydrite cm3_Anorthite cm3_Calcite cm3_Chalcedony cm3_Clinochl-7A cm3_Illite cm3_K-feldspar cm3_Kaolinite cm3_Quartz cm3_Paragonite cm3_Phlogopite cm3_Zoisite cm3_Laumontite'
    pp_coefs = '1 1 1 1 1 1 1 1 1 1 1 1 1 1'
  []
  [pH]
    type = PointValue
    variable = 'pH'
  []
  [kg_solvent_H2O]
    type = PointValue
    variable = 'kg_solvent_H2O'
  []
[]

[Outputs]
  csv = true
[]

Results are shown in Figure 1, Figure 2 and Figure 3 The figures reveal that:

• The small amount of Calcite dissolves entirely within years

• About 25% of the small amount of Anhydrite dissolves rapidly, but after 100 years it returns to its original volume

• Between 1 and 100 years after the simulation commences, all Anorthite, Chalcedony and Clinochl-7A dissolve

• Between 1 and 100 years after the simulation commences, Zoisite, Illite, Paragony and Quartz increase in volume. Laumontite increases slightly too, before returning to its original small value.

• The pH drops to around 4.8 during the course of the first 100 years before rebounding to 5.8.

• The solvent water mass roughly halves during the first 100 years, corresponding to small porosity decreases.

The results depend quite crucially on the kinetic rate constants, which are poorly constrained by experiment. For instance, increasing Laumontite's rate constant by a factor of 10 means that all the water gets absorbed into minerals as they precipitate, resulting in geothermal reservoir with zero porosity.

This concludes the study on the in-situ dynamics of the reservoir. Attention is now turned to the injectate and, finally, reactive-transport modelling.

Figure 1: Change in mineral volume when the in-situ reservoir water is in contact with the minerals of Table 3.

Figure 2: Percentage change in mineral volume when the in-situ reservoir water is in contact with the minerals of Table 3.

Figure 3: Changes in the water volume and chemistry when the in-situ reservoir water is in contact with the minerals of Table 3.

Water3 at 70C

The reactive-transport simulation involves injecting Water3 at 70C and it is useful to determine the free molalities of species in this water. To find these:

1. The initial equilibrium is found at 20degC. This is the temperature at which the bulk composition was measured, and at this temperature most species are supersaturated. However, since measurements were performed in the absence of free minerals, their precipitation must be retarded in some way, so all minerals are prevented from precipitating in the model

2. The pH constraint is removed and the system is raised to 70degC, which is the injection temperature. This causes the pH to drop from 6.2 to 6.1, and only Kaolinite and Illite are supersaturated

An input file that performs this simulation is

# Equilibrium model "Water 3" from "Subtask 2C.4.7 Geochemical Modeling SSimmons-VPatil.pdf".  The steps followed in this input file are:
# 1. The initial equilibrium is found at 20degC.  This is the temperature at which the bulk composition was measured, and at this temperature most species are supersaturated.  However, since measurements were performed in the absence of free minerals, their precipitation must be retarded in some way, so all minerals are prevented from precipitating in the model
# 2. The pH constraint is removed and the system is raised to 70degC, which is the injection temperature.  This causes the pH to drop from 6.2 to 6.1, and only Kaolinite and Illite are supersaturated
# 3. The free molality of the species is measured for use in other models
[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = '../../../../geochemistry/database/moose_geochemdb.json'
    basis_species = 'H2O H+ Na+ K+ Ca++ Mg++ SiO2(aq) Al+++ Cl- SO4-- HCO3-'
    equilibrium_minerals = 'Albite Anhydrite Anorthite Calcite Chalcedony Clinochl-7A Illite K-feldspar Kaolinite Quartz Paragonite Phlogopite Zoisite Laumontite'
    remove_all_extrapolated_secondary_species = true
  []
[]

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = 'Cl-'
  constraint_species = 'H2O H+      Na+     K+      Ca++    Mg++    SiO2(aq) Al+++   Cl-     SO4--   HCO3-'
  constraint_value = '  1.0 6.31E-7 1.32E-4 2.81E-5 7.76E-5 2.88E-5 2.73E-4  3.71E-6 1.41E-5 1.04E-5 3.28E-4'
  constraint_meaning = 'kg_solvent_water activity 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'
  prevent_precipitation = 'Albite Anhydrite Anorthite Calcite Chalcedony Clinochl-7A Illite K-feldspar Kaolinite Quartz Paragonite Phlogopite Zoisite Laumontite'
  initial_temperature = 20
  remove_fixed_activity_name = 'H+'
  remove_fixed_activity_time = 0
  temperature = 70
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  execute_console_output_on = 'final' # only CSV output needed
[]

[Executioner]
  type = Transient
  dt = 1
  end_time = 1
[]

[GlobalParams]
  point = '0 0 0'
  reactor = reactor
[]

[AuxVariables]
  [transported_H2O]
  []
  [transported_H+]
  []
  [transported_Na+]
  []
  [transported_K+]
  []
  [transported_Ca++]
  []
  [transported_Mg++]
  []
  [transported_SiO2]
  []
  [transported_Al+++]
  []
  [transported_Cl-]
  []
  [transported_SO4--]
  []
  [transported_HCO3-]
  []
[]

[AuxKernels]
  [transported_H2O]
    type = GeochemistryQuantityAux
    species = 'H2O'
    variable = transported_H2O
    quantity = transported_moles_in_original_basis
  []
  [transported_H+]
    type = GeochemistryQuantityAux
    species = 'H+'
    variable = transported_H+
    quantity = transported_moles_in_original_basis
  []
  [transported_Na+]
    type = GeochemistryQuantityAux
    species = 'Na+'
    variable = transported_Na+
    quantity = transported_moles_in_original_basis
  []
  [transported_K+]
    type = GeochemistryQuantityAux
    species = 'K+'
    variable = transported_K+
    quantity = transported_moles_in_original_basis
  []
  [transported_Ca++]
    type = GeochemistryQuantityAux
    species = 'Ca++'
    variable = transported_Ca++
    quantity = transported_moles_in_original_basis
  []
  [transported_Mg++]
    type = GeochemistryQuantityAux
    species = 'Mg++'
    variable = transported_Mg++
    quantity = transported_moles_in_original_basis
  []
  [transported_SiO2]
    type = GeochemistryQuantityAux
    species = 'SiO2(aq)'
    variable = transported_SiO2
    quantity = transported_moles_in_original_basis
  []
  [transported_Al+++]
    type = GeochemistryQuantityAux
    species = 'Al+++'
    variable = transported_Al+++
    quantity = transported_moles_in_original_basis
  []
  [transported_Cl-]
    type = GeochemistryQuantityAux
    species = 'Cl-'
    variable = transported_Cl-
    quantity = transported_moles_in_original_basis
  []
  [transported_SO4--]
    type = GeochemistryQuantityAux
    species = 'SO4--'
    variable = transported_SO4--
    quantity = transported_moles_in_original_basis
  []
  [transported_HCO3-]
    type = GeochemistryQuantityAux
    species = 'HCO3-'
    variable = transported_HCO3-
    quantity = transported_moles_in_original_basis
  []
[]

[Postprocessors]
  [temperature]
    type = PointValue
    variable = 'solution_temperature'
  []
  [bulk_Cl]
    type = PointValue
    variable = 'bulk_moles_Cl-'
  []
  [kg_solvent_H2O]
    type = PointValue
    variable = 'kg_solvent_H2O'
  []
  [molal_H+]
    type = PointValue
    variable = 'molal_H+'
  []
  [molal_Na+]
    type = PointValue
    variable = 'molal_Na+'
  []
  [molal_K+]
    type = PointValue
    variable = 'molal_K+'
  []
  [molal_Ca++]
    type = PointValue
    variable = 'molal_Ca++'
  []
  [molal_Mg++]
    type = PointValue
    variable = 'molal_Mg++'
  []
  [molal_SiO2aq]
    type = PointValue
    variable = 'molal_SiO2(aq)'
  []
  [molal_Al+++]
    type = PointValue
    variable = 'molal_Al+++'
  []
  [molal_Cl-]
    type = PointValue
    variable = 'molal_Cl-'
  []
  [molal_SO4--]
    type = PointValue
    variable = 'molal_SO4--'
  []
  [molal_HCO3-]
    type = PointValue
    variable = 'molal_HCO3-'
  []
  [transported_H2O]
    type = PointValue
    variable = transported_H2O
  []
  [transported_H+]
    type = PointValue
    variable = transported_H+
  []
  [transported_Na+]
    type = PointValue
    variable = transported_Na+
  []
  [transported_K+]
    type = PointValue
    variable = transported_K+
  []
  [transported_Ca++]
    type = PointValue
    variable = transported_Ca++
  []
  [transported_Mg++]
    type = PointValue
    variable = transported_Mg++
  []
  [transported_SiO2]
    type = PointValue
    variable = transported_SiO2
  []
  [transported_Al+++]
    type = PointValue
    variable = transported_Al+++
  []
  [transported_Cl-]
    type = PointValue
    variable = transported_Cl-
  []
  [transported_SO4--]
    type = PointValue
    variable = transported_SO4--
  []
  [transported_HCO3-]
    type = PointValue
    variable = transported_HCO3-
  []
[]

[Outputs]
  csv = true
[]

Note that in the TimeDependentReactionSolver:

• The pH is initially set to 6.2 by setting the activity of H+ to

• All minerals are prevented from precipitating

• The initial temperature is 20C

• After the system has equilibrated during the initial setup phase, the activity constraint on H+ is removed at

• The temperature after the initial setup phase is set to 70C

Injectate impacts on the reservoir mineralogy

Water3 at 70C (from the previous section) is mixed with rock with composition of Table 3 in order to assess the possible changes in the reservoir. This is not a reactive transport simulation: it is simply flushing Water3 repeatedly through the rock mineral assemblage. That is, as soon as a mineral dissolves, the aqueous components are swept away and replaced by a new batch of Water3. Or, as soon as a precipitate forms from Water3, it is retained and a new batch of Water3 is introduced, ready to potentially form more precipitate. The rate constants of Table 5 are used.

An input file that performs this simulation is

# Simulation to assess possible changes in the reservoir.  The rock composition from natural_reservoir.i is mixed with the water from water_3.i  Note that the free_concentration values are used from water_3.i and that composition is held fixed throughout this entire simulation.  This models water_3 continually flushing through the rock mineral assemblage: as soon as a mineral dissolves the aqueous components are swept away and replaced by a new batch of water_3; as soon as mineral precipitates more water_3 sweeps into the system providing a limitless source of aqueous components (in set ratios) at 70degC
# The results depend on the kinetic rates used and these are recognised to be poorly constrained by experiment
[UserObjects]
  [rate_Albite]
    type = GeochemistryKineticRate
    kinetic_species_name = Albite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 69.8E3
    one_over_T0 = 0.003354
  []
  [rate_Anhydrite]
    type = GeochemistryKineticRate
    kinetic_species_name = Anhydrite
    intrinsic_rate_constant = 1.0E-7
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 14.3E3
    one_over_T0 = 0.003354
  []
  [rate_Anorthite]
    type = GeochemistryKineticRate
    kinetic_species_name = Anorthite
    intrinsic_rate_constant = 1.0E-13
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 17.8E3
    one_over_T0 = 0.003354
  []
  [rate_Calcite]
    type = GeochemistryKineticRate
    kinetic_species_name = Calcite
    intrinsic_rate_constant = 1.0E-10
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 23.5E3
    one_over_T0 = 0.003354
  []
  [rate_Chalcedony]
    type = GeochemistryKineticRate
    kinetic_species_name = Chalcedony
    intrinsic_rate_constant = 1.0E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 90.1E3
    one_over_T0 = 0.003354
  []
  [rate_Clinochl-7A]
    type = GeochemistryKineticRate
    kinetic_species_name = Clinochl-7A
    intrinsic_rate_constant = 1.0E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 88.0E3
    one_over_T0 = 0.003354
  []
  [rate_Illite]
    type = GeochemistryKineticRate
    kinetic_species_name = Illite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 29E3
    one_over_T0 = 0.003354
  []
  [rate_K-feldspar]
    type = GeochemistryKineticRate
    kinetic_species_name = K-feldspar
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 38E3
    one_over_T0 = 0.003354
  []
  [rate_Kaolinite]
    type = GeochemistryKineticRate
    kinetic_species_name = Kaolinite
    intrinsic_rate_constant = 1E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22.2E3
    one_over_T0 = 0.003354
  []
  [rate_Quartz]
    type = GeochemistryKineticRate
    kinetic_species_name = Quartz
    intrinsic_rate_constant = 1E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 90.1E3
    one_over_T0 = 0.003354
  []
  [rate_Paragonite]
    type = GeochemistryKineticRate
    kinetic_species_name = Paragonite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22E3
    one_over_T0 = 0.003354
  []
  [rate_Phlogopite]
    type = GeochemistryKineticRate
    kinetic_species_name = Phlogopite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22E3
    one_over_T0 = 0.003354
  []
  [rate_Laumontite]
    type = GeochemistryKineticRate
    kinetic_species_name = Laumontite
    intrinsic_rate_constant = 1.0E-15
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 17.8E3
    one_over_T0 = 0.003354
  []
  [rate_Zoisite]
    type = GeochemistryKineticRate
    kinetic_species_name = Zoisite
    intrinsic_rate_constant = 1E-16
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 66.1E3
    one_over_T0 = 0.003354
  []
  [definition]
    type = GeochemicalModelDefinition
    database_file = '../../../../geochemistry/database/moose_geochemdb.json'
    basis_species = 'H2O H+ Na+ K+ Ca++ Mg++ SiO2(aq) Al+++ Cl- SO4-- HCO3-'
    remove_all_extrapolated_secondary_species = true
    kinetic_minerals = 'Albite Anhydrite Anorthite Calcite Chalcedony Clinochl-7A Illite K-feldspar Kaolinite Quartz Paragonite Phlogopite Zoisite Laumontite'
    kinetic_rate_descriptions = 'rate_Albite rate_Anhydrite rate_Anorthite rate_Calcite rate_Chalcedony rate_Clinochl-7A rate_Illite rate_K-feldspar rate_Kaolinite rate_Quartz rate_Paragonite rate_Phlogopite rate_Zoisite rate_Laumontite'
  []
[]

[TimeDependentReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = 'Cl-'
  constraint_species = 'H2O              H+                  Na+                K+                  Ca++                Mg++               SiO2(aq)            Al+++               Cl-                SO4--               HCO3-'
# Following numbers are from water_3_out.csv
  constraint_value = '  0.99999999549877 8.0204734722945e-07 0.0001319920398478 2.8097346859027e-05 7.7328020546464e-05 2.874602030221e-05 0.00027284654762868 4.4715524787497e-12 0.0002253530818877 1.0385772502298e-05 0.00012427759434288'
  constraint_meaning = 'kg_solvent_water free_concentration       free_concentration    free_concentration      free_concentration     free_concentration       free_concentration      free_concentration       bulk_composition free_concentration       free_concentration'
  constraint_unit = '   kg               molal               molal            molal              molal             molal               molal              molal               moles              molal               molal'
  initial_temperature = 70
  temperature = 70
  close_system_at_time = 1E20 # keep the free molalities specified above for all time
  kinetic_species_name = '         Albite             Anorthite          K-feldspar         Quartz             Phlogopite         Paragonite         Calcite            Anhydrite          Chalcedony         Illite             Kaolinite          Clinochl-7A        Zoisite            Laumontite'
  kinetic_species_initial_value = '4.324073236492E+02 4.631370307325E+01 2.685015418378E+02 7.720095013956E+02 1.235192062541E+01 7.545461404965E-01 4.234651808835E-04 4.000485907930E-04 4.407616361072E+00 1.342524904876E+01 1.004823151125E+00 4.728132387707E-01 7.326007326007E-01 4.818116116598E-01'
  kinetic_species_unit = '         moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles'
  evaluate_kinetic_rates_always = true # otherwise will easily "run out" of dissolving species
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  execute_console_output_on = ''
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = 'max(1E6, 0.3 * t)'
  []
  end_time = 4E11
[]

[GlobalParams]
  point = '0 0 0'
[]

[Postprocessors]
  [temperature]
    type = PointValue
    variable = 'solution_temperature'
  []
  [cm3_Albite]
    type = PointValue
    variable = 'free_cm3_Albite'
  []
  [cm3_Anhydrite]
    type = PointValue
    variable = 'free_cm3_Anhydrite'
  []
  [cm3_Anorthite]
    type = PointValue
    variable = 'free_cm3_Anorthite'
  []
  [cm3_Calcite]
    type = PointValue
    variable = 'free_cm3_Calcite'
  []
  [cm3_Chalcedony]
    type = PointValue
    variable = 'free_cm3_Chalcedony'
  []
  [cm3_Clinochl-7A]
    type = PointValue
    variable = 'free_cm3_Clinochl-7A'
  []
  [cm3_Illite]
    type = PointValue
    variable = 'free_cm3_Illite'
  []
  [cm3_K-feldspar]
    type = PointValue
    variable = 'free_cm3_K-feldspar'
  []
  [cm3_Kaolinite]
    type = PointValue
    variable = 'free_cm3_Kaolinite'
  []
  [cm3_Quartz]
    type = PointValue
    variable = 'free_cm3_Quartz'
  []
  [cm3_Paragonite]
    type = PointValue
    variable = 'free_cm3_Paragonite'
  []
  [cm3_Phlogopite]
    type = PointValue
    variable = 'free_cm3_Phlogopite'
  []
  [cm3_Zoisite]
    type = PointValue
    variable = 'free_cm3_Zoisite'
  []
  [cm3_Laumontite]
    type = PointValue
    variable = 'free_cm3_Laumontite'
  []
  [cm3_mineral]
    type = LinearCombinationPostprocessor
    pp_names = 'cm3_Albite cm3_Anhydrite cm3_Anorthite cm3_Calcite cm3_Chalcedony cm3_Clinochl-7A cm3_Illite cm3_K-feldspar cm3_Kaolinite cm3_Quartz cm3_Paragonite cm3_Phlogopite cm3_Zoisite cm3_Laumontite'
    pp_coefs = '1 1 1 1 1 1 1 1 1 1 1 1 1 1'
  []
  [pH]
    type = PointValue
    variable = 'pH'
  []
[]

[Outputs]
  csv = true
[]

The results are shown in Figure 4 and Figure 5, which should be explored keeping the baseline of Figure 1 and Figure 2 in mind.

• The cool injectate causes complete dissolution of the small quantities of Anhydrite

• The large quantities of Anorthite and small quantities of Calcite completely dissolve, but this could also be due to natural reservoir processes

• After 100 years, the cool injectate causes gradual and complete dissolution of the small quantities of Zoisite and Laumontite

• After 100 years, the small quantities of Clinochl-7A gradually dissolve, but this could also be due to natural reservoir processes

• After 100 years, large quantities of Albite begin to dissolve

• After 1000 years, quantities of K-feldspar begin to dissolve

• Small quantities of Kaolinite precipitate, but no other minerals precipitate

The key results are:

• Flushing with Water3 at 70C prevents any appreciable precipitation

• Hence, the dissolution of Anorthite commencing after about 1 year, and the dissolution of Albite and K-feldspar later on lead to significant loss of mineral, and a consequential increase in porosity or fracture aperture

Once again, the results depend quite crucially on the kinetic rate constants, which are poorly constrained by experiment.

Figure 4: Change in mineral volume when Water3 at 70C is flushed through the mineral assemblage of Table 3.

Figure 5: Percentage change in mineral volume when Water3 at 70C is flushed through the mineral assemblage of Table 3.

Injection simulation without geochemistry

The simulation involves injecting cold fluid into a thin, hot subsurface confined aquifer. The PorousFlow module is used to simulate the transport of solutes and temperature. The model parameters are displayed in Table 6.

PorousFlow operates using mass-fractions and not mole numbers, so the species concentrations in the reservoir Table 4 and injectate Table 1 must be converted, as shown in Table 7.

The input file that simulates this situation is (ignore the MultiApps and Transfers for the time being):

# Input file modified from RobPodgorney version
# - 2D instead of 3D with different resolution.  Effectively this means a 1m height of RobPodgorney aquifer is simulated.  RobPodgorney total mass flux is 2.5kg/s meaning 0.25kg/s is appropriate here
# - Celsius instead of Kelvin
# - no use of PorousFlowPointEnthalpySourceFromPostprocessor since that is not yet merged into MOOSE: a DirichletBC is used instead
# - Use of PorousFlowFullySaturated instead of PorousFlowUnsaturated, and the save_component_rate_in feature to record the change in kg of each species at each node for passing to the Geochem simulation
# - MultiApps and Transfers to transfer information between this simulation and the aquifer_geochemistry.i simulation
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 225
    ny = 200
    xmin = -400
    xmax = 500
    ymin = -400
    ymax = 400
  []
  [injection_node]
    input = gen
    type = ExtraNodesetGenerator
    new_boundary = injection_node
    coord = '0 0 0'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [f_H]
    initial_condition = 8.201229858451E-07
  []
  [f_Na]
    initial_condition = 2.281094143525E-03
  []
  [f_K]
    initial_condition = 2.305489507836E-04
  []
  [f_Ca]
    initial_condition = 5.818776782059E-04
  []
  [f_Mg]
    initial_condition = 1.539513498238E-07
  []
  [f_SiO2]
    initial_condition = 2.691822196469E-04
  []
  [f_Al]
    initial_condition = 4.457519474122E-08
  []
  [f_Cl]
    initial_condition = 4.744309776594E-03
  []
  [f_SO4]
    initial_condition = 9.516650880811E-06
  []
  [f_HCO3]
    initial_condition = 5.906126982324E-05
  []
  [porepressure]
    initial_condition = 20E6
  []
  [temperature]
    initial_condition = 220 # degC
    scaling = 1E-6 # fluid enthalpy is roughly 1E6
  []
[]

[BCs]
  [source_temperature]
    type = DirichletBC
    boundary = injection_node
    variable = temperature
    value = 70 # degC
  []
[]

[DiracKernels]
  [inject_H]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 4.790385871045E-08
    variable = f_H
  []
  [inject_Na]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 7.586252963780E-07
    variable = f_Na
  []
  [inject_K]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.746517625125E-07
    variable = f_K
  []
  [inject_Ca]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 7.775129478597E-07
    variable = f_Ca
  []
  [inject_Mg]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 1.749872109005E-07
    variable = f_Mg
  []
  [inject_SiO2]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 4.100547515915E-06
    variable = f_SiO2
  []
  [inject_Al]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.502408592080E-08
    variable = f_Al
  []
  [inject_Cl]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 1.997260386272E-06
    variable = f_Cl
  []
  [inject_SO4]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.497372164191E-07
    variable = f_SO4
  []
  [inject_HCO3]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 5.003150992902E-06
    variable = f_HCO3
  []
  [inject_H2O]
    type = PorousFlowPointSourceFromPostprocessor
    point = ' 0 0 0'
    mass_flux = 2.499865905987E-01
    variable = porepressure
  []

  [produce_H]
    type = PorousFlowPeacemanBorehole
    variable = f_H
    SumQuantityUO = produced_mass_H
    mass_fraction_component = 0
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Na]
    type = PorousFlowPeacemanBorehole
    variable = f_Na
    SumQuantityUO = produced_mass_Na
    mass_fraction_component = 1
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_K]
    type = PorousFlowPeacemanBorehole
    variable = f_K
    SumQuantityUO = produced_mass_K
    mass_fraction_component = 2
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Ca]
    type = PorousFlowPeacemanBorehole
    variable = f_Ca
    SumQuantityUO = produced_mass_Ca
    mass_fraction_component = 3
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Mg]
    type = PorousFlowPeacemanBorehole
    variable = f_Mg
    SumQuantityUO = produced_mass_Mg
    mass_fraction_component = 4
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_SiO2]
    type = PorousFlowPeacemanBorehole
    variable = f_SiO2
    SumQuantityUO = produced_mass_SiO2
    mass_fraction_component = 5
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Al]
    type = PorousFlowPeacemanBorehole
    variable = f_Al
    SumQuantityUO = produced_mass_Al
    mass_fraction_component = 6
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_Cl]
    type = PorousFlowPeacemanBorehole
    variable = f_Cl
    SumQuantityUO = produced_mass_Cl
    mass_fraction_component = 7
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_SO4]
    type = PorousFlowPeacemanBorehole
    variable = f_SO4
    SumQuantityUO = produced_mass_SO4
    mass_fraction_component = 8
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_HCO3]
    type = PorousFlowPeacemanBorehole
    variable = f_HCO3
    SumQuantityUO = produced_mass_HCO3
    mass_fraction_component = 9
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [produce_H2O]
    type = PorousFlowPeacemanBorehole
    variable = porepressure
    SumQuantityUO = produced_mass_H2O
    mass_fraction_component = 10
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    character = 1
  []
  [remove_heat_at_production_well]
    type = PorousFlowPeacemanBorehole
    variable = temperature
    SumQuantityUO = produced_heat
    point_file = production.bh
    line_length = 1
    bottom_p_or_t = 20E6
    unit_weight = '0 0 0'
    use_mobility = true
    use_enthalpy = true
    character = 1
  []
[]

[UserObjects]
  [produced_mass_H]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Na]
    type = PorousFlowSumQuantity
  []
  [produced_mass_K]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Ca]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Mg]
    type = PorousFlowSumQuantity
  []
  [produced_mass_SiO2]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Al]
    type = PorousFlowSumQuantity
  []
  [produced_mass_Cl]
    type = PorousFlowSumQuantity
  []
  [produced_mass_SO4]
    type = PorousFlowSumQuantity
  []
  [produced_mass_HCO3]
    type = PorousFlowSumQuantity
  []
  [produced_mass_H2O]
    type = PorousFlowSumQuantity
  []
  [produced_heat]
    type = PorousFlowSumQuantity
  []
[]

[Postprocessors]
  [heat_extracted]
    type = PorousFlowPlotQuantity
    uo = produced_heat
  []
  [approx_production_temperature]
    type = PointValue
    point = '100 0 0'
    variable = temperature
  []
  [mass_extracted_H]
    type = PorousFlowPlotQuantity
    uo = produced_mass_H
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Na]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Na
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_K]
    type = PorousFlowPlotQuantity
    uo = produced_mass_K
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Ca]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Ca
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Mg]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Mg
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_SiO2]
    type = PorousFlowPlotQuantity
    uo = produced_mass_SiO2
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Al]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Al
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_Cl]
    type = PorousFlowPlotQuantity
    uo = produced_mass_Cl
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_SO4]
    type = PorousFlowPlotQuantity
    uo = produced_mass_SO4
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_HCO3]
    type = PorousFlowPlotQuantity
    uo = produced_mass_HCO3
    execute_on = 'initial timestep_end'
  []
  [mass_extracted_H2O]
    type = PorousFlowPlotQuantity
    uo = produced_mass_H2O
    execute_on = 'initial timestep_end'
  []
  [mass_extracted]
    type = LinearCombinationPostprocessor
    pp_names = 'mass_extracted_H mass_extracted_Na mass_extracted_K mass_extracted_Ca mass_extracted_Mg mass_extracted_SiO2 mass_extracted_Al mass_extracted_Cl mass_extracted_SO4 mass_extracted_HCO3 mass_extracted_H2O'
    pp_coefs = '1 1 1 1 1 1 1 1 1 1 1'
    execute_on = 'initial timestep_end'
  []
  [dt]
    type = TimestepSize
    execute_on = 'timestep_begin'
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    thermal_expansion = 2E-4
    bulk_modulus = 2E9
    viscosity = 1E-3
    density0 = 980
    cv = 4000.0
    cp = 4000.0
    porepressure_coefficient = 0
  []
[]

[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = temperature
  mass_fraction_vars = 'f_H f_Na f_K f_Ca f_Mg f_SiO2 f_Al f_Cl f_SO4 f_HCO3'
  save_component_rate_in = 'rate_H rate_Na rate_K rate_Ca rate_Mg rate_SiO2 rate_Al rate_Cl rate_SO4 rate_HCO3 rate_H2O' # change in kg at every node / dt
  fp = the_simple_fluid
  temperature_unit = Celsius
[]

[AuxVariables]
  [rate_H]
  []
  [rate_Na]
  []
  [rate_K]
  []
  [rate_Ca]
  []
  [rate_Mg]
  []
  [rate_SiO2]
  []
  [rate_Al]
  []
  [rate_Cl]
  []
  [rate_SO4]
  []
  [rate_HCO3]
  []
  [rate_H2O]
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.01
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2.5 0 0  0 2.5 0  0 0 2.5'
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    density = 2750.0
    specific_heat_capacity = 900.0
  []
[]

[Preconditioning]
  active = typically_efficient
  [typically_efficient]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_hypre_type'
    petsc_options_value = ' hypre    boomeramg'
  []
  [strong]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      ilu           NONZERO                   2'
  []
  [probably_too_strong]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 31536000 #1 year

  [TimeStepper]
    type = SolutionTimeAdaptiveDT
    dt = 500
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[MultiApps]
  [react]
    type = TransientMultiApp
    input_files = aquifer_geochemistry.i
    clone_master_mesh = true
    execute_on = 'timestep_end'
  []
[]
[Transfers]
  [changes_due_to_flow]
    type = MultiAppCopyTransfer
    source_variable = 'rate_H rate_Na rate_K rate_Ca rate_Mg rate_SiO2 rate_Al rate_Cl rate_SO4 rate_HCO3 rate_H2O temperature'
    variable = 'pf_rate_H pf_rate_Na pf_rate_K pf_rate_Ca pf_rate_Mg pf_rate_SiO2 pf_rate_Al pf_rate_Cl pf_rate_SO4 pf_rate_HCO3 pf_rate_H2O temperature'
    to_multi_app = react
  []
  [massfrac_from_geochem]
    type = MultiAppCopyTransfer
    source_variable = 'massfrac_H massfrac_Na massfrac_K massfrac_Ca massfrac_Mg massfrac_SiO2 massfrac_Al massfrac_Cl massfrac_SO4 massfrac_HCO3'
    variable = 'f_H f_Na f_K f_Ca f_Mg f_SiO2 f_Al f_Cl f_SO4 f_HCO3'
    from_multi_app = react
  []
[]

Some results concerning the produced fluid are shown in Figure 6, Figure 7 and Figure 8. These figures demonstrate that:

• apart at the very beginning of the simulation, the production rate equals the injection rate;

• fluid break-through occurs quite rapidly — certainly in less than a month after injection commences;

• thermal break-through occurs after about 6 months of injection.

Figure 6: Production rate.

Figure 7: Production fraction of SiO2(aq).

Figure 8: Production temperature.

Some contour plots are shown in Figure 9 and Figure 10. The former demonstrates the injectate fluid-breakthrough after only 1 week of injecting.

Figure 9: Mass-fraction of SiO2(aq) contour after only 1 week of injection. The green dots show the injection and production wells.

Figure 10: Temperature contour after 1 year of injection. The green dots show the injection and production wells.

Reactive transport simulation

The sections above have prepared the ground for a reactive transport simulation. The geochemistry part of the reactive-transport simulation is:

# Simulates geochemistry in the aquifer.  This input file may be run in standalone fashion, which will study the natural kinetically-controlled mineral changes in the same way as natural_reservoir.i.  To simulate the FORGE injection scenario, run the porous_flow.i simulation which couples to this input file using MultiApps.
# This file receives pf_rate_H pf_rate_Na pf_rate_K pf_rate_Ca pf_rate_Mg pf_rate_SiO2 pf_rate_Al pf_rate_Cl pf_rate_SO4 pf_rate_HCO3 pf_rate_H2O and temperature as AuxVariables from porous_flow.i
# The pf_rate quantities are kg/s changes of fluid-component mass at each node, but the geochemistry module expects rates-of-changes of moles at every node.  Secondly, since this input file considers just 1 litre of aqueous solution at every node, the nodal_void_volume is used to convert pf_rate_* into rate_*_per_1l, which is measured in mol/s/1_litre_of_aqueous_solution.
# This file sends massfrac_H massfrac_Na massfrac_K massfrac_Ca massfrac_Mg massfrac_SiO2 massfrac_Al massfrac_Cl massfrac_SO4 massfrac_HCO3 to porous_flow.i.  These are computed from the corresponding transported_* quantities.
# The results depend on the kinetic rates used and these are recognised to be poorly constrained by experiment
[UserObjects]
  [rate_Albite]
    type = GeochemistryKineticRate
    kinetic_species_name = Albite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 69.8E3
    one_over_T0 = 0.003354
  []
  [rate_Anhydrite]
    type = GeochemistryKineticRate
    kinetic_species_name = Anhydrite
    intrinsic_rate_constant = 1.0E-7
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 14.3E3
    one_over_T0 = 0.003354
  []
  [rate_Anorthite]
    type = GeochemistryKineticRate
    kinetic_species_name = Anorthite
    intrinsic_rate_constant = 1.0E-13
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 17.8E3
    one_over_T0 = 0.003354
  []
  [rate_Calcite]
    type = GeochemistryKineticRate
    kinetic_species_name = Calcite
    intrinsic_rate_constant = 1.0E-10
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 23.5E3
    one_over_T0 = 0.003354
  []
  [rate_Chalcedony]
    type = GeochemistryKineticRate
    kinetic_species_name = Chalcedony
    intrinsic_rate_constant = 1.0E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 90.1E3
    one_over_T0 = 0.003354
  []
  [rate_Clinochl-7A]
    type = GeochemistryKineticRate
    kinetic_species_name = Clinochl-7A
    intrinsic_rate_constant = 1.0E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 88.0E3
    one_over_T0 = 0.003354
  []
  [rate_Illite]
    type = GeochemistryKineticRate
    kinetic_species_name = Illite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 29E3
    one_over_T0 = 0.003354
  []
  [rate_K-feldspar]
    type = GeochemistryKineticRate
    kinetic_species_name = K-feldspar
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 38E3
    one_over_T0 = 0.003354
  []
  [rate_Kaolinite]
    type = GeochemistryKineticRate
    kinetic_species_name = Kaolinite
    intrinsic_rate_constant = 1E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22.2E3
    one_over_T0 = 0.003354
  []
  [rate_Quartz]
    type = GeochemistryKineticRate
    kinetic_species_name = Quartz
    intrinsic_rate_constant = 1E-18
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 90.1E3
    one_over_T0 = 0.003354
  []
  [rate_Paragonite]
    type = GeochemistryKineticRate
    kinetic_species_name = Paragonite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22E3
    one_over_T0 = 0.003354
  []
  [rate_Phlogopite]
    type = GeochemistryKineticRate
    kinetic_species_name = Phlogopite
    intrinsic_rate_constant = 1E-17
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 22E3
    one_over_T0 = 0.003354
  []
  [rate_Laumontite]
    type = GeochemistryKineticRate
    kinetic_species_name = Laumontite
    intrinsic_rate_constant = 1.0E-15
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 17.8E3
    one_over_T0 = 0.003354
  []
  [rate_Zoisite]
    type = GeochemistryKineticRate
    kinetic_species_name = Zoisite
    intrinsic_rate_constant = 1E-16
    multiply_by_mass = true
    area_quantity = 10
    activation_energy = 66.1E3
    one_over_T0 = 0.003354
  []
  [definition]
    type = GeochemicalModelDefinition
    database_file = '../../../../geochemistry/database/moose_geochemdb.json'
    basis_species = 'H2O H+ Na+ K+ Ca++ Mg++ SiO2(aq) Al+++ Cl- SO4-- HCO3-'
    remove_all_extrapolated_secondary_species = true
    kinetic_minerals = 'Albite Anhydrite Anorthite Calcite Chalcedony Clinochl-7A Illite K-feldspar Kaolinite Quartz Paragonite Phlogopite Zoisite Laumontite'
    kinetic_rate_descriptions = 'rate_Albite rate_Anhydrite rate_Anorthite rate_Calcite rate_Chalcedony rate_Clinochl-7A rate_Illite rate_K-feldspar rate_Kaolinite rate_Quartz rate_Paragonite rate_Phlogopite rate_Zoisite rate_Laumontite'
  []
  [nodal_void_volume_uo]
    type = NodalVoidVolume
    porosity = porosity
    execute_on = 'initial timestep_end' # "initial" means this is evaluated properly for the first timestep
  []
[]

[SpatialReactionSolver]
  model_definition = definition
  geochemistry_reactor_name = reactor
  charge_balance_species = 'Cl-'
  constraint_species = 'H2O              H+                  Na+              K+                 Ca++              Mg++                SiO2(aq)           Al+++               Cl-                SO4--               HCO3-'
# Following numbers are from water_60_to_220degC_out.csv
  constraint_value = '  1.0006383866109  9.5165072498215e-07 0.100020379171   0.0059389061065    0.011570884507621 4.6626763057447e-06 0.0045110404925255 5.8096968688789e-17 0.13500708594394   6.6523540147676e-05 7.7361407898089e-05'
  constraint_meaning = 'kg_solvent_water free_concentration       free_concentration    free_concentration      free_concentration     free_concentration       free_concentration      free_concentration       bulk_composition free_concentration       free_concentration'
  constraint_unit = '   kg               molal               molal            molal              molal             molal               molal              molal               moles              molal               molal'
  initial_temperature = 220
  temperature = temperature
  kinetic_species_name = '         Albite             Anorthite          K-feldspar         Quartz             Phlogopite         Paragonite         Calcite            Anhydrite          Chalcedony         Illite             Kaolinite          Clinochl-7A        Zoisite            Laumontite'
  kinetic_species_initial_value = '4.324073236492E+02 4.631370307325E+01 2.685015418378E+02 7.720095013956E+02 1.235192062541E+01 7.545461404965E-01 4.234651808835E-04 4.000485907930E-04 4.407616361072E+00 1.342524904876E+01 1.004823151125E+00 4.728132387707E-01 7.326007326007E-01 4.818116116598E-01'
  kinetic_species_unit = '         moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles              moles'
  evaluate_kinetic_rates_always = true # otherwise will easily "run out" of dissolving species
  source_species_names = 'H2O H+ Na+ K+ Ca++ Mg++ SiO2(aq) Al+++ Cl- SO4-- HCO3-'
  source_species_rates = 'rate_H2O_per_1l rate_H_per_1l rate_Na_per_1l rate_K_per_1l rate_Ca_per_1l rate_Mg_per_1l rate_SiO2_per_1l rate_Al_per_1l rate_Cl_per_1l rate_SO4_per_1l rate_HCO3_per_1l'
  ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
  execute_console_output_on = ''
  add_aux_molal = false # save some memory and reduce variables in output exodus
  add_aux_mg_per_kg = false # save some memory and reduce variables in output exodus
  add_aux_free_mg = false # save some memory and reduce variables in output exodus
  add_aux_activity = false # save some memory and reduce variables in output exodus
  add_aux_bulk_moles = false # save some memory and reduce variables in output exodus
  adaptive_timestepping = true
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 15
    ny = 10
    xmin = -100
    xmax = 200
    ymin = -100
    ymax = 100
  []
  [injection_node]
    input = gen
    type = ExtraNodesetGenerator
    new_boundary = injection_node
    coord = '0 0 0'
  []
[]

[Executioner]
  type = Transient
  [TimeStepper]
    type = FunctionDT
    function = 'max(1E6, 0.3 * t)'
  []
  end_time = 4E12
[]

[AuxVariables]
  [temperature]
    initial_condition = 220.0
  []
  [porosity]
    initial_condition = 0.01
  []
  [nodal_void_volume]
  []
  [free_cm3_Kfeldspar] # necessary because of the minus sign in K-feldspar which does not parse correctly in the porosity AuxKernel
  []
  [free_cm3_Clinochl7A] # necessary because of the minus sign in Clinochl-7A which does not parse correctly in the porosity AuxKernel
  []
  [pf_rate_H] # change in H mass (kg/s) at each node provided by the porous-flow simulation
  []
  [pf_rate_Na]
  []
  [pf_rate_K]
  []
  [pf_rate_Ca]
  []
  [pf_rate_Mg]
  []
  [pf_rate_SiO2]
  []
  [pf_rate_Al]
  []
  [pf_rate_Cl]
  []
  [pf_rate_SO4]
  []
  [pf_rate_HCO3]
  []
  [pf_rate_H2O] # change in H2O mass (kg/s) at each node provided by the porous-flow simulation
  []
  [rate_H_per_1l]
  []
  [rate_Na_per_1l]
  []
  [rate_K_per_1l]
  []
  [rate_Ca_per_1l]
  []
  [rate_Mg_per_1l]
  []
  [rate_SiO2_per_1l]
  []
  [rate_Al_per_1l]
  []
  [rate_Cl_per_1l]
  []
  [rate_SO4_per_1l]
  []
  [rate_HCO3_per_1l]
  []
  [rate_H2O_per_1l]
  []
  [transported_H]
  []
  [transported_Na]
  []
  [transported_K]
  []
  [transported_Ca]
  []
  [transported_Mg]
  []
  [transported_SiO2]
  []
  [transported_Al]
  []
  [transported_Cl]
  []
  [transported_SO4]
  []
  [transported_HCO3]
  []
  [transported_H2O]
  []
  [transported_mass]
  []
  [massfrac_H]
  []
  [massfrac_Na]
  []
  [massfrac_K]
  []
  [massfrac_Ca]
  []
  [massfrac_Mg]
  []
  [massfrac_SiO2]
  []
  [massfrac_Al]
  []
  [massfrac_Cl]
  []
  [massfrac_SO4]
  []
  [massfrac_HCO3]
  []
  [massfrac_H2O]
  []
[]

[AuxKernels]
  [free_cm3_Kfeldspar]
    type = GeochemistryQuantityAux
    variable = free_cm3_Kfeldspar
    species = 'K-feldspar'
    quantity = free_cm3
    execute_on = 'timestep_begin timestep_end'
  []
  [free_cm3_Clinochl7A]
    type = GeochemistryQuantityAux
    variable = free_cm3_Clinochl7A
    species = 'Clinochl-7A'
    quantity = free_cm3
    execute_on = 'timestep_begin timestep_end'
  []
  [porosity_auxk]
    type = ParsedAux
    coupled_variables = 'free_cm3_Albite free_cm3_Anhydrite free_cm3_Anorthite free_cm3_Calcite free_cm3_Chalcedony free_cm3_Clinochl7A free_cm3_Illite free_cm3_Kfeldspar free_cm3_Kaolinite free_cm3_Quartz free_cm3_Paragonite free_cm3_Phlogopite free_cm3_Zoisite free_cm3_Laumontite'
    expression = '1000.0 / (1000.0 + free_cm3_Albite + free_cm3_Anhydrite + free_cm3_Anorthite + free_cm3_Calcite + free_cm3_Chalcedony + free_cm3_Clinochl7A + free_cm3_Illite + free_cm3_Kfeldspar + free_cm3_Kaolinite + free_cm3_Quartz + free_cm3_Paragonite + free_cm3_Phlogopite + free_cm3_Zoisite + free_cm3_Laumontite)'
    variable = porosity
    execute_on = 'timestep_end'
  []
  [nodal_void_volume_auxk]
    type = NodalVoidVolumeAux
    variable = nodal_void_volume
    nodal_void_volume_uo = nodal_void_volume_uo
    execute_on = 'initial timestep_end' # "initial" to ensure it is properly evaluated for the first timestep
  []
  [rate_H_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_H nodal_void_volume'
    variable = rate_H_per_1l
    expression = 'pf_rate_H / 1.0079 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_Na_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_Na nodal_void_volume'
    variable = rate_Na_per_1l
    expression = 'pf_rate_Na / 22.9898 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_K_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_K nodal_void_volume'
    variable = rate_K_per_1l
    expression = 'pf_rate_K / 39.0983 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_Ca_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_Ca nodal_void_volume'
    variable = rate_Ca_per_1l
    expression = 'pf_rate_Ca / 40.08 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_Mg_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_Mg nodal_void_volume'
    variable = rate_Mg_per_1l
    expression = 'pf_rate_Mg / 24.305 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_SiO2_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_SiO2 nodal_void_volume'
    variable = rate_SiO2_per_1l
    expression = 'pf_rate_SiO2 / 60.0843 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_Al_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_Al nodal_void_volume'
    variable = rate_Al_per_1l
    expression = 'pf_rate_Al / 26.9815 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_Cl_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_Cl nodal_void_volume'
    variable = rate_Cl_per_1l
    expression = 'pf_rate_Cl / 35.453 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_SO4_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_SO4 nodal_void_volume'
    variable = rate_SO4_per_1l
    expression = 'pf_rate_SO4 / 96.0576 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_HCO3_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_HCO3 nodal_void_volume'
    variable = rate_HCO3_per_1l
    expression = 'pf_rate_HCO3 / 61.0171 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [rate_H2O_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_H2O nodal_void_volume'
    variable = rate_H2O_per_1l
    expression = 'pf_rate_H2O / 18.01801802 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []
  [transported_H_auxk]
    type = GeochemistryQuantityAux
    variable = transported_H
    species = 'H+'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_Na_auxk]
    type = GeochemistryQuantityAux
    variable = transported_Na
    species = 'Na+'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_K_auxk]
    type = GeochemistryQuantityAux
    variable = transported_K
    species = 'K+'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_Ca_auxk]
    type = GeochemistryQuantityAux
    variable = transported_Ca
    species = 'Ca++'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_Mg_auxk]
    type = GeochemistryQuantityAux
    variable = transported_Mg
    species = 'Mg++'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_SiO2_auxk]
    type = GeochemistryQuantityAux
    variable = transported_SiO2
    species = 'SiO2(aq)'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_Al_auxk]
    type = GeochemistryQuantityAux
    variable = transported_Al
    species = 'Al+++'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_Cl_auxk]
    type = GeochemistryQuantityAux
    variable = transported_Cl
    species = 'Cl-'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_SO4_auxk]
    type = GeochemistryQuantityAux
    variable = transported_SO4
    species = 'SO4--'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_HCO3_auxk]
    type = GeochemistryQuantityAux
    variable = transported_HCO3
    species = 'HCO3-'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_H2O_auxk]
    type = GeochemistryQuantityAux
    variable = transported_H2O
    species = 'H2O'
    quantity = transported_moles_in_original_basis
    execute_on = 'timestep_begin'
  []
  [transported_mass_auxk]
    type = ParsedAux
    coupled_variables = ' transported_H transported_Na transported_K transported_Ca transported_Mg transported_SiO2 transported_Al transported_Cl transported_SO4 transported_HCO3 transported_H2O'
    variable = transported_mass
    expression = 'transported_H * 1.0079 + transported_Cl * 35.453 + transported_SO4 * 96.0576 + transported_HCO3 * 61.0171 + transported_SiO2 * 60.0843 + transported_Al * 26.9815 + transported_Ca * 40.08 + transported_Mg * 24.305 + transported_K * 39.0983 + transported_Na * 22.9898 + transported_H2O * 18.01801802'
    execute_on = 'timestep_end'
  []
  [massfrac_H_auxk]
    type = ParsedAux
    coupled_variables = 'transported_H transported_mass'
    variable = massfrac_H
    expression = 'transported_H * 1.0079 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_Na_auxk]
    type = ParsedAux
    coupled_variables = 'transported_Na transported_mass'
    variable = massfrac_Na
    expression = 'transported_Na * 22.9898 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_K_auxk]
    type = ParsedAux
    coupled_variables = 'transported_K transported_mass'
    variable = massfrac_K
    expression = 'transported_K * 39.0983 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_Ca_auxk]
    type = ParsedAux
    coupled_variables = 'transported_Ca transported_mass'
    variable = massfrac_Ca
    expression = 'transported_Ca * 40.08 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_Mg_auxk]
    type = ParsedAux
    coupled_variables = 'transported_Mg transported_mass'
    variable = massfrac_Mg
    expression = 'transported_Mg * 24.305 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_SiO2_auxk]
    type = ParsedAux
    coupled_variables = 'transported_SiO2 transported_mass'
    variable = massfrac_SiO2
    expression = 'transported_SiO2 * 60.0843 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_Al_auxk]
    type = ParsedAux
    coupled_variables = 'transported_Al transported_mass'
    variable = massfrac_Al
    expression = 'transported_Al * 26.9815 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_Cl_auxk]
    type = ParsedAux
    coupled_variables = 'transported_Cl transported_mass'
    variable = massfrac_Cl
    expression = 'transported_Cl * 35.453 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_SO4_auxk]
    type = ParsedAux
    coupled_variables = 'transported_SO4 transported_mass'
    variable = massfrac_SO4
    expression = 'transported_SO4 * 96.0576 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_HCO3_auxk]
    type = ParsedAux
    coupled_variables = 'transported_HCO3 transported_mass'
    variable = massfrac_HCO3
    expression = 'transported_HCO3 * 61.0171 / transported_mass'
    execute_on = 'timestep_end'
  []
  [massfrac_H2O_auxk]
    type = ParsedAux
    coupled_variables = 'transported_H2O transported_mass'
    variable = massfrac_H2O
    expression = 'transported_H2O * 18.01801802 / transported_mass'
    execute_on = 'timestep_end'
  []
[]

[GlobalParams]
  point = '0 0 0'
  reactor = reactor
[]

[Postprocessors]
  [temperature]
    type = PointValue
    variable = 'solution_temperature'
  []
  [porosity]
    type = PointValue
    variable = porosity
  []
  [solution_temperature]
    type = PointValue
    variable = solution_temperature
  []
  [massfrac_H]
    type = PointValue
    variable = massfrac_H
  []
  [massfrac_Na]
    type = PointValue
    variable = massfrac_Na
  []
  [massfrac_K]
    type = PointValue
    variable = massfrac_K
  []
  [massfrac_Ca]
    type = PointValue
    variable = massfrac_Ca
  []
  [massfrac_Mg]
    type = PointValue
    variable = massfrac_Mg
  []
  [massfrac_SiO2]
    type = PointValue
    variable = massfrac_SiO2
  []
  [massfrac_Al]
    type = PointValue
    variable = massfrac_Al
  []
  [massfrac_Cl]
    type = PointValue
    variable = massfrac_Cl
  []
  [massfrac_SO4]
    type = PointValue
    variable = massfrac_SO4
  []
  [massfrac_HCO3]
    type = PointValue
    variable = massfrac_HCO3
  []
  [massfrac_H2O]
    type = PointValue
    variable = massfrac_H2O
  []
  [cm3_Albite]
    type = PointValue
    variable = 'free_cm3_Albite'
  []
  [cm3_Anhydrite]
    type = PointValue
    variable = 'free_cm3_Anhydrite'
  []
  [cm3_Anorthite]
    type = PointValue
    variable = 'free_cm3_Anorthite'
  []
  [cm3_Calcite]
    type = PointValue
    variable = 'free_cm3_Calcite'
  []
  [cm3_Chalcedony]
    type = PointValue
    variable = 'free_cm3_Chalcedony'
  []
  [cm3_Clinochl-7A]
    type = PointValue
    variable = 'free_cm3_Clinochl-7A'
  []
  [cm3_Illite]
    type = PointValue
    variable = 'free_cm3_Illite'
  []
  [cm3_K-feldspar]
    type = PointValue
    variable = 'free_cm3_K-feldspar'
  []
  [cm3_Kaolinite]
    type = PointValue
    variable = 'free_cm3_Kaolinite'
  []
  [cm3_Quartz]
    type = PointValue
    variable = 'free_cm3_Quartz'
  []
  [cm3_Paragonite]
    type = PointValue
    variable = 'free_cm3_Paragonite'
  []
  [cm3_Phlogopite]
    type = PointValue
    variable = 'free_cm3_Phlogopite'
  []
  [cm3_Zoisite]
    type = PointValue
    variable = 'free_cm3_Zoisite'
  []
  [cm3_Laumontite]
    type = PointValue
    variable = 'free_cm3_Laumontite'
  []
  [cm3_mineral]
    type = LinearCombinationPostprocessor
    pp_names = 'cm3_Albite cm3_Anhydrite cm3_Anorthite cm3_Calcite cm3_Chalcedony cm3_Clinochl-7A cm3_Illite cm3_K-feldspar cm3_Kaolinite cm3_Quartz cm3_Paragonite cm3_Phlogopite cm3_Zoisite cm3_Laumontite'
    pp_coefs = '1 1 1 1 1 1 1 1 1 1 1 1 1 1'
  []
  [pH]
    type = PointValue
    variable = 'pH'
  []
[]

[Outputs]
  [exo]
    type = Exodus
    execute_on = final
  []
  csv = true
[]

Exploring this input file reveals:

• about half the file is virtually identical to the natural_reservoir.i simulation studied above that simulates the in-situ kinetics without the influence of an injectate

• the other half of the input file involves the interface with the porous_flow.i simulation.

The PorousFlow simulation (listed in the previous section) controls the flow of aqueous species, and interfaces with the aquifer-geochemistry simulation.

1. The PorousFlow simulation provides the AuxVariables called pf_rate_H, pf_rate_Na, pf_rate_K, etc. These are the changes of mass of each species at each node.

2. The PorousFlow simulation also provides an AuxVariable temperature that specifies the temperature at each node.

3. Since the pf_rate_* variables have units kg.s but the Geochemistry module expects rates-of-changes of moles, a conversion must take place. Secondly, the aquifer-geochemistry simulation considers just 1 litre of aqueous solution at every node, so the nodal_void_volume (amount of fluid at each node) is used in the conversion:

  [rate_H_per_1l_auxk]
    type = ParsedAux
    coupled_variables = 'pf_rate_H nodal_void_volume'
    variable = rate_H_per_1l
    expression = 'pf_rate_H / 1.0079 / nodal_void_volume'
    execute_on = 'timestep_begin'
  []

4. The geochemistry code solves the geochemical system at each node, which depends on the rates of source species (rate_*_per_1l) and the kinetics.

5. The mass-fraction of each species at each node is computed from the transported mole numbers, and sent back to the porous_flow.i simulation in order for the next step of transport:

  [transported_mass_auxk]
    type = ParsedAux
    coupled_variables = ' transported_H transported_Na transported_K transported_Ca transported_Mg transported_SiO2 transported_Al transported_Cl transported_SO4 transported_HCO3 transported_H2O'
    variable = transported_mass
    expression = 'transported_H * 1.0079 + transported_Cl * 35.453 + transported_SO4 * 96.0576 + transported_HCO3 * 61.0171 + transported_SiO2 * 60.0843 + transported_Al * 26.9815 + transported_Ca * 40.08 + transported_Mg * 24.305 + transported_K * 39.0983 + transported_Na * 22.9898 + transported_H2O * 18.01801802'
    execute_on = 'timestep_end'
  []
  [massfrac_H_auxk]
    type = ParsedAux
    coupled_variables = 'transported_H transported_mass'
    variable = massfrac_H
    expression = 'transported_H * 1.0079 / transported_mass'
    execute_on = 'timestep_end'
  []

Some results are shown in Figure 11 to Figure 16. The following may be surmised from these figures.

• The regions 50m from the injection well, and between the two wells experience temperature changes.

• The region with approximately 300m from the wells is impacted by the injectate fluid chemistry. The impacts are more obvious for some species than others.

• The porosity decreases in response to the injectate, but it appears that porosity changes in response to cool injectate more than heated injectate.

Figure 11: Temperature contour after 1 year of injection in the reactive-transport simulation. The green dots show the injection and production wells and the black ring is the 100C contour.

Figure 12: Porosity contour after 1 year of injection in the reactive-transport simulation. The green dots show the injection and production wells and the black ring is the 100C contour.

Figure 13: Contour of free Albite volume after 1 year of injection in the reactive-transport simulation. The green dots show the injection and production wells and the black ring is the 100C contour.

Figure 14: Contour of pH after 1 year of injection in the reactive-transport simulation. The green dots show the injection and production wells and the black ring is the 100C contour.

Figure 15: Contour of SiO2(aq) mass fraction after 1 year of injection in the reactive-transport simulation. The green dots show the injection and production wells and the black ring is the 100C contour.

Figure 16: Contour of SO4– mass fraction after 1 year of injection in the reactive-transport simulation. The green dots show the injection and production wells and the black ring is the 100C contour.

References