Enable retrieving alkalinity from phreeqc2026 engine

[jjstickel]

Summary

Major changes:

• creates infrastructure for Solution.engine.get_alkalinity()
• retrieve alkalinity for the phreeqc2026 engine
• use the engine's alkalinity, if available; otherwise fall-back to alkalinity calculated within pyEQL class

Todos

• review overall approach
• determine whether the changed test conditions are generally appropriate
• also fix pyEQL's alkalinity calculation(?) see alkalinity of seawater is too high for phreeqc engines #458

Checklist

• Google format doc strings added.
• Code linted with ruff. (For guidance in fixing rule violates, see rule list)
• Type annotations included. Check with mypy.
• Tests added for new features/fixes.
• I have run the tests locally and they passed.

Tip: Install pre-commit hooks to auto-check types and linting before every commit:

pip install -U pre-commit
pre-commit install

[jjstickel]

Disclosure: I used Claude to help, but the changes are few and make sense to me.

Sorry about lumping removal of docs/examples/.ipynb_checkpoints in with this PR. I can figure out how to unroll that if the PR is invited to proceed.

[rkingsbury]

Thanks very much @jjstickel ; this is a helpful contribution! I'd like to proceed but with a few structural changes:

1. If not too much trouble, can you split the .gitignore changes into another PR? That is a very helpful change; I'd just prefer to keep it separate
2. The changes to expose ALK to the underlying PHREEQC solution (pyEQL/phreeqc/core.py and solution.py) are very helpful; I'm happy with those as-is
3. I'm not (yet) ready to make get_alkalinity an engine method, so can you please revert the changes in pyEQL/engines.py and pyEQL/solution.py? In general, I we try to use engine-specific methods only for those quantities that do not have an unambiguous definition, and rely on pyEQL.Solution native methods for everything else.

