diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index 3967f1a18..f112cf8cc 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -206,14 +206,15 @@ total system energy.
 The :class:`openmc.deplete.IndependentOperator` uses inner classes subclassed
 from those listed above to perform similar calculations.
 
-The following classes are used to define transfer rates to model continuous
-removal or feed of nuclides during depletion.
+The following classes are used to define external source rates or transfer rates
+to model continuous removal or feed of nuclides during depletion.
 
 .. autosummary::
    :toctree: generated
    :nosignatures:
    :template: myclass.rst
 
+   transfer_rates.ExternalSourceRates
    transfer_rates.TransferRates
 
 Intermediate Classes
diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index ab087b98e..c304b8dc6 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -29,7 +29,7 @@ from .results import Results, _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \
     _SECONDS_PER_DAY, _SECONDS_PER_JULIAN_YEAR
 from .pool import deplete
 from .reaction_rates import ReactionRates
-from .transfer_rates import TransferRates
+from .transfer_rates import TransferRates, ExternalSourceRates
 
 
 __all__ = [
@@ -607,9 +607,14 @@ class Integrator(ABC):
               next time step. Expected to be of the same shape as ``n0``
 
     transfer_rates : openmc.deplete.TransferRates
-        Instance of TransferRates class to perform continuous transfer during depletion
+        Transfer rates for the depletion system used to model continuous
+        removal/feed between materials.
 
         .. versionadded:: 0.14.0
+    external_source_rates : openmc.deplete.ExternalSourceRates
+        External source rates for the depletion system.
+
+        .. versionadded:: 0.15.3
 
     """
 
@@ -686,6 +691,7 @@ class Integrator(ABC):
         self.source_rates = np.asarray(source_rates)
 
         self.transfer_rates = None
+        self.external_source_rates = None
 
         if isinstance(solver, str):
             # Delay importing of cram module, which requires this file
@@ -731,11 +737,11 @@ class Integrator(ABC):
 
         self._solver = func
 
-    def _timed_deplete(self, n, rates, dt, matrix_func=None):
+    def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None):
         start = time.time()
         results = deplete(
-            self._solver, self.chain, n, rates, dt, matrix_func,
-            self.transfer_rates)
+            self._solver, self.chain, n, rates, dt, i, matrix_func,
+            self.transfer_rates, self.external_source_rates)
         return time.time() - start, results
 
     @abstractmethod
@@ -885,13 +891,14 @@ class Integrator(ABC):
         self.operator.finalize()
 
     def add_transfer_rate(
-            self,
-            material: Union[str, int, Material],
-            components: Sequence[str],
-            transfer_rate: float,
-            transfer_rate_units: str = '1/s',
-            destination_material: Optional[Union[str, int, Material]] = None
-        ):
+        self,
+        material: str | int | Material,
+        components: Sequence[str],
+        transfer_rate: float,
+        transfer_rate_units: str = '1/s',
+        timesteps: Sequence[int] | None = None,
+        destination_material: str | int | Material | None = None
+    ):
         """Add transfer rates to depletable material.
 
         Parameters
@@ -905,18 +912,79 @@ class Integrator(ABC):
         transfer_rate : float
             Rate at which elements are transferred. A positive or negative values
             set removal of feed rates, respectively.
-        destination_material : openmc.Material or str or int, Optional
-            Destination material to where nuclides get fed.
         transfer_rate_units : {'1/s', '1/min', '1/h', '1/d', '1/a'}
             Units for values specified in the transfer_rate argument. 's' means
             seconds, 'min' means minutes, 'h' means hours, 'a' means Julian years.
+        timesteps : list of int, optional
+            List of timestep indices where to set external source rates.
+            Defaults to None, which means the external source rate is set for
+            all timesteps.
+        destination_material : openmc.Material or str or int, Optional
+            Destination material to where nuclides get fed.
 
         """
         if self.transfer_rates is None:
-            self.transfer_rates = TransferRates(self.operator, self.operator.model)
+            if hasattr(self.operator, 'model'):
+                materials = self.operator.model.materials
+            elif hasattr(self.operator, 'materials'):
+                materials = self.operator.materials
+            self.transfer_rates = TransferRates(
+                self.operator, materials, len(self.timesteps))
+
+        if self.external_source_rates is not None and destination_material:
+            raise ValueError('Currently is not possible to set a transfer rate '
+                             'with destination matrial in combination with '
+                             'external source rates.')
+
+        self.transfer_rates.set_transfer_rate(
+            material, components, transfer_rate, transfer_rate_units,
+            timesteps, destination_material)
+
+    def add_external_source_rate(
+        self,
+        material: str | int | Material,
+        composition: dict[str, float],
+        rate: float,
+        rate_units: str = 'g/s',
+        timesteps: Sequence[int] | None = None
+    ):
+        """Add external source rates to depletable material.
+
+        Parameters
+        ----------
+        material : openmc.Material or str or int
+            Depletable material
+        composition : dict of str to float
+            External source rate composition vector, where key can be an element
+            or a nuclide and value the corresponding weight percent.
+        rate : float
+            External source rate in units of mass per time. A positive or
+            negative value corresponds to a feed or removal rate, respectively.
+        units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'}
+            Units for values specified in the `rate` argument. 's' for seconds,
+            'min' for minutes, 'h' for hours, 'a' for Julian years.
+        timesteps : list of int, optional
+            List of timestep indices where to set external source rates.
+            Defaults to None, which means the external source rate is set for
+            all timesteps.
+
+        """
+        if self.external_source_rates is None:
+            if hasattr(self.operator, 'model'):
+                materials = self.operator.model.materials
+            elif hasattr(self.operator, 'materials'):
+                materials = self.operator.materials
+            self.external_source_rates = ExternalSourceRates(
+                self.operator, materials, len(self.timesteps))
+
+        if self.transfer_rates is not None and self.transfer_rates.index_transfer:
+            raise ValueError('Currently is not possible to set an external '
+                             'source rate in combination with transfer rates '
+                             'with destination matrial.')
+
+        self.external_source_rates.set_external_source_rate(
+            material, composition, rate, rate_units, timesteps)
 
-        self.transfer_rates.set_transfer_rate(material, components, transfer_rate,
-                                      transfer_rate_units, destination_material)
 
 @add_params
 class SIIntegrator(Integrator):
@@ -1047,10 +1115,10 @@ class SIIntegrator(Integrator):
         return inherited
 
     def integrate(
-            self,
-            output: bool = True,
-            path: PathLike = "depletion_results.h5"
-        ):
+        self,
+        output: bool = True,
+        path: PathLike = "depletion_results.h5"
+    ):
         """Perform the entire depletion process across all steps
 
         Parameters
diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 4998a2c3d..9f234683a 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -703,7 +703,7 @@ class Chain:
         # Return CSC representation instead of DOK
         return matrix.tocsc()
 
-    def form_rr_term(self, tr_rates, mats):
+    def form_rr_term(self, tr_rates, current_timestep, mats):
         """Function to form the transfer rate term matrices.
 
         .. versionadded:: 0.14.0
@@ -712,6 +712,8 @@ class Chain:
         ----------
         tr_rates : openmc.deplete.TransferRates
             Instance of openmc.deplete.TransferRates
+        current_timestep : int
+            Current timestep index
         mats : string or two-tuple of strings
             Two cases are possible:
 
@@ -740,32 +742,74 @@ class Chain:
 
         for i, nuc in enumerate(self.nuclides):
             elm = re.split(r'\d+', nuc.name)[0]
-            # Build transfer terms matrices
+            # Build transfer terms (nuclide transfer only)
             if isinstance(mats, str):
                 mat = mats
-                components = tr_rates.get_components(mat)
+                components = tr_rates.get_components(mat, current_timestep)
+                if not components:
+                    break
                 if elm in components:
-                    matrix[i, i] = sum(tr_rates.get_transfer_rate(mat, elm))
+                    matrix[i, i] = sum(
+                        tr_rates.get_external_rate(mat, elm, current_timestep))
                 elif nuc.name in components:
-                    matrix[i, i] = sum(tr_rates.get_transfer_rate(mat, nuc.name))
+                    matrix[i, i] = sum(
+                        tr_rates.get_external_rate(mat, nuc.name, current_timestep))
                 else:
                     matrix[i, i] = 0.0
-            #Build transfer terms matrices
+
+            # Build transfer terms (transfer from one material into another)
             elif isinstance(mats, tuple):
                 dest_mat, mat = mats
-                if dest_mat in tr_rates.get_destination_material(mat, elm):
-                    dest_mat_idx = tr_rates.get_destination_material(mat, elm).index(dest_mat)
-                    matrix[i, i] = tr_rates.get_transfer_rate(mat, elm)[dest_mat_idx]
-                elif dest_mat in tr_rates.get_destination_material(mat, nuc.name):
-                    dest_mat_idx = tr_rates.get_destination_material(mat, nuc.name).index(dest_mat)
-                    matrix[i, i] = tr_rates.get_transfer_rate(mat, nuc.name)[dest_mat_idx]
+                components = tr_rates.get_components(mat, current_timestep, dest_mat)
+                if elm in components:
+                    matrix[i, i] = tr_rates.get_external_rate(
+                        mat, elm, current_timestep, dest_mat)[0]
+                elif nuc.name in components:
+                    matrix[i, i] = tr_rates.get_external_rate(
+                        mat, nuc.name, current_timestep, dest_mat)[0]
                 else:
                     matrix[i, i] = 0.0
-            #Nothing else is allowed
 
         # Return CSC instead of DOK
         return matrix.tocsc()
 
+    def form_ext_source_term(self, ext_source_rates, current_timestep, mat):
+        """Function to form the external source rate term vectors.
+
+        .. versionadded:: 0.15.3
+
+        Parameters
+        ----------
+        ext_source_rates : openmc.deplete.ExternalSourceRates
+            Instance of openmc.deplete.ExternalSourceRates
+        current_timestep : int
+            Current timestep index
+        mat : string
+            Material id
+
+        Returns
+        -------
+        scipy.sparse.csc_matrix
+            Sparse vector representing external source term.
+
+        """
+        if not ext_source_rates.get_components(mat, current_timestep):
+            return
+        # Use DOK as intermediate representation
+        n = len(self)
+        vector = sp.dok_matrix((n, 1))
+
+        for i, nuc in enumerate(self.nuclides):
+            # Build source term vector
+            if nuc.name in ext_source_rates.get_components(mat, current_timestep):
+                vector[i] = sum(ext_source_rates.get_external_rate(
+                    mat, nuc.name, current_timestep))
+            else:
+                vector[i] = 0.0
+
+        # Return CSC instead of DOK
+        return vector.tocsc()
+
     def get_branch_ratios(self, reaction="(n,gamma)"):
         """Return a dictionary with reaction branching ratios
 
diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py
index 50810c88a..000cb2c41 100644
--- a/openmc/deplete/integrators.py
+++ b/openmc/deplete/integrators.py
@@ -40,8 +40,8 @@ class PredictorIntegrator(Integrator):
             Time in [s] for the entire depletion interval
         source_rate : float
             Power in [W] or source rate in [neutron/sec]
-        _i : int or None
-            Iteration index. Not used
+        _i : int, optional
+            Current iteration count. Not used
 
         Returns
         -------
@@ -54,7 +54,7 @@ class PredictorIntegrator(Integrator):
             with predictor
 
         """
-        proc_time, n_end = self._timed_deplete(n, rates, dt)
+        proc_time, n_end = self._timed_deplete(n, rates, dt, _i)
         return proc_time, [n_end], []
 
 
@@ -106,12 +106,12 @@ class CECMIntegrator(Integrator):
             Eigenvalue and reaction rates from transport simulations
         """
         # deplete across first half of interval
-        time0, n_middle = self._timed_deplete(n, rates, dt / 2)
+        time0, n_middle = self._timed_deplete(n, rates, dt / 2, _i)
         res_middle = self.operator(n_middle, source_rate)
 
         # deplete across entire interval with BOS concentrations,
         # MOS reaction rates
-        time1, n_end = self._timed_deplete(n, res_middle.rates, dt)
+        time1, n_end = self._timed_deplete(n, res_middle.rates, dt, _i)
 
         return time0 + time1, [n_middle, n_end], [res_middle]
 
@@ -172,26 +172,26 @@ class CF4Integrator(Integrator):
         """
         # Step 1: deplete with matrix 1/2*A(y0)
         time1, n_eos1 = self._timed_deplete(
-            n_bos, bos_rates, dt, matrix_func=cf4_f1)
+            n_bos, bos_rates, dt, _i, matrix_func=cf4_f1)
         res1 = self.operator(n_eos1, source_rate)
 
         # Step 2: deplete with matrix 1/2*A(y1)
         time2, n_eos2 = self._timed_deplete(
-            n_bos, res1.rates, dt, matrix_func=cf4_f1)
+            n_bos, res1.rates, dt, _i, matrix_func=cf4_f1)
         res2 = self.operator(n_eos2, source_rate)
 
         # Step 3: deplete with matrix -1/2*A(y0)+A(y2)
         list_rates = list(zip(bos_rates, res2.rates))
         time3, n_eos3 = self._timed_deplete(
-            n_eos1, list_rates, dt, matrix_func=cf4_f2)
+            n_eos1, list_rates, dt, _i, matrix_func=cf4_f2)
         res3 = self.operator(n_eos3, source_rate)
 
         # Step 4: deplete with two matrix exponentials
         list_rates = list(zip(bos_rates, res1.rates, res2.rates, res3.rates))
         time4, n_inter = self._timed_deplete(
-            n_bos, list_rates, dt, matrix_func=cf4_f3)
+            n_bos, list_rates, dt, _i, matrix_func=cf4_f3)
         time5, n_eos5 = self._timed_deplete(
-            n_inter, list_rates, dt, matrix_func=cf4_f4)
+            n_inter, list_rates, dt, _i, matrix_func=cf4_f4)
 
         return (time1 + time2 + time3 + time4 + time5,
                 [n_eos1, n_eos2, n_eos3, n_eos5],
@@ -249,17 +249,17 @@ class CELIIntegrator(Integrator):
             simulation
         """
         # deplete to end using BOS rates
-        proc_time, n_ce = self._timed_deplete(n_bos, rates, dt)
+        proc_time, n_ce = self._timed_deplete(n_bos, rates, dt, _i)
         res_ce = self.operator(n_ce, source_rate)
 
         # deplete using two matrix exponentials
         list_rates = list(zip(rates, res_ce.rates))
 
         time_le1, n_inter = self._timed_deplete(
-            n_bos, list_rates, dt, matrix_func=celi_f1)
+            n_bos, list_rates, dt, _i, matrix_func=celi_f1)
 
         time_le2, n_end = self._timed_deplete(
-            n_inter, list_rates, dt, matrix_func=celi_f2)
+            n_inter, list_rates, dt, _i, matrix_func=celi_f2)
 
         return proc_time + time_le1 + time_le1, [n_ce, n_end], [res_ce]
 
