diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst
index 6e9a7f2ab..915f698b6 100644
--- a/docs/source/io_formats/settings.rst
+++ b/docs/source/io_formats/settings.rst
@@ -902,7 +902,12 @@ attributes/sub-elements:
 
 The ``<surf_source_write>`` element triggers OpenMC to bank particles crossing
 certain surfaces and write out the source bank in a separate file called
-``surface_source.h5``. This element has the following attributes/sub-elements:
+``surface_source.h5``. One or multiple surface IDs and one cell ID can be used
+to select the surfaces of interest. If no surface IDs are declared, every surface
+of the model is eligible to bank particles. In that case, a cell ID (using
+either the ``cell``, ``cellfrom`` or ``cellto`` attributes) can be used to select
+every surface of a specific cell. This element has the following
+attributes/sub-elements:
 
   :surface_ids:
     A list of integers separated by spaces indicating the unique IDs of surfaces
@@ -928,6 +933,34 @@ certain surfaces and write out the source bank in a separate file called
 
     .. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf
 
+  :cell:
+    An integer representing the cell ID used to determine if particles crossing
+    identified surfaces are to be banked. Particles coming from or going to this
+    declared cell will be banked if they cross the identified surfaces.
+
+    *Default*: None
+
+  :cellfrom:
+    An integer representing the cell ID used to determine if particles crossing
+    identified surfaces are to be banked. Particles coming from this declared
+    cell will be banked if they cross the identified surfaces.
+
+    *Default*: None
+
+  :cellto:
+    An integer representing the cell ID used to determine if particles crossing
+    identified surfaces are to be banked. Particles going to this declared cell
+    will be banked if they cross the identified surfaces.
+
+    *Default*: None
+
+.. note:: The ``cell``, ``cellfrom`` and ``cellto`` attributes cannot be
+          used simultaneously.
+
+.. note:: Surfaces with boundary conditions that are not "transmission" or "vacuum"
+          are not eligible to store any particles when using ``cell``, ``cellfrom``
+          or ``cellto`` attributes. It is recommended to use surface IDs instead.
+
 ------------------------------
 ``<survival_biasing>`` Element
 ------------------------------
diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst
index eb02f654c..d73d90cbb 100644
--- a/docs/source/usersguide/settings.rst
+++ b/docs/source/usersguide/settings.rst
@@ -277,6 +277,9 @@ source file can be manually generated with the :func:`openmc.write_source_file`
 function. This is particularly useful for coupling OpenMC with another program
 that generates a source to be used in OpenMC.
 
+Surface Sources
++++++++++++++++
+
 A source file based on particles that cross one or more surfaces can be
 generated during a simulation using the :attr:`Settings.surf_source_write`
 attribute::
@@ -287,7 +290,51 @@ attribute::
   }
 
 In this example, at most 10,000 source particles are stored when particles cross
-surfaces with IDs of 1, 2, or 3.
+surfaces with IDs of 1, 2, or 3. If no surface IDs are declared, particles
+crossing any surface of the model will be banked::
+
+  settings.surf_source_write = {'max_particles': 10000}
+
+A cell ID can also be used to bank particles that are crossing any surface of
+a cell that particles are either coming from or going to::
+
+  settings.surf_source_write = {'cell': 1, 'max_particles': 10000}
+
+In this example, particles that are crossing any surface that bounds cell 1 will
+be banked excluding any surface that does not use a 'transmission' or 'vacuum'
+boundary condition.
+
+.. note:: Surfaces with boundary conditions that are not "transmission" or "vacuum"
+          are not eligible to store any particles when using ``cell``, ``cellfrom``
+          or ``cellto`` attributes. It is recommended to use surface IDs instead.
+
+Surface IDs can be used in combination with a cell ID::
+
+  settings.surf_source_write = {
+      'cell': 1,
+      'surfaces_ids': [1, 2, 3],
+      'max_particles': 10000
+  }
+
+In that case, only particles that are crossing the declared surfaces coming from
+cell 1 or going to cell 1 will be banked. To account specifically for particles
+leaving or entering a given cell, ``cellfrom`` and ``cellto`` are also available
+to respectively account for particles coming from a cell::
+
+  settings.surf_source_write = {
+      'cellfrom': 1,
+      'max_particles': 10000
+  }
+
+or particles going to a cell::
+
+  settings.surf_source_write = {
+      'cellto': 1,
+      'max_particles': 10000
+  }
+
+.. note:: The ``cell``, ``cellfrom`` and ``cellto`` attributes cannot be
+          used simultaneously.
 
 .. _compiled_source:
 
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index e423880f8..6a2e67049 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -74,7 +74,7 @@ public:
   void pht_secondary_particles();
 
   //! Cross a surface and handle boundary conditions
-  void cross_surface();
+  void cross_surface(const Surface& surf);
 
   //! Cross a vacuum boundary condition.
   //
@@ -127,6 +127,8 @@ std::string particle_type_to_str(ParticleType type);
 
 ParticleType str_to_particle_type(std::string str);
 
+void add_surf_source_to_bank(Particle& p, const Surface& surf);
+
 } // namespace openmc
 
 #endif // OPENMC_PARTICLE_H
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 5c765e2e6..164148cce 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -211,9 +211,15 @@ public:
   // resets all coordinate levels for the particle
   void clear()
   {
-    for (auto& level : coord_)
+    for (auto& level : coord_) {
       level.reset();
+    }
     n_coord_ = 1;
+
+    for (auto& cell : cell_last_) {
+      cell = C_NONE;
+    }
+    n_coord_last_ = 1;
   }
 
   // Initialize all internal state from position and direction
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index df4a54a6f..a95f1ced9 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -16,6 +16,14 @@
 
 namespace openmc {
 
+// Type of surface source write
+enum class SSWCellType {
+  None,
+  Both,
+  From,
+  To,
+};
+
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -125,6 +133,10 @@ extern int
   max_history_splits; //!< maximum number of particle splits for weight windows
 extern int64_t max_surface_particles; //!< maximum number of particles to be
                                       //!< banked on surfaces per process
+extern int64_t ssw_cell_id;           //!< Cell id for the surface source
+                                      //!< write setting
+extern SSWCellType ssw_cell_type;     //!< Type of option for the cell
+                                      //!< argument of surface source write
 extern TemperatureMethod
   temperature_method; //!< method for choosing temperatures
 extern double
diff --git a/openmc/settings.py b/openmc/settings.py
index 43f763a19..80d766ab1 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -201,6 +201,15 @@ class Settings:
         :max_particles: Maximum number of particles to be banked on surfaces per
                    process (int)
         :mcpl: Output in the form of an MCPL-file (bool)
+        :cell: Cell ID used to determine if particles crossing identified
+               surfaces are to be banked. Particles coming from or going to this
+               declared cell will be banked (int)
+        :cellfrom: Cell ID used to determine if particles crossing identified
+                   surfaces are to be banked. Particles coming from this
+                   declared cell will be banked (int)
+        :cellto: Cell ID used to determine if particles crossing identified
+                 surfaces are to be banked. Particles going to this declared
+                 cell will be banked (int)
     survival_biasing : bool
         Indicate whether survival biasing is to be used
     tabular_legendre : dict
@@ -693,23 +702,30 @@ class Settings:
 
     @surf_source_write.setter
     def surf_source_write(self, surf_source_write: dict):
-        cv.check_type('surface source writing options', surf_source_write, Mapping)
+        cv.check_type("surface source writing options", surf_source_write, Mapping)
         for key, value in surf_source_write.items():
-            cv.check_value('surface source writing key', key,
-                           ('surface_ids', 'max_particles', 'mcpl'))
-            if key == 'surface_ids':
-                cv.check_type('surface ids for source banking', value,
-                              Iterable, Integral)
+            cv.check_value(
+                "surface source writing key",
+                key,
+                ("surface_ids", "max_particles", "mcpl", "cell", "cellfrom", "cellto"),
+            )
+            if key == "surface_ids":
+                cv.check_type(
+                    "surface ids for source banking", value, Iterable, Integral
+                )
                 for surf_id in value:
-                    cv.check_greater_than('surface id for source banking',
-                                          surf_id, 0)
-            elif key == 'max_particles':
-                cv.check_type('maximum particle banks on surfaces per process',
-                              value, Integral)
-                cv.check_greater_than('maximum particle banks on surfaces per process',
-                                      value, 0)
-            elif key == 'mcpl':
-                cv.check_type('write to an MCPL-format file', value, bool)
+                    cv.check_greater_than("surface id for source banking", surf_id, 0)
+            elif key == "mcpl":
+                cv.check_type("write to an MCPL-format file", value, bool)
+            elif key in ("max_particles", "cell", "cellfrom", "cellto"):
+                name = {
+                    "max_particles": "maximum particle banks on surfaces per process",
+                    "cell": "Cell ID for source banking (from or to)",
+                    "cellfrom": "Cell ID for source banking (from only)",
+                    "cellto": "Cell ID for source banking (to only)",
+                }[key]
+                cv.check_type(name, value, Integral)
+                cv.check_greater_than(name, value, 0)
 
         self._surf_source_write = surf_source_write
 
@@ -1207,16 +1223,18 @@ class Settings:
     def _create_surf_source_write_subelement(self, root):
         if self._surf_source_write:
             element = ET.SubElement(root, "surf_source_write")
-            if 'surface_ids' in self._surf_source_write:
+            if "surface_ids" in self._surf_source_write:
                 subelement = ET.SubElement(element, "surface_ids")
-                subelement.text = ' '.join(
-                    str(x) for x in self._surf_source_write['surface_ids'])
-            if 'max_particles' in self._surf_source_write:
-                subelement = ET.SubElement(element, "max_particles")
-                subelement.text = str(self._surf_source_write['max_particles'])
-            if 'mcpl' in self._surf_source_write:
+                subelement.text = " ".join(
+                    str(x) for x in self._surf_source_write["surface_ids"]
+                )
+            if "mcpl" in self._surf_source_write:
                 subelement = ET.SubElement(element, "mcpl")
-                subelement.text = str(self._surf_source_write['mcpl']).lower()
+                subelement.text = str(self._surf_source_write["mcpl"]).lower()
+            for key in ("max_particles", "cell", "cellfrom", "cellto"):
+                if key in self._surf_source_write:
+                    subelement = ET.SubElement(element, key)
+                    subelement.text = str(self._surf_source_write[key])
 
     def _create_confidence_intervals(self, root):
         if self._confidence_intervals is not None:
@@ -1381,9 +1399,9 @@ class Settings:
             elem.text = str(self._create_fission_neutrons).lower()
 
     def _create_create_delayed_neutrons_subelement(self, root):
-       if self._create_delayed_neutrons is not None:
-           elem = ET.SubElement(root, "create_delayed_neutrons")
-           elem.text = str(self._create_delayed_neutrons).lower()
+        if self._create_delayed_neutrons is not None:
+            elem = ET.SubElement(root, "create_delayed_neutrons")
+            elem.text = str(self._create_delayed_neutrons).lower()
 
     def _create_delayed_photon_scaling_subelement(self, root):
         if self._delayed_photon_scaling is not None:
@@ -1610,17 +1628,18 @@ class Settings:
 
     def _surf_source_write_from_xml_element(self, root):
         elem = root.find('surf_source_write')
-        if elem is not None:
-            for key in ('surface_ids', 'max_particles','mcpl'):
-                value = get_text(elem, key)
-                if value is not None:
-                    if key == 'surface_ids':
-                        value = [int(x) for x in value.split()]
-                    elif key in ('max_particles'):
-                        value = int(value)
-                    elif key == 'mcpl':
-                        value = value in ('true', '1')
-                    self.surf_source_write[key] = value
+        if elem is None:
+            return
+        for key in ('surface_ids', 'max_particles', 'mcpl', 'cell', 'cellto', 'cellfrom'):
+            value = get_text(elem, key)
+            if value is not None:
+                if key == 'surface_ids':
+                    value = [int(x) for x in value.split()]
+                elif key == 'mcpl':
+                    value = value in ('true', '1')
+                elif key in ('max_particles', 'cell', 'cellfrom', 'cellto'):
+                    value = int(value)
+                self.surf_source_write[key] = value
 
     def _confidence_intervals_from_xml_element(self, root):
         text = get_text(root, 'confidence_intervals')
diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index 2f2502f6e..a29a2589f 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -271,8 +271,10 @@ void DAGUniverse::init_geometry()
                                   : dagmc_instance_->id_by_index(2, i + 1);
 
     // set surface source attribute if needed
-    if (contains(settings::source_write_surf_id, s->id_))
+    if (contains(settings::source_write_surf_id, s->id_) ||
+        settings::source_write_surf_id.empty()) {
       s->surf_source_ = true;
+    }
 
     // set BCs
     std::string bc_value =
diff --git a/src/particle.cpp b/src/particle.cpp
index a91113c61..0c6a33b78 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -168,6 +168,12 @@ void Particle::event_calculate_xs()
     // Set birth cell attribute
     if (cell_born() == C_NONE)
       cell_born() = lowest_coord().cell;
+
+    // Initialize last cells from current cell
+    for (int j = 0; j < n_coord(); ++j) {
+      cell_last(j) = coord(j).cell;
+    }
+    n_coord_last() = n_coord();
   }
 
   // Write particle track.
