diff --git a/docs/_static/examples/meshlib/plot_calculate_diff_thumbnail.png b/docs/_static/examples/meshlib/plot_calculate_diff_thumbnail.png
new file mode 100644
index 0000000000000000000000000000000000000000..b87fba84518db946a24ba9c3cc30202187ddd68c
Binary files /dev/null and b/docs/_static/examples/meshlib/plot_calculate_diff_thumbnail.png differ
diff --git a/docs/examples/howto_meshlib/plot_calculate_diff.py b/docs/examples/howto_meshlib/plot_calculate_diff.py
new file mode 100644
index 0000000000000000000000000000000000000000..ea6828acf646ca0162e1450072415062d89256b2
--- /dev/null
+++ b/docs/examples/howto_meshlib/plot_calculate_diff.py
@@ -0,0 +1,159 @@
+"""
+Differences between meshes
+==========================
+
+.. sectionauthor:: Feliks Kiszkurno (Helmholtz Centre for Environmental Research GmbH - UFZ)
+
+This example explains how to use functions from meshlib to calculate differences between meshes.
+"""
+
+# %%
+
+# sphinx_gallery_start_ignore
+
+# sphinx_gallery_thumbnail_path = '_static/examples/meshlib/plot_calculate_diff_thumbnail.png'
+
+# sphinx_gallery_end_ignore
+
+from ogstools.meshlib import difference, difference_matrix, difference_pairwise
+from ogstools.meshplotlib.examples import meshseries_THM_2D
+from ogstools.propertylib import presets
+
+mesh_property = presets.temperature
+
+# %%
+# 0. Introduction
+# ---------------
+# This example shows how to calculate differences between meshes.
+# It is possible to calculate the difference between multiple meshes at the same time. Multiple meshes can be provided either as list or Numpy arrays.
+# 4 ways of calculating the difference are presented here.
+
+# %%
+# 1. Difference between two meshes
+# --------------------------------
+# The simplest case is calculating the difference between two meshes. For this example, we read two different timesteps from a MeshSeries. It is not required that they belong to the same MeshSeries object. As long, as the meshes share the same topology and contain the mesh_property of interest, the difference will work fine.
+mesh1 = meshseries_THM_2D.read(0)
+mesh2 = meshseries_THM_2D.read(-1)
+
+# %% [markdown]
+# The following call will return a mesh containing the difference of the mesh_property between the two provided meshes:
+#
+# .. math::
+#
+#   \text{Mesh}_1 - \text{Mesh}_2
+#
+
+mesh_diff = difference(mesh_property, mesh1, mesh2)
+
+# %% [markdown]
+# This returned object will be a PyVista UnstructuredGrid object:
+
+# %%
+print(f"Type of mesh_diff: {type(mesh_diff)}")
+
+# %% [markdown]
+# 2. Pairwise difference
+# ----------------------
+# In order to calculate pairwise differences, two lists or arrays of equal length have to be provided as the input. They can contain an arbitrary number of different meshes, as long as their length is equal.
+#
+# Consider the two following lists:
+#
+# .. math::
+#
+#   \text{List}_{1} = \begin{bmatrix} A_{1} & B_{1} & C_{1}\\ \end{bmatrix}
+#
+# and
+#
+# .. math::
+#
+#   \text{List}_{2} = \begin{bmatrix} A_{2} & B_{2} & C_{2}\\ \end{bmatrix}
+#
+# The output array will contain following differences between meshes:
+#
+# .. math::
+#
+#   \begin{bmatrix} A_{1}-A_{2} & B_{1}-B_{2} & C_{1}-C_{2}\\ \end{bmatrix}
+#
+# and will have the same shape as the input lists. As in the example below:
+
+meshes_1 = [mesh1] * 3
+meshes_2 = [mesh2] * 3
+
+mesh_diff_pair_wise = difference_pairwise(mesh_property, meshes_1, meshes_2)
+
+# %%
+print(f"Length of mesh_list1: {len(meshes_1)}")
+
+# %%
+print(f"Length of mesh_list2: {len(meshes_2)}")
+
+# %%
+print(f"Shape of mesh_diff_pair_wise: {mesh_diff_pair_wise.shape}")
+
+# %% [markdown]
+# 3. Matrix difference - one array
+# --------------------------------
+# If only one list or array is provided, the differences between every possible pair of objects in the array will be returned.
+#
+# Consider following list:
+#
+# .. math::
+#
+#   \text{List} = \begin{bmatrix} A & B & C\\ \end{bmatrix}
+#
+# The output array will contain following differences between meshes:
+#
+# .. math::
+#
+#   \begin{bmatrix} A-A & B-A & C-A\\ A-B & B-B & C-B \\ A-C & B-C & C-C \\ \end{bmatrix}
+#
+# and will have shape of (len(mesh_list), len(mesh_list)). As in the following example:
+
+mesh_list = [mesh1, mesh2, mesh1, mesh2]
+
+mesh_diff_matrix = difference_matrix(mesh_property, mesh_list)
+
+# %%
+print(f"Length of mesh_list1: {len(mesh_list)}")
+
+# %%
+print(f"Shape of mesh_list1: {mesh_diff_matrix.shape}")
+
+# %% [markdown]
+# 4. Matrix difference - two arrays of different length
+# -----------------------------------------------------
+# Unlike difference_pairwise(), difference_matrix() can accept two lists/arrays of different lengths. As in Section 3, the differences between all possible pairs of meshes in both lists/arrays will be calculated.
+#
+# Consider following lists:
+#
+# .. math::
+#
+#   \text{List}_1 = \begin{bmatrix} A_1 & B_1 & C_1\\ \end{bmatrix}
+#
+# .. math::
+#
+#   \text{List}_2 = \begin{bmatrix} A_2 & B_2 \\ \end{bmatrix}
+#
+# The output array will contain following differences between meshes:
+#
+# .. math::
+#
+#   \begin{bmatrix} A_1-A_2 & A_1-B_2 \\ B_1-A_2 & B_1-B_2 \\ C_1-A_2 & C_1-B_2 \\ \end{bmatrix}
+#
+# and will have a shape of (len(mesh_list_matrix_1), len(mesh_list_matrix_2)). As in the following example:
+
+mesh_list_matrix_1 = [mesh1, mesh2, mesh1]
+mesh_list_matrix_2 = [mesh2, mesh1]
+
+mesh_diff_matrix = difference_matrix(
+    mesh_property, mesh_list_matrix_1, mesh_list_matrix_2
+)
+
+# %%
+print(f"Length of mesh_list_matrix_1: {len(mesh_list_matrix_1)}")
+
+# %%
+print(f"Length of mesh_list_matrix_2: {len(mesh_list_matrix_2)}")
+
+# %%
+print(f"Shape of mesh_diff_matrix: {mesh_diff_matrix.shape}")
diff --git a/docs/examples/howto_meshplotlib/plot_special.py b/docs/examples/howto_meshplotlib/plot_special.py
index 9de4cc7e1378180be0a7350dbe154d67eeae8d25..48ae7d864291993e0355e838ad01f85bb992117c 100644
--- a/docs/examples/howto_meshplotlib/plot_special.py
+++ b/docs/examples/howto_meshplotlib/plot_special.py
@@ -44,5 +44,5 @@ fig = plot_limit(mesh_series, si, "min")
 # Difference between the last and he first timestep:
 
 # %%