@@ -316,20 +316,20 @@ class EPCRK4Integrator(Integrator):
         """
 
         # Step 1: deplete with matrix A(y0) / 2
-        time1, n1 = self._timed_deplete(n, rates, dt, matrix_func=rk4_f1)
+        time1, n1 = self._timed_deplete(n, rates, dt, _i, matrix_func=rk4_f1)
         res1 = self.operator(n1, source_rate)
 
         # Step 2: deplete with matrix A(y1) / 2
-        time2, n2 = self._timed_deplete(n, res1.rates, dt, matrix_func=rk4_f1)
+        time2, n2 = self._timed_deplete(n, res1.rates, dt, _i, matrix_func=rk4_f1)
         res2 = self.operator(n2, source_rate)
 
         # Step 3: deplete with matrix A(y2)
-        time3, n3 = self._timed_deplete(n, res2.rates, dt)
+        time3, n3 = self._timed_deplete(n, res2.rates, dt, _i)
         res3 = self.operator(n3, source_rate)
 
         # Step 4: deplete with matrix built from weighted rates
         list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates))
-        time4, n4 = self._timed_deplete(n, list_rates, dt, matrix_func=rk4_f4)
+        time4, n4 = self._timed_deplete(n, list_rates, dt, _i, matrix_func=rk4_f4)
 
         return (time1 + time2 + time3 + time4, [n1, n2, n3, n4], [res1, res2, res3])
 
@@ -414,9 +414,9 @@ class LEQIIntegrator(Integrator):
             self._prev_rates, bos_res.rates, repeat(prev_dt), repeat(dt)))
 
         time1, n_inter = self._timed_deplete(
-            n_bos, le_inputs, dt, matrix_func=leqi_f1)
+            n_bos, le_inputs, dt, i, matrix_func=leqi_f1)
         time2, n_eos0 = self._timed_deplete(
-            n_inter, le_inputs, dt, matrix_func=leqi_f2)
+            n_inter, le_inputs, dt, i, matrix_func=leqi_f2)
 
         res_inter = self.operator(n_eos0, source_rate)
 
@@ -425,9 +425,9 @@ class LEQIIntegrator(Integrator):
             repeat(prev_dt), repeat(dt)))
 
         time3, n_inter = self._timed_deplete(
-            n_bos, qi_inputs, dt, matrix_func=leqi_f3)
+            n_bos, qi_inputs, dt, i, matrix_func=leqi_f3)
         time4, n_eos1 = self._timed_deplete(
-            n_inter, qi_inputs, dt, matrix_func=leqi_f4)
+            n_inter, qi_inputs, dt, i, matrix_func=leqi_f4)
 
         # store updated rates
         self._prev_rates = copy.deepcopy(bos_res.rates)
@@ -478,7 +478,7 @@ class SICELIIntegrator(SIIntegrator):
             Eigenvalue and reaction rates from intermediate transport
             simulations
         """
