diff --git a/.docs/Notebooks/gridgen_example.py b/.docs/Notebooks/gridgen_example.py
index c02616408d..bcca0f08e1 100644
--- a/.docs/Notebooks/gridgen_example.py
+++ b/.docs/Notebooks/gridgen_example.py
@@ -25,7 +25,6 @@
 
 # +
 import os
-import shutil
 import sys
 from tempfile import TemporaryDirectory
 
diff --git a/autotest/test_gridgen.py b/autotest/test_gridgen.py
index 1326221173..3eedb0cc8d 100644
--- a/autotest/test_gridgen.py
+++ b/autotest/test_gridgen.py
@@ -11,6 +11,7 @@
 from modflow_devtools.misc import has_pkg
 
 import flopy
+from flopy.discretization.unstructuredgrid import UnstructuredGrid
 from flopy.discretization.vertexgrid import VertexGrid
 from flopy.utils.gridgen import Gridgen
 
@@ -26,7 +27,7 @@ def test_ctor_accepts_path_or_string_model_ws(function_tmpdir):
     assert g.model_ws == function_tmpdir
 
 
-def base_grid():
+def get_structured_grid():
     """Get a small version of the first grid in:
     .docs/Notebooks/gridgen_example.py"""
 
@@ -57,64 +58,96 @@ def base_grid():
 
 
 @requires_exe("gridgen")
-def test_add_active_domain(function_tmpdir):
-    bgrid = base_grid()
-    grids = []
+# GRIDGEN seems not to like paths containing "[" or "]", as
+# function_tmpdir does with parametrization, do it manually
+# @pytest.mark.parametrize("grid_type", ["vertex", "unstructured"])
+def test_add_active_domain(function_tmpdir):  # , grid_type):
+    bgrid = get_structured_grid()
 
     # test providing active domain various ways
-    for feature in [
-        [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]],
-        function_tmpdir / "ad0.shp",
-        function_tmpdir / "ad0",
-        "ad0.shp",
-        "ad0",
-    ]:
-        gridgen = Gridgen(bgrid, model_ws=function_tmpdir)
-        gridgen.add_active_domain(
-            feature,
-            range(bgrid.nlay),
-        )
-        gridgen.build()
-        grid = VertexGrid(**gridgen.get_gridprops_vertexgrid())
-        grid.plot()
-        grids.append(grid)
-        # plt.show()
-
-        assert grid.nnodes < bgrid.nnodes
-        assert grid.ncpl != bgrid.ncpl
-        assert all(grid.ncpl == g.ncpl for g in grids)
-        assert all(grid.nnodes == g.nnodes for g in grids)
+    for grid_type in ["vertex", "unstructured"]:
+        grids = []
+        for feature in [
+            [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]],
+            function_tmpdir / "ad0.shp",
+            function_tmpdir / "ad0",
+            "ad0.shp",
+            "ad0",
+        ]:
+            print(
+                "Testing add_active_domain() for",
+                grid_type,
+                "grid with features",
+                feature,
+            )
+            gridgen = Gridgen(bgrid, model_ws=function_tmpdir)
+            gridgen.add_active_domain(
+                feature,
+                range(bgrid.nlay),
+            )
+            gridgen.build()
+            grid = (
+                VertexGrid(**gridgen.get_gridprops_vertexgrid())
+                if grid_type == "vertex"
+                else UnstructuredGrid(
+                    **gridgen.get_gridprops_unstructuredgrid()
+                )
+            )
+            grid.plot()
+            grids.append(grid)
+            # plt.show()
+
+            assert grid.nnodes < bgrid.nnodes
+            assert not np.array_equal(grid.ncpl, bgrid.ncpl)
+            assert all(np.array_equal(grid.ncpl, g.ncpl) for g in grids)
+            assert all(grid.nnodes == g.nnodes for g in grids)
 
 
 @requires_exe("gridgen")
-def test_add_refinement_feature(function_tmpdir):
-    bgrid = base_grid()
-    grids = []
+# GRIDGEN seems not to like paths containing "[" or "]", as
+# function_tmpdir does with parametrization, do it manually
+# @pytest.mark.parametrize("grid_type", ["vertex", "unstructured"])
+def test_add_refinement_feature(function_tmpdir):  # , grid_type):
+    bgrid = get_structured_grid()
 
     # test providing refinement feature various ways
-    for features in [
-        [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]],
-        function_tmpdir / "rf0.shp",
-        function_tmpdir / "rf0",
-        "rf0.shp",
-        "rf0",
-    ]:
-        gridgen = Gridgen(bgrid, model_ws=function_tmpdir)
-        gridgen.add_refinement_features(
-            features,
-            "polygon",
-            1,
-            range(bgrid.nlay),
-        )
-        gridgen.build()
-        grid = VertexGrid(**gridgen.get_gridprops_vertexgrid())
-        grid.plot()
-        # plt.show()
-
-        assert grid.nnodes > bgrid.nnodes
-        assert grid.ncpl != bgrid.ncpl
-        assert all(grid.ncpl == g.ncpl for g in grids)
-        assert all(grid.nnodes == g.nnodes for g in grids)
+    for grid_type in ["vertex", "unstructured"]:
+        grids = []
+        for features in [
+            [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]],
+            function_tmpdir / "rf0.shp",
+            function_tmpdir / "rf0",
+            "rf0.shp",
+            "rf0",
+        ]:
+            print(
+                "Testing add_refinement_feature() for",
+                grid_type,
+                "grid with features",
+                features,
+            )
+            gridgen = Gridgen(bgrid, model_ws=function_tmpdir)
+            gridgen.add_refinement_features(
+                features,
+                "polygon",
+                1,
+                range(bgrid.nlay),
+            )
+            gridgen.build()
+            grid = (
+                VertexGrid(**gridgen.get_gridprops_vertexgrid())
+                if grid_type == "vertex"
+                else UnstructuredGrid(
+                    **gridgen.get_gridprops_unstructuredgrid()
+                )
+            )
+            grid.plot()
+            # plt.show()
+
+            assert grid.nnodes > bgrid.nnodes
+            assert not np.array_equal(grid.ncpl, bgrid.ncpl)
+            assert all(np.array_equal(grid.ncpl, g.ncpl) for g in grids)
+            assert all(grid.nnodes == g.nnodes for g in grids)
 
 
 @pytest.mark.slow