diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 98709fb33ea..b3f862381a9 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -75,14 +75,17 @@ class Region {
   //! \param r The 3D Cartesian coordinate to check.
   //! \param u A direction used to "break ties" the coordinates are very
   //!   close to a surface.
+  //! \param t The time coordinate to check.
+  //! \param speed Particle speed to "break ties" for moving surface.
   //! \param on_surface The signed index of a surface that the coordinate is
   //!   known to be on.  This index takes precedence over surface sense
   //!   calculations.
-  bool contains(Position r, Direction u, int32_t on_surface) const;
+  bool contains(
+    Position r, Direction u, double t, double speed, int32_t on_surface) const;
 
   //! Find the oncoming boundary of this cell.
   std::pair<double, int32_t> distance(
-    Position r, Direction u, int32_t on_surface) const;
+    Position r, Direction u, double t, double speed, int32_t on_surface) const;
 
   //! Get the BoundingBox for this cell.
   BoundingBox bounding_box(int32_t cell_id) const;
@@ -108,7 +111,9 @@ class Region {
 
   //! Determine if a particle is inside the cell for a simple cell (only
   //! intersection operators)
-  bool contains_simple(Position r, Direction u, int32_t on_surface) const;
+  bool contains_simple(
+
+    Position r, Direction u, double t, double speed, int32_t on_surface) const;
 
   //! Determine if a particle is inside the cell for a complex cell.
   //!
@@ -116,6 +121,8 @@ class Region {
   //! if short circuiting can be used. Short cicuiting uses the relative and
   //! absolute depth of parentheses in the expression.
   bool contains_complex(Position r, Direction u, int32_t on_surface) const;
+  bool contains_complex(
+    Position r, Direction u, double t, double speed, int32_t on_surface) const;
 
   //! BoundingBox if the paritcle is in a simple cell.
   BoundingBox bounding_box_simple() const;
@@ -175,14 +182,17 @@ class Cell {
   //! \param r The 3D Cartesian coordinate to check.
   //! \param u A direction used to "break ties" the coordinates are very
   //!   close to a surface.
+  //! \param t The time coordinate to check.
+  //! \param speed Particle speed to "break ties" for moving surface.
   //! \param on_surface The signed index of a surface that the coordinate is
   //!   known to be on.  This index takes precedence over surface sense
   //!   calculations.
-  virtual bool contains(Position r, Direction u, int32_t on_surface) const = 0;
+  virtual bool contains(Position r, Direction u, double t, double speed,
+    int32_t on_surface) const = 0;
 
   //! Find the oncoming boundary of this cell.
   virtual std::pair<double, int32_t> distance(
-    Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0;
+    Position r, Direction u, double t, double speed, int32_t on_surface, GeometryState* p) const = 0;
 
   //! Write all information needed to reconstruct the cell to an HDF5 group.
   //! \param group_id An HDF5 group id.
@@ -371,14 +381,15 @@ class CSGCell : public Cell {
   vector<int32_t> surfaces() const override { return region_.surfaces(); }
 
   std::pair<double, int32_t> distance(Position r, Direction u,
-    int32_t on_surface, GeometryState* p) const override
+    double t, double speed, int32_t on_surface, GeometryState* p) const override
   {
-    return region_.distance(r, u, on_surface);
+    return region_.distance(r, u, t, speed, on_surface);
   }
 
-  bool contains(Position r, Direction u, int32_t on_surface) const override
+  bool contains(Position r, Direction u, double t, double speed,
+    int32_t on_surface) const override
   {
-    return region_.contains(r, u, on_surface);
+    return region_.contains(r, u, t, speed, on_surface);
   }
 
   BoundingBox bounding_box() const override
diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h
index 82ef1b644f3..ec9031af592 100644
--- a/include/openmc/dagmc.h
+++ b/include/openmc/dagmc.h
@@ -46,9 +46,11 @@ class DAGSurface : public Surface {
 
   moab::EntityHandle mesh_handle() const;
 
-  double evaluate(Position r) const override;
+  double evaluate(Position r, double t) const override;
   double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
+  double dot_normal(
+    Position r, Direction u, double t, double speed) const override;
   Direction reflect(
     Position r, Direction u, GeometryState* p = nullptr) const override;
 
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 0f37719b94f..4e97733ce64 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -37,7 +37,7 @@ class Particle : public ParticleData {
   //==========================================================================
   // Methods
 
-  double speed() const;
+  double get_speed() const;
 
   //! create a secondary particle
   //
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index be7fa74837e..508dce216e2 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -316,6 +316,16 @@ class GeometryState {
   Direction& u_local() { return coord_[n_coord_ - 1].u; }
   const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
 
+  // Accessors for time (units are seconds).
+  double& time() { return time_; }
+  const double& time() const { return time_; }
+  double& time_last() { return time_last_; }
+  const double& time_last() const { return time_last_; }
+
+  // Accessors for speed (cm/s).
+  double& speed() { return speed_; }
+  const double& speed() const { return speed_; }
+
   // Surface token for the surface that the particle is currently on
   int& surface() { return surface_; }
   const int& surface() const { return surface_; }
@@ -364,6 +374,11 @@ class GeometryState {
   Position r_last_;         //!< previous coordinates
   Direction u_last_;        //!< previous direction coordinates
 
+  double time_;      //!< time
+  double time_last_; //!< previous time
+
+  double speed_;
+
   int surface_ {
     SURFACE_NONE}; //!< surface token for surface the particle is currently on
 
@@ -437,8 +452,6 @@ class ParticleData : public GeometryState {
   double wgt_born_ {1.0};
   double wgt_ww_born_ {-1.0};
   double mu_;
-  double time_ {0.0};
-  double time_last_ {0.0};
   double wgt_last_ {1.0};
 
   bool fission_ {false};
@@ -561,13 +574,6 @@ class ParticleData : public GeometryState {
   double& mu() { return mu_; }
   const double& mu() const { return mu_; }
 
-  // Tracks the time of a particle as it traverses the problem.
-  // Units are seconds.
-  double& time() { return time_; }
-  const double& time() const { return time_; }
-  double& time_last() { return time_last_; }
-  const double& time_last() const { return time_last_; }
-
   // Particle lifetime
   double& lifetime() { return lifetime_; }
   const double& lifetime() const { return lifetime_; }
diff --git a/include/openmc/plot.h b/include/openmc/plot.h
index 69e801c5fc3..8073a51e83e 100644
--- a/include/openmc/plot.h
+++ b/include/openmc/plot.h
@@ -171,6 +171,7 @@ class SlicePlotBase {
   Position origin_;           //!< Plot origin in geometry
   Position width_;            //!< Plot width in geometry
   PlotBasis basis_;           //!< Plot basis (XY/XZ/YZ)
+  double time_;               //!< Plot time
   array<size_t, 3> pixels_;   //!< Plot size in pixels
   bool slice_color_overlaps_; //!< Show overlapping cells?
   int slice_level_ {-1};      //!< Plot universe level
@@ -223,6 +224,7 @@ T SlicePlotBase::get_map() const
     GeometryState p;
     p.r() = xyz;
     p.u() = dir;
+    p.time() = time_;
     p.coord(0).universe = model::root_universe;
     int level = slice_level_;
     int j {};
diff --git a/include/openmc/source.h b/include/openmc/source.h
index d58195f1089..12099ff25a3 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -91,7 +91,7 @@ class Source {
 
   // Methods for constraints
   void read_constraints(pugi::xml_node node);
-  bool satisfies_spatial_constraints(Position r) const;
+  bool satisfies_spatial_constraints(Position r, double time) const;
   bool satisfies_energy_constraints(double E) const;
   bool satisfies_time_constraints(double time) const;
 
diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 549e8f6abac..a9510f28c73 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -39,6 +39,12 @@ class Surface {
   std::string name_;                 //!< User-defined name
   unique_ptr<BoundaryCondition> bc_; //!< Boundary condition
   bool surf_source_ {false}; //!< Activate source banking for the surface?
+  
+  //!< Moving surface parameters
+  bool moving_ {false};
+  vector<double> moving_time_grid_;      // Time grid points [0.0, ..., INFTY]
+  vector<Position> moving_translations_; // Translations at the time grid points
+  vector<Position> moving_velocities_;   // Velocities within time bins
 
   explicit Surface(pugi::xml_node surf_node);
   Surface();
@@ -49,9 +55,11 @@ class Surface {
   //! \param r The 3D Cartesian coordinate of a point.
   //! \param u A direction used to "break ties" and pick a sense when the
   //!   point is very close to the surface.
+  //! \param t The time for the evaluation
+  //! \param speed The point speed to "break ties" with moving surface.
   //! \return true if the point is on the "positive" side of the surface and
   //!   false otherwise.
-  bool sense(Position r, Direction u) const;
+  bool sense(Position r, Direction u, double t, double speed) const;
 
   //! Determine the direction of a ray reflected from the surface.
   //! \param[in] r The point at which the ray is incident.
@@ -66,23 +74,37 @@ class Surface {
 
   //! Evaluate the equation describing the surface.
   //!
-  //! Surfaces can be described by some function f(x, y, z) = 0.  This member
-  //! function evaluates that mathematical function.
+  //! Surfaces can be described by some function f(x, y, z) = 0. This member
+  //! function evaluates that mathematical function. The time variable t is
+  //! needed for evaluation of moving surfaces.
   //! \param r A 3D Cartesian coordinate.
-  virtual double evaluate(Position r) const = 0;
+  //! \param t The time for the evaluation.
+  virtual double evaluate(Position r, double t) const;
 
   //! Compute the distance between a point and the surface along a ray.
   //! \param r A 3D Cartesian coordinate.
   //! \param u The direction of the ray.
+  //! \param t The time of the moving point (for moving surface).
+  //! \param speed The speed of moving point (for moving surface).
   //! \param coincident A hint to the code that the given point should lie
   //!   exactly on the surface.
-  virtual double distance(Position r, Direction u, bool coincident) const = 0;
+  virtual double distance(Position r, Direction u, double t, double speed, bool coincident) const;
 
   //! Compute the local outward normal direction of the surface.
   //! \param r A 3D Cartesian coordinate.
   //! \return Normal direction
   virtual Direction normal(Position r) const = 0;
 
+  //! Compute the dot product of the local outward normal direction of the
+  //! surface to a geometry coordinate.
+  //! \param r A 3D Cartesian coordinate.
+  //! \param u The direction of the ray.
+  //! \param t The time for the evaluation.
+  //! \param speed The speed of the particle.
+  //! \return The dot product
+  virtual double dot_normal(
+    Position r, Direction u, double t, double speed) const;
+
   //! Write all information needed to reconstruct the surface to an HDF5 group.
   //! \param group_id An HDF5 group id.
   void to_hdf5(hid_t group_id) const;
@@ -100,6 +122,10 @@ class Surface {
   virtual GeometryType geom_type() const { return GeometryType::CSG; }
 
 protected:
+  //! Static surface functions
+  virtual double _evaluate(Position r) const = 0;
+  virtual double _distance(Position r, Direction u, bool coincident) const = 0;
+
   virtual void to_hdf5_inner(hid_t group_id) const = 0;
 };
 
@@ -112,13 +138,15 @@ class Surface {
 class SurfaceXPlane : public Surface {
 public:
   explicit SurfaceXPlane(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
   BoundingBox bounding_box(bool pos_side) const override;
 
   double x0_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -130,13 +158,15 @@ class SurfaceXPlane : public Surface {
 class SurfaceYPlane : public Surface {
 public:
   explicit SurfaceYPlane(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
   BoundingBox bounding_box(bool pos_side) const override;
 
   double y0_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -148,13 +178,15 @@ class SurfaceYPlane : public Surface {
 class SurfaceZPlane : public Surface {
 public:
   explicit SurfaceZPlane(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
   BoundingBox bounding_box(bool pos_side) const override;
 
   double z0_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -166,12 +198,14 @@ class SurfaceZPlane : public Surface {
 class SurfacePlane : public Surface {
 public:
   explicit SurfacePlane(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   double A_, B_, C_, D_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -184,13 +218,15 @@ class SurfacePlane : public Surface {
 class SurfaceXCylinder : public Surface {
 public:
   explicit SurfaceXCylinder(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
   BoundingBox bounding_box(bool pos_side) const override;
 
   double y0_, z0_, radius_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -203,13 +239,15 @@ class SurfaceXCylinder : public Surface {
 class SurfaceYCylinder : public Surface {
 public:
   explicit SurfaceYCylinder(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
   BoundingBox bounding_box(bool pos_side) const override;
 
   double x0_, z0_, radius_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -222,13 +260,15 @@ class SurfaceYCylinder : public Surface {
 class SurfaceZCylinder : public Surface {
 public:
   explicit SurfaceZCylinder(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
   BoundingBox bounding_box(bool pos_side) const override;
 
   double x0_, y0_, radius_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -241,13 +281,15 @@ class SurfaceZCylinder : public Surface {
 class SurfaceSphere : public Surface {
 public:
   explicit SurfaceSphere(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
   BoundingBox bounding_box(bool pos_side) const override;
 
   double x0_, y0_, z0_, radius_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -260,12 +302,14 @@ class SurfaceSphere : public Surface {
 class SurfaceXCone : public Surface {
 public:
   explicit SurfaceXCone(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   double x0_, y0_, z0_, radius_sq_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -278,12 +322,14 @@ class SurfaceXCone : public Surface {
 class SurfaceYCone : public Surface {
 public:
   explicit SurfaceYCone(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   double x0_, y0_, z0_, radius_sq_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -296,12 +342,14 @@ class SurfaceYCone : public Surface {
 class SurfaceZCone : public Surface {
 public:
   explicit SurfaceZCone(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   double x0_, y0_, z0_, radius_sq_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -314,13 +362,15 @@ class SurfaceZCone : public Surface {
 class SurfaceQuadric : public Surface {
 public:
   explicit SurfaceQuadric(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
   double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -332,12 +382,14 @@ class SurfaceQuadric : public Surface {
 class SurfaceXTorus : public Surface {
 public:
   explicit SurfaceXTorus(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   double x0_, y0_, z0_, A_, B_, C_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -349,12 +401,14 @@ class SurfaceXTorus : public Surface {
 class SurfaceYTorus : public Surface {
 public:
   explicit SurfaceYTorus(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   double x0_, y0_, z0_, A_, B_, C_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
@@ -366,12 +420,14 @@ class SurfaceYTorus : public Surface {
 class SurfaceZTorus : public Surface {
 public:
   explicit SurfaceZTorus(pugi::xml_node surf_node);
-  double evaluate(Position r) const override;
-  double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
   void to_hdf5_inner(hid_t group_id) const override;
 
   double x0_, y0_, z0_, A_, B_, C_;
+
+protected:
+  double _evaluate(Position r) const override;
+  double _distance(Position r, Direction u, bool coincident) const override;
 };
 
 //==============================================================================
diff --git a/include/openmc/universe.h b/include/openmc/universe.h
index e8fbacfdc7b..15910b4f39c 100644
--- a/include/openmc/universe.h
+++ b/include/openmc/universe.h
@@ -60,7 +60,8 @@ class UniversePartitioner {
   explicit UniversePartitioner(const Universe& univ);
 
   //! Return the list of cells that could contain the given coordinates.
-  const vector<int32_t>& get_cells(Position r, Direction u) const;
+  const vector<int32_t>& get_cells(
+    Position r, Direction u, double t, double speed) const;
 
 private:
   //! A sorted vector of indices to surfaces that partition the universe
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 55f309bfbd3..566d8562f04 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -909,6 +909,7 @@ def plot(
         width: Sequence[float] | None = None,
         pixels: int | Sequence[int] = 40000,
         basis: str = 'xy',
+        time: float = 0.0,
         color_by: str = 'cell',
         colors: dict | None = None,
         seed: int | None = None,
@@ -1005,6 +1006,7 @@ def plot(
             plot.width = width
             plot.pixels = pixels
             plot.basis = basis
+            plot.time = time
             plot.color_by = color_by
             plot.show_overlaps = show_overlaps
             if overlap_color is not None:
@@ -1105,7 +1107,6 @@ def plot(
             if outline != 'only':
                 axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs)
 
-
         if n_samples:
             # Sample external source particles
             particles = self.sample_external_source(n_samples)
diff --git a/openmc/plots.py b/openmc/plots.py
index 40fcf0c78da..405776c8206 100644
--- a/openmc/plots.py
+++ b/openmc/plots.py
@@ -185,6 +185,8 @@
             the image aspect ratio.
         basis : {'xy', 'xz', 'yz'}
             The basis directions for the plot
+        time : float
+            Simulation time at which geometry is evaluated
         color_by : {'cell', 'material'}
             Indicate whether the plot should be colored by cell or by material
         colors : dict
@@ -673,6 +675,8 @@ class Plot(PlotBase):
         The type of the plot
     basis : {'xy', 'xz', 'yz'}
         The basis directions for the plot
+    time : float
+        Time at which geometry is evaluated
     meshlines : dict
         Dictionary defining type, id, linewidth and color of a mesh to be
         plotted on top of a plot
@@ -685,6 +689,7 @@ def __init__(self, plot_id=None, name=''):
         self._origin = [0., 0., 0.]
         self._type = 'slice'
         self._basis = 'xy'
+        self._time = 0.0
         self._meshlines = None
 
     @property
@@ -725,6 +730,15 @@ def basis(self, basis):
         cv.check_value('plot basis', basis, _BASES)
         self._basis = basis
 
+    @property
+    def time(self):
+        return self._time
+
+    @time.setter
+    def time(self, time):
+        cv.check_type('plot time', time, Real)
+        self._time = time
+
     @property
     def meshlines(self):
         return self._meshlines
@@ -766,6 +780,7 @@ def __repr__(self):
         string += '{: <16}=\t{}\n'.format('\tBasis', self._basis)
         string += '{: <16}=\t{}\n'.format('\tWidth', self._width)
         string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin)
+        string += '{: <16}=\t{}\n'.format('\tTime', self._time)
         string += '{: <16}=\t{}\n'.format('\tPixels', self._pixels)
         string += '{: <16}=\t{}\n'.format('\tColor by', self._color_by)
         string += '{: <16}=\t{}\n'.format('\tBackground', self._background)
@@ -899,6 +914,9 @@ def to_xml_element(self):
         subelement = ET.SubElement(element, "width")
         subelement.text = ' '.join(map(str, self._width))
 
+        subelement = ET.SubElement(element, "time")
+        subelement.text = str(self._time)
+
         if self._colors:
             self._colors_to_xml(element)
 
diff --git a/openmc/surface.py b/openmc/surface.py
index 840c3125d88..a4587efd8fc 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -148,7 +148,13 @@ class Surface(IDManagerMixin, ABC):
         Name of the surface
     type : str
         Type of the surface
-
+    moving : bool
+        Whether or not the surface is moving
+    moving_velocities : iterable of iterable of float
+        3D velocities in which the surface moves
+    moving_durations : iterable of float
+        Durations during which the surface moves with the associated
+        velocities.
     """
 
     next_id = 1
@@ -167,6 +173,10 @@ def __init__(self, surface_id=None, boundary_type='transmission',
         # Value    - coefficient value
         self._coefficients = {}
 
+        self._moving = False
+        self._moving_velocities = np.array([[]])
+        self._moving_durations = np.array([])
+
     def __neg__(self):
         return Halfspace(self, '-')
 
@@ -192,6 +202,24 @@ def __repr__(self):
 
         string += coefficients
 
+        # TODO: print moving durations and velocities
+        if self._moving:
+            '''
+            string += '{0: <20}{1}{2}\n'.format(
+                'moving velocities', '=\t', np.array_repr(self._moving_velocities).replace('\n', ' '))
+            string += '{0: <20}{1}{2}\n'.format('moving durations',
+                                                '=\t', self._moving_durations)
+            '''
+            string += "\tMove parameters\n"
+            N_move = len(self._moving_durations)
+            string += "T        Vx       Vy       Vz\n"
+            for i in range(N_move):
+                string += '{0: <8} {1: <8} {2: <8} {3: <8}\n'.format(
+                    self._moving_durations[i],
+                    self._moving_velocities[i,0],
+                    self._moving_velocities[i,1],
+                    self._moving_velocities[i,2])
+
         return string
 
     @property
@@ -234,6 +262,18 @@ def albedo(self, albedo):
     def coefficients(self):
         return self._coefficients
 
+    @property
+    def moving(self):
+        return self._moving
+
+    @property
+    def moving_velicities(self):
+        return self._moving_velocities
+
+    @property
+    def moving_durations(self):
+        return self._moving_durations
+
     def bounding_box(self, side):
         """Determine an axis-aligned bounding box.
 
@@ -327,6 +367,68 @@ def is_equal(self, other):
 
         return np.allclose(coeffs1, coeffs2, rtol=0., atol=self._atol)
 
+    def move(self, velocities, durations):
+        """Set the surface continuous movement based on the given velocities
+        and durations.
+
+        Parameters
+        ==========
+        velocities : iterable of iterable of float
+            3D velocities in which the surface moves
+        durations : iterable of float
+            Durations during which the surface moves with the associated
+            velocities.
+        """
+        self._moving = True
+        self._moving_velocities = np.array(velocities)
+        self._moving_durations = np.array(durations)
+
+        # Add the final static condition
+        self._moving_velocities = np.append(self._moving_velocities, np.zeros((1,3)), axis=0)
+        self._moving_durations = np.append(self._moving_durations, np.inf)
+
+        # Set moving time grid and translations for evaluation convenience
+        N_move = len(self._moving_durations)
+        self._moving_time_grid = np.zeros(N_move + 1)
+        self._moving_translations = np.zeros((N_move + 1, 3))
+        for n in range(N_move - 1):
+            duration = self._moving_durations[n]
+            velocity = self._moving_velocities[n]
+
+            t_start = self._moving_time_grid[n]
+            self._moving_time_grid[n + 1] = t_start + duration
+
+            trans_start = self._moving_translations[n]
+            self._moving_translations[n + 1] = trans_start + velocity * duration
+        # Manual assignment to avoid Numpy's warning
+        self._moving_time_grid[N_move] = np.inf
+        self._moving_translations[N_move] = self._moving_translations[N_move - 1]
+
+
+    def _adjust_point(self, point, time):
+        """Adjust point based on the surface movement"""
+        if not self._moving:
+            return point
+
+        # Get moving index
+        time_grid = self._moving_time_grid
+        idx = np.searchsorted(time_grid, time, side='right') - 1
+
+        # Get moving translation, velocity, and starting time
+        trans_0 = self._moving_translations[idx]
+        time_0 = self._moving_time_grid[idx]
+        V = self._moving_velocities[idx]
+
+        # Translate the particle
+        t_local = time - time_0
+        x = point[0] - (trans_0[0] + V[0] * t_local)
+        y = point[1] - (trans_0[1] + V[1] * t_local)
+        z = point[2] - (trans_0[2] + V[2] * t_local)
+
+        print(time_grid)
+        print(time, idx)
+        return (x, y, z)
+
     @abstractmethod
     def _get_base_coeffs(self):
         """Return polynomial coefficients representing the implicit surface
@@ -335,19 +437,22 @@ def _get_base_coeffs(self):
         """
 
     @abstractmethod
-    def evaluate(self, point):
-        """Evaluate the surface equation at a given point.
+    def evaluate(self, point, time=0.0):
+        """Evaluate the surface equation at a given point and time.
 
         Parameters
         ----------
         point : 3-tuple of float
             The Cartesian coordinates, :math:`(x',y',z')`, at which the surface
             equation should be evaluated.
+        time : float, optional
+            The time :math:`t'` at which the surface equation should be evaluated
 
         Returns
         -------
         float
             Evaluation of the surface polynomial at point :math:`(x',y',z')`
+            and time :math:`t'`
 
         """
 
@@ -432,6 +537,12 @@ def to_xml_element(self):
         element.set("coeffs", ' '.join([str(self._coefficients.setdefault(key, 0.0))
                                         for key in self._coeff_keys]))
 
+        if self._moving:
+            element.set("moving_velocities", ' '.join(
+                [str(value) for value in self._moving_velocities[:-1]]))
+            element.set("moving_durations", ' '.join(
+                [str(value) for value in self._moving_durations[:-1]]))
+
         return element
 
     @staticmethod
@@ -574,14 +685,16 @@ def bounding_box(self, side):
 
         return BoundingBox(ll, ur)
 
-    def evaluate(self, point):
-        """Evaluate the surface equation at a given point.
+    def evaluate(self, point, time=0.0):
+        """Evaluate the surface equation at a given point and time.
 
         Parameters
         ----------
         point : 3-tuple of float
             The Cartesian coordinates, :math:`(x',y',z')`, at which the surface
             equation should be evaluated.
+        time : float, optional
+            The time :math:`t'` at which the surface equation should be evaluated
 
         Returns
         -------
@@ -590,7 +703,7 @@ def evaluate(self, point):
 
         """
 
-        x, y, z = point
+        x, y, z = self._adjust_point(point, time)
         a, b, c, d = self._get_base_coeffs()
         return a*x + b*y + c*z - d
 
@@ -872,8 +985,9 @@ def __init__(self, x0=0., *args, **kwargs):
     c = SurfaceCoefficient(0.)
     d = x0
 
-    def evaluate(self, point):
-        return point[0] - self.x0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        return x - self.x0
 
 
 class YPlane(PlaneMixin, Surface):
@@ -937,8 +1051,9 @@ def __init__(self, y0=0., *args, **kwargs):
     c = SurfaceCoefficient(0.)
     d = y0
 
-    def evaluate(self, point):
-        return point[1] - self.y0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        return y - self.y0
 
 
 class ZPlane(PlaneMixin, Surface):
@@ -1002,8 +1117,9 @@ def __init__(self, z0=0., *args, **kwargs):
     c = SurfaceCoefficient(1.)
     d = z0
 
-    def evaluate(self, point):
-        return point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        return z - self.z0
 
 
 class QuadricMixin:
@@ -1060,14 +1176,16 @@ def eigh(self, coeffs=None):
         """
         return np.linalg.eigh(self.get_Abc(coeffs=coeffs)[0])
 
-    def evaluate(self, point):
-        """Evaluate the surface equation at a given point.
+    def evaluate(self, point, time=0.0):
+        """Evaluate the surface equation at a given point and time.
 
         Parameters
         ----------
         point : 3-tuple of float
             The Cartesian coordinates, :math:`(x',y',z')`, in [cm] at which the
             surface equation should be evaluated.
+        time : float, optional
+            The time :math:`t'` at which the surface equation should be evaluated
 
         Returns
         -------
@@ -1076,7 +1194,7 @@ def evaluate(self, point):
             Jz' + K = 0`
 
         """
-        x = np.asarray(point)
+        x = np.asarray(self._adjust_point(point, time))
         A, b, c = self.get_Abc()
         return x.T @ A @ x + b.T @ x + c
 
@@ -1452,9 +1570,10 @@ def bounding_box(self, side):
         elif side == '+':
             return BoundingBox.infinite()
 
-    def evaluate(self, point):
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        y -= self.y0
+        z -= self.z0
         return y*y + z*z - self.r**2
 
 
@@ -1550,9 +1669,10 @@ def bounding_box(self, side):
         elif side == '+':
             return BoundingBox.infinite()
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        z -= self.z0
         return x*x + z*z - self.r**2
 
 
@@ -1648,9 +1768,10 @@ def bounding_box(self, side):
         elif side == '+':
             return BoundingBox.infinite()
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
         return x*x + y*y - self.r**2
 
 
@@ -1745,10 +1866,11 @@ def bounding_box(self, side):
         elif side == '+':
             return BoundingBox.infinite()
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
+        z -= self.z0
         return x*x + y*y + z*z - self.r**2
 
 
@@ -2003,10 +2125,11 @@ def _get_base_coeffs(self):
 
         return (a, b, c, d, e, f, g, h, j, k)
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
+        z -= self.z0
         return y*y + z*z - self.r2*x*x
 
 
@@ -2105,10 +2228,11 @@ def _get_base_coeffs(self):
 
         return (a, b, c, d, e, f, g, h, j, k)
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
+        z -= self.z0
         return x*x + z*z - self.r2*y*y
 
 
@@ -2207,10 +2331,11 @@ def _get_base_coeffs(self):
 
         return (a, b, c, d, e, f, g, h, j, k)
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
+        z -= self.z0
         return x*x + y*y - self.r2*z*z
 
 
@@ -2409,10 +2534,11 @@ class XTorus(TorusMixin, Surface):
     """
     _type = 'x-torus'
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
+        z -= self.z0
         a = self.a
         b = self.b
         c = self.c
@@ -2484,10 +2610,11 @@ class YTorus(TorusMixin, Surface):
     """
     _type = 'y-torus'
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
+        z -= self.z0
         a = self.a
         b = self.b
         c = self.c
@@ -2559,10 +2686,11 @@ class ZTorus(TorusMixin, Surface):
 
     _type = 'z-torus'
 
-    def evaluate(self, point):
-        x = point[0] - self.x0
-        y = point[1] - self.y0
-        z = point[2] - self.z0
+    def evaluate(self, point, time=0.0):
+        x, y, z = self._adjust_point(point, time)
+        x -= self.x0
+        y -= self.y0
+        z -= self.z0
         a = self.a
         b = self.b
         c = self.c
diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp
index 7216ac89649..a2fea6ea578 100644
--- a/src/boundary_condition.cpp
+++ b/src/boundary_condition.cpp
@@ -101,22 +101,22 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf)
   Position origin {0, 0, 0};
   Direction u = surf1.normal(origin);
   double d1;
-  double e1 = surf1.evaluate(origin);
+  double e1 = surf1.evaluate(origin, 0.0); // Time = 0.0
   if (e1 > FP_COINCIDENT) {
-    d1 = -surf1.distance(origin, -u, false);
+    d1 = -surf1.distance(origin, -u, 0.0, 0.0, false); // Time = speed = 0.0
   } else if (e1 < -FP_COINCIDENT) {
-    d1 = surf1.distance(origin, u, false);
+    d1 = surf1.distance(origin, u, 0.0, 0.0, false); // Time = speed = 0.0
   } else {
     d1 = 0.0;
   }
 
   // Compute the distance from the second surface to the origin.
   double d2;
-  double e2 = surf2.evaluate(origin);
+  double e2 = surf2.evaluate(origin, 0.0); // Time = 0.0
   if (e2 > FP_COINCIDENT) {
-    d2 = -surf2.distance(origin, -u, false);
+    d2 = -surf2.distance(origin, -u, 0.0, 0.0, false); // Time = speed = 0.0
   } else if (e2 < -FP_COINCIDENT) {
-    d2 = surf2.distance(origin, u, false);
+    d2 = surf2.distance(origin, u, 0.0, 0.0, false); // Time = speed = 0.0
   } else {
     d2 = 0.0;
   }
@@ -216,14 +216,14 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
   }
 
   // Make sure both surfaces intersect the origin
-  if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) {
+  if (std::abs(surf1.evaluate({0, 0, 0}, 0.0)) > FP_COINCIDENT) { // Time = 0.0
     throw std::invalid_argument(fmt::format(
       "Rotational periodic BCs are only "
       "supported for rotations about the origin, but surface {} does not "
       "intersect the origin.",
       surf1.id_));
   }
-  if (std::abs(surf2.evaluate({0, 0, 0})) > FP_COINCIDENT) {
+  if (std::abs(surf2.evaluate({0, 0, 0}, 0.0)) > FP_COINCIDENT) { // Time = 0.0
     throw std::invalid_argument(fmt::format(
       "Rotational periodic BCs are only "
       "supported for rotations about the origin, but surface {} does not "
diff --git a/src/cell.cpp b/src/cell.cpp
index 4b26992299e..5c6c6cec35c 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -775,7 +775,7 @@ std::string Region::str() const
 //==============================================================================
 
 std::pair<double, int32_t> Region::distance(
-  Position r, Direction u, int32_t on_surface) const
+  Position r, Direction u, double t, double speed, int32_t on_surface) const
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
@@ -788,7 +788,7 @@ std::pair<double, int32_t> Region::distance(
     // Calculate the distance to this surface.
     // Note the off-by-one indexing
     bool coincident {std::abs(token) == std::abs(on_surface)};
-    double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)};
+    double d {model::surfaces[abs(token) - 1]->distance(r, u, t, speed, coincident)};
 
     // Check if this distance is the new minimum.
     if (d < min_dist) {
@@ -804,30 +804,45 @@ std::pair<double, int32_t> Region::distance(
 
 //==============================================================================
 
-bool Region::contains(Position r, Direction u, int32_t on_surface) const
+bool Region::contains(
+  Position r, Direction u, double t, double speed, int32_t on_surface) const
 {
   if (simple_) {
-    return contains_simple(r, u, on_surface);
+    return contains_simple(r, u, t, speed, on_surface);
   } else {
-    return contains_complex(r, u, on_surface);
+    return contains_complex(r, u, t, speed, on_surface);
   }
 }
 
 //==============================================================================
 
-bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const
+bool Region::contains_simple(
+  Position r, Direction u, double t, double speed, int32_t on_surface) const
 {
   for (int32_t token : expression_) {
     // Assume that no tokens are operators. Evaluate the sense of particle with
     // respect to the surface and see if the token matches the sense. If the
     // particle's surface attribute is set and matches the token, that
-    // overrides the determination based on sense().
+    // overrides the determination based on sense(). However, if the surface
+    // is moving, we need to factor in the surface-particle relative velocity,
+    // so the surface-sense calculation needs to be done.
+   
+    if (model::surfaces[abs(token) - 1]->moving_) {
+      // Note the off-by-one indexing
+      bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed);
+      if (sense != (token > 0)) {
+        return false;
+      }
+      continue;
+    }
+    // Surface is not moving
+    
     if (token == on_surface) {
     } else if (-token == on_surface) {
       return false;
     } else {
       // Note the off-by-one indexing
-      bool sense = model::surfaces[abs(token) - 1]->sense(r, u);
+      bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed);
       if (sense != (token > 0)) {
         return false;
       }
@@ -838,7 +853,8 @@ bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const
 
 //==============================================================================
 
-bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const
+bool Region::contains_complex(
+  Position r, Direction u, double t, double speed, int32_t on_surface) const
 {
   bool in_cell = true;
   int total_depth = 0;
@@ -851,14 +867,21 @@ bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const
     // If the token is a union or intersection check to
     // short circuit
     if (token < OP_UNION) {
-      if (token == on_surface) {
-        in_cell = true;
-      } else if (-token == on_surface) {
-        in_cell = false;
-      } else {
+      if (model::surfaces[abs(token) - 1]->moving_) {
         // Note the off-by-one indexing
-        bool sense = model::surfaces[abs(token) - 1]->sense(r, u);
+        bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed);
         in_cell = (sense == (token > 0));
+      } else {
+        // Surface is not moving
+        if (token == on_surface) {
+          in_cell = true;
+        } else if (-token == on_surface) {
+          in_cell = false;
+        } else {
+          // Note the off-by-one indexing
+          bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed);
+          in_cell = (sense == (token > 0));
+        }
       }
     } else if ((token == OP_UNION && in_cell == true) ||
                (token == OP_INTERSECTION && in_cell == false)) {
diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index 13436088652..879c1a8a241 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -769,7 +769,7 @@ moab::EntityHandle DAGSurface::mesh_handle() const
   return dagmc_ptr()->entity_by_index(2, dag_index());
 }
 
-double DAGSurface::evaluate(Position r) const
+double DAGSurface::evaluate(Position r, double t) const
 {
   return 0.0;
 }
@@ -800,6 +800,12 @@ Direction DAGSurface::normal(Position r) const
   return dir;
 }
 
+double DAGSurface::dot_normal(
+  Position r, Direction u, double t, double speed) const
+{
+  return u.dot(normal(r));
+}
+
 Direction DAGSurface::reflect(Position r, Direction u, GeometryState* p) const
 {
   assert(p);
diff --git a/src/geometry.cpp b/src/geometry.cpp
index 9e975c00dee..d7d76fe7b1f 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -43,7 +43,8 @@ bool check_cell_overlap(GeometryState& p, bool error)
     // Loop through each cell on this level
     for (auto index_cell : univ.cells_) {
       Cell& c = *model::cells[index_cell];
-      if (c.contains(p.coord(j).r, p.coord(j).u, p.surface())) {
+      if (c.contains(
+            p.coord(j).r, p.coord(j).u, p.time(), p.surface(), p.speed())) {
         if (index_cell != p.coord(j).cell) {
           if (error) {
             fatal_error(
@@ -110,7 +111,7 @@ bool find_cell_inner(
     for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) {
       i_cell = *it;
 
-      // Make sure the search cell is in the same universe.
+      // Make sure the search cell is in the same `universe.
       int i_universe = p.lowest_coord().universe;
       if (model::cells[i_cell]->universe_ != i_universe)
         continue;
@@ -118,8 +119,10 @@ bool find_cell_inner(
       // Check if this cell contains the particle.
       Position r {p.r_local()};
       Direction u {p.u_local()};
+      double t {p.time()};
+      double speed {p.speed()};
       auto surf = p.surface();
-      if (model::cells[i_cell]->contains(r, u, surf)) {
+      if (model::cells[i_cell]->contains(r, u, t, speed, surf)) {
         p.lowest_coord().cell = i_cell;
         found = true;
         break;
@@ -369,10 +372,12 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
     const auto& coord {p.coord(i)};
     const Position& r {coord.r};
     const Direction& u {coord.u};
+    const double t {p.time()};
+    const double speed {p.speed()};
     Cell& c {*model::cells[coord.cell]};
 
     // Find the oncoming surface in this cell and the distance to it.
-    auto surface_distance = c.distance(r, u, p.surface(), &p);
+    auto surface_distance = c.distance(r, u, t, speed, p.surface(), &p);
     d_surf = surface_distance.first;
     level_surf_cross = surface_distance.second;
 
diff --git a/src/particle.cpp b/src/particle.cpp
index 6cfdf2175af..d1bec35239d 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -42,7 +42,7 @@ namespace openmc {
 // Particle implementation
 //==============================================================================
 
-double Particle::speed() const
+double Particle::get_speed() const
 {
   // Determine mass in eV/c^2
   double mass;
@@ -168,6 +168,15 @@ void Particle::event_calculate_xs()
   r_last() = r();
   time_last() = time();
 
+  // Calculate speed.
+  // In MG mode, we can't determine particle speed if cell, and thus material, 
+  // is not known. Just set to zero for now.
+  if (settings::run_CE || !(lowest_coord().cell == C_NONE)) {
+    this->speed() = this->get_speed();
+  } else {
+    this->speed() = 0.0;
+  }
+
   // Reset event variables
   event() = TallyEvent::KILL;
   event_nuclide() = NUCLIDE_NONE;
@@ -192,6 +201,11 @@ void Particle::event_calculate_xs()
       cell_last(j) = coord(j).cell;
     }
     n_coord_last() = n_coord();
+    
+    // Calculate speed for MG mode
+    if (!settings::run_CE) {
+      this->speed() = this->get_speed();
+    }
   }
 
   // Write particle track.
diff --git a/src/particle_data.cpp b/src/particle_data.cpp
index fef8359a356..e81debf82ca 100644
--- a/src/particle_data.cpp
+++ b/src/particle_data.cpp
@@ -64,7 +64,7 @@ void GeometryState::advance_to_boundary_from_void()
 
   for (auto c_i : root_universe->cells_) {
     auto dist =
-      model::cells.at(c_i)->distance(root_coord.r, root_coord.u, 0, this);
+      model::cells.at(c_i)->distance(root_coord.r, root_coord.u, this->time(), this->speed(), 0, this);
     if (dist.first < boundary().distance) {
       boundary().distance = dist.first;
       boundary().surface = dist.second;
diff --git a/src/plot.cpp b/src/plot.cpp
index dbc25b21e4d..a8a18a6fdbb 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -145,6 +145,7 @@ void Plot::print_info() const
 
   // Plot parameters
   fmt::print("Origin: {} {} {}\n", origin_[0], origin_[1], origin_[2]);
+  fmt::print("Time: {}\n", time_);
 
   if (PlotType::slice == type_) {
     fmt::print("Width: {:4} {:4}\n", width_[0], width_[1]);
@@ -721,6 +722,7 @@ Plot::Plot(pugi::xml_node plot_node, PlotType type)
 {
   set_output_path(plot_node);
   set_basis(plot_node);
+  time_ = get_node_array<double>(plot_node, "time")[0];
   set_origin(plot_node);
   set_width(plot_node);
   set_meshlines(plot_node);
diff --git a/src/simulation.cpp b/src/simulation.cpp
index d9a4fd131dc..a919fe597f2 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -580,6 +580,16 @@ void initialize_history(Particle& p, int64_t index_source)
   }
   p.current_work() = index_source;
 
+  // Speed will be needed if the particle happens to be coincident on a moving
+  // surface. In MG mode, we get a chicken-egg situation as we need material to 
+  // get speed, but we need speed to determine cell and thus material. 
+  // Put zero for now.
+  if (settings::run_CE) {
+    p.speed() = p.get_speed();
+  } else {
+    p.speed() = 0.0;
+  }
+
   // set identifier for particle
   p.id() = simulation::work_index[mpi::rank] + index_source;
 
diff --git a/src/source.cpp b/src/source.cpp
index 65873f5a797..c70ca122849 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -189,7 +189,7 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const
       accepted = true;
     } else {
       // Check whether sampled site satisfies constraints
-      accepted = satisfies_spatial_constraints(site.r) &&
+      accepted = satisfies_spatial_constraints(site.r, site.time) &&
                  satisfies_energy_constraints(site.E) &&
                  satisfies_time_constraints(site.time);
       if (!accepted) {
@@ -223,11 +223,12 @@ bool Source::satisfies_time_constraints(double time) const
   return time > time_bounds_.first && time < time_bounds_.second;
 }
 
-bool Source::satisfies_spatial_constraints(Position r) const
+bool Source::satisfies_spatial_constraints(Position r, double time) const
 {
   GeometryState geom_state;
   geom_state.r() = r;
   geom_state.u() = {0.0, 0.0, 1.0};
+  geom_state.time() = time;
 
   // Reject particle if it's not in the geometry at all
   bool found = exhaustive_find_cell(geom_state);
@@ -583,9 +584,30 @@ SourceSite MeshSource::sample(uint64_t* seed) const
   // Sample a mesh element based on the relative strengths
   int32_t element = space_->sample_element_index(seed);
 
-  // Sample the distribution for the specific mesh element; note that the
-  // spatial distribution has been set for each element using MeshElementSpatial
-  return source(element)->sample_with_constraints(seed);
+  SourceSite site;
+  while (true) {
+    // Sample the distribution for the specific mesh element; note that the
+    // spatial distribution has been set for each element using MeshElementSpatial
+    site = source(element)->sample_with_constraints(seed);
+   
+    // Reject time?
+    if (!satisfies_time_constraints(site.time)) { continue; }
+
+    // Sample position and apply rejection on spatial domains
+    Position r;
+    do {
+      r = space_->mesh()->sample_element(element, seed);
+    } while (!this->satisfies_spatial_constraints(r, site.time));
+    // Replace the position
+    site.r = r;
+    
+    // Apply other rejections
+    if (satisfies_energy_constraints(site.E)) {
+      break;
+    }
+  }
+
+  return site;
 }
 
 //==============================================================================
diff --git a/src/surface.cpp b/src/surface.cpp
index ea19514a311..c4d54e6404d 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -15,6 +15,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
+#include "openmc/search.h"
 #include "openmc/settings.h"
 #include "openmc/string_utils.h"
 #include "openmc/xml_interface.h"
@@ -111,20 +112,76 @@ Surface::Surface(pugi::xml_node surf_node)
       bc_->set_albedo(surf_alb);
     }
   }
+  
+  // Not moving?
+  if (!(check_for_node(surf_node, "moving_velocities") ||
+        check_for_node(surf_node, "moving_durations"))) {
+    moving_ = false;
+    return;
+  }
+
+  // Now, set the surface moving parameters
+  moving_ = true;
+
+  // Moving durations
+  auto durations = get_node_array<double>(surf_node, "moving_durations");
+  const int N_move = durations.size() + 1;
+
+  // Moving time grids
+  moving_time_grid_.resize(N_move + 1);
+  moving_time_grid_[0] = 0.0;
+  for (int n = 0; n < N_move - 1; n++) {
+    moving_time_grid_[n + 1] = moving_time_grid_[n] + durations[n];
+  }
+  moving_time_grid_[N_move] = INFTY;
+
+  // Moving velocities
+  moving_velocities_.resize(N_move);
+  std::string velocities_spec = get_node_value(surf_node, "moving_velocities");
+  // Parse
+  std::vector<double> numbers;
+  for (int i = 0; i < velocities_spec.size();) {
+    if (velocities_spec[i] == '-' || std::isdigit(velocities_spec[i])) {
+      int j = i + 1;
+      while (j < velocities_spec.size() && (std::isdigit(velocities_spec[j]) || velocities_spec[j] == '.')) {
+        j++;
+      }
+      numbers.push_back(std::stod(velocities_spec.substr(i, j - i)));
+      i = j;
+    }
+    i++;
+  }
+  // Assign to velocities
+  for (int n = 0; n < N_move - 1; n++) {
+    int idx = 3 * n;
+    moving_velocities_[n][0] = numbers[idx];
+    moving_velocities_[n][1] = numbers[idx + 1];
+    moving_velocities_[n][2] = numbers[idx + 2];
+  }
+  moving_velocities_[N_move - 1] *= 0.0;
+
+  // Moving translations
+  moving_translations_.resize(N_move + 1);
+  moving_translations_[0] *= 0.0;
+  for (int n = 0; n < N_move - 1; n++) {
+    moving_translations_[n + 1] =
+      moving_translations_[n] + moving_velocities_[n] * durations[n];
+  }
+  moving_translations_[N_move] = moving_translations_[N_move - 1];
 }
 
-bool Surface::sense(Position r, Direction u) const
+bool Surface::sense(Position r, Direction u, double t, double speed) const
 {
   // Evaluate the surface equation at the particle's coordinates to determine
   // which side the particle is on.
-  const double f = evaluate(r);
+  const double f = evaluate(r, t);
 
   // Check which side of surface the point is on.
   if (std::abs(f) < FP_COINCIDENT) {
     // Particle may be coincident with this surface. To determine the sense, we
     // look at the direction of the particle relative to the surface normal (by
     // default in the positive direction) via their dot product.
-    return u.dot(normal(r)) > 0.0;
+    return dot_normal(r, u, t, speed) > 0.0;
   }
   return f > 0.0;
 }
@@ -160,6 +217,133 @@ Direction Surface::diffuse_reflect(
   return u / u.norm();
 }
 
+double Surface::evaluate(Position r, double t) const
+{
+  if (!moving_) {
+    return _evaluate(r);
+  }
+  // The surface moves
+
+  // Get moving index
+  int idx =
+    lower_bound_index(moving_time_grid_.begin(), moving_time_grid_.end(), t);
+
+  // Get moving translation, velocity, and starting time
+  Position translation = moving_translations_[idx];
+  double time_0 = moving_time_grid_[idx];
+  Position velocity = moving_velocities_[idx];
+
+  // Move the position relative to the surface movement
+  double t_local = t - time_0;
+  Position r_moved = r - (translation + velocity * t_local);
+
+  // Evaluate the moved position
+  return _evaluate(r_moved);
+}
+
+double Surface::distance(Position r, Direction u, double t, double speed, bool coincident) const
+{
+  if (!moving_) {
+    return _distance(r, u, coincident);
+  }
+  // The surface moves
+  
+  coincident = false;
+
+  // Store the origin coordinate
+  Position r_origin {r};
+  Position u_origin {u};
+  double t_origin {t};
+
+  // Get moving interval index
+  int idx =
+    lower_bound_index(moving_time_grid_.begin(), moving_time_grid_.end(), t);
+
+  // Distance accumulator
+  double distance_total = 0.0;
+
+  // Evaluate the current and the subsequent intervals until intersecting
+  while (idx < moving_velocities_.size()) {
+    // Get moving translation, velocity, and starting time
+    Position translation = moving_translations_[idx];
+    double time_0 = moving_time_grid_[idx];
+    Position velocity = moving_velocities_[idx];
+
+    // Adjust the position based on the surface movement
+    double t_local = t - time_0;
+    r -= (translation + velocity * t_local);
+
+    // Adjust to relative direction
+    u -= velocity / speed;
+
+    // Normalize scaled relative direction
+    double norm = u.norm();
+    u /= norm;
+        
+    // Get distance using the static function based on the adjusted position
+    // and direction
+    double distance = _distance(r, u, coincident) * norm;
+
+    // Rescale relative direction
+    u *= norm;
+        
+    // Intersection within the interval?
+    double t_distance = distance / speed;
+    double t_interval = moving_time_grid_[idx + 1] - t;
+    if (t_distance < t_interval) {
+      // Return the total distance
+      return distance_total + distance;
+    }
+    // Not intersecting. 
+    // Prepare to check the next interval.
+
+    // Accumulate distance
+    distance_total += t_interval * speed;
+
+    // March forward the coordinate
+    r = r_origin + distance_total * u_origin;
+    u = u_origin;
+    t = moving_time_grid_[idx + 1];
+
+    idx++;
+  }
+
+  // No intersection
+  return INFTY;
+}
+
+double Surface::dot_normal(
+  Position r, Direction u, double t, double speed) const
+{
+  if (!moving_) {
+    return u.dot(normal(r));
+  }
+  // The surface moves
+
+  // Get moving index
+  int idx =
+    lower_bound_index(moving_time_grid_.begin(), moving_time_grid_.end(), t);
+
+  // Get moving translation, velocity, and starting time
+  Position translation = moving_translations_[idx];
+  double time_0 = moving_time_grid_[idx];
+  Position velocity = moving_velocities_[idx];
+
+  // Move the position relative to the surface movement
+  double t_local = t - time_0;
+  Position r_moved = r - (translation + velocity * t_local);
+
+  // Get the relative direction
+  Direction u_relative = u - velocity / speed;
+  
+  // Normalize scaled relative direction
+  double norm = u_relative.norm();
+  u_relative /= norm;
+  
+  // Get the dot product
+  return u_relative.dot(normal(r_moved));
+}
+
 void Surface::to_hdf5(hid_t group_id) const
 {
   hid_t surf_group = create_group(group_id, fmt::format("surface {}", id_));
@@ -226,12 +410,12 @@ SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&x0_});
 }
 
-double SurfaceXPlane::evaluate(Position r) const
+double SurfaceXPlane::_evaluate(Position r) const
 {
   return r.x - x0_;
 }
 
-double SurfaceXPlane::distance(Position r, Direction u, bool coincident) const
+double SurfaceXPlane::_distance(Position r, Direction u, bool coincident) const
 {
   return axis_aligned_plane_distance<0>(r, u, coincident, x0_);
 }
@@ -266,12 +450,12 @@ SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&y0_});
 }
 
-double SurfaceYPlane::evaluate(Position r) const
+double SurfaceYPlane::_evaluate(Position r) const
 {
   return r.y - y0_;
 }
 
-double SurfaceYPlane::distance(Position r, Direction u, bool coincident) const
+double SurfaceYPlane::_distance(Position r, Direction u, bool coincident) const
 {
   return axis_aligned_plane_distance<1>(r, u, coincident, y0_);
 }
@@ -306,12 +490,12 @@ SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&z0_});
 }
 
-double SurfaceZPlane::evaluate(Position r) const
+double SurfaceZPlane::_evaluate(Position r) const
 {
   return r.z - z0_;
 }
 
-double SurfaceZPlane::distance(Position r, Direction u, bool coincident) const
+double SurfaceZPlane::_distance(Position r, Direction u, bool coincident) const
 {
   return axis_aligned_plane_distance<2>(r, u, coincident, z0_);
 }
@@ -346,12 +530,12 @@ SurfacePlane::SurfacePlane(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&A_, &B_, &C_, &D_});
 }
 
-double SurfacePlane::evaluate(Position r) const
+double SurfacePlane::_evaluate(Position r) const
 {
   return A_ * r.x + B_ * r.y + C_ * r.z - D_;
 }
 
-double SurfacePlane::distance(Position r, Direction u, bool coincident) const
+double SurfacePlane::_distance(Position r, Direction u, bool coincident) const
 {
   const double f = A_ * r.x + B_ * r.y + C_ * r.z - D_;
   const double projection = A_ * u.x + B_ * u.y + C_ * u.z;
@@ -465,12 +649,12 @@ SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node)
   read_coeffs(surf_node, id_, {&y0_, &z0_, &radius_});
 }
 
-double SurfaceXCylinder::evaluate(Position r) const
+double SurfaceXCylinder::_evaluate(Position r) const
 {
   return axis_aligned_cylinder_evaluate<1, 2>(r, y0_, z0_, radius_);
 }
 
-double SurfaceXCylinder::distance(
+double SurfaceXCylinder::_distance(
   Position r, Direction u, bool coincident) const
 {
   return axis_aligned_cylinder_distance<0, 1, 2>(
@@ -508,12 +692,12 @@ SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node)
   read_coeffs(surf_node, id_, {&x0_, &z0_, &radius_});
 }
 
-double SurfaceYCylinder::evaluate(Position r) const
+double SurfaceYCylinder::_evaluate(Position r) const
 {
   return axis_aligned_cylinder_evaluate<0, 2>(r, x0_, z0_, radius_);
 }
 
-double SurfaceYCylinder::distance(
+double SurfaceYCylinder::_distance(
   Position r, Direction u, bool coincident) const
 {
   return axis_aligned_cylinder_distance<1, 0, 2>(
@@ -552,12 +736,12 @@ SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
   read_coeffs(surf_node, id_, {&x0_, &y0_, &radius_});
 }
 
-double SurfaceZCylinder::evaluate(Position r) const
+double SurfaceZCylinder::_evaluate(Position r) const
 {
   return axis_aligned_cylinder_evaluate<0, 1>(r, x0_, y0_, radius_);
 }
 
-double SurfaceZCylinder::distance(
+double SurfaceZCylinder::_distance(
   Position r, Direction u, bool coincident) const
 {
   return axis_aligned_cylinder_distance<2, 0, 1>(
@@ -595,7 +779,7 @@ SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_});
 }
 
-double SurfaceSphere::evaluate(Position r) const
+double SurfaceSphere::_evaluate(Position r) const
 {
   const double x = r.x - x0_;
   const double y = r.y - y0_;
@@ -603,7 +787,7 @@ double SurfaceSphere::evaluate(Position r) const
   return x * x + y * y + z * z - radius_ * radius_;
 }
 
-double SurfaceSphere::distance(Position r, Direction u, bool coincident) const
+double SurfaceSphere::_distance(Position r, Direction u, bool coincident) const
 {
   const double x = r.x - x0_;
   const double y = r.y - y0_;
@@ -761,12 +945,12 @@ SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_});
 }
 
-double SurfaceXCone::evaluate(Position r) const
+double SurfaceXCone::_evaluate(Position r) const
 {
   return axis_aligned_cone_evaluate<0, 1, 2>(r, x0_, y0_, z0_, radius_sq_);
 }
 
-double SurfaceXCone::distance(Position r, Direction u, bool coincident) const
+double SurfaceXCone::_distance(Position r, Direction u, bool coincident) const
 {
   return axis_aligned_cone_distance<0, 1, 2>(
     r, u, coincident, x0_, y0_, z0_, radius_sq_);
@@ -793,12 +977,12 @@ SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_});
 }
 
-double SurfaceYCone::evaluate(Position r) const
+double SurfaceYCone::_evaluate(Position r) const
 {
   return axis_aligned_cone_evaluate<1, 0, 2>(r, y0_, x0_, z0_, radius_sq_);
 }
 
-double SurfaceYCone::distance(Position r, Direction u, bool coincident) const
+double SurfaceYCone::_distance(Position r, Direction u, bool coincident) const
 {
   return axis_aligned_cone_distance<1, 0, 2>(
     r, u, coincident, y0_, x0_, z0_, radius_sq_);
@@ -825,12 +1009,12 @@ SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) : Surface(surf_node)
   read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_});
 }
 
-double SurfaceZCone::evaluate(Position r) const
+double SurfaceZCone::_evaluate(Position r) const
 {
   return axis_aligned_cone_evaluate<2, 0, 1>(r, z0_, x0_, y0_, radius_sq_);
 }
 
-double SurfaceZCone::distance(Position r, Direction u, bool coincident) const
+double SurfaceZCone::_distance(Position r, Direction u, bool coincident) const
 {
   return axis_aligned_cone_distance<2, 0, 1>(
     r, u, coincident, z0_, x0_, y0_, radius_sq_);
@@ -858,7 +1042,7 @@ SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) : Surface(surf_node)
     surf_node, id_, {&A_, &B_, &C_, &D_, &E_, &F_, &G_, &H_, &J_, &K_});
 }
 
-double SurfaceQuadric::evaluate(Position r) const
+double SurfaceQuadric::_evaluate(Position r) const
 {
   const double x = r.x;
   const double y = r.y;
@@ -867,7 +1051,7 @@ double SurfaceQuadric::evaluate(Position r) const
          z * (C_ * z + F_ * x + J_) + K_;
 }
 
-double SurfaceQuadric::distance(
+double SurfaceQuadric::_distance(
   Position r, Direction ang, bool coincident) const
 {
   const double& x = r.x;
@@ -1024,7 +1208,7 @@ void SurfaceXTorus::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-double SurfaceXTorus::evaluate(Position r) const
+double SurfaceXTorus::_evaluate(Position r) const
 {
   double x = r.x - x0_;
   double y = r.y - y0_;
@@ -1033,7 +1217,7 @@ double SurfaceXTorus::evaluate(Position r) const
          std::pow(std::sqrt(y * y + z * z) - A_, 2) / (C_ * C_) - 1.;
 }
 
-double SurfaceXTorus::distance(Position r, Direction u, bool coincident) const
+double SurfaceXTorus::_distance(Position r, Direction u, bool coincident) const
 {
   double x = r.x - x0_;
   double y = r.y - y0_;
@@ -1077,7 +1261,7 @@ void SurfaceYTorus::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-double SurfaceYTorus::evaluate(Position r) const
+double SurfaceYTorus::_evaluate(Position r) const
 {
   double x = r.x - x0_;
   double y = r.y - y0_;
@@ -1086,7 +1270,7 @@ double SurfaceYTorus::evaluate(Position r) const
          std::pow(std::sqrt(x * x + z * z) - A_, 2) / (C_ * C_) - 1.;
 }
 
-double SurfaceYTorus::distance(Position r, Direction u, bool coincident) const
+double SurfaceYTorus::_distance(Position r, Direction u, bool coincident) const
 {
   double x = r.x - x0_;
   double y = r.y - y0_;
@@ -1130,7 +1314,7 @@ void SurfaceZTorus::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-double SurfaceZTorus::evaluate(Position r) const
+double SurfaceZTorus::_evaluate(Position r) const
 {
   double x = r.x - x0_;
   double y = r.y - y0_;
@@ -1139,7 +1323,7 @@ double SurfaceZTorus::evaluate(Position r) const
          std::pow(std::sqrt(x * x + y * y) - A_, 2) / (C_ * C_) - 1.;
 }
 
-double SurfaceZTorus::distance(Position r, Direction u, bool coincident) const
+double SurfaceZTorus::_distance(Position r, Direction u, bool coincident) const
 {
   double x = r.x - x0_;
   double y = r.y - y0_;
diff --git a/src/universe.cpp b/src/universe.cpp
index b4ef6b4b265..89cbe09dd93 100644
--- a/src/universe.cpp
+++ b/src/universe.cpp
@@ -39,11 +39,14 @@ void Universe::to_hdf5(hid_t universes_group) const
 
 bool Universe::find_cell(GeometryState& p) const
 {
-  const auto& cells {
-    !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
+  const auto& cells {!partitioner_ ? cells_
+                                   : partitioner_->get_cells(p.r_local(),
+                                       p.u_local(), p.time(), p.speed())};
 
   Position r {p.r_local()};
   Position u {p.u_local()};
+  double t {p.time()};
+  double speed {p.speed()};
   auto surf = p.surface();
   int32_t i_univ = p.lowest_coord().universe;
 
@@ -51,7 +54,7 @@ bool Universe::find_cell(GeometryState& p) const
     if (model::cells[i_cell]->universe_ != i_univ)
       continue;
     // Check if this cell contains the particle
-    if (model::cells[i_cell]->contains(r, u, surf)) {
+    if (model::cells[i_cell]->contains(r, u, t, speed, surf)) {
       p.lowest_coord().cell = i_cell;
       return true;
     }
@@ -178,7 +181,7 @@ UniversePartitioner::UniversePartitioner(const Universe& univ)
 }
 
 const vector<int32_t>& UniversePartitioner::get_cells(
-  Position r, Direction u) const
+  Position r, Direction u, double t, double speed) const
 {
   // Perform a binary search for the partition containing the given coordinates.
   int left = 0;
@@ -187,7 +190,7 @@ const vector<int32_t>& UniversePartitioner::get_cells(
   while (true) {
     // Check the sense of the coordinates for the current surface.
     const auto& surf = *model::surfaces[surfs_[middle]];
-    if (surf.sense(r, u)) {
+    if (surf.sense(r, u, t, speed)) {
       // The coordinates lie in the positive halfspace.  Recurse if there are
       // more surfaces to check.  Otherwise, return the cells on the positive
       // side of this surface.