-        proc_time, n_eos = self._timed_deplete(n_bos, bos_rates, dt)
+        proc_time, n_eos = self._timed_deplete(n_bos, bos_rates, dt, _i)
         n_inter = copy.deepcopy(n_eos)
 
         # Begin iteration
@@ -494,9 +494,9 @@ class SICELIIntegrator(SIIntegrator):
 
             list_rates = list(zip(bos_rates, res_bar.rates))
             time1, n_inter = self._timed_deplete(
-                n_bos, list_rates, dt, matrix_func=celi_f1)
+                n_bos, list_rates, dt, _i, matrix_func=celi_f1)
             time2, n_inter = self._timed_deplete(
-                n_inter, list_rates, dt, matrix_func=celi_f2)
+                n_inter, list_rates, dt, _i, matrix_func=celi_f2)
             proc_time += time1 + time2
 
         # end iteration
@@ -560,9 +560,9 @@ class SILEQIIntegrator(SIIntegrator):
         inputs = list(zip(self._prev_rates, bos_rates,
                           repeat(prev_dt), repeat(dt)))
         proc_time, n_inter = self._timed_deplete(
-            n_bos, inputs, dt, matrix_func=leqi_f1)
+            n_bos, inputs, dt, i, matrix_func=leqi_f1)
         time1, n_eos = self._timed_deplete(
-            n_inter, inputs, dt, matrix_func=leqi_f2)
+            n_inter, inputs, dt, i, matrix_func=leqi_f2)
 
         proc_time += time1
         n_inter = copy.deepcopy(n_eos)
@@ -580,9 +580,9 @@ class SILEQIIntegrator(SIIntegrator):
             inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates,
                               repeat(prev_dt), repeat(dt)))
             time1, n_inter = self._timed_deplete(
-                n_bos, inputs, dt, matrix_func=leqi_f3)
+                n_bos, inputs, dt, i, matrix_func=leqi_f3)
             time2, n_inter = self._timed_deplete(
-                n_inter, inputs, dt, matrix_func=leqi_f4)
+                n_inter, inputs, dt, i, matrix_func=leqi_f4)
             proc_time += time1 + time2
 
         return proc_time, [n_eos, n_inter], [res_bar]
diff --git a/openmc/deplete/pool.py b/openmc/deplete/pool.py
index 27ecaa4dd..03b050af3 100644
--- a/openmc/deplete/pool.py
+++ b/openmc/deplete/pool.py
@@ -5,7 +5,7 @@ Provided to avoid some circular imports
 from itertools import repeat, starmap
 from multiprocessing import Pool
 
-from scipy.sparse import bmat
+from scipy.sparse import bmat, hstack, vstack, csc_matrix
 import numpy as np
 
 from openmc.mpi import comm
@@ -40,8 +40,8 @@ def _distribute(items):
             return items[j:j + chunk_size]
         j += chunk_size
 
