Diffusion Coefficient from Phreeqc

src/pyEQL/engines.py

@@ -1128,3 +1128,23 @@ def equilibrate(
         # weights for water. Since water amount is passed to PHREEQC in kg but returned in moles, each
         # call to equilibrate can thus result in a slight change in the Solution mass.
         solution.components[solution.solvent] = orig_solvent_moles
+
+    def get_diffusion_coefficient(self, solution: "solution.Solution", solute: str) -> ureg.Quantity:
+        if (self.ppsol is None) or (solution.components != self._stored_comp):
+            self._destroy_ppsol()
+            self._setup_ppsol(solution)
+
+        # translate the species into keys that phreeqc will understand
+        k = standardize_formula(solute)
+        if "[" in k:
+            spl = k.split("[")
+            el = spl[0]
+            chg = spl[1].split("]")[0]
+            if chg[-1] == "1":
+                chg = chg[0]  # just pass + or -, not +1 / -1
+            k = el + chg
+
+        diffusion_coefficient = self.ppsol.get_diffusion_coefficient(k)
+        if diffusion_coefficient is None:
+            return None
+        return ureg.Quantity(diffusion_coefficient, "meter ** 2 / s")

src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py

@@ -13,7 +13,7 @@
     "TOT['water']",
     "OSMOTIC",
 )
-SPECIES_PROPS = ("MOL", "ACT")
+SPECIES_PROPS = ("MOL", "ACT", "DIFF_C")
 EQ_SPECIES_PROPS = ("SI",)
 
 

src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py

@@ -20,9 +20,11 @@ def _get_calculated_props(self):
 
     def _get_calculated_prop(self, which, species: str | None = None, eq_species: str | None = None):
         if species is not None:
-            return self._calculated_props["species"][species][which]
+            species_props = self._calculated_props["species"].get(species)
+            return species_props[which] if species_props is not None else None
         if eq_species is not None:
-            return self._calculated_props["eq_species"][eq_species][which]
+            eq_species_props = self._calculated_props["eq_species"].get(species)
+            return eq_species_props[which] if eq_species_props is not None else None
         return self._calculated_props[which]
 
     def get_activity(self, species) -> float:
@@ -31,6 +33,9 @@ def get_activity(self, species) -> float:
     def get_molality(self, species) -> float:
         return self._get_calculated_prop("MOL", species=species)
 
+    def get_diffusion_coefficient(self, species) -> float:
+        return self._get_calculated_prop("DIFF_C", species=species)
+
     def get_osmotic_coefficient(self) -> float:
         return self._get_calculated_prop("OSMOTIC")
 

src/pyEQL/solution.py

@@ -2071,7 +2071,7 @@ def _get_property(self, solute: str, name: str) -> Any | None:
             return ureg.Quantity(val)
         return None
 
-    def get_transport_number(self, solute: str) -> Quantity:
+    def get_transport_number(self, solute: str, use_engine: bool = False) -> Quantity:
         r"""Calculate the transport number of the solute in the solution.
 
         Args:
@@ -2114,7 +2114,7 @@ def get_transport_number(self, solute: str) -> Quantity:
             # cancels out
             # using species amounts in mol is equivalent to using concentrations in mol/L
             # since there is only one solution volume, and it's much faster.
-            term = self.get_molar_conductivity(item).magnitude * mol
+            term = self.get_molar_conductivity(item, use_engine=use_engine).magnitude * mol
 
             if item == solute:
                 numerator = term
@@ -2123,7 +2123,7 @@ def get_transport_number(self, solute: str) -> Quantity:
 
         return ureg.Quantity(numerator / denominator, "dimensionless")
 
-    def _get_molar_conductivity(self, solute: str) -> Quantity:
+    def _get_molar_conductivity(self, solute: str, use_engine: bool = False) -> Quantity:
         r"""
         Calculate the molar (equivalent) conductivity for a solute.
 
@@ -2160,7 +2160,7 @@ def _get_molar_conductivity(self, solute: str) -> Quantity:
         See Also:
             :py:meth:`get_diffusion_coefficient`
         """
-        D = self.get_diffusion_coefficient(solute)
+        D = self.get_diffusion_coefficient(solute, use_engine=use_engine)
 
         if D != 0:
             molar_cond = (
@@ -2173,14 +2173,18 @@ def _get_molar_conductivity(self, solute: str) -> Quantity:
 
         return molar_cond.to("mS / cm / (mol/L)")
 
-    def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = True) -> Quantity:
+    def _get_diffusion_coefficient(
+        self, solute: str, activity_correction: bool = True, use_engine: bool = False
+    ) -> Quantity:
         r"""
         Get the **temperature-adjusted** diffusion coefficient of a solute.
 
         Args:
             solute: the solute for which to retrieve the diffusion coefficient.
             activity_correction: If True (default), adjusts the diffusion coefficient for the effects of ionic
                 strength using a model from Ref 2.
+            use_engine: Whether to use the underlying phreeqc engine to determine diffusion coefficient.
+                Only the 'pyeql' engine is supported.
 
         Notes:
             This method is equivalent to self.get_property(solute, "transport.diffusion_coefficient")
