diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py
index 32007dcb4..a6d994c02 100644
--- a/festim/exports/derived_quantities/average_surface.py
+++ b/festim/exports/derived_quantities/average_surface.py
@@ -3,9 +3,41 @@
 
 
 class AverageSurface(SurfaceQuantity):
+    """
+    Computes the average value of a field on a given surface
+    int(f ds) / int (1 * ds)
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the field
+
+    Notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+
+    """
+
     def __init__(self, field, surface) -> None:
         super().__init__(field=field, surface=surface)
-        self.title = "Average {} surface {}".format(self.field, self.surface)
+
+    @property
+    def title(self):
+        quantity_title = f"Average {self.field} surface {self.surface}"
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (K)"
+            else:
+                return quantity_title + " (H m-3)"
+        else:
+            return quantity_title
 
     def compute(self):
         return f.assemble(self.function * self.ds(self.surface)) / f.assemble(
diff --git a/festim/exports/derived_quantities/average_volume.py b/festim/exports/derived_quantities/average_volume.py
index 49372e721..9c4140161 100644
--- a/festim/exports/derived_quantities/average_volume.py
+++ b/festim/exports/derived_quantities/average_volume.py
@@ -3,9 +3,40 @@
 
 
 class AverageVolume(VolumeQuantity):
+    """
+    Computes the average value of a field in a given volume
+    int(f dx) / int (1 * dx)
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the field
+
+    Notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+    """
+
     def __init__(self, field, volume: int) -> None:
-        super().__init__(field, volume)
-        self.title = "Average {} volume {}".format(self.field, self.volume)
+        super().__init__(field=field, volume=volume)
+
+    @property
+    def title(self):
+        quantity_title = f"Average {self.field} volume {self.volume}"
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (K)"
+            else:
+                return quantity_title + " (H m-3)"
+        else:
+            return quantity_title
 
     def compute(self):
         return f.assemble(self.function * self.dx(self.volume)) / f.assemble(
diff --git a/festim/exports/derived_quantities/derived_quantities.py b/festim/exports/derived_quantities/derived_quantities.py
index 0cb0e108d..4664f0ad5 100644
--- a/festim/exports/derived_quantities/derived_quantities.py
+++ b/festim/exports/derived_quantities/derived_quantities.py
@@ -23,6 +23,8 @@ class DerivedQuantities(list):
         nb_iterations_between_exports (int, optional): number of
             iterations between each export. If None, the file will be
             exported at the last timestep. Defaults to None.
+        show_units (bool, optional): will show the units of each
+            derived quantity in the title in export
     """
 
     def __init__(
@@ -31,6 +33,7 @@ def __init__(
         filename: str = None,
         nb_iterations_between_compute: int = 1,
         nb_iterations_between_exports: int = None,
+        show_units=False,
     ) -> None:
         # checks that input is list
         if len(args) == 0:
@@ -43,8 +46,9 @@ def __init__(
         self.filename = filename
         self.nb_iterations_between_compute = nb_iterations_between_compute
         self.nb_iterations_between_exports = nb_iterations_between_exports
+        self.show_units = show_units
 
-        self.data = [self.make_header()]
+        self.data = []
         self.t = []
 
     @property
@@ -108,6 +112,12 @@ def filename(self, value):
     def make_header(self):
         header = ["t(s)"]
         for quantity in self:
+            quantity.show_units = self.show_units
+            if self.show_units is False:
+                warnings.warn(
+                    "The current derived_quantities title style will be deprecated in a future release, please use show_units=True instead",
+                    DeprecationWarning,
+                )
             header.append(quantity.title)
         return header
 
@@ -139,6 +149,11 @@ def compute(self, t):
                 value = quantity.compute(self.volume_markers)
             else:
                 value = quantity.compute()
+
+            # check if first time writing data
+            if len(self.data) == 0:
+                self.data = [self.make_header()]
+
             quantity.data.append(value)
             quantity.t.append(t)
             row.append(value)
diff --git a/festim/exports/derived_quantities/derived_quantity.py b/festim/exports/derived_quantities/derived_quantity.py
index 113ac26ba..489ab3fa6 100644
--- a/festim/exports/derived_quantities/derived_quantity.py
+++ b/festim/exports/derived_quantities/derived_quantity.py
@@ -3,6 +3,8 @@
 
 class DerivedQuantity(Export):
     """
+    Parent class of all derived quantities
+
     Args:
         field (str, int):  the field ("solute", 0, 1, "T", "retention")
     """
@@ -18,16 +20,19 @@ def __init__(self, field) -> None:
         self.Q = None
         self.data = []
         self.t = []
+        self.show_units = False
 
 
 class VolumeQuantity(DerivedQuantity):
-    def __init__(self, field: str or int, volume: int) -> None:
-        """DerivedQuantity relative to a volume
+    """DerivedQuantity relative to a volume
+
+    Args:
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
 
-        Args:
-            field (str, int): the field ("solute", 0, 1, "T", "retention")
-            volume (int): the volume id
-        """
+    """
+
+    def __init__(self, field: str or int, volume: int) -> None:
         super().__init__(field)
         self.volume = volume
 
@@ -44,13 +49,16 @@ def volume(self, value):
 
 
 class SurfaceQuantity(DerivedQuantity):
+    """DerivedQuantity relative to a surface
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    """
+
     def __init__(self, field: str or int, surface: int) -> None:
-        """DerivedQuantity relative to a surface
 
-        Args:
-            field (str, int):  the field ("solute", 0, 1, "T", "retention")
-            surface (int): the surface id
-        """
         super().__init__(field)
         self.surface = surface
 
diff --git a/festim/exports/derived_quantities/hydrogen_flux.py b/festim/exports/derived_quantities/hydrogen_flux.py
index 09f6b9d70..26f39626c 100644
--- a/festim/exports/derived_quantities/hydrogen_flux.py
+++ b/festim/exports/derived_quantities/hydrogen_flux.py
@@ -2,7 +2,25 @@
 
 
 class HydrogenFlux(SurfaceFlux):
-    """Equivalent to SurfaceFlux("solute", ...)"""
+    """
+    Computes the surface flux of hydrogen at a given surface
+
+    Args:
+        surface (int): the surface id
+
+    Attribtutes
+        field (str): the hydrogen solute field
+        surface (int): the surface id
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the hydrogen solute field
+
+    Notes:
+        units are in H/m2/s in 1D, H/m/s in 2D and H/s in 3D domains
+
+    """
 
     def __init__(self, surface) -> None:
         super().__init__(field="solute", surface=surface)
diff --git a/festim/exports/derived_quantities/maximum_surface.py b/festim/exports/derived_quantities/maximum_surface.py
index 35eab5f16..f4b790755 100644
--- a/festim/exports/derived_quantities/maximum_surface.py
+++ b/festim/exports/derived_quantities/maximum_surface.py
@@ -1,20 +1,43 @@
-from festim import DerivedQuantity
+from festim import SurfaceQuantity
 import fenics as f
 import numpy as np
 
 
-class MaximumSurface(DerivedQuantity):
+class MaximumSurface(SurfaceQuantity):
     """
+    Computes the maximum value of a field on a given surface
+
     Args:
-        field (str): the field from which the maximum
-            is computed (ex: "solute", "retention", "T"...)
-        surface (int): the surface id where the maximum is computed
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the field
+
+    Notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+
     """
 
     def __init__(self, field, surface) -> None:
-        super().__init__(field)
-        self.surface = surface
-        self.title = "Maximum {} surface {}".format(self.field, self.surface)
+        super().__init__(field=field, surface=surface)
+
+    @property
+    def title(self):
+        quantity_title = f"Maximum {self.field} surface {self.surface}"
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (K)"
+            else:
+                return quantity_title + " (H m-3)"
+        else:
+            return quantity_title
 
     def compute(self, surface_markers):
         """Maximum of f over subdomains facets marked with self.surface"""
diff --git a/festim/exports/derived_quantities/maximum_volume.py b/festim/exports/derived_quantities/maximum_volume.py
index 1c7c66d75..7dc7f0b18 100644
--- a/festim/exports/derived_quantities/maximum_volume.py
+++ b/festim/exports/derived_quantities/maximum_volume.py
@@ -4,9 +4,40 @@
 
 
 class MaximumVolume(VolumeQuantity):
+    """
+    Computes the maximum value of a field in a given volume
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function for
+            the field
+
+    notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+
+    """
+
     def __init__(self, field, volume) -> None:
-        super().__init__(field, volume=volume)
-        self.title = "Maximum {} volume {}".format(self.field, self.volume)
+        super().__init__(field=field, volume=volume)
+
+    @property
+    def title(self):
+        quantity_title = f"Maximum {self.field} volume {self.volume}"
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (K)"
+            else:
+                return quantity_title + " (H m-3)"
+        else:
+            return quantity_title
 
     def compute(self, volume_markers):
         """Minimum of f over subdomains cells marked with self.volume"""
diff --git a/festim/exports/derived_quantities/minimum_surface.py b/festim/exports/derived_quantities/minimum_surface.py
index 2ea9c7bc5..f782db5f0 100644
--- a/festim/exports/derived_quantities/minimum_surface.py
+++ b/festim/exports/derived_quantities/minimum_surface.py
@@ -1,20 +1,43 @@
-from festim import DerivedQuantity
+from festim import SurfaceQuantity
 import fenics as f
 import numpy as np
 
 
-class MinimumSurface(DerivedQuantity):
+class MinimumSurface(SurfaceQuantity):
     """
+    Computes the minimum value of a field on a given surface
+
     Args:
-        field (str): the field from which the minimum
-            is computed (ex: "solute", "retention", "T"...)
-        surface (int): the surface id where the minimum is computed
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the field
+
+    Notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+
     """
 
     def __init__(self, field, surface) -> None:
-        super().__init__(field)
-        self.surface = surface
-        self.title = "Minimum {} surface {}".format(self.field, self.surface)
+        super().__init__(field=field, surface=surface)
+
+    @property
+    def title(self):
+        quantity_title = f"Minimum {self.field} surface {self.surface}"
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (K)"
+            else:
+                return quantity_title + " (H m-3)"
+        else:
+            return quantity_title
 
     def compute(self, surface_markers):
         """Minimum of f over subdomains facets marked with self.surface"""
diff --git a/festim/exports/derived_quantities/minimum_volume.py b/festim/exports/derived_quantities/minimum_volume.py
index 02b064789..b11e6bb21 100644
--- a/festim/exports/derived_quantities/minimum_volume.py
+++ b/festim/exports/derived_quantities/minimum_volume.py
@@ -4,9 +4,40 @@
 
 
 class MinimumVolume(VolumeQuantity):
+    """
+    Computes the minimum value of a field in a given volume
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function for
+            the field
+
+    notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+
+    """
+
     def __init__(self, field, volume) -> None:
-        super().__init__(field, volume=volume)
-        self.title = "Minimum {} volume {}".format(self.field, self.volume)
+        super().__init__(field=field, volume=volume)
+
+    @property
+    def title(self):
+        quantity_title = f"Minimum {self.field} volume {self.volume}"
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (K)"
+            else:
+                return quantity_title + " (H m-3)"
+        else:
+            return quantity_title
 
     def compute(self, volume_markers):
         """Minimum of f over subdomains cells marked with self.volume"""
diff --git a/festim/exports/derived_quantities/point_value.py b/festim/exports/derived_quantities/point_value.py
index 10b046560..80d5a82e0 100644
--- a/festim/exports/derived_quantities/point_value.py
+++ b/festim/exports/derived_quantities/point_value.py
@@ -2,20 +2,44 @@
 
 
 class PointValue(DerivedQuantity):
-    """DerivedQuantity relative to a point
+    """
+    Computes the value of a field at a given point
 
     Args:
-        field (str, int): the field ("solute", 0, 1, "T", "retention")
-        point (int, float, tuple, list): the point coordinates
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        x (int, float, tuple, list): the point coordinates
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        x (int, float, tuple, list): the point coordinates
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the field
+
+    Notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+
     """
 
     def __init__(self, field: str or int, x: int or float or tuple or list) -> None:
-        super().__init__(field)
+        super().__init__(field=field)
         # make sure x is an iterable
         if not hasattr(x, "__iter__"):
             x = [x]
         self.x = x
-        self.title = "{} value at {}".format(field, x)
+
+    @property
+    def title(self):
+        quantity_title = f"{self.field} value at {self.x}"
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (K)"
+            else:
+                return quantity_title + " (H m-3)"
+        else:
+            return quantity_title
 
     def compute(self):
         """The value at the point"""
diff --git a/festim/exports/derived_quantities/surface_flux.py b/festim/exports/derived_quantities/surface_flux.py
index 227da62e3..7c1d73c21 100644
--- a/festim/exports/derived_quantities/surface_flux.py
+++ b/festim/exports/derived_quantities/surface_flux.py
@@ -4,24 +4,61 @@
 
 
 class SurfaceFlux(SurfaceQuantity):
-    def __init__(self, field, surface) -> None:
-        """
+    """
+    Computes the surface flux of a field at a given surface in cartesian coordinates
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attribtutes
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        export_unit (str): the unit of the derived quantity in the export file
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the field
+
+    Notes:
         Object to compute the flux J of a field u through a surface
-        J = integral(-prop * grad(u) . n ds)
+        J = integral(+prop * grad(u) . n ds)
         where prop is the property of the field (D, thermal conductivity, etc)
         u is the field
         n is the normal vector of the surface
         ds is the surface measure.
+        units are in H/m2/s in 1D, H/m/s in 2D and H/s in 3D domains for hydrogen
+        concentration and W/m2 in 1D, W/m in 2D and W in 3D domains for temperature
 
-        Note: this will probably won't work correctly for axisymmetric meshes.
-        Use only with cartesian coordinates.
+    """
 
-        Args:
-            field (str, int):  the field ("solute", 0, 1, "T", "retention")
-            surface (int): the surface id
-        """
+    def __init__(self, field, surface) -> None:
         super().__init__(field=field, surface=surface)
-        self.title = "Flux surface {}: {}".format(self.surface, self.field)
+
+    @property
+    def export_unit(self):
+        # obtain domain dimension
+        dim = self.function.function_space().mesh().topology().dim()
+        # return unit depending on field and dimension of domain
+        if self.field == "T":
+            return f"W m{dim-3}".replace(" m0", "")
+        else:
+            return f"H m{dim-3} s-1".replace(" m0", "")
+
+    @property
+    def title(self):
+        quantity_title = f"Flux surface {self.surface}: {self.field}"
+        if self.show_units:
+            if self.field == "T":
+                quantity_title = f"Heat flux surface {self.surface}"
+            else:
+                quantity_title = f"{self.field} flux surface {self.surface}"
+
+            return quantity_title + f" ({self.export_unit})"
+
+        else:
+            return quantity_title
 
     @property
     def prop(self):
@@ -69,10 +106,25 @@ class SurfaceFluxCylindrical(SurfaceFlux):
     """
 
     def __init__(self, field, surface, azimuth_range=(0, 2 * np.pi)) -> None:
-        super().__init__(field, surface)
+        super().__init__(field=field, surface=surface)
         self.r = None
         self.azimuth_range = azimuth_range
 
+    @property
+    def title(self):
+        if self.field == "T":
+            quantity_title = f"Heat flux surface {self.surface}"
+        else:
+            quantity_title = f"{self.field} flux surface {self.surface}"
+
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (W)"
+            else:
+                return quantity_title + " (H s-1)"
+        else:
+            return quantity_title
+
     @property
     def azimuth_range(self):
         return self._azimuth_range
@@ -134,11 +186,26 @@ class SurfaceFluxSpherical(SurfaceFlux):
     def __init__(
         self, field, surface, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi)
     ) -> None:
-        super().__init__(field, surface)
+        super().__init__(field=field, surface=surface)
         self.r = None
         self.polar_range = polar_range
         self.azimuth_range = azimuth_range
 
+    @property
+    def title(self):
+        if self.field == "T":
+            quantity_title = f"Heat flux surface {self.surface}"
+        else:
+            quantity_title = f"{self.field} flux surface {self.surface}"
+
+        if self.show_units:
+            if self.field == "T":
+                return quantity_title + " (W)"
+            else:
+                return quantity_title + " (H s-1)"
+        else:
+            return quantity_title
+
     @property
     def polar_range(self):
         return self._polar_range
diff --git a/festim/exports/derived_quantities/thermal_flux.py b/festim/exports/derived_quantities/thermal_flux.py
index 176df3799..8dd4b4032 100644
--- a/festim/exports/derived_quantities/thermal_flux.py
+++ b/festim/exports/derived_quantities/thermal_flux.py
@@ -2,5 +2,25 @@
 
 
 class ThermalFlux(SurfaceFlux):
+    """
+    Computes the surface flux of heat at a given surface
+
+    Args:
+        surface (int): the surface id
+
+    Attribtutes
+        surface (int): the surface id
+        field (str): the temperature field
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the temperature field
+
+    Notes:
+        units are in W/m2 in 1D, W/m in 2D and W in 3D domains
+
+    """
+
     def __init__(self, surface) -> None:
         super().__init__(field="T", surface=surface)
diff --git a/festim/exports/derived_quantities/total_surface.py b/festim/exports/derived_quantities/total_surface.py
index 057d7f9af..649c1308b 100644
--- a/festim/exports/derived_quantities/total_surface.py
+++ b/festim/exports/derived_quantities/total_surface.py
@@ -3,9 +3,50 @@
 
 
 class TotalSurface(SurfaceQuantity):
+    """
+    Computes the total value of a field on a given surface
+    int(f ds)
+
+    Args:
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attribtutes
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        export_unit (str): the unit of the derived quantity for exporting
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the hydrogen solute field
+
+    Notes:
+        units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen
+        concentration and K in 1D, K m in 2D and K m2 in 3D domains for temperature
+
+    """
+
     def __init__(self, field, surface) -> None:
-        super().__init__(field, surface=surface)
-        self.title = "Total {} surface {}".format(self.field, self.surface)
+        super().__init__(field=field, surface=surface)
+
+    @property
+    def export_unit(self):
+        # obtain domain dimension
+        dim = self.function.function_space().mesh().topology().dim()
+        # return unit depending on field and dimension of domain
+        if self.field == "T":
+            return f"K m{dim-1}".replace(" m0", "").replace(" m1", " m")
+        else:
+            return f"H m{dim-3}".replace(" m0", "")
+
+    @property
+    def title(self):
+        quantity_title = f"Total {self.field} surface {self.surface}"
+        if self.show_units:
+            return quantity_title + f" ({self.export_unit})"
+        else:
+            return quantity_title
 
     def compute(self):
         return f.assemble(self.function * self.ds(self.surface))
diff --git a/festim/exports/derived_quantities/total_volume.py b/festim/exports/derived_quantities/total_volume.py
index 497db138b..164c55f54 100644
--- a/festim/exports/derived_quantities/total_volume.py
+++ b/festim/exports/derived_quantities/total_volume.py
@@ -3,9 +3,50 @@
 
 
 class TotalVolume(VolumeQuantity):
+    """
+    Computes the total value of a field in a given volume
+    int(f dx)
+
+    Args:
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+
+    Attribtutes
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        export_unit (str): the unit of the derived quantity for exporting
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the hydrogen solute field
+
+    Notes:
+        units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen
+        concentration and K m in 1D, K m2 in 2D and K m3 in 3D domains for temperature
+
+    """
+
     def __init__(self, field, volume) -> None:
-        super().__init__(field, volume=volume)
-        self.title = "Total {} volume {}".format(self.field, self.volume)
+        super().__init__(field=field, volume=volume)
+
+    @property
+    def export_unit(self):
+        # obtain domain dimension
+        dim = self.function.function_space().mesh().topology().dim()
+        # return unit depending on field and dimension of domain
+        if self.field == "T":
+            return f"K m{dim}".replace("1", "")
+        else:
+            return f"H m{dim-3}".replace(" m0", "")
+
+    @property
+    def title(self):
+        quantity_title = f"Total {self.field} volume {self.volume}"
+        if self.show_units:
+            return quantity_title + f" ({self.export_unit})"
+        else:
+            return quantity_title
 
     def compute(self):
         return f.assemble(self.function * self.dx(self.volume))
diff --git a/festim/exports/exports.py b/festim/exports/exports.py
index a84e80d9e..b79cd2f2b 100644
--- a/festim/exports/exports.py
+++ b/festim/exports/exports.py
@@ -136,6 +136,6 @@ def initialise_derived_quantities(self, dx, ds, materials):
         """
         for export in self:
             if isinstance(export, festim.DerivedQuantities):
-                export.data = [export.make_header()]
+                export.data = []
                 export.assign_measures_to_quantities(dx, ds)
                 export.assign_properties_to_quantities(materials)
diff --git a/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py b/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py
index 1a81e91c5..b6b198df4 100644
--- a/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py
+++ b/test/unit/test_exports/test_derived_quantities/test_derived_quantities.py
@@ -6,6 +6,12 @@
     TotalVolume,
     MaximumVolume,
     MinimumVolume,
+    MinimumSurface,
+    MaximumSurface,
+    AverageSurface,
+    PointValue,
+    SurfaceFluxCylindrical,
+    SurfaceFluxSpherical,
     Materials,
 )
 import fenics as f
@@ -23,7 +29,21 @@ class TestMakeHeader:
     tot_surf_2 = TotalSurface("trap1", 7)
     tot_vol_1 = TotalVolume("trap2", 5)
     min_vol_1 = MinimumVolume("retention", 2)
+    min_vol_2 = MinimumVolume("T", 2)
     max_vol_1 = MaximumVolume("T", 2)
+    max_vol_2 = MaximumVolume("trap2", 2)
+    min_surface_1 = MinimumSurface("solute", 1)
+    min_surface_2 = MinimumSurface("T", 2)
+    max_surface_1 = MaximumSurface("solute", 8)
+    max_surface_2 = MaximumSurface("T", 9)
+    avg_surface_1 = AverageSurface("solute", 4)
+    avg_surface_2 = AverageSurface("T", 6)
+    point_1 = PointValue(field="retention", x=2)
+    point_2 = PointValue(field="T", x=9)
+    cyl_surface_flux_1 = SurfaceFluxCylindrical("solute", 2)
+    cyl_surface_flux_2 = SurfaceFluxCylindrical("T", 3)
+    sph_surface_flux_1 = SurfaceFluxSpherical("solute", 5)
+    sph_surface_flux_2 = SurfaceFluxSpherical("T", 6)
 
     def test_simple(self):
         """
@@ -76,6 +96,54 @@ def test_all_quantities(self):
         ]
         assert header == expected_header
 
+    def test_with_units_simple(self):
+        """Test with quantities that don't require mesh dimension for unit"""
+        my_derv_quant = DerivedQuantities(
+            [
+                self.average_vol_1,
+                self.average_vol_2,
+                self.min_vol_1,
+                self.min_vol_2,
+                self.max_vol_1,
+                self.max_vol_2,
+                self.min_surface_1,
+                self.min_surface_2,
+                self.max_surface_1,
+                self.max_surface_2,
+                self.avg_surface_1,
+                self.avg_surface_2,
+                self.point_1,
+                self.point_2,
+                self.cyl_surface_flux_1,
+                self.cyl_surface_flux_2,
+                self.sph_surface_flux_1,
+                self.sph_surface_flux_2,
+            ],
+            show_units=True,
+        )
+        header = my_derv_quant.make_header()
+        expected_header = ["t(s)"] + [
+            "Average solute volume 3 (H m-3)",
+            "Average T volume 4 (K)",
+            "Minimum retention volume 2 (H m-3)",
+            "Minimum T volume 2 (K)",
+            "Maximum T volume 2 (K)",
+            "Maximum trap2 volume 2 (H m-3)",
+            "Minimum solute surface 1 (H m-3)",
+            "Minimum T surface 2 (K)",
+            "Maximum solute surface 8 (H m-3)",
+            "Maximum T surface 9 (K)",
+            "Average solute surface 4 (H m-3)",
+            "Average T surface 6 (K)",
+            "retention value at [2] (H m-3)",
+            "T value at [9] (K)",
+            "solute flux surface 2 (H s-1)",
+            "Heat flux surface 3 (W)",
+            "solute flux surface 5 (H s-1)",
+            "Heat flux surface 6 (W)",
+        ]
+        assert header == expected_header
+
 
 class TestAssignMeasuresToQuantities:
     """
@@ -211,7 +279,11 @@ def test_simple(self):
 
         my_derv_quant.data = []
         my_derv_quant.compute(t)
-        assert my_derv_quant.data[0] == expected_data
+
+        # title created in compute() method and appended to data
+        # so test line 2 for first data entry
+
+        assert my_derv_quant.data[1] == expected_data
 
     def test_two_quantities(self):
         """Check for the case of two festim.DerivedQuantity objects"""
@@ -232,7 +304,10 @@ def test_two_quantities(self):
         my_derv_quant.data = []
         my_derv_quant.compute(t)
 
-        assert my_derv_quant.data[0] == expected_data
+        # title created in compute() method and appended to data
+        # so test line 2 for first data entry
+
+        assert my_derv_quant.data[1] == expected_data
 
     def test_all_quantities(self):
         """Check for the case of many festim.DerivedQuantity objects"""
@@ -262,7 +337,10 @@ def test_all_quantities(self):
         my_derv_quant.data = []
         my_derv_quant.compute(t)
 
-        assert my_derv_quant.data[0] == expected_data
+        # title created in compute() method and appended to data
+        # so test line 2 for first data entry
+
+        assert my_derv_quant.data[1] == expected_data
 
 
 class TestWrite:
diff --git a/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py b/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py
index 1b1b71181..96aa175e4 100644
--- a/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py
+++ b/test/unit/test_exports/test_derived_quantities/test_hydrogen_flux.py
@@ -1,4 +1,6 @@
 from festim import HydrogenFlux
+from .tools import c_1D, c_2D, c_3D
+import pytest
 
 
 def test_field_is_solute():
@@ -9,3 +11,24 @@ def test_field_is_solute():
 
     my_flux = HydrogenFlux(2)
     assert my_flux.field == "solute"
+
+
+@pytest.mark.parametrize(
+    "function, expected_title",
+    [
+        (c_1D, "solute flux surface 3 (H m-2 s-1)"),
+        (c_2D, "solute flux surface 3 (H m-1 s-1)"),
+        (c_3D, "solute flux surface 3 (H s-1)"),
+    ],
+)
+def test_title_with_units(function, expected_title):
+    my_flux = HydrogenFlux(3)
+    my_flux.function = function
+    my_flux.show_units = True
+
+    assert my_flux.title == expected_title
+
+
+def test_title_without_units():
+    my_flux = HydrogenFlux(4)
+    assert my_flux.title == "Flux surface 4: solute"
diff --git a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
index 55446254e..b6d3323be 100644
--- a/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
+++ b/test/unit/test_exports/test_derived_quantities/test_surface_flux.py
@@ -3,6 +3,8 @@
 import math
 import numpy as np
 import pytest
+from festim import SurfaceFlux, k_B
+from .tools import c_1D, c_2D, c_3D
 
 
 @pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
@@ -264,3 +266,54 @@ def test_soret_raises_error_spherical():
     my_flux = SurfaceFluxSpherical("T", 1)
     with pytest.raises(NotImplementedError):
         my_flux.compute(soret=True)
+
+
+@pytest.mark.parametrize(
+    "function, field, expected_title",
+    [
+        (c_1D, "solute", "solute flux surface 3 (H m-2 s-1)"),
+        (c_1D, "T", "Heat flux surface 3 (W m-2)"),
+        (c_2D, "solute", "solute flux surface 3 (H m-1 s-1)"),
+        (c_2D, "T", "Heat flux surface 3 (W m-1)"),
+        (c_3D, "solute", "solute flux surface 3 (H s-1)"),
+        (c_3D, "T", "Heat flux surface 3 (W)"),
+    ],
+)
+def test_title_with_units(function, field, expected_title):
+    my_flux = SurfaceFlux(field=field, surface=3)
+    my_flux.function = function
+    my_flux.show_units = True
+
+    assert my_flux.title == expected_title
+
+
+def test_cylindrical_flux_title_no_units_solute():
+    """A simple test to check that the title is set correctly in
+    festim.CylindricalSurfaceFlux with a solute field without units"""
+
+    my_h_flux = SurfaceFluxCylindrical("solute", 2)
+    assert my_h_flux.title == "solute flux surface 2"
+
+
+def test_cylindrical_flux_title_no_units_temperature():
+    """A simple test to check that the title is set correctly in
+    festim.CylindricalSurfaceFlux with a T field without units"""
+
+    my_heat_flux = SurfaceFluxCylindrical("T", 4)
+    assert my_heat_flux.title == "Heat flux surface 4"
+
+
+def test_spherical_flux_title_no_units_solute():
+    """A simple test to check that the title is set correctly in
+    festim.SphericalSurfaceFlux with a solute field without units"""
+
+    my_h_flux = SurfaceFluxSpherical("solute", 3)
+    assert my_h_flux.title == "solute flux surface 3"
+
+
+def test_spherical_flux_title_no_units_temperature():
+    """A simple test to check that the title is set correctly in
+    festim.CSphericalSurfaceFlux with a T field without units"""
+
+    my_heat_flux = SurfaceFluxSpherical("T", 5)
+    assert my_heat_flux.title == "Heat flux surface 5"
diff --git a/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py b/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py
index fac0e5f46..92776a1b9 100644
--- a/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py
+++ b/test/unit/test_exports/test_derived_quantities/test_thermal_flux.py
@@ -1,4 +1,6 @@
 from festim import ThermalFlux
+import pytest
+from .tools import c_1D, c_2D, c_3D
 
 
 def test_field_is_T():
@@ -8,3 +10,25 @@ def test_field_is_T():
     """
     my_flux = ThermalFlux(2)
     assert my_flux.field == "T"
+
+
+@pytest.mark.parametrize(
+    "function, expected_title",
+    [
+        (c_1D, "Heat flux surface 5 (W m-2)"),
+        (c_2D, "Heat flux surface 5 (W m-1)"),
+        (c_3D, "Heat flux surface 5 (W)"),
+    ],
+)
+def test_title_with_units(function, expected_title):
+    my_flux = ThermalFlux(5)
+    my_flux.function = function
+    my_flux.show_units = True
+
+    assert my_flux.title == expected_title
+
+
+def test_title_without_units():
+    my_flux = ThermalFlux(5)
+
+    assert my_flux.title == "Flux surface 5: T"
diff --git a/test/unit/test_exports/test_derived_quantities/test_total_surface.py b/test/unit/test_exports/test_derived_quantities/test_total_surface.py
index f5e152595..ea0cd8354 100644
--- a/test/unit/test_exports/test_derived_quantities/test_total_surface.py
+++ b/test/unit/test_exports/test_derived_quantities/test_total_surface.py
@@ -2,6 +2,8 @@
 from festim import TotalSurface
 import fenics as f
 import pytest
+from .tools import c_1D, c_2D, c_3D
+import pytest
 
 
 @pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
@@ -43,3 +45,22 @@ def test_minimum(self):
 
         produced = self.my_total.compute()
         assert produced == expected
+
+
+@pytest.mark.parametrize(
+    "function, field, expected_title",
+    [
+        (c_1D, "solute", "Total solute surface 8 (H m-2)"),
+        (c_1D, "T", "Total T surface 8 (K)"),
+        (c_2D, "solute", "Total solute surface 8 (H m-1)"),
+        (c_2D, "T", "Total T surface 8 (K m)"),
+        (c_3D, "solute", "Total solute surface 8 (H)"),
+        (c_3D, "T", "Total T surface 8 (K m2)"),
+    ],
+)
+def test_title_with_units(function, field, expected_title):
+    my_export = TotalSurface(surface=8, field=field)
+    my_export.function = function
+    my_export.show_units = True
+
+    assert my_export.title == expected_title
diff --git a/test/unit/test_exports/test_derived_quantities/test_total_volume.py b/test/unit/test_exports/test_derived_quantities/test_total_volume.py
index 8f2e3cfa5..5fd3a8cb3 100644
--- a/test/unit/test_exports/test_derived_quantities/test_total_volume.py
+++ b/test/unit/test_exports/test_derived_quantities/test_total_volume.py
@@ -1,6 +1,8 @@
 from festim import TotalVolume
 import fenics as f
 import pytest
+from .tools import c_1D, c_2D, c_3D
+import pytest
 
 
 @pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)])
