diff --git a/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py b/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
index bc46a4af66365fda14a35761f40fde6854447764..b08f0f71e91e7cda5c1c13fc94e679dc42b4d95a 100644
--- a/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
+++ b/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
@@ -19,7 +19,6 @@ from ogstools.meshlib.region import (
 )
 
 # See other examples for different meshing algorithms
-from ogstools.propertylib import ScalarProperty  # For visualization only
 
 # %%
 # Define a simple surface
@@ -36,7 +35,7 @@ mesh = to_region_tetraeder(ls, 40).mesh
 # Visualize the prism mesh
 
 slices = np.reshape(list(mesh.slice_along_axis(n=4, axis="y")), (2, 2))
-fig = mpl.plot(slices, property=ScalarProperty("MaterialIDs"))
+fig = mpl.plot(slices, "MaterialIDs")
 
 
 # %%
diff --git a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
index 263003afa3edd71b90652b33b64e1c081eb301cb..be1b6bc9ba047c02f16fb9570aba8c1d3b805a6e 100644
--- a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
+++ b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
@@ -23,9 +23,7 @@ from ogstools.meshlib.region import (
     to_region_tetraeder,
     to_region_voxel,
 )
-from ogstools.propertylib import ScalarProperty  # For visualization only
 
-material_id = ScalarProperty("MaterialIDs")
 mpl.setup.reset()
 mpl.setup.length.output_unit = "km"
 mpl.setup.figsize = [24, 12]
@@ -62,7 +60,7 @@ tm = to_region_tetraeder(layer_set1, resolution=200).mesh
 
 mesh = pm
 slices = np.reshape(list(mesh.slice_along_axis(n=4, axis="y")), (2, 2))
-fig = mpl.plot(slices, property=material_id)
+fig = mpl.plot(slices, "MaterialIDs")
 for ax, slice in zip(fig.axes, np.ravel(slices)):
     ax.set_title(f"z = {slice.center[2]:.1f} {mpl.setup.length.data_unit}")
 
@@ -78,6 +76,6 @@ names = [
 ]
 
 x_slices = np.reshape([mesh.slice("x") for mesh in meshes], (-1, 1))
-fig = mpl.plot(x_slices, property=material_id)
+fig = mpl.plot(x_slices, "MaterialIDs")
 for ax, name in zip(fig.axes, names):
     ax.set_title(name)
diff --git a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
index 3594b169cacc07895a208f6256f291bd1ff991c8..787e44326db7adfebd2df63a884842d270f789d1 100644
--- a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
+++ b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
@@ -19,14 +19,14 @@ mpl.setup.reset()
 mpl.setup.length.output_unit = "km"
 mpl.setup.material_names = {i + 1: f"Layer {i+1}" for i in range(26)}
 mesh = meshseries_THM_2D.read(1)
-fig = mpl.plot(mesh, property=THM.material_id)
+fig = mpl.plot(mesh, THM.material_id)
 
 # %%
 # Now, let's plot the temperature field (point_data) at the first timestep.
 # The default temperature property from the `propertylib` reads the temperature
 # data as Kelvin and converts them to degrees Celsius.
 
-fig = mpl.plot(mesh, property=THM.temperature)
+fig = mpl.plot(mesh, THM.temperature)
 
 # %%
 # This example has hydraulically deactivated subdomains:
diff --git a/tests/test_meshplotlib.py b/tests/test_meshplotlib.py
index 2326ca94c8a8c5e6b3c2cbe39d91710e798e5a09..b3e79df79b32f573032bc1c2ea5cabc2686e92d0 100644
--- a/tests/test_meshplotlib.py
+++ b/tests/test_meshplotlib.py
@@ -53,8 +53,8 @@ class MeshplotlibTest(unittest.TestCase):
         setup.material_names = {i + 1: f"Layer {i+1}" for i in range(26)}
         meshseries = examples.meshseries_THM_2D
         mesh = meshseries.read(1)
-        plot(mesh, property=THM.material_id)
-        plot(mesh, property=THM.temperature)
+        plot(mesh, THM.material_id)
+        plot(mesh, THM.temperature)
         plot(mesh, ScalarProperty("pressure_active"))
         plot(mesh.threshold((1, 3), "MaterialIDs"), THM.velocity)
         fig = plot(mesh, THM.displacement[0])
@@ -102,7 +102,7 @@ class MeshplotlibTest(unittest.TestCase):
             f"{THIS_DIR}/data/meshplotlib/"
             "2D_single_fracture_HT_2D_single_fracture.xdmf"
         ).read(0)
-        plot(mesh, property=THM.temperature)
+        plot(mesh, THM.temperature)
 
 
 if __name__ == "__main__":