diff --git a/docs/images/depth_optimization.png b/docs/images/depth_optimization.png
index 9178f34d489c2bb1e52ad5d3e8876c543c7dcda9..1a5e88e1e03c7bb9b72a186760ab8f531b10a6fd 100644
Binary files a/docs/images/depth_optimization.png and b/docs/images/depth_optimization.png differ
diff --git a/docs/images/exceedance_probabilities.png b/docs/images/exceedance_probabilities.png
index 09b4b7169022dc34948d3572c666cf10dbbfd00b..6fd558433eee0e2d7eda07377bd90e178d422d0c 100644
Binary files a/docs/images/exceedance_probabilities.png and b/docs/images/exceedance_probabilities.png differ
diff --git a/docs/images/input_maps_1.png b/docs/images/input_maps_1.png
new file mode 100644
index 0000000000000000000000000000000000000000..d7db153913e1bcc0b108f11d5f1da1d7e6772695
Binary files /dev/null and b/docs/images/input_maps_1.png differ
diff --git a/docs/images/input_maps_2.png b/docs/images/input_maps_2.png
new file mode 100644
index 0000000000000000000000000000000000000000..e03cb9855dc425fa1791a028c1431d69bfe6a624
Binary files /dev/null and b/docs/images/input_maps_2.png differ
diff --git a/docs/images/kh.png b/docs/images/kh.png
index 7bf38d058240ef8a6aab2bd4d6ae3c386658116b..5171a21188d2e332848d29e66c99e505001777c6 100644
Binary files a/docs/images/kh.png and b/docs/images/kh.png differ
diff --git a/docs/images/maps.png b/docs/images/maps.png
index ec227cd134019ad10afe44b31c385f15418ccac5..84dc50736af2be840f45112ed99617ae26e6080a 100644
Binary files a/docs/images/maps.png and b/docs/images/maps.png differ
diff --git a/docs/images/npv.png b/docs/images/npv.png
index 70fb00fddf0f88c890c907ff56c4045f14f13994..f6c82b203e5e462eb67198ded2a733a9978c1cfd 100644
Binary files a/docs/images/npv.png and b/docs/images/npv.png differ
diff --git a/docs/images/pos.png b/docs/images/pos.png
index 4cc03539264ff7d707f2f6d276ec063061ceeb61..97ff98d1e9d4d410d039837d7454ccb66940f62c 100644
Binary files a/docs/images/pos.png and b/docs/images/pos.png differ
diff --git a/docs/images/reservoir_distributions.png b/docs/images/reservoir_distributions.png
index cee50f34bc1f333529bf49d61cbf47b1cc526ec5..5f9d802e79ebf7be9980f9af77a08ee625e08523 100644
Binary files a/docs/images/reservoir_distributions.png and b/docs/images/reservoir_distributions.png differ
diff --git a/docs/images/simulation_results.png b/docs/images/simulation_results.png
index 6ac25e63ad56ed1e3625baa873346ac7610c9d10..b2369708bf8e3fc4733cd8076b8184f7a9ceedb8 100644
Binary files a/docs/images/simulation_results.png and b/docs/images/simulation_results.png differ
diff --git a/docs/usage/customised_stochastic_simulations.md b/docs/usage/customised_stochastic_simulations.md
index f788fdd30b4278f4e3878512726ca12d3ac94c59..e6b994e4a93b7e68083b0c8d82a0d677a2c4dd82 100644
--- a/docs/usage/customised_stochastic_simulations.md
+++ b/docs/usage/customised_stochastic_simulations.md
@@ -26,7 +26,7 @@ output_data_path.mkdir(parents=True, exist_ok=True)
 # generate simulation samples across desired reservoir properties
 Nsamples = 1000
 thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples)
-porosity_samples = np.random.uniform(low=0.5, high=0.8, size=Nsamples)
+porosity_samples = np.random.uniform(low=0.2, high=0.3, size=Nsamples)
 ntg_samples = np.random.uniform(low=0.25, high=0.5, size=Nsamples)
 depth_samples = np.random.normal(loc=2000, scale=100, size=Nsamples)
 permeability_samples = np.random.normal(loc=400, scale=100, size=Nsamples)
