diff --git a/asltk/asldata.py b/asltk/asldata.py
index a72a56f..8cc3346 100644
--- a/asltk/asldata.py
+++ b/asltk/asldata.py
@@ -86,6 +86,8 @@ def __init__(
 
         if kwargs.get('m0') is not None:
             average_m0 = kwargs.get('average_m0', False)
+            if self._asl_image:
+                self._asl_image._average_m0 = average_m0
 
             if isinstance(kwargs.get('m0'), str):
                 m0_path = kwargs.get('m0')
diff --git a/asltk/reconstruction/cbf_mapping.py b/asltk/reconstruction/cbf_mapping.py
index 22d72be..8fb703d 100644
--- a/asltk/reconstruction/cbf_mapping.py
+++ b/asltk/reconstruction/cbf_mapping.py
@@ -14,7 +14,7 @@
 from asltk.logging_config import get_logger, log_processing_step
 from asltk.models.signal_dynamic import asl_model_buxton
 from asltk.mri_parameters import MRIParameters
-from asltk.utils.io import ImageIO
+from asltk.utils.io import ImageIO, clone_image
 
 # Global variables to assist multi cpu threading
 cbf_map = None
@@ -406,14 +406,22 @@ def create_map(
         )
 
         # Prepare output maps
-        cbf_map_image = ImageIO(self._asl_data('m0').get_image_path())
-        cbf_map_image.update_image_data(self._cbf_map)
+        base_volume = ImageIO(self._asl_data('m0').get_image_path())
+        cbf_map_image = clone_image(base_volume)
+        cbf_map_image.update_image_data(
+            self._cbf_map, self._asl_data._asl_image._average_m0
+        )
 
-        cbf_map_norm_image = ImageIO(self._asl_data('m0').get_image_path())
-        cbf_map_norm_image.update_image_data(self._cbf_map * (60 * 60 * 1000))
+        cbf_map_norm_image = clone_image(base_volume)
+        cbf_map_norm_image.update_image_data(
+            self._cbf_map * (60 * 60 * 1000),
+            self._asl_data._asl_image._average_m0,
+        )
 
-        att_map_image = ImageIO(self._asl_data('m0').get_image_path())
-        att_map_image.update_image_data(self._att_map)
+        att_map_image = clone_image(base_volume)
+        att_map_image.update_image_data(
+            self._att_map, self._asl_data._asl_image._average_m0
+        )
 
         output_maps = {
             'cbf': cbf_map_image,
diff --git a/asltk/reconstruction/multi_te_mapping.py b/asltk/reconstruction/multi_te_mapping.py
index fc6adc5..9bbb66f 100644
--- a/asltk/reconstruction/multi_te_mapping.py
+++ b/asltk/reconstruction/multi_te_mapping.py
@@ -12,7 +12,7 @@
 from asltk.models.signal_dynamic import asl_model_multi_te
 from asltk.mri_parameters import MRIParameters
 from asltk.reconstruction import CBFMapping
-from asltk.utils.io import ImageIO
+from asltk.utils.io import ImageIO, clone_image
 
 # Global variables to assist multi cpu threading
 cbf_map = None
@@ -388,19 +388,27 @@ def create_map(
             )
 
             # Prepare output maps
-            cbf_map_image = ImageIO(self._asl_data('m0').get_image_path())
-            cbf_map_image.update_image_data(self._cbf_map)
+            base_volume = ImageIO(self._asl_data('m0').get_image_path())
+            cbf_map_image = clone_image(base_volume)
+            cbf_map_image.update_image_data(
+                self._cbf_map, self._asl_data._asl_image._average_m0
+            )
 
-            cbf_map_norm_image = ImageIO(self._asl_data('m0').get_image_path())
+            cbf_map_norm_image = clone_image(base_volume)
             cbf_map_norm_image.update_image_data(
-                self._cbf_map * (60 * 60 * 1000)
+                self._cbf_map * (60 * 60 * 1000),
+                self._asl_data._asl_image._average_m0,
             )
 
-            att_map_image = ImageIO(self._asl_data('m0').get_image_path())
-            att_map_image.update_image_data(self._att_map)
+            att_map_image = clone_image(base_volume)
+            att_map_image.update_image_data(
+                self._att_map, self._asl_data._asl_image._average_m0
+            )
 
-            t1blgm_map_image = ImageIO(self._asl_data('m0').get_image_path())
-            t1blgm_map_image.update_image_data(self._t1blgm_map)
+            t1blgm_map_image = clone_image(base_volume)
+            t1blgm_map_image.update_image_data(
+                self._t1blgm_map, self._asl_data._asl_image._average_m0
+            )
 
             # Create output maps dictionary
             output_maps = {
diff --git a/asltk/reconstruction/t2_mapping.py b/asltk/reconstruction/t2_mapping.py
index 03f96c5..28cb2a0 100644
--- a/asltk/reconstruction/t2_mapping.py
+++ b/asltk/reconstruction/t2_mapping.py
@@ -184,7 +184,7 @@ def create_map(
             )
 
             # Prepare output maps
-            # TODO At the moment, the T2 maps and mean T2 maps are as ImageIO object, however, the Spacing, Dimension are not given as a 4D array. Ceck if can be imported from the m0 image is 3D.
+            # TODO At the moment, the T2 maps and mean T2 maps are as ImageIO object, however, the Spacing, Dimension are not given as a 4D array. Check if can be imported from the m0 image is 3D.
             t2_maps_image = ImageIO(
                 image_array=np.array(
                     [
diff --git a/tests/reconstruction/test_cbf_mapping.py b/tests/reconstruction/test_cbf_mapping.py
index ea7a253..d8db922 100644
--- a/tests/reconstruction/test_cbf_mapping.py
+++ b/tests/reconstruction/test_cbf_mapping.py
@@ -194,3 +194,29 @@ def test_cbf_map_normalized_flag_true_result_cbf_map_rescaled():
     out_norm_array[out_norm_array == 0] = np.nan
     mean_px_value = np.nanmean(out_norm_array)
     assert mean_px_value < 500 and mean_px_value > 50
+
+
+def test_cbf_object_create_map_success_using_average_m0_option(tmp_path):
+    base_m0 = ImageIO(M0).get_as_numpy()
+    m04d = np.stack([base_m0 for _ in range(3)], axis=0)
+    # Save m04d into a temporary .nii.gz file
+    ImageIO(image_array=m04d).save_image(
+        os.path.join(tmp_path, 'm0_4d.nii.gz')
+    )
+
+    asldata4d = ASLData(
+        pcasl=PCASL_MTE,
+        m0=os.path.join(tmp_path, 'm0_4d.nii.gz'),
+        ld_values=[100.0, 100.0, 150.0, 150.0, 400.0, 800.0, 1800.0],
+        pld_values=[170.0, 270.0, 370.0, 520.0, 670.0, 1070.0, 1870.0],
+        average_m0=True,
+    )
+    mte = CBFMapping(asldata4d)
+    mask = ImageIO(M0_BRAIN_MASK)
+
+    mte.set_brain_mask(mask)
+
+    output = mte.create_map()
+
+    assert output['cbf'].get_as_numpy().shape == base_m0.shape
+    assert output['att'].get_as_numpy().shape == base_m0.shape
diff --git a/tests/reconstruction/test_multi_te_mapping.py b/tests/reconstruction/test_multi_te_mapping.py
index 5d1dc22..a325839 100644
--- a/tests/reconstruction/test_multi_te_mapping.py
+++ b/tests/reconstruction/test_multi_te_mapping.py
@@ -209,3 +209,42 @@ def test_multite_asl_object_create_map_using_provided_cbf_att_maps(capfd):
     if re.search('multiTE-ASL', out):
         test_pass = True
     assert test_pass
+
+
+def test_multite_asl_object_create_map_success_using_average_m0_option(
+    tmp_path,
+):
+    base_m0 = ImageIO(M0).get_as_numpy()
+    m04d = np.stack([base_m0 for _ in range(3)], axis=0)
+    # Save m04d into a temporary .nii.gz file
+    ImageIO(image_array=m04d).save_image(
+        os.path.join(tmp_path, 'm0_4d.nii.gz')
+    )
+
+    asldata4d = ASLData(
+        pcasl=PCASL_MTE,
+        m0=os.path.join(tmp_path, 'm0_4d.nii.gz'),
+        ld_values=[100.0, 100.0, 150.0, 150.0, 400.0, 800.0, 1800.0],
+        pld_values=[170.0, 270.0, 370.0, 520.0, 670.0, 1070.0, 1870.0],
+        te_values=[
+            13.56,
+            67.82,
+            122.08,
+            176.33,
+            230.59,
+            284.84,
+            339.100,
+            393.36,
+        ],
+        average_m0=True,
+    )
+    mte = MultiTE_ASLMapping(asldata4d)
+    mask = ImageIO(M0_BRAIN_MASK)
+
+    mte.set_brain_mask(mask)
+
+    output = mte.create_map()
+
+    assert output['cbf'].get_as_numpy().shape == base_m0.shape
+    assert output['att'].get_as_numpy().shape == base_m0.shape
+    assert output['t1blgm'].get_as_numpy().shape == base_m0.shape