From 3cb64b3b13a494945454897aede7989920a8133f Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Wed, 13 May 2026 11:13:33 -0400
Subject: [PATCH 01/20] get_sat_idx

---
 src/pyEQL/phreeqc/solution.py |  5 ++++
 src/pyEQL/solution.py         | 54 +++++++++++++++++++++++++++++++++--
 tests/test_solution.py        |  9 ++++++
 3 files changed, 65 insertions(+), 3 deletions(-)

diff --git a/src/pyEQL/phreeqc/solution.py b/src/pyEQL/phreeqc/solution.py
index fb74aee9..fa3c9a97 100644
--- a/src/pyEQL/phreeqc/solution.py
+++ b/src/pyEQL/phreeqc/solution.py
@@ -55,6 +55,11 @@ def equalize(self, phases: list[str], to_si: list[float], in_phase: list[float])
     def si(self, eq_species) -> float:
         return self._get_calculated_prop("SI", eq_species=eq_species)
 
+    @property
+    def phases(self) -> dict[str, float]:
+        eq_species = self._get_calculated_prop("eq_species") or {}
+        return {phase: props["SI"] for phase, props in eq_species.items()}
+
     """
     The following properties are somewhat redundant, but included in here
     so we can act as a drop-in replacement for PhreeqPython as far as its
diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py
index 7e57225f..b8d46c8c 100644
--- a/src/pyEQL/solution.py
+++ b/src/pyEQL/solution.py
@@ -17,6 +17,8 @@
 from typing import Any, Literal
 
 import numpy as np
+import pandas as pd
+import plotly.express as px
 from maggma.stores import JSONStore, Store
 from monty.dev import deprecated
 from monty.json import MontyDecoder, MSONable
@@ -777,8 +779,8 @@ def charge_balance(self) -> float:
         r"""
         Return the signed charge balance of the solution, positive or negative.
 
-        Return the signed charge balance of the solution, positive or negative. The charge balance represents the net electric charge 
-        of the solution and SHOULD equal zero at all times, but due to numerical errors will usually have a small nonzero value. 
+        Return the signed charge balance of the solution, positive or negative. The charge balance represents the net electric charge
+        of the solution and SHOULD equal zero at all times, but due to numerical errors will usually have a small nonzero value.
         Positive values indicate excess cationic charge, while negative values indivate excess anionic charge. It is calculated according to:
 
         .. math:: CB = \sum_i C_i z_i
@@ -1699,7 +1701,7 @@ def equilibrate(
                 typically not considered due to its low solubility and limited
                 impact on aqueous speciation.
             solids:
-                A list of solids used to achieve liquid–solid equilibrium. Each
+                A list of solids used to achieve liquid-solid equilibrium. Each
                 solid in this list should be the name of a mineral phase present
                 in the Phreeqc database (e.g. "Calcite"). We assume a target
                 saturation index of 0 and an infinite amount of material.
@@ -2784,3 +2786,49 @@ def add_solvent(self, formula: str, amount: str):  # pragma: no cover
         mw = self.get_property(formula, "molecular_weight")
         target_mol = quantity.to("moles", "chem", mw=mw, volume=self.volume, solvent_mass=self.solvent_mass)
         self.components[formula] = target_mol.to("moles").magnitude
+
+    def get_saturation_index(self, phases=None, get_plot=None) -> float:
+        """
+        Calculate the saturation index of a solute in the solution.
+        """
+
+        engine = self.engine
+
+        if not hasattr(engine, "ppsol"):
+            raise NotImplementedError(f"Engine {type(engine).__name__} does not support saturation index calculations.")
+
+        # caching method from Phrqsol
+        if (engine.ppsol is None) or (self.components != engine._stored_comp):
+            engine._destroy_ppsol()
+            engine._setup_ppsol(self)
+
+        ppsol = engine.ppsol
+
+        phases = list(ppsol.phases.keys())
+        eq_species_dict = {phase: ppsol.si(phase) for phase in phases}
+
+        sorted_eq_species_dict = dict(sorted(eq_species_dict.items(), key=lambda item: item[1], reverse=True))
+
+        if get_plot:
+            df = pd.DataFrame(
+                {"species": list(sorted_eq_species_dict.keys()), "si": list(sorted_eq_species_dict.values())}
+            )
+
+            fig = px.bar(
+                df,
+                x="species",
+                y="si",
+                labels={"species": "Mineral Phase", "si": "Saturation Index"},
+                color="si",
+                color_continuous_scale="Mint",
+            )
+
+            fig.update_layout(
+                xaxis_tickangle=-45,
+                template="plotly_white",
+                height=500,
+            )
+
+            fig.show()
+
+        return sorted_eq_species_dict
diff --git a/tests/test_solution.py b/tests/test_solution.py
index 1de69699..eac7371c 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1155,3 +1155,12 @@ def test_should_use_major_salt_molar_volume_to_calculate_solute_volume_when_para
                 "OH-", "size.molar_volume"
             )
         assert solute_volume_without_protons_and_hydroxide.m == expected_solute_volume.m
+
+
+class TestSaturationIndex:
+    @staticmethod
+    @pytest.mark.parametrize("engine", ["native", "phreeqc", "phreeqc2026"])
+    def test_halite_si_positive(engine):
+        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
+        si = solution.get_saturation_index()
+        assert si["Halite"] > 0

From 0a334fd67972e0c73827e925bf35333152acde85 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Wed, 13 May 2026 11:25:46 -0400
Subject: [PATCH 02/20] correct halite naming

---
 tests/test_solution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index eac7371c..f760ada9 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1160,7 +1160,7 @@ def test_should_use_major_salt_molar_volume_to_calculate_solute_volume_when_para
 class TestSaturationIndex:
     @staticmethod
     @pytest.mark.parametrize("engine", ["native", "phreeqc", "phreeqc2026"])
-    def test_halite_si_positive(engine):
+    def test_halite_si(engine):
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
         si = solution.get_saturation_index()
         assert si["Halite"] > 0

From c782684e613c981d9f4d9463c1e219b993e86531 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Wed, 13 May 2026 16:04:05 -0400
Subject: [PATCH 03/20] add monkeypatch for plt

---
 tests/test_solution.py | 32 +++++++++++++++-----------------
 1 file changed, 15 insertions(+), 17 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 27a16c38..fd203a84 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -12,6 +12,7 @@
 from itertools import zip_longest
 
 import numpy as np
+import plotly.graph_objs as go
 import pytest
 import yaml
 from monty.serialization import dumpfn, loadfn
@@ -90,7 +91,7 @@ def s6_Ca():
             ["CO3-2", "6 mM"],  # no contribution to alk or hardness
             ["SO4-2", "60 mM"],  # -120 meq/L
             ["Br-", "20 mM"],  # -20 meq/L