-def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
-            *matrix_args):
+def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
+            transfer_rates=None, external_source_rates=None, *matrix_args):
     """Deplete materials using given reaction rates for a specified time
 
     Parameters
@@ -58,15 +58,21 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
         Reaction rates (from transport operator)
     dt : float
         Time in [s] to deplete for
+    current_timestep : int
+        Current timestep index
     maxtrix_func : callable, optional
         Function to form the depletion matrix after calling ``matrix_func(chain,
         rates, fission_yields)``, where ``fission_yields = {parent: {product:
         yield_frac}}`` Expected to return the depletion matrix required by
         ``func``
     transfer_rates : openmc.deplete.TransferRates, Optional
-        Object to perform continuous reprocessing.
+        Transfer rates for continuous removal/feed.
 
         .. versionadded:: 0.14.0
+    external_source_rates : openmc.deplete.ExternalSourceRates, Optional
+        External source rates for continuous removal/feed.
+
+        .. versionadded:: 0.15.3
     matrix_args: Any, optional
         Additional arguments passed to matrix_func
 
@@ -93,15 +99,17 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
         matrices = map(matrix_func, repeat(chain), rates, fission_yields,
                        *matrix_args)
 
-    if transfer_rates is not None:
+    if (transfer_rates is not None and
+        current_timestep in transfer_rates.external_timesteps):
         # Calculate transfer rate terms as diagonal matrices
         transfers = map(chain.form_rr_term, repeat(transfer_rates),
-                        transfer_rates.local_mats)
+                        repeat(current_timestep), transfer_rates.local_mats)
+
         # Subtract transfer rate terms from Bateman matrices
         matrices = [matrix - transfer for (matrix, transfer) in zip(matrices,
                                                                     transfers)]
 
-        if len(transfer_rates.index_transfer) > 0:
+        if current_timestep in transfer_rates.index_transfer:
             # Gather all on comm.rank 0
             matrices = comm.gather(matrices)
             n = comm.gather(n)
@@ -112,10 +120,12 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
                 n = [n_elm for n_mat in n for  n_elm in n_mat]
 
                 # Calculate transfer rate terms as diagonal matrices
-                transfer_pair = {
-                    mat_pair: chain.form_rr_term(transfer_rates, mat_pair)
-                    for mat_pair in transfer_rates.index_transfer
-                }
+                transfer_pair = {}
+                for mat_pair in transfer_rates.index_transfer[current_timestep]:
+                    transfer_matrix = chain.form_rr_term(transfer_rates,
+                                                         current_timestep,
+                                                         mat_pair)
+                    transfer_pair[mat_pair] = transfer_matrix
 
                 # Combine all matrices together in a single matrix of matrices
                 # to be solved in one go
@@ -129,7 +139,7 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
                         if row == col:
                             # Fill the diagonals with the Bateman matrices
                             cols.append(matrices[row])
-                        elif mat_pair in transfer_rates.index_transfer:
+                        elif mat_pair in transfer_rates.index_transfer[current_timestep]:
                             # Fill the off-diagonals with the transfer pair matrices
                             cols.append(transfer_pair[mat_pair])
                         else:
@@ -155,6 +165,25 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
 
             return n_result
 
+    if (external_source_rates is not None and
+        current_timestep in external_source_rates.external_timesteps):
+        # Calculate external source term vectors
+        sources = map(chain.form_ext_source_term, repeat(external_source_rates),
+                      repeat(current_timestep), external_source_rates.local_mats)
+
+        # stack vector column at the end of the matrix
+        matrices = [
+            hstack([matrix, source])
+            for matrix, source in zip(matrices, sources)
+        ]
+
+        # Add a last row of zeroes to the matrices and append 1 to the last row
+        # of the nuclide vectors
+        for i, matrix in enumerate(matrices):
+            if not np.equal(*matrix.shape):
+                matrices[i] = vstack([matrix, csc_matrix([0]*matrix.shape[1])])
+                n[i] = np.append(n[i], 1.0)
+
     inputs = zip(matrices, n, repeat(dt))
 
     if USE_MULTIPROCESSING:
@@ -163,4 +192,10 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
     else:
         n_result = list(starmap(func, inputs))
 
+    # Remove extra value at the end of the nuclide vectors
+    if (external_source_rates is not None and
+        current_timestep in external_source_rates.external_timesteps):
+        external_source_rates.reformat_nuclide_vectors(n)
+        external_source_rates.reformat_nuclide_vectors(n_result)
+
     return n_result
diff --git a/openmc/deplete/transfer_rates.py b/openmc/deplete/transfer_rates.py
index 01c9b2e35..4c28c7d15 100644
--- a/openmc/deplete/transfer_rates.py
+++ b/openmc/deplete/transfer_rates.py
@@ -1,31 +1,31 @@
+from collections import defaultdict
 from numbers import Real
 import re
+from typing import Sequence
+
+import numpy as np
 
 from openmc.checkvalue import check_type, check_value
 from openmc import Material
-from openmc.data import ELEMENT_SYMBOL
+from openmc.data import ELEMENT_SYMBOL, isotopes, AVOGADRO, atomic_mass
+from .results import _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \
+    _SECONDS_PER_DAY, _SECONDS_PER_JULIAN_YEAR
 
 
-class TransferRates:
-    """Class for defining continuous removals and feeds.
 
-    Molten Salt Reactors (MSRs) benefit from continuous reprocessing,
-    which removes fission products and feeds fresh fuel into the system. MSRs
-    inspired the development of this class.
+class ExternalRates:
+    """External rates class for defining addition terms of depletion equation.
 
-    An instance of this class can be passed directly to an instance of one of
-    the :class:`openmc.deplete.Integrator` classes.
-
-    .. versionadded:: 0.14.0
+    .. versionadded:: 0.15.3
 
     Parameters
     ----------
     operator : openmc.TransportOperator
         Depletion operator
-    model : openmc.Model
-        OpenMC model containing materials and geometry. If using
-        :class:`openmc.deplete.CoupledOperator`, the model must also contain
-        a :class:`opnemc.Settings` object.
+    materials : openmc.Materials
+        OpenMC materials.
+    number_of_timesteps : int
+        Total number of depletion timesteps
 
     Attributes
     ----------
@@ -33,22 +33,25 @@ class TransferRates:
         All burnable material IDs.
     local_mats : list of str
         All burnable material IDs being managed by a single process
-    transfer_rates : dict of str to dict
-        Container of transfer rates, components (elements and/or nuclides) and
-        destination material
-    index_transfer : Set of pair of str
-        Pair of strings needed to build final matrix (destination_material, mat)
+    number_of_timesteps : int
+        Total number of depletion timesteps
+    external_rates : dict of str to dict
+        Container of timesteps, external rates, components (elements and/or
+        nuclides) and optionally destination material
+    external_timesteps : list of int
+        Container of all timesteps indeces with an external rate defined.
     """
 
-    def __init__(self, operator, model):
+    def __init__(self, operator, materials, number_of_timesteps):
 
-        self.materials = model.materials
+        self.materials = materials
         self.burnable_mats = operator.burnable_mats
         self.local_mats = operator.local_mats
+        self.number_of_timesteps = number_of_timesteps
 
-        #initialize transfer rates container dict
-        self.transfer_rates = {mat: {} for mat in self.burnable_mats}
-        self.index_transfer = set()
+        # initialize transfer rates container dict
+        self.external_rates = {mat: defaultdict(list) for mat in self.burnable_mats}
+        self.external_timesteps = []
 
     def _get_material_id(self, val):
         """Helper method for getting material id from Material obj or name.
@@ -71,7 +74,7 @@ class TransferRates:
                 check_value('Material ID', str(val), self.burnable_mats)
             else:
                 check_value('Material name', val,
-                        [mat.name for mat in self.materials if mat.depletable])
+                            [mat.name for mat in self.materials if mat.depletable])
                 val = [mat.id for mat in self.materials if mat.name == val][0]
 
         elif isinstance(val, int):
@@ -79,7 +82,13 @@ class TransferRates:
 
         return str(val)
 
-    def get_transfer_rate(self, material, component):
+    def get_external_rate(
+        self,
+        material: str | int | Material,
+        component: str,
+        timestep: int,
+        destination_material: str | int | Material | None = None
+    ):
         """Return transfer rate for given material and element.
 
         Parameters
@@ -88,62 +97,114 @@ class TransferRates:
             Depletable material
         component : str
             Element or nuclide to get transfer rate value
+        timestep : int
+            Current timestep index
+        destination_material : openmc.Material or str or int, Optional
+            Destination material to where nuclides get fed
 
         Returns
         -------
