diff --git a/docs/source/io_formats/depletion_results.rst b/docs/source/io_formats/depletion_results.rst
index d35e251469..4f92384a0a 100644
--- a/docs/source/io_formats/depletion_results.rst
+++ b/docs/source/io_formats/depletion_results.rst
@@ -12,9 +12,10 @@ The current version of the depletion results file format is 1.0.
              - **version** (*int[2]*) -- Major and minor version of the
                statepoint file format.
 
-:Datasets: - **eigenvalues** (*double[][]*) -- k-eigenvalues at each
+:Datasets: - **eigenvalues** (*double[][][2]*) -- k-eigenvalues at each
              time/stage. This array has shape (number of timesteps, number of
-             stages).
+             stages, value). The last axis contains the eigenvalue and the 
+             associated uncertainty
            - **number** (*double[][][][]*) -- Total number of atoms. This array
              has shape (number of timesteps, number of stages, number of
              materials, number of nuclides).
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index ffbffbc734..21c96b0471 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -525,7 +525,8 @@ class Operator(TransportOperator):
         rates = self.reaction_rates
         rates[:, :, :] = 0.0
 
-        k_combined = openmc.capi.keff()[0]
+        # Get k and uncertainty
+        k_combined = openmc.capi.keff()
 
         # Extract tally bins
         materials = self.burnable_mats
diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
index 6170a44643..b0dc4f3bca 100644
--- a/openmc/deplete/results.py
+++ b/openmc/deplete/results.py
@@ -237,8 +237,8 @@ class Results(object):
                               chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
                               dtype='float64')
 
-        handle.create_dataset("eigenvalues", (1, n_stages),
-                              maxshape=(None, n_stages), dtype='float64')
+        handle.create_dataset("eigenvalues", (1, n_stages, 2),
+                              maxshape=(None, n_stages, 2), dtype='float64')
 
         handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
 
diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py
index 9fb6eec86c..7f33c0573f 100644
--- a/openmc/deplete/results_list.py
+++ b/openmc/deplete/results_list.py
@@ -91,11 +91,13 @@ class ResultsList(list):
         time : numpy.ndarray
             Array of times in [s]
         eigenvalue : numpy.ndarray
-            k-eigenvalue at each time
+            k-eigenvalue at each time. Column 0
+            contains the eigenvalue, while column
+            1 contains the associated uncertainty
 
         """
         time = np.empty_like(self, dtype=float)
-        eigenvalue = np.empty_like(self, dtype=float)
+        eigenvalue = np.empty((len(self), 2), dtype=float)
 
         # Get time/eigenvalue at each point
         for i, result in enumerate(self):
diff --git a/tests/regression_tests/deplete/test_reference.h5 b/tests/regression_tests/deplete/test_reference.h5
index 2f1ef375a9..3ee9d4447a 100644
Binary files a/tests/regression_tests/deplete/test_reference.h5 and b/tests/regression_tests/deplete/test_reference.h5 differ
diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py
index 8cade45b14..40090264ae 100644
--- a/tests/unit_tests/test_deplete_integrator.py
+++ b/tests/unit_tests/test_deplete_integrator.py
@@ -64,8 +64,9 @@ def test_results_save(run_in_tmpdir):
         r1[:] = np.random.rand(2, 2, 2)
 
     # Create global terms
-    eigvl1 = np.random.rand(stages)
-    eigvl2 = np.random.rand(stages)
+    # Col 0: eig, Col 1: uncertainty
+    eigvl1 = np.random.rand(stages, 2)
+    eigvl2 = np.random.rand(stages, 2)
 
     eigvl1 = comm.bcast(eigvl1, root=0)
     eigvl2 = comm.bcast(eigvl2, root=0)
diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py
index 5a066be6e4..c4b7aceba5 100644
--- a/tests/unit_tests/test_deplete_resultslist.py
+++ b/tests/unit_tests/test_deplete_resultslist.py
@@ -44,6 +44,8 @@ def test_get_eigenvalue(res):
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
     k_ref = [1.16984322, 1.19097427, 1.03012572, 1.20045627]
+    u_ref = [0.0375587, 0.0347639, 0.07216021, 0.02839642]
 
     np.testing.assert_allclose(t, t_ref)
-    np.testing.assert_allclose(k, k_ref)
+    np.testing.assert_allclose(k[:, 0], k_ref)
+    np.testing.assert_allclose(k[:, 1], u_ref)