Regarding 3 - you have highlighted an important discrepancy between the pyEQL.Solution alkalinity and the underlying PHREEQC one which I want to instead treat as a bug in get_alkalinity. We recently implemented the weak acid-base definition of alkalinity (see #299 and #398 ) and perhaps that introduced a bug we weren't aware of. OR perhaps it isn't appropriate to rely on the Stumm and Morgan definition as the default (which is our current approach). Maybe it's more appropriate to use the proton condition as default.

I welcome additional thoughts!

[jjstickel]

Thanks for the feedback. I'll work on these as I have time (likely later this week or next week).

[rkingsbury]

FYI, I (actually Claude Code) took care of the .gitignore issue in #460 and I updated your branch here, so that part is all done.

[jjstickel]

3. I'm not (yet) ready to make get_alkalinity an engine method, so can you please revert the changes in pyEQL/engines.py and pyEQL/solution.py? In general, I we try to use engine-specific methods only for those quantities that do not have an unambiguous definition, and rely on pyEQL.Solution native methods for everything else.

I added a keyword argument to select between returning the alkalinity calculated by the engine vs. by pyEQL, and set it to be pyEQL by default.

I work with geochemists who have been using PHREEQC for years and trust it a great deal. It is, and will continue to be, necessary for me in my work to report the alkalinity that PHREEQC calculates. Also, from the current bug in alkalinity, it seems that alkalinity calculations are not always straightforward, and so having a feature to calculate it a couple different ways could be useful. (I am not a geochemist, and so I am not able to suggest how to improve pyEQL's alkalinity calculation without performing some literature research).

Does the keyword approach work for you? If this all gets sorted and everyone is convinced that pyEQL is doing the correct calculation [1], then the keyword argument could be removed in the future.

[1] If the "correct" alkalinity calculation ends up being different from PHREEQC, even by a small fraction, it would be important to document why the calculations are different.

[rkingsbury]

I work with geochemists who have been using PHREEQC for years and trust it a great deal. It is, and will continue to be, necessary for me in my work to report the alkalinity that PHREEQC calculates. Also, from the current bug in alkalinity, it seems that alkalinity calculations are not always straightforward, and so having a feature to calculate it a couple different ways could be useful. (I am not a geochemist, and so I am not able to suggest how to improve pyEQL's alkalinity calculation without performing some literature research).

Yes, totally understand. I'm concerned that there is this big a difference between what PHREEQC vs. pyEQL are returning. Clearly something going on.

Does the keyword approach work for you? If this all gets sorted and everyone is convinced that pyEQL is doing the correct calculation [1], then the keyword argument could be removed in the future.

[1] If the "correct" alkalinity calculation ends up being different from PHREEQC, even by a small fraction, it would be important to document why the calculations are different.

I'm not opposed to a keyword argument if there turns out to be some structural difference in the calculations that's correct, but I still feel we shouldn't implement that (especially at the class level) until the bug is better understood. I've done a bit of digging (see the Issue thread) and this seems related to complexes and complex solutions (hence why it appears after an equilibrate() call).

How about this - let's contain your fix to the body of Solution.alkalinity itself, functionallydoing the same thing as alkalinity_calc='engine' and then falling back to the native pyEQL calcs if that doesn't work. That way we don't have to change the Solution call signature or modify engines.py. So the body of Solution.alkalinity can be something like

    @property
    def alkalinity(self) -> Quantity:
        try:
            if (s.engine.ppsol is None) or (s.components != sengine._stored_comp):
                s.engine._destroy_ppsol()
                s.engine._setup_ppsol(solution)
                alk_eq_per_kgw = s.engine.ppsol.get_alkalinity()
                kgw = s.engine.ppsol.get_kgw()
                vol_L = s.volume.to("L").magnitude
                alk_eq_per_L = alk_eq_per_kgw * kgw / vol_L
                return (ureg.Quantity(alk_eq_per_L, "mol/L") * EQUIV_WT_CACO3).to("mg/L")
        except ValueError:
            pass
       
      ... continue with the existing method as written
        

And you would revert the changes in engines.py and Solution.__init__.

I feel confident that we can find a way to fix the behavior of the built-in alkalinity method, so this will provide a temporary fix in a less intrusive way. What do you think?

[jjstickel]

That's great you were able to fix the internal alkalinity calculation in #462. However, the alkalinity calculated by phreeqc is still different by more than 1% for seawater (114 vs 117 mg/L).[1] I've made changes to this PR to use the engine's alkalinity for an internal check (but removing the code from engines.py, like you requested). You can let me know what you think about keeping this check or something like it.

Separately, is it helpful to keep the line,

self._stored_comp = solution.components.copy()

in pyEQL/engines.py? Claude suggested it as an efficiency improvement, but I don't know the base code well enough to know if it should be kept.

[1] Interestingly, both values also increase by about 1 mg/L after equilibration with air.

[rkingsbury]

That's great you were able to fix the internal alkalinity calculation in #462. However, the alkalinity calculated by phreeqc is still different by more than 1% for seawater (114 vs 117 mg/L).[1]

Interesting. At that magnitude, the difference could be something subtle like a difference in the precision of the molecular weights used in the two respective databases. In any case, once you've finished putting this infrastructure in place to compare the two directly, it'll be much easier to pin down.

I've made changes to this PR to use the engine's alkalinity for an internal check (but removing the code from engines.py, like you requested). You can let me know what you think about keeping this check or something like it.

Thank you for the edits. Yes, I'm OK with this approach at least for now; there's no harm in a warning. I requested several specific changes on the implementation to align with the norms in our repo, but nothing changes at a high level.

Separately, is it helpful to keep the line,

self._stored_comp = solution.components.copy()

in pyEQL/engines.py? Claude suggested it as an efficiency improvement, but I don't know the base code well enough to know if it should be kept.

Yes, I think so! See comment, but ideally I'd prefer it in a dedicated PR. Claude Code should be able to do that very easily.

Finally, you will need to make edits to two places in test/phreeqc/test_phreeqc.py to account for the fact that there is now an extra column (ALK) coming out of the phreeqc output, e.g.

>           assert set(props.keys()) == set(expected[solution_index].keys())
E           assert {'ALK', 'CELL...s', 'species'} == {'CELL_NO', '...s', 'species'}
E             
E             Extra items in the left set:
E             'ALK'
E             
E             Full diff:
E               {
E             +     'ALK',
E                   'CELL_NO',
E                   'OSMOTIC',
E                   "TOT['water']",
E                   'eq_species',
E                   'species',
E               }

and

>       assert phreeqc.get_selected_output_column_count() == 27
E       assert 28 == 27
E        +  where 28 = get_selected_output_column_count()
E        +    where get_selected_output_column_count = <pyEQL.phreeqc.core.Phreeqc object at 0x7f0b355dcb40>.get_selected_output_column_count

[rkingsbury]

Thank you @jjstickel ! I jumped in and updated test/phreeqc/test_phreeqc.py so those tests will pass - all the failures were expected based on the fact that your change expands the output that comes back from PHREEQC itself. That's very "under the hood" stuff so I don't expect you to know how to modify it.

I left in place the one remaining test failure, which I realized is a consequence of the way we detect the engine:

        # check against alkalinity provided by the engine
        if hasattr(self.engine, 'ppsol'):
            try:
                if (self.engine.ppsol is None) or (self.components\
                                                   != self.engine._stored_comp):
                    self.engine._destroy_ppsol()
                    self.engine._setup_ppsol(self)
>               alk_eq_per_kgw = self.engine.ppsol.get_alkalinity()
                                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
E               AttributeError: 'Solution' object has no attribute 'get_alkalinity'

The phreeqc2026 engine exposes get_alkalinity per your changes here, but the legacy phreeqc engine does not, yet both engines have the ppsol attribute. This test is failing b/c it uses the legacy phreeqc engine.

The fix is to change

if hasattr(self.engine, 'ppsol'):

to

if isinstance(self.engine, Phreeqc2026EOS):`

The NativeEOS inherits from Phreeqc2026EOS, so the isinstance check will capture both cases.

If you agree with that change, I'm ready to merge once it lands.

[jjstickel]

Thanks for correcting the tests.

if isinstance(self.engine, Phreeqc2026EOS):

This didn't work, maybe because PhreeqcEOS also inherits from Phreeqc2026EOS?

Instead, I changed the try-except block to also catch the AttributeError. Does that work for you?

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 87.78%. Comparing base (ade469e) to head (845653a).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #459      +/-   ##
==========================================
+ Coverage   87.67%   87.78%   +0.10%     
==========================================
  Files          14       14              
  Lines        1972     1989      +17     
  Branches      344      347       +3     
==========================================
+ Hits         1729     1746      +17     
  Misses        192      192              
  Partials       51       51              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

All changes addressed

[rkingsbury]

Thank you @jjstickel ! Now I can let Claude Code dig into the PHREEQC vs. pyEQL discrepancy.