-        transfer_rate : list of floats
-            Transfer rate values
+        external_rate : list of floats
+            External rate values
 
         """
         material_id = self._get_material_id(material)
         check_type('component', component, str)
-        return [i[0] for i in self.transfer_rates[material_id][component]]
-
-    def get_destination_material(self, material, component):
-        """Return destination material for given material and
-        component, if defined.
-
-        Parameters
-        ----------
-        material : openmc.Material or str or int
-            Depletable material
-        component : str
-            Element or nuclide that gets transferred to another material.
-
-        Returns
-        -------
-        destination_material_id : list of str
-            Depletable material ID to where the element or nuclide gets
-            transferred
-
-        """
-        material_id = self._get_material_id(material)
-        check_type('component', component, str)
-        if component in self.transfer_rates[material_id]:
-            return [i[1] for i in self.transfer_rates[material_id][component]]
+        if destination_material is not None:
+            dest_mat_id = self._get_material_id(destination_material)
+            return [i[1] for i in self.external_rates[material_id][component]
+                    if timestep in i[0] and dest_mat_id == i[2]]
         else:
-            return []
+            return [i[1] for i in self.external_rates[material_id][component]
+                    if timestep in i[0]]
 
-    def get_components(self, material):
-        """Extract removing elements and/or nuclides for a given material
+    def get_components(self, material, timestep, destination_material=None):
+        """Extract removing elements and/or nuclides for a given material at a
+        given timestep
 
         Parameters
         ----------
         material : openmc.Material or str or int
             Depletable material
+        timestep : int
+            Current timestep index
+        destination_material : openmc.Material or str or int, Optional
+            Destination material to where nuclides get fed
 
         Returns
         -------
-        elements : list
-            List of elements and nuclides where transfer rates exist
+        components : list
+            List of elements or nuclides with external rates set at a given
+            timestep
 
         """
         material_id = self._get_material_id(material)
-        if material_id in self.transfer_rates:
-            return self.transfer_rates[material_id].keys()
+        if destination_material is not None:
+            dest_mat_id = self._get_material_id(destination_material)
+        else:
+            dest_mat_id = None
+
+        all_components = []
+        if material_id in self.external_rates:
+            mat_components = self.external_rates[material_id]
+
+            for component in mat_components:
+                if dest_mat_id:
+                    # check for both timestep and destination material ids
+                    if np.isin(timestep, [val[0] for val in mat_components[component]]) and \
+                       np.isin(dest_mat_id, [val[2] for val in mat_components[component]]):
+                        all_components.append(component)
+                else:
+                    # check only for timesteps
+                    if np.isin(timestep, [val[0] for val in mat_components[component]]):
+                        all_components.append(component)
+        return all_components
+
+
+class TransferRates(ExternalRates):
+    """Class for defining continuous removals and feeds.
+
+    Molten Salt Reactors (MSRs) benefit from continuous reprocessing,
+    which removes fission products and feeds fresh fuel into the system. MSRs
+    inspired the development of this class.
+
+    An instance of this class can be passed directly to an instance of one of
+    the :class:`openmc.deplete.Integrator` classes.
+
+    .. versionadded:: 0.14.0
+
+    Parameters
+    ----------
+    operator : openmc.TransportOperator
+        Depletion operator
+    materials : openmc.Materials
+        OpenMC materials.
+    number_of_timesteps : int
+        Total number of depletion timesteps
+
+    Attributes
+    ----------
+    burnable_mats : list of str
+        All burnable material IDs.
+    local_mats : list of str
+        All burnable material IDs being managed by a single process
+    external_rates : dict of str to dict
+        Container of timesteps, transfer rates, components (elements and/or
+        nuclides) and destination material
+    external_timesteps : list of int
+        Container of all timesteps indeces with an external rate defined.
+    index_transfer : Set of pair of str
+        Pair of strings needed to build final matrix (destination_material, mat)
+    """
+
+    def __init__(self, operator, materials, number_of_timesteps):
+        super().__init__(operator, materials, number_of_timesteps)
+        self.index_transfer = defaultdict(list)
+        self.chain_nuclides = [nuc.name for nuc in operator.chain.nuclides]
 
     def set_transfer_rate(self, material, components, transfer_rate,
-                          transfer_rate_units='1/s', destination_material=None):
+                          transfer_rate_units='1/s', timesteps=None,
+                          destination_material=None):
         """Set element and/or nuclide transfer rates in a depletable material.
 
         Parameters
@@ -157,11 +218,14 @@ class TransferRates:
         transfer_rate : float
             Rate at which elements and/or nuclides are transferred. A positive or
             negative value corresponds to a removal or feed rate, respectively.
-        destination_material : openmc.Material or str or int, Optional
-            Destination material to where nuclides get fed.
         transfer_rate_units : {'1/s', '1/min', '1/h', '1/d', '1/a'}
             Units for values specified in the transfer_rate argument. 's' for
             seconds, 'min' for minutes, 'h' for hours, 'a' for Julian years.
+        timesteps : list of int, Optional
+            List of timestep indeces where to set transfer rates.
+            Default to None means the transfer rate is set for all timesteps.
+        destination_material : openmc.Material or str or int, Optional
+            Destination material to where nuclides get fed.
 
         """
         material_id = self._get_material_id(material)
@@ -183,19 +247,25 @@ class TransferRates:
         if transfer_rate_units in ('1/s', '1/sec'):
             unit_conv = 1
         elif transfer_rate_units in ('1/min', '1/minute'):
-            unit_conv = 60
+            unit_conv = _SECONDS_PER_MINUTE
         elif transfer_rate_units in ('1/h', '1/hr', '1/hour'):
-            unit_conv = 60*60
+            unit_conv = _SECONDS_PER_HOUR
         elif transfer_rate_units in ('1/d', '1/day'):
-            unit_conv = 24*60*60
+            unit_conv = _SECONDS_PER_DAY
         elif transfer_rate_units in ('1/a', '1/year'):
-            unit_conv = 365.25*24*60*60
+            unit_conv = _SECONDS_PER_JULIAN_YEAR
         else:
-            raise ValueError('Invalid transfer rate unit '
-                             f'"{transfer_rate_units}"')
+            raise ValueError(f'Invalid transfer rate unit "{transfer_rate_units}"')
+
+        if timesteps is not None:
+            for timestep in timesteps:
+                check_value('timestep', timestep, range(self.number_of_timesteps))
+            timesteps = np.array(timesteps)
+        else:
+            timesteps = np.arange(self.number_of_timesteps)
 
         for component in components:
-            current_components = self.transfer_rates[material_id].keys()
+            current_components = self.external_rates[material_id].keys()
             split_component = re.split(r'\d+', component)
             element = split_component[0]
             if element not in ELEMENT_SYMBOL.values():
@@ -219,11 +289,144 @@ class TransferRates:
                                          f'where element {element} already has '
                                          'a transfer rate.')
 
-            if component in self.transfer_rates[material_id]:
-                self.transfer_rates[material_id][component].append(
-                    (transfer_rate / unit_conv, destination_material_id))
-            else:
-                self.transfer_rates[material_id][component] = [
-                    (transfer_rate / unit_conv, destination_material_id)]
+            self.external_rates[material_id][component].append(
+                (timesteps, transfer_rate/unit_conv, destination_material_id))
+
             if destination_material_id is not None:
