diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 8bbe126df4..c00c80399f 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -329,7 +329,7 @@ class NormalizationHelper(ABC):
             `fission_rates` for :meth:`update`.
         """
 
-    def update(self, fission_rates):
+    def update(self, fission_rates, mat_index=None):
         """Update the normalization based on fission rates (only used for
         energy-based normalization)
 
@@ -339,6 +339,8 @@ class NormalizationHelper(ABC):
             fission reaction rate for each isotope in the specified
             material. Should be ordered corresponding to initial
             ``rate_index`` used in :meth:`prepare`
+        mat_index : int
+            Material index
         """
 
     @property
@@ -523,6 +525,8 @@ class Integrator(ABC):
         Power density of the reactor in [W/gHM]. It is multiplied by
         initial heavy metal inventory to get total power if ``power``
         is not speficied.
+    flux : float or iterable of float, optional
+        Neutron flux in [neur/s-cm^2] for each interval in :attr: `timesteps`
     source_rates : float or iterable of float, optional
         Source rate in [neutron/sec] for each interval in :attr:`timesteps`
 
@@ -597,8 +601,10 @@ class Integrator(ABC):
                 source_rates = power_density * operator.heavy_metal
             else:
                 source_rates = [p*operator.heavy_metal for p in power_density]
+        elif flux is not None:
+            source_rates = flux
         elif source_rates is None:
-            raise ValueError("Either power, power_density, or source_rates must be set")
+            raise ValueError("Either power, power_density, flux, or source_rates must be set")
 
         if not isinstance(source_rates, Iterable):
             # Ensure that rate is single value if that is the case
diff --git a/openmc/deplete/atom_number.py b/openmc/deplete/atom_number.py
index 5430ddc35d..78ceecca66 100644
--- a/openmc/deplete/atom_number.py
+++ b/openmc/deplete/atom_number.py
@@ -126,6 +126,23 @@ class AtomNumber:
         return [nuc for nuc, ind in self.index_nuc.items()
                 if ind < self.n_nuc_burn]
 
+    def get_mat_volume(self, mat):
+        """Return material volume
+
+        Parameters
+        ----------
+        mat : str, int, openmc.Material, or slice
+            Material index.
+
+        Returns
+        -------
+        float
+            Material volume in [cm^3]
+
+        """
+        mat = self._get_mat_index(mat)
+        return self.volume[mat]
+
     def get_atom_density(self, mat, nuc):
         """Return atom density of given material and nuclide
 
diff --git a/openmc/deplete/flux_operator.py b/openmc/deplete/flux_operator.py
index cfcec242a2..9ff0dee5d0 100644
--- a/openmc/deplete/flux_operator.py
+++ b/openmc/deplete/flux_operator.py
@@ -54,6 +54,15 @@ class FluxDepletionOperator(OpenMCOperator):
         values will be pulled from the ``chain_file``.
     prev_results : Results, optional
         Results from a previous depletion calculation.
+    normalization_mode : {"constant-power", "constant-flux"}
+        Indicate how reaction rates should be calculated.
+        ``"constant-power"`` uses the fission Q values from the depletion chain to
+        compute the flux based on the power. ``"constant-flux"`` uses the value stored in `_normalization_helper` as the flux.
+    fission_q : dict, optional
+        Dictionary of nuclides and their fission Q values [eV]. If not given,
+        values will be pulled from the ``chain_file``. Only applicable
+        if ``"normalization_mode" == "constant-power"``.
+
     reduce_chain : bool, optional
         If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
         depletion chain up to ``reduce_chain_level``. Default is False.
@@ -108,6 +117,7 @@ class FluxDepletionOperator(OpenMCOperator):
                  flux,
                  chain_file,
                  keff=None,
+                 normalization_mode = 'constant-flux',
                  fission_q=None,
                  prev_results=None,
                  reduce_chain=False,
@@ -124,7 +134,8 @@ class FluxDepletionOperator(OpenMCOperator):
         self.flux = flux
 
         helper_kwargs = dict()
