Multiple questions on hydrothermal modeling

[Flopi]

Dear Forum,
I am trying to use PHREEQC to simulate the equilbrium of high temperature minerals (in a hydrothermal setting) and I have a few questions about it.

Below is the code I am currently running:

PHASES
   Serpentine(FeIII)
   Fe2Si2O5(OH)4 + 6H+ = 2Fe+3 + 2H4SiO4 + 1H2O
   -analytic -9.80801362E+02 -1.6542353E-01 4.7339634E+04 3.5649533E+02 -2.1360519E+06
   -Vm 105.03

   Cronstedtite-7A #From Carbfix
   Fe2Fe2SiO5(OH)4 + 10 H+ = H4SiO4 + 2 Fe+2 + 2 Fe+3 + 5 H2O
   log_k      16.2603
   -delta_H   -244.266   kJ/mol
   -analytic   -2.3783e2 -7.1026e-2 1.7752e4 8.7147e1 2.7707e2
   -Vm      110.9 # HDN+78


SOLID_SOLUTIONS 1
   Olivine
      -comp Forsterite 6.12
      -comp Fayalite 0.68

   Orthopyroxene
      -comp Enstatite(alpha) 0
      -comp Ferrosilite(alpha) 0

   Clinopyroxene
      -comp Diopside 0
      -comp Hedenbergite 0

   Serpentine
      -comp Lizardite 0
      -comp Greenalite 0
      -comp Cronstedtite-7A 0
      -comp Kaolinite 0

   Brucite_Fe
      -comp Brucite 0
      -comp Fe(OH)2 0

   Talc_Fe
      -comp Talc 0
      -comp Minnesotaite 0

   Chlorite
      -comp Clinochlore 0
      -comp Chamosite(daphnite) 0

   Tremolite
      -comp Tremolite 0
      -comp Ferrotremolite 0
END

SOLUTION 1
    temp      25
   pressure 500
    pH        7.8
    pe        0
    redox     O(-2)/O(0)
    units     mmol/kgw
    density   1
    Al        3.7e-05
    C(4)      2.34
    Ca        10.2
    Cl        546
    Fe(2)     1.5e-06
    K         9.8
    Mg        23
    Na        464
    O(0)      0.25
    Si        0.11
    -water    1 # kg
SAVE SOLUTION 1
END

USE SOLUTION 1
USE SOLID_SOLUTION 1

REACTION_TEMPERATURE 1
    25 300 in 12 steps

EQUILIBRIUM_PHASES
   magnetite 0 0

PRINT
   selected_output true

SELECTED_OUTPUT
   -file Dunite
   -temp
   -solid_solutions    Forsterite   Fayalite Enstatite(alpha) Ferrosilite(alpha) Diopside   Hedenbergite Lizardite Greenalite Cronstedtite-7A Kaolinite   Brucite Fe(OH)2 Talc Minnesotaite Clinochlore Chamosite(daphnite) Tremolite Ferrotremolite
   -equilibrium_phases magnetite


I have 3 questions about it:

- the main issue I have here is that the code crashes at 200°C. It works fine before. However if I remove the EQUILIBRIUM_PHASES with magnetite, it's all good all the way to 300°C. Any idea why magnetite makes the whole thing crash ? Am I using EQUILIBRIUM_PHASES and SOLID_SOLUTIONS right ?

- If I ever wanted to fix the compositional ratio of a solid_solution. Typically force olivine to be Fo90 (i.e 90%mol of forsterite and 10%mol of fayalite) throughout the whole process. This means that olivine can dissolve but should retain always the same proportion of forsterite and fayalite. What would be the best way to implement it ?

- Finally, a bit of unrelated question. What is the syntax in selected_output to obtain the mole quantity of each solid solution (e.g. writing in the .out file the mole quantity of Serpentine).

Thank you for your help,

[dlparkhurst]

(1) You will need to include the database that you are using to be able to replicate the problems with temperature.

(2) To fix a solid solution composition, you would define a PHASE with the given composition. For ideal solid solution, log K would be sum of the log K's using the appropriate mole fractions. Then include the phase in EQUILIBRIUM_PHASES.

(3) You have SELECTED_OUTPUT; -solid_solution, which will print to the selected output file.

You can use the Basic function S_S("Ferrotremolite") in USER_PUNCH, which also goes to the selected output file.

For the output file, by default there is a block for solid solutions that is printed.

[Flopi]

Thanks for the answers.

I am using the Thermoddem database
https://thermoddem.brgm.fr/databases/phreeqc

I also tried with the CarbFix database https://www.carbfix.com/database and it crashed the same way but at different temperatures.

[dlparkhurst]

It is apparently a difficult numerical problem for PHREEQC. I was able to get a solution using some of the KNOBS options:

Code: [Select]

...
END
KNOBS
-step 5
-pe   2
-iter 500
USE SOLUTION 1
USE SOLID_SOLUTION 1
...

This definition causes the numerical method to take smaller steps in the direction of the numerical solution, and, if necessary, more steps to arrive at the solution. I won't promise that this set of parameters will work for all of your simulations, but perhaps it is a start.

[Flopi]

That does sthe work indeed. Thanks a lot !

[Flopi]

As a follow up with my question, I was wondering about the 300°C limit on the code and how "hard" that limit is. I know it's from the critical point of water but I'm wondering if there is a way to circumvent it and go beyond. And potentially, how hard would it be to modify the core code so that it could routinely calculate speciations and water-rock interactions at high T and P (typically up to 600-800°C and 1 to 2 kbar), similarly to what EQ3/6 can do.

[dlparkhurst]

The fit we use for the dielectric constant of water goes negative at about 300 or 350 C. So, better parameterizations of some of the basic water properties are needed to go to higher temperatures. In addition, the temperature dependence of log Ks and Pitzer parameters is spotty, so a much more careful evaluation of the aqueous model is needed.