diff --git a/openmc/plots.py b/openmc/plots.py
index 072a9a319e3..0489345ff31 100644
--- a/openmc/plots.py
+++ b/openmc/plots.py
@@ -624,14 +624,15 @@ def to_xml_element(self):
         return element
 
 
-class Plot(PlotBase):
-    """Definition of a finite region of space to be plotted.
+class SlicePlot(PlotBase):
+    """Definition of a 2D slice plot of the geometry.
 
-    OpenMC is capable of generating two-dimensional slice plots, or
-    three-dimensional voxel or projection plots. Colors that are used in plots can be given as
-    RGB tuples, e.g. (255, 255, 255) would be white, or by a string indicating a
+    Colors that are used in plots can be given as RGB tuples, e.g. 
+    (255, 255, 255) would be white, or by a string indicating a
     valid `SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.
 
+    .. versionadded:: 0.15.1
+
     Parameters
     ----------
     plot_id : int
@@ -646,7 +647,7 @@ class Plot(PlotBase):
     name : str
         Name of the plot
     pixels : Iterable of int
-        Number of pixels to use in each direction
+        Number of pixels to use in each direction (2 values)
     filename : str
         Path to write the plot to
     color_by : {'cell', 'material'}
@@ -669,11 +670,9 @@ class Plot(PlotBase):
     level : int
         Universe depth to plot at
     width : Iterable of float
-        Width of the plot in each basis direction
+        Width of the plot in each basis direction (2 values)
     origin : tuple or list of ndarray
-        Origin (center) of the plot
-    type : {'slice', 'voxel'}
-        The type of the plot
+        Origin (center) of the plot (3 values)
     basis : {'xy', 'xz', 'yz'}
         The basis directions for the plot
     meshlines : dict
@@ -686,10 +685,21 @@ def __init__(self, plot_id=None, name=''):
         super().__init__(plot_id, name)
         self._width = [4.0, 4.0]
         self._origin = [0., 0., 0.]
-        self._type = 'slice'
         self._basis = 'xy'
         self._meshlines = None
 
+    @property
+    def pixels(self):
+        return self._pixels
+
+    @pixels.setter
+    def pixels(self, pixels):
+        cv.check_type('plot pixels', pixels, Iterable, Integral)
+        cv.check_length('plot pixels', pixels, 2, 2)
+        for dim in pixels:
+            cv.check_greater_than('plot pixels', dim, 0)
+        self._pixels = pixels
+
     @property
     def width(self):
         return self._width
@@ -697,7 +707,7 @@ def width(self):
     @width.setter
     def width(self, width):
         cv.check_type('plot width', width, Iterable, Real)
-        cv.check_length('plot width', width, 2, 3)
+        cv.check_length('plot width', width, 2, 2)
         self._width = width
 
     @property
@@ -710,15 +720,6 @@ def origin(self, origin):
         cv.check_length('plot origin', origin, 3)
         self._origin = origin
 
-    @property
-    def type(self):
-        return self._type
-
-    @type.setter
-    def type(self, plottype):
-        cv.check_value('plot type', plottype, ['slice', 'voxel'])
-        self._type = plottype
-
     @property
     def basis(self):
         return self._basis
@@ -761,11 +762,10 @@ def meshlines(self, meshlines):
         self._meshlines = meshlines
 
     def __repr__(self):
-        string = 'Plot\n'
+        string = 'SlicePlot\n'
         string += '{: <16}=\t{}\n'.format('\tID', self._id)
         string += '{: <16}=\t{}\n'.format('\tName', self._name)
         string += '{: <16}=\t{}\n'.format('\tFilename', self._filename)
-        string += '{: <16}=\t{}\n'.format('\tType', self._type)
         string += '{: <16}=\t{}\n'.format('\tBasis', self._basis)
         string += '{: <16}=\t{}\n'.format('\tWidth', self._width)
         string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin)
@@ -881,7 +881,7 @@ def highlight_domains(self, geometry, domains, seed=1,
                 self._colors[domain] = (r, g, b)
 
     def to_xml_element(self):
-        """Return XML representation of the slice/voxel plot
+        """Return XML representation of the slice plot
 
         Returns
         -------
@@ -891,10 +891,8 @@ def to_xml_element(self):
         """
 
         element = super().to_xml_element()
-        element.set("type", self._type)
-
-        if self._type == 'slice':
-            element.set("basis", self._basis)
+        element.set("type", "slice")
+        element.set("basis", self._basis)
 
         subelement = ET.SubElement(element, "origin")
         subelement.text = ' '.join(map(str, self._origin))
@@ -940,8 +938,8 @@ def from_xml_element(cls, elem):
 
         Returns
         -------
-        openmc.Plot
-            Plot object
+        openmc.SlicePlot
+            SlicePlot object
 
         """
         plot_id = int(get_text(elem, "id"))
@@ -950,9 +948,7 @@ def from_xml_element(cls, elem):
         if "filename" in elem.keys():
             plot.filename = get_text(elem, "filename")
         plot.color_by = get_text(elem, "color_by")
-        plot.type = get_text(elem, "type")
-        if plot.type == 'slice':
-            plot.basis = get_text(elem, "basis")
+        plot.basis = get_text(elem, "basis")
 
         plot.origin = tuple(get_elem_list(elem, "origin", float))
         plot.width = tuple(get_elem_list(elem, "width", float))
@@ -1034,9 +1030,246 @@ def to_ipython_image(self, openmc_exec='openmc', cwd='.'):
         # Return produced image
         return _get_plot_image(self, cwd)
 
+
+
+class VoxelPlot(PlotBase):
+    """Definition of a 3D voxel plot of the geometry.
+
+    Colors that are used in plots can be given as RGB tuples, e.g. 
+    (255, 255, 255) would be white, or by a string indicating a
+    valid `SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.
+
+    .. versionadded:: 0.15.1
+
+    Parameters
+    ----------
+    plot_id : int
+        Unique identifier for the plot
+    name : str
+        Name of the plot
+
+    Attributes
+    ----------
+    id : int
+        Unique identifier
+    name : str
+        Name of the plot
+    pixels : Iterable of int
+        Number of pixels to use in each direction (3 values)
+    filename : str
+        Path to write the plot to
+    color_by : {'cell', 'material'}
+        Indicate whether the plot should be colored by cell or by material
+    background : Iterable of int or str
+        Color of the background
+    mask_components : Iterable of openmc.Cell or openmc.Material or int
+        The cells or materials (or corresponding IDs) to mask
+    mask_background : Iterable of int or str
+        Color to apply to all cells/materials listed in mask_components
+    show_overlaps : bool
+        Indicate whether or not overlapping regions are shown
+    overlap_color : Iterable of int or str
+        Color to apply to overlapping regions
+    colors : dict
+        Dictionary indicating that certain cells/materials should be
+        displayed with a particular color. The keys can be of type
+        :class:`~openmc.Cell`, :class:`~openmc.Material`, or int (ID for a
+        cell/material).
+    level : int
+        Universe depth to plot at
+    width : Iterable of float
+        Width of the plot in each dimension (3 values)
+    origin : tuple or list of ndarray
+        Origin (center) of the plot (3 values)
+
+    """
+
+    def __init__(self, plot_id=None, name=''):
+        super().__init__(plot_id, name)
+        self._width = [4.0, 4.0, 4.0]
+        self._origin = [0., 0., 0.]
+        self._pixels = [400, 400, 400]
+
+    @property
+    def pixels(self):
+        return self._pixels
+
+    @pixels.setter
+    def pixels(self, pixels):
+        cv.check_type('plot pixels', pixels, Iterable, Integral)
+        cv.check_length('plot pixels', pixels, 3, 3)
+        for dim in pixels:
+            cv.check_greater_than('plot pixels', dim, 0)
+        self._pixels = pixels
+
+    @property
+    def width(self):
+        return self._width
+
+    @width.setter
+    def width(self, width):
+        cv.check_type('plot width', width, Iterable, Real)
+        cv.check_length('plot width', width, 3, 3)
+        self._width = width
+
+    @property
+    def origin(self):
+        return self._origin
+
+    @origin.setter
+    def origin(self, origin):
+        cv.check_type('plot origin', origin, Iterable, Real)
+        cv.check_length('plot origin', origin, 3)
+        self._origin = origin
+
+    def __repr__(self):
+        string = 'VoxelPlot\n'
+        string += '{: <16}=\t{}\n'.format('\tID', self._id)
+        string += '{: <16}=\t{}\n'.format('\tName', self._name)
+        string += '{: <16}=\t{}\n'.format('\tFilename', self._filename)
+        string += '{: <16}=\t{}\n'.format('\tWidth', self._width)
+        string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin)
+        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)
+        string += '{: <16}=\t{}\n'.format('\tMask components',
+                                          self._mask_components)
+        string += '{: <16}=\t{}\n'.format('\tMask background',
+                                          self._mask_background)
+        string += '{: <16}=\t{}\n'.format('\tOverlap Color',
+                                          self._overlap_color)
+        string += '{: <16}=\t{}\n'.format('\tColors', self._colors)
+        string += '{: <16}=\t{}\n'.format('\tLevel', self._level)
+        return string
+
+    def to_xml_element(self):
+        """Return XML representation of the voxel plot
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing plot data
+
+        """
+
+        element = super().to_xml_element()
+        element.set("type", "voxel")
+
+        subelement = ET.SubElement(element, "origin")
+        subelement.text = ' '.join(map(str, self._origin))
+
+        subelement = ET.SubElement(element, "width")
+        subelement.text = ' '.join(map(str, self._width))
+
+        if self._colors:
+            self._colors_to_xml(element)
+
+        if self._show_overlaps:
+            subelement = ET.SubElement(element, "show_overlaps")
+            subelement.text = "true"
+
+            if self._overlap_color is not None:
+                color = self._overlap_color
+                if isinstance(color, str):
+                    color = _SVG_COLORS[color.lower()]
+                subelement = ET.SubElement(element, "overlap_color")
+                subelement.text = ' '.join(str(x) for x in color)
+
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate plot object from an XML element
+
+        Parameters
+        ----------
+        elem : lxml.etree._Element
+            XML element
+
+        Returns
+        -------
+        openmc.VoxelPlot
+            VoxelPlot object
+
+        """
+        plot_id = int(get_text(elem, "id"))
+        name = get_text(elem, 'name', '')
+        plot = cls(plot_id, name)
+        if "filename" in elem.keys():
+            plot.filename = get_text(elem, "filename")
+        plot.color_by = get_text(elem, "color_by")
+
+        plot.origin = tuple(get_elem_list(elem, "origin", float))
+        plot.width = tuple(get_elem_list(elem, "width", float))
+        plot.pixels = tuple(get_elem_list(elem, "pixels"))
+        background = get_elem_list(elem, "background")
+        if background is not None:
+            plot._background = tuple(background)
+
+        # Set plot colors
+        colors = {}
+        for color_elem in elem.findall("color"):
+            uid = int(get_text(color_elem, "id"))
+            colors[uid] = tuple(get_elem_list(color_elem, "rgb", int))
+        plot.colors = colors
+
+        # Set masking information
+        mask_elem = elem.find("mask")
+        if mask_elem is not None:
+            plot.mask_components = get_elem_list(mask_elem, "components", int)
+            background = get_elem_list(mask_elem, "background", int)
+            if background is not None:
+                plot.mask_background = tuple(background)
+
+        # show overlaps
+        overlap = get_text(elem, "show_overlaps")
+        if overlap is not None:
+            plot.show_overlaps = (overlap in ('true', '1'))
+        overlap_color = get_elem_list(elem, "overlap_color", int)
+        if overlap_color is not None:
+            plot.overlap_color = tuple(overlap_color)
+
+        # Set universe level
+        level = get_text(elem, "level")
+        if level is not None:
+            plot.level = int(level)
+
+        return plot
+
+    def to_ipython_image(self, openmc_exec='openmc', cwd='.'):
+        """Render plot as an image
+
+        This method runs OpenMC in plotting mode to produce a .png file.
+
+        .. versionchanged:: 0.13.0
+            The *convert_exec* argument was removed since OpenMC now produces
+            .png images directly.
+
+        Parameters
+        ----------
+        openmc_exec : str
+            Path to OpenMC executable
+        cwd : str, optional
+            Path to working directory to run in
+
+        Returns
+        -------
+        IPython.display.Image
+            Image generated
+
+        """
+        # Create plots.xml
+        Plots([self]).export_to_xml(cwd)
+
+        # Run OpenMC in geometry plotting mode
+        openmc.plot_geometry(False, openmc_exec, cwd)
+
+        # Return produced image
+        return _get_plot_image(self, cwd)
+
     def to_vtk(self, output: PathLike | None = None,
                openmc_exec: str = 'openmc', cwd: str = '.'):
