diff --git a/src/porespy/filters/_displacement.py b/src/porespy/filters/_displacement.py
index 24bcf5397..0b4c32d49 100644
--- a/src/porespy/filters/_displacement.py
+++ b/src/porespy/filters/_displacement.py
@@ -212,6 +212,7 @@ def find_trapped_clusters(
     im: npt.ArrayLike,
     seq: npt.ArrayLike,
     outlets: npt.ArrayLike,
+    min_size: int = 0,
     conn: Literal["min", "max"] = "min",
     method: Literal["queue", "labels"] = "labels",
 ):
@@ -231,6 +232,9 @@ def find_trapped_clusters(
     outlets : ndarray
         An image the same size as ``im`` with ``True`` indicating outlets
         and ``False`` elsewhere.
+    min_size : scalar
+        The threshold size of clusters.  Clusters with this many voxels or fewer
+        will be ignored.
     conn : str
         Controls the shape of the structuring element used to determine if voxels
         are connected.  Options are:
@@ -288,7 +292,12 @@ def find_trapped_clusters(
     else:
         raise Exception(f"{method} is not a supported method")
 
-    return (seq_temp == -1) * im
+    trapped = (seq_temp == -1) * im
+
+    if min_size > 0:
+        trapped = trim_small_clusters(im=trapped, min_size=min_size)
+
+    return trapped
 
 
 def _find_trapped_clusters_labels(
diff --git a/src/porespy/filters/_lt_methods.py b/src/porespy/filters/_lt_methods.py
index 5ec7b346f..9aab5b3b0 100644
--- a/src/porespy/filters/_lt_methods.py
+++ b/src/porespy/filters/_lt_methods.py
@@ -90,8 +90,8 @@ def porosimetry(
         from porespy.simulations import drainage_dt
         drn = drainage_dt(im=im, inlets=inlets, steps=sizes, smooth=smooth)
     elif method in ['dsi', 'bf']:
-        from porespy.simulations import drainage_dsi
-        drn = drainage_dsi(im=im, inlets=inlets, steps=sizes, smooth=smooth)
+        from porespy.simulations import drainage_bf
+        drn = drainage_bf(im=im, inlets=inlets, steps=sizes, smooth=smooth)
     if method in ['fft', 'conv']:
         from porespy.simulations import drainage_fft
         drn = drainage_fft(im=im, inlets=inlets, steps=sizes, smooth=smooth)
diff --git a/src/porespy/simulations/_drainage.py b/src/porespy/simulations/_drainage.py
index 60a2845b9..fd38a94ae 100644
--- a/src/porespy/simulations/_drainage.py
+++ b/src/porespy/simulations/_drainage.py
@@ -8,7 +8,7 @@
     fftmorphology,
     find_small_clusters,
     find_trapped_clusters,
-    pc_to_satn,
+    seq_to_satn,
     trim_disconnected_voxels,
 )
 from porespy.metrics import pc_map_to_pc_curve
@@ -32,7 +32,7 @@
     "drainage_dt",
     "drainage_fft",
     "drainage_dt_fft",
-    "drainage_dsi",
+    "drainage_bf",
 ]
 
 
@@ -41,7 +41,7 @@
 strel = get_strel()
 
 
-def drainage_dsi(
+def drainage_bf(
     im,
     inlets=None,
     outlets=None,
@@ -477,19 +477,14 @@ def drainage(
     outlets : ndarray, optional
         A boolean image with ``True`` values indicating the outlet locations.
         If this is provided then trapped voxels of wetting phase are found and
-        all the output images are adjusted accordingly. Note that trapping can
-        be assessed as a postprocessing step as well, so if this is not provided
-        trapping can still be considered.
+        all the output images are adjusted accordingly.
     residual : ndarray, optional
         A boolean array indicating the locations of any residual invading
         phase. This is added to the intermediate image prior to trimming
         disconnected clusters, so will create connections to some clusters
         that would otherwise be removed. The residual phase is indicated
         in the capillary pressure map by ``-np.inf`` values, since these voxels
-        are invaded at all applied capillary pressures. Note that the presence of
-        residual non-wetting phase makes it impossible to correctly deal with
-        trapping, so if both `residual` and `outlets` are given then an error is
-        raised.
+        are invaded at all applied capillary pressures.
     steps : int or array_like (default = 25)
         The range of pressures to apply. If an integer is given then the given
         number of steps will be created between the lowest and highest values in
@@ -541,12 +536,6 @@ def drainage(
                    value of capillary pressure (``pc``).
         ========== ============================================================
 
-    Notes
-    -----
-    This algorithm only provides sensible results for gravity stabilized
-    configurations, meaning the more dense fluid is on the bottom. Be sure that
-    ``inlets`` are specified accordingly.
-
     References
     ----------
     .. [4] Chadwick EA, Hammen LH, Schulz VP, Bazylak A, Ioannidis MA, Gostick JT.
@@ -562,19 +551,16 @@ def drainage(
     to view online example.
 
     """
-    if (residual is not None) and (outlets is not None):
-        raise Exception("Trapping cannot be properly assessed if residual present")
-
     im = np.array(im, dtype=bool)
 
-    if dt is None:
-        dt = edt(im)
-
     if outlets is not None:
         outlets = outlets * im
         if np.sum(inlets * outlets):
             raise Exception("Specified inlets and outlets overlap")
 
+    if dt is None:
+        dt = edt(im)
+
     if pc is None:
         pc = 2.0 / dt
     pc[~im] = 0  # Remove any infs or nans from pc computation
@@ -592,16 +578,31 @@ def drainage(
         Ps = np.unique(steps)  # To ensure they are in ascending order
 
     # Initialize empty arrays to accumulate results of each loop
-    nwp_mask = np.zeros_like(im, dtype=bool)
-    im_seq = np.zeros_like(im, dtype=int)
     im_pc = np.zeros_like(im, dtype=float)
     im_size = np.zeros_like(im, dtype=float)
+    im_seq = np.zeros_like(im, dtype=int)
+    trapped = np.zeros_like(im, dtype=bool)
+    if residual is not None:
+        im_seq[residual] = 1
+    if (outlets is not None) and (residual is not None):
+        trapped = find_trapped_clusters(
+            im=im,
+            seq=(im*2.0 - residual*1.0).astype(int),
+            outlets=outlets,
+            min_size=min_size,
+            method="labels" if len(Ps) < 100 else "queue",
+            conn=conn,
+        )
+    im_seq[trapped] = -1
+
+    # Initialize arrays used inside loop
+    nwp_mask = np.zeros_like(im, dtype=bool)
     seeds = np.zeros_like(im, dtype=bool)
 
     desc = inspect.currentframe().f_code.co_name  # Get current func name
-    for step, p in enumerate(tqdm(Ps, desc=desc, **settings.tqdm)):
+    for step, P in enumerate(tqdm(Ps, desc=desc, **settings.tqdm)):
         # Find all locations in image invadable at current pressure
-        invadable = (pc <= p) * im  # Equivalent to erosion
+        invadable = (pc <= P) * im  # Equivalent to erosion
         # Trim locations not connected to the inlets, if given
         if inlets is not None:
             invadable = trim_disconnected_voxels(
@@ -626,11 +627,20 @@ def drainage(
         # Deal with impact of residual, if present
         if residual is not None:
             if np.any(nwp_mask):
-                # Find invadable pixels connected to surviving residual
-                temp = trim_disconnected_voxels(residual, nwp_mask, conn=conn) * ~nwp_mask
+                # Find nwp pixels connected to residual
+                temp = trim_disconnected_voxels(
+                    im=residual,
+                    inlets=nwp_mask,
+                    conn=conn,
+                )
                 if np.any(temp):
                     # Trim invadable pixels not connected to residual
-                    invadable = trim_disconnected_voxels(invadable, temp, conn=conn)
+                    invadable = (pc <= P) * im  # Find full set of invadable again
+                    invadable = trim_disconnected_voxels(
+                        im=invadable,
+                        inlets=temp,
+                        conn=conn,
+                    )
                     coords = np.where(invadable)
                     radii = dt[coords].astype(int)
                     nwp_mask = _insert_disks_at_points_parallel(
@@ -641,10 +651,29 @@ def drainage(
                         smooth=True,
                         overwrite=False,
                     )
+        # This could be nested inside above if-block, but is outdented for clarity
+        if (residual is not None) and (outlets is not None):
+            nwp_mask = trim_disconnected_voxels(
+                im=nwp_mask * ~trapped,
+                inlets=inlets,
+                conn=conn,
+            )
+            trapped += find_trapped_clusters(
+                im=im,
+                seq=((~nwp_mask)*im*2.0 - residual*1.0).astype(int),
+                outlets=outlets,
+                min_size=min_size,
+                method="labels" if len(Ps) < 100 else "queue",
+                conn=conn,
+            )
+            trapped[residual] = False
+            nwp_mask[trapped] = False  # Set nwp in trapped regions to 0
+            im_seq[trapped] = -1
+
         mask = nwp_mask * (im_seq == 0) * im
         if np.any(mask):
             im_seq[mask] = step + 1
-            im_pc[mask] = p
+            im_pc[mask] = P
             if np.size(radii) > 0:
                 im_size[mask] = np.amin(radii)
         # Add new locations to list of invaded locations
@@ -653,41 +682,34 @@ def drainage(
     # Set uninvaded voxels to inf
     im_pc[(im_seq == 0) * im] = np.inf
 
-    # Add residual is given
+    # Add residual if given
     if residual is not None:
         im_pc[residual] = -np.inf
-        im_seq[residual] = 0
+        im_seq[residual] = 1
 
-    # Analyze trapping and adjust computed images accordingly
-    trapped = None  # Initialize trapped to None in case outlets not given
-    if outlets is not None:
+    # Analyze trapping as a post-processing step if no residual
+    if (outlets is not None) and (residual is None):
         trapped = find_trapped_clusters(
             im=im,
             seq=im_seq,
             outlets=outlets,
+            min_size=min_size,
             method="labels" if len(Ps) < 100 else "queue",
             conn=conn,
         )
-        if min_size > 0:
-            temp = find_small_clusters(
-                im=im,
-                trapped=trapped,
-                min_size=min_size,
-                conn=conn,
-            )
-            trapped = temp.im_trapped
         trapped[im_seq == -1] = True
         im_pc[trapped] = np.inf  # Trapped defender only displaced as Pc -> inf
-        if residual is not None:  # Re-add residual to inv
-            im_pc[residual] = -np.inf  # Residual defender is always present
 
-    # Initialize results object
+    im_seq = make_contiguous(im_seq, mode='symmetric')
+    satn = seq_to_satn(seq=im_seq, im=im, mode="drainage")
+
+    # Initialize results object to collect data
     results = Results()
-    results.im_snwp = pc_to_satn(pc=im_pc, im=im, mode="drainage")
+    results.im_snwp = satn
     results.im_seq = im_seq
-    # results.im_seq = pc_to_seq(pc=pc_inv, im=im, mode='drainage')
     results.im_pc = im_pc
     results.im_trapped = trapped
+
     if trapped is not None:
         results.im_seq[trapped] = -1
         results.im_snwp[trapped] = -1
@@ -695,19 +717,25 @@ def drainage(
     im_size[im_pc == np.inf] = np.inf
     im_size[im_pc == -np.inf] = -np.inf
     results.im_size = im_size
-    results.pc, results.snwp = pc_map_to_pc_curve(
+
+    pc_curve = pc_map_to_pc_curve(
         im=im,
         pc=results.im_pc,
         seq=results.im_seq,
         mode="drainage",
     )
+    results.pc = pc_curve.pc
+    results.snwp = pc_curve.snwp
     return results
 
 
+# %%
 if __name__ == "__main__":
     from copy import copy
 
     import matplotlib.pyplot as plt
+    from skimage.segmentation import clear_border
+    import scipy.ndimage as spim
 
     import porespy as ps
     ps.visualization.set_mpl_style()
@@ -717,24 +745,40 @@ def drainage(
     cm.set_over("k")
 
     # %% Run this cell to regenerate the variables in drainage
+    seed = np.random.randint(100000)  # 12129, 61227
     bg = "white"
     plots = True
-    im = ps.generators.blobs(
-        shape=[500, 500],
-        porosity=0.7,
-        blobiness=1.5,
-        seed=16,
+    # im = ps.generators.blobs(
+    #     shape=[500, 500],
+    #     porosity=0.7,
+    #     blobiness=1.5,
+    #     seed=seed,
+    # )
+    im = ~ps.generators.random_spheres(
+        [600, 600],
+        r=15,
+        clearance=15,
+        seed=1,
+        edges='extended',
+        phi=0.2,
     )
     im = ps.filters.fill_invalid_pores(im)
     inlets = np.zeros_like(im)
     inlets[0, :] = True
     outlets = np.zeros_like(im)
     outlets[-1, :] = True
-
-    lt = ps.filters.local_thickness(im)
     dt = edt(im)
-    residual = lt > 25
-    steps = 25
+
+    imb = ps.simulations.imbibition_dt(
+        im=im,
+        inlets=outlets,
+        outlets=inlets,
+    )
+
+    residual = imb.im_seq == -1
+    # residual = clear_border(spim.label(residual)[0]) > 0
+
+    steps = 50
     pc = ps.filters.capillary_transform(
         im=im,
         dt=dt,
@@ -751,7 +795,7 @@ def drainage(
         im=im,
         pc=pc,
         inlets=inlets,
-        steps=30,
+        steps=steps,
         min_size=5,
     )
     drn2 = ps.simulations.drainage(
@@ -759,19 +803,22 @@ def drainage(
         pc=pc,
         inlets=inlets,
         outlets=outlets,
-        steps=30,
+        steps=steps,
     )
     drn3 = ps.simulations.drainage(
         im=im,
         pc=pc,
         inlets=inlets,
         residual=residual,
-        steps=30,
+        steps=steps,
     )
     drn5 = ps.simulations.drainage(
         im=im,
         pc=pc,
-        steps=30,
+        steps=steps,
+        inlets=inlets,
+        outlets=outlets,
+        residual=residual,
     )
 
     # %% Visualize the invasion configurations for each scenario
@@ -782,13 +829,17 @@ def drainage(
         )
         # drn1.im_pc[~im] = -1
         ax["(a)"].imshow(drn1.im_seq / im, origin="lower", cmap=cm, vmin=0)
+        ax["(a)"].axis(False)
 
         vmax = drn2.im_seq.max()
         ax["(b)"].imshow(drn2.im_seq / im, origin="lower", cmap=cm, vmin=0, vmax=vmax)
+        ax["(b)"].axis(False)
 
         ax["(c)"].imshow(drn3.im_seq / im, origin="lower", cmap=cm, vmin=0)
+        ax["(c)"].axis(False)
 
         ax["(d)"].imshow(drn5.im_seq / im, origin="lower", cmap=cm, vmin=0)
+        ax["(d)"].axis(False)
 
         pc, s = ps.metrics.pc_map_to_pc_curve(
             pc=drn1.im_pc, seq=drn1.im_seq, im=im, mode="drainage"
@@ -808,6 +859,6 @@ def drainage(
         pc, s = ps.metrics.pc_map_to_pc_curve(
             pc=drn5.im_pc, seq=drn5.im_seq, im=im, mode="drainage"
         )
-        ax["(e)"].plot(np.log10(pc), s, "m-*", label="local thickness")
+        ax["(e)"].plot(np.log10(pc), s, "m-*", label="drainage w residual + trapping")
 
         ax["(e)"].legend()
diff --git a/src/porespy/simulations/_examples.py b/src/porespy/simulations/_examples.py
new file mode 100644
index 000000000..db625ebac
--- /dev/null
+++ b/src/porespy/simulations/_examples.py
@@ -0,0 +1,75 @@
+import numpy as np
+
+
+if __name__ == "__main__":
+    from copy import copy
+    import porespy as ps
+    from edt import edt
+    import matplotlib.pyplot as plt
+
+    ps.visualization.set_mpl_style()
+    cm = copy(plt.cm.plasma)
+    cm.set_under('grey')
+    # cm.set_over('grey')
+    cm.set_bad('k')
+    conn = 'min'
+    steps = 50
+
+    im = ~ps.generators.random_spheres(
+        [400, 400],
+        r=10,
+        # phi=0.15,
+        clearance=10,
+        seed=16,
+        edges="extended",
+        phi=0.25,
+    )
+    inlets = ps.generators.faces(shape=im.shape, inlet=0)
+    outlets = ps.generators.faces(shape=im.shape, outlet=0)
+    dt = edt(im)
+    pc = ps.filters.capillary_transform(
+        im=im,
+        dt=dt,
+        sigma=0.072,
+        theta=180,
+        g=0,
+        voxel_size=1e-5,
+    )
+
+    drn1 = ps.simulations.drainage(
+        im=im,
+        pc=pc,
+        inlets=inlets,
+        outlets=outlets,
+        steps=steps,
+        conn=conn,
+    )
+
+    imb2 = ps.simulations.imbibition(
+        im=im,
+        pc=pc,
+        inlets=outlets,
+        outlets=inlets,
+        residual=drn1.im_trapped,
+        steps=steps,
+        conn=conn,
+    )
+
+    drn2 = ps.simulations.drainage(
+        im=im,
+        pc=pc,
+        inlets=inlets,
+        outlets=outlets,
+        residual=imb2.im_trapped,
+        steps=steps,
+        conn=conn,
+    )
+
+# %%
+    tmp = np.copy(drn2.im_seq).astype(float)
+    # tmp[drn2.im_trapped] = tmp.max() + 1
+    tmp[~im] = np.nan
+
+    fig, ax = plt.subplots(figsize=[5, 5])
+    ax.imshow(tmp, cmap=cm, vmin=0, origin='lower')
+    ax.axis(False)
diff --git a/src/porespy/simulations/_imbibition.py b/src/porespy/simulations/_imbibition.py
index 2b47c719a..6b95a2f3f 100644
--- a/src/porespy/simulations/_imbibition.py
+++ b/src/porespy/simulations/_imbibition.py
@@ -33,11 +33,11 @@
     'imbibition_dt',
     'imbibition_dt_fft',
     'imbibition_fft',
-    'imbibition_dsi',
+    'imbibition_bf',
 ]
 
 
-def imbibition_dsi(
+def imbibition_bf(
     im,
     inlets=None,
     outlets=None,
@@ -492,6 +492,10 @@ def imbibition(
         An image the same shape as ``im`` with ``True`` values indicating the
         wetting fluid inlet(s).  If ``None`` then the wetting film is able to
         appear anywhere within the domain.
+    outlets : ndarray, optional
+        A boolean image with ``True`` values indicating the outlet locations.
+        If this is provided then trapped voxels of non-wetting phase are found and
+        all the output images are adjusted accordingly.
     residual : ndarray, optional
         A boolean mask the same shape as ``im`` with ``True`` values
         indicating to locations of residual wetting phase.
@@ -553,13 +557,16 @@ def imbibition(
     """
     im = np.array(im, dtype=bool)
 
+    if outlets is not None:
+        outlets = outlets * im
+        if np.sum(inlets * outlets):
+            raise Exception("Specified inlets and outlets overlap")
+
     if dt is None:
         dt = edt(im)
 
     if pc is None:
         pc = 2/dt
-
-    pc = np.copy(pc)
     pc[~im] = 0  # Remove any infs or nans from pc computation
 
     if isinstance(steps, int):
@@ -577,16 +584,30 @@ def imbibition(
     # Initialize empty arrays to accumulate results of each loop
     im_pc = np.zeros_like(im, dtype=float)
     im_seq = np.zeros_like(im, dtype=int)
+    trapped = np.zeros_like(im)
+    if residual is not None:
+        im_seq[residual] = 1
+    if (outlets is not None) and (residual is not None):
+        trapped = find_trapped_clusters(
+            im=im,
+            seq=(im*2.0 - residual*1.0).astype(int),
+            outlets=outlets,
+            min_size=min_size,
+            method="labels" if len(Ps) < 100 else "queue",
+            conn=conn,
+        )
+    im_seq[trapped] = -1
 
     desc = inspect.currentframe().f_code.co_name  # Get current func name
     for step, P in enumerate(tqdm(Ps, desc=desc, **settings.tqdm)):
         # This can be made faster if I find a way to get only seeds on edge, so
         # less spheres need to be drawn
-        invadable = (pc <= P)*im
+        invadable = (pc <= P)*im  # This means 'invadable by non-wetting phase'
         # Using FFT-based erosion to find edges.  When struct is small, this is
         # quite fast so it saves time overall by reducing the number of spheres
         # that need to be inserted.
         edges = (~erode(invadable, r=1, smooth=False, method='conv'))*invadable
+        # edges = invadable
         nwp_mask = np.zeros_like(im, dtype=bool)
         if np.any(edges):
             coords = np.where(edges)
@@ -606,50 +627,68 @@ def imbibition(
                 inlets=inlets,
                 conn=conn,
             )*im
-        if residual is not None:
-            nwp_mask = nwp_mask * ~residual
+
+        # Deal with impact of residual, if present
+        if (residual is not None) and (outlets is not None):
+                nwp_mask = ~trim_disconnected_voxels(
+                    im=~nwp_mask * ~trapped,
+                    inlets=inlets,
+                    conn=conn,
+                )
+                trapped += find_trapped_clusters(
+                    im=im,
+                    seq=(nwp_mask*im*2.0 - residual*1.0).astype(int),
+                    outlets=outlets,
+                    min_size=min_size,
+                    method="labels" if len(Ps) < 100 else "queue",
+                    conn=conn,
+                )
+                trapped[residual] = False
+                nwp_mask[trapped] = True
+                im_seq[trapped] = -1
 
         mask = (nwp_mask == 0) * (im_seq == 0) * im
         if np.any(mask):
-            im_seq[mask] = step
+            im_seq[mask] = step + 1
             im_pc[mask] = P
-    im_seq = make_contiguous(im_seq)
 
-    trapped = None  # Initialize trapped to None in case outlets not given
-    if outlets is not None:
-        if inlets is not None:
-            outlets[inlets] = False  # Ensure outlets do not overlap inlets
+    # Set uninvaded voxels to -inf
+    im_pc[(im_seq == 0) * im] = -np.inf
+
+    # Add residual if given
+    if residual is not None:
+        im_pc[residual] = np.inf
+        im_seq[residual] = 1
+
+    # Check for trapping as a post-processing step if no residual
+    if (outlets is not None) and (residual is None):
         trapped = find_trapped_clusters(
             im=im,
             seq=im_seq,
             outlets=outlets,
+            min_size=min_size,
             method='labels' if len(Ps) < 100 else 'queue',
             conn=conn,
         )
-        if min_size > 0:
-            temp = find_small_clusters(
-                im=im,
-                trapped=trapped,
-                min_size=min_size,
-                conn=conn,
-            )
-            trapped = temp.im_trapped
+        trapped[im_seq == -1] = True
         im_pc[trapped] = -np.inf
-        im_seq[trapped] = -1
-
-    if residual is not None:
-        im_pc[residual] = np.inf
-        im_seq[residual] = 0
 
+    im_seq = make_contiguous(im_seq, mode='symmetric')
     satn = seq_to_satn(im=im, seq=im_seq, mode='imbibition')
-    # Collect data in a Results object
+
+    # Initialize Results object to collect data
     results = Results()
     results.im_snwp = satn
     results.im_seq = im_seq
     results.im_pc = im_pc
     results.im_trapped = trapped
 
-    pc_curve = pc_map_to_pc_curve(pc=im_pc, im=im, seq=im_seq, mode='imbibition')
+    pc_curve = pc_map_to_pc_curve(
+        pc=im_pc,
+        im=im,
+        seq=im_seq,
+        mode='imbibition',
+    )
     results.pc = pc_curve.pc
     results.snwp = pc_curve.snwp
     return results
@@ -710,25 +749,27 @@ def _insert_disks_npoints_nradii_1value_parallel(
     cm = copy(plt.cm.turbo)
     cm.set_under('k')
     cm.set_over('grey')
+    steps = 50
 
     i = np.random.randint(1, 100000)  # bad: 38364, good: 65270, 71698
     i = 50591
     # i = 59477  # Bug in pc curve if lowest point is not 0.99 x min(pc)
     # i = 38364
     print(i)
-    im = ps.generators.blobs([500, 500], porosity=0.65, blobiness=2, seed=i)
+    im = ps.generators.blobs([500, 500], porosity=0.8, blobiness=1.5, seed=i)
     im = ps.filters.fill_invalid_pores(im)
 
     inlets = ps.generators.faces(im.shape, inlet=0)
     outlets = ps.generators.faces(im.shape, outlet=0)
-    lt = ps.filters.local_thickness_dt(im)
-    residual = (lt < 8)*im
     pc = ps.filters.capillary_transform(im=im, voxel_size=1e-4)
 
-    imb1 = imbibition(im=im, pc=pc, inlets=inlets, min_size=1)
-    imb2 = imbibition(im=im, pc=pc, inlets=inlets, outlets=outlets, min_size=1)
-    imb3 = imbibition(im=im, pc=pc, inlets=inlets, residual=residual, min_size=1)
-    imb4 = imbibition(im=im, pc=pc, inlets=inlets, outlets=outlets, residual=residual, min_size=1)
+    drn = ps.simulations.drainage(im=im, pc=pc, inlets=inlets, outlets=outlets, steps=steps)
+    residual = drn.im_trapped
+
+    imb1 = imbibition(im=im, pc=pc, inlets=inlets, steps=steps, min_size=5)
+    imb2 = imbibition(im=im, pc=pc, inlets=inlets, outlets=outlets, steps=steps, min_size=5)
+    imb3 = imbibition(im=im, pc=pc, inlets=inlets, residual=residual, steps=steps, min_size=5)
+    imb4 = imbibition(im=im, pc=pc, inlets=inlets, outlets=outlets, residual=residual, steps=steps, min_size=5)
 
     # %%
 
diff --git a/test/unit/test_displacement_refs.py b/test/unit/test_displacement_refs.py
index 32662a0ef..3c10b2ef4 100644
--- a/test/unit/test_displacement_refs.py
+++ b/test/unit/test_displacement_refs.py
@@ -29,7 +29,7 @@ def test_drainage_2D_no_trapping_smooth(self):
             inlets=self.inlets,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im,
             inlets=self.inlets,
             smooth=smooth,
@@ -40,7 +40,7 @@ def test_drainage_2D_no_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -71,7 +71,7 @@ def test_drainage_2D_no_trapping_not_smooth(self):
             inlets=self.inlets,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im,
             inlets=self.inlets,
             smooth=smooth,
@@ -82,7 +82,7 @@ def test_drainage_2D_no_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -113,7 +113,7 @@ def test_drainage_3D_no_trapping_not_smooth(self):
             inlets=self.inlets3D,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             smooth=smooth,
@@ -124,7 +124,7 @@ def test_drainage_3D_no_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -155,7 +155,7 @@ def test_drainage_3D_no_trapping_smooth(self):
             inlets=self.inlets3D,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             smooth=smooth,
@@ -166,7 +166,7 @@ def test_drainage_3D_no_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -199,7 +199,7 @@ def test_drainage_2D_w_trapping_smooth(self):
             outlets=self.outlets,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im,
             inlets=self.inlets,
             outlets=self.outlets,
@@ -212,7 +212,7 @@ def test_drainage_2D_w_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -245,7 +245,7 @@ def test_drainage_2D_w_trapping_not_smooth(self):
             outlets=self.outlets,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im,
             inlets=self.inlets,
             outlets=self.outlets,
@@ -258,7 +258,7 @@ def test_drainage_2D_w_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -291,7 +291,7 @@ def test_drainage_3D_w_trapping_not_smooth(self):
             outlets=self.outlets3D,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             outlets=self.outlets3D,
@@ -304,7 +304,7 @@ def test_drainage_3D_w_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -337,7 +337,7 @@ def test_drainage_3D_w_trapping_smooth(self):
             outlets=self.outlets3D,
             smooth=smooth,
         )
-        drn['dsi'] = ps.simulations.drainage_dsi(
+        drn['bf'] = ps.simulations.drainage_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             outlets=self.outlets3D,
@@ -350,7 +350,7 @@ def test_drainage_3D_w_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = drn[a].im_seq == drn[b].im_seq
         assert np.all(tmp)
         tmp = drn[a].im_size == drn[b].im_size
@@ -371,7 +371,7 @@ def test_drainage_3D_w_trapping_smooth(self):
     def test_imbibition_2D_no_trapping_smooth(self):
         imb = {}
         smooth = True
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im,
             inlets=self.inlets,
             smooth=smooth,
@@ -392,7 +392,7 @@ def test_imbibition_2D_no_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size
@@ -413,7 +413,7 @@ def test_imbibition_2D_no_trapping_smooth(self):
     def test_imbibition_2D_no_trapping_not_smooth(self):
         imb = {}
         smooth = False
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im,
             inlets=self.inlets,
             smooth=smooth,
@@ -434,7 +434,7 @@ def test_imbibition_2D_no_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size
@@ -455,7 +455,7 @@ def test_imbibition_2D_no_trapping_not_smooth(self):
     def test_imbibition_3D_no_trapping_smooth(self):
         imb = {}
         smooth = True
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             smooth=smooth,
@@ -476,7 +476,7 @@ def test_imbibition_3D_no_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size
@@ -497,7 +497,7 @@ def test_imbibition_3D_no_trapping_smooth(self):
     def test_imbibition_3D_no_trapping_not_smooth(self):
         imb = {}
         smooth = False
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             smooth=smooth,
@@ -518,7 +518,7 @@ def test_imbibition_3D_no_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size
@@ -539,7 +539,7 @@ def test_imbibition_3D_no_trapping_not_smooth(self):
     def test_imbibition_2D_w_trapping_smooth(self):
         imb = {}
         smooth = True
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im,
             inlets=self.inlets,
             outlets=self.outlets,
@@ -564,7 +564,7 @@ def test_imbibition_2D_w_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size
@@ -585,7 +585,7 @@ def test_imbibition_2D_w_trapping_smooth(self):
     def test_imbibition_2D_w_trapping_not_smooth(self):
         imb = {}
         smooth = False
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im,
             inlets=self.inlets,
             outlets=self.outlets,
@@ -610,7 +610,7 @@ def test_imbibition_2D_w_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size
@@ -631,7 +631,7 @@ def test_imbibition_2D_w_trapping_not_smooth(self):
     def test_imbibition_3D_w_trapping_smooth(self):
         imb = {}
         smooth = True
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             outlets=self.outlets3D,
@@ -656,7 +656,7 @@ def test_imbibition_3D_w_trapping_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size
@@ -677,7 +677,7 @@ def test_imbibition_3D_w_trapping_smooth(self):
     def test_imbibition_3D_w_trapping_not_smooth(self):
         imb = {}
         smooth = False
-        imb['dsi'] = ps.simulations.imbibition_dsi(
+        imb['bf'] = ps.simulations.imbibition_bf(
             im=self.im3D,
             inlets=self.inlets3D,
             outlets=self.outlets3D,
@@ -702,7 +702,7 @@ def test_imbibition_3D_w_trapping_not_smooth(self):
             smooth=smooth,
         )
 
-        a, b = 'dt', 'dsi'
+        a, b = 'dt', 'bf'
         tmp = imb[a].im_seq == imb[b].im_seq
         assert np.all(tmp)
         tmp = imb[a].im_size == imb[b].im_size