-        ],  
+        ],
         volume="1 L",
         balance_charge="Ca+2",
     )
@@ -98,29 +99,29 @@ def s6_Ca():
 
 @pytest.fixture
 def s7():
-    # unstability solution in specific redox and pH combinations
+    # instability solution in specific redox and pH combinations
     return Solution(
         [
             ["Na+", "100 mM"],  # 100 meq/L
             ["Cl-", "100 mM"],  # -100 meq/L
         ],
         volume="1 L",
-        pH = 20,
-        pE = 1,
+        pH=20,
+        pE=1,
     )
 
 
 @pytest.fixture
 def s8():
-    # unstability solution in specific redox and pH combinations
+    # instability solution in specific redox and pH combinations
     return Solution(
         [
             ["Na+", "100 mM"],  # 100 meq/L
             ["Cl-", "100 mM"],  # -100 meq/L
         ],
         volume="1 L",
-        pH = 0,
-        pE = -10,
+        pH=0,
+        pE=-10,
     )
 
 
@@ -368,26 +369,20 @@ def test_charge_balance(s3, s5, s5_pH, s6, s6_Ca):
 
     with pytest.raises(ValueError, match=r"Charge balancing species Zr\[\+4\] was not found"):
         s = Solution({"Na+": "2 mM", "Cl-": "2 mM"}, balance_charge="Zr[+4]")
-        
+
 
 def test_water_stability_oxidizing(s7, caplog):
     with caplog.at_level(logging.WARNING, logger=s7.logger.name):
         s7._check_water_stability()
 
-    assert any(
-        "Oxygen evolution may occur" in message
-        for message in caplog.messages
-    )
+    assert any("Oxygen evolution may occur" in message for message in caplog.messages)
 
 
 def test_water_stability_reducing(s8, caplog):
     with caplog.at_level(logging.WARNING, logger=s8.logger.name):
         s8._check_water_stability()
 
-    assert any(
-        "Hydrogen evolution may occur" in r.message
-        for r in caplog.records
-    )
+    assert any("Hydrogen evolution may occur" in r.message for r in caplog.records)
 
 
 def test_alkalinity_hardness(s3, s5, s6):
@@ -1208,7 +1203,10 @@ def test_should_use_major_salt_molar_volume_to_calculate_solute_volume_when_para
 class TestSaturationIndex:
     @staticmethod
     @pytest.mark.parametrize("engine", ["native", "phreeqc", "phreeqc2026"])
-    def test_halite_si(engine):
+    def test_halite_si(engine, monkeypatch):
+        monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
         si = solution.get_saturation_index()
         assert si["Halite"] > 0
+        si_plot = solution.get_saturation_index(get_plot=True)
+        assert isinstance(si_plot, dict)

From f17b807b2bf91c5310619b22ed92c5bd232b0ca6 Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Wed, 13 May 2026 18:19:48 -0400
Subject: [PATCH 04/20] Update test_solution.py

---
 tests/test_solution.py | 49 ++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 47 insertions(+), 2 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index fd203a84..5352a868 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1198,7 +1198,7 @@ def test_should_use_major_salt_molar_volume_to_calculate_solute_volume_when_para
                 "OH-", "size.molar_volume"
             )
         assert solute_volume_without_protons_and_hydroxide.m == expected_solute_volume.m
-
+        
 
 class TestSaturationIndex:
     @staticmethod
@@ -1207,6 +1207,51 @@ def test_halite_si(engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
         si = solution.get_saturation_index()
-        assert si["Halite"] > 0
+        assert si["Halite"] > 0.01
         si_plot = solution.get_saturation_index(get_plot=True)
         assert isinstance(si_plot, dict)
+
+    def test_halite_si_under(self, monkeypatch):
+        solution = Solution({"Na+": "0.001 mol/L", "Cl-": "0.001 mol/L"}, engine="phreeqc")
+        si = solution.get_saturation_index()
+        assert si["Halite"] < -0.01
+
+    def test_saturation_index_sorted(self, monkeypatch):
+        monkeypatch.setattr(go.Figure, "show", lambda self: None)
+        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine="native")
+        si = solution.get_saturation_index(get_plot=True)
+        values = list(si.values())
+        assert values == sorted(values, reverse=True)
+
+    def test_halite_si_trend(self, monkeypatch):
+        monkeypatch.setattr(go.Figure, "show", lambda self: None)
+        solution_low = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine="native")
+        solution_high = Solution({"Na+": "10 mol/L", "Cl-": "10 mol/L"}, engine="native")
+        si_low = solution_low.get_saturation_index()
+        si_high = solution_high.get_saturation_index()
+        assert si_high["Halite"] > si_low["Halite"]
+
+    def test_halite_si_matches_phreeqc(self, monkeypatch):
+        monkeypatch.setattr(go.Figure, "show", lambda self: None)
+        composition = {"Na+": "10 mol/L", "K+": "10 mol/L","Cl-": "10 mol/L"}
+        native_solution = Solution(composition, engine="native")
+        phreeqc_solution = Solution(composition, engine="phreeqc")
+        native_si = native_solution.get_saturation_index()
+        phreeqc_si = phreeqc_solution.get_saturation_index()
+        assert pytest.approx(native_si["Halite"],rel=1e-3,abs=1e-3) == phreeqc_si["Halite"]
+
+    def test_calcite_si_matches_phreeqc(self, monkeypatch):
+        monkeypatch.setattr(go.Figure, "show", lambda self: None)
+        composition = {"Ca2+": "2 mmol/L","CO3-2": "2 mmol/L","H+": "7.0 pH"}
+        native = Solution(composition, engine="native").get_saturation_index()
+        phreeqc = Solution(composition, engine="phreeqc").get_saturation_index()
+        assert pytest.approx(native["Calcite"],rel=1e-3,abs=1e-3) == phreeqc["Calcite"]
+    
+    def test_multi_mineral_si(self, monkeypatch):
+        monkeypatch.setattr(go.Figure, "show", lambda self: None)
+        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mmol/L", "Cl-": "10 mol/L"}, engine="native")
+        si = solution.get_saturation_index()
+        assert isinstance(si, dict)
+        assert "Halite" in si
+        assert "Calcite" in si
+        assert len(si) >= 2

From 7a634840ca4cc9759d1dc7e602faf31d3908736d Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Thu, 14 May 2026 15:27:39 -0400
Subject: [PATCH 05/20] Update test_solution.py

