diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index ff024fd969b6d6c85b69bca8c0569c628997c2b2..1f0e390d385ccf1883110932731304de569b3892 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -328,8 +328,7 @@ def get_combined_levels(
 def resolve_property(property: Union[Property, str]) -> Property:
     if isinstance(property, Property):
         return property
-    prop = THM.find_property(property)
-    return prop if prop else Property(property)
+    return THM.find_property(property)
 
 
 def _plot_on_figure(
diff --git a/ogstools/propertylib/property_collection.py b/ogstools/propertylib/property_collection.py
index 3a9fa96114ef514cc43bc023bb9b08dcab35dadf..62d013d0276922fa4433fc1db238f2bb5d68d816 100644
--- a/ogstools/propertylib/property_collection.py
+++ b/ogstools/propertylib/property_collection.py
@@ -43,7 +43,7 @@ class PropertyCollection:
 
     def find_property(
         self, output_name: str, dim: Opt[Literal[2, 3]] = None
-    ) -> Opt[Property]:
+    ) -> Property:
         """Return predefined property with given output_name."""
         for prop in self.get_properties(dim):
             if prop.output_name == output_name:
@@ -52,4 +52,4 @@ class PropertyCollection:
         for prop in self.get_properties(dim):
             if prop.data_name == output_name:
                 return prop
-        return None
+        return Property(output_name)
diff --git a/ogstools/studies/convergence/convergence_study.py b/ogstools/studies/convergence/convergence_study.py
new file mode 100644
index 0000000000000000000000000000000000000000..d1df704adfb4690195c9c681d32860a5cad35bb9
--- /dev/null
+++ b/ogstools/studies/convergence/convergence_study.py
@@ -0,0 +1,92 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.15.2
+#   kernelspec:
+#     display_name: .venv
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+# # Convergence study
+#
+# This script performs a convergence study and generates plots to analyze the
+# convergence of numerical simulations.
+# Below are the custom parameters for this report.
+
+# %% tags=["parameters"]
+mesh_paths: "pv.DataSet" = None
+property_name: str = ""
+
+# %% [markdown]
+# Import required modules and setup plots.
+
+# %%
+import pyvista as pv  # noqa: E402
+
+import ogstools.meshplotlib as mpl  # noqa: E402
+from ogstools.propertylib import THM  # noqa: E402
+from ogstools.studies.convergence import (  # noqa: E402
+    convergence_metrics,
+    plot_convergence,
+    plot_convergence_errors,
+    richardson_extrapolation,
+)
+
+mpl.setup.reset()
+mpl.setup.show_element_edges = True
+mpl.setup.ax_aspect_ratio = 1
+
+# %% [markdown]
+# Read the meshes, get Property object from property name, define topology and calculate Richardson extrapolation.
+
+# %%
+meshes = [pv.read(mesh_path) for mesh_path in mesh_paths]
+mesh_property = THM.find_property(property_name)
+mesh_property.output_unit = ""
+mesh_property.data_unit = ""
+topology = meshes[-3]
+richardson = richardson_extrapolation(meshes, mesh_property, topology)
+
+# %% [markdown]
+# Plotting the grid convergence.
+
+# %%
+fig = mpl.plot(richardson, "grid_convergence")
+
+# %% [markdown]
+# Plotting the 3 finest discretizations.
+
+# %%
+fig = mpl.plot(meshes[-3:], mesh_property)
+
+# %% [markdown]
+# Plotting the Richardson extrapolation.
+
+# %%
+fig = mpl.plot(richardson, mesh_property)
+
+# %% [markdown]
+# Table of convergence metrics.
+
+# %%
+mpl.core.plt.rcdefaults()
+metrics = convergence_metrics(meshes, richardson, mesh_property)
+metrics.style.format("{:,.4g}").hide()
+
+# %% [markdown]
+# Absolute convergence metrics.
+
+# %%
+fig = plot_convergence(metrics, mesh_property)
+
+# %% [markdown]
+# Relative errors in loglog-scale.
+
+# %%
+fig = plot_convergence_errors(metrics)
diff --git a/ogstools/studies/convergence/generate_report.py b/ogstools/studies/convergence/generate_report.py
index 01efea0049214bba72c8b3942f397b2f58177a9f..7ca772f5807e373cec7b62a39a9c9b90c1f5dd24 100644
--- a/ogstools/studies/convergence/generate_report.py
+++ b/ogstools/studies/convergence/generate_report.py
@@ -2,16 +2,21 @@
 
 from pathlib import Path
 
+import jupytext
 import papermill as pm
 
 from ogstools.studies.convergence.examples import mesh_paths
 
 parent = Path(__file__).resolve().parent
 
+nb = jupytext.read("convergence_study.py")
+# TODO: make temporary
+jupytext.write(nb, "convergence_study.ipynb", fmt="py:percent")
+
 property_name = "pressure"
 params = {"mesh_paths": mesh_paths, "property_name": property_name}
 pm.execute_notebook(
-    input_path=str(parent) + "/template.ipynb",
+    input_path=str(parent) + "/convergence_study.ipynb",
     output_path="convergence_study_" + property_name + ".ipynb",
     parameters=params,
 )
diff --git a/pyproject.toml b/pyproject.toml
index d80f375ef3319d15715e91c538b531c8ef2dac21..6af386e61115ff9048657d3dba72660a4c8e8bd9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,6 +12,7 @@ readme = "README.md"
 requires-python = '>=3.9'
 dependencies = [
   "h5py>=3.9.0",
+  "jupytext>=1.15.2",
   "meshio==5.3.4",
   "matplotlib>=3.7.1",
   "pyvista>=0.39.1",
@@ -20,6 +21,7 @@ dependencies = [
   "Pint>=0.22",
   "Pillow>=9.5.0",
   "pandas>=2.0.3",
+  "papermill>=2.4.0",
   "ogs>=6.4.5.dev0",
 ]