@@ -2218,8 +2222,20 @@ def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = Tr
             pyEQL.activity_correction._debye_parameter_activity
 
         """
-        D = self.get_property(solute, "transport.diffusion_coefficient")
         rform = standardize_formula(solute)
+
+        if use_engine:
+            assert self._engine == "pyeql", "Only supported for pyeql engine"
+            D = self.engine.get_diffusion_coefficient(self, solute)
+            if D is None or D.magnitude == 0:
+                self.logger.warning(
+                    f"Diffusion coefficient not found for species {rform}. Using default value of "
+                    f"{self.default_diffusion_coeff} m**2/s."
+                )
+                return ureg.Quantity(self.default_diffusion_coeff, "m**2/s")
+            return D
+
+        D = self.get_property(solute, "transport.diffusion_coefficient")
         if D is None or D.magnitude == 0:
             self.logger.warning(
                 f"Diffusion coefficient not found for species {rform}. Using default value of "

tests/phreeqc/test_phreeqc.py

@@ -406,7 +406,7 @@ def test_add_solution_input():
         5 PUNCH CELL_NO, TOT['water'], OSMOTIC, EOL_NOTAB$
         10 t = SYS("aq", count, name$, type$, moles)
         20 FOR i = 1 to count
-        30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i))
+        30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)), DIFF_C(name$(i))
         40 NEXT i
         50 PUNCH EOL$
         60 p = SYS("phases", count, name$, type$, moles)
@@ -454,9 +454,9 @@ def test_add_solution_output():
     # heading + soln 0 data + soln 1 data (regardless of whether we have -headings in USER_PUNCH or not)
     assert phreeqc.get_selected_output_row_count() == 3
     # <cell_no>, <tot_water>, <osmotic>, "\n",
-    # +[<name>, <molality>, <activity>] repeated for each (4) species
+    # +[<name>, <molality>, <activity>, <diff_c>] repeated for each (4) species
     # +"\n" + [<equilibrium_species>, <si>] repeated for each (3) equilibrium species
-    assert phreeqc.get_selected_output_column_count() == 23
+    assert phreeqc.get_selected_output_column_count() == 27
 
 
 def test_kgw():

tests/test_solution_pyeql_engine.py

@@ -0,0 +1,92 @@
+"""
+pyEQL volume and concentration methods test suite
+=================================================
+
+This file contains tests for the volume and concentration-related methods
+used by pyEQL's Solution class
+"""
+
+import numpy as np
+import pytest
+
+from pyEQL import Solution
+
+
+@pytest.fixture
+def s1():
+    return Solution(volume="2 L", engine="pyeql")
+
+
+@pytest.fixture
+def s2():
+    return Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], volume="2 L", engine="pyeql")
+
+
+def test_diffusion_transport(s1, s2):
+    # test ionic strength adjustment
+    # assert s1.get_diffusion_coefficient("H+", use_engine=True) > s2.get_diffusion_coefficient("H+", use_engine=True)
+    # for Na+, d=122, a1=1.52, a2=3.7, A=1.173802/2.303 at 25 DegC, B = 3.2843078+10
+    factor = np.exp(
+        -1.52
+        * 1.173802
+        / 2.303
+        * 1
+        * np.sqrt(s2.ionic_strength.magnitude)
+        / (1 + 3.2843078e10 * np.sqrt(s2.ionic_strength.magnitude) * 3.7 / (1 + s2.ionic_strength.magnitude**0.75))
+    )
+    assert np.isclose(
+        factor * s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude,
+        s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude,
+        atol=5e-11,
+    )
+    s_dilute = Solution({"Na+": "1 mmol/L", "Cl-": "1 mmol/L"}, engine="pyeql")
+    assert np.isclose(
+        s_dilute.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude,
+        1.334e-9,
+        atol=1e-11,
+    )
+    assert np.isclose(s_dilute.get_transport_number("Na+", use_engine=True), 0.396, atol=1e-3)
+    assert np.isclose(s_dilute.get_transport_number("Cl-", use_engine=True), 0.604, atol=1e-3)
+
+    # test setting a default value
+    s2.default_diffusion_coeff = 0
+    assert s2.get_diffusion_coefficient("Cs+", use_engine=True).magnitude == 0
+    s2.default_diffusion_coeff = 1e-9
+    assert s2.get_diffusion_coefficient("Cs+", activity_correction=False, use_engine=True).magnitude == 1e-9
+    s2.default_diffusion_coeff = 0
+    assert s2.get_diffusion_coefficient("Cs+", activity_correction=True, use_engine=True).magnitude < 1e-9
+    d25 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude
+    nu25 = s2.water_substance.nu
+    s2.temperature = "40 degC"
+    d40 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude
+    nu40 = s2.water_substance.nu
+    assert np.isclose(
+        d40,
+        d25 * np.exp(122 / (273.15 + 40) - 122 / 298.15) * (nu25 / nu40),
+        atol=5e-10,
+    )
+
+    # test correction factors for concentration, as per Appelo 2017 Fig 5
+    D1 = (
+        Solution({"Na+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql")
+        .get_diffusion_coefficient("Na+", use_engine=True)
+        .magnitude
+    )
+    D2 = (
+        Solution({"Na+": "1.7 mol/kg", "Cl-": "1.7 mol/kg"}, engine="pyeql")
+        .get_diffusion_coefficient("Na+", use_engine=True)
+        .magnitude
+    )
+    assert np.isclose(D2 / D1, 1, atol=1e-2)
+
+    D1 = (
+        Solution({"K+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql")
+        .get_diffusion_coefficient("K+", use_engine=True)
+        .magnitude
+    )
+    D2 = (
+        Solution({"K+": "0.5 mol/kg", "Cl-": "0.5 mol/kg"}, engine="pyeql")
+        .get_diffusion_coefficient("K+", use_engine=True)
+        .magnitude
+    )
+    assert np.isclose(D2 / D1, 1, atol=1e-2)