How can I export the different pH-Fe species combinations?

[Juan Sebastian]

Hello Mr. Parkhurst,

I understand that questions related to PhreePlot must be directed to David Kinninburg, but I believe you might be able to assist me with this, so I will try it. I have a script that exports a distribution of species diagram for Fe, but I also want to export the different pH-Fe species combinations that I got across the entire plot (in a.sel file). In the USER_PUNCH section, I used an example you provided to export all species and their log molalities in two columns, which has been very useful. However, I'm having trouble exporting the specific information I described.

Hopefully, you can help me, as you always do.

This is my code:

SPECIATION
  jobTitle                         "Fe Monterrey"
  calculationType                     species
  calculationMethod                    1
  mainSpecies                         Fe
  xmin                                 -4.0
  xmax                                 -0.1
  resolution                           100
 
PLOT
  plotTitle                            "Fe speciation"
  customXcolumn                         2
  pxmax                                      4
  minimumYValueForPlotting        5.0
  legendTitle                          "Fe species"
 
 CHEMISTRY

include 'speciesvsph.inc'

SELECTED_OUTPUT
-reset  false
USER_PUNCH
-headings Species_names     Log_Molality
10 t = SYS("aq", count, name$, type$, moles)
20 FOR i = 1 to count
30 PUNCH PAD$(name$(i),12), STR_F$(LM(name$(i)),6, 2), EOL_NOTAB$
40 NEXT i
50 END

PHASES
Fix_H+
H+ = H+
log_k 0.0

SOLUTION 26
pH 4
pe 10.350000
temp 42.800000
-units mg/L
Na 402.943233
K 1038.472893
Mg 28848.877736
Ca 550.174358
Alkalinity 0.000000
S(6) 160666.600000
Cl 236.666667
Al 7648.121721
As 41.176146
B 19.064162
Ba 0.016402
Cd 0.067811
Cu 0.573103
Fe 2011.838793
Li 68.412078
Mn 81.797433
Ni 1.162858
P 32.229078
Pb 1.916603
Rb 0.000000
S 0.000000
Se 0.000000
Si 63.320155
Sr 40.798466
Zn 13.500483
SAVE solution 26
END

USE solution 26
EQUILIBRIUM_PHASES   
   Fix_H+ <x_axis>  HCl 150 
  -force_equality true
   O2(g)      -28.11  0.1
   -force_equality true
   CO2(g)     -3.5    0.1
   -force_equality true
 
END

[Juan Sebastian]

Thank you very much! David Kinninburg mentioned that you reached out to him to help us with this question. He has already helped me resolve some doubts related to this. Thanks again! I?m going to paste David Kinninburg?s response here in case it?s useful to someone else:


Species plots have a special format explained in the speciesvsph.inc
file that you are including. This includes a USER_PUNCH block.

You are redefining this (which you should not) so delete

#SELECTED_OUTPUT
#-reset  false
#USER_PUNCH
#-headings Species_names     Log_Molality
#10 t = SYS("aq", count, name$, type$, moles)
#20 FOR i = 1 to count
#30 PUNCH PAD$(name$(i),12), STR_F$(LM(name$(i)),6, 2), EOL_NOTAB$
#40 NEXT i
#50 END

and it should then give you a species plot.

Then look at the .out file which gives pH, then pairs of species name,
species %age, ... in decreasing abundance.
The .pts file simplifies this to a simple species-concn table suitable
for importing into spreadsheets.

To get the molal concentrations that I think you want, you should modify
the .inc file by making it PUNCH concentrations not %ages. See my
species_concnvsph.inc  file.

Then add this to the directory with your ppi file (or the system
directory):

include species_concnvsph.inc

NB to get 'all T' to work, change PRINT;-reset T in the include file.
Also look at the .all for combined phreeqc.out output.

[David Kinniburgh]

I'll add two points to my initial response to Juan pasted above:

(i) running the code that Juan posted stopped abruptly during the plotting phase with a message:

*** Error: Allocation error at 398 53224684 xg,yg,curve2column.

This is a bit cryptic but arises because PhreePlot is trying to create an array of size 30800 x 25380432! This is just too large and plainly wrong. It arises from the overwriting of the correct USER_PUNCH block in speciesvsph.inc with the one defined by Juan - which has the wrong format - as described above. So either delete this one and modify the speciesvsph.inc block to your needs, or define a second SELECTED_OUTPUT/USER_PUNCH pair, e.g. SELECTED_OUTPUT 2/USER_PUNCH 2 to write the new data to, probably adding -file to direct to a file (PhreePlot always uses 1 as its principal output). Also add

Code: [Select]

selectedOutputFile TRUE
to the PhreePlot section to force the file to actually be written.

(ii) for completeness, to get what Juan wanted, the revised species_concnvsph.inc file contains this:

Code: [Select]

USER_PUNCH
  -start
