diff --git a/include/core/HY_Features.hpp b/include/core/HY_Features.hpp
index 7f1333c5fe..73a11fe587 100644
--- a/include/core/HY_Features.hpp
+++ b/include/core/HY_Features.hpp
@@ -190,6 +190,11 @@ namespace hy_features {
           std::cout<<"Catchment topology is dendritic."<<std::endl;
         }
 
+        /**
+         * Release formulations
+        */
+        void release_formulations();
+
         /**
          * @brief Destroy the hy features object
          * 
diff --git a/include/core/catchment/HY_CatchmentArea.hpp b/include/core/catchment/HY_CatchmentArea.hpp
index 603780f7a1..d1e13ec228 100644
--- a/include/core/catchment/HY_CatchmentArea.hpp
+++ b/include/core/catchment/HY_CatchmentArea.hpp
@@ -17,15 +17,17 @@ class HY_CatchmentArea : public HY_CatchmentRealization, public GM_Object
     HY_CatchmentArea();
     HY_CatchmentArea(std::shared_ptr<data_access::GenericDataProvider> forcing, utils::StreamHandler output_stream);
     //HY_CatchmentArea(forcing_params forcing_config, utils::StreamHandler output_stream); //TODO not sure I like this pattern
-    void set_output_stream(std::string file_path){output = utils::FileStreamHandler(file_path.c_str());}
-    void write_output(std::string out){ output<<out; }
+    void set_output_stream(std::string file_path){
+        output = std::make_unique<utils::FileStreamHandler>(file_path.c_str());
+    }
+    void write_output(std::string out){ *output<<out; }
     virtual ~HY_CatchmentArea();
 
     protected:
 
     polygon_t bounds;
 
-    utils::StreamHandler output;
+    std::unique_ptr<utils::StreamHandler> output;
 
     private:
 };
diff --git a/include/realizations/catchment/Bmi_Multi_Formulation.hpp b/include/realizations/catchment/Bmi_Multi_Formulation.hpp
index 3ed5c17749..9df11fc0e9 100644
--- a/include/realizations/catchment/Bmi_Multi_Formulation.hpp
+++ b/include/realizations/catchment/Bmi_Multi_Formulation.hpp
@@ -659,7 +659,7 @@ namespace realization {
         template<class T>
         std::shared_ptr<T> init_nested_module(int mod_index, std::string identifier, geojson::PropertyMap properties) {
             std::shared_ptr<data_access::GenericDataProvider> wfp = std::make_shared<data_access::WrappedDataProvider>(this);
-            std::shared_ptr<T> mod = std::make_shared<T>(identifier, wfp, output);
+            std::shared_ptr<T> mod = std::make_shared<T>(identifier, wfp, *output);
 
             // Since this is a nested formulation, support usage of the '{{id}}' syntax for init config file paths.
             Catchment_Formulation::config_pattern_substitution(properties, BMI_REALIZATION_CFG_PARAM_REQ__INIT_CONFIG,
diff --git a/include/realizations/catchment/Formulation_Manager.hpp b/include/realizations/catchment/Formulation_Manager.hpp
index a61c68833b..998062bdf1 100644
--- a/include/realizations/catchment/Formulation_Manager.hpp
+++ b/include/realizations/catchment/Formulation_Manager.hpp
@@ -27,6 +27,9 @@ namespace realization {
 
         public:
 
+            void clear(){
+                formulations.clear();            
+            }
             std::shared_ptr<Simulation_Time> Simulation_Time_Object;
 
             Formulation_Manager(std::stringstream &data) {
diff --git a/include/utilities/FileStreamHandler.hpp b/include/utilities/FileStreamHandler.hpp
index 1e4108fbaa..1ee36804f9 100644
--- a/include/utilities/FileStreamHandler.hpp
+++ b/include/utilities/FileStreamHandler.hpp
@@ -14,7 +14,13 @@ namespace utils
                 stream->open(path, std::ios::trunc);
                 output_stream = stream;
             }
-            virtual ~FileStreamHandler(){}
+
+            virtual ~FileStreamHandler() override {
+                auto tmp = std::static_pointer_cast<std::ofstream>(output_stream);
+                if(tmp){
+                    tmp->close();
+                }
+            }
     };
 
 }
diff --git a/include/utilities/StreamHandler.hpp b/include/utilities/StreamHandler.hpp
index 3f2b45f370..5c1fdf2db8 100644
--- a/include/utilities/StreamHandler.hpp
+++ b/include/utilities/StreamHandler.hpp
@@ -39,7 +39,7 @@ namespace utils
 
             /** Deconstructor for a StreamHandler */
 
-            ~StreamHandler()
+            virtual ~StreamHandler()
             {
 
             }
diff --git a/src/NGen.cpp b/src/NGen.cpp
index 74817737db..e80134cf77 100644
--- a/src/NGen.cpp
+++ b/src/NGen.cpp
@@ -539,6 +539,17 @@ int main(int argc, char *argv[]) {
 
     } //done time
 
