diff --git a/ChemistryLib/CreatePhreeqcIO.cpp b/ChemistryLib/CreatePhreeqcIO.cpp
index 9100dc63a6f636b5ddd9ca0cbdbeee9b2eec389a..61b542d0f3e6dcf0ba41713b13a47cafca5905d5 100644
--- a/ChemistryLib/CreatePhreeqcIO.cpp
+++ b/ChemistryLib/CreatePhreeqcIO.cpp
@@ -7,8 +7,6 @@
  *
  */
 
-#include <fstream>
-
 #include "BaseLib/ConfigTreeUtil.h"
 #include "BaseLib/Error.h"
 #include "CreateOutput.h"
diff --git a/ChemistryLib/Output.cpp b/ChemistryLib/Output.cpp
index 64593e54754665ea5c905a0806d55a370e4eba08..a904efb1bb2981091da66ac5246f6ad72303f4ad 100644
--- a/ChemistryLib/Output.cpp
+++ b/ChemistryLib/Output.cpp
@@ -7,69 +7,69 @@
  *
  */
 
-#include <fstream>
+#include <ostream>
 
 #include "Output.h"
 
 namespace ChemistryLib
 {
-std::ofstream& operator<<(std::ofstream& out,
-                          BasicOutputSetups const& basic_output_setups)
+std::ostream& operator<<(std::ostream& os,
+                         BasicOutputSetups const& basic_output_setups)
 {
-    out << "-file " << basic_output_setups.output_file << "\n";
-    out << "-high_precision " << std::boolalpha
-        << basic_output_setups.use_high_precision << "\n";
-    out << "-simulation " << std::boolalpha
-        << basic_output_setups.display_simulation_id << "\n";
-    out << "-state " << std::boolalpha << basic_output_setups.display_state
-        << "\n";
-    out << "-distance " << std::boolalpha
-        << basic_output_setups.display_distance << "\n";
-    out << "-time " << std::boolalpha
-        << basic_output_setups.display_current_time << "\n";
-    out << "-step " << std::boolalpha << basic_output_setups.display_time_step
-        << "\n";
+    os << "-file " << basic_output_setups.output_file << "\n";
+    os << "-high_precision " << std::boolalpha
+       << basic_output_setups.use_high_precision << "\n";
+    os << "-simulation " << std::boolalpha
+       << basic_output_setups.display_simulation_id << "\n";
+    os << "-state " << std::boolalpha << basic_output_setups.display_state
+       << "\n";
+    os << "-distance " << std::boolalpha << basic_output_setups.display_distance
+       << "\n";
+    os << "-time " << std::boolalpha << basic_output_setups.display_current_time
+       << "\n";
+    os << "-step " << std::boolalpha << basic_output_setups.display_time_step
+       << "\n";
 
-    return out;
+    return os;
 }
 
-std::ofstream& operator<<(std::ofstream& out, Output const& output)
+std::ostream& operator<<(std::ostream& os, Output const& output)
 {
-    out << output.basic_output_setups;
+    os << output.basic_output_setups;
 
     auto const component_items =
         output.getOutputItemsByItemType(ItemType::Component);
-    out << "-totals";
+    os << "-totals";
     for (auto const& component_item : component_items)
     {
-        out << " " << component_item.name;
+        os << " " << component_item.name;
     }
-    out << "\n";
+    os << "\n";
 
     auto const equilibrium_phase_items =
         output.getOutputItemsByItemType(ItemType::EquilibriumPhase);
     if (!equilibrium_phase_items.empty())
     {
-        out << "-equilibrium_phases";
+        os << "-equilibrium_phases";
         for (auto const& equilibrium_phase_item : equilibrium_phase_items)
         {
-            out << " " << equilibrium_phase_item.name;
+            os << " " << equilibrium_phase_item.name;
         }
-        out << "\n";
+        os << "\n";
     }
 
     auto const kinetic_reactant_items =
         output.getOutputItemsByItemType(ItemType::KineticReactant);
     if (!kinetic_reactant_items.empty())
     {
-        out << "-kinetic_reactants";
+        os << "-kinetic_reactants";
         for (auto const& kinetic_reactant_item : kinetic_reactant_items)
         {
-            out << " " << kinetic_reactant_item.name;
+            os << " " << kinetic_reactant_item.name;
         }
-        out << "\n";
+        os << "\n";
     }
 
-    return out;
+    return os;
 }
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/Output.h b/ChemistryLib/Output.h
index cc5b1f2ea576b66e4114b7147005df78374ca2d8..1cd4730672966ad4de351e832f0f5d3277254601 100644
--- a/ChemistryLib/Output.h
+++ b/ChemistryLib/Output.h
@@ -37,8 +37,8 @@ public:
                display_distance + display_current_time + display_time_step;
     }
 