-diff_mesh = difference(mesh_series.read(-1), mesh_series.read(0), si)
+diff_mesh = difference(si, mesh_series.read(-1), mesh_series.read(0))
 fig = plot(diff_mesh, si.delta)
diff --git a/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py b/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
index 4e9030a901d1eeb5d8d2ea761dbe4942eb650060..f3e5b8a020c1f6f2ceb8c6dbbea22be44e9a7b53 100644
--- a/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
+++ b/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
@@ -50,7 +50,7 @@ ax[0].set_title(r"$T(\mathrm{t}_{0})$")
 plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
 ax[1].set_title(r"$T(\mathrm{t}_{end})$")
 diff_mesh = difference(
-    meshseries.read(0), meshseries.read(1), presets.temperature
+    presets.temperature, meshseries.read(1), meshseries.read(0)
 )
 plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
 ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
diff --git a/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py b/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
index ba05f7276df2d26e5162cbf2d36afceb4ccdca79..26773da7ac9d81b1781e4917dea35fde0ab39042 100644
--- a/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
+++ b/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
@@ -54,8 +54,8 @@ fig = plot([mesh_a, mesh_b], temperature)
 fig, ax = plt.subplots(2, 2, figsize=(40, 20), sharex=True, sharey=True)
 plot(mesh_a, temperature, fig=fig, ax=ax[0][0])
 plot(mesh_b, temperature, fig=fig, ax=ax[1][0])
