Equilibrate concentration & species wrong using PHREEQC in pyEQL when specified with pure elements concentration

[abhardwaj73]

When specifying [M] without the charge of the metallic cation:

• It doesn't assign it to the respective cation [M+n] as done in PHREEQC
• Causes elements to be double counted

Example:

comp = {"Na": "313.449122807018 mg/L", "Ca": "146.356944444444 mg/L", "Mg": "10.0854166666667 mg/L"}
solution = Solution(comp, balance_charge='auto')
solution.equilibrate()
species_concentrations = solution.list_concentrations(unit='mol/L')
print(species_concentrations)

Output

Na(aq):  313.4491 mg / l
Na[+1]:  313.4491 mg / l
Ca(aq):  146.3569 mg / l
Ca[+2]:  146.3568 mg / l
Mg[+2]:  10.0854 mg / l
Mg(aq):  10.0854 mg / l
OH[-1]:  0.0019 mg / l
CaOH[+1]: 0.0002 mg / l

This output confirms that ionic species like [Na+1] and neutral species Na(aq) are both present. This shows that there could be double counting of elements.

Comparing it with the case of providing charge in the concentration:

comp = {"Na+": "313.449122807018 mg/L", "Ca++": "146.356944444444 mg/L", "Mg++": "10.0854166666667 mg/L"}
solution = Solution(comp, balance_charge='auto')
solution.equilibrate()
species_concentrations = solution.list_concentrations(unit='mol/L')
print(species_concentrations)

Output

OH[-1]:  370.1403 mg / l
Na[+1]:  313.3702 mg / l
Ca[+2]:  146.3199 mg / l
Mg[+2]:  10.0829 mg / l
H[+1]:   0.0001 mg / l
CaOH[+1]:  0.0002 mg / l

[vineetbansal]

@rkingsbury - I've been diving into the code and looked into this a bit, and here's what seems to be happening in a simpler case of:

solution = Solution({"Na": "0.001 mg/L"}, balance_charge='auto', engine='native')

• On initialization, solution has components = 'H2O(aq)', 'OH[-1]', 'H[+1]', 'Na(aq)'
• After equilibrate(), it adds 4 more (after processing through phreeqc) = 'Na[+1]', 'NaOH(aq)', 'O2(aq)', 'H2(aq)'
• phreeqc by itself comes up with species 'H+', 'H2', 'H2O', 'Na+', 'NaOH', 'O2', 'OH-', i.e. it accounts for all species that the solution class was coming up by itself, while correctly assigning valence +1 to Na.

So it seems to me that for the NativeEOS class, its add_solute method can pass its incoming formula (str) argument to the underlying phreeqc engine and get its components that way (still passing it through the FormulaDict class to standardize the formulas). Further invocations of add_solute will keep adding keys to its components.

Let me know if this seems like a reasonable approach or if I'm missing something.

[rkingsbury]

So it seems to me that for the NativeEOS class, its add_solute method can pass its incoming formula (str) argument to the underlying phreeqc engine and get its components that way (still passing it through the FormulaDict class to standardize the formulas). Further invocations of add_solute will keep adding keys to its components.

Thanks for investigating @vineetbansal . Yes, this is consistent with how (I think) the code works now and how it can continue to work. The only problem is that when the phreeqc engine turns Na(aq) into Na[+1] and then hands that composition back to pyEQL, pyEQL will think that Na(aq) is a missing species - one that needs to be added back in. And so the final composition will include both Na(aq) and Na[+1] (with identical amounts) which is not correct:

{'H2O(aq)': 55.34457595832858, 'OH[-1]': 9.6102268597005e-08, 'H[+1]': 5.260466010841792e-08, 'Na(aq)': 4.349760921132654e-08, 'Na[+1]': 4.349760848858708e-08, 'NaOH(aq)': 7.227393962327349e-16, 'O2(aq)': 1.00841670337401e-24, 'H2(aq)': 7.911616444921677e-35}

