Discrepancy between SIs and SIs from 'get_selected_output'

[jeremyfirst22]

Hello,

There seems to be a discrepancy between the saturation indices returned by Solution.si(), and those returned by PhreeqPython.ip.get_selected_output.

If I run:

from phreeqpython import PhreeqPython
pp = PhreeqPython("/hpcdata/fijmjl/iphreeqc-3.7.3-15968/database/phreeqc.dat", debug=True)
soln = pp.add_solution(        
    {    
        "units": "ppm",
        "Mg": 516, 
        "Mn(2)": 2.6, 
    }, 
)
minerals = ["Hausmannite", "Manganite", "Pyrochroite", "Pyrolusite"]
{k: soln.si(k) for k in minerals}

I get the following SIs:

SOLUTION 0
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 0
END 

{'Hausmannite': -18.90189294986606,
 'Manganite': -8.963964316622022,
 'Pyrochroite': -5.823964316622019,
 'Pyrolusite': -18.003964316622024}

But using the selected output array:

from phreeqpython import PhreeqPython
pp = PhreeqPython("/hpcdata/fijmjl/iphreeqc-3.7.3-15968/database/phreeqc.dat", debug=True)
phases = soln.phases.keys()
pp.ip.run_string(
"""
SOLUTION 1
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 2

SELECTED_OUTPUT 1
-file selected_output.sel
-reset false
-saturation_indices Hausmannite Manganite Pyrochroite Pyrolusite
END
"""
)
{k: v for k, v in zip(*pp.ip.get_selected_output_array())}

I get different SIs for the minerals:

SOLUTION 1
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 2

SELECTED_OUTPUT 1
-file selected_output.sel
-reset false
-saturation_indices Hausmannite Manganite Pyrochroite Pyrolusite
END

{'si_Hausmannite': -10.902521584503347,
 'si_Manganite': -4.964278633940662,
 'si_Pyrochroite': -5.824278633940661,
 'si_Pyrolusite': -10.004435792599985}

Worth noting that the latter results match what I get from running the phreeqc desktop application. Any ideas on what might be going on?

Thank you in advance for the help!

[AbelHeinsbroek]

Hi Jeremy,

I think I've solved this issue in the fix_si branch of phreeqpython, can you test and confirm?

Regards,

Abel

[jeremyfirst22]

Hi Abel,

Thank you for the quick response! I am still getting slightly different results in the 'fix_si' branch, but I am not sure if they are significantly different or not. For the first code snippet above, I get:

SOLUTION 0
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 0
END 

{'Hausmannite': -10.901892948191104,
 'Manganite': -4.963964316063699,
 'Pyrochroite': -5.8239643160637,
 'Pyrolusite': -10.003964316063701}

While for the second snippet, I get:

SOLUTION 1
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 2

SELECTED_OUTPUT 1
-file selected_output.sel
-reset false
-saturation_indices Hausmannite Manganite Pyrochroite Pyrolusite
END

{'si_Hausmannite': -10.902521582829458,
 'si_Manganite': -4.964278633382879,
 'si_Pyrochroite': -5.824278633382876,
 'si_Pyrolusite': -10.004435792042472}

These are significantly closer than previously, but I figure the Python code is directly fetching the SIs from the IPhreeqc instance, so I would expect them to be exact. Is that not the case?

-- Jeremy

[AbelHeinsbroek]

Hi Jeremy,

The way (I)PhreeqC works is; run a single simulation, print/output all relevant data, close simulation and store only the solution input data to be used in a subsequent simulation.

To make PhreeqPython's object-oriented approach work I've had to modify the inner workings of IPHREEQC to store certain properties, such as pH, speciation, specific conductance, etc after each simulation run, so that they can be requested 'on the fly' instead of having to specify the requested output beforehand as you do in your example.

The code-base of IPhreeqC is a bit obtuse in some places, with loads of global variables that get accessed and changes from lots of different files. I think I made a small error in storing the saturation indices during simulation run, but will need some time to figure out what causes this subtle difference, as they should be exactly the same.

The SI calculation for example is derived from this piece of code: https://github.com/Vitens/VIPhreeqc/blob/master/src/phreeqcpp/print.cpp#L1180