-        helper_kwargs = {'fission_yield_opts': fission_yield_opts}
+        helper_kwargs = {'normalization_mode': normalization_mode,
+                         'fission_yield_opts': fission_yield_opts}
 
         super().__init__(
             materials,
@@ -152,7 +163,7 @@ class FluxDepletionOperator(OpenMCOperator):
         volume : float
             Volume of the material being depleted in [cm^3]
         nuclides : dict of str to float
-            Dictionary with nuclide names as keys and nuclide concentrations as
+            ,Dictionary with nuclide names as keys and nuclide concentrations as
             values. Nuclide concentration units are [atom/cm^3].
         micro_xs : pandas.DataFrame
             DataFrame with nuclides names as index and microscopic cross section
@@ -164,9 +175,15 @@ class FluxDepletionOperator(OpenMCOperator):
         keff : 2-tuple of float, optional
            keff eigenvalue and uncertainty from transport calculation.
            Default is None.
+        normalization_mode : {"constant-power", "constant-flux"}
+            Indicate how reaction rates should be calculated.
+            ``"constant-power"`` uses the fission Q values from the depletion
+            chain to compute the flux based on the power. ``"constant-flux"``
+            uses the value stored in `_normalization_helper` as the flux.
         fission_q : dict, optional
-            Dictionary of nuclides and their fission Q values [eV]. If not given,
-            values will be pulled from the ``chain_file``.
+            Dictionary of nuclides and their fission Q values [eV]. If not
+            given, values will be pulled from the ``chain_file``. Only
+            applicable if ``"normalization_mode" == "constant-power"``.
         prev_results : Results, optional
             Results from a previous depletion calculation.
         reduce_chain : bool, optional
@@ -215,6 +232,49 @@ class FluxDepletionOperator(OpenMCOperator):
         """Finds nuclides with cross section data"""
         return set(cross_sections.index)
 
+    class _FluxDepletionNormalizationHelper(ChainFissionHelper):
+        """Class for calculating one-group flux
+        based on a power.
+
+        flux = Power / X, where X = volume * sum_i(Q_i * fission_micro_xs_i * density_i)
+
+        Parameters
+        ----------
+        op : openmc.deplete.FluxDepletionOperator
+            Reference to the object encapsulate _FluxDepletionNormalizationHelper.
+            We pass this so we don't have to duplicate the ``number`` object.
+
+        """
+
+        def __init__(self, op):
+            self._op = op
+            rates = self.reaction_rates
+            self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()}
+            super().__init__()
+
+        def update(self, fission_rates, mat_index=None):
+        """Update 'energy' produced with fission rates in a material. What this
+        actually calculates is the quantity X.
+
+        Parameters
+        ----------
+        fission_rates : numpy.ndarray
+            fission reaction rate for each isotope in the specified
+            material. Should be ordered corresponding to initial
+            ``rate_index`` used in :meth:`prepare`
+        mat_index : int
+            Material index
+
+        """
+            volume = self._op.number.get_mat_volume(mat_index)
+            densities = np.empty(shape(fission_rates))
+            for i_nuc in nuc_ind_map:
+                nuc = self.nuc_ind_map[i_nuc]
+                densities[i_nuc] = self._op.number.get_atom_density(mat_index, nuc)
+            fission_rates = fission_rates * volume * densities
+
+            super.update(fission_rates)
+
     class _FluxDepletionRateHelper(ReactionRateHelper):
         """Class for generating one-group reaction rates with flux and
         one-group cross sections.
@@ -246,11 +306,11 @@ class FluxDepletionOperator(OpenMCOperator):
 
         def __init__(self, n_nuc, n_react, op, nuc_ind_map. rxn_ind_map)
             super().__init__(n_nuc, n_react)
-            self._op = op
-            rates = self.reaction_rates
-            # Get classes to assit working with tallies
+            rates = op.reaction_rates
+
             self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()}
             self.rxn_ind_map = {ind: rxn for rxn, ind in rates.index_rx.items()}
+            self._op = op
 
         def generate_tallies(self, materials, scores):
             """Unused in this case"""
@@ -269,12 +329,14 @@ class FluxDepletionOperator(OpenMCOperator):
                 Ordering of reactions
             """
             self._results_cache.fill(0.0)
+
+            volume = self._op.number.get_mat_volume(mat_id)
             for i_nuc, i_react in product(nuc_index, react_index):
                 nuc = self.nuc_ind_map[i_nuc]
                 rxn = self.rxn_ind_map[i_react]
                 density = self._op.number.get_atom_density(mat_id, nuc)
                 self._results_cache[i_nuc,
-                                    i_react] = self._op.cross_sections[rxn][nuc] * density
+                                    i_react] = self._op.cross_sections[rxn][nuc] * density * volume
 
             return self._results_cache
 
@@ -286,15 +348,17 @@ class FluxDepletionOperator(OpenMCOperator):
         helper_kwargs : dict
             Keyword arguments for helper classes
 
-
         """
 
+        normalization_mode = helper_kwargs['normalization_mode']
         fission_yield_opts = helper_kwargs.get('fission_yield_opts', {})
 
         self._rate_helper = self._FluxDepletionRateHelper(
             self.reaction_rates.n_nuc, self.reaction_rates.n_react, self)
-
-        self._normalization_helper = SourceRateHelper()
+        if normalization_mode == "constant-power":
+            self._normalization_helper = self.FluxDepletionNormalizationHelper()
+        else:
+            self._normalization_helper = SourceRateHelper()
 
         # Select and create fission yield helper
         fission_helper = ConstantFissionYieldHelper
diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py
index 33bd182d7d..0aaa3b08ab 100644
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@@ -423,7 +423,7 @@ class ChainFissionHelper(EnergyNormalizationHelper):
 
         self._fission_q_vector = fission_qs
 
-    def update(self, fission_rates):
+    def update(self, fission_rates, mat_index=None):
         """Update energy produced with fission rates in a material
 
         Parameters
@@ -432,6 +432,8 @@ class ChainFissionHelper(EnergyNormalizationHelper):
             fission reaction rate for each isotope in the specified
             material. Should be ordered corresponding to initial
             ``rate_index`` used in :meth:`prepare`
+        mat_index : int
+            Unused
         """
         self._energy += dot(fission_rates, self._fission_q_vector)
 
diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py
index 0ea1f093c3..aafa1adf9c 100644
--- a/openmc/deplete/openmc_operator.py
+++ b/openmc/deplete/openmc_operator.py
@@ -527,7 +527,7 @@ class OpenMCOperator(TransportOperator):
 
             # Accumulate energy from fission
             if fission_ind is not None:
-                self._normalization_helper.update(tally_rates[:, fission_ind])
+                self._normalization_helper.update(tally_rates[:, fission_ind], mat_index=mat_index)
 
             # Divide by total number and store
             rates[i] = self._rate_helper.divide_by_adens(number)