diff --git a/magnet_on_bentheimer.py b/magnet_on_bentheimer.py
deleted file mode 100644
index fc132b56e..000000000
--- a/magnet_on_bentheimer.py
+++ /dev/null
@@ -1,178 +0,0 @@
-import pandas as pd
-import porespy as ps
-import openpnm as op
-import numpy as np
-import scipy.ndimage as spim
-from edt import edt
-from skimage.morphology import square, cube
-
-
-# %% Set simulation parameters
-np.random.seed(10)
-export = False
-fill_blind_pores = True
-padding = None
-l_max=None
-bw = None
-endpoints=True
-name = 'Bentheimer'
-res = 5.82e-6
-Kx, Ky, Kz, Kavg = [], [], [], []
-porosity = []
-
-# %% Import and prepare image
-f = r"C:\Users\jeff\OneDrive - University of Waterloo\Research\Elidek - Bentheimer\bentheimer_binary.npz"
-imb = np.load(f)['im'][:500, :500, :500]
-imb = ps.filters.fill_blind_pores(imb, conn=6, surface=True)
-imb = ps.filters.trim_floating_solid(imb, conn=6, surface=True)
-im = imb
-porosity.append(np.sum(im)/np.product(im.shape)*100)
-print(f'porosity: {np.sum(im)/np.product(im.shape)*100}')
-
-# %% Get distance transform
-dt = edt(im, parallel=8)
-
-# %% Get Skeleton
-# sk = ps.networks.skeletonize_magnet2(im)
-sk, im = ps.networks.skeleton(im, padding=None, surface=True, parallel=False)
-labels, N = spim.label(sk, structure=ps.tools.ps_rect(3, 3))
-values, counts = np.unique(labels, return_counts=True)
-if N > 1:
-    print(f"There are {N} separated segments in the skeleton")
-    print(f"Their sizes are {counts[1:]}")
-
-# %% Find junctions and end points on skeleton
-juncs, endpts = ps.networks.analyze_skeleton(sk)
-pores, throats = ps.networks.partition_skeleton(sk, juncs + endpts, dt)
-new_juncs, pores, new_throats =\
-    ps.networks.find_throat_junctions(im=im, pores=pores, throats=throats, dt=dt)
-# pores, throats = ps.networks.partition_skeleton(sk, juncs + new_juncs + endpts, dt=dt)
-
-# net = ps.networks.sk_to_network(pores, throats, dt)
-
-# %% insert pores at junction points
-pts = (pores + new_juncs) > 0
-fbd = ps.networks.insert_pore_bodies(sk, dt, pts)
-# convert spheres to network dictionary
-net = ps.networks.spheres_to_network(
-    sk=sk, dt=dt, fbd=fbd, voxel_size=res, boundary_width=0)
-
-
-# %% import network to openpnm
-net = op.io.network_from_porespy(net)
-net['pore.diameter'] = net['pore.radius']*2
-net['throat.diameter'] = net['throat.radius']*2
-
-# check network health
-h = op.utils.check_network_health(net)
-op.topotools.trim(net, pores=np.append(h['disconnected_pores'], h['isolated_pores']))
-h = op.utils.check_network_health(net)
-
-# add geometry models
-geo = op.models.collections.geometry.cubes_and_cuboids.copy()
-del geo['pore.diameter'], geo['throat.diameter']
-net['pore.diameter'] = net['pore.diameter'].copy()
-net['throat.diameter'] = net['throat.diameter'].copy()
-net.add_model_collection(geo)
-net.regenerate_models()
-
-# add phase
-phase = op.phase.Phase(network=net)
-phase['pore.viscosity'] = 1e-3
-
-# physics
-phys = op.models.collections.physics.basic.copy()
-phase.add_model_collection(phys)
-phase.regenerate_models()
-
-# stokes flow simulation to estimate permeability
-Pin = 5e-6
-Pout = 0
-A = (im.shape[0]*im.shape[1]) * res**2
-L = im.shape[1] * res
-mu = phase['pore.viscosity'].max()
-
-# label boundary pores
-xmin = net.coords[:, 0] <= net.coords[:, 0].max()*0.02
-xmax = net.coords[:, 0] >= net.coords[:, 0].max()*0.98
-flow_x = op.algorithms.StokesFlow(network=net, phase=phase)
-flow_x.set_value_BC(pores=xmin, values=Pin)
-flow_x.set_value_BC(pores=xmax, values=Pout)
-flow_x.run()
-
-Q_x = flow_x.rate(pores=xmin, mode='group')[0]
-K_x = Q_x * L * mu / (A * (Pin - Pout))/0.98e-12*1000
-print(f'K_x is: {K_x:.2f} mD')
-
-ymin = net.coords[:, 1] <= net.coords[:, 1].max()*0.02
-ymax = net.coords[:, 1] >= net.coords[:, 1].max()*0.98
-flow_y = op.algorithms.StokesFlow(network=net, phase=phase)
-flow_y.set_value_BC(pores=ymin, values=Pin)
-flow_y.set_value_BC(pores=ymax, values=Pout)
-flow_y.run()
-
-Q_y = flow_y.rate(pores=ymin, mode='group')[0]
-K_y = Q_y * L * mu / (A * (Pin - Pout))/0.98e-12*1000
-print(f'K_y is: {K_y:.2f} mD')
-
-zmin = net.coords[:, 2] <= net.coords[:, 2].max()*0.02
-zmax = net.coords[:, 2] >= net.coords[:, 2].max()*0.98
-flow_z = op.algorithms.StokesFlow(network=net, phase=phase)
-flow_z.set_value_BC(pores=zmin, values=Pin)
-flow_z.set_value_BC(pores=zmax, values=Pout)
-flow_z.run()
-
-Q_z = flow_z.rate(pores=zmin, mode='group')[0]
-K_z = Q_z * L * mu / (A * (Pin - Pout))/0.98e-12*1000
-print(f'K_z is: {K_z:.2f} mD')
-
-K = np.average([K_x, K_y, K_z])
-print(f'K is: {K:.2f} mD')
-
-# number of pore vs. skeleton clusters in network
-from scipy.sparse import csgraph as csg
-am = net.create_adjacency_matrix(fmt='coo', triu=True)
-Np, cluster_num = csg.connected_components(am, directed=False)
-print('Pore clusters:', Np)
-# number of skeleton pieces
-b = square(3) if im.ndim == 2 else cube(3)
-_, Ns = spim.label(input=sk.astype('bool'), structure=b)
-print('Skeleton clusters:', Ns)
-
-print('The expected permeability is 2490 mD')
-
-# append data
-Kx.append(K_x)
-Ky.append(K_y)
-Kz.append(K_z)
-Kavg.append(K)
-
-
-# export
-if export:
-    ps.io.to_stl(~sk, name + '-sk')
-    ps.io.to_stl(im, name + '-im')
-    net['pore.coords'] = net['pore.coords']/res
-    net['pore.diameter'] = net['pore.diameter']/res
-    net['throat.radius'] = net['throat.radius']/res
-    net['pore.coords'] += 10  # res
-    proj = net.project
-    op.io.project_to_xdmf(proj, filename=name + '-network')
-
-# create and export data frame
-if 0:
-    data = {
-        "Name": name,
-        "Resolution": res,
-        "boundary_width": bw,
-        "l_max": l_max,
-        "padding": str(padding),
-        "endpoints": endpoints,
-        "porosity": porosity,
-        "K_x": Kx,
-        "K_y": Ky,
-        "K_z": Kz,
-        "K": Kavg
-    }
-    df = pd.DataFrame(data=data)
-    df.to_csv('MAGNET Permeability.csv', index=False)
diff --git a/src/porespy/filters/_funcs.py b/src/porespy/filters/_funcs.py
index 9f9430cf4..f3d579246 100644
--- a/src/porespy/filters/_funcs.py
+++ b/src/porespy/filters/_funcs.py
@@ -688,21 +688,21 @@ def local_thickness(
         then parallel processing does not occur.  `2` is equivalent to
         `[2, 2, 2]` for a 3D image.  The number of cores used is specified in
         `porespy.settings.ncores` and defaults to all cores.
-        
+
     parallel_kw : dict
         Dictionary containing the settings for parallelization by chunking. The
         optional settings include divs (scalar or list of scalars,
         default = 1), overlap (scalar or list of scalars, unused in this func),
         and cores (scalar, default is all available cores).
-        
+
         Divs is the number of times to divide the image for parallel
         processing. If `1` then parallel processing does not occur. `2` is
         equivalent to `[2, 2, 2]` for a 3D image.
-        
+
         Overlap is the amount of overlap to apply between chunks. In this
         function, the overlap is controlled by the distance transform and
         cannot be altered in any way.
-        
+
         Cores is the number of cores that will be used to parallel process all
         domains. If ``None`` then all cores will be used but user can specify
         any integer values to control the memory usage. Setting value to 1 will
@@ -865,21 +865,21 @@ def porosimetry(
         equivalent to `[2, 2, 2]` for a 3D image.  The number of cores
         used is specified in `porespy.settings.ncores` and defaults to
         all cores.
-    
+
     parallel_kw : dict
         Dictionary containing the settings for parallelization by chunking. The
         optional settings include divs (scalar or list of scalars,
         default = 1), overlap (scalar or list of scalars, unused in this func),
         and cores (scalar, default is all available cores).
-        
+
         Divs is the number of times to divide the image for parallel
         processing. If `1` then parallel processing does not occur. `2` is
         equivalent to `[2, 2, 2]` for a 3D image.
-        
+
         Overlap is the amount of overlap to apply between chunks. In this
         function, the overlap is controlled by the distance transform and
         cannot be altered in any way.
-        
+
         Cores is the number of cores that will be used to parallel process all
         domains. If ``None`` then all cores will be used but user can specify
         any integer values to control the memory usage. Setting value to 1 will
@@ -1221,23 +1221,23 @@ def chunked_func(
     func : function handle
         The function which should be applied to each chunk, such as
         `spipy.ndimage.binary_dilation`.
-        
+
     parallel_kw : dict
         Dictionary containing the settings for parallelization by chunking. The
         optional settings include divs (scalar or list of scalars,
         default = [2, 2, 2]), overlap (scalar or list of scalars, optional),
         and cores (scalar, default is all available cores).
-        
+
         Divs is the number of times to divide the image for parallel
         processing. If `1` then parallel processing does not occur. `2` is
         equivalent to `[2, 2, 2]` for a 3D image.
-        
+
         Overlap is the amount of overlap to include when dividing up the image.
         This value will almost always be the size (i.e. raduis) of the
         structuring element. If not specified then the amount of overlap
         is inferred from the size of the structuring element, in which
         case the `strel_arg` must be specified.
-        
+
         Cores is the number of cores that will be used to parallel process all
         domains. If ``None`` then all cores will be used but user can specify
         any integer values to control the memory usage. Setting value to 1 will
diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py
index f4d29c7d6..1094e1369 100644
--- a/src/porespy/metrics/_funcs.py
+++ b/src/porespy/metrics/_funcs.py
@@ -5,15 +5,16 @@
 import scipy.ndimage as spim
 import scipy.spatial as sptl
 import scipy.stats as spst
-from skimage.morphology import skeletonize
 from numba import njit
 from scipy import fft as sp_ft
 from skimage.measure import regionprops
+from skimage.morphology import skeletonize, ball, disk, square, cube
 from porespy import settings
 from porespy.filters import (
     local_thickness,
     pc_to_seq,
     fill_closed_pores,
+    find_invalid_pores,
 )
 from porespy.tools import (
     Results,
@@ -21,7 +22,8 @@
     extend_slice,
     get_tqdm,
     get_slices_slabs,
-    get_edt
+    get_edt,
+    ps_round,
 )
 
 
@@ -30,13 +32,16 @@
     "chord_counts",
     "chord_length_distribution",
     "find_h",
+    "is_percolating",
     "lineal_path_distribution",
     "pore_size_distribution",
     "radial_density_distribution",
     "porosity",
+    "find_porosity_threshold",
     "porosity_profile",
     "satn_profile",
     "two_point_correlation",
+    "percolating_porosity",
     "phase_fraction",
     "pc_curve",
     "pc_map_to_pc_curve",
@@ -47,6 +52,169 @@
 edt = get_edt()
 tqdm = get_tqdm()
 logger = logging.getLogger(__name__)
+strel = {2: {'min': disk(1), 'max': square(3)}, 3: {'min': ball(1), 'max': cube(3)}}
+
+
+def is_percolating(im, axis=None, inlets=None, outlets=None, conn='min'):
+    R"""
+    Determines if a percolating path exists across the domain (in the specified
+    direction) or between given inlets and outlets.
+
+    Parameters
+    ----------
+    im : ndarray
+        Image of the void space with `True` indicating void space.
+    axis : int
+        The axis along which percolation is checked. If `None` (default) then
+        percolation is checked in all dimensions.
+    conn : str
+        Can be either `'min'` or `'max'` and controls the shape of the structuring
+        element used to determine voxel connectivity.  The default if `'min'` which
+        imposes the strictest criteria, so that voxels must share a face to be
+        considered connected.
+
+    Returns
+    -------
+    percolating : bool or list of bools
+        A boolean value indicating if the domain percolates in the given direction.
+        If `axis=None` then all directions are checked and the result is returned
+        as a list like `[True, False, True]` indicating that the domain percolates
+        in the `x` and `z` directions, but not `y`.
+
+    """
+    if (inlets is not None) and (outlets is not None):
+        pass
+    elif axis is not None:
+        im = np.swapaxes(im, 0, axis) == 1
+        inlets = np.zeros_like(im)
+        inlets[0, ...] = True
+        inlets *= im
+        outlets = np.zeros_like(im)
+        outlets[-1, ...] = True
+        outlets *= im
+    else:
+        ans = []
+        for ax in range(im.ndim):
+            ans.append(is_percolating(im, axis=ax, conn=conn))
+        return ans
+
+    labels, N = spim.label(im, structure=strel[im.ndim][conn])
+    a = np.unique(labels[inlets])
+    a = a[a > 0]
+    b = np.unique(labels[outlets])
+    b = b[b > 0]
+    hits = np.isin(a, b)
+    return np.any(hits)
+
+
+def find_porosity_threshold(im, axis=0, conn='min'):
+    r"""
+    Finds the porosity of the image at the percolation threshold
+
+    This function progressively dilates the solid and reports the porosity at the
+    step just before there are no percolating paths (in the specified direction)
+
+    Parameters
+    ----------
+    im : ndarray
+        Image of the void space with `True` indicating void space
+    axis : int
+        The axis along which percolation is checked
+    conn : str
+        Can be either `'min'` or `'max'` and controls the shape of the structuring
+        element used to determine voxel connectivity.  The default if `'min'` which
+        imposes the strictest criteria, so that voxels must share a face to be
+        considered connected.
+
+    Returns
+    -------
+    results
+        A Results object with the following attributes:
+
+        ================ ===========================================================
+        Attribute        Description
+        ================ ===========================================================
+        eps_orig         The total porosity of the original image, including closed
+                         and surface pores
+        eps_orig_perc    The percolating porosity of the original image (i.e. with
+                         closed and surface pores filled)
+        eps_thresh       The total porosity of the image after eroding the void
+                         space results in no percolating paths
+        eps_thresh_perc  The percolating porosity of the eroded image (with closed
+                         and surface pores filled)
+        R                The threshold to apply to the distance transform to
+                         obtain the percolating image (i.e. im = dt >= R)
+        ================ ===========================================================
+    """
+    if axis is None:
+        raise Exception('axis must be specified')
+
+    def _check_percolation(dt, R, step, axis, conn):
+        while True:
+            im2 = dt >= R
+            check = np.array(is_percolating(im2, axis=axis, conn=conn))
+            if not np.all(check):
+                break
+            R += step
+        return R
+
+    dt = edt(im)
+
+    R = _check_percolation(dt, R=1, step=10, axis=axis, conn=conn)
+    R = _check_percolation(dt, R=max(1, R-10), step=4, axis=axis, conn=conn)
+    R = _check_percolation(dt, R=max(1, R-4), step=1, axis=axis, conn=conn)
+
+    im2 = dt >= (R - 1)
+    eps_thresh_total = porosity(im2)
+    eps_thresh_perc = percolating_porosity(im2, axis=axis)
+
+    from porespy.tools import Results
+    r = Results()
+    r.eps_orig = porosity(im)
+    r.eps_orig_perc = percolating_porosity(im, axis=axis)
+    r.eps_thresh = eps_thresh_total
+    r.eps_thresh_perc = eps_thresh_perc
+    r.R = R - 1
+    return r
+
+
+def percolating_porosity(im, axis=0, inlets=None, outlets=None, conn='min'):
+    r"""
+    Finds volume fraction of void space which belongs to percolating paths
+    across the domain in the direction specified.
+
+    Parameters
+    ----------
+    im : ndarray
+        Image of the void space with `True` indicating void space
+    axis : int
+        The axis along which percolation is checked
+    conn : str
+        Can be either `'min'` or `'max'` and controls the shape of the structuring
+        element used to determine voxel connectivity.  The default if `'min'` which
+        imposes the strictest criteria, so that voxels must share a face to be
+        considered connected.
+    inlets, outlets : ndarray
+        Boolean arrays indicating the locations of the inlets and outlets. These
+        are useful if the domain is not cubic or if special inlet and outlet
+        locations are desired.
+    """
+    se = strel[im.ndim][conn]
+    if (inlets is None) and (outlets is None):
+        im2 = np.swapaxes(im, 0, axis)
+        inlets = np.zeros_like(im2, dtype=bool)
+        inlets[0, ...] = True
+        outlets = np.zeros_like(im2, dtype=bool)
+        outlets[-1, ...] = True
+    labels, N = spim.label(im2, structure=se)
+    a = np.unique(labels*inlets)
+    a = a[a > 0]
+    b = np.unique(labels*outlets)
+    b = b[b > 0]
+    hits = np.intersect1d(a, b)
+    im3 = np.isin(labels, hits)
+    eps = porosity(im3)
+    return eps
 
 
 def boxcount(im, bins=10):
diff --git a/test/unit/test_metrics.py b/test/unit/test_metrics.py
index 9f16da3a3..686186273 100644
--- a/test/unit/test_metrics.py
+++ b/test/unit/test_metrics.py
@@ -514,6 +514,34 @@ def test_bond_number(self):
             source='dt', method='max')
         assert np.around(bo, decimals=5) == 0.79461
 
+    def test_is_percolating(self):
+        im = np.ones([20, 20], dtype=bool)
+        im[:, 10:12] = False
+        im[10, 10] = True
+        im[11, 11] = True
+        assert ps.metrics.is_percolating(im=im, axis=0)
+        assert not ps.metrics.is_percolating(im=im, axis=1, conn='min')
+        assert ps.metrics.is_percolating(im=im, axis=1, conn='max')
+        assert ps.metrics.is_percolating(im=im, axis=None) == [True, False]
+        assert ps.metrics.is_percolating(im=im, axis=None, conn='max') == [True, True]
+
+    def test_pecolating_porosity(self):
+        im = np.ones([20, 20], dtype=bool)
+        im[:10, 10] = False
+        im[10:, 11] = False
+        e = im.sum()/im.size
+        assert ps.metrics.percolating_porosity(im, axis=0, conn='min') == e
+        assert ps.metrics.percolating_porosity(im, axis=1, conn='min') == 0
+        assert ps.metrics.percolating_porosity(im, axis=0, conn='max') == e
+        assert ps.metrics.percolating_porosity(im, axis=1, conn='max') == e
+
+    def test_find_porosity_threshold(self):
+        im = ps.generators.lattice_spheres([31, 51], r=15, offset=5, spacing=30)
+        ep = ps.metrics.find_porosity_threshold(im, axis=0)
+        assert ep.eps_thresh == 0.4484503478810879
+        ep = ps.metrics.find_porosity_threshold(im, axis=1)
+        assert ep.eps_thresh == 0.05439595192915876
+
 
 if __name__ == '__main__':
     t = MetricsTest()