diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index dc9370ab5..8731a9a13 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -57,6 +57,16 @@ The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also
 have some knowledge of how nuclides transmute and decay. This is handled by the
 :class:`Chain` class.
 
+The :class:`IndependentOperator` class requires a set of fluxes and microscopic
+cross sections. The following function can be used to generate this information:
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myfunction.rst
+
+   get_microxs_and_flux
+
 Minimal Example
 ---------------
 
diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst
index 439e30305..261900ce6 100644
--- a/docs/source/usersguide/depletion.rst
+++ b/docs/source/usersguide/depletion.rst
@@ -197,43 +197,54 @@ across all material instances.
 Transport-independent depletion
 ===============================
 
-.. warning::
+This category of operator uses multigroup microscopic cross sections along with
+multigroup flux spectra to obtain transmutation reaction rates. The cross
+sections are pre-calculated, so there is no need for direct coupling between a
+transport-independent operator and a transport solver. The :mod:`openmc.deplete`
+module offers a single transport-independent operator,
+:class:`~openmc.deplete.IndependentOperator`, and only one operator is needed
+since, in theory, any transport code could calculate the multigroup microscopic
+cross sections. The :class:`~openmc.deplete.IndependentOperator` class has two
+constructors. The default constructor requires a :class:`openmc.Materials`
+instance, a list of multigroup flux arrays, and a list of
+:class:`~openmc.deplete.MicroXS` instances containing multigroup microscopic
+cross sections in units of barns. This might look like the following::
 
-   This feature is still under heavy development and has yet to be rigorously
-   verified. API changes and feature additions are possible and likely in
-   the near future.
-
-This category of operator uses one-group microscopic cross sections to obtain
-transmutation reaction rates. The cross sections are pre-calculated, so there is
-no need for direct coupling between a transport-independent operator and a
-transport solver. The :mod:`openmc.deplete` module offers a single
-transport-independent operator, :class:`~openmc.deplete.IndependentOperator`,
-and only one operator is needed since, in theory, any transport code could
-calcuate the one-group microscopic cross sections.
-
-The :class:`~openmc.deplete.IndependentOperator` class has two constructors.
-The default constructor requires a :class:`openmc.Materials` instance, a
-:class:`~openmc.deplete.MicroXS` instance containing one-group microscoic cross
-sections in units of barns, and a path to a depletion chain file::
-
-    materials = openmc.Materials()
+    materials = openmc.Materials([m1, m2, m3])
     ...
 
-    # load in the microscopic cross sections
-    micro_xs = openmc.deplete.MicroXS.from_csv(micro_xs_path)
+    # Assign fluxes (generated from any code)
+    flux_m1 = numpy.array([...])
+    flux_m2 = numpy.array([...])
+    flux_m3 = numpy.array([...])
+    fluxes = [flux_m1, flux_m2, flux_m3]
 
-    op = openmc.deplete.IndependentOperator(materials, micro_xs, chain_file)
+    # Assign microscopic cross sections
+    micro_m1 = openmc.deplete.MicroXS.from_csv('xs_m1.csv')
+    micro_m2 = openmc.deplete.MicroXS.from_csv('xs_m2.csv')
+    micro_m3 = openmc.deplete.MicroXS.from_csv('xs_m3.csv')
+    micros = [micro_m1, micro_m2, micro_m3]
+
+    # Create operator
+    op = openmc.deplete.IndependentOperator(materials, fluxes, micros)
+
+For more details on the :class:`~openmc.deplete.MicroXS` class, including how to
+use OpenMC's transport solver to generate microscopic cross sections and fluxes
+for use with :class:`~openmc.deplete.IndependentOperator`, see :ref:`micros`.
 
 .. note::
 
-   The same statements from :ref:`coupled-depletion` about which
-   materials are depleted and the requirement for depletable materials to have
-   a specified volume also apply here.
+   The same statements from :ref:`coupled-depletion` about which materials are
+   depleted and the requirement for depletable materials to have a specified
+   volume also apply here.
 
 An alternate constructor,
 :meth:`~openmc.deplete.IndependentOperator.from_nuclides`, accepts a volume and
 dictionary of nuclide concentrations in place of the :class:`openmc.Materials`