@@ -264,16 +270,16 @@ void Particle::event_advance()
 
 void Particle::event_cross_surface()
 {
-  // Set surface that particle is on and adjust coordinate levels
-  surface() = boundary().surface_index;
-  n_coord() = boundary().coord_level;
-
   // Saving previous cell data
   for (int j = 0; j < n_coord(); ++j) {
     cell_last(j) = coord(j).cell;
   }
   n_coord_last() = n_coord();
 
+  // Set surface that particle is on and adjust coordinate levels
+  surface() = boundary().surface_index;
+  n_coord() = boundary().coord_level;
+
   if (boundary().lattice_translation[0] != 0 ||
       boundary().lattice_translation[1] != 0 ||
       boundary().lattice_translation[2] != 0) {
@@ -284,7 +290,17 @@ void Particle::event_cross_surface()
     event() = TallyEvent::LATTICE;
   } else {
     // Particle crosses surface
-    cross_surface();
+    // TODO: off-by-one
+    const auto& surf {model::surfaces[std::abs(surface()) - 1].get()};
+    // If BC, add particle to surface source before crossing surface
+    if (surf->surf_source_ && surf->bc_) {
+      add_surf_source_to_bank(*this, *surf);
+    }
+    cross_surface(*surf);
+    // If no BC, add particle to surface source after crossing surface
+    if (surf->surf_source_ && !surf->bc_) {
+      add_surf_source_to_bank(*this, *surf);
+    }
     if (settings::weight_window_checkpoint_surface) {
       apply_weight_windows(*this);
     }
@@ -418,6 +434,12 @@ void Particle::event_revive_from_secondary()
         // Set birth cell attribute
         if (cell_born() == C_NONE)
           cell_born() = lowest_coord().cell;
+
+        // Initialize last cells from current cell
+        for (int j = 0; j < n_coord(); ++j) {
+          cell_last(j) = coord(j).cell;
+        }
+        n_coord_last() = n_coord();
       }
       pht_secondary_particles();
     }
@@ -502,40 +524,22 @@ void Particle::pht_secondary_particles()
   }
 }
 
-void Particle::cross_surface()
+void Particle::cross_surface(const Surface& surf)
 {
-  int i_surface = std::abs(surface());
-  // TODO: off-by-one
-  const auto& surf {model::surfaces[i_surface - 1].get()};
-  if (settings::verbosity >= 10 || trace()) {
-    write_message(1, "    Crossing surface {}", surf->id_);
-  }
 
-  if (surf->surf_source_ && simulation::current_batch > settings::n_inactive &&
-      !simulation::surf_source_bank.full()) {
-    SourceSite site;
-    site.r = r();
-    site.u = u();
-    site.E = E();
-    site.time = time();
-    site.wgt = wgt();
-    site.delayed_group = delayed_group();
-    site.surf_id = surf->id_;
-    site.particle = type();
-    site.parent_id = id();
-    site.progeny_id = n_progeny();
-    int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
+  if (settings::verbosity >= 10 || trace()) {
+    write_message(1, "    Crossing surface {}", surf.id_);
   }
 
 // if we're crossing a CSG surface, make sure the DAG history is reset
 #ifdef DAGMC
-  if (surf->geom_type_ == GeometryType::CSG)
+  if (surf.geom_type_ == GeometryType::CSG)
     history().reset();
 #endif
 
   // Handle any applicable boundary conditions.
-  if (surf->bc_ && settings::run_mode != RunMode::PLOTTING) {
-    surf->bc_->handle_particle(*this, *surf);
+  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING) {
+    surf.bc_->handle_particle(*this, surf);
     return;
   }
 
@@ -544,10 +548,10 @@ void Particle::cross_surface()
 
 #ifdef DAGMC
   // in DAGMC, we know what the next cell should be
-  if (surf->geom_type_ == GeometryType::DAG) {
-    int32_t i_cell =
-      next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) -
-      1;
+  if (surf.geom_type_ == GeometryType::DAG) {
+    int32_t i_cell = next_cell(std::abs(surface()), cell_last(n_coord() - 1),
+                       lowest_coord().universe) -
+                     1;
     // save material and temp
     material_last() = material();
     sqrtkT_last() = sqrtkT();
@@ -561,8 +565,9 @@ void Particle::cross_surface()
 #endif
 
   bool verbose = settings::verbosity >= 10 || trace();
-  if (neighbor_list_find_cell(*this, verbose))
+  if (neighbor_list_find_cell(*this, verbose)) {
     return;
+  }
 
   // ==========================================================================
   // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
@@ -586,7 +591,7 @@ void Particle::cross_surface()
 
     if (!exhaustive_find_cell(*this, verbose)) {
       mark_as_lost("After particle " + std::to_string(id()) +
-                   " crossed surface " + std::to_string(surf->id_) +
+                   " crossed surface " + std::to_string(surf.id_) +
                    " it could not be located in any cell and it did not leak.");
       return;
     }
@@ -650,7 +655,7 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
   u() = new_u;
 
   // Reassign particle's cell and surface
-  coord(0).cell = cell_last(n_coord_last() - 1);
+  coord(0).cell = cell_last(0);
   surface() = -surface();
 
   // If a reflective surface is coincident with a lattice or universe
@@ -873,4 +878,90 @@ ParticleType str_to_particle_type(std::string str)
   }
 }
 
+void add_surf_source_to_bank(Particle& p, const Surface& surf)
+{
+  if (simulation::current_batch <= settings::n_inactive ||
+      simulation::surf_source_bank.full()) {
+    return;
+  }
+
+  // If a cell/cellfrom/cellto parameter is defined
+  if (settings::ssw_cell_id != C_NONE) {
+
+    // Retrieve cell index and storage type
+    int cell_idx = model::cell_map[settings::ssw_cell_id];
+
+    if (surf.bc_) {
+      // Leave if cellto with vacuum boundary condition
+      if (surf.bc_->type() == "vacuum" &&
+          settings::ssw_cell_type == SSWCellType::To) {
+        return;
+      }
+
+      // Leave if other boundary condition than vacuum
+      if (surf.bc_->type() != "vacuum") {
+        return;
+      }
+    }
+
+    // Check if the cell of interest has been exited
+    bool exited = false;
+    for (int i = 0; i < p.n_coord_last(); ++i) {
+      if (p.cell_last(i) == cell_idx) {
+        exited = true;
+      }
+    }
+
+    // Check if the cell of interest has been entered
+    bool entered = false;
+    for (int i = 0; i < p.n_coord(); ++i) {
+      if (p.coord(i).cell == cell_idx) {
+        entered = true;
+      }
+    }
+
+    // Vacuum boundary conditions: return if cell is not exited
+    if (surf.bc_) {
+      if (surf.bc_->type() == "vacuum" && !exited) {
+        return;
+      }
+    } else {
+
+      // If we both enter and exit the cell of interest
+      if (entered && exited) {
+        return;
+      }
+
+      // If we did not enter nor exit the cell of interest
+      if (!entered && !exited) {
+        return;
+      }
+
+      // If cellfrom and the cell before crossing is not the cell of
+      // interest
+      if (settings::ssw_cell_type == SSWCellType::From && !exited) {
+        return;
+      }
+
+      // If cellto and the cell after crossing is not the cell of interest
+      if (settings::ssw_cell_type == SSWCellType::To && !entered) {
+        return;
+      }
+    }
+  }
+
+  SourceSite site;
+  site.r = p.r();
+  site.u = p.u();
+  site.E = p.E();
+  site.time = p.time();
+  site.wgt = p.wgt();
+  site.delayed_group = p.delayed_group();
+  site.surf_id = surf.id_;
+  site.particle = p.type();
+  site.parent_id = p.id();
+  site.progeny_id = p.n_progeny();
+  int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
+}
+
 } // namespace openmc
diff --git a/src/settings.cpp b/src/settings.cpp
index f5cd1b7e9..6c3226b04 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -119,6 +119,8 @@ std::unordered_set<int> sourcepoint_batch;
 std::unordered_set<int> statepoint_batch;
 std::unordered_set<int> source_write_surf_id;
 int64_t max_surface_particles;
+int64_t ssw_cell_id {C_NONE};
+SSWCellType ssw_cell_type {SSWCellType::None};
 TemperatureMethod temperature_method {TemperatureMethod::NEAREST};
 double temperature_tolerance {10.0};
 double temperature_default {293.6};
@@ -747,7 +749,9 @@ void read_settings_xml(pugi::xml_node root)
     // Get surface source write node
     xml_node node_ssw = root.child("surf_source_write");
 
-    // Determine surface ids at which crossing particles are to be banked
+    // Determine surface ids at which crossing particles are to be banked.
+    // If no surfaces are specified, all surfaces in the model will be used
+    // to bank source points.
     if (check_for_node(node_ssw, "surface_ids")) {
       auto temp = get_node_array<int>(node_ssw, "surface_ids");
       for (const auto& b : temp) {
@@ -759,7 +763,12 @@ void read_settings_xml(pugi::xml_node root)
     if (check_for_node(node_ssw, "max_particles")) {
       max_surface_particles =
         std::stoll(get_node_value(node_ssw, "max_particles"));
+    } else {
+      fatal_error("A maximum number of particles needs to be specified "
+                  "using the 'max_particles' parameter to store surface "
+                  "source points.");
     }
+
     if (check_for_node(node_ssw, "mcpl")) {
       surf_mcpl_write = get_node_value_bool(node_ssw, "mcpl");
 
@@ -769,6 +778,27 @@ void read_settings_xml(pugi::xml_node root)
                     "surface source files.");
       }
     }
+    // Get cell information
+    if (check_for_node(node_ssw, "cell")) {
+      ssw_cell_id = std::stoll(get_node_value(node_ssw, "cell"));
+      ssw_cell_type = SSWCellType::Both;
+    }
+    if (check_for_node(node_ssw, "cellfrom")) {
+      if (ssw_cell_id != C_NONE) {
+        fatal_error(
+          "'cell', 'cellfrom' and 'cellto' cannot be used at the same time.");
+      }
+      ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellfrom"));
+      ssw_cell_type = SSWCellType::From;
+    }
+    if (check_for_node(node_ssw, "cellto")) {
+      if (ssw_cell_id != C_NONE) {
+        fatal_error(
+          "'cell', 'cellfrom' and 'cellto' cannot be used at the same time.");
+      }
+      ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellto"));
+      ssw_cell_type = SSWCellType::To;
+    }
   }
 
   // If source is not separate and is to be written out in the statepoint file,