The reason this behavior exists in pyEQL is that if we start with a Solution where there is a species / element not present in PHREEQC's databases, the phreeqc engine will just silently throw it away. We want to make sure that doesn't happen (basically, guaranteeing mass balance before and after equilibrate. Here are two illustrative cases.

First case - a species (ReO4[-2]) is not in PHREEQC's database, but the element is

This case is actually still a problem, because here, PHREEQC transforms ReO4[-2] into ReO4[-1], but our pyEQL logic retains the original ReO4[-2] as well. We should fix this.

>>> s1 = Solution({"ReO4-2": "0.001 mg/L"}, balance_charge='auto', engine='native', log_level='info')
[09/25/25 08:37:30] WARNING  [2025-09-25 08:37:30,745] [ WARNING] --- Volume of solute ReO4[-2] will be   engines.py:653
                             ignored because partial molar volume data are not available.
                             (engines.py:653)
>>> s1.equilibrate()
[09/25/25 08:37:39] WARNING  [2025-09-25 08:37:39,111] [ WARNING] --- After equilibration, the amounts of engines.py:680
                             species {'ReO4[-2]'} were not modified by PHREEQC. These species are likely
                             absent from its database. (engines.py:680)
>>> s1.components
{'H2O(aq)': 55.34457595832858, 'H[+1]': 9.974194146715052e-08, 'OH[-1]': 9.609888987846696e-08, 'ReO4[-2]': 3.996729076923446e-09, 'ReO4[-1]': 3.996729076923443e-09, 'O2(aq)': 1.0084167095921508e-24, 'H2(aq)': 7.911616483047811e-35}

Second case - entire element is not in PHREEQC database

Here, the element Rh is simply not present at all. If we relied only on PHREEQC to return the components, we would find that the Rh+3 had been discarded. So here, the existing pyEQL logic does what it's supposed to do.

>>> s1 = Solution({"Rh+3": "0.001 mg/L"}, balance_charge='auto', engine='native', log_level='info')
[09/25/25 08:42:46] WARNING  [2025-09-25 08:42:46,538] [ WARNING] --- Volume of solute Rh[+3] will be     engines.py:653
                             ignored because partial molar volume data are not available.
                             (engines.py:653)
                    INFO     [2025-09-25 08:42:46,540] [    INFO] --- Solution is not electroneutral    solution.py:2286
                             (C.B. = 2.9152960726103063e-08 eq/L). (solution.py:2286)
                    INFO     [2025-09-25 08:42:46,541] [    INFO] --- Adjusting OH[-1] to compensate.   solution.py:2292
                             (solution.py:2292)
                    WARNING  [2025-09-25 08:42:46,544] [ WARNING] --- Volume of solute Rh[+3] will be     engines.py:653
                             ignored because partial molar volume data are not available.
                             (engines.py:653)
>>> s1.components
{'H2O(aq)': 55.34457595832858, 'OH[-1]': 1.2915296072610306e-07, 'H[+1]': 1e-07, 'Rh[+3]': 9.717653575367691e-09}
>>> s1.equilibrate()
[09/25/25 08:42:49] WARNING  [2025-09-25 08:42:49,125] [ WARNING] --- After equilibration, the amounts of engines.py:680
                             species {'Rh[+3]'} were not modified by PHREEQC. These species are likely
                             absent from its database. (engines.py:680)
                    INFO     [2025-09-25 08:42:49,129] [    INFO] --- Solution is not electroneutral    solution.py:2286
                             (C.B. = 3.2796550812162786e-08 eq/L). (solution.py:2286)
                    INFO     [2025-09-25 08:42:49,130] [    INFO] --- Adjusting OH[-1] to compensate.   solution.py:2292
                             (solution.py:2292)
                    WARNING  [2025-09-25 08:42:49,138] [ WARNING] --- Volume of solute Rh[+3] will be     engines.py:653
                             ignored because partial molar volume data are not available.
                             (engines.py:653)
                    WARNING  [2025-09-25 08:42:49,140] [ WARNING] --- Volume of solute O2(aq) will be     engines.py:653
                             ignored because partial molar volume data are not available.
                             (engines.py:653)
                    WARNING  [2025-09-25 08:42:49,142] [ WARNING] --- Volume of solute H2(aq) will be     engines.py:653
                             ignored because partial molar volume data are not available.
                             (engines.py:653)
>>> s1.components
{'H2O(aq)': 55.34457595832858, 'OH[-1]': 1.289102123028436e-07, 'H[+1]': 9.975725157654929e-08, 'Rh[+3]': 9.717653575367691e-09, 'O2(aq)': 1.008416691965259e-24, 'H2(aq)': 7.911616291116781e-35}

[rkingsbury]

The relevant logic is here. Basically, we update the Solution composition using values from the phreeqc instance rather than overwrite it.

Arguably, a more robust solution would be to just use the phreeqc composition as-is, but then do a check that the total concentrations of all elements (e.g., get_total_amount) in the original composition have not changed. That might impose a performance penalty compared to the current behavior though; which can be an issue if pyEQL is used as part of a transport model or something.

[rkingsbury]

Fixed by #282