-instance::
+instance. Note that while the normal constructor allows multiple materials to be
+depleted with a single operator, the
+:meth:`~openmc.deplete.IndependentOperator.from_nuclides` classmethod only works
+for a single material::
 
     nuclides = {'U234': 8.92e18,
                 'U235': 9.98e20,
@@ -244,6 +255,7 @@ instance::
     volume = 0.5
     op = openmc.deplete.IndependentOperator.from_nuclides(volume,
                                                           nuclides,
+                                                          flux,
                                                           micro_xs,
                                                           chain_file,
                                                           nuc_units='atom/cm3')
@@ -253,18 +265,21 @@ transport-depletion calculation and follow the same steps from there.
 
 .. note::
 
-   Ideally, one-group cross section data should be available for every
-   reaction in the depletion chain. If cross section data is not present for
-   a nuclide in the depletion chain with at least one reaction, that reaction
-   will not be simulated.
+   Ideally, multigroup cross section data should be available for every reaction
+   in the depletion chain. If cross section data is not present for a nuclide in
+   the depletion chain with at least one reaction, that reaction will not be
+   simulated.
+
+.. _micros:
 
 Loading and Generating Microscopic Cross Sections
 -------------------------------------------------
 
-As mentioned earlier, any transport code could be used to calculate one-group
-microscopic cross sections. The :mod:`openmc.deplete` module provides the
-:class:`~openmc.deplete.MicroXS` class, which contains methods to read in
-pre-calculated cross sections from a ``.csv`` file or from data arrays::
+As mentioned above, any transport code could be used to calculate multigroup
+microscopic cross sections and fluxes. The :mod:`openmc.deplete` module provides
+the :class:`~openmc.deplete.MicroXS` class, which can either be instantiated
+from pre-calculated cross sections in a ``.csv`` file or from data arrays
+directly::
 
     micro_xs = MicroXS.from_csv(micro_xs_path)
 
@@ -273,37 +288,31 @@ pre-calculated cross sections from a ``.csv`` file or from data arrays::
     data = np.array([[0.1, 0.2],
                      [0.3, 0.4],
                      [0.01, 0.5]])
-    micro_xs = MicroXS.from_array(nuclides, reactions, data)
+    micro_xs = MicroXS(data, nuclides, reactions)
 
 .. important::
 
-   Both :meth:`~openmc.deplete.MicroXS.from_csv()` and
-   :meth:`~openmc.deplete.MicroXS.from_array()` assume the cross section values
-   provided are in barns by defualt, but have no way of verifying this. Make
-   sure your cross sections are in the correct units before passing to a
+   The cross section values are assumed to be in units of barns. Make sure your
+   cross sections are in the correct units before passing to a
    :class:`~openmc.deplete.IndependentOperator` object.
 
-The :class:`~openmc.deplete.MicroXS` class also contains a method to generate one-group microscopic cross sections using OpenMC's transport solver. The
-:meth:`~openmc.deplete.MicroXS.from_model()` method will produce a
-:class:`~openmc.deplete.MicroXS` instance with microscopic cross section data in
-units of barns::
+Additionally, a convenience function,
+:func:`~openmc.deplete.get_microxs_and_flux`, can provide the needed fluxes and
+cross sections using OpenMC's transport solver::
 
-    import openmc
+    model = openmc.Model()
+    ...
 
-    model = openmc.Model.from_xml()
+    fluxes, micros = openmc.deplete.get_microxs_and_flux(model, materials)
 
-    micro_xs = openmc.deplete.MicroXS.from_model(model,
-                                                 model.materials[0],
-                                                 chain_file)
-
-If you are running :meth:`~openmc.deplete.MicroXS.from_model()` on a cluster
+If you are running :func:`~openmc.deplete.get_microxs_and_flux` on a cluster
 where temporary files are created on a local filesystem that is not shared
 across nodes, you'll need to set an environment variable pointing to a local
 directoy so that each MPI process knows where to store output files used to
 calculate the microscopic cross sections. In order of priority, they are
-:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in
-further details can read the documentation for the `tempfile <https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir>`_ module.
-
+:envvar:`TMPDIR`. :envvar:`TEMP`, and :envvar:`TMP`. Users interested in further
+details can read the documentation for the `tempfile
+<https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir>`_ module.
 
 Caveats
 -------
@@ -325,24 +334,30 @@ normalizing reaction rates:
    the time integrator is a flux, and obtains the reaction rates by multiplying
    the cross sections by the ``source-rate``.
 2. ``fission-q`` normalization, which uses the ``power`` or ``power_density``
-   provided by the time integrator to obtain reaction rates by computing a value
-   for the flux based on this power. The equation we use for this calculation is
+   provided by the time integrator to obtain normalized reaction rates by
+   computing a normalization factor as the ratio of the user-specified power to
+   the "observed" power based on fission reaction rates. The equation for the
+   normalization factor is
 
    .. math::
       :label: fission-q
 
-      \phi = \frac{P}{\sum\limits_i (Q_i \sigma^f_i N_i)}
+      f = \frac{P}{\sum\limits_m \sum\limits_i \left(Q_i N_{i,m} \sum\limits_g
+      \sigma^f_{i,g,m} \phi_{g,m} \right)}
 
    where :math:`P` is the power, :math:`Q_i` is the fission Q value for nuclide
-   :math:`i`, :math:`\sigma_i^f` is the microscopic fission cross section for
-   nuclide :math:`i`, and :math:`N_i` is the number of atoms of nuclide
-   :math:`i`. This equation makes the same assumptions and issues as discussed
-   in :ref:`energy-deposition`. Unfortunately, the proposed solution in that
-   section does not apply here since we are decoupled from transport code.
-   However, there is a method to converge to a more accurate value for flux by
-   using substeps during time integration. `This paper
-   <https://doi.org/10.1016/j.anucene.2016.05.031>`_ provides a good discussion
-   of this method.
+   :math:`i`, :math:`\sigma_{i,g,m}^f` is the microscopic fission cross section
+   for nuclide :math:`i` in energy group :math:`g` for material :math:`m`,
+   :math:`\phi_{g,m}` is the neutron flux in group :math:`g` for material
+   :math:`m`, and :math:`N_{i,m}` is the number of atoms of nuclide :math:`i`
+   for material :math:`m`. Reaction rates are then multiplied by :math:`f` so
+   that the total fission power matches :math:`P`. This equation makes the same
+   assumptions and issues as discussed in :ref:`energy-deposition`.
+   Unfortunately, the proposed solution in that section does not apply here
+   since we are decoupled from transport code. However, there is a method to
+   converge to a more accurate value for flux by using substeps during time
+   integration. `This paper <https://doi.org/10.1016/j.anucene.2016.05.031>`_
+   provides a good discussion of this method.
 
 .. warning::
 
@@ -359,14 +374,14 @@ useful for running many different cases of a particular scenario. A
 transport-independent depletion simulation using ``fission-q`` normalization
 will sum the fission energy values across all materials into :math:`Q_i` in
 Equation :math:numref:`fission-q`, and Equation :math:numref:`fission-q`
-provides the flux we use to calculate the reaction rates in each material.
+provides the normalization factor applied to reaction rates in each material.
 This can be useful for running a scenario with multiple depletable materials
 that are part of the same reactor. This behavior may change in the future.
 
 Time integration
 ~~~~~~~~~~~~~~~~
 
-The values of the one-group microscopic cross sections passed to
+The values of the microscopic cross sections passed to
 :class:`openmc.deplete.IndependentOperator` are fixed for the entire depletion
 simulation. This implicit assumption may produce inaccurate results for certain
 scenarios.
diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py
index 9a1fc6a3f..1c9d96f67 100644
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@@ -196,13 +196,13 @@ class DirectReactionRateHelper(ReactionRateHelper):
         """
         self._rate_tally_means_cache = None
 
-    def get_material_rates(self, mat_id, nuc_index, rx_index):
+    def get_material_rates(self, mat_index, nuc_index, rx_index):
         """Return an array of reaction rates for a material
 
         Parameters
         ----------
-        mat_id : int
-            Unique ID for the requested material
+        mat_index : int
+            Index for the material
         nuc_index : iterable of int
             Index for each nuclide in :attr:`nuclides` in the
             desired reaction rate matrix
@@ -216,7 +216,7 @@ class DirectReactionRateHelper(ReactionRateHelper):
             reaction rates in this material
         """
         self._results_cache.fill(0.0)
-        full_tally_res = self.rate_tally_means[mat_id]
+        full_tally_res = self.rate_tally_means[mat_index]
         for i_tally, (i_nuc, i_rx) in enumerate(product(nuc_index, rx_index)):
             self._results_cache[i_nuc, i_rx] = full_tally_res[i_tally]
 
diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py
index 757278406..aa71942df 100644
--- a/openmc/deplete/independent_operator.py
+++ b/openmc/deplete/independent_operator.py
@@ -5,8 +5,10 @@ transport solver by using user-provided one-group cross sections.
 
 """
 
+from __future__ import annotations
+from collections.abc import Iterable
 import copy
-from itertools import product
+from typing import List, Set
 
 import numpy as np
 from uncertainties import ufloat
@@ -37,14 +39,20 @@ class IndependentOperator(OpenMCOperator):
 
     .. versionadded:: 0.13.1
 
+    .. versionchanged:: 0.13.4
+        Arguments updated to include list of fluxes and microscopic cross
+        sections.
+
     Parameters
     ----------
     materials : openmc.Materials
         Materials to deplete.
-    micro_xs : MicroXS
-        One-group microscopic cross sections in [b]. If the
-        :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will
-        be run.
+    fluxes : list of numpy.ndarray
+        Flux in each group in [n-cm/src] for each domain
+    micros : list of MicroXS
+        Cross sections in [b] for each domain. If the
+        :class:`~openmc.deplete.MicroXS` object is empty, a decay-only
+        calculation will be run.
     chain_file : str
         Path to the depletion chain XML file. Defaults to
         ``openmc.config['chain_file']``.
@@ -109,7 +117,8 @@ class IndependentOperator(OpenMCOperator):
 
     def __init__(self,
                  materials,
-                 micro_xs,
+                 fluxes,
+                 micros,
                  chain_file=None,
                  keff=None,
                  normalization_mode='fission-q',
@@ -120,7 +129,7 @@ class IndependentOperator(OpenMCOperator):
                  fission_yield_opts=None):
         # Validate micro-xs parameters
         check_type('materials', materials, openmc.Materials)
-        check_type('micro_xs', micro_xs, MicroXS)
+        check_type('micros', micros, Iterable, MicroXS)
         if keff is not None:
             check_type('keff', keff, tuple, float)
             keff = ufloat(*keff)
@@ -132,10 +141,15 @@ class IndependentOperator(OpenMCOperator):
         helper_kwargs = {'normalization_mode': normalization_mode,
                          'fission_yield_opts': fission_yield_opts}
 
-        cross_sections = micro_xs * 1e-24
+        # Sort fluxes and micros in same order that materials get sorted
+        index_sort = np.argsort([mat.id for mat in materials])
+        fluxes = [fluxes[i] for i in index_sort]
+        micros = [micros[i] for i in index_sort]
+
+        self.fluxes = fluxes
         super().__init__(
             materials,
-            cross_sections,
+            micros,
             chain_file,
             prev_results,
             fission_q=fission_q,
@@ -145,6 +159,7 @@ class IndependentOperator(OpenMCOperator):
 
     @classmethod
     def from_nuclides(cls, volume, nuclides,
+                      flux,
                       micro_xs,
                       chain_file=None,
                       nuc_units='atom/b-cm',
@@ -163,10 +178,11 @@ class IndependentOperator(OpenMCOperator):
         nuclides : dict of str to float
             Dictionary with nuclide names as keys and nuclide concentrations as
             values.
+        flux : numpy.ndarray
+            Flux in each group in [n-cm/src]
         micro_xs : MicroXS
-            One-group microscopic cross sections in [b]. If the
-            :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation
-            will be run.
+            Cross sections in [b]. If the :class:`~openmc.deplete.MicroXS`
+            object is empty, a decay-only calculation will be run.
         chain_file : str, optional
             Path to the depletion chain XML file. Defaults to
             ``openmc.config['chain_file']``.
@@ -203,8 +219,11 @@ class IndependentOperator(OpenMCOperator):
         """
         check_type('nuclides', nuclides, dict, str)
         materials = cls._consolidate_nuclides_to_material(nuclides, nuc_units, volume)
+        fluxes = [flux]
+        micros = [micro_xs]
         return cls(materials,
-                   micro_xs,
+                   fluxes,
+                   micros,
                    chain_file,
                    keff=keff,
                    normalization_mode=normalization_mode,
@@ -256,9 +275,9 @@ class IndependentOperator(OpenMCOperator):
                 self.prev_res.append(new_res)
 
 
-    def _get_nuclides_with_data(self, cross_sections):
+    def _get_nuclides_with_data(self, cross_sections: List[MicroXS]) -> Set[str]:
         """Finds nuclides with cross section data"""
-        return set(cross_sections.index)
+        return set(cross_sections[0].nuclides)
 
     class _IndependentRateHelper(ReactionRateHelper):
         """Class for generating one-group reaction rates with flux and
@@ -285,7 +304,7 @@ class IndependentOperator(OpenMCOperator):
 
         """
 
-        def __init__(self, op):
+        def __init__(self, op: IndependentOperator):
             rates = op.reaction_rates
             super().__init__(rates.n_nuc, rates.n_react)
 
@@ -301,13 +320,13 @@ class IndependentOperator(OpenMCOperator):
             """Unused in this case"""
             pass
 
-        def get_material_rates(self, mat_id, nuc_index, react_index):
+        def get_material_rates(self, mat_index, nuc_index, react_index):
             """Return 2D array of [nuclide, reaction] reaction rates
 
             Parameters
             ----------
-            mat_id : int
-                Unique ID for the requested material
+            mat_index : int
+                Index for the material
             nuc_index : list of str
                 Ordering of desired nuclides
             react_index : list of str
@@ -315,20 +334,18 @@ class IndependentOperator(OpenMCOperator):
             """
             self._results_cache.fill(0.0)
 
-            # Get volume in units of [b-cm]
-            volume_b_cm = 1e24 * self._op.number.get_mat_volume(mat_id)
+            # Get flux and microscopic cross sections from operator
+            flux = self._op.fluxes[mat_index]
+            xs = self._op.cross_sections[mat_index]
 
-            for i_nuc, i_react in product(nuc_index, react_index):
+            for i_nuc in nuc_index:
                 nuc = self.nuc_ind_map[i_nuc]
-                rx = self.rx_ind_map[i_react]
+                for i_rx in react_index:
+                    rx = self.rx_ind_map[i_rx]
 
-                # OK, this is kind of weird, but we multiply by volume in [b-cm]
-                # only because OpenMCOperator._calculate_reaction_rates has to
-                # divide it out later. It might make more sense to account for
-                # the source rate (flux) here rather than in the normalization
-                # helper.
-                self._results_cache[i_nuc, i_react] = \
-                    self._op.cross_sections[rx][nuc] * volume_b_cm
+                    # Determine reaction rate by multiplying xs in [b] by flux
+                    # in [n-cm/src] to give [(reactions/src)*b-cm/atom]
+                    self._results_cache[i_nuc, i_rx] = (xs[nuc, rx] * flux).sum()
 
             return self._results_cache
 
diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py
index 050389cc2..59285e213 100644
--- a/openmc/deplete/microxs.py
+++ b/openmc/deplete/microxs.py
@@ -1,15 +1,17 @@
 """MicroXS module
 
-A pandas.DataFrame storing microscopic cross section data with
-nuclide names as row indices and reaction names as column indices.
+A class for storing microscopic cross section data that can be used with the
+IndependentOperator class for depletion.
 """
 
+from __future__ import annotations
 import tempfile
+from typing import List, Tuple, Iterable, Optional, Union
 
-from pandas import DataFrame, read_csv, Series
+import pandas as pd
 import numpy as np
 
-from openmc.checkvalue import check_type, check_value, check_iterable_type
+from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike
 from openmc.exceptions import DataError
 from openmc import StatePoint
 import openmc
@@ -20,168 +22,178 @@ _valid_rxns = list(REACTIONS)
 _valid_rxns.append('fission')
 
 
-class MicroXS(DataFrame):
+def get_microxs_and_flux(
+        model: openmc.Model,
+        domains,
+        nuclides: Optional[Iterable[str]] = None,
+        reactions: Optional[Iterable[str]] = None,
+        energies: Optional[Union[Iterable[float], str]] = None,
+        chain_file: Optional[PathLike] = None,
+        run_kwargs=None
+    ) -> Tuple[List[np.ndarray], List[MicroXS]]:
+    """Generate a microscopic cross sections and flux from a Model
+
+    .. versionadded:: 0.13.4
+
+    Parameters
+    ----------
+    model : openmc.Model
+        OpenMC model object. Must contain geometry, materials, and settings.
+    domains : list of openmc.Material or openmc.Cell or openmc.Universe
+        Domains in which to tally reaction rates.
+    nuclides : list of str
+        Nuclides to get cross sections for. If not specified, all burnable
+        nuclides from the depletion chain file are used.
+    reactions : list of str
+        Reactions to get cross sections for. If not specified, all neutron
+        reactions listed in the depletion chain file are used.
+    energies : iterable of float or str
+        Energy group boundaries in [eV] or the name of the group structure
+    chain_file : str, optional
+        Path to the depletion chain XML file that will be used in depletion
+        simulation. Used to determine cross sections for materials not
+        present in the inital composition. Defaults to
+        ``openmc.config['chain_file']``.
+    run_kwargs : dict, optional
+        Keyword arguments passed to :meth:`openmc.Model.run`
+
+    Returns
+    -------
+    list of numpy.ndarray
+        Flux in each group in [n-cm/src] for each domain
+    list of MicroXS
+        Cross section data in [b] for each domain
+
+    """
+    # Save any original tallies on the model
+    original_tallies = model.tallies
+
+    # Determine what reactions and nuclides are available in chain
+    if chain_file is None:
+        chain_file = openmc.config.get('chain_file')
+        if chain_file is None:
+            raise DataError(
+                "No depletion chain specified and could not find depletion "
+                "chain in openmc.config['chain_file']"
+            )
+    chain = Chain.from_xml(chain_file)
+    if reactions is None:
+        reactions = chain.reactions
+    if not nuclides:
+        cross_sections = _find_cross_sections(model)
+        nuclides_with_data = _get_nuclides_with_data(cross_sections)
+        nuclides = [nuc.name for nuc in chain.nuclides
+                    if nuc.name in nuclides_with_data]
+
+    # Set up the reaction rate and flux tallies
+    if energies is None:
+        energy_filter = openmc.EnergyFilter([0.0, 100.0e6])
+    elif isinstance(energies, str):
+        energy_filter = openmc.EnergyFilter.from_group_structure(energies)
+    else:
+        energy_filter = openmc.EnergyFilter(energies)
+    if isinstance(domains[0], openmc.Material):
+        domain_filter = openmc.MaterialFilter(domains)
+    elif isinstance(domains[0], openmc.Cell):
+        domain_filter = openmc.CellFilter(domains)
+    elif isinstance(domains[0], openmc.Universe):
+        domain_filter = openmc.UniverseFilter(domains)
+    else:
+        raise ValueError(f"Unsupported domain type: {type(domains[0])}")
+
+    rr_tally = openmc.Tally(name='MicroXS RR')
+    rr_tally.filters = [domain_filter, energy_filter]
+    rr_tally.nuclides = nuclides
+    rr_tally.multiply_density = False
+    rr_tally.scores = reactions
+
+    flux_tally = openmc.Tally(name='MicroXS flux')
+    flux_tally.filters = [domain_filter, energy_filter]
+    flux_tally.scores = ['flux']
+    model.tallies = openmc.Tallies([rr_tally, flux_tally])
+
+    # create temporary run
+    with tempfile.TemporaryDirectory() as temp_dir:
+        if run_kwargs is None:
+            run_kwargs = {}
+        else:
+            run_kwargs = dict(run_kwargs)
+        run_kwargs.setdefault('cwd', temp_dir)
+        statepoint_path = model.run(**run_kwargs)
+
+        with StatePoint(statepoint_path) as sp:
+            rr_tally = sp.tallies[rr_tally.id]
+            rr_tally._read_results()
+            flux_tally = sp.tallies[flux_tally.id]
+            flux_tally._read_results()
+
+    # Get reaction rates and flux values
+    reaction_rates = rr_tally.get_reshaped_data()  # (domains, groups, nuclides, reactions)
+    flux = flux_tally.get_reshaped_data()  # (domains, groups, 1, 1)
+
+    # Make energy groups last dimension
+    reaction_rates = np.moveaxis(reaction_rates, 1, -1)  # (domains, nuclides, reactions, groups)
+    flux = np.moveaxis(flux, 1, -1)  # (domains, 1, 1, groups)
+
+    # Divide RR by flux to get microscopic cross sections
+    xs = np.empty_like(reaction_rates) # (domains, nuclides, reactions, groups)
+    d, _, _, g = np.nonzero(flux)
+    xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
+
+    # Reset tallies
+    model.tallies = original_tallies
+
+    # Create lists where each item corresponds to one domain
+    fluxes = list(flux.squeeze((1, 2)))
+    micros = [MicroXS(xs_i, nuclides, reactions) for xs_i in xs]
+    return fluxes, micros
+
+
+class MicroXS:
     """Microscopic cross section data for use in transport-independent depletion.
 
     .. versionadded:: 0.13.1
 
+    .. versionchanged:: 0.13.4
+        Class was heavily refactored and no longer subclasses pandas.DataFrame.
+
+    Parameters
+    ----------
+    data : numpy.ndarray of floats
+        3D array containing microscopic cross section values for each
+        nuclide, reaction, and energy group. Cross section values are assumed to
+        be in [b], and indexed by [nuclide, reaction, energy group]
+    nuclides : list of str
+        List of nuclide symbols for that have data for at least one
+        reaction.
+    reactions : list of str
+        List of reactions. All reactions must match those in
+        :data:`openmc.deplete.chain.REACTIONS`
+
     """
-
-    @classmethod
-    def from_model(cls,
-                   model,
-                   domain,
-                   nuclides=None,
-                   reactions=None,
-                   chain_file=None,
-                   energy_bounds=(0, 20e6),
-                   run_kwargs=None):
-        """Generate a one-group cross-section dataframe using OpenMC.
-
-        Note that the ``openmc`` executable must be compiled.
-
-        Parameters
-        ----------
-        model : openmc.Model
-            OpenMC model object. Must contain geometry, materials, and settings.
-        domain : openmc.Material or openmc.Cell or openmc.Universe
-            Domain in which to tally reaction rates.
-        nuclides : list of str
-            Nuclides to get cross sections for. If not specified, all burnable
-            nuclides from the depletion chain file are used.
-        reactions : list of str
-            Reactions to get cross sections for. If not specified, all neutron
-            reactions listed in the depletion chain file are used.
-        chain_file : str, optional
-            Path to the depletion chain XML file that will be used in depletion
-            simulation. Used to determine cross sections for materials not
-            present in the inital composition. Defaults to
-            ``openmc.config['chain_file']``.
-        energy_bound : 2-tuple of float, optional
-            Bounds for the energy group.
-        run_kwargs : dict, optional
-            Keyword arguments passed to :meth:`openmc.model.Model.run`
-
-        Returns
-        -------
-        MicroXS
-            Cross section data in [b]
-
-        """
-        # Save any original tallies on the model
-        original_tallies = model.tallies
-
-        # Determine what reactions and nuclides are available in chain
-        if chain_file is None:
-            chain_file = openmc.config.get('chain_file')
-            if chain_file is None:
-                raise DataError(
-                    "No depletion chain specified and could not find depletion "
-                    "chain in openmc.config['chain_file']"
-                )
-        chain = Chain.from_xml(chain_file)
-        if reactions is None:
-            reactions = chain.reactions
-        if not nuclides:
-            cross_sections = _find_cross_sections(model)
-            nuclides_with_data = _get_nuclides_with_data(cross_sections)
-            nuclides = [nuc.name for nuc in chain.nuclides
-                        if nuc.name in nuclides_with_data]
-
-        # Set up the reaction rate and flux tallies
-        energy_filter = openmc.EnergyFilter(energy_bounds)
-        if isinstance(domain, openmc.Material):
-            domain_filter = openmc.MaterialFilter([domain])
-        elif isinstance(domain, openmc.Cell):
-            domain_filter = openmc.CellFilter([domain])
-        elif isinstance(domain, openmc.Universe):
-            domain_filter = openmc.UniverseFilter([domain])
-        else:
-            raise ValueError(f"Unsupported domain type: {type(domain)}")
-
-        rr_tally = openmc.Tally(name='MicroXS RR')
-        rr_tally.filters = [domain_filter, energy_filter]
-        rr_tally.nuclides = nuclides
-        rr_tally.multiply_density = False
-        rr_tally.scores = reactions
-
-        flux_tally = openmc.Tally(name='MicroXS flux')
-        flux_tally.filters = [domain_filter, energy_filter]
-        flux_tally.scores = ['flux']
-        tallies = openmc.Tallies([rr_tally, flux_tally])
-
-        model.tallies = tallies
-
-        # create temporary run
-        with tempfile.TemporaryDirectory() as temp_dir:
-            if run_kwargs is None:
-                run_kwargs = {}
-            run_kwargs.setdefault('cwd', temp_dir)
-            statepoint_path = model.run(**run_kwargs)
-
-            with StatePoint(statepoint_path) as sp:
-                rr_tally = sp.tallies[rr_tally.id]
-                rr_tally._read_results()
-                flux_tally = sp.tallies[flux_tally.id]
-                flux_tally._read_results()
-
-        # Get reaction rates and flux values
-        reaction_rates = rr_tally.mean.sum(axis=0)  # (nuclides, reactions)
-        flux = flux_tally.mean[0, 0, 0]
-
-        # Divide RR by flux to get microscopic cross sections
-        xs = reaction_rates / flux
-
-        # Build Series objects
-        series = {}
-        for i, rx in enumerate(reactions):
-            series[rx] = Series(xs[..., i], index=rr_tally.nuclides)
-
-        # Revert to the original tallies and materials
-        model.tallies = original_tallies
-
-        return cls(series).rename_axis('nuclide')
-
-    @classmethod
-    def from_array(cls, nuclides, reactions, data):
-        """
-        Creates a ``MicroXS`` object from arrays.
-
-        Parameters
-        ----------
-        nuclides : list of str
-            List of nuclide symbols for that have data for at least one
-            reaction.
-        reactions : list of str
-            List of reactions. All reactions must match those in
-            :data:`openmc.deplete.chain.REACTIONS`
-        data : ndarray of floats
-            Array containing one-group microscopic cross section values for
-            each nuclide and reaction. Cross section values are assumed to be
-            in [b].
-
-        Returns
-        -------
-        MicroXS
-        """
-
+    def __init__(self, data: np.ndarray, nuclides: List[str], reactions: List[str]):
         # Validate inputs
-        if data.shape != (len(nuclides), len(reactions)):
+        if data.shape[:2] != (len(nuclides), len(reactions)):
             raise ValueError(
                 f'Nuclides list of length {len(nuclides)} and '
                 f'reactions array of length {len(reactions)} do not '
                 f'match dimensions of data array of shape {data.shape}')
+        check_iterable_type('nuclides', nuclides, str)
+        check_iterable_type('reactions', reactions, str)
+        check_type('data', data, np.ndarray, expected_iter_type=float)
+        for reaction in reactions:
+            check_value('reactions', reaction, _valid_rxns)
 
-        cls._validate_micro_xs_inputs(
-            nuclides, reactions, data)
-        micro_xs = cls(index=nuclides, columns=reactions, data=data)
+        self.data = data
+        self.nuclides = nuclides
+        self.reactions = reactions
 
-        return micro_xs
+    # TODO: Add a classmethod for generating MicroXS directly from cross section
+    # data using a known flux spectrum
 
     @classmethod
     def from_csv(cls, csv_file, **kwargs):
-        """
-        Load a ``MicroXS`` object from a ``.csv`` file.
+        """Load data from a comma-separated values (csv) file.
 
         Parameters
         ----------
@@ -199,17 +211,37 @@ class MicroXS(DataFrame):
         if 'float_precision' not in kwargs:
             kwargs['float_precision'] = 'round_trip'
 
-        micro_xs = cls(read_csv(csv_file, index_col=0, **kwargs))
+        df = pd.read_csv(csv_file, **kwargs)
+        df.set_index(['nuclides', 'reactions', 'groups'], inplace=True)
+        nuclides = list(df.index.unique(level='nuclides'))
+        reactions = list(df.index.unique(level='reactions'))
+        groups = list(df.index.unique(level='groups'))
+        shape = (len(nuclides), len(reactions), len(groups))
+        data = df.values.reshape(shape)
+        return cls(data, nuclides, reactions)
 
-        cls._validate_micro_xs_inputs(list(micro_xs.index),
-                                      list(micro_xs.columns),
-                                      micro_xs.to_numpy())
-        return micro_xs
+    def __getitem__(self, index):
+        nuc, rx = index
+        i_nuc = self.nuclides.index(nuc)
+        i_rx = self.reactions.index(rx)
+        return self.data[i_nuc, i_rx]
+
+    def to_csv(self, *args, **kwargs):
+        """Write data to a comma-separated values (csv) file
+
+        Parameters
+        ----------
+        *args
+            Positional arguments passed to :meth:`pandas.DataFrame.to_csv`
+        **kwargs
+            Keyword arguments passed to :meth:`pandas.DataFrame.to_csv`
+
+        """
+        groups = self.data.shape[2]
+        multi_index = pd.MultiIndex.from_product(
+            [self.nuclides, self.reactions, range(1, groups + 1)],
+            names=['nuclides', 'reactions', 'groups']
+        )
+        df = pd.DataFrame({'xs': self.data.flatten()}, index=multi_index)
+        df.to_csv(*args, **kwargs)
 
-    @staticmethod
-    def _validate_micro_xs_inputs(nuclides, reactions, data):
-        check_iterable_type('nuclides', nuclides, str)
-        check_iterable_type('reactions', reactions, str)
-        check_type('data', data, np.ndarray, expected_iter_type=float)
-        for reaction in reactions:
-            check_value('reactions', reaction, _valid_rxns)
diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py
index d3b9ef911..319d98701 100644
--- a/openmc/deplete/openmc_operator.py
+++ b/openmc/deplete/openmc_operator.py
@@ -6,6 +6,7 @@ transport-independent transport operators.
 """
 
 from abc import abstractmethod
+from typing import List, Tuple, Dict
 
 import numpy as np
 
@@ -31,9 +32,9 @@ class OpenMCOperator(TransportOperator):
     ----------
     materials : openmc.Materials
         List of all materials in the model
-    cross_sections : str or pandas.DataFrame
-        Path to continuous energy cross section library, or object containing
-        one-group cross-sections.
+    cross_sections : str or list of MicroXS
+        Path to continuous energy cross section library, or list of objects
+        containing cross sections.
     chain_file : str, optional
         Path to the depletion chain XML file. Defaults to
         openmc.config['chain_file'].
@@ -60,9 +61,9 @@ class OpenMCOperator(TransportOperator):
     ----------
     materials : openmc.Materials
         All materials present in the model
-    cross_sections : str or MicroXS
-            Path to continuous energy cross section library, or object
-            containing one-group cross-sections.
+    cross_sections : str or list of MicroXS
+        Path to continuous energy cross section library, or list of objects
+        containing cross sections.
     output_dir : pathlib.Path
         Path to output directory to save results.
     round_number : bool
@@ -164,7 +165,7 @@ class OpenMCOperator(TransportOperator):
         """Assign distribmats for each burnable material"""
         pass
 
-    def _get_burnable_mats(self):
+    def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]:
         """Determine depletable materials, volumes, and nuclides
 
         Returns
diff --git a/openmc/deplete/reaction_rates.py b/openmc/deplete/reaction_rates.py
index b806dfad8..4aaf829a9 100644
--- a/openmc/deplete/reaction_rates.py
+++ b/openmc/deplete/reaction_rates.py
@@ -2,6 +2,8 @@
 
 An ndarray to store reaction rates with string, integer, or slice indexing.
 """
+from typing import Dict
+
 import numpy as np
 
 
@@ -51,6 +53,10 @@ class ReactionRates(np.ndarray):
     # the __array_finalize__ method (discussed here:
     # https://docs.scipy.org/doc/numpy/user/basics.subclassing.html)
 
+    index_mat: Dict[str, int]
+    index_nuc: Dict[str, int]
+    index_rx: Dict[str, int]
+
     def __new__(cls, local_mats, nuclides, reactions, from_results=False):
         # Create appropriately-sized zeroed-out ndarray
         shape = (len(local_mats), len(nuclides), len(reactions))
diff --git a/tests/micro_xs_simple.csv b/tests/micro_xs_simple.csv
index 71981c718..146896aa2 100644
--- a/tests/micro_xs_simple.csv
+++ b/tests/micro_xs_simple.csv
@@ -1,13 +1,25 @@
-nuclide,"(n,gamma)",fission
-U234,22.231989822002454,0.4962074466374984
-U235,10.479008971197121,48.41787337164606
-U238,0.8673334105437321,0.1046788058876236
-U236,8.651710446071224,0.31948392400019293
-O16,7.497851000107522e-05,0.0
-O17,0.0004079227797153372,0.0
-I135,6.842395323713929,0.0
-Xe135,227463.8642699061,0.0
-Xe136,0.023178960347535887,0.0
-Cs135,2.1721665580713623,0.0
-Gd157,12786.099392370175,0.0
-Gd156,3.4006085445846983,0.0
+nuclides,reactions,groups,xs
+U234,"(n,gamma)",1,22.23198982200245
+U234,fission,1,0.4962074466374984
+U235,"(n,gamma)",1,10.47900897119712
+U235,fission,1,48.41787337164606
+U238,"(n,gamma)",1,0.8673334105437321
+U238,fission,1,0.1046788058876236
+U236,"(n,gamma)",1,8.651710446071224
+U236,fission,1,0.3194839240001929
+O16,"(n,gamma)",1,7.497851000107522e-05
+O16,fission,1,0.0
+O17,"(n,gamma)",1,0.0004079227797153
+O17,fission,1,0.0
+I135,"(n,gamma)",1,6.842395323713929
+I135,fission,1,0.0
+Xe135,"(n,gamma)",1,227463.8642699061
+Xe135,fission,1,0.0
+Xe136,"(n,gamma)",1,0.0231789603475358
+Xe136,fission,1,0.0
+Cs135,"(n,gamma)",1,2.1721665580713623
+Cs135,fission,1,0.0
+Gd157,"(n,gamma)",1,12786.099392370175
+Gd157,fission,1,0.0
+Gd156,"(n,gamma)",1,3.4006085445846983
+Gd156,fission,1,0.0
diff --git a/tests/regression_tests/deplete_no_transport/test.py b/tests/regression_tests/deplete_no_transport/test.py
index 64769c0b0..54c29cd5b 100644
--- a/tests/regression_tests/deplete_no_transport/test.py
+++ b/tests/regression_tests/deplete_no_transport/test.py
@@ -36,13 +36,16 @@ def chain_file():
     return Path(__file__).parents[2] / 'chain_simple.xml'
 
 
-@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, flux", [
-    (True, True, 'source-rate', None, 1164719970082145.0),
-    (False, True, 'source-rate', None, 1164719970082145.0),
+neutron_per_cm2_sec = 1164719970082145.0
+
+
+@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, source_rate", [
+    (True, True, 'source-rate', None, 1.0),
+    (False, True, 'source-rate', None, 1.0),
     (True, True, 'fission-q', 174, None),
     (False, True, 'fission-q', 174, None),
-    (True, False, 'source-rate', None, 1164719970082145.0),
-    (False, False, 'source-rate', None, 1164719970082145.0),
+    (True, False, 'source-rate', None, 1.0),
+    (False, False, 'source-rate', None, 1.0),
     (True, False, 'fission-q', 174, None),
     (False, False, 'fission-q', 174, None)])
 def test_against_self(run_in_tmpdir,
@@ -53,7 +56,7 @@ def test_against_self(run_in_tmpdir,
                       from_nuclides,
                       normalization_mode,
                       power,
-                      flux):
+                      source_rate):
     """Transport free system test suite.
 
     Runs an OpenMC transport-free depletion calculation and verifies
@@ -61,8 +64,10 @@ def test_against_self(run_in_tmpdir,
 
     """
     # Create operator