10 main$ = "<mainspecies>"
20 PUNCH "pH", -LA("H+")
30 totel = SYS(main$, n, n$, t$, moles)
35 IF (totel <=0) THEN 120
36 h2o = TOT("water")
40 FOR i = 1 TO n
50   IF(MID$(n$(i), 1, 2) = "H(" OR MID$(n$(i), 1, 2) = "O(" OR n$(i) = "O" OR n$(i) = "H") THEN 100
60   IF(t$(i) = "equi") THEN 70 ELSE 90
70   IF (instr(n$(i),"(g)")=0) THEN PUNCH n$(i)+"(s)", moles(i)/h2o ELSE PUNCH n$(i), moles(i)/h2o
80   GOTO 100
90   PUNCH n$(i), moles(i)/h2o
100 NEXT i
120 END


[David Kinniburgh]

There is now another way to create the species concentration-pH table that Juan wanted. This makes use of the relatively new optional sixth parameter to the SYS() command. Setting this to 1 changes the output ordering from abundance-based to alphabetic. This then provides the means to get a fixed order output which is a necessary precursor to producing such tables and plots.

This means that you can now use PhreePlot's 'custom' type rather than its 'species' type. This makes the code more transparent and amenable to customization. There are three steps involved: (i) generate the list of species involved to be used in the USER_PUNCH;-headings block and for plotting with PhreePlot's 'lines' setting; (ii) produce the output table; (iii) generate the plot.

The following PhreePlot code does this. There are two simulations: the first is the preloop simulation that runs a single instance of the speciation at arbitrary pH and PUNCH'es the species lists for both the -heading and 'lines' settings. This output automatically generate tags that can be used in subsequent simulations. The second simulation is the main loop and iterates 101 times from from pH 2 to 12. It uses the generated <heading$> tag to produce the 'out' file. Finally the plot is made from the 'out' file with the lines to be plotted specified by the generated <aqlist$> tag and subject to the given minimumYValueForPlotting setting.

Note: this relatively simple example is to demonstrate the mechanics only. In practice, you would want to consider other species, including mineral species, that may be involved, e.g. here by including additional EQUILIBRIUM_PHASES, manually or automatically.

Code: [Select]

# PhreePlot demo to show how to make a spreadsheet-like file with aq species concn vs pH for Fe.
# Accumulated tabular output is sent to the 'out' file.
# Also plots the results.

# For PhreePlot, see [url]https://www.phreeplot.org[/url].

# N.B. This uses a 'custom' calculationType rather than a 'species' type.

SPECIATION
  calculationType          custom
  calculationMethod        1
  xmin                     2
  xmax                     12
  resolution               101
  customxColumn            pH

PLOT
  ytitle                   "species concn (mol/kgw)"
  lines                    <aqlist$>
  minimumYValueForPlotting 1e-5

CHEMISTRY

######################################################################

# SIMULATION 1 - preloop
# execute once to generate the tags containing the species for the -heading and 'lines' setting

PHASES
Fix_H+
H+ = H+
log_k 0.0

SOLUTION 1
-units mol/kgw
pH 1.5
Na 0.1
Cl 0.13
Fe 1e-3

SELECTED_OUTPUT 1
  -reset                false

USER_PUNCH 1
  -headings aqlist$ heading$  # these automatically make character tags
10 aqlist$ = ''
20 main$ = "Fe"
30 IF(STEP_NO = 0) THEN 130
40 REM Now concatenate the species names
50 REM N.B. the optional sixth parameter (1) means that the species output is in a fixed alphabetic order
60 totel = SYS(main$, n, n$, t$, moles, 1)
70 FOR i = 1 TO n
80   aqlist$ =  aqlist$ + " " + n$(i)
90 NEXT i
100 PUNCH aqlist$
110 heading$ = "pH" + aqlist$
120 PUNCH heading$
130 END

EQUILIBRIUM_PHASES 1 
   Fix_H+ -2  NaOH
  -force_equality true
   O2(g) -0.67
  -force_equality true
END

######################################################################

# SIMULATION 2 - last so by default, main loop
# loop on this to generate the data

USE solution 1

SELECTED_OUTPUT 1
  -reset false

USER_PUNCH 1
    -heading <heading$>
10 main$ = "Fe"
20 IF(STEP_NO = 0) THEN 90
30 PUNCH <x_axis>
40 REM Now output the concentrations in alphabetic order
50 totel = SYS(main$, n, n$, t$, moles, 1)
60 FOR i = 1 TO n
70   PUNCH MOL(n$(i)) # remember for mass balance, species may be polynuclear
80 NEXT i
90 END

EQUILIBRIUM_PHASES 1 
   Fix_H+ -<x_axis>  NaOH
  -force_equality true
   O2(g) -0.67
  -force_equality true
END


[Juan Sebastian]

Thank you very much for this thorough explanation! I was able to run my code using your examples.

This is exactly what I was looking for!