diff --git a/datalab/datalab_session/data_operations/median.py b/datalab/datalab_session/data_operations/median.py
index 5ac97fb..fd552fe 100644
--- a/datalab/datalab_session/data_operations/median.py
+++ b/datalab/datalab_session/data_operations/median.py
@@ -57,7 +57,9 @@ def operate(self):
         # using the numpy library's median method
         median = np.median(stacked_data, axis=2)
 
-        output = create_output(self.cache_key, median)
+        self.set_operation_progress(0.80)
+
+        output = create_output(self.cache_key, median, comment=f'Product of Datalab Median on files {", ".join([image["basename"] for image in input])}')
 
         self.set_output(output)
         log.info(f'Median output: {self.get_output()}')
diff --git a/datalab/datalab_session/data_operations/normalization.py b/datalab/datalab_session/data_operations/normalization.py
index f842244..9698117 100644
--- a/datalab/datalab_session/data_operations/normalization.py
+++ b/datalab/datalab_session/data_operations/normalization.py
@@ -52,7 +52,7 @@ def operate(self):
             median = np.median(image)
             normalized_image = image / median
 
-            output = create_output(self.cache_key, normalized_image, index=index)
+            output = create_output(self.cache_key, normalized_image, index=index, comment=f'Product of Datalab Normalization on file {input[index]["basename"]}')
             output_files.append(output)
 
             self.set_operation_progress(self.get_operation_progress() + .40 * (index + 1) / len(input))
diff --git a/datalab/datalab_session/data_operations/rgb_stack.py b/datalab/datalab_session/data_operations/rgb_stack.py
index 46e238c..1610bd1 100644
--- a/datalab/datalab_session/data_operations/rgb_stack.py
+++ b/datalab/datalab_session/data_operations/rgb_stack.py
@@ -78,7 +78,8 @@ def operate(self):
 
         self.set_operation_progress(0.8)
         
-        output =  create_output(self.cache_key, stacked_data, large_jpg=large_jpg_path, small_jpg=small_jpg_path)
+        rgb_comment = f'Product of Datalab RGB Stack on files {", ".join([image["basename"] for image in rgb_input_list])}'
+        output =  create_output(self.cache_key, stacked_data, large_jpg=large_jpg_path, small_jpg=small_jpg_path, comment=rgb_comment)
 
         self.set_operation_progress(1.0)
         self.set_output(output)
diff --git a/datalab/datalab_session/data_operations/stacking.py b/datalab/datalab_session/data_operations/stacking.py
index 096bb0e..741b411 100644
--- a/datalab/datalab_session/data_operations/stacking.py
+++ b/datalab/datalab_session/data_operations/stacking.py
@@ -63,7 +63,8 @@ def operate(self):
         
         self.set_operation_progress(0.8)
 
-        output = create_output(self.cache_key, stacked_sum)
+        stacking_comment = f'Product of Datalab Stacking. Stack of {", ".join([image["basename"] for image in input_files])}'
+        output = create_output(self.cache_key, stacked_sum, comment=stacking_comment)
 
         self.set_output(output)
         log.info(f'Stacked output: {self.get_output()}')
diff --git a/datalab/datalab_session/data_operations/subtraction.py b/datalab/datalab_session/data_operations/subtraction.py
index 6a41291..5c42ff9 100644
--- a/datalab/datalab_session/data_operations/subtraction.py
+++ b/datalab/datalab_session/data_operations/subtraction.py
@@ -75,7 +75,8 @@ def operate(self):
 
             difference_array = np.subtract(input_image, subtraction_image)
 
-            outputs.append(create_output(self.cache_key, difference_array, index=index))
+            subtraction_comment = f'Product of Datalab Subtraction of {subtraction_file_input[0]["basename"]} subtracted from {input_files[index]["basename"]}'
+            outputs.append(create_output(self.cache_key, difference_array, index=index, comment=subtraction_comment))
 
             self.set_operation_progress(self.get_operation_progress() + .50 * (index + 1) / len(input_files))
 
diff --git a/datalab/datalab_session/file_utils.py b/datalab/datalab_session/file_utils.py
index 164b677..0c04094 100644
--- a/datalab/datalab_session/file_utils.py
+++ b/datalab/datalab_session/file_utils.py
@@ -27,13 +27,14 @@ def get_hdu(path: str, extension: str = 'SCI') -> list[fits.HDUList]:
 def get_fits_dimensions(fits_file, extension: str = 'SCI') -> tuple:
   return fits.open(fits_file)[extension].shape
 
-def create_fits(key: str, image_arr: np.ndarray) -> str:
+def create_fits(key: str, image_arr: np.ndarray, comment=None) -> str:
   """
   Creates a fits file with the given key and image array
   Returns the the path to the fits_file
   """
 
   header = fits.Header([('KEY', key)])
+  header.add_comment(comment) if comment else None
   primary_hdu = fits.PrimaryHDU(header=header)
   image_hdu = fits.ImageHDU(data=image_arr, name='SCI')
 
@@ -109,7 +110,7 @@ def scale_points(height_1: int, width_1: int, height_2: int, width_2: int, x_poi
 
   return x_points, y_points
 
-def create_output(cache_key, np_array=None, fits_file=None, large_jpg=None, small_jpg=None, index=None):
+def create_output(cache_key, np_array=None, fits_file=None, large_jpg=None, small_jpg=None, index=None, comment=None):
   """
   A more automated way of creating output for a dev
   Dev can specify just a cache_key and np array and the function will create the fits and jpgs
@@ -117,7 +118,7 @@ def create_output(cache_key, np_array=None, fits_file=None, large_jpg=None, smal
   """
 
   if np_array is not None and fits_file is None:
-    fits_file = create_fits(cache_key, np_array)
+    fits_file = create_fits(cache_key, np_array, comment)
 
   if not large_jpg or not small_jpg:
     large_jpg, small_jpg = create_jpgs(cache_key, fits_file)