diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 14d991949461768776552bd881ea2867f460db16..cde24306b386c0fd1e53530258a256231dfbff5e 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -212,12 +212,14 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
         os << "SOLUTION " << chemical_system_id + 1 << "\n";
         os << aqueous_solution << "\n";
 
-        auto const& equilibrium_phases =
-            phreeqc_io._equilibrium_phases[chemical_system_id];
+        auto const& equilibrium_phases = phreeqc_io._equilibrium_phases;
         if (!equilibrium_phases.empty())
         {
             os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
-            os << equilibrium_phases << "\n";
+            for (auto const& equilibrium_phase : equilibrium_phases)
+            {
+                equilibrium_phase.print(os, chemical_system_id);
+            }
         }
 
         auto const& kinetic_reactants = phreeqc_io._kinetic_reactants;
diff --git a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp
index 378eb2a4cdbc1779015a01c08e0f243f2f00c706..52253ab3c5c221ac11dcebefa1af93c298b8b453 100644
--- a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp
+++ b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp
@@ -13,15 +13,10 @@
 
 namespace ChemistryLib
 {
-std::ostream& operator<<(std::ostream& os,
-                         EquilibriumPhase const& equilibrium_phase)
+void EquilibriumPhase::print(std::ostream& os,
+                             std::size_t const chemical_system_id) const
 {
-    os << equilibrium_phase.name;
-
-    os << " " << equilibrium_phase.saturation_index;
-
-    os << " " << equilibrium_phase.amount << "\n";
-
-    return os;
+    os << name << " " << saturation_index << " "
+       << (*amount)[chemical_system_id] << "\n";
 }
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h
index 51e2eec934eee1af5ff0a1e8b83b9e39c342ba17..19fe180834bc34bdda5bd579439dc4914fc7c15b 100644
--- a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h
+++ b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h
@@ -34,8 +34,7 @@ struct EquilibriumPhase
     {
     }
 
-    friend std::ostream& operator<<(std::ostream& os,
-                                    EquilibriumPhase const& equilibrium_phase);
+    void print(std::ostream& os, std::size_t const chemical_system_id) const;
 
     std::string const name;
     MeshLib::PropertyVector<double>* amount;