@@ -42,3 +44,22 @@ def test_minimum(self):
 
         produced = self.my_total.compute()
         assert produced == expected
+
+
+@pytest.mark.parametrize(
+    "function, field, expected_title",
+    [
+        (c_1D, "solute", "Total solute volume 2 (H m-2)"),
+        (c_1D, "T", "Total T volume 2 (K m)"),
+        (c_2D, "solute", "Total solute volume 2 (H m-1)"),
+        (c_2D, "T", "Total T volume 2 (K m2)"),
+        (c_3D, "solute", "Total solute volume 2 (H)"),
+        (c_3D, "T", "Total T volume 2 (K m3)"),
+    ],
+)
+def test_title_with_units(function, field, expected_title):
+    my_export = TotalVolume(volume=2, field=field)
+    my_export.function = function
+    my_export.show_units = True
+
+    assert my_export.title == expected_title
diff --git a/test/unit/test_exports/test_derived_quantities/tools.py b/test/unit/test_exports/test_derived_quantities/tools.py
new file mode 100644
index 000000000..8f9ce3089
--- /dev/null
+++ b/test/unit/test_exports/test_derived_quantities/tools.py
@@ -0,0 +1,13 @@
+import fenics as f
+
+mesh_1D = f.UnitIntervalMesh(10)
+V_1D = f.FunctionSpace(mesh_1D, "P", 1)
+c_1D = f.interpolate(f.Expression("x[0]", degree=1), V_1D)
+
+mesh_2D = f.UnitSquareMesh(10, 10)
+V_2D = f.FunctionSpace(mesh_2D, "P", 1)
+c_2D = f.interpolate(f.Expression("x[0]", degree=1), V_2D)
+
+mesh_3D = f.UnitCubeMesh(10, 10, 10)
+V_3D = f.FunctionSpace(mesh_3D, "P", 1)
+c_3D = f.interpolate(f.Expression("x[0]", degree=1), V_3D)