diff --git a/include/bmi_serialization.h b/include/bmi_serialization.h
index 2ceb763..b1e776c 100644
--- a/include/bmi_serialization.h
+++ b/include/bmi_serialization.h
@@ -8,7 +8,7 @@ extern "C" {
 #include "bmi.h"
 
 const int serialize_topmodel(Bmi* bmi);
-const int deserialize_topmodel(Bmi* bmi, const char* buffer);
+const int deserialize_topmodel(Bmi* bmi, char* buffer);
 
 #ifdef __cplusplus
 }
diff --git a/include/topmodel.h b/include/topmodel.h
index b6333da..6fe0306 100755
--- a/include/topmodel.h
+++ b/include/topmodel.h
@@ -26,28 +26,28 @@
 
 /*** Function/subroutine prototypes ***/
 extern void convert_dist_to_histords(const double * const dist_from_outlet, const int num_channels,
-					const double * const chv, const double * const rv, const double dt, double* const tch);
+					const double chv, const double rv, const double dt, double* const tch);
 
 extern void calc_time_delay_histogram(int num_channels, double area,
-		double* tch, double *cum_dist_area_with_dist,
+		const double *tch, const double *cum_dist_area_with_dist,
 		int *num_time_delay_histo_ords, int *num_delay,	double **time_delay_histogram);
 
 extern void init_water_balance(
-				int num_topodex_values, double dt, double *sr0, 
-				double *szm, double *Q0, double *t0, double tl,
+				int num_topodex_values, double dt, double sr0, 
+				double szm, double Q0, double t0, double tl,
 				double **stor_unsat_zone, double *szq,
 				double **deficit_local, double **deficit_root_zone, 
 				double *sbar, double *bal);
 
