diff --git a/poisson/discrete/finite_volume.py b/poisson/discrete/finite_volume.py
index ead6170bfb3456ffb15638645b396dc31da4e877..30df399d984ea291669dc0faededc2dbd6c52067 100644
--- a/poisson/discrete/finite_volume.py
+++ b/poisson/discrete/finite_volume.py
@@ -14,7 +14,12 @@ import warnings
 import itertools
 from collections import namedtuple
 
-from scipy.spatial import qhull
+from scipy.spatial import ConvexHull
+# TODO: remove the options once the code depends on scipy > 1.8
+try:
+    from scipy.spatial.qhull import _Qhull
+except ImportError:
+    from scipy.spatial._qhull import _Qhull
 
 import numpy as np
 
@@ -275,7 +280,7 @@ class Voronoi(FiniteVolume):
             -----------
                 vor is a voronoi object comming from:
                     scipy.spatial.Voronoi,
-                    scipy.spatial.qhull_Qhull.get_voronoi_diagram,
+                    scipy.spatial.qhull._Qhull.get_voronoi_diagram,
                     poisson.poisson.VoronoiMesh._make_voronoi,
 
                 check: boolean. True to test each VoronoiMesh property
@@ -430,8 +435,7 @@ class Voronoi(FiniteVolume):
             qhull_options = b"Qbb Qc Qx Qz"
 
         # Initialize Qhull object
-        qhull_obj = qhull._Qhull(b"v", grid, qhull_options,
-                                 **kwargs)
+        qhull_obj = _Qhull(b"v", grid, qhull_options, **kwargs)
 
         # FOr incremental processing, e.g. add points use qhull._QhullUser
         voronoi = namedtuple('voronoi_prop', ['vertices', 'ridge_points',
@@ -824,7 +828,7 @@ class Voronoi(FiniteVolume):
             regs = [[self.vertices[reg]
                      for reg in self.region_vertices[i]]
                     for i in index_of_regions]
-            volume = np.array([qhull.ConvexHull(reg).volume for reg in regs])
+            volume = np.array([ConvexHull(reg).volume for reg in regs])
             return volume