+    flux = neutron_per_cm2_sec * fuel.volume
     op = _create_operator(from_nuclides,
                           fuel,
+                          flux,
                           micro_xs,
                           chain_file,
                           normalization_mode)
@@ -75,12 +80,12 @@ def test_against_self(run_in_tmpdir,
     openmc.deplete.PredictorIntegrator(op,
                                        dt,
                                        power=power,
-                                       source_rates=flux,
+                                       source_rates=source_rate,
                                        timestep_units='s').integrate()
 
     # Get path to test and reference results
     path_test = op.output_dir / 'depletion_results.h5'
-    if flux is not None:
+    if power is None:
         ref_path = 'test_reference_source_rate.h5'
     else:
         ref_path = 'test_reference_fission_q.h5'
@@ -99,8 +104,8 @@ def test_against_self(run_in_tmpdir,
     _assert_same_mats(res_test, res_ref)
 
     tol = 1.0e-14
-    assert_atoms_equal(res_test, res_ref, tol)
-    assert_reaction_rates_equal(res_test, res_ref, tol)
+    assert_atoms_equal(res_ref, res_test, tol)
+    assert_reaction_rates_equal(res_ref, res_test, tol)
 
 
 @pytest.mark.parametrize("multiproc, dt, time_units, time_type, atom_tol, rx_tol ", [
@@ -123,7 +128,8 @@ def test_against_coupled(run_in_tmpdir,
                          atom_tol,
                          rx_tol):
     # Create operator
-    op = _create_operator(False, fuel, micro_xs, chain_file, 'fission-q')
+    flux = neutron_per_cm2_sec * fuel.volume
+    op = _create_operator(False, fuel, flux, micro_xs, chain_file, 'fission-q')
 
     # Power and timesteps
     dt = [dt]  # single step
@@ -151,12 +157,13 @@ def test_against_coupled(run_in_tmpdir,
     # Assert same mats
     _assert_same_mats(res_test, res_ref)
 
-    assert_atoms_equal(res_test, res_ref, atom_tol)
-    assert_reaction_rates_equal(res_test, res_ref, rx_tol)
+    assert_atoms_equal(res_ref, res_test, atom_tol)
+    assert_reaction_rates_equal(res_ref, res_test, rx_tol)
 
 
 def _create_operator(from_nuclides,
                      fuel,
+                     flux,
                      micro_xs,
                      chain_file,
                      normalization_mode):
@@ -167,13 +174,15 @@ def _create_operator(from_nuclides,
 
         op = IndependentOperator.from_nuclides(fuel.volume,
                                                nuclides,
+                                               flux,
                                                micro_xs,
                                                chain_file,
                                                normalization_mode=normalization_mode)
 
     else:
         op = IndependentOperator(openmc.Materials([fuel]),
-                                 micro_xs,
+                                 [flux],
+                                 [micro_xs],
                                  chain_file,
                                  normalization_mode=normalization_mode)
 
diff --git a/tests/regression_tests/microxs/test.py b/tests/regression_tests/microxs/test.py
index 3d67ee8a3..f4b3649fd 100644
--- a/tests/regression_tests/microxs/test.py
+++ b/tests/regression_tests/microxs/test.py
@@ -4,7 +4,7 @@ from pathlib import Path
 import numpy as np
 import pytest
 import openmc
-from openmc.deplete import MicroXS
+from openmc.deplete import MicroXS, get_microxs_and_flux
 
 from tests.regression_tests import config
 
@@ -47,13 +47,13 @@ def model():
 
 
 def test_from_model(model):
-    fuel = model.materials[0]
+    domains = model.materials[:1]
     nuclides = ['U234', 'U235', 'U238', 'U236', 'O16', 'O17', 'I135', 'Xe135',
                 'Xe136', 'Cs135', 'Gd157', 'Gd156']
-    test_xs = MicroXS.from_model(model, fuel, nuclides, chain_file=CHAIN_FILE)
+    _, test_xs = get_microxs_and_flux(model, domains, nuclides, chain_file=CHAIN_FILE)
     if config['update']:
-        test_xs.to_csv('test_reference.csv')
+        test_xs[0].to_csv('test_reference.csv')
 
     ref_xs = MicroXS.from_csv('test_reference.csv')
 
-    np.testing.assert_allclose(test_xs, ref_xs, rtol=1e-11)
+    np.testing.assert_allclose(test_xs[0].data, ref_xs.data, rtol=1e-11)
diff --git a/tests/regression_tests/microxs/test_reference.csv b/tests/regression_tests/microxs/test_reference.csv
index 0fc41862d..3e679f980 100644
--- a/tests/regression_tests/microxs/test_reference.csv
+++ b/tests/regression_tests/microxs/test_reference.csv
@@ -1,13 +1,25 @@
-nuclide,"(n,gamma)",fission
-U234,21.418670317831197,0.5014588470882195
-U235,10.343944102483244,47.46718472611891
-U238,0.8741166723597251,0.10829568455139126
-U236,9.083486784689326,0.3325287927011428
-O16,7.548646353912453e-05,0.0
-O17,0.0004018486221310307,0.0
-I135,6.6912565089429235,0.0
-Xe135,223998.64185667288,0.0
-Xe136,0.022934362666193576,0.0
-Cs135,2.28453952223533,0.0
-Gd157,12582.079620036275,0.0
-Gd156,2.9421127515332417,0.0
+nuclides,reactions,groups,xs
+U234,"(n,gamma)",1,21.418670317831076
+U234,fission,1,0.5014588470882162
+U235,"(n,gamma)",1,10.343944102483215
+U235,fission,1,47.46718472611895
+U238,"(n,gamma)",1,0.8741166723597229
+U238,fission,1,0.10829568455139067
+U236,"(n,gamma)",1,9.08348678468935
+U236,fission,1,0.3325287927011424
+O16,"(n,gamma)",1,7.548646353912426e-05
+O16,fission,1,0.0
+O17,"(n,gamma)",1,0.00040184862213103105
+O17,fission,1,0.0
+I135,"(n,gamma)",1,6.691256508942912
+I135,fission,1,0.0
+Xe135,"(n,gamma)",1,223998.6418566729
+Xe135,fission,1,0.0
+Xe136,"(n,gamma)",1,0.022934362666193503
+Xe136,fission,1,0.0
+Cs135,"(n,gamma)",1,2.2845395222353204
+Cs135,fission,1,0.0
+Gd157,"(n,gamma)",1,12582.07962003624
+Gd157,fission,1,0.0
+Gd156,"(n,gamma)",1,2.942112751533234
+Gd156,fission,1,0.0
diff --git a/tests/unit_tests/test_deplete_decay_products.py b/tests/unit_tests/test_deplete_decay_products.py
index 34ca1b067..751a96ca6 100644
--- a/tests/unit_tests/test_deplete_decay_products.py
+++ b/tests/unit_tests/test_deplete_decay_products.py
@@ -1,4 +1,5 @@
 import openmc.deplete
+import numpy as np
 import pytest
 
 
@@ -16,11 +17,14 @@ def test_deplete_decay_products(run_in_tmpdir):
 </depletion_chain>
         """)
 
+    # Create MicroXS object with no cross sections
+    micro_xs = openmc.deplete.MicroXS(np.empty((0, 0)), [], [])
+
     # Create depletion operator with no reactions
-    micro_xs = openmc.deplete.MicroXS()
     op = openmc.deplete.IndependentOperator.from_nuclides(
         volume=1.0,
         nuclides={'Li5': 1.0},
+        flux=0.0,
         micro_xs=micro_xs,
         chain_file='test_chain.xml',
         normalization_mode='source-rate'
diff --git a/tests/unit_tests/test_deplete_independent_operator.py b/tests/unit_tests/test_deplete_independent_operator.py
index f6cc7e200..813ac06de 100644
--- a/tests/unit_tests/test_deplete_independent_operator.py
+++ b/tests/unit_tests/test_deplete_independent_operator.py
@@ -4,14 +4,13 @@
 
 from pathlib import Path
 
-import pytest
 from openmc.deplete import IndependentOperator, MicroXS
 from openmc import Material, Materials
-import numpy as np
 
 CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
 ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv"
 
+
 def test_operator_init():
     """The test uses a temporary dummy chain. This file will be removed
     at the end of the test, and only contains a depletion_chain node."""
@@ -22,15 +21,18 @@ def test_operator_init():
                 'U236': 4.5724195495061115e+18,
                 'O16': 4.639065406771322e+22,
                 'O17': 1.7588724018066158e+19}
+    flux = 1.0
     micro_xs = MicroXS.from_csv(ONE_GROUP_XS)
     IndependentOperator.from_nuclides(
-        volume, nuclides, micro_xs, CHAIN_PATH, nuc_units='atom/cm3')
+        volume, nuclides, flux, micro_xs, CHAIN_PATH, nuc_units='atom/cm3')
 
     fuel = Material(name="uo2")
     fuel.add_element("U", 1, percent_type="ao", enrichment=4.25)
     fuel.add_element("O", 2)
     fuel.set_density("g/cc", 10.4)
-    fuel.depletable=True
+    fuel.depletable = True
     fuel.volume = 1
     materials = Materials([fuel])
-    IndependentOperator(materials, micro_xs, CHAIN_PATH)
+    fluxes = [1.0]
+    micros = [micro_xs]
+    IndependentOperator(materials, fluxes, micros, CHAIN_PATH)
diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py
index 557082eda..582c58465 100644
--- a/tests/unit_tests/test_deplete_microxs.py
+++ b/tests/unit_tests/test_deplete_microxs.py
@@ -43,17 +43,18 @@ def test_from_array():
                      [0., 0.],
                      [0., 0.1],
                      [0., 0.1]])
+    data.shape = (12, 2, 1)
 
-    MicroXS.from_array(nuclides, reactions, data)
+    MicroXS(data, nuclides, reactions)
     with pytest.raises(ValueError, match=r'Nuclides list of length \d* and '
                        r'reactions array of length \d* do not '
-                       r'match dimensions of data array of shape \(\d*\,d*\)'):
-        MicroXS.from_array(nuclides, reactions, data[:, 0])
+                       r'match dimensions of data array of shape \(\d*\, \d*\)'):
+        MicroXS(data[:, 0], nuclides, reactions)
 
 def test_csv():
     ref_xs = MicroXS.from_csv(ONE_GROUP_XS)
     ref_xs.to_csv('temp_xs.csv')
     temp_xs = MicroXS.from_csv('temp_xs.csv')
-    assert np.all(ref_xs == temp_xs)
+    assert np.all(ref_xs.data == temp_xs.data)
     remove('temp_xs.csv')