---
 tests/test_solution.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 5352a868..0017b6aa 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1234,24 +1234,24 @@ def test_halite_si_trend(self, monkeypatch):
     def test_halite_si_matches_phreeqc(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         composition = {"Na+": "10 mol/L", "K+": "10 mol/L","Cl-": "10 mol/L"}
-        native_solution = Solution(composition, engine="native")
+        phreeqc2026_solution = Solution(composition, engine="phreeqc2026")
         phreeqc_solution = Solution(composition, engine="phreeqc")
-        native_si = native_solution.get_saturation_index()
+        phreeqc2026_si = phreeqc2026_solution.get_saturation_index()
         phreeqc_si = phreeqc_solution.get_saturation_index()
-        assert pytest.approx(native_si["Halite"],rel=1e-3,abs=1e-3) == phreeqc_si["Halite"]
+        assert pytest.approx(phreeqc2026_si["Halite"],rel=1e-3,abs=1e-3) == phreeqc_si["Halite"]
 
     def test_calcite_si_matches_phreeqc(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         composition = {"Ca2+": "2 mmol/L","CO3-2": "2 mmol/L","H+": "7.0 pH"}
-        native = Solution(composition, engine="native").get_saturation_index()
+        phreeqc2026 = Solution(composition, engine="phreeqc2026").get_saturation_index()
         phreeqc = Solution(composition, engine="phreeqc").get_saturation_index()
-        assert pytest.approx(native["Calcite"],rel=1e-3,abs=1e-3) == phreeqc["Calcite"]
+        assert pytest.approx(phreeqc2026["Calcite"],rel=1e-3,abs=1e-3) == phreeqc["Calcite"]
     
     def test_multi_mineral_si(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mmol/L", "Cl-": "10 mol/L"}, engine="native")
+        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mol/L", "Cl-": "10 mol/L"}, pH=11, engine="native")
+        solution.equilibrate(gases={"CO2": -0.5})
         si = solution.get_saturation_index()
         assert isinstance(si, dict)
-        assert "Halite" in si
-        assert "Calcite" in si
-        assert len(si) >= 2
+        assert si["Halite"] > 2.5
+        assert si["Calcite"] > 0.03

From 65a1384355f4de3dc2e1d6c0ada7d44fda365992 Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Thu, 14 May 2026 16:00:42 -0400
Subject: [PATCH 06/20] Update test_solution.py

---
 tests/test_solution.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 0017b6aa..852d1913 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1242,10 +1242,10 @@ def test_halite_si_matches_phreeqc(self, monkeypatch):
 
     def test_calcite_si_matches_phreeqc(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        composition = {"Ca2+": "2 mmol/L","CO3-2": "2 mmol/L","H+": "7.0 pH"}
-        phreeqc2026 = Solution(composition, engine="phreeqc2026").get_saturation_index()
-        phreeqc = Solution(composition, engine="phreeqc").get_saturation_index()
-        assert pytest.approx(phreeqc2026["Calcite"],rel=1e-3,abs=1e-3) == phreeqc["Calcite"]
+        composition = {"Ca2+": "2 mmol/L","CO3-2": "2 mmol/L","H+": "10**(-10.3) mol/L"}
+        phreeqc2026_si = Solution(composition, engine="phreeqc2026").get_saturation_index()
+        phreeqc_si = Solution(composition, engine="phreeqc").get_saturation_index()
+        assert pytest.approx(phreeqc2026_si["Calcite"],rel=1e-3,abs=1e-3) == phreeqc_si["Calcite"]
     
     def test_multi_mineral_si(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)

From 5871b373b1ac87f3265f71a90c2bfb85e380e2da Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Thu, 14 May 2026 16:02:41 -0400
Subject: [PATCH 07/20] Update test_solution.py

---
 tests/test_solution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 852d1913..7f484565 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1253,5 +1253,5 @@ def test_multi_mineral_si(self, monkeypatch):
         solution.equilibrate(gases={"CO2": -0.5})
         si = solution.get_saturation_index()
         assert isinstance(si, dict)
-        assert si["Halite"] > 2.5
+        assert si["Halite"] > 0.39
         assert si["Calcite"] > 0.03

From 30964be2d5560f412e2b6230de89c2412b9f54e8 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Fri, 15 May 2026 13:27:13 -0400
Subject: [PATCH 08/20] added corrections

---
 src/pyEQL/solution.py | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py
index 5b329aa2..9663f385 100644
--- a/src/pyEQL/solution.py
+++ b/src/pyEQL/solution.py
@@ -17,8 +17,6 @@
 from typing import Any, Literal
 
 import numpy as np
-import pandas as pd
-import plotly.express as px
 from maggma.stores import JSONStore, Store
 from monty.dev import deprecated
 from monty.json import MontyDecoder, MSONable
@@ -2437,12 +2435,14 @@ def _adjust_charge_balance(self, atol=1e-8) -> None:
                 return
 
     def _check_water_stability(self, tol=1e-6) -> None:
-        """Helper method to adjust the thermodynamic stability of the Solution."""     
+        """Helper method to adjust the thermodynamic stability of the Solution."""
         temp = self.temperature.to("K")
         E0_O2 = 1.229 * ureg.V
 
         lower_limit = -float(self.pH)
-        upper_limit = (ureg.faraday_constant * E0_O2 / (2.303 * ureg.R * temp)).to_base_units().magnitude - float(self.pH)
+        upper_limit = (ureg.faraday_constant * E0_O2 / (2.303 * ureg.R * temp)).to_base_units().magnitude - float(
+            self.pH
+        )
 
         if self.pE < lower_limit - tol:
             msg = (
@@ -2847,9 +2847,15 @@ def add_solvent(self, formula: str, amount: str):  # pragma: no cover
         target_mol = quantity.to("moles", "chem", mw=mw, volume=self.volume, solvent_mass=self.solvent_mass)
         self.components[formula] = target_mol.to("moles").magnitude
 
-    def get_saturation_index(self, phases=None, get_plot=None) -> float:
+    def get_saturation_index(self, get_plot=None) -> dict:
         """
-        Calculate the saturation index of a solute in the solution.
+        Calculate the saturation index of a solute in the solution based on the current engine.
+        Args:
+            get_plot (bool, optional):
+                If True, displays an interactive bar plot of saturation indices sorted from most oversaturated to least. Defaults to None (no plot).
+        Returns:
+            dict:
+                A dictionary with mineral phase names as keys and their saturation index values as values, sorted in descending order (most oversaturated to least oversaturated).
         """
 
         engine = self.engine
@@ -2870,6 +2876,9 @@ def get_saturation_index(self, phases=None, get_plot=None) -> float:
         sorted_eq_species_dict = dict(sorted(eq_species_dict.items(), key=lambda item: item[1], reverse=True))
 
         if get_plot:
+            import pandas as pd  # noqa: PLC0415
+            import plotly.express as px  # noqa: PLC0415
+
             df = pd.DataFrame(
                 {"species": list(sorted_eq_species_dict.keys()), "si": list(sorted_eq_species_dict.values())}
             )

From 114a569942203243e675dbab76a927a3b40e2c22 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Fri, 15 May 2026 13:27:43 -0400
Subject: [PATCH 09/20] added corrections

---
 ...tutorial_charge_balancing-checkpoint.ipynb | 292 ++++++++++++++++++
 1 file changed, 292 insertions(+)
 create mode 100644 docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb

diff --git a/docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb b/docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb
new file mode 100644
index 00000000..eda00cc6
--- /dev/null
+++ b/docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb
@@ -0,0 +1,292 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "76c3ccc7",
+   "metadata": {},
+   "source": [
+    "# `pyEQL` Tutorial: Charge balancing\n",
+    "\n",
+    "![pyEQL Logo](../../pyeql-logo.png)\n",
+    "\n",
+    "`pyEQL` is an open-source `python` library for solution chemistry calculations and ion properties developed by the [Kingsbury Lab](https://www.kingsburylab.org/) at Princeton University.\n",
+    "\n",
+    "[Documentation](https://pyeql.readthedocs.io/en/latest/) | [How to Install](https://pyeql.readthedocs.io/en/latest/installation.html) | [GitHub](https://github.com/rkingsbury/pyEQL) "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3d9ba4d0",
+   "metadata": {},
+   "source": [
+    "This tutorial introduces strategies for handling charge balance in a `Solution` object."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "dd840142",
+   "metadata": {},
+   "source": [
+    "### Charge balancing strategies\n",
+    "- Approach 0: `None (default)` in which case no charge balancing will be performed\n",
+    "- Approach 1: `auto` which will use the majority cation or anion (i.e., that with the largest concentration) as needed\n",
+    "- Approach 2: `pH` which will adjust the solution pH to balance charge\n",
+    "- Approach 3: `target_ion` which will use a specified cation or anion to balance charge\n",
+    "<!-- - `pE` (not currently implemented) which will adjust the redox equilibrium to balance charge, or the name of a dissolved species e.g. 'Ca+2' or 'Cl-' that will be added/subtracted to balance charge. -->\n",
+    "\n",
+    "### Important note between `balance_charge` and `charge_balance`:\n",
+    "- The `balance_charge` *kwarg* controls how charge balancing is performed.\n",
+    "- The `charge_balance` *property* returns the net electric charge of the solution."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5509e64a",
+   "metadata": {},
+   "source": [
+    "### Import `Solution` object"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "ae238c71",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/st2591/.conda/envs/st2591_tiger3_dev/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
+      "  from .autonotebook import tqdm as notebook_tqdm\n"
+     ]
+    }
+   ],
+   "source": [
+    "from pyEQL import Solution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "91c165bb",
+   "metadata": {},
+   "source": [
+    "### Example of a non-charge balanced solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "4702be4f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ion_dict = {\"Na+\": \"0.1 mol/L\", \"K+\": \"0.2 mol/L\", \"SO4-2\": \"0.50 mol/L\"}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ee5d0ae2",
+   "metadata": {},
+   "source": [
+    "### **Approach 0:** `balance_charge` = `None (default)`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "abed484d",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Solution charge balance: -7.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "sol = Solution(solutes=ion_dict, pH=7)\n",
+    "print(f\"Solution charge balance: {sol.charge_balance}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "387ae5c1",
+   "metadata": {},
+   "source": [
+    "### **Approach 1:** `balance_charge` = \"auto\"\n",
+    "\n",
+    "The `balance_charge = \"auto\"` will use the majority cation or anion (i.e., that with the largest concentration), in this case `K[+1]`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "e0b10d4e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Balance charge: auto, Solution charge balance: -5.551115288561905e-17\n"
+     ]
+    }
+   ],
+   "source": [
+    "sol1 = Solution(solutes=ion_dict, pH=7, balance_charge=\"auto\")\n",
+    "print(f\"Balance charge: {sol1.balance_charge}, Solution charge balance: {sol1.charge_balance}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "32164496",
+   "metadata": {},
+   "source": [
+    "### **Approach 2:** `balance_charge` = \"pH\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "d787c7bd",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Adjusted pH: -0.8450980400142569, Solution charge balance: 3.4778541082882235e-16\n"
+     ]
+    }
+   ],
+   "source": [
+    "sol2 = Solution(solutes=ion_dict, pH=7, balance_charge=\"pH\")\n",
+    "print(f\"Adjusted pH: {sol2.pH}, Solution charge balance: {sol2.charge_balance}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d078528c",
+   "metadata": {},
+   "source": [
+    "### **Approach 3:** `balance_charge` = \"target ion\"\n",
+    "Charge balance can also be done by designating a particular \"target ion\", in this case `Mg[+2]`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "d633c3c9",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Charge balancing species Mg[+2] was not found in the solution!. Species {'OH[-1]', 'SO4[-2]', 'Na[+1]', 'K[+1]', 'H[+1]'} were found.\n"
+     ]
+    }
+   ],
+   "source": [
+    "try:\n",
+    "    sol3 = Solution(solutes=ion_dict, pH=7, balance_charge=\"Mg+2\")\n",
+    "    print(f\"Solution charge balance: {sol3.charge_balance}\")\n",
+    "except Exception as e:\n",
+    "    print(e)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e29fd23b",
+   "metadata": {},
+   "source": [
+    "Notice: To use `balance_charge`, the target ion must be present in the `Solution` object. \n",
+    "We add a trace concentration (0.5 mol/L)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "ead247ce",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Balance charge: Mg[+2], Solution charge balance: 5.421010862427522e-20\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/st2591/pyEQL_shaun/pyEQL/src/pyEQL/solution.py:498: RuntimeWarning: invalid value encountered in log10\n",
+      "  return float(-1 * np.log10(amt))\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Add target ion in Solution object with 0.5 mol/L\n",
+    "ion_dict[\"Mg+2\"] = \"0.5 mol/L\"\n",
+    "try:\n",
+    "    sol3 = Solution(solutes=ion_dict, pH=7, balance_charge=\"Mg+2\")\n",
+    "    print(f\"Balance charge: {sol3.balance_charge}, Solution charge balance: {sol3.charge_balance}\")\n",
+    "except Exception as e:\n",
+    "    print(e)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d9ee787f-ce70-4a04-9f9b-41c8345332e0",
+   "metadata": {},
+   "source": [
+    "Notice: Equilibration or `equilibrate()` may adjust the solution’s charge balance, regardless of the charge‑balancing method selected."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "784b1cf7-cc92-48f7-a02e-a2767b82bd92",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol3 = Solution(solutes=ion_dict, pH=7, balance_charge=\"Mg+2\")\n",
+    "sol3.equilibrate()\n",
+    "print(f\"Balance charge: {sol3.balance_charge}, Solution charge balance: {sol3.charge_balance}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ad6b1b29-3337-4cc7-87fb-c91631851dc1",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "st2591_tiger3_dev [~/.conda/envs/st2591_tiger3_dev/]",
+   "language": "python",
+   "name": "conda_st2591_tiger3_dev"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From 9baf96c6e9c33075855194b97bd3ad767fbed53d Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Fri, 15 May 2026 13:39:57 -0400
Subject: [PATCH 10/20] added engine parameterize test

---
 tests/test_solution.py | 34 ++++++++++++++++++----------------
 1 file changed, 18 insertions(+), 16 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 7f484565..7d32fb4e 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1198,11 +1198,11 @@ def test_should_use_major_salt_molar_volume_to_calculate_solute_volume_when_para
                 "OH-", "size.molar_volume"
             )
         assert solute_volume_without_protons_and_hydroxide.m == expected_solute_volume.m
-        
 
+
+@pytest.mark.parametrize("engine", ["native", "phreeqc", "phreeqc2026"])
 class TestSaturationIndex:
     @staticmethod
-    @pytest.mark.parametrize("engine", ["native", "phreeqc", "phreeqc2026"])
     def test_halite_si(engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
@@ -1211,45 +1211,47 @@ def test_halite_si(engine, monkeypatch):
         si_plot = solution.get_saturation_index(get_plot=True)
         assert isinstance(si_plot, dict)
 
-    def test_halite_si_under(self, monkeypatch):
-        solution = Solution({"Na+": "0.001 mol/L", "Cl-": "0.001 mol/L"}, engine="phreeqc")
+    def test_halite_si_under(self, engine, monkeypatch):
+        solution = Solution({"Na+": "0.001 mol/L", "Cl-": "0.001 mol/L"}, engine=engine)
         si = solution.get_saturation_index()
         assert si["Halite"] < -0.01
 
-    def test_saturation_index_sorted(self, monkeypatch):
+    def test_saturation_index_sorted(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine="native")
+        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
         si = solution.get_saturation_index(get_plot=True)
         values = list(si.values())
         assert values == sorted(values, reverse=True)
 
-    def test_halite_si_trend(self, monkeypatch):
+    def test_halite_si_trend(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution_low = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine="native")
-        solution_high = Solution({"Na+": "10 mol/L", "Cl-": "10 mol/L"}, engine="native")
+        solution_low = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine=engine)
+        solution_high = Solution({"Na+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
         si_low = solution_low.get_saturation_index()
         si_high = solution_high.get_saturation_index()
         assert si_high["Halite"] > si_low["Halite"]
 
     def test_halite_si_matches_phreeqc(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        composition = {"Na+": "10 mol/L", "K+": "10 mol/L","Cl-": "10 mol/L"}
+        composition = {"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}
         phreeqc2026_solution = Solution(composition, engine="phreeqc2026")
         phreeqc_solution = Solution(composition, engine="phreeqc")
         phreeqc2026_si = phreeqc2026_solution.get_saturation_index()
         phreeqc_si = phreeqc_solution.get_saturation_index()
-        assert pytest.approx(phreeqc2026_si["Halite"],rel=1e-3,abs=1e-3) == phreeqc_si["Halite"]
+        assert pytest.approx(phreeqc2026_si["Halite"], rel=1e-3, abs=1e-3) == phreeqc_si["Halite"]
 
     def test_calcite_si_matches_phreeqc(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        composition = {"Ca2+": "2 mmol/L","CO3-2": "2 mmol/L","H+": "10**(-10.3) mol/L"}
+        composition = {"Ca2+": "2 mmol/L", "CO3-2": "2 mmol/L", "H+": "10**(-10.3) mol/L"}
         phreeqc2026_si = Solution(composition, engine="phreeqc2026").get_saturation_index()
         phreeqc_si = Solution(composition, engine="phreeqc").get_saturation_index()
-        assert pytest.approx(phreeqc2026_si["Calcite"],rel=1e-3,abs=1e-3) == phreeqc_si["Calcite"]
-    
-    def test_multi_mineral_si(self, monkeypatch):
+        assert pytest.approx(phreeqc2026_si["Calcite"], rel=1e-3, abs=1e-3) == phreeqc_si["Calcite"]
+
+    def test_multi_mineral_si(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mol/L", "Cl-": "10 mol/L"}, pH=11, engine="native")
+        solution = Solution(
+            {"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mol/L", "Cl-": "10 mol/L"}, pH=11, engine=engine
+        )
         solution.equilibrate(gases={"CO2": -0.5})
         si = solution.get_saturation_index()
         assert isinstance(si, dict)

From 303e0d833b7a304eb2edce0a3cf3ad4f20a4d130 Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Sat, 16 May 2026 08:54:50 -0400
Subject: [PATCH 11/20] Update test_solution.py

---
 tests/test_solution.py | 39 ++++++++++++++++++++-------------------
 1 file changed, 20 insertions(+), 19 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 7d32fb4e..5f2b3172 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1203,7 +1203,7 @@ def test_should_use_major_salt_molar_volume_to_calculate_solute_volume_when_para
 @pytest.mark.parametrize("engine", ["native", "phreeqc", "phreeqc2026"])
 class TestSaturationIndex:
     @staticmethod
-    def test_halite_si(engine, monkeypatch):
+    def test_halite_si_over(engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
         si = solution.get_saturation_index()
@@ -1216,6 +1216,11 @@ def test_halite_si_under(self, engine, monkeypatch):
         si = solution.get_saturation_index()
         assert si["Halite"] < -0.01
 
+    def test_halite_si_near(self, engine, monkeypatch):
+        solution = Solution({"Na+": "6 mol/L", "Cl-": "6 mol/L"}, engine=engine)
+        si = solution.get_saturation_index()
+        assert -0.1 < si["Halite"] < 0.1
+
     def test_saturation_index_sorted(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
@@ -1225,21 +1230,12 @@ def test_saturation_index_sorted(self, engine, monkeypatch):
 
     def test_halite_si_trend(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution_low = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine=engine)
-        solution_high = Solution({"Na+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
+        solution_low = Solution({"Na+": "0.4 mol/L", "Cl-": "0.4 mol/L"}, engine=engine)
+        solution_high = Solution({"Na+": "0.75 mol/L", "Cl-": "0.75 mol/L"}, engine=engine)
         si_low = solution_low.get_saturation_index()
         si_high = solution_high.get_saturation_index()
         assert si_high["Halite"] > si_low["Halite"]
 
-    def test_halite_si_matches_phreeqc(self, monkeypatch):
-        monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        composition = {"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}
-        phreeqc2026_solution = Solution(composition, engine="phreeqc2026")
-        phreeqc_solution = Solution(composition, engine="phreeqc")
-        phreeqc2026_si = phreeqc2026_solution.get_saturation_index()
-        phreeqc_si = phreeqc_solution.get_saturation_index()
-        assert pytest.approx(phreeqc2026_si["Halite"], rel=1e-3, abs=1e-3) == phreeqc_si["Halite"]
-
     def test_calcite_si_matches_phreeqc(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         composition = {"Ca2+": "2 mmol/L", "CO3-2": "2 mmol/L", "H+": "10**(-10.3) mol/L"}
@@ -1247,13 +1243,18 @@ def test_calcite_si_matches_phreeqc(self, monkeypatch):
         phreeqc_si = Solution(composition, engine="phreeqc").get_saturation_index()
         assert pytest.approx(phreeqc2026_si["Calcite"], rel=1e-3, abs=1e-3) == phreeqc_si["Calcite"]
 
-    def test_multi_mineral_si(self, engine, monkeypatch):
+    @pytest.mark.parametrize("engine", ["native"])
+    def test_saturation_index_ideal_not_supported(self, monkeypatch):
+        solution = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine="ideal",)
+        with pytest.raises(NotImplementedError):
+            solution.get_saturation_index()
+
+    @pytest.mark.skip(reason="temporarily disabled")
+    def test_multi_mineral_si(self, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution = Solution(
-            {"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mol/L", "Cl-": "10 mol/L"}, pH=11, engine=engine
-        )
-        solution.equilibrate(gases={"CO2": -0.5})
+        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mmol/L", "Cl-": "10 mol/L"}, engine="native")
         si = solution.get_saturation_index()
         assert isinstance(si, dict)
-        assert si["Halite"] > 0.39
-        assert si["Calcite"] > 0.03
+        assert "Halite" in si
+        assert "Calcite" in si
+        assert len(si) >= 2

From 8abffd45a99984f09ffcba0c6d9fcca9289db08c Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Sat, 16 May 2026 09:13:45 -0400
Subject: [PATCH 12/20] Update test_solution.py

---
 tests/test_solution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 5f2b3172..c203c3b0 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1236,7 +1236,7 @@ def test_halite_si_trend(self, engine, monkeypatch):
         si_high = solution_high.get_saturation_index()
         assert si_high["Halite"] > si_low["Halite"]
 
-    def test_calcite_si_matches_phreeqc(self, monkeypatch):
+    def test_calcite_si_matches_phreeqc(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         composition = {"Ca2+": "2 mmol/L", "CO3-2": "2 mmol/L", "H+": "10**(-10.3) mol/L"}
         phreeqc2026_si = Solution(composition, engine="phreeqc2026").get_saturation_index()

From 693e730ec2cd049e2ef25b21f56d55598af4e535 Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Sat, 16 May 2026 09:32:15 -0400
Subject: [PATCH 13/20] Update test_solution.py

---
 tests/test_solution.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index c203c3b0..3e2dabb8 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1243,9 +1243,8 @@ def test_calcite_si_matches_phreeqc(self, engine, monkeypatch):
         phreeqc_si = Solution(composition, engine="phreeqc").get_saturation_index()
         assert pytest.approx(phreeqc2026_si["Calcite"], rel=1e-3, abs=1e-3) == phreeqc_si["Calcite"]
 
-    @pytest.mark.parametrize("engine", ["native"])
-    def test_saturation_index_ideal_not_supported(self, monkeypatch):
-        solution = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine="ideal",)
+    def test_saturation_index_ideal_not_supported(self, engine, monkeypatch):
+        solution = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine="ideal")
         with pytest.raises(NotImplementedError):
             solution.get_saturation_index()
 

From 562e3077dd470ca0347a1e49d8178eaf0aa55011 Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Sat, 16 May 2026 09:42:56 -0400
Subject: [PATCH 14/20] Update test_solution.py

---
 tests/test_solution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 3e2dabb8..7c6ac509 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1249,7 +1249,7 @@ def test_saturation_index_ideal_not_supported(self, engine, monkeypatch):
             solution.get_saturation_index()
 
     @pytest.mark.skip(reason="temporarily disabled")
-    def test_multi_mineral_si(self, monkeypatch):
+    def test_multi_mineral_si(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mmol/L", "Cl-": "10 mol/L"}, engine="native")
         si = solution.get_saturation_index()

From ecc71cdfb196897765c5827bcc46deb8dc2525f3 Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Sat, 16 May 2026 09:55:23 -0400
Subject: [PATCH 15/20] Update test_solution.py

---
 tests/test_solution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index 7c6ac509..fbe4ec2a 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1219,7 +1219,7 @@ def test_halite_si_under(self, engine, monkeypatch):
     def test_halite_si_near(self, engine, monkeypatch):
         solution = Solution({"Na+": "6 mol/L", "Cl-": "6 mol/L"}, engine=engine)
         si = solution.get_saturation_index()
-        assert -0.1 < si["Halite"] < 0.1
+        assert -1 < si["Halite"] < 1
 
     def test_saturation_index_sorted(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)

From 3d4364f4426848f017f3242cdba9e6b1bf71dac4 Mon Sep 17 00:00:00 2001
From: YitongPan1 <119685432+YitongPan1@users.noreply.github.com>
Date: Mon, 18 May 2026 06:14:58 -0500
Subject: [PATCH 16/20] Update test_solution.py

---
 tests/test_solution.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index fbe4ec2a..56d5b4b3 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1249,7 +1249,7 @@ def test_saturation_index_ideal_not_supported(self, engine, monkeypatch):
             solution.get_saturation_index()
 
     @pytest.mark.skip(reason="temporarily disabled")
-    def test_multi_mineral_si(self, engine, monkeypatch):
+    def test_multi_equilibrate_si(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mmol/L", "Cl-": "10 mol/L"}, engine="native")
         si = solution.get_saturation_index()
@@ -1257,3 +1257,6 @@ def test_multi_mineral_si(self, engine, monkeypatch):
         assert "Halite" in si
         assert "Calcite" in si
         assert len(si) >= 2
+        solution.equilibrate()
+        assert -1 < si["Halite"] < 1
+        assert -1 < si["Calcite"] < 1

From 8bbf078ee978652def410d56b99c8d7f6d104382 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Thu, 28 May 2026 22:10:57 -0400
Subject: [PATCH 17/20] solution.py correction

---
 src/pyEQL/solution.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py
index 9663f385..8788f5b9 100644
--- a/src/pyEQL/solution.py
+++ b/src/pyEQL/solution.py
@@ -2848,8 +2848,18 @@ def add_solvent(self, formula: str, amount: str):  # pragma: no cover
         self.components[formula] = target_mol.to("moles").magnitude
 
     def get_saturation_index(self, get_plot=None) -> dict:
-        """
-        Calculate the saturation index of a solute in the solution based on the current engine.
+        r"""
+        Calculate the saturation index of a solute in the solution.
+        Notes:
+            The saturation index (:math:`\mathrm{SI}`) is defined as log10(IAP/Ksp), where IAP is the ion activity product and Ksp is the solubility product constant.
+            This method calculates the saturation index based on the active engine and database from `__init__`. The interpretation of the saturation index values is as follows:
+
+            - :math:`\mathrm{SI} < 0`: The solution is **undersaturated**. The solid tends to dissolve if present.
+
+            - :math:`\mathrm{SI} = 0`: The solution is **at saturation equilibrium**. Therefore, at the saturation limit, the SI is zero.
+
+            - :math:`\mathrm{SI} > 0`: The solution is **supersaturated**. Precipitation is thermodynamically favored, although kinetic factors may delay or prevent it.
+
         Args:
             get_plot (bool, optional):
                 If True, displays an interactive bar plot of saturation indices sorted from most oversaturated to least. Defaults to None (no plot).

From b0ade6dcf911bcca45b7389f64535361b158e347 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Fri, 29 May 2026 12:41:30 -0400
Subject: [PATCH 18/20] small corrections

---
 tests/test_solution.py | 39 +++++++++++++++++----------------------
 1 file changed, 17 insertions(+), 22 deletions(-)

diff --git a/tests/test_solution.py b/tests/test_solution.py
index c8602a3d..07be19e1 100644
--- a/tests/test_solution.py
+++ b/tests/test_solution.py
@@ -1252,6 +1252,8 @@ def test_halite_si_over(engine, monkeypatch):
         assert si["Halite"] > 0.01
         si_plot = solution.get_saturation_index(get_plot=True)
         assert isinstance(si_plot, dict)
+        values = list(si.values())
+        assert values == sorted(values, reverse=True)
 
     def test_halite_si_under(self, engine, monkeypatch):
         solution = Solution({"Na+": "0.001 mol/L", "Cl-": "0.001 mol/L"}, engine=engine)
@@ -1261,40 +1263,33 @@ def test_halite_si_under(self, engine, monkeypatch):
     def test_halite_si_near(self, engine, monkeypatch):
         solution = Solution({"Na+": "6 mol/L", "Cl-": "6 mol/L"}, engine=engine)
         si = solution.get_saturation_index()
-        assert -1 < si["Halite"] < 1
-
-    def test_saturation_index_sorted(self, engine, monkeypatch):
-        monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution = Solution({"Na+": "10 mol/L", "K+": "10 mol/L", "Cl-": "10 mol/L"}, engine=engine)
-        si = solution.get_saturation_index(get_plot=True)
-        values = list(si.values())
-        assert values == sorted(values, reverse=True)
-
-    def test_halite_si_trend(self, engine, monkeypatch):
-        monkeypatch.setattr(go.Figure, "show", lambda self: None)
-        solution_low = Solution({"Na+": "0.4 mol/L", "Cl-": "0.4 mol/L"}, engine=engine)
-        solution_high = Solution({"Na+": "0.75 mol/L", "Cl-": "0.75 mol/L"}, engine=engine)
-        si_low = solution_low.get_saturation_index()
-        si_high = solution_high.get_saturation_index()
-        assert si_high["Halite"] > si_low["Halite"]
+        assert -0.5 < si["Halite"] < 0.0
 
     def test_calcite_si_matches_phreeqc(self, engine, monkeypatch):
+        from pyEQL.engines import Phreeqc2026EOS, PhreeqcEOS  # noqa: PLC0415
+
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
+        phreeqc_eos = PhreeqcEOS(phreeqc_db="phreeqc.dat")
+        phreeqc2026_eos = Phreeqc2026EOS(phreeqc_db="phreeqc.dat")
         composition = {"Ca2+": "2 mmol/L", "CO3-2": "2 mmol/L", "H+": "10**(-10.3) mol/L"}
-        phreeqc2026_si = Solution(composition, engine="phreeqc2026").get_saturation_index()
-        phreeqc_si = Solution(composition, engine="phreeqc").get_saturation_index()
+
+        phreeqc_si = Solution(composition, engine=phreeqc_eos).get_saturation_index()
+        phreeqc2026_si = Solution(composition, engine=phreeqc2026_eos).get_saturation_index()
         assert pytest.approx(phreeqc2026_si["Calcite"], rel=1e-3, abs=1e-3) == phreeqc_si["Calcite"]
+        assert phreeqc_si["Calcite"] == pytest.approx(2.14, rel=1e-2, abs=1e-2)
+        assert phreeqc2026_si["Calcite"] == pytest.approx(2.14, rel=1e-2, abs=1e-2)
 
     def test_saturation_index_ideal_not_supported(self, engine, monkeypatch):
         solution = Solution({"Na+": "1 mol/L", "Cl-": "1 mol/L"}, engine="ideal")
         with pytest.raises(NotImplementedError):
             solution.get_saturation_index()
 
-    @pytest.mark.skip(reason="temporarily disabled")
+    # @pytest.mark.skip(reason="temporarily disabled")
     def test_multi_equilibrate_si(self, engine, monkeypatch):
         monkeypatch.setattr(go.Figure, "show", lambda self: None)
         solution = Solution(
-            {"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "2 mmol/L", "Cl-": "10 mol/L"}, engine="native"
+            {"Na+": "10 mol/L", "K+": "10 mol/L", "Ca2+": "0.05 mol/L", "Cl-": "10 mol/L", "CO3-2": "0.05 mol/L"},
+            engine="native",
         )
         si = solution.get_saturation_index()
         assert isinstance(si, dict)
@@ -1302,5 +1297,5 @@ def test_multi_equilibrate_si(self, engine, monkeypatch):
         assert "Calcite" in si
         assert len(si) >= 2
         solution.equilibrate()
-        assert -1 < si["Halite"] < 1
-        assert -1 < si["Calcite"] < 1
+        assert 0 < si["Halite"] < 1
+        assert 0 < si["Calcite"] < 1

From 2f5e5fa93eb934b0b2130e894c3f90c5af9bec24 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Fri, 29 May 2026 12:43:28 -0400
Subject: [PATCH 19/20] removed checkpoints

---
 ...tutorial_charge_balancing-checkpoint.ipynb | 292 ------------------
 1 file changed, 292 deletions(-)
 delete mode 100644 docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb

diff --git a/docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb b/docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb
deleted file mode 100644
index eda00cc6..00000000
--- a/docs/examples/.ipynb_checkpoints/pyeql_tutorial_charge_balancing-checkpoint.ipynb
+++ /dev/null
@@ -1,292 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "id": "76c3ccc7",
-   "metadata": {},
-   "source": [
-    "# `pyEQL` Tutorial: Charge balancing\n",
-    "\n",
-    "![pyEQL Logo](../../pyeql-logo.png)\n",
-    "\n",
-    "`pyEQL` is an open-source `python` library for solution chemistry calculations and ion properties developed by the [Kingsbury Lab](https://www.kingsburylab.org/) at Princeton University.\n",
-    "\n",
-    "[Documentation](https://pyeql.readthedocs.io/en/latest/) | [How to Install](https://pyeql.readthedocs.io/en/latest/installation.html) | [GitHub](https://github.com/rkingsbury/pyEQL) "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "3d9ba4d0",
-   "metadata": {},
-   "source": [
-    "This tutorial introduces strategies for handling charge balance in a `Solution` object."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "dd840142",
-   "metadata": {},
-   "source": [
-    "### Charge balancing strategies\n",
-    "- Approach 0: `None (default)` in which case no charge balancing will be performed\n",
-    "- Approach 1: `auto` which will use the majority cation or anion (i.e., that with the largest concentration) as needed\n",
-    "- Approach 2: `pH` which will adjust the solution pH to balance charge\n",
-    "- Approach 3: `target_ion` which will use a specified cation or anion to balance charge\n",
-    "<!-- - `pE` (not currently implemented) which will adjust the redox equilibrium to balance charge, or the name of a dissolved species e.g. 'Ca+2' or 'Cl-' that will be added/subtracted to balance charge. -->\n",
-    "\n",
-    "### Important note between `balance_charge` and `charge_balance`:\n",
-    "- The `balance_charge` *kwarg* controls how charge balancing is performed.\n",
-    "- The `charge_balance` *property* returns the net electric charge of the solution."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "5509e64a",
-   "metadata": {},
-   "source": [
-    "### Import `Solution` object"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "id": "ae238c71",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/home/st2591/.conda/envs/st2591_tiger3_dev/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
-      "  from .autonotebook import tqdm as notebook_tqdm\n"
-     ]
-    }
-   ],
-   "source": [
-    "from pyEQL import Solution"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "91c165bb",
-   "metadata": {},
-   "source": [
-    "### Example of a non-charge balanced solution."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "id": "4702be4f",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "ion_dict = {\"Na+\": \"0.1 mol/L\", \"K+\": \"0.2 mol/L\", \"SO4-2\": \"0.50 mol/L\"}"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "ee5d0ae2",
-   "metadata": {},
-   "source": [
-    "### **Approach 0:** `balance_charge` = `None (default)`"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "id": "abed484d",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Solution charge balance: -7.0\n"
-     ]
-    }
-   ],
-   "source": [
-    "sol = Solution(solutes=ion_dict, pH=7)\n",
-    "print(f\"Solution charge balance: {sol.charge_balance}\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "387ae5c1",
-   "metadata": {},
-   "source": [
-    "### **Approach 1:** `balance_charge` = \"auto\"\n",
-    "\n",
-    "The `balance_charge = \"auto\"` will use the majority cation or anion (i.e., that with the largest concentration), in this case `K[+1]`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "id": "e0b10d4e",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Balance charge: auto, Solution charge balance: -5.551115288561905e-17\n"
-     ]
-    }
-   ],
-   "source": [
-    "sol1 = Solution(solutes=ion_dict, pH=7, balance_charge=\"auto\")\n",
-    "print(f\"Balance charge: {sol1.balance_charge}, Solution charge balance: {sol1.charge_balance}\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "32164496",
-   "metadata": {},
-   "source": [
-    "### **Approach 2:** `balance_charge` = \"pH\""
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "id": "d787c7bd",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Adjusted pH: -0.8450980400142569, Solution charge balance: 3.4778541082882235e-16\n"
-     ]
-    }
-   ],
-   "source": [
-    "sol2 = Solution(solutes=ion_dict, pH=7, balance_charge=\"pH\")\n",
-    "print(f\"Adjusted pH: {sol2.pH}, Solution charge balance: {sol2.charge_balance}\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "d078528c",
-   "metadata": {},
-   "source": [
-    "### **Approach 3:** `balance_charge` = \"target ion\"\n",
-    "Charge balance can also be done by designating a particular \"target ion\", in this case `Mg[+2]`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 6,
-   "id": "d633c3c9",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Charge balancing species Mg[+2] was not found in the solution!. Species {'OH[-1]', 'SO4[-2]', 'Na[+1]', 'K[+1]', 'H[+1]'} were found.\n"
-     ]
-    }
-   ],
-   "source": [
-    "try:\n",
-    "    sol3 = Solution(solutes=ion_dict, pH=7, balance_charge=\"Mg+2\")\n",
-    "    print(f\"Solution charge balance: {sol3.charge_balance}\")\n",
-    "except Exception as e:\n",
-    "    print(e)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "e29fd23b",
-   "metadata": {},
-   "source": [
-    "Notice: To use `balance_charge`, the target ion must be present in the `Solution` object. \n",
-    "We add a trace concentration (0.5 mol/L)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 7,
-   "id": "ead247ce",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Balance charge: Mg[+2], Solution charge balance: 5.421010862427522e-20\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/home/st2591/pyEQL_shaun/pyEQL/src/pyEQL/solution.py:498: RuntimeWarning: invalid value encountered in log10\n",
-      "  return float(-1 * np.log10(amt))\n"
-     ]
-    }
-   ],
-   "source": [
-    "# Add target ion in Solution object with 0.5 mol/L\n",
-    "ion_dict[\"Mg+2\"] = \"0.5 mol/L\"\n",
-    "try:\n",
-    "    sol3 = Solution(solutes=ion_dict, pH=7, balance_charge=\"Mg+2\")\n",
-    "    print(f\"Balance charge: {sol3.balance_charge}, Solution charge balance: {sol3.charge_balance}\")\n",
-    "except Exception as e:\n",
-    "    print(e)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "d9ee787f-ce70-4a04-9f9b-41c8345332e0",
-   "metadata": {},
-   "source": [
-    "Notice: Equilibration or `equilibrate()` may adjust the solution’s charge balance, regardless of the charge‑balancing method selected."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "784b1cf7-cc92-48f7-a02e-a2767b82bd92",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "sol3 = Solution(solutes=ion_dict, pH=7, balance_charge=\"Mg+2\")\n",
-    "sol3.equilibrate()\n",
-    "print(f\"Balance charge: {sol3.balance_charge}, Solution charge balance: {sol3.charge_balance}\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "ad6b1b29-3337-4cc7-87fb-c91631851dc1",
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "st2591_tiger3_dev [~/.conda/envs/st2591_tiger3_dev/]",
-   "language": "python",
-   "name": "conda_st2591_tiger3_dev"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.12.2"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}

From 6cd4529dcae6d49613c5916bde069d14a113d2c9 Mon Sep 17 00:00:00 2001
From: Sui Xiong Tay <st2591@tiger3.princeton.edu>
Date: Fri, 29 May 2026 13:47:24 -0400
Subject: [PATCH 20/20] trigger ci