-        """Render plot as an voxel image
+        """Render plot as a voxel image
 
         This method runs OpenMC in plotting mode to produce a .vti file.
 
@@ -1057,10 +1290,6 @@ def to_vtk(self, output: PathLike | None = None,
             Path of the .vti file produced
 
         """
-        if self.type != 'voxel':
-            raise ValueError(
-                'Generating a VTK file only works for voxel plots')
-
         # Create plots.xml
         Plots([self]).export_to_xml(cwd)
 
@@ -1080,6 +1309,286 @@ def to_vtk(self, output: PathLike | None = None,
         return voxel_to_vtk(h5_voxel_file, output)
 
 
+class Plot(SlicePlot):
+    """Legacy Plot class for backward compatibility.
+
+    .. deprecated:: 0.15.1
+        Use :class:`SlicePlot` for 2D slice plots or :class:`VoxelPlot` for 3D voxel plots.
+
+    """
+
+    def __init__(self, plot_id=None, name=''):
+        import warnings
+        warnings.warn(
+            "The Plot class is deprecated. Use SlicePlot for 2D slice plots "
+            "or VoxelPlot for 3D voxel plots.",
+            DeprecationWarning,
+            stacklevel=2
+        )
+        super().__init__(plot_id, name)
+        self._type = 'slice'
+
+    @property
+    def type(self):
+        return self._type
+
+    @type.setter
+    def type(self, plottype):
+        cv.check_value('plot type', plottype, ['slice', 'voxel'])
+        self._type = plottype
+
+    @property
+    def width(self):
+        return self._width
+
+    @width.setter
+    def width(self, width):
+        cv.check_type('plot width', width, Iterable, Real)
+        cv.check_length('plot width', width, 2, 3)
+        self._width = width
+
+    @property
+    def pixels(self):
+        return self._pixels
+
+    @pixels.setter
+    def pixels(self, pixels):
+        cv.check_type('plot pixels', pixels, Iterable, Integral)
+        cv.check_length('plot pixels', pixels, 2, 3)
+        for dim in pixels:
+            cv.check_greater_than('plot pixels', dim, 0)
+        self._pixels = pixels
+
+    def to_xml_element(self):
+        """Return XML representation of the plot
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing plot data
+
+        """
+        if self._type == 'voxel':
+            # Convert to VoxelPlot for proper XML generation
+            voxel_plot = VoxelPlot(self.id, self.name)
+            voxel_plot._width = self._width
+            voxel_plot._origin = self._origin
+            voxel_plot._pixels = self._pixels
+            voxel_plot._filename = self._filename
+            voxel_plot._color_by = self._color_by
+            voxel_plot._background = self._background
+            voxel_plot._mask_components = self._mask_components
+            voxel_plot._mask_background = self._mask_background
+            voxel_plot._show_overlaps = self._show_overlaps
+            voxel_plot._overlap_color = self._overlap_color
+            voxel_plot._colors = self._colors
+            voxel_plot._level = self._level
+            return voxel_plot.to_xml_element()
+        else:
+            # Use parent SlicePlot implementation
+            return super().to_xml_element()
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate plot object from an XML element
+
+        Parameters
+        ----------
+        elem : lxml.etree._Element
+            XML element
+
+        Returns
+        -------
+        openmc.Plot
+            Plot object
+
+        """
+        plot_type = get_text(elem, 'type')
+        if plot_type == 'voxel':
+            # Create a Plot but with voxel type
+            plot_id = int(get_text(elem, "id"))
+            name = get_text(elem, 'name', '')
+            plot = cls(plot_id, name)
+            plot._type = 'voxel'
+            if "filename" in elem.keys():
+                plot.filename = get_text(elem, "filename")
+            plot.color_by = get_text(elem, "color_by")
+            
+            plot.origin = tuple(get_elem_list(elem, "origin", float))
+            plot.width = tuple(get_elem_list(elem, "width", float))
+            plot.pixels = tuple(get_elem_list(elem, "pixels"))
+            background = get_elem_list(elem, "background")
+            if background is not None:
+                plot._background = tuple(background)
+            
+            # Set plot colors
+            colors = {}
+            for color_elem in elem.findall("color"):
+                uid = int(get_text(color_elem, "id"))
+                colors[uid] = tuple(get_elem_list(color_elem, "rgb", int))
+            plot.colors = colors
+            
+            # Set masking information
+            mask_elem = elem.find("mask")
+            if mask_elem is not None:
+                plot.mask_components = get_elem_list(mask_elem, "components", int)
+                background = get_elem_list(mask_elem, "background", int)
+                if background is not None:
+                    plot.mask_background = tuple(background)
+            
+            # show overlaps
+            overlap = get_text(elem, "show_overlaps")
+            if overlap is not None:
+                plot.show_overlaps = (overlap in ('true', '1'))
+            overlap_color = get_elem_list(elem, "overlap_color", int)
+            if overlap_color is not None:
+                plot.overlap_color = tuple(overlap_color)
+            
+            # Set universe level
+            level = get_text(elem, "level")
+            if level is not None:
+                plot.level = int(level)
+            
+            return plot
+        else:
+            # Use SlicePlot.from_xml_element but return as Plot
+            plot_id = int(get_text(elem, "id"))
+            name = get_text(elem, 'name', '')
+            plot = cls(plot_id, name)
+            if "filename" in elem.keys():
+                plot.filename = get_text(elem, "filename")
+            plot.color_by = get_text(elem, "color_by")
+            plot.basis = get_text(elem, "basis")
+            
+            plot.origin = tuple(get_elem_list(elem, "origin", float))
+            plot.width = tuple(get_elem_list(elem, "width", float))
+            plot.pixels = tuple(get_elem_list(elem, "pixels"))
+            background = get_elem_list(elem, "background")
+            if background is not None:
+                plot._background = tuple(background)
+            
+            # Set plot colors
+            colors = {}
+            for color_elem in elem.findall("color"):
+                uid = int(get_text(color_elem, "id"))
+                colors[uid] = tuple(get_elem_list(color_elem, "rgb", int))
+            plot.colors = colors
+            
+            # Set masking information
+            mask_elem = elem.find("mask")
+            if mask_elem is not None:
+                plot.mask_components = get_elem_list(mask_elem, "components", int)
+                background = get_elem_list(mask_elem, "background", int)
+                if background is not None:
+                    plot.mask_background = tuple(background)
+            
+            # show overlaps
+            overlap = get_text(elem, "show_overlaps")
+            if overlap is not None:
+                plot.show_overlaps = (overlap in ('true', '1'))
+            overlap_color = get_elem_list(elem, "overlap_color", int)
+            if overlap_color is not None:
+                plot.overlap_color = tuple(overlap_color)
+            
+            # Set universe level
+            level = get_text(elem, "level")
+            if level is not None:
+                plot.level = int(level)
+            
+            # Set meshlines
+            mesh_elem = elem.find("meshlines")
+            if mesh_elem is not None:
+                meshlines = {'type': get_text(mesh_elem, 'meshtype')}
+                mesh_id = get_text(mesh_elem, 'id')
+                if mesh_id is not None:
+                    meshlines['id'] = int(mesh_id)
+                linewidth = get_text(mesh_elem, 'linewidth')
+                if linewidth is not None:
+                    meshlines['linewidth'] = int(linewidth)
+                color = get_elem_list(mesh_elem, 'color', int)
+                if color is not None:
+                    meshlines['color'] = tuple(color)
+                plot.meshlines = meshlines
+            
+            return plot
+
+    @classmethod 
+    def from_geometry(cls, geometry,
+                      basis: str = 'xy',
+                      slice_coord: float = 0.):
+        """Generate plot from a geometry object
+        
+        Parameters
+        ----------
+        geometry : openmc.Geometry
+            Geometry object to create plot from
+        basis : {'xy', 'xz', 'yz'}
+            The basis directions
+        slice_coord : float
+            The position of the slice
+            
+        Returns
+        -------
+        openmc.Plot
+            Plot object
+            
+        """
+        import warnings
+        # Suppress deprecation warning when called as class method
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", DeprecationWarning)
+            plot = super().from_geometry(geometry, basis, slice_coord)
+            # Convert from SlicePlot to Plot
+            new_plot = cls.__new__(cls)
+            new_plot.__dict__.update(plot.__dict__)
+            new_plot._type = 'slice'
+            return new_plot
+
+    def to_vtk(self, output: PathLike | None = None,
+               openmc_exec: str = 'openmc', cwd: str = '.'):
+        """Render plot as a voxel image
+
+        This method runs OpenMC in plotting mode to produce a .vti file.
+
+        .. versionadded:: 0.14.0
+
+        Parameters
+        ----------
+        output : path-like
+            Path of the output .vti file produced
+        openmc_exec : str
+            Path to OpenMC executable
+        cwd : str, optional
+            Path to working directory to run in
+
+        Returns
+        -------
+        Path
+            Path of the .vti file produced
+
+        """
+        if self.type != 'voxel':
+            raise ValueError(
+                'Generating a VTK file only works for voxel plots')
+        
+        # Convert to VoxelPlot and call its to_vtk method
+        voxel_plot = VoxelPlot(self.id, self.name)
+        voxel_plot._width = self._width
+        voxel_plot._origin = self._origin
+        voxel_plot._pixels = self._pixels
+        voxel_plot._filename = self._filename
+        voxel_plot._color_by = self._color_by
+        voxel_plot._background = self._background
+        voxel_plot._mask_components = self._mask_components
+        voxel_plot._mask_background = self._mask_background
+        voxel_plot._show_overlaps = self._show_overlaps
+        voxel_plot._overlap_color = self._overlap_color
+        voxel_plot._colors = self._colors
+        voxel_plot._level = self._level
+        
+        return voxel_plot.to_vtk(output, openmc_exec, cwd)
+
+
 class RayTracePlot(PlotBase):
     """Definition of a camera's view of OpenMC geometry
 
@@ -1901,8 +2410,13 @@ def from_xml_element(cls, elem):
                 plots.append(WireframeRayTracePlot.from_xml_element(e))
             elif plot_type == 'solid_raytrace':
                 plots.append(SolidRayTracePlot.from_xml_element(e))
-            elif plot_type in ('slice', 'voxel'):
-                plots.append(Plot.from_xml_element(e))
+            elif plot_type == 'slice':
+                plots.append(SlicePlot.from_xml_element(e))
+            elif plot_type == 'voxel':
+                plots.append(VoxelPlot.from_xml_element(e))
+            elif plot_type is None:
+                # For backward compatibility, assume slice if no type specified
+                plots.append(SlicePlot.from_xml_element(e))
             else:
                 raise ValueError("Unknown plot type: {}".format(plot_type))
         return plots
diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp
index b8e78071003..2439d27837a 100644
--- a/src/mcpl_interface.cpp
+++ b/src/mcpl_interface.cpp
@@ -67,6 +67,8 @@ using mcpl_hdr_set_srcname_fpt = void (*)(
 using mcpl_add_particle_fpt = void (*)(
   mcpl_outfile_t* outfile_handle, const mcpl_particle_repr_t* particle);
 using mcpl_close_outfile_fpt = void (*)(mcpl_outfile_t* outfile_handle);
+using mcpl_hdr_add_data_fpt = void (*)(mcpl_outfile_t* outfile_handle,
+  const char* key, uint32_t datalength, const char* data);
 using mcpl_hdr_add_stat_sum_fpt = void (*)(
   mcpl_outfile_t* outfile_handle, const char* key, double value);
 
@@ -112,6 +114,7 @@ struct McplApi {
   mcpl_hdr_set_srcname_fpt hdr_set_srcname;
   mcpl_add_particle_fpt add_particle;
   mcpl_close_outfile_fpt close_outfile;
+  mcpl_hdr_add_data_fpt hdr_add_data;
   mcpl_hdr_add_stat_sum_fpt hdr_add_stat_sum;
 
   explicit McplApi(LibraryHandleType lib_handle)
@@ -151,6 +154,15 @@ struct McplApi {
     close_outfile = reinterpret_cast<mcpl_close_outfile_fpt>(
       load_symbol_platform("mcpl_close_outfile"));
 
+    // Try to load mcpl_hdr_add_data (available in MCPL >= 2.1.0)
+    // Set to nullptr if not available for graceful fallback
+    try {
+      hdr_add_data = reinterpret_cast<mcpl_hdr_add_data_fpt>(
+        load_symbol_platform("mcpl_hdr_add_data"));
+    } catch (const std::runtime_error&) {
+      hdr_add_data = nullptr;
+    }
+
     // Try to load mcpl_hdr_add_stat_sum (available in MCPL >= 2.1.0)
     // Set to nullptr if not available for graceful fallback
     try {
diff --git a/tests/unit_tests/test_plots.py b/tests/unit_tests/test_plots.py
index fad574ee697..3f01985a0ca 100644
--- a/tests/unit_tests/test_plots.py
+++ b/tests/unit_tests/test_plots.py
@@ -80,8 +80,7 @@ def test_voxel_plot(run_in_tmpdir):
     geometry.export_to_xml()
     materials = openmc.Materials()
     materials.export_to_xml()
-    vox_plot = openmc.Plot()
-    vox_plot.type = 'voxel'
+    vox_plot = openmc.VoxelPlot()
     vox_plot.id = 12
     vox_plot.width = (1500., 1500., 1500.)
     vox_plot.pixels = (200, 200, 200)
@@ -97,8 +96,9 @@ def test_voxel_plot(run_in_tmpdir):
     assert Path('h5_voxel_plot.h5').is_file()
     assert Path('another_test_voxel_plot.vti').is_file()
 
-    slice_plot = openmc.Plot()
-    with pytest.raises(ValueError):
+    # SlicePlot should not have to_vtk method
+    slice_plot = openmc.SlicePlot()
+    with pytest.raises(AttributeError):
         slice_plot.to_vtk('shimmy.vti')
 
 
@@ -160,7 +160,7 @@ def test_from_geometry():
 
 
 def test_highlight_domains():
-    plot = openmc.Plot()
+    plot = openmc.SlicePlot()
     plot.color_by = 'material'
     plots = openmc.Plots([plot])
 
@@ -200,11 +200,11 @@ def test_to_xml_element_proj(myprojectionplot):
 
 
 def test_plots(run_in_tmpdir):
-    p1 = openmc.Plot(name='plot1')
+    p1 = openmc.SlicePlot(name='plot1')
     p1.origin = (5., 5., 5.)
     p1.colors = {10: (255, 100, 0)}
     p1.mask_components = [2, 4, 6]
-    p2 = openmc.Plot(name='plot2')
+    p2 = openmc.SlicePlot(name='plot2')
     p2.origin = (-3., -3., -3.)
     plots = openmc.Plots([p1, p2])
     assert len(plots) == 2
@@ -213,7 +213,7 @@ def test_plots(run_in_tmpdir):
     plots = openmc.Plots([p1, p2, p3])
     assert len(plots) == 3
 
-    p4 = openmc.Plot(name='plot4')
+    p4 = openmc.VoxelPlot(name='plot4')
     plots.append(p4)
     assert len(plots) == 4
 
@@ -288,10 +288,9 @@ def test_phong_plot_roundtrip():
 def test_plot_directory(run_in_tmpdir):
     pwr_pin = openmc.examples.pwr_pin_cell()
 
-    # create a standard plot, expected to work
-    plot = openmc.Plot()
+    # create a standard slice plot, expected to work
+    plot = openmc.SlicePlot()
     plot.filename = 'plot_1'
-    plot.type = 'slice'
     plot.pixels = (10, 10)
     plot.color_by = 'material'
     plot.width = (100., 100.)
diff --git a/tests/unit_tests/test_slice_voxel_plots.py b/tests/unit_tests/test_slice_voxel_plots.py
new file mode 100644
index 00000000000..a97ae002ef4
--- /dev/null
+++ b/tests/unit_tests/test_slice_voxel_plots.py
@@ -0,0 +1,327 @@
+"""Tests for SlicePlot and VoxelPlot classes
+
+This module tests the functionality of the new SlicePlot and VoxelPlot
+classes that replace the legacy Plot class.
+"""
+import warnings
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+import openmc
+from openmc.plots import _SVG_COLORS
+
+
+def test_slice_plot_initialization():
+    """Test SlicePlot initialization with defaults"""
+    plot = openmc.SlicePlot()
+    assert plot.width == [4.0, 4.0]
+    assert plot.pixels == [400, 400]
+    assert plot.basis == 'xy'
+    assert plot.origin == [0., 0., 0.]
+
+
+def test_slice_plot_width_validation():
+    """Test that SlicePlot only accepts 2 values for width"""
+    plot = openmc.SlicePlot()
+    
+    # Should accept 2 values
+    plot.width = [10.0, 20.0]
+    assert plot.width == [10.0, 20.0]
+    
+    # Should reject 1 value
+    with pytest.raises(ValueError, match='must be of length "2"'):
+        plot.width = [10.0]
+    
+    # Should reject 3 values
+    with pytest.raises(ValueError, match='must be of length "2"'):
+        plot.width = [10.0, 20.0, 30.0]
+
+
+def test_slice_plot_pixels_validation():
+    """Test that SlicePlot only accepts 2 values for pixels"""
+    plot = openmc.SlicePlot()
+    
+    # Should accept 2 values
+    plot.pixels = [100, 200]
+    assert plot.pixels == [100, 200]
+    
+    # Should reject 1 value
+    with pytest.raises(ValueError, match='must be of length "2"'):
+        plot.pixels = [100]
+    
+    # Should reject 3 values
+    with pytest.raises(ValueError, match='must be of length "2"'):
+        plot.pixels = [100, 200, 300]
+
+
+def test_slice_plot_basis():
+    """Test that SlicePlot has basis attribute"""
+    plot = openmc.SlicePlot()
+    
+    # Test all valid basis values
+    for basis in ['xy', 'xz', 'yz']:
+        plot.basis = basis
+        assert plot.basis == basis
+    
+    # Test invalid basis
+    with pytest.raises(ValueError):
+        plot.basis = 'invalid'
+
+
+def test_slice_plot_meshlines():
+    """Test that SlicePlot has meshlines attribute"""
+    plot = openmc.SlicePlot()
+    
+    meshlines = {
+        'type': 'tally',
+        'id': 1,
+        'linewidth': 2,
+        'color': (255, 0, 0)
+    }
+    plot.meshlines = meshlines
+    assert plot.meshlines == meshlines
+
+
+def test_slice_plot_xml_roundtrip():
+    """Test SlicePlot XML serialization and deserialization"""
+    plot = openmc.SlicePlot(name='test_slice')
+    plot.width = [15.0, 25.0]
+    plot.pixels = [150, 250]
+    plot.basis = 'xz'
+    plot.origin = [1.0, 2.0, 3.0]
+    plot.color_by = 'material'
+    plot.filename = 'test_plot'
+    
+    # Convert to XML and back
+    elem = plot.to_xml_element()
+    new_plot = openmc.SlicePlot.from_xml_element(elem)
+    
+    # Check all attributes preserved
+    assert new_plot.name == plot.name
+    assert new_plot.width == pytest.approx(plot.width)
+    assert new_plot.pixels == tuple(plot.pixels)
+    assert new_plot.basis == plot.basis
+    assert new_plot.origin == pytest.approx(plot.origin)
+    assert new_plot.color_by == plot.color_by
+    assert new_plot.filename == plot.filename
+
+
+def test_slice_plot_from_geometry():
+    """Test creating SlicePlot from geometry"""
+    # Create simple geometry
+    s = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    c = openmc.Cell(region=-s)
+    univ = openmc.Universe(cells=[c])
+    geom = openmc.Geometry(univ)
+    
+    # Test all basis options
+    for basis in ['xy', 'xz', 'yz']:
+        plot = openmc.SlicePlot.from_geometry(geom, basis=basis)
+        assert plot.basis == basis
+        assert plot.width == pytest.approx([20.0, 20.0])
+        assert plot.origin == pytest.approx([0.0, 0.0, 0.0])
+
+
+def test_voxel_plot_initialization():
+    """Test VoxelPlot initialization with defaults"""
+    plot = openmc.VoxelPlot()
+    assert plot.width == [4.0, 4.0, 4.0]
+    assert plot.pixels == [400, 400, 400]
+    assert plot.origin == [0., 0., 0.]
+
+
+def test_voxel_plot_width_validation():
+    """Test that VoxelPlot only accepts 3 values for width"""
+    plot = openmc.VoxelPlot()
+    
+    # Should accept 3 values
+    plot.width = [10.0, 20.0, 30.0]
+    assert plot.width == [10.0, 20.0, 30.0]
+    
+    # Should reject 2 values
+    with pytest.raises(ValueError, match='must be of length "3"'):
+        plot.width = [10.0, 20.0]
+    
+    # Should reject 1 value
+    with pytest.raises(ValueError, match='must be of length "3"'):
+        plot.width = [10.0]
+
+
+def test_voxel_plot_pixels_validation():
+    """Test that VoxelPlot only accepts 3 values for pixels"""
+    plot = openmc.VoxelPlot()
+    
+    # Should accept 3 values
+    plot.pixels = [100, 200, 300]
+    assert plot.pixels == [100, 200, 300]
+    
+    # Should reject 2 values
+    with pytest.raises(ValueError, match='must be of length "3"'):
+        plot.pixels = [100, 200]
+    
+    # Should reject 1 value
+    with pytest.raises(ValueError, match='must be of length "3"'):
+        plot.pixels = [100]
+
+
+def test_voxel_plot_no_basis():
+    """Test that VoxelPlot does not have basis attribute"""
+    plot = openmc.VoxelPlot()
+    assert not hasattr(plot, 'basis')
+    assert not hasattr(plot, '_basis')
+
+
+def test_voxel_plot_no_meshlines():
+    """Test that VoxelPlot does not have meshlines attribute"""
+    plot = openmc.VoxelPlot()
+    assert not hasattr(plot, 'meshlines')
+    assert not hasattr(plot, '_meshlines')
+
+
+def test_voxel_plot_has_to_vtk():
+    """Test that VoxelPlot has to_vtk method"""
+    plot = openmc.VoxelPlot()
+    assert hasattr(plot, 'to_vtk')
+    assert callable(plot.to_vtk)
+
+
+def test_voxel_plot_xml_roundtrip():
+    """Test VoxelPlot XML serialization and deserialization"""
+    plot = openmc.VoxelPlot(name='test_voxel')
+    plot.width = [10.0, 20.0, 30.0]
+    plot.pixels = [100, 200, 300]
+    plot.origin = [1.0, 2.0, 3.0]
+    plot.color_by = 'cell'
+    plot.filename = 'voxel_plot'
+    
+    # Convert to XML and back
+    elem = plot.to_xml_element()
+    new_plot = openmc.VoxelPlot.from_xml_element(elem)
+    
+    # Check all attributes preserved
+    assert new_plot.name == plot.name
+    assert new_plot.width == pytest.approx(plot.width)
+    assert new_plot.pixels == tuple(plot.pixels)
+    assert new_plot.origin == pytest.approx(plot.origin)
+    assert new_plot.color_by == plot.color_by
+    assert new_plot.filename == plot.filename
+
+
+def test_plot_deprecation_warning():
+    """Test that Plot class raises deprecation warning"""
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter("always")
+        plot = openmc.Plot()
+        
+        assert len(w) == 1
+        assert issubclass(w[0].category, DeprecationWarning)
+        assert "deprecated" in str(w[0].message).lower()
+
+
+def test_plot_slice_compatibility():
+    """Test Plot class with slice type"""
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        plot = openmc.Plot()
+        
+        plot.type = 'slice'
+        plot.width = [10.0, 20.0]
+        plot.pixels = [100, 200]
+        plot.basis = 'yz'
+        
+        assert plot.type == 'slice'
+        assert plot.width == [10.0, 20.0]
+        assert plot.pixels == [100, 200]
+        assert plot.basis == 'yz'
+
+
+def test_plot_voxel_compatibility():
+    """Test Plot class with voxel type"""
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        plot = openmc.Plot()
+        
+        plot.type = 'voxel'
+        plot.width = [10.0, 20.0, 30.0]
+        plot.pixels = [100, 200, 300]
+        
+        assert plot.type == 'voxel'
+        assert plot.width == [10.0, 20.0, 30.0]
+        assert plot.pixels == [100, 200, 300]
+
+
+def test_plot_xml_roundtrip_slice():
+    """Test XML roundtrip for Plot with slice type"""
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        plot = openmc.Plot(name='legacy_slice')
+        plot.type = 'slice'
+        plot.width = [5.0, 10.0]
+        plot.pixels = [50, 100]
+        plot.basis = 'xz'
+        
+        elem = plot.to_xml_element()
+        new_plot = openmc.Plot.from_xml_element(elem)
+        
+        assert new_plot.type == plot.type
+        assert new_plot.width == pytest.approx(plot.width)
+        assert new_plot.pixels == tuple(plot.pixels)
+        assert new_plot.basis == plot.basis
+
+
+def test_plot_xml_roundtrip_voxel():
+    """Test XML roundtrip for Plot with voxel type"""
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        plot = openmc.Plot(name='legacy_voxel')
+        plot.type = 'voxel'
+        plot.width = [5.0, 10.0, 15.0]
+        plot.pixels = [50, 100, 150]
+        
+        elem = plot.to_xml_element()
+        new_plot = openmc.Plot.from_xml_element(elem)
+        
+        assert new_plot.type == plot.type
+        assert new_plot.width == pytest.approx(plot.width)
+        assert new_plot.pixels == tuple(plot.pixels)
+
+
+def test_plots_collection_mixed_types():
+    """Test Plots collection with different plot types"""
+    slice_plot = openmc.SlicePlot(name='slice')
+    voxel_plot = openmc.VoxelPlot(name='voxel')
+    wireframe_plot = openmc.WireframeRayTracePlot(name='wireframe')
+    
+    plots = openmc.Plots([slice_plot, voxel_plot, wireframe_plot])
+    
+    assert len(plots) == 3
+    assert isinstance(plots[0], openmc.SlicePlot)
+    assert isinstance(plots[1], openmc.VoxelPlot)
+    assert isinstance(plots[2], openmc.WireframeRayTracePlot)
+
+
+def test_plots_collection_xml_roundtrip(run_in_tmpdir):
+    """Test XML export and import with new plot types"""
+    s1 = openmc.SlicePlot(name='slice1')
+    s1.width = [10.0, 20.0]
+    s1.basis = 'xz'
+    
+    v1 = openmc.VoxelPlot(name='voxel1')
+    v1.width = [10.0, 20.0, 30.0]
+    
+    plots = openmc.Plots([s1, v1])
+    plots.export_to_xml()
+    
+    # Read back
+    new_plots = openmc.Plots.from_xml()
+    
+    assert len(new_plots) == 2
+    assert isinstance(new_plots[0], openmc.SlicePlot)
+    assert isinstance(new_plots[1], openmc.VoxelPlot)
+    assert new_plots[0].name == 'slice1'
+    assert new_plots[1].name == 'voxel1'
+    assert new_plots[0].basis == 'xz'
+    assert new_plots[0].width == pytest.approx([10.0, 20.0])
+    assert new_plots[1].width == pytest.approx([10.0, 20.0, 30.0])
\ No newline at end of file