diff --git a/.docs/Notebooks/horizontal_flow_barriers_from_linework.py b/.docs/Notebooks/horizontal_flow_barriers_from_linework.py
new file mode 100644
index 000000000..6db2ed3cc
--- /dev/null
+++ b/.docs/Notebooks/horizontal_flow_barriers_from_linework.py
@@ -0,0 +1,305 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# # Creating HFB input data from shapefile linework
+#
+# This notebook shows examples of how to create a horizontal flow barrier package (HFB) from shapefile information in FloPy. FloPy has support for creating HFB inputs from LineString like vector data for all discretization types (`DIS`, `DISV`, `DISU`).
+#
+# This notebook starts off by generating two simple models:
+#    - a structured (DIS) model
+#    - a vertex (DISV) model
+#
+# And then shows how to generate HFB input for each of these models.
+
+# In[1]:
+
+
+from tempfile import TemporaryDirectory
+
+import geopandas as gpd
+import matplotlib.pyplot as plt
+import numpy as np
+from shapely.geometry import LineString
+
+import flopy
+from flopy.utils import make_hfb_array
+
+temp_dir = TemporaryDirectory()
+workspace = temp_dir.name
+
+
+# ### Generate a Structured (DIS) model
+
+# In[2]:
+
+
+def simple_structured_model():
+    """
+    Method to generate a 10x10 structured grid model
+    """
+    lx = 100
+    ly = 100
+    nlay = 1
+    nrow = 10
+    ncol = 10
+    delc = np.full((nrow,), ly / nrow)
+    delr = np.full((ncol,), ly / ncol)
+    top = np.full((nrow, ncol), 10)
+    botm = np.zeros((nlay, nrow, ncol))
+    idomain = np.ones(botm.shape, dtype=int)
+
+    sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct")
+    ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model")
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delc=delc,
+        delr=delr,
+        top=top,
+        botm=botm,
+        idomain=idomain,
+    )
+    npf = flopy.mf6.ModflowGwfnpf(gwf)
+    sto = flopy.mf6.ModflowGwfsto(gwf)
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=np.full((nlay, nrow, ncol), 9))
+
+    return sim
+
+
+# ### Generate a voronoi (DISV) model
+
+# In[3]:
+
+
+def simple_vertex_model(workspace):
+    """
+    Method to generate a voroni vertex grid model
+    """
+    from flopy.utils.triangle import Triangle
+    from flopy.utils.voronoi import VoronoiGrid
+
+    geom = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)]
+    tri = Triangle(angle=30, model_ws=workspace)
+    tri.add_polygon(geom)
+    tri.add_region((5, 5), 0, maximum_area=40)
+    tri.build()
+
+    vor = VoronoiGrid(tri)
+    pkg_props = vor.get_disv_gridprops()
+    ncpl = pkg_props["ncpl"]
+    nlay = 1
+    top = np.full((ncpl,), 10)
+    botm = np.zeros((nlay, ncpl))
+    idomain = np.ones(botm.shape, dtype=int)
+
+    sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct")
+    ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model")
+    disv = flopy.mf6.ModflowGwfdisv(
+        gwf, nlay=1, top=top, botm=botm, idomain=idomain, **pkg_props
+    )
+    npf = flopy.mf6.ModflowGwfnpf(gwf)
+    sto = flopy.mf6.ModflowGwfsto(gwf)
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=np.full((nlay, ncpl), 9))
+
+    return sim
+
+
+# ### Create HFB inputs for a Structured (DIS) model
+
+# The `make_hfb_array` utility is used to intersect `LineString` like data with `modelgrid` instances to produce HFB data. The `make_hfb_array` utility accepts two parameters:
+#    - `modelgrid` : a flopy Grid instance
+#    - `geom` : a geospatial LineString like object that can be: an iterable of points, GeoJSON, Shapely.geometry.LineString, or a shapefile.Shape object.
+
+# First define a fault line and visualize it with the model
+
+# In[4]:
+
+
+# define a fault line
+fault = [(0, 4), (50, 55), (100, 55)]
+
+# build simulation
+sim = simple_structured_model()
+gwf = sim.get_model()
+modelgrid = gwf.modelgrid
+
+fig, ax = plt.subplots(figsize=(5, 5))
+pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax)
+pmv.plot_grid()
+plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--")
+# now create hfb data with the `make_hfb_array` utility
+
+# In[5]:
+
+
+hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+hfbs[0:5]
+
+
+# notice that the `hydchr` parameter is not filled out; this is left up to the user
+
+# In[6]:
+
+
+hfbs["hydchr"] = 1e-05
+
+
+# Now we can build an HFB package and compare it to the fault line that was defined above
+
+# In[7]:
+
+
+hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfbs})
+
+
+# In[8]:
+
+
+# plot up the results
+fig, ax = plt.subplots(figsize=(5, 5))
+pmv = flopy.plot.PlotMapView(modelgrid=modelgrid)
+pmv.plot_grid()
+pmv.plot_bc(package=hfb, color="b", lw=2)
+plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--")
+# ### Create HFB inputs for a Vertex (DISV) model
+
+# The `make_hfb_array` utility is used to intersect `LineString` like data with `modelgrid` instances to produce HFB data. The `make_hfb_array` utility accepts two parameters:
+#    - `modelgrid` : a flopy Grid instance
+#    - `geom` : a geospatial LineString like object that can be: an iterable of points, GeoJSON, Shapely.geometry.LineString, or a shapefile.Shape object.
+
+# First define a fault line and visualize it with the model
+
+# In[9]:
+
+
+# define a fault line
+fault = [(0, 4), (50, 55), (100, 55)]
+
+# build simulation
+sim = simple_vertex_model(workspace)
+gwf = sim.get_model()
+modelgrid = gwf.modelgrid
+
+fig, ax = plt.subplots(figsize=(5, 5))
+pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax)
+pmv.plot_grid()
+plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--")
+# now create hfb data with the `make_hfb_array` utility
+
+# In[10]:
+
+
+hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+hfbs[0:5]
+
+
+# notice that the `hydchr` parameter is not filled out; this is left up to the user
+
+# In[11]:
+
+
+hfbs["hydchr"] = 3e-06
+
+
+# Now we can build an HFB package and compare it to the fault line that was defined above
+
+# In[12]:
+
+
+hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfbs})
+
+
+# In[13]:
+
+
+# plot up the results
+fig, ax = plt.subplots(figsize=(5, 5))
+pmv = flopy.plot.PlotMapView(modelgrid=modelgrid)
+pmv.plot_grid()
+pmv.plot_bc(package=hfb, color="b", lw=2)
+plt.plot(np.array(fault).T[0], np.array(fault).T[1], "r--")
+# ### Working with multiple fault geometries
+#
+# The `make_hfb_array` only accepts a single LineString geometry at a time. This was done by design for a couple reasons.
+#    1) This gives the user an opportunity to set the `hydchr` parameter to a unique value for each fault line
+#    2) This also gives the user an opportunity to save the boundary condition data (cellid1, cellid2) to an external file for each fault line that can be useful during model calibration with automated software (e.g., PEST++)
+#
+# As a result of these design descisions recarray data will need to be concatenated when working with multiple LineString geometries.
+
+# For this example, a geodataframe of two fault lines is created and then used to build a HFB package
+
+# In[14]:
+
+
+geom1 = LineString([(0, 4), (50, 44)])
+geom2 = LineString([(50, 44), (88, 100)])
+d = {"name": ["fault1", "fault2"], "geometry": [geom1, geom2]}
+gdf = gpd.GeoDataFrame(d)
+gdf
+
+
+# Visualize the faults on the modelgrid
+
+# In[15]:
+
+
+# build simulation
+sim = simple_structured_model()
+gwf = sim.get_model()
+modelgrid = gwf.modelgrid
+
+fig, ax = plt.subplots(figsize=(5, 5))
+pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax)
+pmv.plot_grid()
+gdf.geometry.plot(color="r", linestyle="--", ax=ax)
+# now create hfb data with the `make_hfb_array` utility and concatenate it
+
+# In[16]:
+
+
+hydchrs = {"fault1": 1e-05, "fault2": 1e-06}
+
+hfb_ras = []
+for name, geom in zip(gdf.name, gdf.geometry):
+    hfbs = make_hfb_array(modelgrid, geom)
+    hfbs["hydchr"] = hydchrs[name]
+    hfb_ras.append(hfbs)
+
+hfbs = np.concatenate(hfb_ras)
+hfbs = hfbs.view(np.recarray)
+hfbs
+
+
+# Now create the hfb package from this data and visualize it for comparison with the fault lines
+
+# In[17]:
+
+
+hfb = flopy.mf6.ModflowGwfhfb(gwf, stress_period_data={0: hfbs})
+
+
+# In[18]:
+
+
+# plot up the results
+fig, ax = plt.subplots(figsize=(5, 5))
+pmv = flopy.plot.PlotMapView(modelgrid=modelgrid)
+pmv.plot_grid()
+pmv.plot_bc(package=hfb, color="b", lw=2)
+gdf.geometry.plot(color="r", linestyle="--", ax=ax)
+# In[19]:
+
+
+try:
+    temp_dir.cleanup()
+except PermissionError:
+    # can occur on windows: https://docs.python.org/3/library/tempfile.html#tempfile.TemporaryDirectory
+    pass
diff --git a/autotest/test_hfb_util.py b/autotest/test_hfb_util.py
new file mode 100644
index 000000000..31f9958df
--- /dev/null
+++ b/autotest/test_hfb_util.py
@@ -0,0 +1,439 @@
+import numpy as np
+from modflow_devtools.markers import requires_exe
+
+import flopy
+from flopy.utils.hfb_util import make_hfb_array
+from flopy.utils.triangle import Triangle
+from flopy.utils.voronoi import VoronoiGrid
+
+
+def structured_sim():
+    lx = 100
+    ly = 100
+    nlay = 1
+    nrow = 10
+    ncol = 10
+    delc = np.full((nrow,), ly / nrow)
+    delr = np.full((ncol,), lx / ncol)
+    top = np.full((nrow, ncol), 10)
+    botm = np.zeros((nlay, nrow, ncol))
+    idomain = np.ones(botm.shape, dtype=int)
+
+    sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct")
+    ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model")
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delc=delc,
+        delr=delr,
+        top=top,
+        botm=botm,
+        idomain=idomain,
+    )
+
+    return sim
+
+
+def vertex_sim(path):
+    geom = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)]
+    tri = Triangle(angle=30, model_ws=path)
+    tri.add_polygon(geom)
+    tri.add_region((5, 5), 0, maximum_area=40)
+    tri.build()
+
+    vor = VoronoiGrid(tri)
+    pkg_props = vor.get_disv_gridprops()
+    ncpl = pkg_props["ncpl"]
+    nlay = 1
+    top = np.full((ncpl,), 10)
+    botm = np.zeros((nlay, ncpl))
+    idomain = np.ones(botm.shape, dtype=int)
+
+    sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct")
+    ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model")
+    disv = flopy.mf6.ModflowGwfdisv(
+        gwf, nlay=1, top=top, botm=botm, idomain=idomain, **pkg_props
+    )
+
+    return sim
+
+
+def unstructured_sim(path):
+    geom = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)]
+    tri = Triangle(angle=30, model_ws=path)
+    tri.add_polygon(geom)
+    tri.add_region((5, 5), 0, maximum_area=40)
+    tri.build()
+
+    vor = VoronoiGrid(tri)
+    grid_props = vor.get_disu6_gridprops()
+
+    top = np.full((grid_props["nodes"],), 10)
+    botm = np.zeros((grid_props["nodes"]))
+    idomain = np.ones(botm.shape, dtype=int)
+
+    sim = flopy.mf6.MFSimulation(sim_ws="tmp_struct")
+    ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")
+    tdis = flopy.mf6.ModflowTdis(sim)
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="hfb_model")
+    disu = flopy.mf6.ModflowGwfdisu(
+        gwf, top=top, bot=botm, idomain=idomain, **grid_props
+    )
+
+    return sim
+
+
+def test_simple_structured():
+    validation = [
+        (0, 0, 8, 0, 1, 8),
+        (0, 0, 8, 0, 0, 9),
+        (0, 1, 7, 0, 2, 7),
+        (0, 1, 7, 0, 1, 8),
+        (0, 2, 6, 0, 3, 6),
+        (0, 2, 6, 0, 2, 7),
+        (0, 3, 5, 0, 4, 5),
+        (0, 3, 5, 0, 3, 6),
+        (0, 4, 4, 0, 4, 5),
+        (0, 4, 4, 0, 5, 4),
+        (0, 5, 3, 0, 5, 4),
+        (0, 5, 3, 0, 6, 3),
+        (0, 6, 2, 0, 6, 3),
+        (0, 6, 2, 0, 7, 2),
+        (0, 7, 1, 0, 7, 2),
+        (0, 7, 1, 0, 8, 1),
+        (0, 8, 0, 0, 8, 1),
+        (0, 8, 0, 0, 9, 0),
+    ]
+
+    fault = [(0, 4), (94.0, 100)]
+    sim = structured_sim()
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+    for row in hfbs:
+        cid1 = row.cellid1
+        cid2 = row.cellid2
+        x = sorted([cid1, cid2])
+        test_set = x[0] + x[1]
+        if test_set not in validation:
+            raise AssertionError(
+                f"HFB Line {x[0]} and {x[1]} are outside of validation data set"
+            )
+
+    if len(hfbs) != len(validation):
+        raise AssertionError(
+            f"HFB data length {len(hfbs)} not equal to validation {len(validation)}"
+        )
+
+
+@requires_exe("triangle")
+def test_simple_vertex(function_tmpdir):
+    validation = [
+        (0, 53, 0, 92),
+        (0, 65, 0, 162),
+        (0, 64, 0, 65),
+        (0, 69, 0, 71),
+        (0, 64, 0, 72),
+        (0, 45, 0, 77),
+        (0, 53, 0, 77),
+        (0, 78, 0, 79),
+        (0, 76, 0, 78),
+        (0, 78, 0, 162),
+        (0, 79, 0, 82),
+        (0, 11, 0, 79),
+        (0, 69, 0, 80),
+        (0, 45, 0, 80),
+        (0, 11, 0, 84),
+        (0, 69, 0, 84),
+        (0, 94, 0, 114),
+        (0, 4, 0, 94),
+        (0, 4, 0, 96),
+        (0, 97, 0, 101),
+        (0, 97, 0, 99),
+        (0, 98, 0, 99),
+        (0, 96, 0, 98),
+        (0, 108, 0, 130),
+        (0, 108, 0, 109),
+        (0, 101, 0, 108),
+        (0, 92, 0, 114),
+        (0, 107, 0, 130),
+        (0, 125, 0, 130),
+        (0, 125, 0, 132),
+        (0, 149, 0, 171),
+        (0, 133, 0, 149),
+        (0, 150, 0, 151),
+        (0, 151, 0, 171),
+        (0, 150, 0, 154),
+        (0, 125, 0, 169),
+        (0, 119, 0, 169),
+        (0, 149, 0, 169),
+    ]
+
+    fault = [(0, 4), (94.0, 100)]
+    sim = vertex_sim(function_tmpdir)
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+    for row in hfbs:
+        cid1 = row.cellid1
+        cid2 = row.cellid2
+        x = sorted([cid1, cid2])
+        test_set = x[0] + x[1]
+        if test_set not in validation:
+            raise AssertionError(
+                f"HFB Line {x[0]} and {x[1]} are outside of validation data set"
+            )
+
+    if len(hfbs) != len(validation):
+        raise AssertionError(
+            f"HFB data length {len(hfbs)} not equal to validation {len(validation)}"
+        )
+
+
+@requires_exe("triangle")
+def test_simple_unstructured(function_tmpdir):
+    validation = [
+        (53, 92),
+        (65, 162),
+        (64, 65),
+        (69, 71),
+        (64, 72),
+        (45, 77),
+        (53, 77),
+        (78, 79),
+        (76, 78),
+        (78, 162),
+        (79, 82),
+        (11, 79),
+        (69, 80),
+        (45, 80),
+        (11, 84),
+        (69, 84),
+        (94, 114),
+        (4, 94),
+        (4, 96),
+        (97, 101),
+        (97, 99),
+        (98, 99),
+        (96, 98),
+        (108, 130),
+        (108, 109),
+        (101, 108),
+        (92, 114),
+        (107, 130),
+        (125, 130),
+        (125, 132),
+        (149, 171),
+        (133, 149),
+        (150, 151),
+        (151, 171),
+        (150, 154),
+        (125, 169),
+        (119, 169),
+        (149, 169),
+    ]
+    fault = [(0, 4), (94.0, 100)]
+    sim = unstructured_sim(function_tmpdir)
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+    for row in hfbs:
+        cid1 = row.cellid1
+        cid2 = row.cellid2
+        x = sorted([cid1, cid2])
+        test_set = x[0] + x[1]
+        if test_set not in validation:
+            raise AssertionError(
+                f"HFB Line {x[0]} and {x[1]} are outside of validation data set"
+            )
+
+    if len(hfbs) != len(validation):
+        raise AssertionError(
+            f"HFB data length {len(hfbs)} not equal to validation {len(validation)}"
+        )
+
+
+def test_multisegment_structured():
+    validation = [
+        (0, 0, 8, 0, 1, 8),
+        (0, 0, 8, 0, 0, 9),
+        (0, 1, 7, 0, 1, 8),
+        (0, 1, 7, 0, 2, 7),
+        (0, 2, 6, 0, 2, 7),
+        (0, 3, 6, 0, 4, 6),
+        (0, 3, 6, 0, 3, 7),
+        (0, 4, 5, 0, 4, 6),
+        (0, 4, 5, 0, 5, 5),
+        (0, 5, 4, 0, 5, 5),
+        (0, 6, 3, 0, 7, 3),
+        (0, 6, 3, 0, 6, 4),
+        (0, 5, 4, 0, 6, 4),
+        (0, 7, 1, 0, 8, 1),
+        (0, 6, 2, 0, 7, 2),
+        (0, 7, 1, 0, 7, 2),
+        (0, 8, 0, 0, 9, 0),
+        (0, 8, 0, 0, 8, 1),
+    ]
+
+    fault = [[0, 11], [55, 45], [89, 100]]
+    sim = structured_sim()
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+    for row in hfbs:
+        cid1 = row.cellid1
+        cid2 = row.cellid2
+        x = sorted([cid1, cid2])
+        test_set = x[0] + x[1]
+        if test_set not in validation:
+            raise AssertionError(
+                f"HFB Line {x[0]} and {x[1]} are outside of validation data set"
+            )
+
+    if len(hfbs) != len(validation):
+        raise AssertionError(
+            f"HFB data length {len(hfbs)} not equal to validation {len(validation)}"
+        )
+
+
+@requires_exe("triangle")
+def test_multisegment_vertex(function_tmpdir):
+    validation = [
+        (0, 11, 0, 84),
+        (0, 11, 0, 79),
+        (0, 40, 0, 103),
+        (0, 4, 0, 40),
+        (0, 40, 0, 114),
+        (0, 40, 0, 53),
+        (0, 45, 0, 53),
+        (0, 45, 0, 77),
+        (0, 45, 0, 80),
+        (0, 29, 0, 53),
+        (0, 64, 0, 72),
+        (0, 64, 0, 65),
+        (0, 69, 0, 80),
+        (0, 69, 0, 71),
+        (0, 69, 0, 84),
+        (0, 78, 0, 79),
+        (0, 76, 0, 78),
+        (0, 78, 0, 162),
+        (0, 79, 0, 82),
+        (0, 97, 0, 110),
+        (0, 97, 0, 102),
+        (0, 98, 0, 102),
+        (0, 98, 0, 103),
+        (0, 4, 0, 103),
+        (0, 107, 0, 130),
+        (0, 107, 0, 108),
+        (0, 108, 0, 110),
+        (0, 125, 0, 132),
+        (0, 125, 0, 130),
+        (0, 133, 0, 149),
+        (0, 150, 0, 151),
+        (0, 150, 0, 154),
+        (0, 65, 0, 162),
+        (0, 125, 0, 169),
+        (0, 119, 0, 169),
+        (0, 149, 0, 169),
+        (0, 149, 0, 171),
+        (0, 151, 0, 171),
+    ]
+
+    fault = [[0, 11], [55, 45], [89, 100]]
+    sim = vertex_sim(function_tmpdir)
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+    for row in hfbs:
+        cid1 = row.cellid1
+        cid2 = row.cellid2
+        x = sorted([cid1, cid2])
+        test_set = x[0] + x[1]
+        if test_set not in validation:
+            raise AssertionError(
+                f"HFB Line {x[0]} and {x[1]} are outside of validation data set"
+            )
+
+    if len(hfbs) != len(validation):
+        raise AssertionError(
+            f"HFB data length {len(hfbs)} not equal to validation {len(validation)}"
+        )
+
+
+def test_colinear_hfb():
+    validation = [
+        (0, 0, 4, 0, 0, 5),
+        (0, 1, 4, 0, 1, 5),
+        (0, 2, 4, 0, 2, 5),
+        (0, 3, 4, 0, 3, 5),
+        (0, 4, 4, 0, 4, 5),
+        (0, 5, 4, 0, 5, 5),
+        (0, 6, 4, 0, 6, 5),
+        (0, 7, 4, 0, 7, 5),
+        (0, 8, 4, 0, 8, 5),
+        (0, 9, 4, 0, 9, 5),
+    ]
+
+    fault = [(50, 0), (50, 100)]
+    sim = structured_sim()
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+    for row in hfbs:
+        cid1 = row.cellid1
+        cid2 = row.cellid2
+        x = sorted([cid1, cid2])
+        test_set = x[0] + x[1]
+        if test_set not in validation:
+            raise AssertionError(
+                f"HFB Line {x[0]} and {x[1]} are outside of validation data set"
+            )
+
+    if len(hfbs) != len(validation):
+        raise AssertionError(
+            f"HFB data length {len(hfbs)} not equal to validation {len(validation)}"
+        )
+
+
+def test_split_cell_hfb():
+    validation = [
+        (0, 4, 0, 0, 5, 0),
+        (0, 4, 1, 0, 5, 1),
+        (0, 4, 2, 0, 5, 2),
+        (0, 4, 3, 0, 5, 3),
+        (0, 4, 4, 0, 5, 4),
+        (0, 4, 5, 0, 5, 5),
+        (0, 4, 6, 0, 5, 6),
+        (0, 4, 7, 0, 5, 7),
+        (0, 4, 8, 0, 5, 8),
+        (0, 4, 9, 0, 5, 9),
+    ]
+
+    fault = [(0, 55), (100, 55)]
+    sim = structured_sim()
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    hfbs = flopy.utils.make_hfb_array(modelgrid, fault)
+    for row in hfbs:
+        cid1 = row.cellid1
+        cid2 = row.cellid2
+        x = sorted([cid1, cid2])
+        test_set = x[0] + x[1]
+        if test_set not in validation:
+            raise AssertionError(
+                f"HFB Line {x[0]} and {x[1]} are outside of validation data set"
+            )
+
+    if len(hfbs) != len(validation):
+        raise AssertionError(
+            f"HFB data length {len(hfbs)} not equal to validation {len(validation)}"
+        )
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 4c0b74e50..0a50597c9 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -411,6 +411,29 @@ def laycbd(self):
         else:
             return self._laycbd
 
