All phases of the database in the simulation.

[PrestoConFuoco]

I have a lot of phases in my database, and sometimes i need to consider them all during simulation. Is there any such mechanism in PHREEQC? I don't want to write all phasenames by hand.
Thanks.

[dlparkhurst]

No, there is not an easy way to list all the phases.

Here is a hack that will work if you have an initial simulation that contains all of the elements that will be encountered. In this case, subsequent simulations (MIX) list SIs in the USER_PRINT data block.  You can change the code to write a USER_PUNCH data block instead of USER_PRINT.

Code: [Select]

TITLE Example 1.--Add uranium and speciate seawater.
SOLUTION 1  SEAWATER FROM NORDSTROM AND OTHERS (1979)
        units   ppm
        pH      8.22
        pe      8.451
        density 1.023
        temp    25.0
        redox   O(0)/O(-2)
        Ca              412.3
        Mg              1291.8
        Na              10768.0
        K               399.1
        Fe              0.002
        Mn              0.0002  pe
        Si              4.28
        Cl              19353.0
        Alkalinity      141.682 as HCO3
        S(6)            2712.0
        O(0)            1.0     O2(g) -0.7
SELECTED_OUTPUT 2
-file phase_list
-reset false
USER_PUNCH 2
10 PUNCH "USER_PRINT", EOL$
20 t = SYS("phases",  count , name$ , type$ , moles )
25 PUNCH str_f$(n*10,4,0), "PRINT '", pad("Phase",20), "',", "'   SI   '", EOL$
30 for n = 1 to count
40 si$ = "SI('" + name$(n) + "')"
50 PUNCH str_f$(n*10,4,0), "PRINT '", pad(name$(n),20), "'", si$, EOL$
60 next n
70 PUNCH "END", EOL$
END
SELECTED_OUTPUT 2
-active false
END
INCLUDE$ phase_list
MIX 1
1 1
END


[John Mahoney]

Your question really gets to the heart of geochemical modeling, which is how does one define activities of the components in a model?  The primary function of a geochemical modeling program is to calculate the activity of a basis set of components for a user defined solution.  It can do that in a number of ways, but I will not get too deeply into that subject here.  Also remember that activity is not the same as a total concentration, but if all the activities are defined then the concentrations can be calculated based upon the speciation reactions.

I have two main concerns, and then some thoughts on possible model setups if you decide to continue:
1)   Database problems in general.  This is discussed for all users and may not apply if you truly have an reliable database.
2)   Significance and ramifications of fixing phases.
3)   Options to setup such a model. 

1 - I would be careful if you intend to “use” all of the phases in your database.  I assume you mean to include them as possible phases in the EQUILIBRIUM_PHASES Keyword block.  Using a format like:

EQUILIBRIUM_PHASES
Phase_name 0 0
Phase2 0 0
And so on.

Unless you have personally vetted all the data in your database, I suspect that there will be numerous issues among the phases and their reactions.   If you are using values from the different databases provided with PHREEQC you will probably have problems.  Several of the databases provided with PHREEQC contain various errors and they have not been corrected.  This would apply mainly to the WATEQ4F.dat and the MINTEQ databases provided with the program.   The LLNL.dat database may be a little better as a lot of that comes out or SPRCRT92, but I would still be careful.

Because of the way many phases have been added to these databases, there are going to be inconsistencies between the solution (speciation-based) model and the assumptions (another model really) use to define the solubility product constants for the phases in the database.   Claims have been made that these databases are internally consistent, but I do not believe that.  In addition to the inconsistencies between phases and speciation reactions, I know that several reactions and phases in the WATEQ4F.dat database are incorrect.   

I would suggest you read the chapter by Rolf Grauer called “Solubility Limitations: An Old Timer’s View”   in Modelling in Aquatic Chemistry for a realistic assessment of databases ; here is the link. 

https://www.oecd-nea.org/dbtdb/pubs/book-pdf/131-152.pdf

The original German text is followed by an English Translation.   It is an irreverent discussion about database problems.   Although this was written over 20 years ago, sadly things are only slightly better.    But at least the newer ANDRA ThermoChemie TDB database appears to indicate some improvement, but I have not gone over this new database in any detailed level to confirm this statement. 

2 - Whenever a phase is selected, whether explicitly by the modeler by including it in an EQUILIBRIUM_PHASE block and including a molar concentration, or by including it as just a possible phase (let’s call that an implicitly fixed phase, no initial concentration, but it may form if it is oversaturated), which the program may select, you are making an important decision about the system being modeled.  Every time a phase is added to a model, one is “playing god” as is implies the modeler knows more about certain conditions, you are fixing an activity in the model.  The most obvious example is fixing a set of pH values to prepare a speciation diagram over a range of pH values.  Computationally, the same thing happens when you fix a solid or gas phase in a model.  Every phase that is fixed lessens the degrees of freedom.  Given the way PHREEQC sets up a problem for the model you appear to be suggesting, the program will select the most stable set of phases based upon the number of components in the solution.  Thus, PHREEQC will pick Hematite over Ferrihydrite as an EQUILIBRIUM_PHASE, and it will pick quartz over all of the other silica polymorphs.  If this is what is really going on, then that is fine, but usually the less stable phases form first under low temperature conditions.  For simplicity let’s start by fixing the pH.  Then if you have precipitated an iron phase, let’s say ferrihydrite, the activity of Fe+3 is defined.  Now if you have enough arsenic to precipitate scorodite the activity of arsenic is now fixed. 

This type of process follows for other phases also.  Fixing a gas phase such as CO2(g) will control the activity of H2CO3 under all pH conditions and therefore the activity of CO3-2 is controlled by the pH and the H2CO3 activity.   So now if calcite forms then the Ca+2 activity is defined through the solubility product for calcite and the previously calculated activity of CO3-2. 
If these stable phases are added into the model, then component activities will decrease and more and more phases will become unsaturated.  As this proceeds PHREEQC will have picked out the most stable phases, one for each component and then it will stop.  This is another example of how the Phase Rule applies to these models.

Because PHREEQC parses the database and lists only the phases based upon the components you have in the solution you will only have the same number of phases as components in your solution.  This is the same reason that the list of saturation indices is limited to phases that can be defined by the components used to define the solution.   This gets into a concept called saturation sufficiency, which says you can not calculate a saturation index if you do not have the appropriate data.   That is the quick explanation of saturation sufficiency or at least my take on the subject, there is more to it than my quick explanation.  Refer to Plummer, Parkhurst, Thorstenson,  1983,  Development of reaction models for ground-water systems, Geochimica et Cosmochimica Acta, vol 47, pp. 665-686 for a better discussion on this topic. 

3- If you still want to list all the phases there are two options.  The first one is to run a solution model and copy the saturation indices block and paste that into EXCEL and clean it up, and then use that in the EQUILIBRIUM _PHASES block for another run.  Or you can do it using a text editor, such as TEXTPAD, that has a block edit function.    A third option is you can use the following USER_PRINT prepared by David Parkhurst about five years ago.  As it is set up now it will list all the phases that have saturation indices greater than -2.  This will reduce the number of phases you have to consider.  The SI value can be changed in line 30.   You can then block this out and clean it up and put it back in the input file. 

Code: [Select]

USER_PRINT
-START
5 PRINT "TEMPERATURE      :", TC,"°C"
10 x = SYS("phases", count, name$, type$, si_values)
20 FOR i = 1 to count
30  IF si_values(i) <-2 then goto 50
40  PRINT PAD(name$(i),20), si_values(i), Phase_Formula(name$(i))
50 NEXT i


There is a lot more I could say about this topic,  but that is enough for now.