-extern void init_discharge_array(int stand_alone, int *num_delay, double *Q0, double area, 
-			int *num_time_delay_histo_ords, double **time_delay_histogram,
+extern void init_discharge_array(int stand_alone, int num_delay, double Q0, double area, 
+			int num_time_delay_histo_ords, double **time_delay_histogram,
                         double **Q);
 
 extern int init(FILE *in_param_fptr, FILE *output_fptr, char *subcat, int stand_alone,
 	      int num_channels, int num_topodex_values, int yes_print_output,
 	      double area, double **time_delay_histogram,
 	      double *cum_dist_area_with_dist, double dt, 
-        double tl, double *dist_from_outlet, 
+        double tl, const double *dist_from_outlet, 
 	      int *num_time_delay_histo_ords,int *num_delay,
 	      double *szm, double *t0, double *chv, double *rv, double *td, double *srmax, 
         double *Q0,double *sr0, int *infex, double *xk0, double *hf, double *dth,
@@ -61,11 +61,11 @@ extern int inputs(FILE *input_fptr, int *nstep, double *dt, double **rain,
 extern void topmod(FILE *output_fptr, int nstep, int num_topodex_values,
                 int yes_print_output,int infex, double dt, double szm,
 	        double *stor_unsat_zone, double *deficit_root_zone,
-                double *deficit_local, double *pe, double *rain,double xk0,double hf, 
-                double *dist_area_lnaotb, double tl, double *lnaotb, double td,
+                double *deficit_local, const double *pe, const double *rain,double xk0,double hf, 
+                const double *dist_area_lnaotb, double tl, const double *lnaotb, double td,
                 double srmax, double *contrib_area, double szq, double *Qout, 
                 int num_time_delay_histo_ords,double *Q,
-                double *time_delay_histogram,char *subcat,double *bal,
+                const double *time_delay_histogram,
                 double *sbar,int num_delay, int current_time_step, int stand_alone,
                 double *sump, double *sumae, double *sumq, double *sumrz, double *sumuz,
                 double *quz, double *qb, double *qof, double *p, double *ep);
diff --git a/include/vecbuf.hpp b/include/vecbuf.hpp
index 6db982c..a20bc70 100644
--- a/include/vecbuf.hpp
+++ b/include/vecbuf.hpp
@@ -28,6 +28,9 @@ class vecbuf : public std::basic_streambuf<CharT> {
     // Forwarder for std::vector::clear()
     constexpr void clear() { vector_.clear(); }
 
+    // Forwarder for std::vector::resize(size)
+    constexpr void resize(size_type size) { vector_.resize(size); }
+
     // Forwarder for std::vector::reserve
     constexpr void reserve(size_type capacity) { vector_.reserve(capacity); setp_from_vector(); }
 
@@ -35,7 +38,7 @@ class vecbuf : public std::basic_streambuf<CharT> {
     constexpr void reserve_additional(size_type additional_capacity) { reserve(size() + additional_capacity); }
 
     // Forwarder for std::vector::data
-    constexpr const value_type* data() const { return vector_.data(); }
+    constexpr value_type* data() { return vector_.data(); }
 
     // Forwarder for std::vector::size
     constexpr size_type size() const { return vector_.size(); }
@@ -112,4 +115,11 @@ class vecbuf : public std::basic_streambuf<CharT> {
 
 };
 
+class membuf : public std::streambuf {
+public:
+    membuf(char *begin, size_t size) {
+        this->setg(begin, begin, begin + size);
+    }
+};
+
 #endif
diff --git a/src/bmi_serialization.cpp b/src/bmi_serialization.cpp
index 03a3448..4daa17d 100644
--- a/src/bmi_serialization.cpp
+++ b/src/bmi_serialization.cpp
@@ -32,6 +32,12 @@ class TopmodelSerializer {
 template<class Archive>
 void TopmodelSerializer::serialize(Archive& ar, const unsigned int version) {
     topmodel_model* model = this->model;
+    if (model->stand_alone == TRUE) {
+        // the number of timesteps makes hindcasting nigh imposible when stand alone
+        auto error = "Topmodel serialization is not currently implemented when running stand alone.";
+        Log(SEVERE, error);
+        throw std::runtime_error(error);
+    }
     ar & model->current_time_step;
 
     // data summed between runs
@@ -62,12 +68,37 @@ void TopmodelSerializer::serialize(Archive& ar, const unsigned int version) {
     ar & boost::serialization::make_array(
         model->deficit_local, num_topodex_values
     );
+
+    // nsteps will always be 1 for non-stand-alone models
     ar & boost::serialization::make_array(
-        model->contrib_area, num_topodex_values
+        model->contrib_area, model->nstep + 1
     );
+
+    // copy the current sizes to detect changes, then archive the model value
+    int num_time_delay_histo_ords = model->num_time_delay_histo_ords;
+    ar & model->num_time_delay_histo_ords;
+    int num_delay = model->num_delay;
+    ar & model->num_delay;
+    size_t num_Q = model->num_delay + model->num_time_delay_histo_ords + 1;
+    if (Archive::is_loading::value) {
+        // if loading and array size has changed, reallocate
+        if (num_time_delay_histo_ords != model->num_time_delay_histo_ords) {
+            if (model->time_delay_histogram != NULL)
+                free(model->time_delay_histogram);
+            model->time_delay_histogram = (double *)malloc(
+                (model->num_time_delay_histo_ords + 1) * sizeof(double)
+            );
+        }
+        if (num_delay != model->num_delay || num_time_delay_histo_ords != model->num_time_delay_histo_ords) {
+            if (model->Q != NULL)
+                free(model->Q);
+            model->Q = (double *)malloc(num_Q * sizeof(double));
+        }
+    }
     ar & boost::serialization::make_array(
-        model->Q, model->num_time_delay_histo_ords + 1
+        model->time_delay_histogram, model->num_time_delay_histo_ords + 1
     );
+    ar & boost::serialization::make_array(model->Q, num_Q);
 }
 
 
@@ -99,15 +130,17 @@ const int serialize_topmodel(Bmi* bmi) {
         free(model->serialized);
     }
     // set size and allocate memory
-    model->serialized_length = stream.size();
-    model->serialized = (char*)malloc(sizeof(char) * model->serialized_length);
+    uint64_t serialized_size = stream.size();
+    model->serialized_length = serialized_size + sizeof(uint64_t);
+    model->serialized = (char*)malloc(model->serialized_length);
     // make sure memory could be allocated
     if (model->serialized == NULL) {
         model->serialized_length = 0;
         return BMI_FAILURE;
     }
     // copy stream data to new allocation
-    memcpy(model->serialized, stream.data(), model->serialized_length);
+    memcpy(model->serialized, &serialized_size, sizeof(uint64_t));
+    memcpy(model->serialized + sizeof(uint64_t), stream.data(), serialized_size);
     return BMI_SUCCESS;
 }
 
@@ -118,9 +151,13 @@ const int serialize_topmodel(Bmi* bmi) {
   * @param buffer Start of data that wil be read as previously serialized state
   * @return int signifiying whether the serialization process completed successfully.
   */
-const int deserialize_topmodel(Bmi* bmi, const char* buffer) {
+const int deserialize_topmodel(Bmi* bmi, char* buffer) {
     TopmodelSerializer serializer(bmi);
-    std::istringstream stream(buffer);
+    // copy size of data out of header
+    uint64_t size;
+    memcpy(&size, buffer, sizeof(uint64_t));
+    // create stream from data after header
+    membuf stream(buffer + sizeof(uint64_t), size);
     boost::archive::binary_iarchive archive(stream);
     try {
         archive >> serializer;
diff --git a/src/bmi_topmodel.c b/src/bmi_topmodel.c
index 54f3256..f3fb362 100755
--- a/src/bmi_topmodel.c
+++ b/src/bmi_topmodel.c
@@ -481,8 +481,6 @@ static int Update(Bmi *self) {
         topmodel->num_time_delay_histo_ords,
         topmodel->Q,
         topmodel->time_delay_histogram,
-        topmodel->subcat,
-        &topmodel->bal,
         &topmodel->sbar,
         topmodel->num_delay,
         topmodel->current_time_step,
@@ -652,6 +650,9 @@ static int Get_var_type(Bmi *self, const char *name, char *type) {
     } else if (strcmp(name, "serialization_free") == 0) {
         strncpy(type, "int", BMI_MAX_TYPE_NAME);
         return BMI_SUCCESS;
+    } else if (strcmp(name, "reset_time") == 0) {
+        strncpy(type, "double", BMI_MAX_TYPE_NAME);
+        return BMI_SUCCESS;
     }
     // If we get here, it means the variable name wasn't recognized
     type[0] = '\0';
@@ -776,7 +777,10 @@ static int Get_var_nbytes(Bmi *self, const char *name, int *nbytes) {
     }
     // special cases for save state
     if (item_count < 1) {
-        if (strcmp(name, "serialization_create") == 0 || strcmp(name, "serialization_size") == 0 || strcmp(name, "serialization_free") == 0) {
+        if (strcmp(name, "serialization_create") == 0
+            || strcmp(name, "serialization_size") == 0
+            || strcmp(name, "serialization_free") == 0
+            || strcmp(name, "reset_time") == 0) {
             item_count = 1;
         } else if (strcmp(name, "serialization_state") == 0) {
             topmodel_model* model = (topmodel_model*)self->data;
@@ -1063,6 +1067,11 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         } else {
             return BMI_FAILURE;
         }
+    } else if (strcmp(name, "reset_time") == 0) {
+        topmodel_model* model = (topmodel_model *)self->data;
+        // current_time_step is mainly used for indexing into config data, so should be safe to reset and nothing else
+        model->current_time_step = 0;
+        return BMI_SUCCESS;
     }
 
     if (self->get_value_ptr(self, name, &dest) == BMI_FAILURE)
@@ -1168,8 +1177,8 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         convert_dist_to_histords(
             topmodel->dist_from_outlet,
             topmodel->num_channels,
-            &topmodel->chv,
-            &topmodel->rv,
+            topmodel->chv,
+            topmodel->rv,
             topmodel->dt,
             tch
         );
@@ -1188,10 +1197,10 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         // Reinitialise discharge array
         init_discharge_array(
             topmodel->stand_alone,
-            &topmodel->num_delay,
-            &topmodel->Q0,
+            topmodel->num_delay,
+            topmodel->Q0,
             topmodel->area,
-            &topmodel->num_time_delay_histo_ords,
+            topmodel->num_time_delay_histo_ords,
             &topmodel->time_delay_histogram,
             &topmodel->Q
         );
@@ -1211,10 +1220,10 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         init_water_balance(
             topmodel->num_topodex_values,
             topmodel->dt,
-            &topmodel->sr0,
-            &topmodel->szm,
-            &topmodel->Q0,
-            &topmodel->t0,
+            topmodel->sr0,
+            topmodel->szm,
+            topmodel->Q0,
+            topmodel->t0,
             topmodel->tl,
             &topmodel->stor_unsat_zone,
             &topmodel->szq,
diff --git a/src/topmodel.c b/src/topmodel.c
index 335931b..4eeea89 100755
--- a/src/topmodel.c
+++ b/src/topmodel.c
@@ -119,13 +119,13 @@ extern void topmod(
     double *stor_unsat_zone,
     double *deficit_root_zone,
     double *deficit_local,
-    double *pe,
-    double *rain,
+    const double *pe,
+    const double *rain,
     double xk0,
     double hf,
-    double *dist_area_lnaotb,
+    const double *dist_area_lnaotb,
     double tl,
-    double *lnaotb,
+    const double *lnaotb,
     double td,
     double srmax,
     double *contrib_area,
@@ -133,9 +133,7 @@ extern void topmod(
     double *Qout,
     int num_time_delay_histo_ords,
     double *Q,
-    double *time_delay_histogram,
-    char *subcat,
-    double *bal,
+    const double *time_delay_histogram,
     double *sbar,
     int num_delay,
     int current_time_step,
@@ -689,8 +687,8 @@ extern int tread(
 extern void convert_dist_to_histords(
     const double *const dist_from_outlet,
     const int num_channels,
-    const double *const chv,
-    const double *const rv,
+    const double chv,
+    const double rv,
     const double dt,
     double *const tch
 ) {
@@ -705,8 +703,8 @@ extern void convert_dist_to_histords(
     double chvdt, rvdt;
     int j;
 
-    chvdt = *chv * dt; // distance water travels in one timestep within channel
-    rvdt  = *rv * dt; // distance water travels as overland flow in one timestep
+    chvdt = chv * dt; // distance water travels in one timestep within channel
+    rvdt  = rv * dt; // distance water travels as overland flow in one timestep
 
     tch[1] = dist_from_outlet[1] / chvdt;
     for (j = 2; j <= num_channels; j++) {
@@ -741,8 +739,8 @@ extern void convert_dist_to_histords(
 extern void calc_time_delay_histogram(
     int num_channels,
     double area,
-    double *tch,
-    double *cum_dist_area_with_dist,
+    const double *tch,
+    const double *cum_dist_area_with_dist,
     int *num_time_delay_histo_ords,
     int *num_delay,
     double **time_delay_histogram
@@ -861,10 +859,10 @@ extern void calc_time_delay_histogram(
  */
 extern void init_discharge_array(
     int stand_alone,
-    int *num_delay,
-    double *Q0,
+    int num_delay,
+    double Q0,
     double area,
-    int *num_time_delay_histo_ords,
+    int num_time_delay_histo_ords,
     double **time_delay_histogram,
     double **Q
 ) {
@@ -876,7 +874,7 @@ extern void init_discharge_array(
             *Q = NULL;
         }
         //*Q = calloc(*num_delay + *num_time_delay_histo_ords + 1, sizeof(double));
-        d_alloc(Q, *num_delay + *num_time_delay_histo_ords);
+        d_alloc(Q, num_delay + num_time_delay_histo_ords);
     }
 
     // declare local variables
@@ -885,14 +883,14 @@ extern void init_discharge_array(
 
     sum = 0.0;
 
-    for (i = 1; i <= (*num_delay); i++) {
-        (*Q)[i] += (*Q0) * area;
+    for (i = 1; i <= (num_delay); i++) {
+        (*Q)[i] += Q0 * area;
     }
 
-    for (i = 1; i <= (*num_time_delay_histo_ords); i++) {
+    for (i = 1; i <= num_time_delay_histo_ords; i++) {
         sum += (*time_delay_histogram)[i];
-        in = (*num_delay) + i;
-        (*Q)[in] += (*Q0) * (area - sum);
+        in = num_delay + i;
+        (*Q)[in] += Q0 * (area - sum);
     };
 
     return;
@@ -931,10 +929,10 @@ extern void init_discharge_array(
 extern void init_water_balance(
     int num_topodex_values,
     double dt,
-    double *sr0,
-    double *szm,
-    double *Q0,
-    double *t0,
+    double sr0,
+    double szm,
+    double Q0,
+    double t0,
     double tl,
     double **stor_unsat_zone,
     double *szq,
@@ -960,20 +958,20 @@ extern void init_water_balance(
     // document the assumption and the requirement for caller to size these arrays to
     // num_topodex_values
 
-    t0dt = (*t0) + log(dt); /* was ALOG - specific log function in fortran*/
+    t0dt = t0 + log(dt); /* was ALOG - specific log function in fortran*/
 
     /*  Calculate SZQ parameter */
     (*szq) = exp(t0dt - tl);
 
     for (ia = 1; ia <= num_topodex_values; ia++) {
         (*stor_unsat_zone)[ia]   = 0.0;
-        (*deficit_root_zone)[ia] = (*sr0);
+        (*deficit_root_zone)[ia] = sr0;
     }
 
-    (*sbar) = -(*szm) * log((*Q0) / (*szq));
+    (*sbar) = (-szm) * log(Q0 / (*szq));
 
     //  Initialise water balance.  BAL is positive for storage
-    (*bal) = -(*sbar) - (*sr0);
+    (*bal) = -(*sbar) - sr0;
 
     return;
 }
@@ -1072,7 +1070,7 @@ extern int init(
     double *cum_dist_area_with_dist,
     double dt,
     double tl,
-    double *dist_from_outlet,
+    const double *dist_from_outlet,
     int *num_time_delay_histo_ords,
     int *num_delay,
     double *szm,
@@ -1149,7 +1147,7 @@ extern int init(
 
     // NJF num_channels is the value provided (SHOULD COME FROM TREAD)
     // Convert distance/area form to time delay histogram ordinates
-    convert_dist_to_histords(dist_from_outlet, num_channels, chv, rv, dt, tch);
+    convert_dist_to_histords(dist_from_outlet, num_channels, *chv, *rv, dt, tch);
 
     // calculate the time_delay_histogram
     calc_time_delay_histogram(
@@ -1165,10 +1163,10 @@ extern int init(
     // Reinitialise discharge array
     init_discharge_array(
         stand_alone,
-        num_delay,
-        Q0,
+        *num_delay,
+        *Q0,
         area,
-        num_time_delay_histo_ords,
+        *num_time_delay_histo_ords,
         time_delay_histogram,
         Q
     );
@@ -1177,10 +1175,10 @@ extern int init(
     init_water_balance(
         num_topodex_values,
         dt,
-        sr0,
-        szm,
-        Q0,
-        t0,
+        *sr0,
+        *szm,
+        *Q0,
+        *t0,
         tl,
         stor_unsat_zone,
         szq,