diff --git a/openmc/model/model.py b/openmc/model/model.py
index 6e4c1c5856d..b4427a4cf59 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -1160,6 +1160,13 @@ def plot(
         y_min = (origin[y] - 0.5*width[1]) * axis_scaling_factor[axis_units]
         y_max = (origin[y] + 0.5*width[1]) * axis_scaling_factor[axis_units]
 
+        # Determine whether any materials contains macroscopic data and if so,
+        # set energy mode accordingly
+        for mat in self.geometry.get_all_materials().values():
+            if mat._macroscopic is not None:
+                self.settings.energy_mode = 'multi-group'
+                break
+
         # Get ID map from the C API
         id_map = self.id_map(
             origin=origin,
@@ -1179,8 +1186,8 @@ def plot(
 
         # Convert ID map to RGB image
         img = id_map_to_rgb(
-            id_map=id_map, 
-            color_by=color_by, 
+            id_map=id_map,
+            color_by=color_by,
             colors=colors,
             overlap_color=overlap_color
         )
@@ -1217,7 +1224,7 @@ def plot(
                 extent=(x_min, x_max, y_min, y_max),
                 **contour_kwargs
             )
-            
+
             # If only showing outline, set the axis limits and aspect explicitly
             if outline == 'only':
                 axes.set_xlim(x_min, x_max)
diff --git a/tests/unit_tests/test_universe.py b/tests/unit_tests/test_universe.py
index efe8552a648..411a00ca77a 100644
--- a/tests/unit_tests/test_universe.py
+++ b/tests/unit_tests/test_universe.py
@@ -104,6 +104,37 @@ def test_plot(run_in_tmpdir, sphere_model):
     plt.close('all')
 
 
+def test_mg_plot(run_in_tmpdir):
+    # Create a simple universe with macroscopic data
+    h2o_data = openmc.Macroscopic('LWTR')
+    water = openmc.Material(name='Water')
+    water.set_density('macro', 1.0)
+    water.add_macroscopic(h2o_data)
+    sph = openmc.Sphere(r=10, boundary_type="vacuum")
+
+    # Create MGXS library and export to HDF5
+    groups = openmc.mgxs.EnergyGroups([1e-5, 20.0e6])
+    h2o_xsdata = openmc.XSdata('LWTR', groups)
+    h2o_xsdata.order = 0
+    h2o_xsdata.set_total([1.0])
+    h2o_xsdata.set_absorption([0.5])
+    scatter_matrix = np.array([[[0.5]]])
+    scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
+    h2o_xsdata.set_scatter_matrix(scatter_matrix)
+    mg_library = openmc.MGXSLibrary(groups)
+    mg_library.add_xsdatas([h2o_xsdata])
+    mg_library.export_to_hdf5('mgxs.h5')
+
+    # Set MG cross sections in config and plot
+    with openmc.config.patch('mg_cross_sections', 'mgxs.h5'):
+        (-sph).plot(width=(200, 200), basis='yz', color_by='cell')
+        (-sph).plot(width=(200, 200), basis='yz', color_by='material')
+
+    # Close plots to avoid warning
+    import matplotlib.pyplot as plt
+    plt.close('all')
+
+
 def test_get_nuclides(uo2):
     c = openmc.Cell(fill=uo2)
     univ = openmc.Universe(cells=[c])