alkalinity of seawater is too high for phreeqc engines

[jjstickel]

I think there is a bug with the alkalinity calculation in pyEQL.Solution. Example:

In [1]: import pyEQL as pyeql

In [2]: sol = pyeql.Solution.from_preset("seawater", engine="phreeqc2026", balance_charge="pH")
/Users/jonathan/Library/Python/3.12/lib/python/site-packages/pyEQL/solution.py:339: UserWarning: After initialization, the calculated solution pH of 8.220 does not match the specified pH of 8.099. This might be a result of erroneous input (e.g., specifying both pH and H+), or it can happen during from_dict / from_preset if you use a different engine than the one which generated the original dict.
  warnings.warn(

In [3]: sol.alkalinity # this looks OK
Out[3]: <Quantity(117.4375715148465, 'milligram / liter')>

In [4]: sol.equilibrate(atmosphere=True)

In [5]: sol.alkalinity # too high!
Out[5]: <Quantity(753.3507804556849, 'milligram / liter')>

Otherwise, the solution and its properties seem reasonable, e.g., pH = 8.2.

As a workaround, I'd like to retrieve the alkalinity calculated by the phreeqc2026 (or phreeqc) engine, but I could not find a way to do that. I do see some discussion of retrieving alkalinity and other properties in issue #310, but it is not clear what are the low-level calls to do this.

[rkingsbury]

Hmmm, thanks for flagging this @jjstickel . @SuixiongTay , can you investigate when you have a few minutes? I wonder if this is related to #398 .

@jjstickel , the underlying PHREEQC object is query-able (although it's generally intended for use only by developers), but I believe the correct code should be something like

sol.engine.ppsol.alkalinity

The engine.ppsol refers to the in-memory PHREEQC solution, and it's the entity that the pyEQL.Solution pulls properties from.

[jjstickel]

sol.engine.ppsol.alkalinity
Unfortunately, this attribute is not available in the current version of pyEQL. I submitted a PR #459 to bring in alkalinity from the phreeqc2026 engine.

[rkingsbury]

Digging further here, it's interesting that the discrepancy is not nearly as large with the native engine (default), which uses llnl.dat rather than phreeqc.dat. This results in a little bit of MgCl[+1] and CaCl[+1] that you don't get from the phreeqc2026 engine.

>>> s=Solution.from_preset('seawater', engine='native', balance_charge='pH')
>>> s.alkalinity
<Quantity(117.4375715148465, 'milligram / liter')>
>>> s.equilibrate()
>>> s.alkalinity
<Quantity(138.58643312800692, 'milligram / liter')>
>>> from pprint import pprint
>>> pprint(s.cations, sort_dicts=False)
{'Na[+1]': 0.4433714025177017,
 'Mg[+2]': 0.03989283543128503,
 'K[+1]': 0.00993963969389008,
 'Ca[+2]': 0.00912890816370233,
 'MgCl[+1]': 0.0048322059815775215,
 'CaCl[+1]': 0.0002349720175303199,
 'MgHCO3[+1]': 0.00017589217605201544,
 'Sr[+2]': 7.81912802211947e-05,
 'MgB(OH)4[+1]': 5.3366260877350515e-05,
 'CaHCO3[+1]': 3.2143754439079596e-05,
 'MgF[+1]': 1.5231616415800963e-05,
 'CaB(OH)4[+1]': 8.033054971734781e-06,
 'SrCl[+1]': 4.78669172386515e-06,
 'CaF[+1]': 5.74606853281073e-07,
 'CaOH[+1]': 1.0658588072772545e-07,
 'H[+1]': 5.442064470280126e-09,
 'SrF[+1]': 1.2010306904232672e-09,
 'SrOH[+1]': 2.8133147113171034e-10,
 'MgOH[+1]': 4.963854324274498e-12}
>>> pprint(s.anions, sort_dicts=False)
{'Cl[-1]': 0.5217477571135759,
 'SO4[-2]': 0.013070267924265954,
 'NaSO4[-1]': 0.006580128353233588,
 'HCO3[-1]': 0.0011945573924713067,
 'Br[-1]': 0.0008321477568366732,
 'KSO4[-1]': 0.0001686344960494305,
 'BO2[-1]': 5.204325165701761e-05,
 'F[-1]': 5.1346905890040505e-05,
 'CO3[-2]': 4.25613922170462e-05,
 'NaCO3[-1]': 1.265396061973774e-05,
 'OH[-1]': 3.302315125178726e-06,
 'HSO4[-1]': 1.4741764891422578e-09,
 'BH2(OF)2[-1]': 8.697358994517376e-15,
 'HF2[-1]': 2.7223126226384054e-15,
 'B2H5O6[-1]': 5.536726401768357e-18,
 'BOHF3[-1]': 5.557698941026352e-21,
 'BrO[-1]': 4.758462532332134e-24,
 'ClO[-1]': 1.9451356506408476e-25,
 'BF4[-1]': 5.518561133906346e-29,
 'Br3[-1]': 2.73868612789012e-29}

Comparing phreeqc2026 engine (so, identical speciation engine, but different activity coefficient methods):

>>> s=Solution.from_preset('seawater', engine='phreeqc2026', balance_charge='pH')
>>> s.alkalinity
<Quantity(117.4375715148465, 'milligram / liter')>
>>> s.equilibrate()
>>> s.alkalinity
<Quantity(750.1329716570397, 'milligram / liter')>
>>> pprint(s.cations, sort_dicts=False)
{'Na[+1]': 0.454076694579413,
 'Mg[+2]': 0.047568706901239065,
 'K[+1]': 0.009994141997974333,
 'Ca[+2]': 0.00956909884097742,
 'MgHCO3[+1]': 0.0002241702871149978,
 'Sr[+2]': 8.39764935415523e-05,
 'MgF[+1]': 2.990307255052956e-05,
 'MgOH[+1]': 2.002642507765606e-05,
 'CaHCO3[+1]': 1.1767367923649085e-05,
 'SrHCO3[+1]': 4.1161804789948005e-07,
 'CaOH[+1]': 1.519796167500615e-07,
 'H[+1]': 4.562072948936468e-09,
 'SrOH[+1]': 4.5220049834386825e-10,
 'CaHSO4[+1]': 2.231148160876166e-11}
>>> pprint(s.anions, sort_dicts=False)
{'Cl[-1]': 0.5442578559819381,
 'NaSO4[-1]': 0.013317073194106448,
 'SO4[-2]': 0.00921392866557388,
 'HCO3[-1]': 0.0012967454181522632,
 'Br[-1]': 0.0008395244921424564,
 'KSO4[-1]': 0.00018089820687298627,
 'Mg(SO4)2[-2]': 0.00012150409668995436,
 'BH2O3[-1]': 8.52909405909086e-05,
 'CO3[-2]': 5.651972067316964e-05,
 'F[-1]': 3.4861809834616796e-05,
 'OH[-1]': 4.67668745098709e-06,
 'BH3O3F[-1]': 4.368815670373759e-09,
 'HSO4[-1]': 7.702228593072517e-10,
 'BH2(OF)2[-1]': 3.532024945674748e-14,
 'HF2[-1]': 1.221783332859699e-14,
 'BOHF3[-1]': 2.922024164532674e-21,
 'BF4[-1]': 8.981395973376415e-28}

Note that I get simliar (but not identical) results if I disable the Stumm-Morgan "conservative ions" definition and rely solely on the weak acid/base definition.

TL;DR I believe this is a weakness in our implementation of get_alkalinity (and, arguably, the S-M definition) related to how to we handle complexes. Things like MgCl[+1] and CaCl[+1] don't get counted, and since the S-M definition sums over all the free cations and anions, the resulting alkalinity for speciated solutions is a bit of a wildcard - it depends on what fraction of the total cations vs. the total anions are complexed. It might be possible to work around this using total concentrations but I'd need to think about it more.

[rkingsbury]

@jjstickel the linked PR should fix the alkalinity discrepancy. As I suspected, pyEQL implementation of the S-M definition did not properly account for speciation (it dates to before speciation was available in the library). The PR modified it to use total element concentrations, resulting in a stable alkalinity before and after equilibration.

I'd still like to merge your changes from #459 that make ALK available via Solution.engine.ppsol, whenever you have time to narrow the scope.

[rkingsbury]

And thanks again for bringing this to my attention! Definitely a serious flaw.

After merging the PR, I'll plan to do a point release with the bug fix in case you need to be able to install from PyPi.

[jjstickel]

Sounds good. I'll try to spend some time on this issue this weekend.

[rkingsbury]

Following up on the remaining alkalinity discrepancy noted in the #459 discussion (pyEQL's Solution.alkalinity vs the phreeqc2026 engine's PhrqSol.get_alkalinity()). The two turn out to compute nearly the same quantity by two different routes, and the gap is fully explained by an exact charge-balance identity.

The two definitions

• PHREEQC ALK = titration (proton-condition) alkalinity of the speciated solution.
• Solution.alkalinity = conservative (Stumm) alkalinity from analytical totals: Σz·(conservative cations) − Σ|z|·(strong-acid anions).

By charge balance these are related exactly:

PHREEQC_ALK = conservative_alk(complete)  −  charge_balance

pyEQL differed from PHREEQC by exactly two omitted terms:

1. Fluoride (F⁻) was missing from the conservative strong-acid-anion set → pyEQL too high by [F⁻].
2. The residual charge imbalance of the input was not subtracted → pyEQL too high by the input's charge-balance residual.

Seawater (default preset)

	mg/L CaCO₃
pyEQL Solution.alkalinity	117.48
PHREEQC get_alkalinity()	106.97
gap	10.52
↳ fluoride term	3.39
↳ charge-imbalance term	7.13

PHREEQC_ALK == conservative_eq − F⁻ − charge_balance closes to 0.0000 mg/L.

Cross-checks (identity holds everywhere)

• Seawater with balance_charge="pH" → gap collapses to 3.39 mg/L, purely fluoride — this is the "≈114 vs 117" observation from the PR.
• Balanced NaHCO₃ (no F, charge-balanced) → the two agree to 0.03 mg/L, confirming they compute the same quantity when neither omitted term is active.
• Badly imbalanced fluoride-free brine → the entire ~1900 mg/L gap is the charge-balance term; the identity still closes to zero.

Resolution

• Fluoride is unambiguously a strong-acid anion and is now added to the conservative set in fix(alkalinity): count fluoride as a conservative strong-acid anion (#458) #467.
• The charge-imbalance term is left as-is by design: subtracting it would redefine alkalinity for non-electroneutral inputs, whereas the conservative and titration definitions already coincide for any charge-balanced solution. The docstring now notes that engine agreement should be evaluated on a charge-balanced solution.

[rkingsbury]

@jjstickel see Claude's diagnosis of the PHREEQC <-> pyEQL discrepancy above. I've opened #467 to partially address. I'm not yet sure what to make of the charge balance piece, and whether to incorporate that into the definition or not.