-                self.index_transfer.add((destination_material_id, material_id))
+                for timestep in timesteps:
+                    self.index_transfer[timestep].append(
+                        (destination_material_id, material_id))
+
+            self.external_timesteps = np.unique(np.concatenate(
+                    [self.external_timesteps, timesteps]))
+
+
+class ExternalSourceRates(ExternalRates):
+    """Class for defining external source rates.
+
+    An instance of this class can be passed directly to an instance of one of
+    the :class:`openmc.deplete.Integrator` classes.
+
+    .. versionadded:: 0.15.3
+
+    Parameters
+    ----------
+    operator : openmc.TransportOperator
+        Depletion operator
+    materials : openmc.Materials
+        OpenMC materials.
+    number_of_timesteps : int
+        Total number of depletion timesteps
+
+    Attributes
+    ----------
+    burnable_mats : list of str
+        All burnable material IDs.
+    local_mats : list of str
+        All burnable material IDs being managed by a single process
+    external_timesteps : list of int
+        Container of all timesteps indeces with an external rate defined.
+    external_rates : dict of str to dict
+        Container of timesteps external source rates, and components
+        (elements and/or nuclides)
+    """
+
+    def reformat_nuclide_vectors(self, vectors):
+        """Remove last element of nuclide vector that was added for handling
+        external source rates by the depletion solver.
+
+        Parameters
+        ----------
+        vectors : list of array
+            List of nuclides vector to reformat
+
+        """
+        for mat_index, i in enumerate(self.local_mats):
+            if self.external_rates[i]:
+                vectors[mat_index] = vectors[mat_index][:-1]
+
+    def set_external_source_rate(
+        self,
+        material: str | int | Material,
+        composition: dict[str, float],
+        rate: float,
+        rate_units: str = 'g/s',
+        timesteps: Sequence[int] | None = None
+    ):
+        """Set element and/or nuclide composition vector external source rates
+        to a depletable material.
+
+        Parameters
+        ----------
+        material : openmc.Material or str or int
+            Depletable material
+        composition : dict of str to float
+            External source rate composition vector, where key can be an element
+            or a nuclide and value the corresponding weight percent.
+        rate : float
+            External source rate in units of mass per time. A positive or
+            negative value corresponds to a feed or removal rate, respectively.
+        units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'}
+            Units for values specified in the `rate` argument. 's' for seconds,
+            'min' for minutes, 'h' for hours, 'a' for Julian years.
+        timesteps : list of int, optional
+            List of timestep indices where to set external source rates. Default
+            to None, which means the external source rate is set for all
+            timesteps.
+
+        """
+
+        material_id = self._get_material_id(material)
+        check_type('rate', rate, Real)
+        check_type('composition', composition, dict, str)
+
+        if rate_units in ('g/s', 'g/sec'):
+            unit_conv = 1
+        elif rate_units in ('g/min', 'g/minute'):
+            unit_conv = _SECONDS_PER_MINUTE
+        elif rate_units in ('g/h', 'g/hr', 'g/hour'):
+            unit_conv = _SECONDS_PER_HOUR
+        elif rate_units in ('g/d', 'g/day'):
+            unit_conv = _SECONDS_PER_DAY
+        elif rate_units in ('g/a', 'g/year'):
+            unit_conv = _SECONDS_PER_JULIAN_YEAR
+        else:
+            raise ValueError(f'Invalid external source rate unit "{rate_units}"')
+
+        if timesteps is not None:
+            for timestep in timesteps:
+                check_value('timestep', timestep, range(self.number_of_timesteps))
+            timesteps = np.asarray(timesteps)
+        else:
+            timesteps = np.arange(self.number_of_timesteps)
+
+        components = composition.keys()
+        percents = composition.values()
+        norm_percents = [float(i) / sum(percents) for i in percents]
+
+        atoms_per_nuc = {}
+        for component, percent in zip(components, norm_percents):
+            split_component = re.split(r'\d+', component)
+            element = split_component[0]
+            if element not in ELEMENT_SYMBOL.values():
+                raise ValueError(f'{component} is not a valid nuclide or element.')
+
+            if len(split_component) == 1:
+                if not isotopes(component):
+                    raise ValueError(f'Cannot add element {component} '
+                                     'as it is not naturally abundant. '
+                                     'Specify a nuclide vector instead.')
+                for nuc, frac in isotopes(component):
+                    atoms_per_nuc[nuc] = (rate / atomic_mass(nuc) * AVOGADRO *
+                                          frac * percent / unit_conv)
+
+            else:
+                atoms_per_nuc[component] = (rate / atomic_mass(component) *
+                                            AVOGADRO * percent / unit_conv)
+
+        for nuc, val in atoms_per_nuc.items():
+            self.external_rates[material_id][nuc].append((timesteps, val, None))
+
+        self.external_timesteps = np.unique(np.concatenate(
+            [self.external_timesteps, timesteps]
+        ))
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_ext_source.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_ext_source.h5
new file mode 100644
index 000000000..aa09e1bcf
Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_ext_source.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_feed.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_feed.h5
index b9be76345..284d5cbe4 100644
Binary files a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_feed.h5 and b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_feed.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_removal.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_removal.h5
index 4b24bed52..271af2410 100644
Binary files a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_removal.h5 and b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_removal.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_transfer.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_transfer.h5
index 51173f778..25365f6f3 100644
Binary files a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_transfer.h5 and b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_transfer.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_feed.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_feed.h5
index d53adc7d8..d080b4470 100644
Binary files a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_feed.h5 and b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_feed.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_removal.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_removal.h5
index 921755e15..36f89e009 100644
Binary files a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_removal.h5 and b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_removal.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_ext_source.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_ext_source.h5
new file mode 100644
index 000000000..f3b1b171e
Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_ext_source.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_transfer.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_transfer.h5
index 469149f67..b0ba99740 100644
Binary files a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_transfer.h5 and b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_transfer.h5 differ
diff --git a/tests/regression_tests/deplete_with_transfer_rates/test.py b/tests/regression_tests/deplete_with_transfer_rates/test.py
index 10d60866a..4ad009d22 100644
--- a/tests/regression_tests/deplete_with_transfer_rates/test.py
+++ b/tests/regression_tests/deplete_with_transfer_rates/test.py
@@ -1,8 +1,7 @@
-""" TransferRates depletion test suite """
+""" ExternalRates depletion test suite """
 
 from pathlib import Path
 import shutil
-import sys
 
 import numpy as np
 import pytest
@@ -11,7 +10,7 @@ import openmc.deplete
 from openmc.deplete import CoupledOperator
 
 from tests.regression_tests import config, assert_reaction_rates_equal, \
-    assert_atoms_equal, assert_same_mats
+    assert_atoms_equal
 
 
 @pytest.fixture
@@ -44,9 +43,11 @@ def model():
     settings.particles = 100
     settings.inactive = 0
     settings.batches = 10
+    settings.seed = 1
 
     return openmc.Model(geometry, materials, settings)
 
