diff --git a/autotest/test_export.py b/autotest/test_export.py
index 4736260ea..f6c050633 100644
--- a/autotest/test_export.py
+++ b/autotest/test_export.py
@@ -1788,51 +1788,51 @@ def test_vtk_export_disv1_model(function_tmpdir):
         idomain=np.ones((nlay, nrow, ncol)),
     )
 
-    with pytest.deprecated_call():
-        from flopy.utils.cvfdutil import gridlist_to_disv_gridprops
-
-        gridprops = gridlist_to_disv_gridprops([mg])
-        gridprops["top"] = 0
-        gridprops["botm"] = np.zeros((nlay, nrow * ncol), dtype=float) - 1
-        gridprops["nlay"] = nlay
-
-        disv = ModflowGwfdisv(gwf, **gridprops)
-        ic = ModflowGwfic(gwf, strt=10)
-        npf = ModflowGwfnpf(gwf)
-
-        # Export model without specifying packages_names parameter
-        # create the vtk output
-        gwf = sim.get_model()
-        vtkobj = Vtk(gwf, binary=False)
-        vtkobj.add_model(gwf)
-        f = function_tmpdir / "gwf.vtk"
-        vtkobj.write(f)
-
-        # load the output using the vtk standard library
-        gridreader = vtkUnstructuredGridReader()
-        gridreader.SetFileName(str(f))
-        gridreader.Update()
-        grid = gridreader.GetOutput()
-
-        # get the points
-        vtk_points = grid.GetPoints()
-        vtk_points = vtk_points.GetData()
-        vtk_points = vtk_to_numpy(vtk_points)
-
-        # get cell locations (ia format of point to cell relationship)
-        cell_locations = vtk_to_numpy(grid.GetCellLocationsArray())
-        cell_locations_answer = np.array([0, 8, 16, 24, 32, 40, 48, 56, 64])
-        print(f"Found cell locations {cell_locations} in vtk file.")
-        print(f"Expecting cell locations {cell_locations_answer}")
-        errmsg = "vtk cell locations do not match expected result."
-        assert np.allclose(cell_locations, cell_locations_answer), errmsg
-
-        cell_types = vtk_to_numpy(grid.GetCellTypesArray())
-        cell_types_answer = np.array(9 * [42])
-        print(f"Found cell types {cell_types} in vtk file.")
-        print(f"Expecting cell types {cell_types_answer}")
-        errmsg = "vtk cell types do not match expected result."
-        assert np.allclose(cell_types, cell_types_answer), errmsg
+    from flopy.utils.cvfdutil import get_disv_gridprops, gridlist_to_verts
+
+    verts, iverts = gridlist_to_verts([mg])
+    gridprops = get_disv_gridprops(verts, iverts)
+    gridprops["top"] = 0
+    gridprops["botm"] = np.zeros((nlay, nrow * ncol), dtype=float) - 1
+    gridprops["nlay"] = nlay
+
+    disv = ModflowGwfdisv(gwf, **gridprops)
+    ic = ModflowGwfic(gwf, strt=10)
+    npf = ModflowGwfnpf(gwf)
+
+    # Export model without specifying packages_names parameter
+    # create the vtk output
+    gwf = sim.get_model()
+    vtkobj = Vtk(gwf, binary=False)
+    vtkobj.add_model(gwf)
+    f = function_tmpdir / "gwf.vtk"
+    vtkobj.write(f)
+
+    # load the output using the vtk standard library
+    gridreader = vtkUnstructuredGridReader()
+    gridreader.SetFileName(str(f))
+    gridreader.Update()
+    grid = gridreader.GetOutput()
+
+    # get the points
+    vtk_points = grid.GetPoints()
+    vtk_points = vtk_points.GetData()
+    vtk_points = vtk_to_numpy(vtk_points)
+
+    # get cell locations (ia format of point to cell relationship)
+    cell_locations = vtk_to_numpy(grid.GetCellLocationsArray())
+    cell_locations_answer = np.array([0, 8, 16, 24, 32, 40, 48, 56, 64])
+    print(f"Found cell locations {cell_locations} in vtk file.")
+    print(f"Expecting cell locations {cell_locations_answer}")
+    errmsg = "vtk cell locations do not match expected result."
+    assert np.allclose(cell_locations, cell_locations_answer), errmsg
+
+    cell_types = vtk_to_numpy(grid.GetCellTypesArray())
+    cell_types_answer = np.array(9 * [42])
+    print(f"Found cell types {cell_types} in vtk file.")
+    print(f"Expecting cell types {cell_types_answer}")
+    errmsg = "vtk cell types do not match expected result."
+    assert np.allclose(cell_types, cell_types_answer), errmsg
 
 
 @pytest.mark.mf6