+  //Close down any file handles we have
+  for(auto& f : nexus_outfiles){
+    f.second.close();
+  }
+  //std::cout<<"Clearning manger with "<<manager.use_count()<<" references\n";
+  //Need to release all formulation references, those stored in hy_features object and the
+  //references used by the formulation manager
+  //this will ensure all formulation destructors are called and resources released.
+  features.release_formulations();
+  manager->clear();
+
 #if NGEN_MPI_ACTIVE
     MPI_Barrier(MPI_COMM_WORLD);
 #endif
diff --git a/src/core/HY_Features.cpp b/src/core/HY_Features.cpp
index 91e59cb487..86aa7bdeb4 100644
--- a/src/core/HY_Features.cpp
+++ b/src/core/HY_Features.cpp
@@ -16,6 +16,11 @@ HY_Features::HY_Features( geojson::GeoJSON catchments, std::string* link_key, st
     HY_Features(network::Network(catchments, link_key), formulations, catchments){  
 }
 
+void HY_Features::release_formulations(){
+  formulations->clear();
+  _catchments.clear();
+}
+
 HY_Features::HY_Features(network::Network network, std::shared_ptr<Formulation_Manager> formulations, geojson::GeoJSON fabric)
   :network(network), formulations(formulations)
 {
diff --git a/src/core/catchment/HY_CatchmentArea.cpp b/src/core/catchment/HY_CatchmentArea.cpp
index 6fd22037d6..e1ffac1348 100644
--- a/src/core/catchment/HY_CatchmentArea.cpp
+++ b/src/core/catchment/HY_CatchmentArea.cpp
@@ -5,7 +5,7 @@ HY_CatchmentArea::HY_CatchmentArea()
     //ctor
 }
 
-HY_CatchmentArea::HY_CatchmentArea(std::shared_ptr<data_access::GenericDataProvider> forcing, utils::StreamHandler output_stream) : HY_CatchmentRealization(forcing), output(output_stream) { }
+HY_CatchmentArea::HY_CatchmentArea(std::shared_ptr<data_access::GenericDataProvider> forcing, utils::StreamHandler output_stream) : HY_CatchmentRealization(forcing), output(std::make_unique<utils::StreamHandler>(output_stream)) { }
 
 HY_CatchmentArea::~HY_CatchmentArea()
 {
diff --git a/src/realizations/catchment/Bmi_C_Formulation.cpp b/src/realizations/catchment/Bmi_C_Formulation.cpp
index 04a0106795..a3323b4212 100644
--- a/src/realizations/catchment/Bmi_C_Formulation.cpp
+++ b/src/realizations/catchment/Bmi_C_Formulation.cpp
@@ -32,7 +32,7 @@ std::shared_ptr<Bmi_C_Adapter> Bmi_C_Formulation::construct_model(const geojson:
                     get_allow_model_exceed_end_time(),
                     is_bmi_model_time_step_fixed(),
                     reg_func,
-                    output);
+                    *output);
 }
 
 std::string Bmi_C_Formulation::get_output_header_line(std::string delimiter) {
diff --git a/src/realizations/catchment/Bmi_Cpp_Formulation.cpp b/src/realizations/catchment/Bmi_Cpp_Formulation.cpp
index 24419b2c75..bb35b93794 100644
--- a/src/realizations/catchment/Bmi_Cpp_Formulation.cpp
+++ b/src/realizations/catchment/Bmi_Cpp_Formulation.cpp
@@ -42,7 +42,7 @@ std::shared_ptr<Bmi_Cpp_Adapter> Bmi_Cpp_Formulation::construct_model(const geoj
                     is_bmi_model_time_step_fixed(),
                     model_create_fname,
                     model_destroy_fname,
-                    output);
+                    *output);
 }
 
 std::string Bmi_Cpp_Formulation::get_output_header_line(std::string delimiter) {
diff --git a/src/realizations/catchment/Bmi_Fortran_Formulation.cpp b/src/realizations/catchment/Bmi_Fortran_Formulation.cpp
index bd6101824f..5439185858 100644
--- a/src/realizations/catchment/Bmi_Fortran_Formulation.cpp
+++ b/src/realizations/catchment/Bmi_Fortran_Formulation.cpp
@@ -38,7 +38,7 @@ std::shared_ptr<Bmi_Fortran_Adapter> Bmi_Fortran_Formulation::construct_model(co
             get_allow_model_exceed_end_time(),
             is_bmi_model_time_step_fixed(),
             reg_func,
-            output);
+            *output);
 }
 
 std::string Bmi_Fortran_Formulation::get_formulation_type() {
diff --git a/src/realizations/catchment/Bmi_Py_Formulation.cpp b/src/realizations/catchment/Bmi_Py_Formulation.cpp
index cc69a80a32..a347719ed1 100644
--- a/src/realizations/catchment/Bmi_Py_Formulation.cpp
+++ b/src/realizations/catchment/Bmi_Py_Formulation.cpp
@@ -28,7 +28,7 @@ shared_ptr<Bmi_Py_Adapter> Bmi_Py_Formulation::construct_model(const geojson::Pr
                     python_type_name,
                     get_allow_model_exceed_end_time(),
                     is_bmi_model_time_step_fixed(),
-                    output);
+                    *output);
 }
 
 time_t realization::Bmi_Py_Formulation::convert_model_time(const double &model_time) {