-diff_ab = difference(mesh_b, mesh_a, temperature)
-diff_ba = difference(mesh_b, mesh_a, temperature)
+diff_ab = difference(temperature, mesh_b, mesh_a)
+diff_ba = difference(temperature, mesh_b, mesh_a)
 plot(diff_ab, temperature.delta, fig=fig, ax=ax[0][1])
 plot(diff_ba, temperature.delta, fig=fig, ax=ax[1][1])
 fig.tight_layout()
diff --git a/ogstools/meshlib/__init__.py b/ogstools/meshlib/__init__.py
index abcff3d357b6461ae7ab383e69863efb77870b3a..ba12045609a504fd7240f48250d93ddad7e8a385 100644
--- a/ogstools/meshlib/__init__.py
+++ b/ogstools/meshlib/__init__.py
@@ -1,7 +1,7 @@
 from .boundary import Boundary, Layer, LocationFrame, Raster
 from .boundary_set import LayerSet
 from .boundary_subset import Surface
-from .data_processing import difference
+from .data_processing import difference, difference_matrix, difference_pairwise
 from .gmsh_meshing import cuboid, rect
 from .mesh_series import MeshSeries
 
@@ -14,6 +14,8 @@ __all__ = [
     "MeshSeries",
     "LayerSet",
     "difference",
+    "difference_pairwise",
+    "difference_matrix",
     "rect",
     "cuboid",
 ]
diff --git a/ogstools/meshlib/data_processing.py b/ogstools/meshlib/data_processing.py
index b547566f2bce87aabebdd604caa3daa20f888a06..61f359d4ee389f19daf7e7a3bfae8020d4ecfcbb 100644
--- a/ogstools/meshlib/data_processing.py
+++ b/ogstools/meshlib/data_processing.py
@@ -1,24 +1,86 @@
-from typing import Union
+from itertools import product
+from typing import Optional, Union
 
+import numpy as np
 import pyvista as pv
 
 from ogstools.propertylib import Property, presets
 
 
 def difference(
+    mesh_property: Union[Property, str],
     mesh1: pv.UnstructuredGrid,
     mesh2: pv.UnstructuredGrid,
-    mesh_property: Union[Property, str],
 ) -> pv.UnstructuredGrid:
     """
     Compute the difference between two meshes based on a specified property.
 
+    :param mesh_property: The property to of interest.
     :param mesh1: The first mesh to be subtracted from.
     :param mesh2: The second mesh whose data is subtracted from the first mesh.
-    :param mesh_property: The property to of interest.
     :returns: A new mesh representing the difference between mesh1 and mesh2.
     """
     mesh_property = presets.get_preset(mesh_property, mesh1)
     diff_mesh = mesh1.copy(deep=True)
     diff_mesh[mesh_property.data_name] -= mesh2[mesh_property.data_name]
     return diff_mesh