@@ -2228,7 +2228,7 @@ def test_mf6_chd_shapefile_export_unstructured(function_tmpdir, use_pandas, spar
 
 def disv_sim(name, tmpdir):
     from flopy.discretization import StructuredGrid
-    from flopy.utils.cvfdutil import gridlist_to_disv_gridprops
+    from flopy.utils.cvfdutil import get_disv_gridprops, gridlist_to_verts
 
     nlay, nrow, ncol = 3, 3, 3
     mg = StructuredGrid(
@@ -2239,12 +2239,12 @@ def disv_sim(name, tmpdir):
         idomain=np.ones((nlay, nrow, ncol)),
     )
 
-    with pytest.deprecated_call():
-        gridprops = gridlist_to_disv_gridprops([mg])
-        gridprops["top"] = 0
-        ncpl = gridprops["ncpl"]
-        gridprops["botm"] = np.zeros((nlay, ncpl), dtype=float) - 1
-        gridprops["nlay"] = nlay
+    verts, iverts = gridlist_to_verts([mg])
+    gridprops = get_disv_gridprops(verts, iverts)
+    gridprops["top"] = 0
+    ncpl = gridprops["ncpl"]
+    gridprops["botm"] = np.zeros((nlay, ncpl), dtype=float) - 1
+    gridprops["nlay"] = nlay
 
     sim = MFSimulation(sim_name=name, sim_ws=tmpdir, exe_name="mf6")
     tdis = ModflowTdis(sim)
diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index 6f1211e60..c390f9b15 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -23,9 +23,9 @@
 from flopy.utils.cvfdutil import (
     area_of_polygon,
     centroid_of_polygon,
-    gridlist_to_disv_gridprops,
     to_cvfd,
 )
+from flopy.utils.lgrutil import Lgr
 from flopy.utils.triangle import Triangle
 from flopy.utils.voronoi import VoronoiGrid
 
@@ -1015,9 +1015,10 @@ def test_tocvfd3():
     delc = 100.0 * np.ones(nrow)
     tp = np.zeros((nrow, ncol))
     bt = -100.0 * np.ones((nlay, nrow, ncol))
-    idomain = np.ones((nlay, nrow, ncol))
-    idomain[:, 2:5, 2:5] = 0
-    sg1 = StructuredGrid(delr=delr, delc=delc, top=tp, botm=bt, idomain=idomain)
+    # Create refine_mask: cells with value 0 will be refined
+    refine_mask = np.ones((nlay, nrow, ncol))
+    refine_mask[:, 2:5, 2:5] = 0
+    sg1 = StructuredGrid(delr=delr, delc=delc, top=tp, botm=bt)
     # inner grid
     nlay = 1
     nrow = ncol = 9
@@ -1036,26 +1037,27 @@ def test_tocvfd3():
         idomain=idomain,
     )
 
-    with pytest.deprecated_call():
-        gridprops = gridlist_to_disv_gridprops([sg1, sg2])
-        assert "ncpl" in gridprops
-        assert "nvert" in gridprops
-        assert "vertices" in gridprops
-        assert "cell2d" in gridprops
-
-        ncpl = gridprops["ncpl"]
-        nvert = gridprops["nvert"]
-        vertices = gridprops["vertices"]
-        cell2d = gridprops["cell2d"]
-        assert ncpl == 121
-        assert nvert == 148
-        assert len(vertices) == nvert
-        assert len(cell2d) == 121
-
-        # spot check information for cell 28 (zero based)
-        answer = [28, 250.0, 150.0, 7, 38, 142, 143, 45, 46, 44, 38]
-        for i, j in zip(cell2d[28], answer):
-            assert i == j, f"{i} not equal {j}"
+    # Use Lgr class to create DISV grid from nested grids
+    lgr = Lgr.from_parent_grid(sg1, refine_mask, ncpp=3, ncppl=1)
+    gridprops = lgr.to_disv_gridprops()
+    assert "ncpl" in gridprops
+    assert "nvert" in gridprops
+    assert "vertices" in gridprops
+    assert "cell2d" in gridprops
+
+    ncpl = gridprops["ncpl"]
+    nvert = gridprops["nvert"]
+    vertices = gridprops["vertices"]
+    cell2d = gridprops["cell2d"]
+    assert ncpl == 121
+    assert nvert == 148  # Lgr includes hanging vertices and deduplicates vertices
+    assert len(vertices) == nvert
+    assert len(cell2d) == 121
+
+    # spot check information for cell 28 (zero based)
+    answer = [28, 250.0, 150.0, 6, 42, 142, 143, 43, 51, 50]
+    for i, j in zip(cell2d[28], answer):
+        assert i == j, f"{i} not equal {j}"
 
 
 @requires_pkg("shapely")
