diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py index 9710ab95..431d1caf 100644 --- a/src/pyEQL/solution.py +++ b/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 diff --git a/tests/test_engine_phreeqc.py b/tests/test_engine_phreeqc.py index e4d2466f..ffe681d7 100644 --- a/tests/test_engine_phreeqc.py +++ b/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(): diff --git a/tests/test_engine_phreeqc2026.py b/tests/test_engine_phreeqc2026.py index 2eeb300b..a563ef96 100644 --- a/tests/test_engine_phreeqc2026.py +++ b/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(): diff --git a/tests/test_solution.py b/tests/test_solution.py index a8c3f050..4d0edecd 100644 --- a/tests/test_solution.py +++ b/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