diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h
index 7ee660d8bd..c7f6413b09 100644
--- a/include/openmc/volume_calc.h
+++ b/include/openmc/volume_calc.h
@@ -8,6 +8,7 @@
 
 #include <string>
 #include <vector>
+#include <gsl/gsl>
 
 namespace openmc {
 
@@ -23,6 +24,31 @@ public:
     std::vector<int> nuclides; //!< Index of nuclides
     std::vector<double> atoms; //!< Number of atoms for each nuclide
     std::vector<double> uncertainty; //!< Uncertainty on number of atoms
+    size_t num_samples;
+
+    Result& operator +=( const Result& other) {
+      Expects(volume.size() == other.volume.size());
+      Expects(atoms.size() == atoms.size());
+
+      size_t total_samples = num_samples + other.num_samples;
+
+      for (int i = 0; i < volume.size(); i++) {
+        // average volume results
+        volume[0] = (num_samples * volume[0] + other.num_samples * other.volume[0]) / total_samples;
+        // propagate error
+        volume[1] = std::sqrt(num_samples * volume[1] *volume[1] + other.num_samples * other.volume[1] * other.volume[1]) / total_samples;
+      }
+
+      for (int i = 0; i < atoms.size(); i++) {
+        atoms[i] = (num_samples * atoms[i] + other.num_samples * other.atoms[i]) / total_samples;
+        uncertainty[i] = std::sqrt(num_samples * uncertainty[i] * uncertainty[i] + other.num_samples * other.uncertainty[i] * other.uncertainty[i]) / total_samples;
+      }
+
+      num_samples = total_samples;
+      
+      return *this;
+
+    }
   }; // Results for a single domain
 
   // Constructors
@@ -36,6 +62,8 @@ public:
   //! \return Vector of results for each user-specified domain
   std::vector<Result> execute() const;
 
+  std::vector<Result> _execute(size_t seed_offset = 0) const;
+
   //! \brief Write volume calculation results to HDF5 file
   //
   //! \param[in] filename  Path to HDF5 file to write
@@ -45,6 +73,7 @@ public:
   // Data members
   int domain_type_; //!< Type of domain (cell, material, etc.)
   int n_samples_; //!< Number of samples to use
+  int seed_offset_;
   Position lower_left_; //!< Lower-left position of bounding box
   Position upper_right_; //!< Upper-right position of bounding box
   std::vector<int> domain_ids_; //!< IDs of domains to find volumes of
diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp
index 6c1c48932f..bd6ada96a7 100644
--- a/src/volume_calc.cpp
+++ b/src/volume_calc.cpp
@@ -71,7 +71,40 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node)
 
 }
 
-std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
+std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const {
+
+  std::vector<VolumeCalculation::Result> results;
+  size_t offset = 0;
+
+  results = _execute(offset);
+  offset += n_samples_;
+
+  double max_err = -INFTY;
+  for (int i = 0; i < results.size(); i++) {
+    max_err = std::max(max_err, results[i].volume[1]);
+  }
+
+  double error_limit = 1E-05;
+  int iters = 1;
+  while (max_err > error_limit) {
+    std::cout << "Iter " << iters++ << std::endl;
+    std::vector<VolumeCalculation::Result> tmp = _execute(offset);
+    max_err = -INFTY;
+    for (int i = 0; i < results.size(); i++) {
+      auto& result = results[i];
+      result += tmp[i];
+      max_err = std::max(max_err, result.volume[1]);
+    }
+    
+    offset += n_samples_;
+
+    std::cout << "Max error: " << max_err << std::endl;
+  }
+
+  return results;
+}
+
+std::vector<VolumeCalculation::Result> VolumeCalculation::_execute(size_t seed_offset) const
 {
   // Shared data that is collected from all threads
   int n = domain_ids_.size();
@@ -102,7 +135,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
     // Sample locations and count hits
     #pragma omp for
     for (int i = i_start; i < i_end; i++) {
-      set_particle_seed(i);
+      set_particle_seed(seed_offset + i);
 
       p.n_coord_ = 1;
       Position xi {prn(), prn(), prn()};
@@ -248,6 +281,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
       result.volume[0] = static_cast<double>(total_hits) / n_samples_ * volume_sample;
       result.volume[1] = std::sqrt(result.volume[0]
         * (volume_sample - result.volume[0]) / n_samples_);
+      result.num_samples = n_samples_;
 
       for (int j = 0; j < n_nuc; ++j) {
         // Determine total number of atoms. At this point, we have values in