diff --git a/autotest/test_lgrutil.py b/autotest/test_lgrutil.py
index f84f071b6..9ff50c419 100644
--- a/autotest/test_lgrutil.py
+++ b/autotest/test_lgrutil.py
@@ -220,6 +220,63 @@ def test_lgr_hanging_vertices():
     b[2] = -3 * dz
     assert np.allclose(gridprops["botm"], b)
 
+    # Verify equivalence with legacy gridlist_to_verts approach
+    # This ensures the Lgr approach produces the same results.
+    # Note: Use single layer for comparison since
+    # gridlist_to_verts counts each layer separately
+    from flopy.utils.cvfdutil import get_disv_gridprops, gridlist_to_verts
+
+    # Create single-layer grids for legacy comparison
+    nlay_compare = 1
+    topp_compare = topp[nrowp // 2, ncolp // 2] * np.ones((nrowp, ncolp))
+    botmp_compare = botmp[0:1, :, :]  # Just first layer
+    idomainp_compare = idomainp[0:1, :, :]
+
+    # Child grid (single layer)
+    sg_child = StructuredGrid(
+        delr=np.full(ncpp, dx / ncpp),
+        delc=np.full(ncpp, dy / ncpp),
+        top=np.full((ncpp, ncpp), topp_compare[nrowp // 2, ncolp // 2]),
+        botm=botmp_compare[:, nrowp // 2 : nrowp // 2 + 1, ncolp // 2 : ncolp // 2 + 1]
+        .repeat(ncpp, axis=1)
+        .repeat(ncpp, axis=2),
+        idomain=np.ones((nlay_compare, ncpp, ncpp), dtype=int),
+        xoff=xoffp + dx * (ncolp // 2),
+        yoff=yoffp + dy * (nrowp // 2),
+    )
+
+    # Parent grid (single layer)
+    sg_parent_legacy = StructuredGrid(
+        delr=delrp,
+        delc=delcp,
+        top=topp_compare,
+        botm=botmp_compare,
+        idomain=idomainp_compare,
+        xoff=xoffp,
+        yoff=yoffp,
+    )
+
+    # Generate gridprops using legacy approach
+    verts_legacy, iverts_legacy = gridlist_to_verts([sg_parent_legacy, sg_child])
+    gridprops_legacy = get_disv_gridprops(verts_legacy, iverts_legacy)
+
+    # Verify ncpl and nvert match
+    assert gridprops["ncpl"] == gridprops_legacy["ncpl"]
+    assert gridprops["nvert"] == gridprops_legacy["nvert"]
+
+    # Verify vertex positions match (may be in different order)
+    verts_lgr_set = {(round(v[1], 9), round(v[2], 9)) for v in gridprops["vertices"]}
+    verts_legacy_set = {
+        (round(v[1], 9), round(v[2], 9)) for v in gridprops_legacy["vertices"]
+    }
+    assert verts_lgr_set == verts_legacy_set
+
+    # Verify cell centroids match
+    for icpl in range(gridprops["ncpl"]):
+        xc_lgr, yc_lgr = gridprops["cell2d"][icpl][1:3]
+        xc_legacy, yc_legacy = gridprops_legacy["cell2d"][icpl][1:3]
+        assert np.allclose([xc_lgr, yc_lgr], [xc_legacy, yc_legacy])
+
 
 def test_lgr_from_parent_grid():
     # Create a parent grid with center cells marked for refinement
@@ -461,347 +518,3 @@ def test_simple_regular_grid_deprecation():
     # Verify backward compatibility attributes
     assert grid.xorigin == xorigin
     assert grid.yorigin == yorigin
-
-
-def test_lgr_matches_legacy_gridlist_to_disv_gridprops():
-    """
-    Lgr.to_disv_gridprops() output should match
-    legacy gridlist_to_disv_gridprops() output.
-    """
-    # TODO: this will need to be removed soon, and the underlying
-    # implementation inlined, when we remove this function.
-    from flopy.utils.cvfdutil import gridlist_to_disv_gridprops
-
-    # 7x7 parent grid with 3x3 refinement in center
-    nlay = 1
-    nrow_p, ncol_p = 7, 7
-    delr_p = delc_p = 100.0 * np.ones(7)
-    top_p = np.zeros((nrow_p, ncol_p))
-    botm_p = -100.0 * np.ones((nlay, nrow_p, ncol_p))
-
-    # legacy approach
-    idomain_p = np.ones((nlay, nrow_p, ncol_p), dtype=int)
-    idomain_p[:, 2:5, 2:5] = 0  # Center 3x3 inactive
-    sg_parent_legacy = StructuredGrid(
-        delc=delc_p,
-        delr=delr_p,
-        top=top_p,
-        botm=botm_p,
-        idomain=idomain_p,
-        xoff=0.0,
-        yoff=0.0,
-    )
-    nrow_c, ncol_c = 9, 9
-    delr_c = delc_c = 100.0 / 3.0 * np.ones(9)
-    top_c = np.zeros((nrow_c, ncol_c))
-    botm_c = -100.0 * np.ones((nlay, nrow_c, ncol_c))
-    idomain_c = np.ones((nlay, nrow_c, ncol_c), dtype=int)
-    sg_child = StructuredGrid(
-        delc=delc_c,
-        delr=delr_c,
-        top=top_c,
-        botm=botm_c,
-        idomain=idomain_c,
-        xoff=200.0,
-        yoff=200.0,
-    )
-    gridprops_legacy = gridlist_to_disv_gridprops([sg_parent_legacy, sg_child])
-
-    # lgr approach
-    parent_grid = StructuredGrid(
-        delc=delc_p,
-        delr=delr_p,
-        top=top_p,
-        botm=botm_p,
-        xoff=0.0,
-        yoff=0.0,
-    )
-    refine_mask = np.ones((nlay, nrow_p, ncol_p))
-    refine_mask[:, 2:5, 2:5] = 0
-    lgr = Lgr.from_parent_grid(parent_grid, refine_mask, ncpp=3)
-    gridprops_lgr = lgr.to_disv_gridprops()
-
-    # Note: legacy gridlist_to_disv_gridprops doesn't return nlay, top, botm
-    # Only check common keys
-    assert gridprops_legacy["ncpl"] == gridprops_lgr["ncpl"]
-    assert gridprops_legacy["nvert"] == gridprops_lgr["nvert"]
-
-    # Verify vertex positions match (may be in different order)
-    verts_legacy_set = {
-        (round(v[1], 9), round(v[2], 9)) for v in gridprops_legacy["vertices"]
-    }
-    verts_lgr_set = {
-        (round(v[1], 9), round(v[2], 9)) for v in gridprops_lgr["vertices"]
-    }
-    assert verts_legacy_set == verts_lgr_set
-
-    for icpl in range(gridprops_legacy["ncpl"]):
-        nverts_legacy = gridprops_legacy["cell2d"][icpl][3]
-        nverts_lgr = gridprops_lgr["cell2d"][icpl][3]
-        assert nverts_legacy >= 3
-        assert nverts_lgr >= 3
-
-    for icpl in range(gridprops_legacy["ncpl"]):
-        xc_legacy, yc_legacy = gridprops_legacy["cell2d"][icpl][1:3]
-        xc_lgr, yc_lgr = gridprops_lgr["cell2d"][icpl][1:3]
-        assert np.allclose([xc_legacy, yc_legacy], [xc_lgr, yc_lgr])
-
-    # Note: top and botm are only in Lgr gridprops,
-    # not in legacy gridlist_to_disv_gridprops.
-    # Just verify they exist in the Lgr version
-    assert "nlay" in gridprops_lgr
-    assert "top" in gridprops_lgr
-    assert "botm" in gridprops_lgr
-    assert gridprops_lgr["nlay"] == nlay
-    assert gridprops_lgr["top"].shape == (gridprops_lgr["ncpl"],)
-    assert gridprops_lgr["botm"].shape == (nlay, gridprops_lgr["ncpl"])
-
-
-def test_hanging_vertex_algorithm_comparison():
-    """
-    Compare the two hanging vertex placement algorithms in detail:
-    1. to_cvfd() - reactive geometric algorithm (legacy)
-    2. Lgr.find_hanging_vertices() - proactive structured algorithm (new)
-
-    This test creates a simple parent-child grid and analyzes the differences
-    in vertex placement, connectivity, and grid structure. Intended mainly to
-    be run manually to inspect the printed output.
-    """
-    # Create a simple 3x3 parent grid with 3x refinement in center cell
-    nlay = 1
-    nrow_p = 3
-    ncol_p = 3
-    delr_p = 100.0
-    delc_p = 100.0
-    top_p = np.zeros((nrow_p, ncol_p))
-    botm_p = np.array([np.full((nrow_p, ncol_p), -10.0)])
-
-    # Setup for both approaches
-    # Legacy approach: use idomain to mark refinement region as inactive
-    idomain_legacy = np.ones((nlay, nrow_p, ncol_p), dtype=int)
-    idomain_legacy[:, 1, 1] = 0  # Center cell inactive
-
-    sg_parent_legacy = StructuredGrid(
-        delr=np.full(ncol_p, delr_p),
-        delc=np.full(nrow_p, delc_p),
-        top=top_p,
-        botm=botm_p,
-        idomain=idomain_legacy,
-        xoff=0.0,
-        yoff=0.0,
-    )
-
-    # Child grid at center position
-    ncpp = 3
-    top_c = np.zeros((ncpp, ncpp))
-    botm_c = np.array([np.full((ncpp, ncpp), -10.0)])
-    idomain_child = np.ones((nlay, ncpp, ncpp), dtype=int)
-    sg_child = StructuredGrid(
-        delr=np.full(ncpp, delr_p / ncpp),
-        delc=np.full(ncpp, delc_p / ncpp),
-        top=top_c,
-        botm=botm_c,
-        idomain=idomain_child,
-        xoff=100.0,
-        yoff=100.0,
-    )
-
-    # New approach: use refine_mask
-    sg_parent_new = StructuredGrid(
-        delr=np.full(ncol_p, delr_p),
-        delc=np.full(nrow_p, delc_p),
-        top=top_p.copy(),
-        botm=botm_p.copy(),
-        xoff=0.0,
-        yoff=0.0,
-    )
-    refine_mask = np.ones((nlay, nrow_p, ncol_p))
-    refine_mask[:, 1, 1] = 0
-
-    # Generate grids with both algorithms
-    # Legacy: to_cvfd() algorithm
-    verts_legacy, iverts_legacy = gridlist_to_verts([sg_parent_legacy, sg_child])
-    gridprops_legacy = get_disv_gridprops(verts_legacy, iverts_legacy)
-
-    # New: Lgr.find_hanging_vertices() algorithm
-    lgr = Lgr.from_parent_grid(sg_parent_new, refine_mask, ncpp=ncpp)
-    gridprops_new = lgr.to_disv_gridprops()
-
-    # Compare basic metrics
-    print("\n" + "=" * 80)
-    print("HANGING VERTEX ALGORITHM COMPARISON")
-    print("=" * 80)
-    print("\nGrid Configuration:")
-    print(f"  Parent: {nrow_p}x{ncol_p} cells, {delr_p}x{delc_p} spacing")
-    print(f"  Child: {ncpp}x{ncpp} refinement at center cell (1,1)")
-    print(
-        f"  Total cells: {nrow_p * ncol_p - 1} "
-        f"parent + {ncpp * ncpp} "
-        f"child = {nrow_p * ncol_p - 1 + ncpp * ncpp}"
-    )
-
-    print(
-        f"\n{'Metric':<30} "
-        f"{'Legacy (to_cvfd)':<20} "
-        f"{'New (Lgr)':<20} "
-        f"{'Difference':<20}"
-    )
-    print("-" * 90)
-    print(
-        f"{'Total vertices':<30} "
-        f"{gridprops_legacy['nvert']:<20} "
-        f"{gridprops_new['nvert']:<20} "
-        f"{gridprops_new['nvert'] - gridprops_legacy['nvert']:<20}"
-    )
-    print(
-        f"{'Total cells':<30} "
-        f"{gridprops_legacy['ncpl']:<20} "
-        f"{gridprops_new['ncpl']:<20} "
-        f"{gridprops_new['ncpl'] - gridprops_legacy['ncpl']:<20}"
-    )
-
-    # Analyze vertex differences
-    verts_legacy_dict = {
-        i: (v[1], v[2]) for i, v in enumerate(gridprops_legacy["vertices"])
-    }
-    verts_new_dict = {i: (v[1], v[2]) for i, v in enumerate(gridprops_new["vertices"])}
-
-    # Find unique vertices in each approach
-    verts_legacy_set = {
-        (round(x, 6), round(y, 6)) for x, y in verts_legacy_dict.values()
-    }
-    verts_new_set = {(round(x, 6), round(y, 6)) for x, y in verts_new_dict.values()}
-
-    only_in_legacy = verts_legacy_set - verts_new_set
-    only_in_new = verts_new_set - verts_legacy_set
-    common = verts_legacy_set & verts_new_set
-
-    print(f"\n{'Vertex Analysis':<30} {'Count':<20}")
-    print("-" * 50)
-    print(f"{'Common vertices':<30} {len(common):<20}")
-    print(f"{'Only in legacy':<30} {len(only_in_legacy):<20}")
-    print(f"{'Only in new':<30} {len(only_in_new):<20}")
-
-    if only_in_legacy:
-        print("\nVertices only in legacy (to_cvfd):")
-        for x, y in sorted(only_in_legacy):
-            print(f"  ({x:8.2f}, {y:8.2f})")
-
-    if only_in_new:
-        print("\nVertices only in new (Lgr):")
-        for x, y in sorted(only_in_new):
-            print(f"  ({x:8.2f}, {y:8.2f})")
-
-    # Analyze cell vertex counts
-    legacy_cell_nverts = [cell[3] for cell in gridprops_legacy["cell2d"]]
-    new_cell_nverts = [cell[3] for cell in gridprops_new["cell2d"]]
-
-    print(f"\n{'Cell Vertex Counts':<30} {'Legacy':<20} {'New':<20}")
-    print("-" * 70)
-    print(
-        f"{'Min vertices per cell':<30} "
-        f"{min(legacy_cell_nverts):<20} "
-        f"{min(new_cell_nverts):<20}"
-    )
-    print(
-        f"{'Max vertices per cell':<30} "
-        f"{max(legacy_cell_nverts):<20} "
-        f"{max(new_cell_nverts):<20}"
-    )
-    print(
-        f"{'Avg vertices per cell':<30} "
-        f"{np.mean(legacy_cell_nverts):<20.2f} "
-        f"{np.mean(new_cell_nverts):<20.2f}"
-    )
-
-    # Find cells with different vertex counts
-    print("\nCells with different vertex counts:")
-    diff_count = 0
-    for icell in range(min(len(legacy_cell_nverts), len(new_cell_nverts))):
-        if legacy_cell_nverts[icell] != new_cell_nverts[icell]:
-            print(
-                f"  Cell {icell}: legacy={legacy_cell_nverts[icell]}, "
-                f"new={new_cell_nverts[icell]}"
-            )
-            diff_count += 1
-    if diff_count == 0:
-        print("  None - all cells have same vertex counts")
-
-    # Identify hanging vertices by analyzing parent-child boundary cells
-    print(f"\n{'Hanging Vertex Analysis':<50}")
-    print("-" * 50)
-
-    # For each parent cell at the boundary, check vertex counts
-    print("\nParent boundary cells (adjacent to refinement):")
-    boundary_parents = [
-        (0, "top-left"),
-        (1, "top-center"),
-        (2, "top-right"),
-        (3, "middle-left"),
-        (5, "middle-right"),
-        (6, "bottom-left"),
-        (7, "bottom-center"),
-        (8, "bottom-right"),
-    ]
-
-    for icell, label in boundary_parents:
-        if icell < len(legacy_cell_nverts):
-            legacy_nv = legacy_cell_nverts[icell]
-            new_nv = new_cell_nverts[icell]
-            status = "DIFFERENT" if legacy_nv != new_nv else "same"
-            print(
-                f"  Cell {icell} ({label:15s}): "
-                f"legacy={legacy_nv}, new={new_nv} [{status}]"
-            )
-
-    # Detailed vertex analysis - check for duplicates
-    print(f"\n{'Duplicate Vertex Analysis':<50}")
-    print("-" * 50)
-
-    # Count occurrences of each coordinate in new approach
-    new_coord_counts = {}
-    for i, (x, y) in verts_new_dict.items():
-        coord = (round(x, 6), round(y, 6))
-        if coord not in new_coord_counts:
-            new_coord_counts[coord] = []
-        new_coord_counts[coord].append(i)
-
-    duplicates_new = {
-        coord: indices
-        for coord, indices in new_coord_counts.items()
-        if len(indices) > 1
-    }
-    if duplicates_new:
-        print(f"Found {len(duplicates_new)} duplicate coordinates in new approach:")
-        for coord, indices in sorted(duplicates_new.items()):
-            print(f"  {coord}: vertex indices {indices}")
-    else:
-        print("No duplicate coordinates in new approach")
-
-    # Count occurrences in legacy approach
-    legacy_coord_counts = {}
-    for i, (x, y) in verts_legacy_dict.items():
-        coord = (round(x, 6), round(y, 6))
-        if coord not in legacy_coord_counts:
-            legacy_coord_counts[coord] = []
-        legacy_coord_counts[coord].append(i)
-
-    duplicates_legacy = {
-        coord: indices
-        for coord, indices in legacy_coord_counts.items()
-        if len(indices) > 1
-    }
-    if duplicates_legacy:
-        print(
-            f"Found {len(duplicates_legacy)} duplicate coordinates in legacy approach:"
-        )
-        for coord, indices in sorted(duplicates_legacy.items()):
-            print(f"  {coord}: vertex indices {indices}")
-    else:
-        print("No duplicate coordinates in legacy approach")
-
-    assert (
-        gridprops_legacy["nvert"] == gridprops_new["nvert"]
-        and len(only_in_legacy) == 0
-        and len(only_in_new) == 0
-    )
diff --git a/autotest/test_postprocessing.py b/autotest/test_postprocessing.py
index 9da2d66f5..0c668077d 100644
--- a/autotest/test_postprocessing.py
+++ b/autotest/test_postprocessing.py
@@ -20,7 +20,6 @@
 from flopy.modflow import Modflow, ModflowDis, ModflowLpf, ModflowUpw
 from flopy.plot import PlotMapView
 from flopy.utils import get_transmissivities
-from flopy.utils.cvfdutil import gridlist_to_disv_gridprops
 from flopy.utils.gridutil import get_disv_kwargs
 from flopy.utils.postprocessing import (
     get_gradients,
diff --git a/flopy/utils/cvfdutil.py b/flopy/utils/cvfdutil.py
index ae7e3cc31..a6b720eee 100644
--- a/flopy/utils/cvfdutil.py
+++ b/flopy/utils/cvfdutil.py
@@ -407,42 +407,3 @@ def get_disv_gridprops(verts, iverts, xcyc=None):
     gridprops["vertices"] = vertices
     gridprops["cell2d"] = cell2d
     return gridprops
-
-
-def gridlist_to_disv_gridprops(gridlist):
-    """
-
-    Take a list of flopy structured model grids and convert them into a
-    dictionary that can be passed into the modflow6 disv package.  Cells from a
-    child grid will patched in to make a single set of vertices.  Cells will
-    be numbered according to consecutive numbering of active cells in the
-    grid list.
-
-    This function is deprecated in 3.8 and will be removed in 3.9.  Use the
-    functionality in flopy.utils.cvfdutil.Lgr() to create a DISV mesh for a
-    nested grid.
-
-    Parameters
-    ----------
-    gridlist : list
-        List of flopy.discretization.modelgrid.  Must be of type structured
-        grids
-
-    Returns
-    -------
-    gridprops : dict
-        Dictionary containing entries that can be passed directly into the
-        modflow6 disv package.
-
-    """
-    warnings.warn(
-        "the gridlist_to_disv_gridprops function is deprecated and will be "
-        "removed in version 3.9. For nested grids, use flopy.utils.lgrutil.Lgr() "
-        "instead. For converting a single structured grid to DISV, use "
-        "get_disv_gridprops(grid.verts, grid.iverts) instead.",
-        PendingDeprecationWarning,
-    )
-
-    verts, iverts = gridlist_to_verts(gridlist)
-    gridprops = get_disv_gridprops(verts, iverts)
-    return gridprops