pyEQL alkalinity incorrect for pure carbonate solutions

[SuixiongTay]

Summary

An enhancement can be made on how pyEQL handles alkalinity. The alkalinity() property considers only a limited set of ions. It excludes important alkalinity contributors like:

• Carbonate species: [HCO3-], [CO3-2]
• Phosphorus species: [HPO4−2], [PO4−3], [H2PO4-]
• Possible sulfur species: [HS−], [S2−], [HSO4-]

References

• Related to Adding tests for solid-liquid equilibrium in equilibrate()  #295
• Stumm & Morgan (1996) Aquatic Chemistry, 3rd ed.

[vineetbansal]

The xfail test that demonstrates this is in tests:

@pytest.mark.xfail(strict=True, reason="alkalinity discrepancy needs to be investigated")
def test_alkalinity():
    solution = Solution({"CO2(aq)": "0.001 mol/L"}, pH=7, volume="1 L", engine="phreeqc")
    solution.equilibrate()
    alk = solution.alkalinity
    # Total alkalinity
    HCO3 = solution.get_amount("HCO3-", "mg/L").magnitude
    CO3 = solution.get_amount("CO3-2", "mg/L").magnitude
    OH = solution.get_amount("OH-", "mg/L").magnitude
    H = solution.get_amount("H+", "mg/L").magnitude
    # Alkalinity calculated from the excess of negative charges from weak acids
    total_alk = HCO3 + 2 * CO3 + OH - H
    assert alk.to("mg/L").magnitude == pytest.approx(total_alk, abs=0.001)

[rkingsbury]

Thanks for flagging this @SuixiongTay . pyEQL uses the "alternate" definition of alkalinity explained in Stumm & Morgan (see comment in the method docstring, which is based on electroneutrality of conservative cations/anions. So it's not accurate to say pyEQL "excludes" weak acid species from alkalinity; it's just that it employs an alternative route to getting them.

I've just re-read both definitions, and I think the only situation in which this "alternative" definition does not work is the case you've identified - where there is nothing in the water except CO2 from the atmosphere. So, I think you've identified an edge case.

One possible fix would be to replace the current definition with the other definition, which is based on a listing of all possible weak acid/base species in solution. My hesitation about this approach is that I think it would be harder to enumerate all possible weak acids/bases than it is to enumerate all the possible conservative cations/anions (ions that do not participate in weak acid-base equilibria).

An alternative, more limited fix would be to just catch the specific case you've identified - if none of the listed base_cations or acid_anions is present in the Solution, then calculate the alkalinity directly from the carbonate species as you've done in your example.

As an aside - on reviewing the code for alkalinity I see a few opportunities to improve its numerical efficiency, so if we decide to make changes to it, let's discuss that aspect as well.