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 89b86447..3a743dc8 100644
--- a/src/pyEQL/solution.py
+++ b/src/pyEQL/solution.py
@@ -2917,3 +2917,68 @@ 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, get_plot=None) -> dict:
+        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).
+        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
+
+        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:
+            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())}
+            )
+
+            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 808edd30..07be19e1 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
@@ -1239,3 +1240,62 @@ 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
+    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()
+        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)
+        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.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"}
+
+        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")
+    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+": "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)
+        assert "Halite" in si
+        assert "Calcite" in si
+        assert len(si) >= 2
+        solution.equilibrate()
+        assert 0 < si["Halite"] < 1
+        assert 0 < si["Calcite"] < 1