diff --git a/src/surface.cpp b/src/surface.cpp
index 12fef070e..50ef2a128 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -63,7 +63,8 @@ Surface::Surface(pugi::xml_node surf_node)
 {
   if (check_for_node(surf_node, "id")) {
     id_ = std::stoi(get_node_value(surf_node, "id"));
-    if (contains(settings::source_write_surf_id, id_)) {
+    if (contains(settings::source_write_surf_id, id_) ||
+        settings::source_write_surf_id.empty()) {
       surf_source_ = true;
     }
   } else {
diff --git a/tests/regression_tests/filter_cellfrom/__init__.py b/tests/regression_tests/filter_cellfrom/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/regression_tests/filter_cellfrom/inputs_true.dat b/tests/regression_tests/filter_cellfrom/inputs_true.dat
new file mode 100644
index 000000000..8a37e8286
--- /dev/null
+++ b/tests/regression_tests/filter_cellfrom/inputs_true.dat
@@ -0,0 +1,151 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="10.97"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+    <material id="3">
+      <density units="g/cm3" value="0.001225"/>
+      <nuclide ao="0.1999242" name="O16"/>
+      <nuclide ao="7.58e-05" name="O17"/>
+      <nuclide ao="0.7970696" name="N14"/>
+      <nuclide ao="0.0029304" name="N15"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="3" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-2.05" id="4" type="x-plane"/>
+    <surface coeffs="2.05" id="5" type="x-plane"/>
+    <surface coeffs="-2.05" id="6" type="y-plane"/>
+    <surface coeffs="2.05" id="7" type="y-plane"/>
+    <surface coeffs="-2.05" id="8" type="z-plane"/>
+    <surface coeffs="2.05" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-6.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="6.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-6.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="6.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-6.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="6.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>2000</particles>
+    <batches>15</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <seed>1</seed>
+  </settings>
+  <tallies>
+    <filter id="5" type="cellfrom">
+      <bins>1</bins>
+    </filter>
+    <filter id="1" type="cell">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="cell">
+      <bins>2</bins>
+    </filter>
+    <filter id="3" type="cell">
+      <bins>3</bins>
+    </filter>
+    <filter id="4" type="cell">
+      <bins>4</bins>
+    </filter>
+    <filter id="6" type="cellfrom">
+      <bins>2</bins>
+    </filter>
+    <filter id="7" type="cellfrom">
+      <bins>3</bins>
+    </filter>
+    <filter id="8" type="cellfrom">
+      <bins>4</bins>
+    </filter>
+    <tally id="1" name="total from 1 in 1">
+      <filters>5 1</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="2" name="total from 1 in 2">
+      <filters>5 2</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="3" name="total from 1 in 3">
+      <filters>5 3</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="4" name="total from 1 in 4">
+      <filters>5 4</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="5" name="total from 2 in 1">
+      <filters>6 1</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="6" name="total from 2 in 2">
+      <filters>6 2</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="7" name="total from 2 in 3">
+      <filters>6 3</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="8" name="total from 2 in 4">
+      <filters>6 4</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="9" name="total from 3 in 1">
+      <filters>7 1</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="10" name="total from 3 in 2">
+      <filters>7 2</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="11" name="total from 3 in 3">
+      <filters>7 3</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="12" name="total from 3 in 4">
+      <filters>7 4</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="13" name="total from 4 in 1">
+      <filters>8 1</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="14" name="total from 4 in 2">
+      <filters>8 2</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="15" name="total from 4 in 3">
+      <filters>8 3</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="16" name="total from 4 in 4">
+      <filters>8 4</filters>
+      <scores>total</scores>
+    </tally>
+    <tally id="17" name="total">
+      <scores>total</scores>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/filter_cellfrom/results_true.dat b/tests/regression_tests/filter_cellfrom/results_true.dat
new file mode 100644
index 000000000..5dd43e1f3
--- /dev/null
+++ b/tests/regression_tests/filter_cellfrom/results_true.dat
@@ -0,0 +1,53 @@
+k-combined:
+9.640806E-02 2.655206E-03
+tally 1:
+6.025999E+00
+3.635317E+00
+tally 2:
+4.523606E-04
+2.051522E-08
+tally 3:
+6.026452E+00
+3.635862E+00
+tally 4:
+0.000000E+00
+0.000000E+00
+tally 5:
+1.530903E+00
+2.357048E-01
+tally 6:
+4.992646E-05
+2.600391E-10
+tally 7:
+1.530953E+00
+2.357201E-01
+tally 8:
+1.889115E+01
+3.574727E+01
+tally 9:
+7.556902E+00
+5.717680E+00
+tally 10:
+5.022871E-04
+2.531352E-08
+tally 11:
+7.557405E+00
+5.718440E+00
+tally 12:
+1.889115E+01
+3.574727E+01
+tally 13:
+0.000000E+00
+0.000000E+00
+tally 14:
+2.663407E-04
+7.161725E-09
+tally 15:
+2.663407E-04
+7.161725E-09
+tally 16:
+8.025710E+01
+6.459305E+02
+tally 17:
+1.067059E+02
+1.140939E+03
diff --git a/tests/regression_tests/filter_cellfrom/test.py b/tests/regression_tests/filter_cellfrom/test.py
new file mode 100644
index 000000000..8fb7509cb
--- /dev/null
+++ b/tests/regression_tests/filter_cellfrom/test.py
@@ -0,0 +1,315 @@
+"""This test ensures that the CellFromFilter works correctly even if the level of
+coordinates (number of encapsulated universes) is different in the cell from
+where the particle originates compared to the cell where the particle is going.
+
+A matrix of reaction rates based on where the particle is coming from and
+where it goes to is calculated and compared to the total reaction rate of the problem.
+The components of this matrix are also compared to other components using symmetric
+properties.
+
+TODO:
+
+- Test with a lattice,
+- Test with mesh,
+- Test with reflective boundary conditions,
+- Test with periodic boundary conditions.
+
+"""
+
+from numpy.testing import assert_allclose, assert_equal
+import numpy as np
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+RTOL = 1.0e-7
+ATOL = 0.0
+
+
+@pytest.fixture
+def model():
+    """Cylindrical core contained in a first box which is contained in a larger box.
+    A lower universe is used to describe the interior of the first box which
+    contains the core and its surrounding space."""
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    fuel = openmc.Material()
+    fuel.add_nuclide("U234", 0.0004524)
+    fuel.add_nuclide("U235", 0.0506068)
+    fuel.add_nuclide("U238", 0.9487090)
+    fuel.add_nuclide("U236", 0.0002318)
+    fuel.add_nuclide("O16", 2.0)
+    fuel.set_density("g/cm3", 10.97)
+
+    water = openmc.Material()
+    water.add_nuclide("H1", 2.0)
+    water.add_nuclide("O16", 1.0)
+    water.set_density("g/cm3", 1.0)
+
+    air = openmc.Material()
+    air.add_element("O", 0.2)
+    air.add_element("N", 0.8)
+    air.set_density("g/cm3", 0.001225)
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    # -----------------------------------------------------------------------------
+    # Cylindrical core
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    core_radius = 2.0
+    core_height = 4.0
+
+    # Surfaces
+    core_cylinder = openmc.ZCylinder(r=core_radius)
+    core_lower_plane = openmc.ZPlane(z0=-core_height / 2.0)
+    core_upper_plane = openmc.ZPlane(z0=core_height / 2.0)
+
+    # Region
+    core_region = -core_cylinder & +core_lower_plane & -core_upper_plane
+
+    # Cells
+    core = openmc.Cell(fill=fuel, region=core_region)
+    outside_core_region = +core_cylinder | -core_lower_plane | +core_upper_plane
+    outside_core = openmc.Cell(fill=air, region=outside_core_region)
+
+    # Universe
+    inside_box1_universe = openmc.Universe(cells=[core, outside_core])
+
+    # -----------------------------------------------------------------------------
+    # Box 1
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box1_size = 4.1
+
+    # Surfaces
+    box1_rpp = openmc.model.RectangularParallelepiped(
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+    )
+
+    # Cell
+    box1 = openmc.Cell(fill=inside_box1_universe, region=-box1_rpp)
+
+    # -----------------------------------------------------------------------------
+    # Box 2
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box2_size = 12
+
+    # Surfaces
+    box2_rpp = openmc.model.RectangularParallelepiped(
+        -box2_size / 2.0, box2_size / 2.0,
+        -box2_size / 2.0, box2_size / 2.0,
+        -box2_size / 2.0, box2_size / 2.0,
+        boundary_type="vacuum"
+    )
+
+    # Cell
+    box2 = openmc.Cell(fill=water, region=-box2_rpp & +box1_rpp)
+
+    # Register geometry
+    model.geometry = openmc.Geometry([box1, box2])
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 2000
+    model.settings.batches = 15
+    model.settings.inactive = 5
+    model.settings.seed = 1
+
+    bounds = [
+        -core_radius,
+        -core_radius,
+        -core_height / 2.0,
+        core_radius,
+        core_radius,
+        core_height / 2.0,
+    ]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    # =============================================================================
+    # Tallies
+    # =============================================================================
+
+    in_core_filter = openmc.CellFilter([core])
+    in_outside_core_filter = openmc.CellFilter([outside_core])
+    in_box1_filter = openmc.CellFilter([box1])
+    in_box2_filter = openmc.CellFilter([box2])
+
+    from_core_filter = openmc.CellFromFilter([core])
+    from_outside_core_filter = openmc.CellFromFilter([outside_core])
+    from_box1_filter = openmc.CellFromFilter([box1])
+    from_box2_filter = openmc.CellFromFilter([box2])
+
+    t1_1 = openmc.Tally(name="total from 1 in 1")
+    t1_1.filters = [from_core_filter, in_core_filter]
+    t1_1.scores = ["total"]
+
+    t1_2 = openmc.Tally(name="total from 1 in 2")
+    t1_2.filters = [from_core_filter, in_outside_core_filter]
+    t1_2.scores = ["total"]
+
+    t1_3 = openmc.Tally(name="total from 1 in 3")
+    t1_3.filters = [from_core_filter, in_box1_filter]
+    t1_3.scores = ["total"]
+
+    t1_4 = openmc.Tally(name="total from 1 in 4")
+    t1_4.filters = [from_core_filter, in_box2_filter]
+    t1_4.scores = ["total"]
+
+    t2_1 = openmc.Tally(name="total from 2 in 1")
+    t2_1.filters = [from_outside_core_filter, in_core_filter]
+    t2_1.scores = ["total"]
+
+    t2_2 = openmc.Tally(name="total from 2 in 2")
+    t2_2.filters = [from_outside_core_filter, in_outside_core_filter]
+    t2_2.scores = ["total"]
+
+    t2_3 = openmc.Tally(name="total from 2 in 3")
+    t2_3.filters = [from_outside_core_filter, in_box1_filter]
+    t2_3.scores = ["total"]
+
+    t2_4 = openmc.Tally(name="total from 2 in 4")
+    t2_4.filters = [from_outside_core_filter, in_box2_filter]
+    t2_4.scores = ["total"]
+
+    t3_1 = openmc.Tally(name="total from 3 in 1")
+    t3_1.filters = [from_box1_filter, in_core_filter]
+    t3_1.scores = ["total"]
+
+    t3_2 = openmc.Tally(name="total from 3 in 2")
+    t3_2.filters = [from_box1_filter, in_outside_core_filter]
+    t3_2.scores = ["total"]
+
+    t3_3 = openmc.Tally(name="total from 3 in 3")
+    t3_3.filters = [from_box1_filter, in_box1_filter]
+    t3_3.scores = ["total"]
+
+    t3_4 = openmc.Tally(name="total from 3 in 4")
+    t3_4.filters = [from_box1_filter, in_box2_filter]
+    t3_4.scores = ["total"]
+
+    t4_1 = openmc.Tally(name="total from 4 in 1")
+    t4_1.filters = [from_box2_filter, in_core_filter]
+    t4_1.scores = ["total"]
+
+    t4_2 = openmc.Tally(name="total from 4 in 2")
+    t4_2.filters = [from_box2_filter, in_outside_core_filter]
+    t4_2.scores = ["total"]
+
+    t4_3 = openmc.Tally(name="total from 4 in 3")
+    t4_3.filters = [from_box2_filter, in_box1_filter]
+    t4_3.scores = ["total"]
+
+    t4_4 = openmc.Tally(name="total from 4 in 4")
+    t4_4.filters = [from_box2_filter, in_box2_filter]
+    t4_4.scores = ["total"]
+
+    tglobal = openmc.Tally(name="total")
+    tglobal.scores = ["total"]
+
+    model.tallies += [
+        t1_1,
+        t1_2,
+        t1_3,
+        t1_4,
+        t2_1,
+        t2_2,
+        t2_3,
+        t2_4,
+        t3_1,
+        t3_2,
+        t3_3,
+        t3_4,
+        t4_1,
+        t4_2,
+        t4_3,
+        t4_4,
+        tglobal,
+    ]
+    return model
+
+
+class CellFromFilterTest(PyAPITestHarness):
+
+    def _compare_results(self):
+        """Additional unit tests on the tally results to check
+        consistency of CellFromFilter."""
+        with openmc.StatePoint(self.statepoint_name) as sp:
+
+            t1_1 = sp.get_tally(name="total from 1 in 1").mean
+            t1_2 = sp.get_tally(name="total from 1 in 2").mean
+            t1_3 = sp.get_tally(name="total from 1 in 3").mean
+            t1_4 = sp.get_tally(name="total from 1 in 4").mean
+
+            t2_1 = sp.get_tally(name="total from 2 in 1").mean
+            t2_2 = sp.get_tally(name="total from 2 in 2").mean
+            t2_3 = sp.get_tally(name="total from 2 in 3").mean
+            t2_4 = sp.get_tally(name="total from 2 in 4").mean
+
+            t3_1 = sp.get_tally(name="total from 3 in 1").mean
+            t3_2 = sp.get_tally(name="total from 3 in 2").mean
+            t3_3 = sp.get_tally(name="total from 3 in 3").mean
+            t3_4 = sp.get_tally(name="total from 3 in 4").mean
+
+            t4_1 = sp.get_tally(name="total from 4 in 1").mean
+            t4_2 = sp.get_tally(name="total from 4 in 2").mean
+            t4_3 = sp.get_tally(name="total from 4 in 3").mean
+            t4_4 = sp.get_tally(name="total from 4 in 4").mean
+
+            tglobal = sp.get_tally(name="total").mean
+
+            # From 1 and 2 is equivalent to from 3
+            assert_allclose(t1_1 + t2_1, t3_1, rtol=RTOL, atol=ATOL)
+            assert_allclose(t1_2 + t2_2, t3_2, rtol=RTOL, atol=ATOL)
+            assert_allclose(t1_3 + t2_3, t3_3, rtol=RTOL, atol=ATOL)
+            assert_allclose(t1_4 + t2_4, t3_4, rtol=RTOL, atol=ATOL)
+
+            # In 1 and 2 equivalent to in 3
+            assert_allclose(t1_1 + t1_2, t1_3, rtol=RTOL, atol=ATOL)
+            assert_allclose(t2_1 + t2_2, t2_3, rtol=RTOL, atol=ATOL)
+            assert_allclose(t3_1 + t3_2, t3_3, rtol=RTOL, atol=ATOL)
+            assert_allclose(t4_1 + t4_2, t4_3, rtol=RTOL, atol=ATOL)
+
+            # Comparison to global from 3
+            assert_allclose(t3_3 + t3_4 + t4_3 + t4_4, tglobal, rtol=RTOL, atol=ATOL)
+
+            # Comparison to global from 1 and 2
+            t_from_1_wo_3 = t1_1 + t1_2 + t1_4
+            t_from_2_wo_3 = t2_1 + t2_2 + t2_4
+            t_from_4_wo_3 = t4_1 + t4_2 + t4_4
+            assert_allclose(
+                t_from_1_wo_3 + t_from_2_wo_3 + t_from_4_wo_3,
+                tglobal,
+                rtol=RTOL,
+                atol=ATOL,
+            )
+
+            # 1 cannot contribute to 4 and 4 cannot contribute to 1 by symmetry
+            assert_equal(t1_4, np.zeros_like(t1_4))
+            assert_equal(t4_1, np.zeros_like(t4_1))
+
+        return super()._compare_results()
+
+
+def test_filter_cellfrom(model):
+    harness = CellFromFilterTest("statepoint.15.h5", model)
+    harness.main()
diff --git a/tests/regression_tests/surface_source/test.py b/tests/regression_tests/surface_source/test.py
index a25b40064..81bba0c94 100644
--- a/tests/regression_tests/surface_source/test.py
+++ b/tests/regression_tests/surface_source/test.py
@@ -118,15 +118,13 @@ class SurfaceSourceTestHarness(PyAPITestHarness):
 
 
 @pytest.mark.surf_source_op('write')
-def test_surface_source_write(model):
+def test_surface_source_write(model, monkeypatch):
     # Test result is based on 1 MPI process
-    np = config['mpi_np']
-    config['mpi_np'] = '1'
+    monkeypatch.setitem(config, "mpi_np", "1")
     harness = SurfaceSourceTestHarness('statepoint.10.h5',
                                        model,
                                        'inputs_true_write.dat')
     harness.main()
-    config['mpi_np'] = np
 
 
 @pytest.mark.surf_source_op('read')
diff --git a/tests/regression_tests/surface_source_write/__init__.py b/tests/regression_tests/surface_source_write/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/regression_tests/surface_source_write/_visualize.py b/tests/regression_tests/surface_source_write/_visualize.py
new file mode 100644
index 000000000..73340cae0
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/_visualize.py
@@ -0,0 +1,66 @@
+"""Helper script to visualize the surface_source.h5 files created with this test.
+"""
+
+import h5py
+import matplotlib.pyplot as plt
+
+
+if __name__ == "__main__":
+
+    # Select an option
+    # "show": 3D visualization using matplotlib
+    # "savefig": 2D representation using matplotlib and storing the fig under plot_2d.png
+    option = "show"
+    # option = "savefig"
+
+    # Select the case from its folder name
+    folder = "case-20"
+
+    # Reading the surface source file
+    with h5py.File(f"{folder}/surface_source_true.h5", "r") as fp:
+        source_bank = fp["source_bank"][()]
+        r_xs = source_bank['r']['x']
+        r_ys = source_bank['r']['y']
+        r_zs = source_bank['r']['z']
+
+    print("Size of the source bank: ", len(source_bank))
+
+    # Select data range to visualize
+    idx_1 = 0
+    idx_2 = -1
+
+    # Show 3D representation
+    if option == "show":
+
+        fig = plt.figure(figsize=(10, 10))
+        ax1 = fig.add_subplot(projection="3d", proj_type="ortho")
+        ax1.scatter(r_xs[idx_1:idx_2], r_ys[idx_1:idx_2], r_zs[idx_1:idx_2], marker=".")
+        ax1.view_init(0, 0)
+        ax1.xaxis.set_ticklabels([])
+        ax1.set_ylabel("y-axis [cm]")
+        ax1.set_zlabel("z-axis [cm]")
+        ax1.set_aspect("equal", "box")
+
+        plt.show()
+
+    # Save 2D representations
+    elif option == "savefig":
+
+        fig = plt.figure(figsize=(14, 5))
+        ax1 = fig.add_subplot(121, projection="3d", proj_type="ortho")
+        ax1.scatter(r_xs[idx_1:idx_2], r_ys[idx_1:idx_2], r_zs[idx_1:idx_2], marker=".")
+        ax1.view_init(0, 0)
+        ax1.xaxis.set_ticklabels([])
+        ax1.set_ylabel("y-axis [cm]")
+        ax1.set_zlabel("z-axis [cm]")
+        ax1.set_aspect("equal", "box")
+
+        ax2 = fig.add_subplot(122, projection="3d", proj_type="ortho")
+        ax2.scatter(r_xs[idx_1:idx_2], r_ys[idx_1:idx_2], r_zs[idx_1:idx_2], marker=".")
+        ax2.view_init(90, -90)
+        ax2.zaxis.set_ticklabels([])
+        ax2.set_xlabel("x-axis [cm]")
+        ax2.set_ylabel("y-axis [cm]")
+        ax2.set_aspect("equal", "box")
+
+        plt.savefig("plot_2d.png")
diff --git a/tests/regression_tests/surface_source_write/case-01/inputs_true.dat b/tests/regression_tests/surface_source_write/case-01/inputs_true.dat
new file mode 100644
index 000000000..ecca3a9f8
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-01/inputs_true.dat
@@ -0,0 +1,54 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-01/results_true.dat b/tests/regression_tests/surface_source_write/case-01/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-01/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-01/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-01/surface_source_true.h5
new file mode 100644
index 000000000..82dca4d03
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-01/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-02/inputs_true.dat b/tests/regression_tests/surface_source_write/case-02/inputs_true.dat
new file mode 100644
index 000000000..8be841ad6
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-02/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>8</surface_ids>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-02/results_true.dat b/tests/regression_tests/surface_source_write/case-02/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-02/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-02/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-02/surface_source_true.h5
new file mode 100644
index 000000000..872efa070
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-02/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-03/inputs_true.dat b/tests/regression_tests/surface_source_write/case-03/inputs_true.dat
new file mode 100644
index 000000000..38faa2a78
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-03/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-03/results_true.dat b/tests/regression_tests/surface_source_write/case-03/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-03/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-03/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-03/surface_source_true.h5
new file mode 100644
index 000000000..04a6ac9a3
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-03/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-04/inputs_true.dat b/tests/regression_tests/surface_source_write/case-04/inputs_true.dat
new file mode 100644
index 000000000..99a5a2899
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-04/inputs_true.dat
@@ -0,0 +1,56 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>2</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-04/results_true.dat b/tests/regression_tests/surface_source_write/case-04/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-04/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-04/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-04/surface_source_true.h5
new file mode 100644
index 000000000..75c309eb4
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-04/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-05/inputs_true.dat b/tests/regression_tests/surface_source_write/case-05/inputs_true.dat
new file mode 100644
index 000000000..49eedb7cf
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-05/inputs_true.dat
@@ -0,0 +1,56 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>3</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-05/results_true.dat b/tests/regression_tests/surface_source_write/case-05/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-05/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-05/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-05/surface_source_true.h5
new file mode 100644
index 000000000..1a2f54d76
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-05/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-06/inputs_true.dat b/tests/regression_tests/surface_source_write/case-06/inputs_true.dat
new file mode 100644
index 000000000..cdb9726d6
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-06/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cell>2</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-06/results_true.dat b/tests/regression_tests/surface_source_write/case-06/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-06/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-06/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-06/surface_source_true.h5
new file mode 100644
index 000000000..bef135e5a
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-06/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-07/inputs_true.dat b/tests/regression_tests/surface_source_write/case-07/inputs_true.dat
new file mode 100644
index 000000000..ff73358ca
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-07/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cell>3</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-07/results_true.dat b/tests/regression_tests/surface_source_write/case-07/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-07/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-07/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-07/surface_source_true.h5
new file mode 100644
index 000000000..f84a5c4c5
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-07/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-08/inputs_true.dat b/tests/regression_tests/surface_source_write/case-08/inputs_true.dat
new file mode 100644
index 000000000..8bda9795d
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-08/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cellfrom>2</cellfrom>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-08/results_true.dat b/tests/regression_tests/surface_source_write/case-08/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-08/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-08/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-08/surface_source_true.h5
new file mode 100644
index 000000000..1fd302a7c
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-08/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-09/inputs_true.dat b/tests/regression_tests/surface_source_write/case-09/inputs_true.dat
new file mode 100644
index 000000000..4e893fffa
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-09/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cellto>2</cellto>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-09/results_true.dat b/tests/regression_tests/surface_source_write/case-09/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-09/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-09/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-09/surface_source_true.h5
new file mode 100644
index 000000000..af832f060
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-09/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-10/inputs_true.dat b/tests/regression_tests/surface_source_write/case-10/inputs_true.dat
new file mode 100644
index 000000000..2d8f1ce64
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-10/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cellfrom>3</cellfrom>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-10/results_true.dat b/tests/regression_tests/surface_source_write/case-10/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-10/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-10/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-10/surface_source_true.h5
new file mode 100644
index 000000000..1512d8692
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-10/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-11/inputs_true.dat b/tests/regression_tests/surface_source_write/case-11/inputs_true.dat
new file mode 100644
index 000000000..307dec547
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-11/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cellto>3</cellto>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-11/results_true.dat b/tests/regression_tests/surface_source_write/case-11/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-11/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-11/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-11/surface_source_true.h5
new file mode 100644
index 000000000..0edb06c88
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-11/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-12/inputs_true.dat b/tests/regression_tests/surface_source_write/case-12/inputs_true.dat
new file mode 100644
index 000000000..087024695
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-12/inputs_true.dat
@@ -0,0 +1,48 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-12/results_true.dat b/tests/regression_tests/surface_source_write/case-12/results_true.dat
new file mode 100644
index 000000000..d793a7e42
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-12/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+4.752377E-02 6.773395E-03
diff --git a/tests/regression_tests/surface_source_write/case-12/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-12/surface_source_true.h5
new file mode 100644
index 000000000..ea81a8250
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-12/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-13/inputs_true.dat b/tests/regression_tests/surface_source_write/case-13/inputs_true.dat
new file mode 100644
index 000000000..a3421298b
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-13/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>3</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-13/results_true.dat b/tests/regression_tests/surface_source_write/case-13/results_true.dat
new file mode 100644
index 000000000..d793a7e42
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-13/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+4.752377E-02 6.773395E-03
diff --git a/tests/regression_tests/surface_source_write/case-13/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-13/surface_source_true.h5
new file mode 100644
index 000000000..08bd80985
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-13/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-14/inputs_true.dat b/tests/regression_tests/surface_source_write/case-14/inputs_true.dat
new file mode 100644
index 000000000..b4c141b5d
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-14/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cellfrom>3</cellfrom>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-14/results_true.dat b/tests/regression_tests/surface_source_write/case-14/results_true.dat
new file mode 100644
index 000000000..d793a7e42
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-14/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+4.752377E-02 6.773395E-03
diff --git a/tests/regression_tests/surface_source_write/case-14/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-14/surface_source_true.h5
new file mode 100644
index 000000000..50cddcc70
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-14/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-15/inputs_true.dat b/tests/regression_tests/surface_source_write/case-15/inputs_true.dat
new file mode 100644
index 000000000..9244d839d
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-15/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cellto>3</cellto>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-15/results_true.dat b/tests/regression_tests/surface_source_write/case-15/results_true.dat
new file mode 100644
index 000000000..d793a7e42
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-15/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+4.752377E-02 6.773395E-03
diff --git a/tests/regression_tests/surface_source_write/case-15/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-15/surface_source_true.h5
new file mode 100644
index 000000000..3e15017a1
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-15/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-16/inputs_true.dat b/tests/regression_tests/surface_source_write/case-16/inputs_true.dat
new file mode 100644
index 000000000..b3139928a
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-16/inputs_true.dat
@@ -0,0 +1,48 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-16/results_true.dat b/tests/regression_tests/surface_source_write/case-16/results_true.dat
new file mode 100644
index 000000000..a97be70ca
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-16/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+1.496403E+00 3.891283E-02
diff --git a/tests/regression_tests/surface_source_write/case-16/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-16/surface_source_true.h5
new file mode 100644
index 000000000..49f67e820
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-16/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-17/inputs_true.dat b/tests/regression_tests/surface_source_write/case-17/inputs_true.dat
new file mode 100644
index 000000000..8fcb01e45
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-17/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>3</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-17/results_true.dat b/tests/regression_tests/surface_source_write/case-17/results_true.dat
new file mode 100644
index 000000000..a97be70ca
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-17/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+1.496403E+00 3.891283E-02
diff --git a/tests/regression_tests/surface_source_write/case-17/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-17/surface_source_true.h5
new file mode 100644
index 000000000..35878eeca
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-17/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-18/inputs_true.dat b/tests/regression_tests/surface_source_write/case-18/inputs_true.dat
new file mode 100644
index 000000000..80d6253a1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-18/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cellfrom>3</cellfrom>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-18/results_true.dat b/tests/regression_tests/surface_source_write/case-18/results_true.dat
new file mode 100644
index 000000000..a97be70ca
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-18/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+1.496403E+00 3.891283E-02
diff --git a/tests/regression_tests/surface_source_write/case-18/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-18/surface_source_true.h5
new file mode 100644
index 000000000..36da9e7a7
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-18/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-19/inputs_true.dat b/tests/regression_tests/surface_source_write/case-19/inputs_true.dat
new file mode 100644
index 000000000..5355675fe
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-19/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cellto>3</cellto>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-19/results_true.dat b/tests/regression_tests/surface_source_write/case-19/results_true.dat
new file mode 100644
index 000000000..a97be70ca
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-19/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+1.496403E+00 3.891283E-02
diff --git a/tests/regression_tests/surface_source_write/case-19/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-19/surface_source_true.h5
new file mode 100644
index 000000000..f4261451d
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-19/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-20/inputs_true.dat b/tests/regression_tests/surface_source_write/case-20/inputs_true.dat
new file mode 100644
index 000000000..44f0d72c7
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-20/inputs_true.dat
@@ -0,0 +1,48 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="4 -5 6 -7 8 -9" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="periodic" coeffs="-3.0" id="4" type="z-plane"/>
+    <surface boundary="periodic" coeffs="3.0" id="5" type="z-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="6" type="x-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="7" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="8" type="y-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="9" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4</surface_ids>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-20/results_true.dat b/tests/regression_tests/surface_source_write/case-20/results_true.dat
new file mode 100644
index 000000000..7ccce7cc3
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-20/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+1.149925E+00 2.542255E-01
diff --git a/tests/regression_tests/surface_source_write/case-20/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-20/surface_source_true.h5
new file mode 100644
index 000000000..fa8a4e628
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-20/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-21/inputs_true.dat b/tests/regression_tests/surface_source_write/case-21/inputs_true.dat
new file mode 100644
index 000000000..893b8b5b9
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-21/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="4 -5 6 -7 8 -9" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="periodic" coeffs="-3.0" id="4" type="z-plane"/>
+    <surface boundary="periodic" coeffs="3.0" id="5" type="z-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="6" type="x-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="7" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-3.0" id="8" type="y-plane"/>
+    <surface boundary="reflective" coeffs="3.0" id="9" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>3</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-21/results_true.dat b/tests/regression_tests/surface_source_write/case-21/results_true.dat
new file mode 100644
index 000000000..7ccce7cc3
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-21/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+1.149925E+00 2.542255E-01
diff --git a/tests/regression_tests/surface_source_write/case-21/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-21/surface_source_true.h5
new file mode 100644
index 000000000..ed87ce849
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-21/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-a01/inputs_true.dat b/tests/regression_tests/surface_source_write/case-a01/inputs_true.dat
new file mode 100644
index 000000000..a65c9f770
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-a01/inputs_true.dat
@@ -0,0 +1,56 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>1 2 3</surface_ids>
+      <max_particles>200</max_particles>
+      <cellfrom>2</cellfrom>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-a01/results_true.dat b/tests/regression_tests/surface_source_write/case-a01/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-a01/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-a01/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-a01/surface_source_true.h5
new file mode 100644
index 000000000..4a8623567
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-a01/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d01/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d01/inputs_true.dat
new file mode 100644
index 000000000..b6157ba9b
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d01/inputs_true.dat
@@ -0,0 +1,33 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d01/results_true.dat b/tests/regression_tests/surface_source_write/case-d01/results_true.dat
new file mode 100644
index 000000000..7fb415cb1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d01/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.263843E-01 3.193920E-02
diff --git a/tests/regression_tests/surface_source_write/case-d01/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d01/surface_source_true.h5
new file mode 100644
index 000000000..8f8cd604d
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d01/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d02/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d02/inputs_true.dat
new file mode 100644
index 000000000..6877634f1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d02/inputs_true.dat
@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>1</surface_ids>
+      <max_particles>300</max_particles>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d02/results_true.dat b/tests/regression_tests/surface_source_write/case-d02/results_true.dat
new file mode 100644
index 000000000..7fb415cb1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d02/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.263843E-01 3.193920E-02
diff --git a/tests/regression_tests/surface_source_write/case-d02/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d02/surface_source_true.h5
new file mode 100644
index 000000000..3abc471e0
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d02/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d03/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d03/inputs_true.dat
new file mode 100644
index 000000000..ddcde89e4
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d03/inputs_true.dat
@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cell>2</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d03/results_true.dat b/tests/regression_tests/surface_source_write/case-d03/results_true.dat
new file mode 100644
index 000000000..7fb415cb1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d03/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.263843E-01 3.193920E-02
diff --git a/tests/regression_tests/surface_source_write/case-d03/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d03/surface_source_true.h5
new file mode 100644
index 000000000..8435e8f17
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d03/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d04/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d04/inputs_true.dat
new file mode 100644
index 000000000..7695ecc79
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d04/inputs_true.dat
@@ -0,0 +1,35 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>1</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>2</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d04/results_true.dat b/tests/regression_tests/surface_source_write/case-d04/results_true.dat
new file mode 100644
index 000000000..7fb415cb1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d04/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.263843E-01 3.193920E-02
diff --git a/tests/regression_tests/surface_source_write/case-d04/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d04/surface_source_true.h5
new file mode 100644
index 000000000..6c45d92e6
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d04/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d05/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d05/inputs_true.dat
new file mode 100644
index 000000000..865633fbb
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d05/inputs_true.dat
@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cellfrom>2</cellfrom>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d05/results_true.dat b/tests/regression_tests/surface_source_write/case-d05/results_true.dat
new file mode 100644
index 000000000..7fb415cb1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d05/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.263843E-01 3.193920E-02
diff --git a/tests/regression_tests/surface_source_write/case-d05/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d05/surface_source_true.h5
new file mode 100644
index 000000000..24b04a862
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d05/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d06/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d06/inputs_true.dat
new file mode 100644
index 000000000..422670736
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d06/inputs_true.dat
@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cellto>2</cellto>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d06/results_true.dat b/tests/regression_tests/surface_source_write/case-d06/results_true.dat
new file mode 100644
index 000000000..7fb415cb1
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d06/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.263843E-01 3.193920E-02
diff --git a/tests/regression_tests/surface_source_write/case-d06/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d06/surface_source_true.h5
new file mode 100644
index 000000000..3eb95bef9
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d06/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d07/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d07/inputs_true.dat
new file mode 100644
index 000000000..77498b64d
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d07/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell fill="1" id="8" region="101 -102 103 -104 105 -106" universe="2"/>
+    <cell id="9" material="41" region="107 -108 109 -110 111 -112 (-101 | 102 | -103 | 104 | -105 | 106)" universe="2"/>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+    <surface coeffs="-22.0" id="101" type="z-plane"/>
+    <surface coeffs="22.0" id="102" type="z-plane"/>
+    <surface coeffs="-22.0" id="103" type="x-plane"/>
+    <surface coeffs="22.0" id="104" type="x-plane"/>
+    <surface coeffs="-22.0" id="105" type="y-plane"/>
+    <surface coeffs="22.0" id="106" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-24.0" id="107" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="24.0" id="108" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-24.0" id="109" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="24.0" id="110" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-24.0" id="111" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="24.0" id="112" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>101 102 103 104 105 106</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>7</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d07/results_true.dat b/tests/regression_tests/surface_source_write/case-d07/results_true.dat
new file mode 100644
index 000000000..63ea1a64d
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d07/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.756086E-01 4.639209E-02
diff --git a/tests/regression_tests/surface_source_write/case-d07/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d07/surface_source_true.h5
new file mode 100644
index 000000000..0332f2dcf
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d07/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-d08/inputs_true.dat b/tests/regression_tests/surface_source_write/case-d08/inputs_true.dat
new file mode 100644
index 000000000..561a38337
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d08/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="40" name="no-void fuel">
+      <density units="g/cc" value="11"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="41" name="water">
+      <density units="g/cc" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell fill="1" id="8" region="101 -102 103 -104 105 -106" universe="2"/>
+    <cell id="9" material="41" region="107 -108 109 -110 111 -112 (-101 | 102 | -103 | 104 | -105 | 106)" universe="2"/>
+    <dagmc_universe filename="../../dagmc/legacy/dagmc.h5m" id="1"/>
+    <surface coeffs="-22.0" id="101" type="z-plane"/>
+    <surface coeffs="22.0" id="102" type="z-plane"/>
+    <surface coeffs="-22.0" id="103" type="x-plane"/>
+    <surface coeffs="22.0" id="104" type="x-plane"/>
+    <surface coeffs="-22.0" id="105" type="y-plane"/>
+    <surface coeffs="22.0" id="106" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-24.0" id="107" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="24.0" id="108" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-24.0" id="109" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="24.0" id="110" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-24.0" id="111" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="24.0" id="112" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-4 -4 -20 4 4 20</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>101 102 103 104 105 106</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>8</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-d08/results_true.dat b/tests/regression_tests/surface_source_write/case-d08/results_true.dat
new file mode 100644
index 000000000..63ea1a64d
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-d08/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+9.756086E-01 4.639209E-02
diff --git a/tests/regression_tests/surface_source_write/case-d08/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-d08/surface_source_true.h5
new file mode 100644
index 000000000..252559204
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-d08/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-e01/inputs_true.dat b/tests/regression_tests/surface_source_write/case-e01/inputs_true.dat
new file mode 100644
index 000000000..99a5a2899
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-e01/inputs_true.dat
@@ -0,0 +1,56 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>2</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-e01/results_true.dat b/tests/regression_tests/surface_source_write/case-e01/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-e01/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-e01/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-e01/surface_source_true.h5
new file mode 100644
index 000000000..3047519e8
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-e01/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-e02/inputs_true.dat b/tests/regression_tests/surface_source_write/case-e02/inputs_true.dat
new file mode 100644
index 000000000..ff73358ca
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-e02/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <cell id="4" material="2" region="-11 10 -13 12 -15 14 (5 | -4 | 7 | -6 | 9 | -8)" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface coeffs="-3.0" id="4" type="x-plane"/>
+    <surface coeffs="3.0" id="5" type="x-plane"/>
+    <surface coeffs="-3.0" id="6" type="y-plane"/>
+    <surface coeffs="3.0" id="7" type="y-plane"/>
+    <surface coeffs="-3.0" id="8" type="z-plane"/>
+    <surface coeffs="3.0" id="9" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="10" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="11" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="12" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="13" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-4.0" id="14" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="4.0" id="15" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <max_particles>300</max_particles>
+      <cell>3</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-e02/results_true.dat b/tests/regression_tests/surface_source_write/case-e02/results_true.dat
new file mode 100644
index 000000000..8979eb554
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-e02/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+5.169123E-02 9.144814E-03
diff --git a/tests/regression_tests/surface_source_write/case-e02/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-e02/surface_source_true.h5
new file mode 100644
index 000000000..806152853
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-e02/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/case-e03/inputs_true.dat b/tests/regression_tests/surface_source_write/case-e03/inputs_true.dat
new file mode 100644
index 000000000..a3421298b
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-e03/inputs_true.dat
@@ -0,0 +1,49 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="11.0"/>
+      <nuclide ao="0.0004524" name="U234"/>
+      <nuclide ao="0.0506068" name="U235"/>
+      <nuclide ao="0.948709" name="U238"/>
+      <nuclide ao="0.0002318" name="U236"/>
+      <nuclide ao="2.0" name="O16"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="2.0" name="H1"/>
+      <nuclide ao="1.0" name="O16"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1 2 -3" universe="1"/>
+    <cell id="2" material="2" region="1 | -2 | 3" universe="1"/>
+    <cell fill="1" id="3" region="-5 4 -7 6 -9 8" universe="2"/>
+    <surface coeffs="0.0 0.0 2.0" id="1" type="z-cylinder"/>
+    <surface coeffs="-2.0" id="2" type="z-plane"/>
+    <surface coeffs="2.0" id="3" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="4" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="5" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="6" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="7" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="-3.0" id="8" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="3.0" id="9" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>5</batches>
+    <inactive>1</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="fission">
+        <parameters>-2.0 -2.0 -2.0 2.0 2.0 2.0</parameters>
+      </space>
+    </source>
+    <surf_source_write>
+      <surface_ids>4 5 6 7 8 9</surface_ids>
+      <max_particles>300</max_particles>
+      <cell>3</cell>
+    </surf_source_write>
+    <seed>1</seed>
+  </settings>
+  </model>
diff --git a/tests/regression_tests/surface_source_write/case-e03/results_true.dat b/tests/regression_tests/surface_source_write/case-e03/results_true.dat
new file mode 100644
index 000000000..d793a7e42
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/case-e03/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+4.752377E-02 6.773395E-03
diff --git a/tests/regression_tests/surface_source_write/case-e03/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-e03/surface_source_true.h5
new file mode 100644
index 000000000..becc9254a
Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-e03/surface_source_true.h5 differ
diff --git a/tests/regression_tests/surface_source_write/test.py b/tests/regression_tests/surface_source_write/test.py
new file mode 100644
index 000000000..f949d8085
--- /dev/null
+++ b/tests/regression_tests/surface_source_write/test.py
@@ -0,0 +1,1126 @@
+"""Test the 'surface_source_write' setting.
+
+Results
+-------
+
+All results are generated using only 1 MPI process.
+
+All results are generated using 1 thread except for "test_consistency_low_realization_number".
+This specific test verifies that when the number of realization (i.e., point being candidate
+to be stored) is lower than the capacity, results are reproducible even with multiple
+threads (i.e., there is no potential thread competition that would produce different
+results in that case).
+
+All results are generated using the history-based mode except for cases e01 to e03.
+
+All results are visually verified using the '_visualize.py' script in the regression test folder.
+
+OpenMC models
+-------------
+
+Four OpenMC models with CSG-only geometries are used to cover the transmission, vacuum,
+reflective and periodic Boundary Conditions (BC):
+
+- model_1: cylindrical core in 2 boxes (vacuum and transmission BC),
+- model_2: cylindrical core in 1 box (vacuum BC),
+- model_3: cylindrical core in 1 box (reflective BC),
+- model_4: cylindrical core in 1 box (periodic BC).
+
+Two models including DAGMC geometries are also used, based on the mesh file 'dagmc.h5m'
+available from tests/regression_tests/dagmc/legacy:
+
+- model_dagmc_1: model adapted from tests/regression_tests/dagmc/legacy,
+- model_dagmc_2: model_dagmc_1 contained in two CSG boxes to introduce multiple level of coordinates.
+
+Test cases
+----------
+
+Test cases using CSG-only geometries:
+
+========  =======  =========  =========================  =====  ===================================
+Folder    Model    Surface    Cell                       BC*    Expected particles
+========  =======  =========  =========================  =====  ===================================
+case-01   model_1  No         No                         T+V    Particles crossing any surface in
+                                                                the model
+case-02   model_1  1          No                         T      Particles crossing this surface
+                                                                only
+case-03   model_1  Multiple   No                         T      Particles crossing the declared
+                                                                surfaces
+case-04   model_1  Multiple   cell (lower universe)      T      Particles crossing the declared
+                                                                surfaces that come from or are
+                                                                coming to the cell
+case-05   model_1  Multiple   cell (root universe)       T      Particles crossing the declared
+                                                                surfaces that come from or are
+                                                                coming to the cell
+case-06   model_1  No         cell (lower universe)      T      Particles crossing any surface that
+                                                                come from or are coming to the cell
+case-07   model_1  No         cell (root universe)       T      Particles crossing any surface that
+                                                                come from or are coming to the cell
+case-08   model_1  No         cellfrom (lower universe)  T      Particles crossing any surface that
+                                                                come from the cell
+case-09   model_1  No         cellto (lower universe)    T      Particles crossing any surface that
+                                                                are coming to the cell
+case-10   model_1  No         cellfrom (root universe)   T      Particles crossing any surface that
+                                                                come from the cell
+case-11   model_1  No         cellto (root universe)     T      Particles crossing any surface that
+                                                                are coming to the cell
+case-12   model_2  Multiple   No                         V      Particles crossing the declared
+                                                                surfaces
+case-13   model_2  Multiple   cell (root universe)       V      Particles crossing any surface that
+                                                                come from or are coming to the cell
+case-14   model_2  Multiple   cellfrom (root universe)   V      Particles crossing any surface that
+                                                                are coming to the cell
+case-15   model_2  Multiple   cellto (root universe)     V      None
+case-16   model_3  Multiple   No                         R      Particles crossing the declared
+                                                                surfaces
+case-17   model_3  Multiple   cell (root universe)       R      None
+case-18   model_3  Multiple   cellfrom (root universe)   R      None
+case-19   model_3  Multiple   cellto (root universe)     R      None
+case-20   model_4  1          No                         P+R    Particles crossing the declared
+                                                                periodic surface
+case-21   model_4  1          cell (root universe)       P+R    None
+========  =======  =========  =========================  =====  ===================================
+
+*: BC stands for Boundary Conditions, T for Transmission, R for Reflective, and V for Vacuum.
+
+An additional case, called 'case-a01', is used to check that the results are comparable when
+the number of threads is set to 2 if the number of realization is lower than the capacity.
+
+Cases e01 to e03 are the event-based cases corresponding to the history-based cases 04, 07 and 13,
+respectively.
+
+Test cases using DAGMC geometries:
+
+========  =============  =========  =====================  =====  ===================================
+Folder    Model          Surface    Cell                   BC*    Expected particles
+========  =============  =========  =====================  =====  ===================================
+case-d01  model_dagmc_1  No         No                     T+V    Particles crossing any surface in
+                                                                  the model
+case-d02  model_dagmc_1  1          No                     T      Particles crossing this surface
+                                                                  only
+case-d03  model_dagmc_1  No         cell                   T      Particles crossing any surface that
+                                                                  come from or are coming to the cell
+case-d04  model_dagmc_1  1          cell                   T      Particles crossing the declared
+                                                                  surface that come from or are
+                                                                  coming to the cell
+case-d05  model_dagmc_1  No         cellfrom               T      Particles crossing any surface that
+                                                                  come from the cell
+case-d06  model_dagmc_1  No         cellto                 T      Particles crossing any surface that
+                                                                  are coming to the cell
+case-d07  model_dagmc_2  Multiple   cell (lower universe)  T      Particles crossing the declared
+                                                                  surfaces that come from or are
+                                                                  coming to the cell
+case-d08  model_dagmc_2  Multiple   cell (root universe)   T      Particles crossing the declared
+                                                                  surfaces that come from or are
+                                                                  coming to the cell
+========  =============  =========  =====================  =====  ===================================
+
+*: BC stands for Boundary Conditions, T for Transmission, and V for Vacuum.
+
+Notes:
+
+- The test cases list is non-exhaustive compared to the number of possible combinations.
+  Test cases have been selected based on use and internal code logic.
+- Cases 08 to 11 are testing that the feature still works even if the level of coordinates
+  before and after crossing a surface is different,
+- Tests on boundary conditions are not performed on DAGMC models as the logic is shared
+  with CSG-only models,
+- Cases that should return an error are tested in the 'test_exceptions' unit test
+  from 'unit_tests/surface_source_write/test.py'.
+
+TODO:
+
+- Test with a lattice.
+
+"""
+
+import os
+import shutil
+from pathlib import Path
+
+import h5py
+import numpy as np
+import openmc
+import openmc.lib
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+from tests.regression_tests import config
+
+
+@pytest.fixture(scope="function")
+def single_thread(monkeypatch):
+    """Set the number of OMP threads to 1 for the test."""
+    monkeypatch.setenv("OMP_NUM_THREADS", "1")
+
+
+@pytest.fixture(scope="function")
+def two_threads(monkeypatch):
+    """Set the number of OMP threads to 2 for the test."""
+    monkeypatch.setenv("OMP_NUM_THREADS", "2")
+
+
+@pytest.fixture(scope="function")
+def single_process(monkeypatch):
+    """Set the number of MPI process to 1 for the test."""
+    monkeypatch.setitem(config, "mpi_np", "1")
+
+
+@pytest.fixture(scope="module")
+def model_1():
+    """Cylindrical core contained in a first box which is contained in a larger box.
+    A lower universe is used to describe the interior of the first box which
+    contains the core and its surrounding space.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    fuel = openmc.Material()
+    fuel.add_nuclide("U234", 0.0004524)
+    fuel.add_nuclide("U235", 0.0506068)
+    fuel.add_nuclide("U238", 0.9487090)
+    fuel.add_nuclide("U236", 0.0002318)
+    fuel.add_nuclide("O16", 2.0)
+    fuel.set_density("g/cm3", 11.0)
+
+    water = openmc.Material()
+    water.add_nuclide("H1", 2.0)
+    water.add_nuclide("O16", 1.0)
+    water.set_density("g/cm3", 1.0)
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    # -----------------------------------------------------------------------------
+    # Cylindrical core
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    core_radius = 2.0
+    core_height = 4.0
+
+    # Surfaces
+    core_cylinder = openmc.ZCylinder(r=core_radius)
+    core_lower_plane = openmc.ZPlane(-core_height / 2.0)
+    core_upper_plane = openmc.ZPlane(core_height / 2.0)
+
+    # Region
+    core_region = -core_cylinder & +core_lower_plane & -core_upper_plane
+
+    # Cells
+    core = openmc.Cell(fill=fuel, region=core_region)
+    outside_core_region = +core_cylinder | -core_lower_plane | +core_upper_plane
+    outside_core = openmc.Cell(fill=water, region=outside_core_region)
+
+    # Universe
+    inside_box1_universe = openmc.Universe(cells=[core, outside_core])
+
+    # -----------------------------------------------------------------------------
+    # Box 1
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box1_size = 6.0
+
+    # Surfaces
+    box1_rpp = openmc.model.RectangularParallelepiped(
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+    )
+
+    # Cell
+    box1 = openmc.Cell(fill=inside_box1_universe, region=-box1_rpp)
+
+    # -----------------------------------------------------------------------------
+    # Box 2
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box2_size = 8
+
+    # Surfaces
+    box2_rpp = openmc.model.RectangularParallelepiped(
+        -box2_size / 2.0, box2_size / 2.0,
+        -box2_size / 2.0, box2_size / 2.0,
+        -box2_size / 2.0, box2_size / 2.0,
+        boundary_type="vacuum"
+    )
+
+    # Cell
+    box2 = openmc.Cell(fill=water, region=-box2_rpp & +box1_rpp)
+
+    # Register geometry
+    model.geometry = openmc.Geometry([box1, box2])
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 100
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    bounds = [
+        -core_radius,
+        -core_radius,
+        -core_height / 2.0,
+        core_radius,
+        core_radius,
+        core_height / 2.0,
+    ]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    return model
+
+
+@pytest.fixture
+def model_2():
+    """Cylindrical core contained in a box.
+    A lower universe is used to describe the interior of the box which
+    contains the core and its surrounding space.
+
+    The box is defined with vacuum boundary conditions.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    fuel = openmc.Material()
+    fuel.add_nuclide("U234", 0.0004524)
+    fuel.add_nuclide("U235", 0.0506068)
+    fuel.add_nuclide("U238", 0.9487090)
+    fuel.add_nuclide("U236", 0.0002318)
+    fuel.add_nuclide("O16", 2.0)
+    fuel.set_density("g/cm3", 11.0)
+
+    water = openmc.Material()
+    water.add_nuclide("H1", 2.0)
+    water.add_nuclide("O16", 1.0)
+    water.set_density("g/cm3", 1.0)
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    # -----------------------------------------------------------------------------
+    # Cylindrical core
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    core_radius = 2.0
+    core_height = 4.0
+
+    # Surfaces
+    core_cylinder = openmc.ZCylinder(r=core_radius)
+    core_lower_plane = openmc.ZPlane(-core_height / 2.0)
+    core_upper_plane = openmc.ZPlane(core_height / 2.0)
+
+    # Region
+    core_region = -core_cylinder & +core_lower_plane & -core_upper_plane
+
+    # Cells
+    core = openmc.Cell(fill=fuel, region=core_region)
+    outside_core_region = +core_cylinder | -core_lower_plane | +core_upper_plane
+    outside_core = openmc.Cell(fill=water, region=outside_core_region)
+
+    # Universe
+    inside_box1_universe = openmc.Universe(cells=[core, outside_core])
+
+    # -----------------------------------------------------------------------------
+    # Box 1
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box1_size = 6.0
+
+    # Surfaces
+    box1_rpp = openmc.model.RectangularParallelepiped(
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+        boundary_type="vacuum"
+    )
+
+    # Cell
+    box1 = openmc.Cell(fill=inside_box1_universe, region=-box1_rpp)
+
+    # Register geometry
+    model.geometry = openmc.Geometry([box1])
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 100
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    bounds = [
+        -core_radius,
+        -core_radius,
+        -core_height / 2.0,
+        core_radius,
+        core_radius,
+        core_height / 2.0,
+    ]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    return model
+
+
+@pytest.fixture
+def model_3():
+    """Cylindrical core contained in a box.
+    A lower universe is used to describe the interior of the box which
+    contains the core and its surrounding space.
+
+    The box is defined with reflective boundary conditions.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    fuel = openmc.Material()
+    fuel.add_nuclide("U234", 0.0004524)
+    fuel.add_nuclide("U235", 0.0506068)
+    fuel.add_nuclide("U238", 0.9487090)
+    fuel.add_nuclide("U236", 0.0002318)
+    fuel.add_nuclide("O16", 2.0)
+    fuel.set_density("g/cm3", 11.0)
+
+    water = openmc.Material()
+    water.add_nuclide("H1", 2.0)
+    water.add_nuclide("O16", 1.0)
+    water.set_density("g/cm3", 1.0)
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    # -----------------------------------------------------------------------------
+    # Cylindrical core
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    core_radius = 2.0
+    core_height = 4.0
+
+    # Surfaces
+    core_cylinder = openmc.ZCylinder(r=core_radius)
+    core_lower_plane = openmc.ZPlane(-core_height / 2.0)
+    core_upper_plane = openmc.ZPlane(core_height / 2.0)
+
+    # Region
+    core_region = -core_cylinder & +core_lower_plane & -core_upper_plane
+
+    # Cells
+    core = openmc.Cell(fill=fuel, region=core_region)
+    outside_core_region = +core_cylinder | -core_lower_plane | +core_upper_plane
+    outside_core = openmc.Cell(fill=water, region=outside_core_region)
+
+    # Universe
+    inside_box1_universe = openmc.Universe(cells=[core, outside_core])
+
+    # -----------------------------------------------------------------------------
+    # Box 1
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box1_size = 6.0
+
+    # Surfaces
+    box1_rpp = openmc.model.RectangularParallelepiped(
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+        -box1_size / 2.0, box1_size / 2.0,
+        boundary_type="reflective"
+    )
+
+    # Cell
+    box1 = openmc.Cell(fill=inside_box1_universe, region=-box1_rpp)
+
+    # Register geometry
+    model.geometry = openmc.Geometry([box1])
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 100
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    bounds = [
+        -core_radius,
+        -core_radius,
+        -core_height / 2.0,
+        core_radius,
+        core_radius,
+        core_height / 2.0,
+    ]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    return model
+
+
+@pytest.fixture
+def model_4():
+    """Cylindrical core contained in a box.
+    A lower universe is used to describe the interior of the box which
+    contains the core and its surrounding space.
+
+    The box is defined with a pair of periodic boundary with reflective
+    boundaries.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    fuel = openmc.Material()
+    fuel.add_nuclide("U234", 0.0004524)
+    fuel.add_nuclide("U235", 0.0506068)
+    fuel.add_nuclide("U238", 0.9487090)
+    fuel.add_nuclide("U236", 0.0002318)
+    fuel.add_nuclide("O16", 2.0)
+    fuel.set_density("g/cm3", 11.0)
+
+    water = openmc.Material()
+    water.add_nuclide("H1", 2.0)
+    water.add_nuclide("O16", 1.0)
+    water.set_density("g/cm3", 1.0)
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    # -----------------------------------------------------------------------------
+    # Cylindrical core
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    core_radius = 2.0
+    core_height = 4.0
+
+    # Surfaces
+    core_cylinder = openmc.ZCylinder(r=core_radius)
+    core_lower_plane = openmc.ZPlane(-core_height / 2.0)
+    core_upper_plane = openmc.ZPlane(core_height / 2.0)
+
+    # Region
+    core_region = -core_cylinder & +core_lower_plane & -core_upper_plane
+
+    # Cells
+    core = openmc.Cell(fill=fuel, region=core_region)
+    outside_core_region = +core_cylinder | -core_lower_plane | +core_upper_plane
+    outside_core = openmc.Cell(fill=water, region=outside_core_region)
+
+    # Universe
+    inside_box1_universe = openmc.Universe(cells=[core, outside_core])
+
+    # -----------------------------------------------------------------------------
+    # Box 1
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box1_size = 6.0
+
+    # Surfaces
+    box1_lower_plane = openmc.ZPlane(-box1_size / 2.0, boundary_type="periodic")
+    box1_upper_plane = openmc.ZPlane(box1_size / 2.0, boundary_type="periodic")
+    box1_left_plane = openmc.XPlane(-box1_size / 2.0, boundary_type="reflective")
+    box1_right_plane = openmc.XPlane(box1_size / 2.0, boundary_type="reflective")
+    box1_rear_plane = openmc.YPlane(-box1_size / 2.0, boundary_type="reflective")
+    box1_front_plane = openmc.YPlane(box1_size / 2.0, boundary_type="reflective")
+
+    # Region
+    box1_region = (
+        +box1_lower_plane
+        & -box1_upper_plane
+        & +box1_left_plane
+        & -box1_right_plane
+        & +box1_rear_plane
+        & -box1_front_plane
+    )
+
+    # Cell
+    box1 = openmc.Cell(fill=inside_box1_universe, region=box1_region)
+
+    # Register geometry
+    model.geometry = openmc.Geometry([box1])
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 100
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    bounds = [
+        -core_radius,
+        -core_radius,
+        -core_height / 2.0,
+        core_radius,
+        core_radius,
+        core_height / 2.0,
+    ]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    return model
+
+
+def return_surface_source_data(filepath):
+    """Read a surface source file and return a sorted array composed
+    of flatten arrays of source data for each surface source point.
+
+    TODO:
+
+    - use read_source_file from source.py instead. Or a dedicated function
+      to produce sorted list of source points for a given file.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the surface source file
+
+    Returns
+    -------
+    data : np.array
+        Sorted array composed of flatten arrays of source data for
+        each surface source point
+
+    """
+    data = []
+    keys = []
+
+    # Read source file
+    with h5py.File(filepath, "r") as f:
+        for point in f["source_bank"]:
+            r = point["r"]
+            u = point["u"]
+            e = point["E"]
+            time = point["time"]
+            wgt = point["wgt"]
+            delayed_group = point["delayed_group"]
+            surf_id = point["surf_id"]
+            particle = point["particle"]
+
+            key = (
+                f"{r[0]:.10e} {r[1]:.10e} {r[2]:.10e} {u[0]:.10e} {u[1]:.10e} {u[2]:.10e}"
+                f"{e:.10e} {time:.10e} {wgt:.10e} {delayed_group} {surf_id} {particle}"
+            )
+
+            keys.append(key)
+
+            values = [*r, *u, e, time, wgt, delayed_group, surf_id, particle]
+            assert len(values) == 12
+            data.append(values)
+
+    data = np.array(data)
+    keys = np.array(keys)
+    sorted_idx = np.argsort(keys)
+
+    return data[sorted_idx]
+
+
+class SurfaceSourceWriteTestHarness(PyAPITestHarness):
+    def __init__(self, statepoint_name, model=None, inputs_true=None, workdir=None):
+        super().__init__(statepoint_name, model, inputs_true)
+        self.workdir = workdir
+
+    def _test_output_created(self):
+        """Make sure surface_source.h5 has also been created."""
+        super()._test_output_created()
+        if self._model.settings.surf_source_write:
+            assert os.path.exists(
+                "surface_source.h5"
+            ), "Surface source file has not been created."
+
+    def _compare_output(self):
+        """Compare surface_source.h5 files."""
+        if self._model.settings.surf_source_write:
+            source_true = return_surface_source_data("surface_source_true.h5")
+            source_test = return_surface_source_data("surface_source.h5")
+            np.testing.assert_allclose(source_true, source_test, rtol=1e-07)
+
+    def main(self):
+        """Accept commandline arguments and either run or update tests."""
+        if config["build_inputs"]:
+            self.build_inputs()
+        elif config["update"]:
+            self.update_results()
+        else:
+            self.execute_test()
+
+    def build_inputs(self):
+        """Build inputs."""
+        base_dir = os.getcwd()
+        try:
+            os.chdir(self.workdir)
+            self._build_inputs()
+        finally:
+            os.chdir(base_dir)
+
+    def execute_test(self):
+        """Build inputs, run OpenMC, and verify correct results."""
+        base_dir = os.getcwd()
+        try:
+            os.chdir(self.workdir)
+            self._build_inputs()
+            inputs = self._get_inputs()
+            self._write_inputs(inputs)
+            self._compare_inputs()
+            self._run_openmc()
+            self._test_output_created()
+            self._compare_output()
+            results = self._get_results()
+            self._write_results(results)
+            self._compare_results()
+        finally:
+            self._cleanup()
+            os.chdir(base_dir)
+
+    def update_results(self):
+        """Update results_true.dat and inputs_true.dat"""
+        base_dir = os.getcwd()
+        try:
+            os.chdir(self.workdir)
+            self._build_inputs()
+            inputs = self._get_inputs()
+            self._write_inputs(inputs)
+            self._overwrite_inputs()
+            self._run_openmc()
+            self._test_output_created()
+            results = self._get_results()
+            self._write_results(results)
+            self._overwrite_results()
+        finally:
+            self._cleanup()
+            os.chdir(base_dir)
+
+    def _overwrite_results(self):
+        """Also add the 'surface_source.h5' file during overwriting."""
+        super()._overwrite_results()
+        if os.path.exists("surface_source.h5"):
+            shutil.copyfile("surface_source.h5", "surface_source_true.h5")
+
+    def _cleanup(self):
+        """Also remove the 'surface_source.h5' file while cleaning."""
+        super()._cleanup()
+        fs = "surface_source.h5"
+        if os.path.exists(fs):
+            os.remove(fs)
+
+
+@pytest.mark.skipif(config["event"] is True, reason="Results from history-based mode.")
+@pytest.mark.parametrize(
+    "folder, model_name, parameter",
+    [
+        ("case-01", "model_1", {"max_particles": 300}),
+        ("case-02", "model_1", {"max_particles": 300, "surface_ids": [8]}),
+        (
+            "case-03",
+            "model_1",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9]},
+        ),
+        (
+            "case-04",
+            "model_1",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cell": 2},
+        ),
+        (
+            "case-05",
+            "model_1",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cell": 3},
+        ),
+        ("case-06", "model_1", {"max_particles": 300, "cell": 2}),
+        ("case-07", "model_1", {"max_particles": 300, "cell": 3}),
+        ("case-08", "model_1", {"max_particles": 300, "cellfrom": 2}),
+        ("case-09", "model_1", {"max_particles": 300, "cellto": 2}),
+        ("case-10", "model_1", {"max_particles": 300, "cellfrom": 3}),
+        ("case-11", "model_1", {"max_particles": 300, "cellto": 3}),
+        (
+            "case-12",
+            "model_2",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9]},
+        ),
+        (
+            "case-13",
+            "model_2",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cell": 3},
+        ),
+        (
+            "case-14",
+            "model_2",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cellfrom": 3},
+        ),
+        (
+            "case-15",
+            "model_2",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cellto": 3},
+        ),
+        (
+            "case-16",
+            "model_3",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9]},
+        ),
+        (
+            "case-17",
+            "model_3",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cell": 3},
+        ),
+        (
+            "case-18",
+            "model_3",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cellfrom": 3},
+        ),
+        (
+            "case-19",
+            "model_3",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cellto": 3},
+        ),
+        (
+            "case-20",
+            "model_4",
+            {"max_particles": 300, "surface_ids": [4]},
+        ),
+        (
+            "case-21",
+            "model_4",
+            {"max_particles": 300, "surface_ids": [4], "cell": 3},
+        ),
+    ],
+)
+def test_surface_source_cell_history_based(
+    folder, model_name, parameter, single_thread, single_process, request
+):
+    """Test on history-based results for CSG-only geometries."""
+    assert os.environ["OMP_NUM_THREADS"] == "1"
+    assert config["mpi_np"] == "1"
+    model = request.getfixturevalue(model_name)
+    model.settings.surf_source_write = parameter
+    harness = SurfaceSourceWriteTestHarness(
+        "statepoint.5.h5", model=model, workdir=folder
+    )
+    harness.main()
+
+
+@pytest.mark.skipif(config["event"] is True, reason="Results from history-based mode.")
+def test_consistency_low_realization_number(model_1, two_threads, single_process):
+    """The objective is to test that the results produced, in a case where
+    the number of potential realization (particle storage) is low
+    compared to the capacity of storage, are still consistent.
+
+    This configuration ensures that the competition between threads does not
+    occur and that the content of the source file created can be compared.
+
+    """
+    assert os.environ["OMP_NUM_THREADS"] == "2"
+    assert config["mpi_np"] == "1"
+    model_1.settings.surf_source_write = {
+        "max_particles": 200,
+        "surface_ids": [1, 2, 3],
+        "cellfrom": 2,
+    }
+    harness = SurfaceSourceWriteTestHarness(
+        "statepoint.5.h5", model=model_1, workdir="case-a01"
+    )
+    harness.main()
+
+
+@pytest.mark.skipif(config["event"] is False, reason="Results from event-based mode.")
+@pytest.mark.parametrize(
+    "folder, model_name, parameter",
+    [
+        (
+            "case-e01",
+            "model_1",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cell": 2},
+        ),
+        ("case-e02", "model_1", {"max_particles": 300, "cell": 3}),
+        (
+            "case-e03",
+            "model_2",
+            {"max_particles": 300, "surface_ids": [4, 5, 6, 7, 8, 9], "cell": 3},
+        ),
+    ],
+)
+def test_surface_source_cell_event_based(
+    folder, model_name, parameter, single_thread, single_process, request
+):
+    """Test on event-based results for CSG-only geometries."""
+    assert os.environ["OMP_NUM_THREADS"] == "1"
+    assert config["mpi_np"] == "1"
+    model = request.getfixturevalue(model_name)
+    model.settings.surf_source_write = parameter
+    harness = SurfaceSourceWriteTestHarness(
+        "statepoint.5.h5", model=model, workdir=folder
+    )
+    harness.main()
+
+
+@pytest.fixture(scope="module")
+def model_dagmc_1():
+    """Model based on the mesh file 'dagmc.h5m' available from
+    tests/regression_tests/dagmc/legacy.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    u235 = openmc.Material(name="no-void fuel")
+    u235.add_nuclide("U235", 1.0, "ao")
+    u235.set_density("g/cc", 11)
+    u235.id = 40
+
+    water = openmc.Material(name="water")
+    water.add_nuclide("H1", 2.0, "ao")
+    water.add_nuclide("O16", 1.0, "ao")
+    water.set_density("g/cc", 1.0)
+    water.add_s_alpha_beta("c_H_in_H2O")
+    water.id = 41
+
+    materials = openmc.Materials([u235, water])
+    model.materials = materials
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    dagmc_univ = openmc.DAGMCUniverse(Path("../../dagmc/legacy/dagmc.h5m"))
+    model.geometry = openmc.Geometry(dagmc_univ)
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 100
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    source_box = openmc.stats.Box([-4, -4, -20], [4, 4, 20], only_fissionable=True)
+    model.settings.source = openmc.IndependentSource(space=source_box)
+
+    return model
+
+
+@pytest.fixture(scope="module")
+def model_dagmc_2():
+    """Model based on the mesh file 'dagmc.h5m' available from
+    tests/regression_tests/dagmc/legacy.
+
+    This model corresponds to the model_dagmc_1 contained in two boxes to introduce
+    multiple level of coordinates from CSG geometry.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    u235 = openmc.Material(name="no-void fuel")
+    u235.add_nuclide("U235", 1.0, "ao")
+    u235.set_density("g/cc", 11)
+    u235.id = 40
+
+    water = openmc.Material(name="water")
+    water.add_nuclide("H1", 2.0, "ao")
+    water.add_nuclide("O16", 1.0, "ao")
+    water.set_density("g/cc", 1.0)
+    water.add_s_alpha_beta("c_H_in_H2O")
+    water.id = 41
+
+    materials = openmc.Materials([u235, water])
+    model.materials = materials
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+
+    dagmc_univ = openmc.DAGMCUniverse(Path("../../dagmc/legacy/dagmc.h5m"))
+
+    # -----------------------------------------------------------------------------
+    # Box 1
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box1_size = 44
+
+    # Surfaces
+    box1_lower_plane = openmc.ZPlane(-box1_size / 2.0, surface_id=101)
+    box1_upper_plane = openmc.ZPlane(box1_size / 2.0, surface_id=102)
+    box1_left_plane = openmc.XPlane(-box1_size / 2.0, surface_id=103)
+    box1_right_plane = openmc.XPlane(box1_size / 2.0, surface_id=104)
+    box1_rear_plane = openmc.YPlane(-box1_size / 2.0, surface_id=105)
+    box1_front_plane = openmc.YPlane(box1_size / 2.0, surface_id=106)
+
+    # Region
+    box1_region = (
+        +box1_lower_plane
+        & -box1_upper_plane
+        & +box1_left_plane
+        & -box1_right_plane
+        & +box1_rear_plane
+        & -box1_front_plane
+    )
+
+    # Cell
+    box1 = openmc.Cell(fill=dagmc_univ, region=box1_region, cell_id=8)
+
+    # -----------------------------------------------------------------------------
+    # Box 2
+    # -----------------------------------------------------------------------------
+
+    # Parameters
+    box2_size = 48
+
+    # Surfaces
+    box2_lower_plane = openmc.ZPlane(
+        -box2_size / 2.0, boundary_type="vacuum", surface_id=107
+    )
+    box2_upper_plane = openmc.ZPlane(
+        box2_size / 2.0, boundary_type="vacuum", surface_id=108
+    )
+    box2_left_plane = openmc.XPlane(
+        -box2_size / 2.0, boundary_type="vacuum", surface_id=109
+    )
+    box2_right_plane = openmc.XPlane(
+        box2_size / 2.0, boundary_type="vacuum", surface_id=110
+    )
+    box2_rear_plane = openmc.YPlane(
+        -box2_size / 2.0, boundary_type="vacuum", surface_id=111
+    )
+    box2_front_plane = openmc.YPlane(
+        box2_size / 2.0, boundary_type="vacuum", surface_id=112
+    )
+
+    # Region
+    inside_box2 = (
+        +box2_lower_plane
+        & -box2_upper_plane
+        & +box2_left_plane
+        & -box2_right_plane
+        & +box2_rear_plane
+        & -box2_front_plane
+    )
+    outside_box1 = (
+        -box1_lower_plane
+        | +box1_upper_plane
+        | -box1_left_plane
+        | +box1_right_plane
+        | -box1_rear_plane
+        | +box1_front_plane
+    )
+
+    box2_region = inside_box2 & outside_box1
+
+    # Cell
+    box2 = openmc.Cell(fill=water, region=box2_region, cell_id=9)
+
+    # Register geometry
+    model.geometry = openmc.Geometry([box1, box2])
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 100
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    source_box = openmc.stats.Box([-4, -4, -20], [4, 4, 20], only_fissionable=True)
+    model.settings.source = openmc.IndependentSource(space=source_box)
+
+    return model
+
+
+@pytest.mark.skipif(
+    not openmc.lib._dagmc_enabled(), reason="DAGMC CAD geometry is not enabled."
+)
+@pytest.mark.skipif(config["event"] is True, reason="Results from history-based mode.")
+@pytest.mark.parametrize(
+    "folder, model_name, parameter",
+    [
+        ("case-d01", "model_dagmc_1", {"max_particles": 300}),
+        ("case-d02", "model_dagmc_1", {"max_particles": 300, "surface_ids": [1]}),
+        ("case-d03", "model_dagmc_1", {"max_particles": 300, "cell": 2}),
+        (
+            "case-d04",
+            "model_dagmc_1",
+            {"max_particles": 300, "surface_ids": [1], "cell": 2},
+        ),
+        ("case-d05", "model_dagmc_1", {"max_particles": 300, "cellfrom": 2}),
+        ("case-d06", "model_dagmc_1", {"max_particles": 300, "cellto": 2}),
+        (
+            "case-d07",
+            "model_dagmc_2",
+            {
+                "max_particles": 300,
+                "surface_ids": [101, 102, 103, 104, 105, 106],
+                "cell": 7,
+            },
+        ),
+        (
+            "case-d08",
+            "model_dagmc_2",
+            {
+                "max_particles": 300,
+                "surface_ids": [101, 102, 103, 104, 105, 106],
+                "cell": 8,
+            },
+        ),
+    ],
+)
+def test_surface_source_cell_dagmc(
+    folder, model_name, parameter, single_thread, single_process, request
+):
+    """Test on models with DAGMC geometries."""
+    assert os.environ["OMP_NUM_THREADS"] == "1"
+    assert config["mpi_np"] == "1"
+    model = request.getfixturevalue(model_name)
+    model.settings.surf_source_write = parameter
+    harness = SurfaceSourceWriteTestHarness(
+        "statepoint.5.h5", model=model, workdir=folder
+    )
+    harness.main()
diff --git a/tests/unit_tests/test_surface_source_write.py b/tests/unit_tests/test_surface_source_write.py
new file mode 100644
index 000000000..472363073
--- /dev/null
+++ b/tests/unit_tests/test_surface_source_write.py
@@ -0,0 +1,248 @@
+"""Test the 'surf_source_write' setting used to store particles that cross
+surfaces in a file for a given simulation."""
+
+from pathlib import Path
+
+import openmc
+import openmc.lib
+import pytest
+import h5py
+import numpy as np
+
+
+@pytest.fixture(scope="module")
+def geometry():
+    """Simple hydrogen sphere geometry"""
+    openmc.reset_auto_ids()
+    material = openmc.Material(name="H1")
+    material.add_element("H", 1.0)
+    sphere = openmc.Sphere(r=1.0, boundary_type="vacuum")
+    cell = openmc.Cell(region=-sphere, fill=material)
+    return openmc.Geometry([cell])
+
+
+@pytest.mark.parametrize(
+    "parameter",
+    [
+        {"max_particles": 200},
+        {"max_particles": 200, "cell": 1},
+        {"max_particles": 200, "cellto": 1},
+        {"max_particles": 200, "cellfrom": 1},
+        {"max_particles": 200, "surface_ids": [2]},
+        {"max_particles": 200, "surface_ids": [2], "cell": 1},
+        {"max_particles": 200, "surface_ids": [2], "cellto": 1},
+        {"max_particles": 200, "surface_ids": [2], "cellfrom": 1},
+    ],
+)
+def test_xml_serialization(parameter, run_in_tmpdir):
+    """Check that the different use cases can be written and read in XML."""
+    settings = openmc.Settings()
+    settings.surf_source_write = parameter
+    settings.export_to_xml()
+
+    read_settings = openmc.Settings.from_xml()
+    assert read_settings.surf_source_write == parameter
+
+
+ERROR_MSG_1 = (
+    "A maximum number of particles needs to be specified "
+    "using the 'max_particles' parameter to store surface "
+    "source points."
+)
+ERROR_MSG_2 = "'cell', 'cellfrom' and 'cellto' cannot be used at the same time."
+
+
+@pytest.mark.parametrize(
+    "parameter, error",
+    [
+        ({"cell": 1}, ERROR_MSG_1),
+        ({"max_particles": 200, "cell": 1, "cellto": 1}, ERROR_MSG_2),
+        ({"max_particles": 200, "cell": 1, "cellfrom": 1}, ERROR_MSG_2),
+        ({"max_particles": 200, "cellto": 1, "cellfrom": 1}, ERROR_MSG_2),
+        ({"max_particles": 200, "cell": 1, "cellto": 1, "cellfrom": 1}, ERROR_MSG_2),
+    ],
+)
+def test_exceptions(parameter, error, run_in_tmpdir, geometry):
+    """Test parameters configuration that should return an error."""
+    settings = openmc.Settings(run_mode="fixed source", batches=5, particles=100)
+    settings.surf_source_write = parameter
+    model = openmc.Model(geometry=geometry, settings=settings)
+    with pytest.raises(RuntimeError, match=error):
+        model.run()
+
+
+@pytest.fixture(scope="module")
+def model():
+    """Simple hydrogen sphere divided in two hemispheres
+    by a z-plane to form 2 cells."""
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # Material
+    material = openmc.Material(name="H1")
+    material.add_element("H", 1.0)
+
+    # Geometry
+    radius = 1.0
+    sphere = openmc.Sphere(r=radius, boundary_type="reflective")
+    plane = openmc.ZPlane(0.0)
+    cell_1 = openmc.Cell(region=-sphere & -plane, fill=material)
+    cell_2 = openmc.Cell(region=-sphere & +plane, fill=material)
+    root = openmc.Universe(cells=[cell_1, cell_2])
+    model.geometry = openmc.Geometry(root)
+
+    # Settings
+    model.settings = openmc.Settings()
+    model.settings.run_mode = "fixed source"
+    model.settings.particles = 100
+    model.settings.batches = 3
+    model.settings.seed = 1
+
+    bounds = [-radius, -radius, -radius, radius, radius, radius]
+    distribution = openmc.stats.Box(bounds[:3], bounds[3:])
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    return model
+
+
+@pytest.mark.parametrize(
+    "parameter",
+    [
+        {"max_particles": 200, "cellto": 2, "surface_ids": [2]},
+        {"max_particles": 200, "cellfrom": 2, "surface_ids": [2]},
+    ],
+)
+def test_particle_direction(parameter, run_in_tmpdir, model):
+    """Test the direction of particles with the 'cellfrom' and 'cellto' parameters
+    on a simple model with only one surface of interest.
+
+    Cell 2 is the upper hemisphere and surface 2 is the plane dividing the sphere
+    into two hemispheres.
+
+    """
+    model.settings.surf_source_write = parameter
+    model.run()
+    with h5py.File("surface_source.h5", "r") as f:
+        source = f["source_bank"]
+
+        assert len(source) == 200
+
+        # We want to verify that the dot product of the surface's normal vector
+        # and the direction of the particle is either positive or negative
+        # depending on cellfrom or cellto. In this case, it is equivalent
+        # to just compare the z component of the direction of the particle.
+        for point in source:
+            if "cellto" in parameter.keys():
+                assert point["u"]["z"] > 0.0
+            elif "cellfrom" in parameter.keys():
+                assert point["u"]["z"] < 0.0
+            else:
+                assert False
+
+
+@pytest.fixture
+def model_dagmc(request):
+    """Model based on the mesh file 'dagmc.h5m' available from
+    tests/regression_tests/dagmc/legacy.
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # =============================================================================
+    # Materials
+    # =============================================================================
+
+    u235 = openmc.Material(name="no-void fuel")
+    u235.add_nuclide("U235", 1.0, "ao")
+    u235.set_density("g/cc", 11)
+    u235.id = 40
+
+    water = openmc.Material(name="water")
+    water.add_nuclide("H1", 2.0, "ao")
+    water.add_nuclide("O16", 1.0, "ao")
+    water.set_density("g/cc", 1.0)
+    water.add_s_alpha_beta("c_H_in_H2O")
+    water.id = 41
+
+    model.materials = openmc.Materials([u235, water])
+
+    # =============================================================================
+    # Geometry
+    # =============================================================================
+    dagmc_path = Path(request.fspath).parent / "../regression_tests/dagmc/legacy/dagmc.h5m"
+    dagmc_univ = openmc.DAGMCUniverse(dagmc_path)
+    model.geometry = openmc.Geometry(dagmc_univ)
+
+    # =============================================================================
+    # Settings
+    # =============================================================================
+
+    model.settings = openmc.Settings()
+    model.settings.particles = 300
+    model.settings.batches = 5
+    model.settings.inactive = 1
+    model.settings.seed = 1
+
+    source_box = openmc.stats.Box([-4, -4, -20], [4, 4, 20])
+    model.settings.source = openmc.IndependentSource(space=source_box)
+
+    return model
+
+
+@pytest.mark.skipif(
+    not openmc.lib._dagmc_enabled(), reason="DAGMC CAD geometry is not enabled."
+)
+@pytest.mark.parametrize(
+    "parameter",
+    [
+        {"max_particles": 200, "cellto": 1},
+        {"max_particles": 200, "cellfrom": 1},
+    ],
+)
+def test_particle_direction_dagmc(parameter, run_in_tmpdir, model_dagmc):
+    """Test the direction of particles with the 'cellfrom' and 'cellto' parameters
+    on a DAGMC model."""
+    model_dagmc.settings.surf_source_write = parameter
+    model_dagmc.run()
+
+    r = 7.0
+    h = 20.0
+
+    with h5py.File("surface_source.h5", "r") as f:
+        source = f["source_bank"]
+
+        assert len(source) == 200
+
+        for point in source:
+
+            x, y, z = point["r"]
+            ux, uy, uz = point["u"]
+
+            # If the point is on the upper or lower circle
+            if np.allclose(abs(z), h):
+                # If the point is also on the cylindrical surface
+                if np.allclose(np.sqrt(x**2 + y**2), r):
+                    if "cellfrom" in parameter.keys():
+                        assert (uz * z > 0) or (ux * x + uy * y > 0)
+                    elif "cellto" in parameter.keys():
+                        assert (uz * z < 0) or (ux * x + uy * y < 0)
+                    else:
+                        assert False
+                # If the point is not on the cylindrical surface
+                else:
+                    if "cellfrom" in parameter.keys():
+                        assert uz * z > 0
+                    elif "cellto" in parameter.keys():
+                        assert uz * z < 0
+                    else:
+                        assert False
+            # If the point is not on the upper or lower circle,
+            # meaning it is on the cylindrical surface
+            else:
+                if "cellfrom" in parameter.keys():
+                    assert ux * x + uy * y > 0
+                elif "cellto" in parameter.keys():
+                    assert ux * x + uy * y < 0
+                else:
+                    assert False
diff --git a/tests/unit_tests/test_tally_multiply_density.py b/tests/unit_tests/test_tally_multiply_density.py
index 66ef41e0d..552d76bd5 100644
--- a/tests/unit_tests/test_tally_multiply_density.py
+++ b/tests/unit_tests/test_tally_multiply_density.py
@@ -3,7 +3,7 @@ import openmc
 import pytest
 
 
-def test_micro_macro_compare():
+def test_micro_macro_compare(run_in_tmpdir):
     # Create simple sphere model with H1 and H2
     mat = openmc.Material()
     mat.add_components({'H1': 1.0, 'H2': 1.0})
diff --git a/tests/unit_tests/weightwindows/test_ww_gen.py b/tests/unit_tests/weightwindows/test_ww_gen.py
index 6fc3747b1..555421461 100644
--- a/tests/unit_tests/weightwindows/test_ww_gen.py
+++ b/tests/unit_tests/weightwindows/test_ww_gen.py
@@ -101,7 +101,7 @@ def labels(params):
 
 
 @pytest.mark.parametrize("filters", test_cases, ids=labels)
-def test_ww_gen(filters, model):
+def test_ww_gen(filters, run_in_tmpdir, model):
 
     tally = openmc.Tally()
     tally.filters = list(filters)