-    friend std::ofstream& operator<<(
-        std::ofstream& out, BasicOutputSetups const& basic_output_setups);
+    friend std::ostream& operator<<(
+        std::ostream& out, BasicOutputSetups const& basic_output_setups);
 
     std::string const output_file;
 
@@ -97,7 +97,7 @@ struct Output
         return matching_items;
     }
 
-    friend std::ofstream& operator<<(std::ofstream& out, Output const& output);
+    friend std::ostream& operator<<(std::ostream& os, Output const& output);
 
     BasicOutputSetups const basic_output_setups;
     std::vector<OutputItem> const accepted_items;
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index a07db023d84969827a87935b800ba8e54d846c31..75ed46671e957152cddfbeb614c08ffca742dc2c 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -9,7 +9,7 @@
 
 #include <boost/algorithm/string.hpp>
 #include <cmath>
-#include <fstream>
+#include <iostream>
 
 #include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTreeUtil.h"
@@ -19,6 +19,18 @@
 
 namespace ChemistryLib
 {
+namespace
+{
+template <typename DataBlock>
+std::ostream& operator<<(std::ostream& os,
+                         std::vector<DataBlock> const& data_blocks)
+{
+    std::copy(data_blocks.begin(), data_blocks.end(),
+              std::ostream_iterator<DataBlock>(os));
+    return os;
+}
+}  // namespace
+
 void PhreeqcIO::doWaterChemistryCalculation(
     std::vector<GlobalVector*>& process_solutions, double const dt)
 {
@@ -134,16 +146,16 @@ void PhreeqcIO::writeInputsToFile()
     out.close();
 }
 
-std::ofstream& operator<<(std::ofstream& out, PhreeqcIO const& phreeqc_io)
+std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
 {
-    out << "SELECTED_OUTPUT" << "\n";
-    out << *phreeqc_io._output << "\n";
+    os << "SELECTED_OUTPUT" << "\n";
+    os << *phreeqc_io._output << "\n";
 
     auto const& reaction_rates = phreeqc_io._reaction_rates;
     if (!reaction_rates.empty())
     {
-        out << "RATES" << "\n";
-        out << reaction_rates << "\n";
+        os << "RATES" << "\n";
+        os << reaction_rates << "\n";
     }
 
     std::size_t const num_chemical_systems =
@@ -154,30 +166,30 @@ std::ofstream& operator<<(std::ofstream& out, PhreeqcIO const& phreeqc_io)
     {
         auto const& aqueous_solution =
             phreeqc_io._aqueous_solutions[chemical_system_id];
-        out << "SOLUTION " << chemical_system_id + 1 << "\n";
-        out << aqueous_solution << "\n";
+        os << "SOLUTION " << chemical_system_id + 1 << "\n";
+        os << aqueous_solution << "\n";
 
         auto const& equilibrium_phases =
             phreeqc_io._equilibrium_phases[chemical_system_id];
         if (!equilibrium_phases.empty())
         {
-            out << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
-            out << equilibrium_phases << "\n";
+            os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
+            os << equilibrium_phases << "\n";
         }
 
         auto const& kinetic_reactants =
             phreeqc_io._kinetic_reactants[chemical_system_id];
         if (!kinetic_reactants.empty())
         {
-            out << "KINETICS " << chemical_system_id + 1 << "\n";
-            out << kinetic_reactants;
-            out << "-steps " << phreeqc_io._dt << "\n" << "\n";
+            os << "KINETICS " << chemical_system_id + 1 << "\n";
+            os << kinetic_reactants;
+            os << "-steps " << phreeqc_io._dt << "\n" << "\n";
         }
 
-        out << "END" << "\n" << "\n";
+        os << "END" << "\n" << "\n";
     }
 
-    return out;
+    return os;
 }
 
 void PhreeqcIO::execute()
@@ -231,7 +243,7 @@ void PhreeqcIO::readOutputsFromFile()
     in.close();
 }
 
