alkalinity: add definition based on weak acids and bases

src/pyEQL/solution.py

@@ -834,8 +834,19 @@ def alkalinity(self) -> Quantity:
 
             Where :math:`C_{B}` and :math:`C_{A}` are conservative cations and anions, respectively
             (i.e. ions that do not participate in acid-base reactions), and :math:`z_{i}` is their signed charge.
-            In this method, the set of conservative cations is all Group I and Group II cations, and the
-            conservative anions are all the anions of strong acids.
+            When conservative cations (Group I and II cations) or strong base anions are present, the alkalinity is calculated according to[stm]_
+
+            .. math::   Alk = \sum_{i} z_{i} C_{B} + \sum_{i} z_{i} C_{A}
+
+            Where :math:`C_{B}` and :math:`C_{A}` are conservative cations and strong base anions, respectively  (i.e. ions that do not participate in acid-base reactions), and :math:`z_{i}` is their signed charge.
+
+            Alternatively, if those species are not present, then alkalinity is calculated based on the concentrations of weak acid and base species according to [stm]_
+
+            .. math:: Alk = -\sum_{i} z_{i} C_{i}
+
+            Where :math:`C_i` is the molar concentration of species i, and :math:`z_i` is its charge.
+
+            The summation should extend over all weak inorganic species that can participate in acid-base reactions. In this method, we consider HCO3[-1], CO3[-2], H2PO4[-1], HPO4[-2], PO4[-3], HS[-1], S[-2], H3SiO4[-1], H2SiO4[-2], B(OH)4[-1], NH3(aq), OH[-1], and H[+1] as the relevant weak acid/base species, while organics are excluded.
 
         References:
             .. [stm] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, pp 165. Wiley Interscience, 1996.
@@ -857,14 +868,44 @@ def alkalinity(self) -> Quantity:
             "Ba[+2]",
             "Ra[+2]",
         }
-        acid_anions = {"Cl[-1]", "Br[-1]", "I[-1]", "SO4[-2]", "NO3[-1]", "ClO4[-1]", "ClO3[-1]"}
+        acid_anions = {
+            "Cl[-1]",
+            "Br[-1]",
+            "I[-1]",
+            "SO4[-2]",
+            "NO3[-1]",
+            "ClO4[-1]",
+            "ClO3[-1]",
+        }
+
+        weak_species = {
+            "HCO3[-1]",
+            "CO3[-2]",
+            "H2PO4[-1]",
+            "HPO4[-2]",
+            "PO4[-3]",
+            "HS[-1]",
+            "S[-2]",
+            "H3SiO4[-1]",
+            "H2SiO4[-2]",
+            "B(OH)4[-1]",
+            "NH3(aq)",
+            "OH[-1]",
+            "H[+1]",
+        }  # Note that organics are excluded
+
+        conservative_species = base_cations.union(acid_anions)
+        # check presence of conservative cations or strong base anions
+        conservative_def = any(item in conservative_species for item in self.components)
 
         for item in self.components:
-            if item in base_cations.union(acid_anions):
-                z = self.get_property(item, "charge")
-                alkalinity += self.get_amount(item, "mol/L") * z
+            if item in conservative_species:
+                # Conservative cations and strong base anions
+                alkalinity += self.get_amount(item, "eq/L")
+            elif item in weak_species and not conservative_def:
+                # Weak acid/base species, exclude organics
+                alkalinity += self.get_amount(item, "eq/L") * (-1)
 
-        # convert the alkalinity to mg/L as CaCO3
         return (alkalinity * EQUIV_WT_CACO3).to("mg/L")
 
     @property

tests/test_engine_phreeqc.py

@@ -373,19 +373,22 @@ def test_equilibrate_logC_pH_carbonate_13():
     assert np.isclose(solution.get_amount("CO3-2", "mol/kg").magnitude, 9.979e-4, atol=1e-5)
 
 
-@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()
 
-    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
+    HCO3 = solution.get_amount("HCO3-", "mol/L").magnitude
+    CO3 = solution.get_amount("CO3-2", "mol/L").magnitude
+    OH = solution.get_amount("OH-", "mol/L").magnitude
+    H = solution.get_amount("H+", "mol/L").magnitude
 
     # Alkalinity calculated from the excess of negative charges from weak acids
     calculated_alk = HCO3 + 2 * CO3 + OH - H
