Repeated calls to equilibrate() with 'auto' charge balancing cause pH to drift (phreeqc2026 engine)

[mdaeron]

This may be a new bug or a resurgence of #141, but I'm getting drifting pH when repeatedly equilibrating the same solution. I'm on a Mac M1 with pyEQL version a169325 (installed directly from github though uv:

from pyEQL import Solution
from pylab import *

initial_conditions = {
	'CO2(aq)': 1.8e-2,
	'H+': 1e-7,
}

S = Solution(
	initial_conditions,
	engine = 'phreeqc2026',
)

x = list(range(200))
y = []
for k in range(200):
	S.equilibrate()
	y.append(S.pH)

figure(figsize = (4,2.5))
subplots_adjust(bottom = 0.2)
plot(x,y)
xlabel('S.equilibrate() iterations')
ylabel('pH')
savefig('pH_drift.png')

yields:

[rkingsbury]

Thanks for reporting @mdaeron . I suspect this is a result of a lack of charge balancing. Unlike PHREEQC, pyEQL does not do automatic charge balancing by default, but you can tell it to do so via the balance_charge kwarg. So, can you let us know what happens if you modify your script to

S = Solution(
	initial_conditions,
	engine = 'phreeqc2026',
        balance_charge='auto'
)

?

[rkingsbury]

Tagging @YitongPan1 and @SuixiongTay as an FYI

[mdaeron]

The problem disappears if I add balance_charge = 'pH' to the Solution() call. As a newcomer I must say it was unexpected that this is not the default behavior. No big deal, though.

[mdaeron]

Ooh, but when I use balance_charge = 'auto' I get funky results:

[rkingsbury]

OK, I'm glad to hear (mostly)! I appreciate the feedback - the charge balancing logic and equilibrate method are both relatively new, so we can work on improving the documentation around those.

I'm intriuged by the oscillatory behavior you got with auto. I'm going to re-title and re-open this issue so we can investigate that. I think it's simply a case of pyEQL "overshooting" the charge balance slightly and somehow shifting the equilibrium, but we should investigate.

[rkingsbury]

Just a note to say this issue still exists as of the latest release (v1.4.0); however now if I set balance_charge='auto' I get a series of the following errors:

