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.