-std::ifstream& operator>>(std::ifstream& in, PhreeqcIO& phreeqc_io)
+std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
 {
     // Skip the headline
     in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index 62d80eb31722707029d17d999bb007d36f2ea97c..8a856afcfc401ed8446e47b8e0fff46402fcf354 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -64,10 +64,10 @@ public:
 
     void readOutputsFromFile();
 
-    friend std::ofstream& operator<<(std::ofstream& out,
-                                     PhreeqcIO const& phreeqc_io);
+    friend std::ostream& operator<<(std::ostream& os,
+                                    PhreeqcIO const& phreeqc_io);
 
-    friend std::ifstream& operator>>(std::ifstream& in, PhreeqcIO& phreeqc_io);
+    friend std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io);
 
     std::string const _phreeqc_input_file;
 
diff --git a/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp b/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
index 932f3f6882b2da4078e5a33a9a07aef6be757c47..ccedd5ceb48519cc1e43627ed4b00e976852df4c 100644
--- a/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
+++ b/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
@@ -7,44 +7,44 @@
  *
  */
 
-#include <fstream>
+#include <ostream>
 
 #include "AqueousSolution.h"
 
 namespace ChemistryLib
 {
-std::ofstream& operator<<(std::ofstream& out,
-                          AqueousSolution const& aqueous_solution)
+std::ostream& operator<<(std::ostream& os,
+                         AqueousSolution const& aqueous_solution)
 {
-    out << "temp " << aqueous_solution.temperature << "\n";
+    os << "temp " << aqueous_solution.temperature << "\n";
 
-    out << "pressure " << aqueous_solution.pressure << "\n";
+    os << "pressure " << aqueous_solution.pressure << "\n";
 
     switch (aqueous_solution.means_of_adjusting_charge)
     {
         case MeansOfAdjustingCharge::pH:
-            out << "pH " << aqueous_solution.pH << " charge"
-                << "\n";
-            out << "pe " << aqueous_solution.pe << "\n";
+            os << "pH " << aqueous_solution.pH << " charge"
+               << "\n";
+            os << "pe " << aqueous_solution.pe << "\n";
             break;
         case MeansOfAdjustingCharge::pe:
-            out << "pH " << aqueous_solution.pH << "\n";
-            out << "pe " << aqueous_solution.pe << " charge"
-                << "\n";
+            os << "pH " << aqueous_solution.pH << "\n";
+            os << "pe " << aqueous_solution.pe << " charge"
+               << "\n";
             break;
         case MeansOfAdjustingCharge::Unspecified:
-            out << "pH " << aqueous_solution.pH << "\n";
-            out << "pe " << aqueous_solution.pe << "\n";
+            os << "pH " << aqueous_solution.pH << "\n";
+            os << "pe " << aqueous_solution.pe << "\n";
             break;
     }
 
-    out << "units mol/kgw\n";
+    os << "units mol/kgw\n";
 
     for (auto const& component : aqueous_solution.components)
     {
-        out << component.name << " " << component.amount << "\n";
+        os << component.name << " " << component.amount << "\n";
     }
 
-    return out;
+    return os;
 }
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/AqueousSolution.h b/ChemistryLib/PhreeqcIOData/AqueousSolution.h
index 403be34c2d0f19efcd0b5fa490cb06938c6b7076..1caaf0909691285e27dd570761990876466b80ff 100644
--- a/ChemistryLib/PhreeqcIOData/AqueousSolution.h
+++ b/ChemistryLib/PhreeqcIOData/AqueousSolution.h
@@ -48,8 +48,8 @@ struct AqueousSolution
     {
     }
 
-    friend std::ofstream& operator<<(std::ofstream& out,
-                                     AqueousSolution const& aqueous_solution);
+    friend std::ostream& operator<<(std::ostream& os,
+                                    AqueousSolution const& aqueous_solution);
 
     double temperature;
     double pressure;
diff --git a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp
index 243263a9af322fefd641d0d84d2add06aa3682e4..378eb2a4cdbc1779015a01c08e0f243f2f00c706 100644
--- a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp
+++ b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.cpp
@@ -7,24 +7,21 @@
  *
  */
 
-#include <fstream>
+#include <ostream>
 
 #include "EquilibriumPhase.h"
 
 namespace ChemistryLib
 {
-std::ofstream& operator<<(
-    std::ofstream& out, std::vector<EquilibriumPhase> const& equilibrium_phases)
+std::ostream& operator<<(std::ostream& os,
+                         EquilibriumPhase const& equilibrium_phase)
 {
-    for (auto const& equilibrium_phase : equilibrium_phases)
-    {
-        out << equilibrium_phase.name;
+    os << equilibrium_phase.name;
 
-        out << " " << equilibrium_phase.saturation_index;
+    os << " " << equilibrium_phase.saturation_index;
 
-        out << " " << equilibrium_phase.amount << "\n";
-    }
+    os << " " << equilibrium_phase.amount << "\n";
 
-    return out;
+    return os;
 }
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h
index b6eb94408aa0b31e12520b61a55835385b75950c..7851377236efa1172fc25e6c6b59f56d5b6fb59d 100644
--- a/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h
+++ b/ChemistryLib/PhreeqcIOData/EquilibriumPhase.h
@@ -33,9 +33,8 @@ struct EquilibriumPhase
     {
     }
 
-    friend std::ofstream& operator<<(
-        std::ofstream& out,
-        std::vector<EquilibriumPhase> const& equilibrium_phases);
+    friend std::ostream& operator<<(std::ostream& os,
+                                    EquilibriumPhase const& equilibrium_phase);
 
     std::string const name;
     double amount;
diff --git a/ChemistryLib/PhreeqcIOData/KineticReactant.cpp b/ChemistryLib/PhreeqcIOData/KineticReactant.cpp
index 58ca69348eb478d71d4d4e80ea188396d6085b74..8f2ab34fd521a6c1ffd6a5daafc6772d40393348 100644
--- a/ChemistryLib/PhreeqcIOData/KineticReactant.cpp
+++ b/ChemistryLib/PhreeqcIOData/KineticReactant.cpp
@@ -7,37 +7,34 @@
  *
  */
 
-#include <fstream>
+#include <ostream>
 
 #include "KineticReactant.h"
 
 namespace ChemistryLib
 {
-std::ofstream& operator<<(std::ofstream& out,
-                          std::vector<KineticReactant> const& kinetic_reactants)
+std::ostream& operator<<(std::ostream& os,
+                         KineticReactant const& kinetic_reactant)
 {
-    for (auto const& kinetic_reactant : kinetic_reactants)
-    {
-        out << kinetic_reactant.name << "\n";
+    os << kinetic_reactant.name << "\n";
 
-        if (!kinetic_reactant.chemical_formula.empty())
-        {
-            out << "-formula " << kinetic_reactant.chemical_formula << "\n";
-        }
+    if (!kinetic_reactant.chemical_formula.empty())
+    {
+        os << "-formula " << kinetic_reactant.chemical_formula << "\n";
+    }
 
-        out << "-m  " << kinetic_reactant.amount << "\n";
+    os << "-m  " << kinetic_reactant.amount << "\n";
 
-        if (!kinetic_reactant.parameters.empty())
+    if (!kinetic_reactant.parameters.empty())
+    {
+        os << "-parms";
+        for (auto const& parameter : kinetic_reactant.parameters)
         {
-            out << "-parms";
-            for (auto const& parameter : kinetic_reactant.parameters)
-            {
-                out << " " << parameter;
-            }
-            out << "\n";
+            os << " " << parameter;
         }
+        os << "\n";
     }
 
-    return out;
+    return os;
 }
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/KineticReactant.h b/ChemistryLib/PhreeqcIOData/KineticReactant.h
index 351c8e45aedac82fe1ef722fc385d3d69db30efa..ed8a983c5569220d11cc2ac02971e04257974019 100644
--- a/ChemistryLib/PhreeqcIOData/KineticReactant.h
+++ b/ChemistryLib/PhreeqcIOData/KineticReactant.h
@@ -31,9 +31,8 @@ struct KineticReactant
     {
     }
 
-    friend std::ofstream& operator<<(
-        std::ofstream& out,
-        std::vector<KineticReactant> const& kinetic_reactants);
+    friend std::ostream& operator<<(std::ostream& os,
+                                    KineticReactant const& kinetic_reactant);
 
     std::string const name;
     std::string const chemical_formula;
diff --git a/ChemistryLib/PhreeqcIOData/ReactionRate.cpp b/ChemistryLib/PhreeqcIOData/ReactionRate.cpp
index 4de068ca057dbfe41170e3c5fd4d2edbe4e00375..af613636fe6075eba5bf011945075ef2e2240751 100644
--- a/ChemistryLib/PhreeqcIOData/ReactionRate.cpp
+++ b/ChemistryLib/PhreeqcIOData/ReactionRate.cpp
@@ -7,29 +7,24 @@
  *
  */
 
-#include <fstream>
+#include <ostream>
 
 #include "ReactionRate.h"
 
 namespace ChemistryLib
 {
-std::ofstream& operator<<(std::ofstream& out,
-                          std::vector<ReactionRate> const& reaction_rates)
+std::ostream& operator<<(std::ostream& os, ReactionRate const& reaction_rate)
 {
-    for (auto const& reaction_rate : reaction_rates)
+    os << reaction_rate.kinetic_reactant << "\n";
+    os << "-start" << "\n";
+    int line_number = 1;
+    for (auto const& expression_statement : reaction_rate.expression_statements)
     {
-        out << reaction_rate.kinetic_reactant << "\n";
-        out << "-start" << "\n";
-        int line_number = 1;
-        for (auto const& expression_statement :
-             reaction_rate.expression_statements)
-        {
-            out << line_number << " " << expression_statement << "\n";
-            ++line_number;
-        }
-        out << "-end" << "\n";
+        os << line_number << " " << expression_statement << "\n";
+        ++line_number;
     }
+    os << "-end" << "\n";
 
-    return out;
+    return os;
 }
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/ReactionRate.h b/ChemistryLib/PhreeqcIOData/ReactionRate.h
index e894488c8b4c4a4a750f4709dd1b4688a440fa29..fe4bcf9a52cdd9ed0b18fab4cafbdfb47f4941cf 100644
--- a/ChemistryLib/PhreeqcIOData/ReactionRate.h
+++ b/ChemistryLib/PhreeqcIOData/ReactionRate.h
@@ -25,8 +25,8 @@ struct ReactionRate
     {
     }
 
-    friend std::ofstream& operator<<(
-        std::ofstream& out, std::vector<ReactionRate> const& reaction_rate);
+    friend std::ostream& operator<<(std::ostream& os,
+                                    ReactionRate const& reaction_rate);
 
     std::string const kinetic_reactant;
     std::vector<std::string> const expression_statements;