+
+
+def difference_pairwise(
+    mesh_property: Union[Property, str],
+    meshes_1: Union[list, np.ndarray],
+    meshes_2: Union[list, np.ndarray],
+) -> np.ndarray:
+    """
+    Compute pairwise difference between meshes from two lists/arrays
+    (they have to be of the same length).
+
+    :param mesh_property: The property to of interest.
+    :param meshes_1: The first list/array of meshes to be subtracted from.
+    :param meshes_2: The second list/array of meshes whose data is subtracted
+                     from the first list/array of meshes - meshes_1.
+    :returns: An array of differences between meshes from meshes_1 and meshes_2.
+    """
+    meshes_1 = np.asarray(meshes_1).flatten()
+    meshes_2 = np.asarray(meshes_2).flatten()
+    if len(meshes_1) != len(meshes_2):
+        msg = "Mismatch in length of provided lists/arrays. \
+              Their length has to be identical to calculate pairwise \
+              difference. Did you intend to use difference_matrix()?"
+        raise RuntimeError(msg)
+    diff_mesh = [
+        difference(mesh_property, m1, m2) for m1, m2 in zip(meshes_1, meshes_2)
+    ]
+    return np.array(diff_mesh)
+
+
+def difference_matrix(
+    mesh_property: Union[Property, str],
+    meshes_1: Union[list, np.ndarray],
+    meshes_2: Optional[Union[list, np.ndarray]] = None,
+) -> np.ndarray:
+    """
+    Compute the difference between all combinations of two meshes
+    from one or two arrays based on a specified property.
+
+    :param mesh_property: The property to of interest.
+    :param meshes_1: The first list/array of meshes to be subtracted from.
+    :param meshes_2: The second list/array of meshes, it is subtracted from
+                     the first list/array of meshes - meshes_1 (optional).
+    :returns: An array of differences between meshes from meshes_1 and meshes_2.
+    """
+    if not isinstance(meshes_1, (list, np.ndarray)):
+        msg = "mesh1 is neither of type list nor np.ndarray"  # type: ignore[unreachable]
+        raise TypeError(msg)
+    meshes_1 = np.asarray(meshes_1).flatten()
+    if not isinstance(meshes_2, (list, np.ndarray)) and meshes_2 is not None:
+        msg = "mesh2 is neither of type list nor np.ndarray."  # type: ignore[unreachable]
+        raise TypeError(msg)
+    if meshes_2 is None:
+        meshes_2 = meshes_1.copy()
+    meshes_2 = np.asarray(meshes_2).flatten()
+    diff_mesh = [
+        difference(mesh_property, m1, m2)
+        for m1, m2 in product(meshes_1, meshes_2)
+    ]
+    return np.array(diff_mesh).reshape((len(meshes_1), len(meshes_2)))
diff --git a/ogstools/studies/templates/convergence_study.py b/ogstools/studies/templates/convergence_study.py
index 30ffe046c2d27ed88c20a8d0f6a706c82043a84b..4cee10d9e7dd16e2447ae19da203961609e78b86 100644
--- a/ogstools/studies/templates/convergence_study.py
+++ b/ogstools/studies/templates/convergence_study.py
@@ -77,7 +77,7 @@ fig = meshplotlib.plot(richardson, mesh_property)
 data_key = mesh_property.data_name
 if reference_solution_path is None:
     diff_mesh = meshlib.difference(
-        richardson, topology.sample(meshes[-1]), mesh_property
+        mesh_property, richardson, topology.sample(meshes[-1])
     )
     fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
 else:
@@ -85,7 +85,7 @@ else:
         meshlib.MeshSeries(reference_solution_path).read_closest(timevalue)
     )
     diff_mesh = meshlib.difference(
-        reference_solution, richardson, mesh_property
+        mesh_property, reference_solution, richardson
     )
     fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
 
diff --git a/tests/test_meshlib.py b/tests/test_meshlib.py
index ffc9e7f2944ac97ecfadb12623b09e2c622123ee..2f46c458aaf1d256f460fb34e0556a3fd9c16136 100644
--- a/tests/test_meshlib.py
+++ b/tests/test_meshlib.py
@@ -3,8 +3,17 @@
 import unittest
 
 import numpy as np
+from pyvista import UnstructuredGrid
 
-from ogstools.meshlib import MeshSeries, examples
+from ogstools.meshlib import (
+    MeshSeries,
+    difference,
+    difference_matrix,
+    difference_pairwise,
+    examples,
+)
+from ogstools.meshplotlib import examples as examples_mpl
+from ogstools.propertylib import presets
 
 
 class UtilsTest(unittest.TestCase):
@@ -45,3 +54,66 @@ class UtilsTest(unittest.TestCase):
         for method in ["nearest", "linear", None]:
             values = mesh_series.probe(points, "temperature", method)
             self.assertTrue(not np.any(np.isnan(values)))