diff --git a/docs/usage/customized_props.md b/docs/usage/customized_props.md
index 7f28563d71e0174df76b556a4a6c595b63ad2f9d..0496e23b56d002c781b1b840f08b63822ed7f74b 100644
--- a/docs/usage/customized_props.md
+++ b/docs/usage/customized_props.md
@@ -45,7 +45,7 @@ import xarray as xr
 input_data = xr.Dataset({
     "thickness": 300,
     "ntg": 0.5,
-    "porosity": 0.5,
+    "porosity": 0.2,
     "depth": 2000,
     "permeability": 300,
 })
@@ -72,7 +72,7 @@ input_data = xr.Dataset({
     "thickness_mean": 300,
     "thickness_sd": 50,
     "ntg": 0.5,
-    "porosity": 0.5,
+    "porosity": 0.2,
     "depth": 5000,
     "ln_permeability_mean": 5,
     "ln_permeability_sd": 0.5,
diff --git a/docs/usage/depth_optimization.md b/docs/usage/depth_optimization.md
index 6fb94bf7140cb69f8a3a127a96a1d4b259940c3e..30105684632791c26942512cf1199fc585f54328 100644
--- a/docs/usage/depth_optimization.md
+++ b/docs/usage/depth_optimization.md
@@ -28,34 +28,43 @@ from matplotlib import pyplot as plt
 from pythermogis import calculate_doublet_performance
 
 # output directory to write output to
-output_data_path = Path(__file__).parent / "resources" / "test_output" / "example9"
+output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example9"
 output_data_path.mkdir(parents=True, exist_ok=True)
 
 # define the poro-depth and the poro-perm relationships:
-# porosity = a * exp(-b * depth), with depth in km
-poro_a = 28.97
-poro_b = 0.21
-poro_std = 1  # 1% uncertainty, if this goes larger it blows up massively...
-
-# ln(permeability) = a * poro + b
-perm_a = 0.334
-perm_b = -1.655
+# porosity = (poro_0-poro_base) * exp(-poro_k * depth) + poro_base, with depth in km
+# values below correspond to the Sirte Basin
+poro_0 = 45
+poro_k = 0.4
+poro_b = 4.0
+poro_std = 1  # 0.01% uncertainty, if this goes larger it blows up massively...
+
+# ln(permeability) =
+# porperm: [-0.0092, 0.76, -6.7]
+perm_a = -0.0092
+perm_b = 0.76
+perm_c = -6.7
 
 # number of depths and number of samples per depth
 n_samples = 500
-n_depths = 10
+n_depths = 20
+
 
-depths = np.linspace(1e3, 5e3, n_depths)
+depth_min = 1e3
+depth_max = 3e3
+depths = np.linspace(depth_min, depth_max, n_depths)
 thickness = 100
 ntg = 1.0
 
-porosity = poro_a * np.exp(-poro_b * (depths * 1e-3))
+porosity = (poro_0-poro_b) * np.exp(-poro_k * (depths * 1e-3)) + poro_b
 
 # expand porosity to have n_samples per depth, with uncertainty in porosity
 porosity = np.random.normal(porosity[:, None], poro_std, size=(n_depths, n_samples))
 porosity = np.clip(porosity, 0.0, 100.0)
 
-permeability = np.exp(perm_a * porosity + perm_b)
+permeability = np.exp(perm_a * porosity**2 + perm_b * porosity +  perm_c)
+# clip to minimum transmissivity of 1Dm
+permeability = np.clip(permeability, 1000/thickness, 1e10)
 
 reservoir_properties = xr.Dataset(
     {
@@ -70,13 +79,18 @@ reservoir_properties = xr.Dataset(
     },
 )
 
+
+
 # For every sample, run a doublet simulation
 simulations = calculate_doublet_performance(
     reservoir_properties,
     chunk_size=100,
 )
+
+
 simulations_stoch = simulations.quantile([0.1,0.5,0.9], dim="samples")
 
+
 # make a plot showing there is a sweet spot between Doublet power, Cost of energy
 # and Transmissivity
 fig, ax = plt.subplots(figsize=(8, 6))
@@ -86,9 +100,11 @@ p10 = simulations_stoch.sel(quantile=0.9)
 p50 = simulations_stoch.sel(quantile=0.5)
 p90 = simulations_stoch.sel(quantile=0.1)
 
+
 # UTC
 ax2.fill_betweenx(p50.depth, p10.utc, p90.utc, color="lightgrey", alpha=0.5)
-ax2.plot(p50.utc, p50.depth, label="UTC [€/kWh]", color="lightgrey")
+ax2.plot(p50.utc, p50.depth, label="UTC [€cts/kWh]", color="lightgrey")
+ax2.set_xlim(0,50)
 
 # power
 ax2.fill_betweenx(p50.depth, p10.power, p90.power, color="black", alpha=0.5)
@@ -110,7 +126,7 @@ ax.plot(
 )
 
 ax.set_ylabel("Depth [m]")
-ax2.set_xlabel("Power [MW] & UTC [€/kWh]")
+ax2.set_xlabel("Power [MW] & UTC [€cts/kWh]")
 ax.set_xlabel("Transmissivity [Dm]", color="tab:green")
 ax.set_ylim([depths[-1], depths[0]])
 
diff --git a/docs/usage/deterministic_doublet.md b/docs/usage/deterministic_doublet.md
index fc8070ea5dfbdd101094e417b05e79f57cf64028..092a095bf0f0933bc9065b021aac1f4d0cd57e7a 100644
--- a/docs/usage/deterministic_doublet.md
+++ b/docs/usage/deterministic_doublet.md
@@ -21,7 +21,7 @@ import xarray as xr
 input_data = xr.Dataset({
     "thickness": 300,
     "ntg": 0.5,
-    "porosity": 0.5,
+    "porosity": 0.2,
     "depth": 2000,
     "permeability": 300,
 })
diff --git a/docs/usage/maprun_analysis.md b/docs/usage/maprun_analysis.md
index 1f5579330cf6b73cbbf2cc389b068aa62b883e99..7e4fb38c632ebb908311368f6211e563b72d4136 100644
--- a/docs/usage/maprun_analysis.md
+++ b/docs/usage/maprun_analysis.md
@@ -41,10 +41,10 @@ from pathlib import Path
 from os import path
 
 # the location of the input data: the data can be found in the resources/example_data directory of the repo
-input_data_path = Path(path.dirname(__file__), "resources") / "test_input" / "example_data"
+input_data_path = Path(__file__).parent.parent / "resources" / "test_input" / "example_data"
 
 # create a directory to write the output files to
-output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "example_data"
+output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example6"
 output_data_path.mkdir(parents=True, exist_ok=True)
 
 # if set to True then simulation is always run, otherwise pre-calculated results are read (if available)
@@ -60,6 +60,23 @@ reservoir_properties["depth"] = resample_grid(read_grid(input_data_path / "ROSL_
 reservoir_properties["ln_permeability_mean"] = np.log(resample_grid(read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"), new_cellsize=new_cellsize))
 reservoir_properties["ln_permeability_sd"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"), new_cellsize=new_cellsize)
 
+
+# output inputs
+variables_to_plot = ["depth", "thickness_mean", "thickness_sd"]
+fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=1, figsize=(7, 5*len(variables_to_plot)), sharex=True, sharey=True)
+for i, variable in enumerate(variables_to_plot):
+    plot_grid(reservoir_properties[variable], axes=axes[i], add_netherlands_shapefile=True)  # See documentation on plot_grid in pygridsio, you can also provide your own shapefile
+plt.tight_layout()  # ensure there is enough spacing
+plt.savefig(output_data_path / "input_maps_1.png")
+
+variables_to_plot = [ "porosity", "ln_permeability_mean", "ln_permeability_sd"]
+fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=1, figsize=(7, 5*len(variables_to_plot)), sharex=True, sharey=True)
+for i, variable in enumerate(variables_to_plot):
+    plot_grid(reservoir_properties[variable], axes=axes[i], add_netherlands_shapefile=True)  # See documentation on plot_grid in pygridsio, you can also provide your own shapefile
+plt.tight_layout()  # ensure there is enough spacing
+plt.savefig(output_data_path / "input_maps_2.png")
+
+
 # simulate
 # if doublet simulation has already been run then read in results, or run the simulation and write results out
 if (output_data_path / "output_results.nc").exists and not run_simulation:
@@ -68,6 +85,9 @@ else:
     results = calculate_doublet_performance_stochastic(reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90])
     results.to_netcdf(output_data_path / "output_results.nc")  # write doublet simulation results to file
 
+    
+
+
 # plot results as maps
 variables_to_plot = ["transmissivity", "power", "pres", "npv"]
 p_values_to_plot = [10, 50, 90]
@@ -103,6 +123,13 @@ plt.savefig(output_data_path / "pos.png")
 ```
 
 ---
+The figure of input maps generated by the above code will look like this:
+
+| ![Maps](../images/input_maps_1.png)    |  ![Maps](../images/input_maps_2.png) |
+|---|---|
+
+
+
 
 The figure generated by the above code will look like this:
 
diff --git a/docs/usage/portfoliorun_analysis.md b/docs/usage/portfoliorun_analysis.md
index 02cefca2d388757f8e79e04c5183c4df092bb323..e4516cbe67ef28bc9d9abdbe0a0275bd985a92d2 100644
--- a/docs/usage/portfoliorun_analysis.md
+++ b/docs/usage/portfoliorun_analysis.md
@@ -41,8 +41,8 @@ from pathlib import Path
 from os import path
 import numpy as np
 
-input_data_path = Path(path.dirname(__file__), "resources") / "test_input" / "example_data"
-output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "example_data"
+input_data_path = Path(__file__).parent.parent / "resources" / "test_input" / "example_data"
+output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example7"
 output_data_path.mkdir(parents=True, exist_ok=True)
 
 # read in reservoir property maps
diff --git a/docs/usage/stochastic_doublet.md b/docs/usage/stochastic_doublet.md
index 45fb276b8164de6d4f71282cf0b468903a5c83f7..5d5730f09694ea73b2151d8700f298007accbb92 100644
--- a/docs/usage/stochastic_doublet.md
+++ b/docs/usage/stochastic_doublet.md
@@ -30,7 +30,7 @@ input_data = xr.Dataset({
     "thickness_mean": 300,
     "thickness_sd": 50,
     "ntg": 0.5,
-    "porosity": 0.5,
+    "porosity": 0.2,
     "depth": 2000,
     "ln_permeability_mean": 4,
     "ln_permeability_sd": 0.5,
@@ -53,14 +53,14 @@ import xarray as xr
 from matplotlib import pyplot as plt
 from pathlib import Path
 from os import path
-output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "example_data"
+output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example5"
 output_data_path.mkdir(parents=True, exist_ok=True)
 
 input_data = xr.Dataset({
     "thickness_mean": 300,
     "thickness_sd": 50,
     "ntg": 0.5,
-    "porosity": 0.5,
+    "porosity": 0.2,
     "depth": 2000,
     "ln_permeability_mean": 5,
     "ln_permeability_sd": 0.5,
diff --git a/tests/documentation/resources/test_output/example9/depth_optimization.png b/tests/documentation/resources/test_output/example9/depth_optimization.png
deleted file mode 100644
index 9178f34d489c2bb1e52ad5d3e8876c543c7dcda9..0000000000000000000000000000000000000000
Binary files a/tests/documentation/resources/test_output/example9/depth_optimization.png and /dev/null differ
diff --git a/tests/documentation/test_doc_examples.py b/tests/documentation/test_doc_examples.py
index 3fe5f3c66cd31cb9b7c477b517bbd05e766824ac..896739baeb77c4c2b46f5389810c50ea5fe1cea3 100644
--- a/tests/documentation/test_doc_examples.py
+++ b/tests/documentation/test_doc_examples.py
@@ -1,3 +1,5 @@
+import pytest
+
 def test_example_1():
     from pythermogis import calculate_doublet_performance
     import xarray as xr
@@ -5,7 +7,7 @@ def test_example_1():
     input_data = xr.Dataset({
         "thickness": 300,
         "ntg": 0.5,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "depth": 2000,
         "permeability": 300,
     })
@@ -21,7 +23,7 @@ def test_example_2():
     input_data = xr.Dataset({
         "thickness": (("x", "y"), np.array([[300, 300], [200, 200]])),
         "ntg": (("x", "y"), np.array([[0.5, 0.5], [0.25, 0.25]])),
-        "porosity": (("x", "y"), np.array([[0.5, 0.5], [0.75, 0.7]])),
+        "porosity": (("x", "y"), np.array([[0.2, 0.2], [0.3, 0.3]])),
         "depth": (("x", "y"), np.array([[5000, 5000], [4500, 4500]])),
         "permeability": (("x", "y"), np.array([[300, 300], [450, 450]])),
     }, coords={"x": [0, 1], "y": [10, 20]})
@@ -36,7 +38,7 @@ def test_example_3():
     input_data = xr.Dataset({
         "thickness": 300,
         "ntg": 0.5,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "depth": 2000,
         "permeability": 300,
     })
@@ -53,7 +55,7 @@ def test_example_4():
         "thickness_mean": 300,
         "thickness_sd": 50,
         "ntg": 0.5,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "depth": 2000,
         "ln_permeability_mean": 4,
         "ln_permeability_sd": 0.5,
@@ -68,14 +70,14 @@ def test_example_5():
     from matplotlib import pyplot as plt
     from pathlib import Path
     from os import path
-    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example_data"
+    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example5"
     output_data_path.mkdir(parents=True, exist_ok=True)
 
     input_data = xr.Dataset({
         "thickness_mean": 300,
         "thickness_sd": 50,
         "ntg": 0.5,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "depth": 2000,
         "ln_permeability_mean": 5,
         "ln_permeability_sd": 0.5,
@@ -98,6 +100,7 @@ def test_example_5():
     axe.set_ylabel("p-value [%]")
     plt.savefig(output_data_path / "kh.png")
 
+@pytest.mark.skip("This test is computationally expensive, skip on the pipeline")
 def test_example_6():
     from pythermogis import calculate_doublet_performance_stochastic, calculate_pos
     from pygridsio import read_grid, plot_grid, resample_grid
@@ -111,14 +114,14 @@ def test_example_6():
     input_data_path = Path(__file__).parent.parent / "resources" / "test_input" / "example_data"
 
     # create a directory to write the output files to
-    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example_data"
+    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example6"
     output_data_path.mkdir(parents=True, exist_ok=True)
 
     # if set to True then simulation is always run, otherwise pre-calculated results are read (if available)
     run_simulation = True
 
     # grids can be in .nc, .asc, .zmap or .tif format: see pygridsio package
-    new_cellsize = 20000  # in m, this sets the resolution of the model; to speed up calcualtions you can increase the cellsize
+    new_cellsize = 5000  # in m, this sets the resolution of the model; to speed up calcualtions you can increase the cellsize
     reservoir_properties = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean")
     reservoir_properties["thickness_sd"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize)
     reservoir_properties["ntg"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize)
@@ -127,6 +130,22 @@ def test_example_6():
     reservoir_properties["ln_permeability_mean"] = np.log(resample_grid(read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"), new_cellsize=new_cellsize))
     reservoir_properties["ln_permeability_sd"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"), new_cellsize=new_cellsize)
 
+    # output inputs
+    variables_to_plot = ["depth", "thickness_mean", "thickness_sd"]
+    fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=1, figsize=(7, 5*len(variables_to_plot)), sharex=True, sharey=True)
+    for i, variable in enumerate(variables_to_plot):
+        plot_grid(reservoir_properties[variable], axes=axes[i], add_netherlands_shapefile=True)  # See documentation on plot_grid in pygridsio, you can also provide your own shapefile
+    plt.tight_layout()  # ensure there is enough spacing
+    plt.savefig(output_data_path / "input_maps_1.png")
+
+    variables_to_plot = [ "porosity", "ln_permeability_mean", "ln_permeability_sd"]
+    fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=1, figsize=(7, 5*len(variables_to_plot)), sharex=True, sharey=True)
+    for i, variable in enumerate(variables_to_plot):
+        plot_grid(reservoir_properties[variable], axes=axes[i], add_netherlands_shapefile=True)  # See documentation on plot_grid in pygridsio, you can also provide your own shapefile
+    plt.tight_layout()  # ensure there is enough spacing
+    plt.savefig(output_data_path / "input_maps_2.png")
+
+
     # simulate
     # if doublet simulation has already been run then read in results, or run the simulation and write results out
     if (output_data_path / "output_results.nc").exists and not run_simulation:
@@ -138,7 +157,7 @@ def test_example_6():
     # plot results as maps
     variables_to_plot = ["transmissivity","power", "pres", "npv"]
     p_values_to_plot = [10, 50, 90]
-    fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=len(p_values_to_plot), figsize=(5*len(variables_to_plot), 5*len(p_values_to_plot)), sharex=True, sharey=True)
+    fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=len(p_values_to_plot), figsize=(3*len(variables_to_plot), 5*len(p_values_to_plot)), sharex=True, sharey=True)
     for j, p_value in enumerate(p_values_to_plot):
         results_p_value = results.sel(p_value=p_value)
         for i, variable in enumerate(variables_to_plot):
@@ -178,7 +197,7 @@ def test_example_7():
     from pathlib import Path
 
     input_data_path = Path(__file__).parent.parent / "resources" / "test_input" / "example_data"
-    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example_data"
+    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example7"
     output_data_path.mkdir(parents=True, exist_ok=True)
 
     # read in reservoir property maps
@@ -265,7 +284,7 @@ def test_example_7():
     plt.tight_layout()  # ensure there is enough spacing
     plt.savefig(output_data_path / "npv.png")
 
-def test_example8():
+def test_example_8():
     from pythermogis import calculate_doublet_performance, calculate_pos, plot_exceedance
     import xarray as xr
     from matplotlib import pyplot as plt
@@ -275,13 +294,13 @@ def test_example8():
 
 
     # output directory to write output to
-    output_data_path =  Path(__file__).parent.parent / "resources" / "test_output" / "exceedance_example"
+    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example8"
     output_data_path.mkdir(parents=True, exist_ok=True)
 
     # generate simulation samples across desired reservoir properties
     Nsamples = 100
     thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples)
-    porosity_samples = np.random.uniform(low=0.5, high=0.8, size=Nsamples)
+    porosity_samples = np.random.uniform(low=0.2, high=0.2, size=Nsamples)
     ntg_samples = np.random.uniform(low=0.25, high=0.5, size=Nsamples)
     depth_samples = np.random.uniform(low=1800, high=2200,  size=Nsamples)
     permeability_samples = np.random.uniform(low=200, high=800, size=Nsamples)
@@ -331,13 +350,14 @@ def test_example8():
     plot_exceedance(results["utc"], ax=axes[0,1])
     plot_exceedance(results["power"], ax=axes[1,0])
     plot_exceedance(results["npv"], ax=axes[1,1])
-    axes[1, 1].axhline(results["pos"], ls="--", c="tab:orange", label=f"POS: {results["pos"]:.1f}%")   # add the probability of success
+    axes[1, 1].axhline(results["pos"], ls="--", c="tab:orange", label=f"POS: {results['pos']:.1f}%")   # add the probability of success
     axes[1, 1].axvline(0, ls="--", c="tab:orange")  # add the probability of success
     axes[1, 1].legend()
     plt.savefig(output_data_path / "exceedance_probabilities.png")
 
 
-def test_example9():
+@pytest.mark.skip("This test is computationally expensive, skip on the pipeline")
+def test_example_9():
     from pathlib import Path
 
     import numpy as np
@@ -347,34 +367,43 @@ def test_example9():
     from pythermogis import calculate_doublet_performance
 
     # output directory to write output to
-    output_data_path = Path(__file__).parent / "resources" / "test_output" / "example9"
+    output_data_path = Path(__file__).parent.parent / "resources" / "test_output" / "example9"
     output_data_path.mkdir(parents=True, exist_ok=True)
 
     # define the poro-depth and the poro-perm relationships:
-    # porosity = a * exp(-b * depth), with depth in km
-    poro_a = 28.97
-    poro_b = 0.21
-    poro_std = 1  # 1% uncertainty, if this goes larger it blows up massively...
-
-    # ln(permeability) = a * poro + b
-    perm_a = 0.334
-    perm_b = -1.655
+    # porosity = (poro_0-poro_base) * exp(-poro_k * depth) + poro_base, with depth in km
+    # values below correspond to the Sirte Basin
+    poro_0 = 45
+    poro_k = 0.4
+    poro_b = 4.0
+    poro_std = 1  # 0.01% uncertainty, if this goes larger it blows up massively...
+
+    # ln(permeability) =
+    # porperm: [-0.0092, 0.76, -6.7]
+    perm_a = -0.0092
+    perm_b = 0.76
+    perm_c = -6.7
 
     # number of depths and number of samples per depth
     n_samples = 500
-    n_depths = 10
+    n_depths = 20
+
 
-    depths = np.linspace(1e3, 5e3, n_depths)
+    depth_min = 1e3
+    depth_max = 3e3
+    depths = np.linspace(depth_min, depth_max, n_depths)
     thickness = 100
     ntg = 1.0
 
-    porosity = poro_a * np.exp(-poro_b * (depths * 1e-3))
+    porosity = (poro_0-poro_b) * np.exp(-poro_k * (depths * 1e-3)) + poro_b
 
     # expand porosity to have n_samples per depth, with uncertainty in porosity
     porosity = np.random.normal(porosity[:, None], poro_std, size=(n_depths, n_samples))
     porosity = np.clip(porosity, 0.0, 100.0)
 
-    permeability = np.exp(perm_a * porosity + perm_b)
+    permeability = np.exp(perm_a * porosity**2 + perm_b * porosity +  perm_c)
+    # clip to minimum transmissivity of 1Dm
+    permeability = np.clip(permeability, 1000/thickness, 1e10)
 
     reservoir_properties = xr.Dataset(
         {
@@ -389,13 +418,18 @@ def test_example9():
         },
     )
 
+
+
     # For every sample, run a doublet simulation
     simulations = calculate_doublet_performance(
         reservoir_properties,
         chunk_size=100,
     )
+
+
     simulations_stoch = simulations.quantile([0.1,0.5,0.9], dim="samples")
 
+
     # make a plot showing there is a sweet spot between Doublet power, Cost of energy
     # and Transmissivity
     fig, ax = plt.subplots(figsize=(8, 6))
@@ -405,9 +439,11 @@ def test_example9():
     p50 = simulations_stoch.sel(quantile=0.5)
     p90 = simulations_stoch.sel(quantile=0.1)
 
+
     # UTC
     ax2.fill_betweenx(p50.depth, p10.utc, p90.utc, color="lightgrey", alpha=0.5)
-    ax2.plot(p50.utc, p50.depth, label="UTC [€/kWh]", color="lightgrey")
+    ax2.plot(p50.utc, p50.depth, label="UTC [€cts/kWh]", color="lightgrey")
+    ax2.set_xlim(0,50)
 
     # power
     ax2.fill_betweenx(p50.depth, p10.power, p90.power, color="black", alpha=0.5)
@@ -429,7 +465,7 @@ def test_example9():
     )
 
     ax.set_ylabel("Depth [m]")
-    ax2.set_xlabel("Power [MW] & UTC [€/kWh]")
+    ax2.set_xlabel("Power [MW] & UTC [€cts/kWh]")
     ax.set_xlabel("Transmissivity [Dm]", color="tab:green")
     ax.set_ylim([depths[-1], depths[0]])
 
@@ -442,4 +478,3 @@ def test_example9():
 
     plt.savefig(output_data_path / "depth_optimization.png")
 
-
diff --git a/tests/java/test_ThermoGISDoublet_Benchmark.py b/tests/java/test_ThermoGISDoublet_Benchmark.py
index dca7e23a34e282f4b6e2ab958c331985a984b898..9dac69befb6774c144829fc719966d3b9b9c9b24 100644
--- a/tests/java/test_ThermoGISDoublet_Benchmark.py
+++ b/tests/java/test_ThermoGISDoublet_Benchmark.py
@@ -1,5 +1,7 @@
 from unittest import TestCase
 
+import pytest
+
 from pythermogis import *
 
 class ThermoGISDoubletBenchmark(TestCase):
@@ -219,6 +221,61 @@ class ThermoGISDoubletBenchmark(TestCase):
         self.assertTrue(np.isclose(10.401010755009017, results.utc(), 0.001))
         self.assertTrue(np.isclose(16.50359210062243, results.capex(), 0.001))
 
+    @pytest.mark.skip("This test requires a fix in the java core. ignore until fix is pushed")
+    def test_calculateDoubletPerformance_directheatHP_App(self):
+        """
+        This is a copy of a test from the Java ThermoGisDoubletTest.java script; to validate that this python implementation of the ThermoGIS Doublet
+        returns the same values.
+        """
+        # Arrange
+        start_jvm()
+        Logger = JClass("logging.Logger")
+        Mockito = JClass("org.mockito.Mockito")
+        RNG = JClass("tno.geoenergy.stochastic.RandomNumberGenerator")
+        ThermoGISDoublet = JClass("thermogis.calc.utc.doublet.ThermoGisDoublet")
+        DoubletInput = JClass("thermogis.calc.utc.doublet.records.DoubletInput")
+        ViscosityMode = JClass("tno.geoenergy.ViscosityMode")
+
+        utc_properties = (self.setup_template_utc_properties_builder()
+                          .setOpexPerPower(100)
+                          .setOpexBase(0)
+                          .setHpDirectHeatInputTemp(80)
+                          #.setHpApplicationMode(True)
+                          .setUseHeatPump(True)
+                          .setDhReturnTemp(50)
+                          .setViscosityMode(ViscosityMode.KESTIN)
+                          .build())
+
+        thickness = 100
+        permeability = 175
+        transmissivity = thickness * permeability
+        temperature = 50
+
+        input = DoubletInput(
+            -999.0, # unknowninput
+            thickness,
+            transmissivity,
+            0.0, # transmissivityWithNtg
+            1.0, # ntg
+            2000.0, # depth
+            0.0,# porosity
+            temperature,
+            None, # ates input
+        )
+        doublet = ThermoGISDoublet(Mockito.mock(Logger), RNG(0), utc_properties)
+
+        # Act
+        results = doublet.calculateDoubletPerformance(input)
+
+        # Assert
+        self.assertTrue(np.isclose(17499.99940, transmissivity, 0.001))
+        self.assertTrue(np.isclose(163.99771118164062, results.flow(), 1))
+        self.assertTrue(np.isclose(60, results.pres(), 0.001))
+        self.assertTrue(np.isclose(4.97566556930542, results.power(), 0.001))
+        self.assertTrue(np.isclose(989.2578125, results.welld(), 0.001))
+        self.assertTrue(np.isclose(4.383212200392486, results.cop(), 0.001))
+        self.assertTrue(np.isclose(10.401010755009017, results.utc(), 0.001))
+        self.assertTrue(np.isclose(16.50359210062243, results.capex(), 0.001))
 
     def test_calculateDoubletPerformance_ORC(self):
         """
diff --git a/tests/simulation/test_deterministic_doublet.py b/tests/simulation/test_deterministic_doublet.py
index 4aedecc0c8608580e2c5ef1b23ee1a68a9d7137d..be6caaea06bb7d903eb70ef349508d1dc62c6d0f 100644
--- a/tests/simulation/test_deterministic_doublet.py
+++ b/tests/simulation/test_deterministic_doublet.py
@@ -7,7 +7,7 @@ def test_deterministic_doublet():
 
     reservoir_properties = xr.Dataset({
         "thickness": 10,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "ntg": 0.5,
         "depth": 2000,
         "permeability": 200
@@ -21,7 +21,7 @@ def test_deterministic_doublet_mask_value():
 
     reservoir_properties = xr.Dataset({
         "thickness": 10,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "ntg": 0.5,
         "depth": 2000,
         "permeability": 200,
@@ -37,7 +37,7 @@ def test_deterministic_doublet_temperature_provided():
 
     reservoir_properties = xr.Dataset({
         "thickness": 10,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "ntg": 0.5,
         "depth": 2000,
         "temperature": 89,
@@ -53,7 +53,7 @@ def test_deterministic_doublet_transmissivity_provided():
 
     reservoir_properties = xr.Dataset({
         "thickness": 10,
-        "porosity": 0.5,
+        "porosity": 0.2,
         "ntg": 0.5,
         "depth": 2000,
         "temperature": 89,