+    @property
+    def cell_area(self):
+        """
+        Use shoelace algorithm for non-self-intersecting polygons to
+        calculate area.
+
+        Returns
+        -------
+            area : np.ndarray
+                numpy array of cell areas in L^2
+        """
+        from ..plot.plotutil import UnstructuredPlotUtilities
+
+        xverts, yverts = self.cross_section_vertices
+        xverts, yverts = UnstructuredPlotUtilities.irregular_shape_patch(xverts, yverts)
+        area_x2 = np.zeros((1, len(xverts)))
+        for i in range(xverts.shape[-1]):
+            # calculate the determinant of each line in polygon
+            area_x2 += xverts[:, i - 1] * yverts[:, i] - yverts[:, i - 1] * xverts[:, i]
+
+        area = np.abs(area_x2 / 2.0)
+        return np.ravel(area)
+
     @property
     def cell_thickness(self):
         """
diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py
index bd15a59b8..0da0e3d10 100644
--- a/flopy/utils/__init__.py
+++ b/flopy/utils/__init__.py
@@ -30,6 +30,7 @@
 
 get_modflow = get_modflow_module.run_main
 from .gridintersect import GridIntersect
+from .hfb_util import make_hfb_array
 from .mflistfile import (
     Mf6ListBudget,
     MfListBudget,
diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py
index f8c9f46c8..1304fe41d 100644
--- a/flopy/utils/gridintersect.py
+++ b/flopy/utils/gridintersect.py
@@ -106,6 +106,8 @@ def __init__(self, mfgrid, rtree=True, local=False):
             self.geoms, self.cellids = self._rect_grid_to_geoms_cellids()
         elif self.mfgrid.grid_type == "vertex":
             self.geoms, self.cellids = self._vtx_grid_to_geoms_cellids()
+        elif self.mfgrid.grid_type == "unstructured":
+            self.geoms, self.cellids = self._usg_grid_to_geoms_cellids()
         else:
             raise NotImplementedError(
                 f"Grid type {self.mfgrid.grid_type} not supported"
@@ -375,7 +377,29 @@ def _usg_grid_to_geoms_cellids(self):
         cellids : array_like
             array of cellids
         """