+
+    def test_diff_two_meshes(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        mesh1 = meshseries.read(0)
+        mesh2 = meshseries.read(-1)
+        mesh_diff = difference(mesh_property, mesh1, mesh2)
+        self.assertTrue(isinstance(mesh_diff, UnstructuredGrid))
+
+    def test_diff_pairwise(self):
+        n = 5
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = [meshseries.read(0)] * n
+        meshes2 = [meshseries.read(-1)] * n
+        meshes_diff = difference_pairwise(mesh_property, meshes1, meshes2)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray) and len(meshes_diff) == n
+        )
+
+    def test_diff_matrix_single_list(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = [meshseries.read(0), meshseries.read(-1)]
+        meshes_diff = difference_matrix(mesh_property, meshes1)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes1))
+        )
+
+    def test_diff_matrix_single_numpy(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = np.array([meshseries.read(0), meshseries.read(-1)])
+        meshes_diff = difference_matrix(mesh_property, meshes1)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes1))
+        )
+
+    def test_diff_matrix_unequal_list(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = [meshseries.read(0), meshseries.read(-1)]
+        meshes2 = [meshseries.read(0), meshseries.read(-1), meshseries.read(-1)]
+        meshes_diff = difference_matrix(mesh_property, meshes1, meshes2)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes2))
+        )
+
+    def test_diff_matrix_unequal_numpy(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = np.array([meshseries.read(0), meshseries.read(-1)])
+        meshes2 = np.array(
+            [meshseries.read(0), meshseries.read(-1), meshseries.read(-1)]
+        )
+        meshes_diff = difference_matrix(mesh_property, meshes1, meshes2)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes2))
+        )
diff --git a/tests/test_meshplotlib.py b/tests/test_meshplotlib.py
index bde3133f30d73190058d88313ce7907bcd5a996c..98f7271bf8ac1e791df0f29b5717b3de2e4f602a 100644
--- a/tests/test_meshplotlib.py
+++ b/tests/test_meshplotlib.py
@@ -135,7 +135,7 @@ class MeshplotlibTest(unittest.TestCase):
     def test_diff_plots(self):
         """Test creation of difference plots."""
         meshseries = examples.meshseries_CT_2D
-        plot(difference(meshseries.read(0), meshseries.read(1), "Si"), "Si")
+        plot(difference("Si", meshseries.read(0), meshseries.read(1)), "Si")
 
     def test_user_defined_ax(self):
         """Test creating plot with subfigures and user provided ax"""
@@ -146,7 +146,7 @@ class MeshplotlibTest(unittest.TestCase):
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
         ax[1].set_title(r"$T(\mathrm{t}_{end})$")
         diff_mesh = difference(
-            meshseries.read(0), meshseries.read(1), presets.temperature
+            presets.temperature, meshseries.read(0), meshseries.read(1)
         )
         plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
         ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
@@ -197,8 +197,8 @@ class MeshplotlibTest(unittest.TestCase):
         ax = ax.flatten()
         plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0])
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
-        diff_ab = difference(mesh_a, mesh_b, presets.temperature)
-        diff_ba = difference(mesh_b, mesh_a, presets.temperature)
+        diff_ab = difference(presets.temperature, mesh_a, mesh_b)
+        diff_ba = difference(presets.temperature, mesh_b, mesh_a)
         plot(diff_ab, presets.temperature.delta, fig=fig, ax=ax[2])
         plot(diff_ba, presets.temperature.delta, fig=fig, ax=ax[3])
         fig.tight_layout()
@@ -211,8 +211,8 @@ class MeshplotlibTest(unittest.TestCase):
         fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)
         plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0][0])
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1][0])
-        diff_ab = difference(mesh_a, mesh_b, presets.temperature)
-        diff_ba = difference(mesh_b, mesh_a, presets.temperature)
+        diff_ab = difference(presets.temperature, mesh_a, mesh_b)
+        diff_ba = difference(presets.temperature, mesh_b, mesh_a)
         plot(diff_ab, presets.temperature.delta, fig=fig, ax=ax[0][1])
         plot(diff_ba, presets.temperature.delta, fig=fig, ax=ax[1][1])
         ax = label_spatial_axes(ax, np.array([0, 1]))