diff --git a/docs/changelog.md b/docs/changelog.md
index 24a8d1cafe6aab804cfa7ccb68d6c8e487684cd5..05f4a84e9e2fbcdaa72a5a42b2ace1ee574050cf 100644
--- a/docs/changelog.md
+++ b/docs/changelog.md
@@ -13,6 +13,8 @@ We also recognise that having to replace the Jar file everytime we update it is
 Fixing a bug caused by the update in the Jar file. calculateDoubletPerformance in the Java source code returns Null if one of the optimizers failed. 
 Pythermogis was still expecting a DoubletOutput class with all values set to -9999 in case of optimizer failure. 
 
-Additionally:
-    - two functions have been added to print out the environment variables THERMOGIS_JAR and JAVA_HOME. This enables users to debug if their are issues with their installation
+Another bug was fixed, this occured when sampling for transmissivity with high or low p-values which caused an indexing error. 
+
+Additional Functionality:
+    - two functions have been added to print out the environment variables THERMOGIS_JAR and JAVA_HOME. This enables users to debug if their are issues with their installation.
 
diff --git a/pixi.lock b/pixi.lock
index b3ca5554899697e34e7eff5bac6c361dee0ddc9b..8de774a2a77285438613385afa57934bedd4a668 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -4830,7 +4830,7 @@ packages:
 - pypi: ./
   name: pythermogis
   version: 1.2.1
-  sha256: 0ae9336fdba95931a8e5b2b0ce4a61acd3bcd3c3d685b2292b11e436d43245b0
+  sha256: adc786bc376cdca6369d5cb58bd4bc7597c9c00a620cbeef48d591239c4ad966
   requires_dist:
   - jpype1>=1.5.2,<2
   - xarray==2024.9.0.*
@@ -4839,6 +4839,7 @@ packages:
   - build>=1.2.2.post1,<2
   - pygridsio>=1.0,<2.0
   - python-dotenv>=1.2.1,<2
+  - dask>=2025.5.1,<2026
   requires_python: '>=3.11'
   editable: true
 - conda: https://conda.anaconda.org/conda-forge/linux-64/python-3.13.3-hf636f53_101_cp313.conda
diff --git a/pyproject.toml b/pyproject.toml
index 0d02ebf473306311ffbd3c57a294d9a5e8da7f23..60532c27e59b82d0e5af179e79d84ba6ed08eaa2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -21,7 +21,8 @@ dependencies = [
     "pytz>=2024.1,<2025",
     "build>=1.2.2.post1,<2",
     "pygridsio>=1.0,<2.0",
-    "python-dotenv>=1.2.1,<2"
+    "python-dotenv>=1.2.1,<2",
+    "dask>=2025.5.1,<2026"
 ]
 
 
diff --git a/src/pythermogis/transmissivity/calculate_thick_perm_trans.py b/src/pythermogis/transmissivity/calculate_thick_perm_trans.py
index aa485e7c230112f5b9149ef0388e66f43fdd6025..9b3bd30b506352900cfe1d49334c4d155f009c90 100644
--- a/src/pythermogis/transmissivity/calculate_thick_perm_trans.py
+++ b/src/pythermogis/transmissivity/calculate_thick_perm_trans.py
@@ -49,7 +49,7 @@ def generate_thickness_permeability_transmissivity_for_pvalues(thickness_mean: f
     # Sampling method for transmissivity
     transmissivity_samples = np.sort(np.exp(ln_permeability_dist.rvs(nSamples) + np.log(thickness_samples)))
 
-    sample_indexes = np.array(p_value_fractions * nSamples)
+    sample_indexes = np.array(p_value_fractions * (nSamples-1))
     transmissivity_pvalues_sampled = transmissivity_samples[sample_indexes.astype(int)]
 
     return thickness_pvalues, permeability_pvalues, transmissivity_pvalues_sampled
\ No newline at end of file
diff --git a/tests/test_doublet_speed.py b/tests/test_doublet_speed.py
index 36be0e513366e375dc6c38fae7c426740c70226f..b9da7a4804528acf37c5e57b9b485349f55d69ca 100644
--- a/tests/test_doublet_speed.py
+++ b/tests/test_doublet_speed.py
@@ -21,7 +21,7 @@ class PyThermoGIS(TestCase):
 
     def test_doublet_speed(self):
         input_grids = self.read_input_grids()
-        p_values = [10,50,90]
+        p_values = [99]
 
         data = input_grids["thickness_mean"].data.flatten()
         non_nan_data = data[~np.isnan(data)]
@@ -29,13 +29,13 @@ class PyThermoGIS(TestCase):
 
         # Run calculation across all dimensions of input_grids, and all provided P_values: The Equivalent calculation with 10 cores takes 43 seconds in the Java code
         start = timeit.default_timer()
-        calculate_doublet_performance_stochastic(input_grids, chunk_size=100, rng_seed=123, p_values=p_values)
+        calculate_doublet_performance_stochastic(input_grids, chunk_size=5, rng_seed=123, p_values=p_values)
         time_elapsed = timeit.default_timer() - start
         print(f"Python calculation took {time_elapsed:.1f} seconds.")
 
 
     def read_input_grids(self):
-        new_cellsize=5000 # in m
+        new_cellsize=10000 # in m
         input_grids = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean")
         input_grids["thickness_sd"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize)
         input_grids["ntg"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize)
diff --git a/tests/test_transmissivity_calculation.py b/tests/test_transmissivity_calculation.py
index aaaf24d3d51168dd113aa2268ce3e9119679e553..ac4e08cffda4c228e0a5e69390895639e439c44b 100644
--- a/tests/test_transmissivity_calculation.py
+++ b/tests/test_transmissivity_calculation.py
@@ -23,7 +23,7 @@ class PyThermoGIS(TestCase):
         ln_perm_mean = 5
         ln_perm_sd = 0.5
         nSamples = 100000
-        p_values_list = [10, 20, 30, 40, 50, 60, 70, 80, 90]
+        p_values_list = [1, 10, 20, 30, 40, 50, 60, 70, 80, 90]
 
         #
         #   Python