-        raise NotImplementedError()
+        warnings.warn(
+            "UnstructuredGrid intersection is experimental", category=UserWarning
+        )
+        shapely = import_optional_dependency("shapely")
+        if self.local:
+            geoms = [
+                shapely.polygons(
+                    list(
+                        zip(
+                            *self.mfgrid.get_local_coords(
+                                *np.array(self.mfgrid.get_cell_vertices(node)).T
+                            )
+                        )
+                    )
+                )
+                for node in range(self.mfgrid.nnodes)
+            ]
+        else:
+            geoms = [
+                shapely.polygons(self.mfgrid.get_cell_vertices(node))
+                for node in range(self.mfgrid.nnodes)
+            ]
+        return np.array(geoms), np.arange(self.mfgrid.nnodes)
 
     def _vtx_grid_to_geoms_cellids(self):
         """internal method, return shapely polygons and cellids for vertex
diff --git a/flopy/utils/hfb_util.py b/flopy/utils/hfb_util.py
new file mode 100644
index 000000000..97d6c8398
--- /dev/null
+++ b/flopy/utils/hfb_util.py
@@ -0,0 +1,366 @@
+import numpy as np
+import pandas as pd
+
+from .geometry import distance
+from .geospatial_utils import GeoSpatialUtil
+from .gridintersect import GridIntersect
+
+
+def _min_distance_index(x0, x1, y0, y1):
+    """
+    Method to get the index of the minimum distance between points
+
+    Parameters
+    ----------
+    x0 : np.ndarray
+        array of x coordinates
+    x1 : float
+        x coordinate of a point
+    y0 : np.ndarray
+        array of y coordinates
+    y1 : float
+        y coordinate of a point
+
+    Returns
+    -------
+        int : index location
+    """
+    dist = distance(x0, y0, x1, y1)
+    idx = np.where(dist == np.min(dist))[0]
+    return idx
+
+
+def _edge_length_lut(xyverts):
+    """
+    Method to create a look-up table based on index of
+    grid cell edge lengths
+
+    Parameters
+    ----------
+    xyverts: np.ndarray
+        numpy array of xy coordinate pairs
+
+    Returns
+    -------
+        dict
+    """
+    xy0 = xyverts.T[:, :-1]
+    xy1 = xyverts.T[:, 1:]
+    dist = distance(xy0[0], xy0[1], xy1[0], xy1[1])
+    dist_lu = {}
+    for v, d in enumerate(dist):
+        if v < len(dist) - 1:
+            dist_lu[(v, v + 1)] = d
+        else:
+            dist_lu[(0, v)] = d
+    return dist_lu
+
+
+def _minimize_hfb_deviance(idxs, xyverts, pts):
+    """
+    Method to minimumize the distance of the HFB fault trace from
+    the intersection points of the fault and model grid cell. Used
+    primarily for tie-breakers where there is not a clear "routing"
+    option based on vertex to vertex edge length distance
+
+    Parameters
+    ----------
+    idxs : iterable
+        list, tuple, or numpy array of polygon vertex indices for test hfb path
+    xyverts : np.array
+        numpy array of (x,y) polygon vertices
+    pts : np.array
+        numpy array of (x,y) points of intersection between line and polygon
+
+    Returns
+    -------
+        float : sum of minimum distances between edge centroids across a potential
+        fault routing option and the intersection points between the fault line
+        and grid node
+    """
+    lcs = []
+    xyverts = xyverts.T
+    pts = pts.T
+    for ix in range(1, len(idxs)):
+        ixx = (idxs[ix - 1], idxs[ix])
+        xc = np.mean(xyverts[0, ixx])
+        yc = np.mean(xyverts[1, ixx])
+        lcs.append([xc, yc])
+
+    mins = []
+    for xc, yc in lcs:
+        dists = distance(pts[0], pts[1], xc, yc)
+        mins.append(np.min(dists))
+
+    return np.sum(mins)
+
+
+def _perturb_intersection_coords(idxs, xyverts, ipt, epsilon=1e-06):
+    """
+    Method to perturb an intersection location by a small amount which is used
+    to handle colinear intersection and intersection at the midpoint of a cell edge.
+
+    Parameters
+    ----------
+    idxs : list
+        list of vertex indices
+    xyverts : np.array
+        numpy array of x,y vertices for a cell
+    ipt : iterable
+        intersection point x,y coordinate
+    epsilon : float
+        perturbation value
+
+    Returns
+    -------
+        ipt : list holding x, y coordinate pair of the perturbed intersection
+    """
+    xyverts = xyverts.T
+    xverts = xyverts[0][idxs]
+    yverts = xyverts[1][idxs]
+
+    if (xverts[0] - xverts[1]) == 0:
+        # vertical line
+        new_vrt = [xverts[0], np.mean(yverts) - epsilon]
+    elif (yverts[0] - yverts[1]) == 0:
+        # horizontal line
+        new_vrt = [np.mean(xverts) - epsilon, yverts[0]]
+    else:
+        # need to adjust across the line
+        m = (yverts[1] - yverts[0]) / [xverts[1] - xverts[0]]
+        if xverts[0] != 0:
+            vidx = 0
+        else:
+            vidx = 1
+        b = yverts[vidx] / (m * xverts[vidx])
+        cx = ipt[0] - epsilon
+        cy = m * cx + b
+        new_vrt = [cx, cy]
+
+    ipt = new_vrt
+    return ipt
+
+
+def _edge_neighbors(modelgrid):
+    """
+    Method to get a dictionary of unique node edges (by ivert) and the nodes the
+    edge is shared between
+
+    Parameters
+    ----------
+    modelgrid : flopy.discretization.Grid object
+
+    Returns
+    -------
+        dict : dictionary of {edge iverts : [nodes]}
+    """
+    node_num = 0
+    geoms = []
+    node_nums = []
+
+    for poly in modelgrid.iverts:
+        poly = [int(i) for i in poly]
+        if poly[0] == poly[-1]:
+            poly = poly[:-1]
+        for v in range(len(poly)):
+            geoms.append(tuple(sorted([poly[v - 1], poly[v]])))
+        node_nums += [node_num] * len(poly)
+        node_num += 1
+
+    edge_nodes = {}
+    for i, item in enumerate(geoms):
+        if item not in edge_nodes:
+            edge_nodes[item] = {
+                node_nums[i],
+            }
+        else:
+            edge_nodes[item].add(node_nums[i])
+
+    return edge_nodes
+
+
+def make_hfb_array(modelgrid, geom):
+    """
+    Method to make a HFB recarray from geospatial linestring information.
+    Note: this method was developed to only accept a single linestring at a
+    time to give the user an opportunity to save unique fault information
+    for calibration purposes.
+
+    Parameters
+    ----------
+    modelgrid : flopy.discretization.Grid object
+        FloPy StructuredGrid, VertexGrid, and UnstructuredGrid are supported
+    geom : geospatial object
+        geom parameter is a geospatial LineString representation.
+        Shapely.geometry.LineString, GeoJson, List of vertices,
+        shapefile.shape types are supported.
+
+    Returns
+    -------
+        np.recarray
+        numpy recarray of cellid1, cellid2, and hydchr. The hydchr field must
+        be set by the user, is set to NaN, and is provided for convenience.
+    """
+    gu = GeoSpatialUtil(geom, "LineString")
+    if gu.shapetype.lower() != "linestring":
+        raise AssertionError(
+            f"{gu.shapetype} is not supported, only LineStrings are supported"
+        )
+
+    geom = gu.points
+
+    if modelgrid.idomain is not None:
+        idomain = modelgrid.idomain.ravel()
+    else:
+        idomain = np.ones((modelgrid.nnodes,), dtype=int)
+    edge_set = _edge_neighbors(modelgrid)
+    vert_ivert = {tuple(i): cnt for cnt, i in enumerate(modelgrid.verts)}
+    iverts = modelgrid.iverts
+    xverts, yverts = modelgrid.cross_section_vertices
+    verts = []
+    for node, xv in enumerate(xverts):
+        yv = yverts[node]
+        vrts = list(zip(xv, yv))
+        if len(vrts) == len(iverts[node]):
+            if iverts[node][0] != iverts[node][-1]:
+                vrts.append(vrts[0])
+        verts.append(np.array(vrts))
+
+    ixs = GridIntersect(modelgrid)
+    result = ixs.intersect(geom, shapetype="LineString")
+    result_adj = []
+    for record in result:
+        node = record.cellid
+        ixshp = record.ixshapes
+        # todo: numpy array this and Transpose
+        coords = np.array(ixshp.coords.xy).T
+        x0, y0 = coords[0, 0], coords[0, 1]
+        x1, y1 = coords[-1, 0], coords[-1, 1]
+
+        xycell = verts[node]
+        xcell = xycell.T[0][:-1]
+        ycell = xycell.T[1][:-1]
+
+        vidx0 = _min_distance_index(xcell, x0, ycell, y0)
+        if len(vidx0) > 1:
+            # perturb line by small epsilon
+            coords[0] = _perturb_intersection_coords(vidx0, xycell, coords[0])
+            x0, y0 = coords[0, 0], coords[0, 1]
+            vidx0 = _min_distance_index(xcell, x0, ycell, y0)
+
+        vidx0 = vidx0[0]
+
+        vidx1 = _min_distance_index(xcell, x1, ycell, y1)
+        if len(vidx1) > 1:
+            # perturb line by small epsilon
+            coords[-1] = _perturb_intersection_coords(vidx1, xycell, coords[-1])
+            x1, y1 = coords[-1, 0], coords[-1, 1]
+            vidx1 = _min_distance_index(xcell, x1, ycell, y1)
+
+        vidx1 = vidx1[0]
+
+        if vidx0 == vidx1:
+            continue
+
+        tmp = tuple(sorted([vidx0, vidx1]))
+
+        if tmp[1] - tmp[0] > 1:
+            nvert = len(xycell) - 1
+            elens = _edge_length_lut(xycell)
+            # construct line routing options
+            o1 = list(range(tmp[0], tmp[1] + 1))
+            o2 = list(range(tmp[1], nvert)) + list(range(0, tmp[0] + 1))
+            # calculate routing distance
+            d1 = np.sum(
+                [elens[tuple(sorted([o1[ix - 1], o1[ix]]))] for ix in range(1, len(o1))]
+            )
+            d2 = np.sum(
+                [elens[tuple(sorted([o2[ix - 1], o2[ix]]))] for ix in range(1, len(o2))]
+            )
+
+            # evaluate distance and break ties if necessary
+            if d1 < d2:
+                tmp = o1
+            elif d2 < d1:
+                tmp = o2
+            else:
+                om1 = _minimize_hfb_deviance(o1, xycell, coords)
+                om2 = _minimize_hfb_deviance(o2, xycell, coords)
+                if om1 <= om2:
+                    tmp = o1
+                else:
+                    tmp = o2
+
+        edges = []
+        for ix in range(1, len(tmp)):
+            eix0 = tmp[ix - 1]
+            eix1 = tmp[ix]
+
+            xyv0 = tuple(xycell[eix0])
+            xyv1 = tuple(xycell[eix1])
+
+            iv0 = vert_ivert[xyv0]
+            iv1 = vert_ivert[xyv1]
+
+            edges.append(tuple(sorted([iv0, iv1])))
+
+        hfb_neighs = []
+        for edge in edges:
+            for n in edge_set[edge]:
+                if n == node:
+                    continue
+                hfb_neighs.append(n)
+
+        res = [int(record.cellid), hfb_neighs]
+        result_adj.append(res)
+
+    hfb_data = []
+    visited = []
+    if modelgrid.nlay is not None and modelgrid.grid_type != "unstructured":
+        for lay in range(modelgrid.nlay):
+            ncpl_adj = lay * modelgrid.ncpl
+            for cid0, hfb_neighs in result_adj:
+                if not idomain[cid0 + ncpl_adj]:
+                    continue
+
+                if modelgrid.grid_type == "structured":
+                    cellid0 = modelgrid.get_lrc(cid0 + ncpl_adj)[0]
+                else:
+                    cellid0 = (lay, cid0)
+
+                for cid1 in hfb_neighs:
+                    if not idomain[cid1 + ncpl_adj]:
+                        continue
+
+                    if modelgrid.grid_type == "structured":
+                        cellid1 = modelgrid.get_lrc(cid1 + ncpl_adj)[0]
+                    else:
+                        cellid1 = (lay, cid1)
+
+                    if cellid0 + cellid1 in visited:
+                        continue
+                    elif cellid1 + cellid0 in visited:
+                        continue
+
+                    visited.append(cellid0 + cellid1)
+                    hfb_data.append((cellid0, cellid1))
+    else:
+        for cellid0, hfb_neighs in result_adj:
+            if not idomain[cellid0]:
+                continue
+
+            for cellid1 in hfb_neighs:
+                if not idomain[cellid1]:
+                    continue
+
+                if (cellid0, cellid1) in visited:
+                    continue
+                elif (cellid1, cellid0) in visited:
+                    continue
+
+                visited.append((cellid0, cellid1))
+                hfb_data.append(((cellid0,), (cellid1,)))
+
+    df = pd.DataFrame(hfb_data, columns=["cellid1", "cellid2"])
+    df["hydchr"] = np.nan
+    return df.to_records(index=False)
diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py
index c9085e2c4..e2ef04f49 100644
--- a/flopy/utils/voronoi.py
+++ b/flopy/utils/voronoi.py
@@ -4,7 +4,7 @@
 import numpy as np
 
 from .cvfdutil import get_disv_gridprops
-from .geometry import point_in_polygon
+from .geometry import distance, point_in_polygon
 from .utl_import import import_optional_dependency
 
 
@@ -309,8 +309,92 @@ def get_disu5_gridprops(self):
         raise NotImplementedError(msg)
 
     def get_disu6_gridprops(self):
-        msg = "This method is not implemented yet."
-        raise NotImplementedError(msg)
+        """
+        Get a dictionary of arguments that can be passed in to the
+        flopy.mf6.ModflowGwfdisu class
+
+        Returns
+        -------
+        disu_gridprops : dict
+            Dictionary of arguments that can be unpacked into the
+            flopy.mf6.ModflowGwfdisu constructor
+
+        """
+        from ..discretization import UnstructuredGrid
+
+        gridprops = self.get_gridprops_unstructuredgrid()
+        ugrid = UnstructuredGrid(**gridprops)
+
+        xcenters = ugrid.xcellcenters
+        ycenters = ugrid.ycellcenters
+        xvertices = ugrid.xvertices
+        yvertices = ugrid.yvertices
+        iverts = ugrid.iverts
+        neighbors = ugrid.neighbors()
+        iac = []
+        ja = []
+        ihc = []
+        cl12 = []
+        hwva = []
+        angldegx = []
+        for node, nnodes in neighbors.items():
+            iac.append(len(nnodes) + 1)
+            ja.extend([node] + nnodes)
+            ihc.extend(
+                [
+                    1,
+                ]
+                * iac[-1]
+            )
+            # cell center to cell center distance...
+            cl12.append(0)
+            hwva.append(0)
+            angldegx.append(0)
+
+            xc = xcenters[node]
+            yc = ycenters[node]
+            nxc = xcenters[nnodes]
+            nyc = ycenters[nnodes]
+            dists = distance(nxc, nyc, xc, yc)
+            cl12.extend(dists)
+
+            angles = np.arctan2(nxc - xc, nyc - yc) * 180.0 / np.pi
+            angles = np.where(angles < 0, angles + 360, angles)
+            angldegx.extend(angles)
+
+            ivrts = iverts[node]
+            xverts = xvertices[node]
+            yverts = yvertices[node]
+            for n in nnodes:
+                nivrts = iverts[n]
+                common = set(ivrts) & set(nivrts)
+                idxs = [ivrts.index(i) for i in common]
+                hwv = distance(
+                    xverts[idxs[0]], yverts[idxs[0]], xverts[idxs[1]], yverts[idxs[1]]
+                )
+                hwva.append(hwv)
+
+        cell2d = []
+        for ix, ivrts in enumerate(iverts):
+            rec = [ix, xcenters[ix], ycenters[ix], len(ivrts)] + ivrts
+            cell2d.append(rec)
+
+        gridprops["iac"] = iac
+        gridprops["ja"] = ja
+        gridprops["ihc"] = ihc
+        gridprops["cl12"] = cl12
+        gridprops["hwva"] = hwva
+        gridprops["angldegx"] = angldegx
+        gridprops["area"] = ugrid.cell_area
+        gridprops["nodes"] = len(iac)
+        gridprops["nja"] = len(ja)
+        gridprops["nvert"] = len(gridprops["vertices"])
+        gridprops["cell2d"] = cell2d
+        gridprops.pop("ncpl")
+        gridprops.pop("iverts")
+        gridprops.pop("xcenters")
+        gridprops.pop("ycenters")
+        return gridprops
 
     def get_gridprops_vertexgrid(self):
         """