ERROR    [2026-02-17 05:33:08,341] [   ERROR] --- There is not enough OH[-1]   [solution.py](file:///home/ryan/miniforge3/envs/feynman/code/pyEQL/src/pyEQL/solution.py):[2389](file:///home/ryan/miniforge3/envs/feynman/code/pyEQL/src/pyEQL/solution.py#2389)
                             present to balance the charge. Try a different species.                               
                             (solution.py:2389)

Yet the pH still decreases as in your original post - no oscillations

On the other hand, with balance_charge=pH, the pH steadily increases, albeit at a slower rate

I get similar behavior with both engines

[rkingsbury]

Update: I believe the origin of this bug is that repeated calls to equilibrate cause the composition of the Solution to change slightly. In the native engine, this also causes the volume of the solution to change slightly. Perhaps more importantly, it seems to change the concentrations enough to eventually cause PHREEQC's equilibrium calcs to fail.

With the native engine, where volume depends on composition:

>>> from pyEQL import Solution
>>> s1 = Solution({"Na+": "0.8 mol", "Cl-": "0.4 mol", "SO4-2": "0.2 mol"}, pH=2, volume="0.1 L")
>>> for k in range(25):
...     s1.equilibrate()
...     print(s1.volume)
...
0.10869897975316335 l
0.10869959166109756 l
0.10870007393411929 l
0.10870045286961588 l
0.10870074986721853 l
0.10870098217771883 l
0.10870116359968529 l
0.10870130510099198 l
0.10870141535566738 l
0.10870150119623463 l
0.10870156798769554 l
0.1087016199322902 l
0.10870166031510985 l
0.1087016917003662 l
0.10870171608720794 l
0.10870173503280028 l
0.10870174974916792 l
0.10870176117916569 l
0.10870177005593672 l
0.10870177694936366 l
0.10870178230230763 l
0.10870178645884954 l
0.10870178968628015 l
0.10870179219221188 l
0.10870179413788728 l

With the phreeqc2026 engine, where solute volume is (currently) assumed to be zero - Note how much the ionic strength changes after a few calls.

>>> s1 = Solution({"Na+": "0.8 mol", "Cl-": "0.4 mol", "SO4-2": "0.2 mol"}, pH=2, volume="0.1 L", engine='phreeqc2026')
>>> for k in range(25):
...     print(s1.ionic_strength)
...     s1.equilibrate()
...
10.034621882974731 mol/kg
6.262738287439613 mol/kg
6.265875939510499 mol/kg
6.2703137665261846 mol/kg
6.276589871546618 mol/kg
6.285464424164794 mol/kg
6.298010707463637 mol/kg
6.315742970530024 mol/kg
6.340795564166596 mol/kg
6.376173327829149 mol/kg
6.426100854442841 mol/kg
6.496509388057287 mol/kg
6.595717395658866 mol/kg
6.73538978752815 mol/kg
6.931911653403657 mol/kg
7.208402342950813 mol/kg
7.597749937386956 mol/kg
8.147300714357069 mol/kg
8.926246479739316 mol/kg
10.03734921753428 mol/kg
11.63513233602117 mol/kg
13.951778508798968 mol/kg
[05/28/26 15:53:29] ERROR    [2026-05-28 15:53:29,480] [   ERROR] --- PHREEQC returned a total Element H  engines.py:407
                             concentration of 13.303613370540704 mol, which differs from the original
                             concentration of 12.623743301856896. This should never occur and indicates
                             an error in the PHREEQC database or calculation. (engines.py:407)
17.33000559819732 mol/kg
                    ERROR    [2026-05-28 15:53:29,488] [   ERROR] --- PHREEQC returned a total Element H  engines.py:407
                             concentration of 14.295725056094058 mol, which differs from the original
                             concentration of 13.303613370540704. This should never occur and indicates
                             an error in the PHREEQC database or calculation. (engines.py:407)
22.269202826464376 mol/kg
Traceback (most recent call last):
  File "/home/ryan/mambaforge/envs/feynman/code/pyEQL/src/pyEQL/engines.py", line 267, in _setup_ppsol
    ppsol = self.pp.add_solution(PHRQSol(d))
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/ryan/mambaforge/envs/feynman/code/pyEQL/src/pyEQL/phreeqc/core.py", line 140, in add_solution
    self()
  File "/home/ryan/mambaforge/envs/feynman/code/pyEQL/src/pyEQL/phreeqc/core.py", line 67, in __call__
    self.run_string(self._str)
RuntimeError: ERROR:               A(H2O) Activity of water has not converged.  Residual: 2.192802e-04

ERROR: Model failed to converge for initial solution  0.


During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 3, in <module>
  File "/home/ryan/mambaforge/envs/feynman/code/pyEQL/src/pyEQL/solution.py", line 1762, in equilibrate
    self.engine.equilibrate(self, atmosphere=atmosphere, solids=solids, gases=gases, **kwargs)
  File "/home/ryan/mambaforge/envs/feynman/code/pyEQL/src/pyEQL/engines.py", line 314, in equilibrate
    self._setup_ppsol(solution)
  File "/home/ryan/mambaforge/envs/feynman/code/pyEQL/src/pyEQL/engines.py", line 270, in _setup_ppsol
    raise ValueError(
ValueError: There is a problem with your input. The error message received from  phreeqc is:

 ERROR:               A(H2O) Activity of water has not converged.       Residual: 2.192802e-04

ERROR: Model failed to converge for initial solution  0.

 Check your input arguments, especially the composition dictionary, and try again.

[rkingsbury]

^ @vineetbansal can you take a look at my last post? I'm wondering if there is some mismatch in numerical precision when the composition gets handed back and forth between pyEQL native code and pyEQL-phreeqc that causes small but accumulative changes in the composition.

In principle, once we call equilibrate, repeated calls should no longer change anything.

[rkingsbury]

Note that I observed similar behavior with the legacy wrapper as well, so it's possible this is down to something internal in PHREEQC.