-    assert solution.alkalinity.to("mg/L").magnitude == pytest.approx(calculated_alk, abs=0.001)
+    # Convert alkalinity from mol/L to mg/L as CaCO3
+    EQUIV_WT_CACO3 = 100.09 / 2  # g/eq
+    calculated_alk_mg_L = calculated_alk * EQUIV_WT_CACO3 * 1000
+
+    assert solution.alkalinity.magnitude == pytest.approx(calculated_alk_mg_L, abs=1e-8)
 
 
 def test_equilibrate_2L():

tests/test_engine_phreeqc2026.py

@@ -378,19 +378,22 @@ def test_equilibrate_logC_pH_carbonate_13():
     assert np.isclose(solution.get_amount("CO3-2", "mol/kg").magnitude, 9.979e-4, atol=1e-5)
 
 
-@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="phreeqc2026")
     solution.equilibrate()
 
-    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
+    HCO3 = solution.get_amount("HCO3-", "mol/L").magnitude
+    CO3 = solution.get_amount("CO3-2", "mol/L").magnitude
+    OH = solution.get_amount("OH-", "mol/L").magnitude
+    H = solution.get_amount("H+", "mol/L").magnitude
 
     # Alkalinity calculated from the excess of negative charges from weak acids
     calculated_alk = HCO3 + 2 * CO3 + OH - H
-    assert solution.alkalinity.to("mg/L").magnitude == pytest.approx(calculated_alk, abs=0.001)
+    # Convert alkalinity from mol/L to mg/L as CaCO3
+    EQUIV_WT_CACO3 = 100.09 / 2  # g/eq
+    calculated_alk_mg_L = calculated_alk * EQUIV_WT_CACO3 * 1000
+
+    assert solution.alkalinity.magnitude == pytest.approx(calculated_alk_mg_L, abs=1e-8)
 
 
 def test_equilibrate_2L():

tests/test_solution.py

@@ -124,6 +124,36 @@ def s8():
     )
 
 
+@pytest.fixture
+def s9():
+    # weak acid for alkalinity with pH variation
+    return Solution(
+        [
+            ["H3PO4(aq)", "1 mM"],  # no contribution to alk or hardness
+            ["H2PO4-", "1 mM"],  # -1 meq/L
+            ["HPO4-2", "1 mM"],  # -2 meq/L
+        ],
+        volume="1 L",
+        pH=3,
+    )
+
+
+@pytest.fixture
+def s10():
+    # both conservative cation and weak acid for alkalinity with pH variation
+    return Solution(
+        [
+            ["Ca+2", "1 mM"],  # 2 meq/L
+            ["Na+", "1 mM"],  # 1 meq/L
+            ["H3SiO4-", "1 mM"],  # -1 meq/L
+            ["H2SiO4-2", "1 mM"],  # -2 meq/L
+            ["SiO2(aq)", "1 mM"],  # no contribution to alk or hardness
+        ],
+        volume="1 L",
+        pH=4,
+    )
+
+
 def test_empty_solution():
     # create an empty solution
     s1 = Solution(database=None)
@@ -384,7 +414,7 @@ def test_water_stability_reducing(s8, caplog):
     assert any("Hydrogen evolution may occur" in r.message for r in caplog.records)
 
 
-def test_alkalinity_hardness(s3, s5, s6):
+def test_alkalinity_hardness(s3, s5, s6, s9, s10):
     assert np.isclose(s3.hardness, 0)
     assert np.isclose(s3.alkalinity, 0)
 
@@ -394,6 +424,12 @@ def test_alkalinity_hardness(s3, s5, s6):
     assert np.isclose(s6.alkalinity.magnitude, -5900, rtol=0.005)
     assert np.isclose(s6.hardness.magnitude, 600, rtol=0.005)
 
+    assert np.isclose(s9.alkalinity.magnitude, 100.09, rtol=0.005)
+    assert np.isclose(s9.hardness.magnitude, 0, rtol=0.005)
+
+    assert np.isclose(s10.alkalinity.magnitude, 150.135, rtol=0.005)
+    assert np.isclose(s10.hardness.magnitude, 100.09, rtol=0.005)
+
 
 def test_pressure_temperature(s5):
     orig_V = s5.volume