+
 @pytest.mark.parametrize("rate, dest_mat, power, ref_result", [
     (1e-5, None, 0.0, 'no_depletion_only_removal'),
     (-1e-5, None, 0.0, 'no_depletion_only_feed'),
@@ -83,6 +84,40 @@ def test_transfer_rates(run_in_tmpdir, model, rate, dest_mat, power, ref_result)
     res_ref = openmc.deplete.Results(path_reference)
     res_test = openmc.deplete.Results(path_test)
 
-    assert_same_mats(res_ref, res_test)
-    assert_atoms_equal(res_ref, res_test, 1e-6)
+    assert_atoms_equal(res_ref, res_test)
+    assert_reaction_rates_equal(res_ref, res_test)
+
+
+@pytest.mark.parametrize("rate, power, ref_result", [
+    (1e-1, 0.0, 'no_depletion_with_ext_source'),
+    (1e-1, 174., 'depletion_with_ext_source'),
+])
+def test_external_source_rates(run_in_tmpdir, model, rate, power, ref_result):
+    """Tests external_rates depletion class with external source rates"""
+
+    chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
+
+    external_source_vector = {'U': 1}
+
+    op = CoupledOperator(model, chain_file)
+    op.round_number = True
+    integrator = openmc.deplete.PredictorIntegrator(
+        op, [1], power, timestep_units='d')
+    integrator.add_external_source_rate('f', external_source_vector, rate)
+    integrator.integrate()
+
+    # Get path to test and reference results
+    path_test = op.output_dir / 'depletion_results.h5'
+    path_reference = Path(__file__).with_name(f'ref_{ref_result}.h5')
+
+    # If updating results, do so and return
+    if config['update']:
+        shutil.copyfile(str(path_test), str(path_reference))
+        return
+
+    # Load the reference/test results
+    res_ref = openmc.deplete.Results(path_reference)
+    res_test = openmc.deplete.Results(path_test)
+
+    assert_atoms_equal(res_ref, res_test)
     assert_reaction_rates_equal(res_ref, res_test)
diff --git a/tests/unit_tests/test_deplete_external_source_rates.py b/tests/unit_tests/test_deplete_external_source_rates.py
new file mode 100644
index 000000000..a8cf9dde5
--- /dev/null
+++ b/tests/unit_tests/test_deplete_external_source_rates.py
@@ -0,0 +1,170 @@
+""" Tests for ExternalSourceRates class """
+
+from pathlib import Path
+
+import pytest
+import numpy as np
+import re
+
+import openmc
+from openmc.data import AVOGADRO, atomic_mass
+from openmc.deplete import CoupledOperator
+from openmc.deplete.transfer_rates import ExternalSourceRates
+from openmc.deplete.abc import (_SECONDS_PER_MINUTE, _SECONDS_PER_HOUR,
+                                _SECONDS_PER_DAY, _SECONDS_PER_JULIAN_YEAR)
+
+CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
+
+@pytest.fixture
+def model():
+    openmc.reset_auto_ids()
+    f = openmc.Material(name="f")
+    f.add_element("U", 1, enrichment=4.25)
+    f.add_element("O", 2)
+    f.set_density("g/cm3", 10.4)
+
+    w = openmc.Material(name="w")
+    w.add_element("O", 1)
+    w.add_element("H", 2)
+    w.set_density("g/cm3", 1.0)
+    w.depletable = True
+
+    # material just to test multiple destination material
+    h = openmc.Material(name="h")
+    h.add_element("He", 1)
+    h.set_density("g/cm3", 1.78e-4)
+
+    radii = [0.42, 0.45]
+    f.volume = np.pi * radii[0] ** 2
+    w.volume = np.pi * (radii[1]**2 - radii[0]**2)
+    h.volume = 1
+    materials = openmc.Materials([f, w, h])
+
+    surf_f = openmc.Sphere(r=radii[0])
+    surf_w = openmc.Sphere(r=radii[1], boundary_type='vacuum')
+    surf_h = openmc.Sphere(x0=10, r=1, boundary_type='vacuum')
+    cell_f = openmc.Cell(fill=f, region=-surf_f)
+    cell_w = openmc.Cell(fill=w, region=+surf_f & -surf_w)
+    cell_h = openmc.Cell(fill=h, region=-surf_h)
+    geometry = openmc.Geometry([cell_f, cell_w, cell_h])
+
+    settings = openmc.Settings()
+    settings.particles = 1000
+    settings.inactive = 10
+    settings.batches = 50
+
+    return openmc.Model(geometry, materials, settings)
+
+
+@pytest.mark.parametrize(
+"case_name, external_source_vectors, external_source_rate, timesteps", [
+    ('elements', [{'U': 0.9, 'Xe': 0.1}], 1, None),
+    ('nuclides', [{'I135': 0.1, 'Gd156': 0.3, 'Gd157': 0.6}], 1, None),
+    ('nuclides_elements', [{'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01, 'U': 0.8,
+                           'Xe': 0.08}], 1, None),
+    ('elements_nuclides', [{'U': 0.78, 'Xe': 0.1, 'I135': 0.01, 'Gd156': 0.1,
+                           'Gd157': 0.01}], 1, None),
+    ('multiple_vectors', [{'U': 1.}, {'Xe': 1}], 1, None),
+    ('timesteps', [{'U': 0.9, 'Xe': 0.1}], 1, [1]),
+    ('rates_invalid_1', [{'Gb': 1.}], 1, None),
+    ('rates_invalid_2', [{'Pu': 1.}], 1, None)
+    ])
+def test_get_set(model, case_name, external_source_vectors, external_source_rate,
+                 timesteps):
+    """Tests the get/set methods"""
+
+    op = CoupledOperator(model, CHAIN_PATH)
+    number_of_timesteps = 2
+    transfer = ExternalSourceRates(op, model.materials, number_of_timesteps)
+
+    if timesteps is None:
+        timesteps = np.arange(number_of_timesteps)
+
+    # Test by Openmc material, material name and material id
+    material= [m for m in model.materials if m.depletable][0]
+
+    for material_input in [material, material.name, material.id]:
+        for external_source_vector in external_source_vectors:
+            if case_name == 'rates_invalid_1':
+                with pytest.raises(ValueError, match='Gb is not a valid '
+                                   'nuclide or element.'):
+                    transfer.set_external_source_rate(material_input,
+                                               external_source_vector,
+                                               external_source_rate)
+            elif case_name == 'rates_invalid_2':
+                with pytest.raises(ValueError, match='Cannot add element Pu'):
+                    transfer.set_external_source_rate(material_input,
+                                               external_source_vector,
+                                               external_source_rate)
+            else:
+                transfer.set_external_source_rate(material_input,
+                                               external_source_vector,
+                                               external_source_rate,
+                                               timesteps=timesteps)
+                for component, percent in external_source_vector.items():
+                    split_component = re.split(r'\d+', component)
+                    if len(split_component) == 1:
+                        for nuc, frac in openmc.data.isotopes(component):
+                            val = external_source_rate * percent * frac * \
+                                AVOGADRO / atomic_mass(nuc)
+                            assert transfer.get_external_rate(
+                                material_input, nuc, timesteps)[0] == pytest.approx(val)
+                    else:
+                        val = external_source_rate * percent * AVOGADRO / atomic_mass(component)
+                        assert transfer.get_external_rate(
+                            material_input, component, timesteps)[0] == pytest.approx(val)
+
+                assert np.all(transfer.external_timesteps == timesteps)
+
+
+@pytest.mark.parametrize("units, unit_conv", [
+    ('g/s', 1),
+    ('g/sec', 1),
+    ('g/min', _SECONDS_PER_MINUTE),
+    ('g/minute', _SECONDS_PER_MINUTE),
+    ('g/h', _SECONDS_PER_HOUR),
+    ('g/hr', _SECONDS_PER_HOUR),
+    ('g/hour', _SECONDS_PER_HOUR),
+    ('g/d', _SECONDS_PER_DAY),
+    ('g/day', _SECONDS_PER_DAY),
+    ('g/a', _SECONDS_PER_JULIAN_YEAR),
+    ('g/year', _SECONDS_PER_JULIAN_YEAR),
+    ])
+def test_units(units, unit_conv, model):
+    """ Units testing"""
+    # create external rate Xe
+    components = ['Xe135', 'U235']
+    external_source_rate = 1.0
+    number_of_timesteps = 2
+    op = CoupledOperator(model, CHAIN_PATH)
+    transfer = ExternalSourceRates(op, model.materials, number_of_timesteps)
+    timesteps = np.arange(number_of_timesteps)
+
+    for component in components:
+        rate = external_source_rate * unit_conv * atomic_mass(component) / AVOGADRO
+        transfer.set_external_source_rate('f', {component: 1}, rate, rate_units=units)
+        assert transfer.get_external_rate(
+            'f', component, timesteps)[0] == pytest.approx(external_source_rate)
+
+
+def test_external_source(run_in_tmpdir, model):
+    """Tests external source depletion class without neither reaction rates nor
+    decay but only external source rates"""
+    # create transfer rate for U
+    vector = {'U235': 1}
+    external_source = 10  # grams
+    op = CoupledOperator(model, CHAIN_PATH)
+    integrator = openmc.deplete.PredictorIntegrator(
+        op, [1, 1], 0.0, timestep_units = 'd')
+    integrator.add_external_source_rate('f', vector, external_source/(24*3600))
+    integrator.integrate()
+
+    # Get number of U238 atoms from results
+    results = openmc.deplete.Results('depletion_results.h5')
+    _, atoms = results.get_atoms(model.materials[0], "U235")
+
+    # Ensure number of atoms equal external source
+    assert atoms[1] - atoms[0] == pytest.approx(
+        external_source * AVOGADRO / atomic_mass('U235'))
+    assert atoms[2] - atoms[1] == pytest.approx(
+        external_source * AVOGADRO / atomic_mass('U235'))
diff --git a/tests/unit_tests/test_deplete_transfer_rates.py b/tests/unit_tests/test_deplete_transfer_rates.py
index 140777cd6..a3228e9fb 100644
--- a/tests/unit_tests/test_deplete_transfer_rates.py
+++ b/tests/unit_tests/test_deplete_transfer_rates.py
@@ -16,6 +16,7 @@ CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
 
 @pytest.fixture
 def model():
+    openmc.reset_auto_ids()
     f = openmc.Material(name="f")
     f.add_element("U", 1, percent_type="ao", enrichment=4.25)
     f.add_element("O", 2)
@@ -54,32 +55,35 @@ def model():
 
     return openmc.Model(geometry, materials, settings)
 
-@pytest.mark.parametrize("case_name, transfer_rates", [
-    ('elements', {'U': 0.01, 'Xe': 0.1}),
-    ('nuclides', {'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01}),
+@pytest.mark.parametrize("case_name, transfer_rates, timesteps", [
+    ('elements', {'U': 0.01, 'Xe': 0.1}, None),
+    ('nuclides', {'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01}, None),
     ('nuclides_elements', {'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01, 'U': 0.01,
-                           'Xe': 0.1}),
+                           'Xe': 0.1}, None),
     ('elements_nuclides', {'U': 0.01, 'Xe': 0.1, 'I135': 0.01, 'Gd156': 0.1,
-                           'Gd157': 0.01}),
+                           'Gd157': 0.01}, None),
     ('multiple_transfer', {'U': 0.01, 'Xe': 0.1, 'I135': 0.01, 'Gd156': 0.1,
-                           'Gd157': 0.01}),
-    ('rates_invalid_1', {'Gd': 0.01, 'Gd157': 0.01, 'Gd156': 0.01}),
-    ('rates_invalid_2', {'Gd156': 0.01, 'Gd157': 0.01, 'Gd': 0.01}),
-    ('rates_invalid_3', {'Gb156': 0.01}),
-    ('rates_invalid_4', {'Gb': 0.01})
+                           'Gd157': 0.01}, None),
+    ('timesteps', {'U': 0.01, 'Xe': 0.1}, [1]),
+    ('rates_invalid_1', {'Gd': 0.01, 'Gd157': 0.01, 'Gd156': 0.01}, None),
+    ('rates_invalid_2', {'Gd156': 0.01, 'Gd157': 0.01, 'Gd': 0.01}, None),
+    ('rates_invalid_3', {'Gb156': 0.01}, None),
+    ('rates_invalid_4', {'Gb': 0.01}, None)
     ])
-def test_get_set(model, case_name, transfer_rates):
+def test_get_set(model, case_name, transfer_rates, timesteps):
     """Tests the get/set methods"""
-
-    openmc.reset_auto_ids()
     op = CoupledOperator(model, CHAIN_PATH)
-    transfer = TransferRates(op, model)
+    number_of_timesteps = 2
+    transfer = TransferRates(op, model.materials, number_of_timesteps)
+
+    if timesteps is None:
+        timesteps = np.arange(number_of_timesteps)
 
     # Test by Openmc material, material name and material id
     material, dest_material, dest_material2 = [m for m in model.materials
                                                if m.depletable]
     for material_input in [material, material.name, material.id]:
-        for dest_material_input in [dest_material, dest_material.name,
+        for dest_material_input in [None, dest_material, dest_material.name,
                                     dest_material.id]:
             if case_name == 'rates_invalid_1':
                 with pytest.raises(ValueError, match='Cannot add transfer '
@@ -117,14 +121,21 @@ def test_get_set(model, case_name, transfer_rates):
                 for component, transfer_rate in transfer_rates.items():
                     transfer.set_transfer_rate(material_input, [component],
                                                transfer_rate,
+                                               timesteps=timesteps,
                                                destination_material=\
                                                dest_material_input)
-                    assert transfer.get_transfer_rate(
-                        material_input, component)[0] == transfer_rate
-                    assert transfer.get_destination_material(
-                        material_input, component)[0] == str(dest_material.id)
-                assert transfer.get_components(material_input) == \
-                    transfer_rates.keys()
+                    assert transfer.get_external_rate(
+                        material_input, component, timesteps,
+                        dest_material_input)[0] == transfer_rate
+                    assert np.all(transfer.external_timesteps == timesteps)
+
+                if timesteps is not None:
+                    for timestep in timesteps:
+                        assert transfer.get_components(material_input, timestep,
+                            dest_material_input) == list(transfer_rates.keys())
+                else:
+                    assert transfer.get_components(material_input, timesteps,
+                            dest_material_input) == list(transfer_rates.keys())
 
             if case_name == 'multiple_transfer':
                 for dest2_material_input in [dest_material2, dest_material2.name,
@@ -135,10 +146,8 @@ def test_get_set(model, case_name, transfer_rates):
                                                destination_material=\
                                                dest2_material_input)
                         for id, dest_mat in zip([0,1],[dest_material,dest_material2]):
-                            assert transfer.get_transfer_rate(
-                                material_input, component)[id] == transfer_rate
-                            assert transfer.get_destination_material(
-                                material_input, component)[id] == str(dest_mat.id)
+                            assert transfer.get_external_rate(
+                                material_input, component, timesteps)[0] == transfer_rate
 
 @pytest.mark.parametrize("transfer_rate_units, unit_conv", [
     ('1/s', 1),
@@ -158,13 +167,15 @@ def test_units(transfer_rate_units, unit_conv, model):
     # create transfer rate Xe
     components = ['Xe', 'U235']
     transfer_rate = 1e-5
+    number_of_timesteps = 2
     op = CoupledOperator(model, CHAIN_PATH)
-    transfer = TransferRates(op, model)
+    transfer = TransferRates(op, model.materials, number_of_timesteps)
 
     for component in components:
         transfer.set_transfer_rate('f', [component], transfer_rate * unit_conv,
                                    transfer_rate_units=transfer_rate_units)
-        assert transfer.get_transfer_rate('f', component)[0] == transfer_rate
+        for timestep in range(transfer.number_of_timesteps):
+            assert transfer.get_external_rate('f', component, timestep)[0] == transfer_rate
 
 
 def test_transfer(run_in_tmpdir, model):