diff --git a/autotest/test_gridgen.py b/autotest/test_gridgen.py
index c76bb5dd7d..e43ad0fba8 100644
--- a/autotest/test_gridgen.py
+++ b/autotest/test_gridgen.py
@@ -1,6 +1,6 @@
 import os
 from pathlib import Path
-from shutil import copyfile, which
+from shutil import which
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -16,7 +16,7 @@
 
 
 @requires_exe("gridgen")
-def test_ctor_accepts_path_or_string(function_tmpdir):
+def test_ctor_accepts_path_or_string_model_ws(function_tmpdir):
     grid = GridCases().structured_small()
 
     g = Gridgen(grid, model_ws=function_tmpdir)
@@ -27,11 +27,14 @@ def test_ctor_accepts_path_or_string(function_tmpdir):
 
 
 def base_grid():
+    """Get a small version of the first grid in:
+    .docs/Notebooks/gridgen_example.py"""
+
     Lx = 100.0
     Ly = 100.0
     nlay = 2
-    nrow = 51
-    ncol = 51
+    nrow = 11
+    ncol = 11
     delr = Lx / ncol
     delc = Ly / nrow
     h0 = 10
@@ -54,37 +57,60 @@ def base_grid():
 
 
 @requires_exe("gridgen")
-def test_add_refinement_feature_accepts_filename(function_tmpdir):
-    # same as 1st grid in .docs/Notebooks/gridgen_example.py
-    grid = base_grid()
-
-    # create initial refined grid (writing refinement
-    # feature shapefile to the workspace)
-    g1 = Gridgen(grid, model_ws=function_tmpdir)
-    g1.add_refinement_features(
+def test_add_active_domain(function_tmpdir):
+    bgrid = base_grid()
+    grids = []
+
+    # test providing active domain various ways
+    for feature in [
         [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]],
-        "polygon",
-        1,
-        range(grid.nlay),
-    )
-    g1.build()
-    g1 = VertexGrid(**g1.get_gridprops_vertexgrid())
-    # g1.plot()
-    # g1.show()
-
-    # recreate gridgen object, using shapefile from the
-    # first run
-    g2 = Gridgen(grid, model_ws=function_tmpdir)
-    g2.add_refinement_features("rf0.shp", "polygon", 1, range(grid.nlay))
-    g2.build()
-    g2 = VertexGrid(**g2.get_gridprops_vertexgrid())
-    # g2.plot()
-    # plt.show()
-
-    assert g1.nnodes > grid.nnodes
-    assert g1.ncpl != grid.ncpl
-    assert g1.ncpl == g2.ncpl
-    assert g1.nnodes == g2.nnodes
+        function_tmpdir / "ad0.shp",
+        "ad0.shp",
+    ]:
+        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)
+
+
+@requires_exe("gridgen")
+def test_add_refinement_feature(function_tmpdir):
+    bgrid = base_grid()
+    grids = []
+
+    # test providing refinement feature various ways
+    for features in [
+        [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]],
+        function_tmpdir / "rf0.shp",
+        "rf0.shp",
+    ]:
+        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)
 
 
 @pytest.mark.slow
diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py
index bccebd6fe4..d32b469a1d 100644
--- a/flopy/utils/gridgen.py
+++ b/flopy/utils/gridgen.py
@@ -6,6 +6,8 @@
 
 import numpy as np
 
+from flopy.utils.flopy_io import relpath_safe
+
 from ..discretization import StructuredGrid
 from ..export.shapefile_utils import shp2recarray
 from ..mbase import resolve_exe
@@ -57,7 +59,7 @@ def features_to_shapefile(
     featuretype : str
         Must be 'point', 'line', 'linestring', or 'polygon'
     filename : str or PathLike
-        Path of the shapefile to write
+        The shapefile to write (extension is optional)
 
     Returns
     -------
@@ -363,15 +365,15 @@ def add_active_domain(self, feature, layers):
         """
         Parameters
         ----------
-        feature : str or list
+        feature : str, path-like or array-like
             feature can be:
-                 a string containing the name of a polygon
-                 a list of polygons
-                 flopy.utils.geometry.Collection object of Polygons
-                 shapely.geometry.Collection object of Polygons
-                 geojson.GeometryCollection object of Polygons
-                 list of shapefile.Shape objects
-                 shapefile.Shapes object
+                a shapefile name (str) or Pathlike
+                a list of polygons
+                a flopy.utils.geometry.Collection object of Polygons
+                a shapely.geometry.Collection object of Polygons
+                a geojson.GeometryCollection object of Polygons
+                a list of shapefile.Shape objects
+                a shapefile.Shapes object
         layers : list
             A list of layers (zero based) for which this active domain
             applies.
@@ -385,37 +387,45 @@ def add_active_domain(self, feature, layers):
         self.nodes = 0
         self.nja = 0
 
-        # Create shapefile or set shapefile to feature
-        adname = f"ad{len(self._addict)}"
-        if isinstance(feature, list):
-            # Create a shapefile
-            adname_w_path = os.path.join(self.model_ws, adname)
-            features_to_shapefile(feature, "polygon", adname_w_path)
-            shapefile = adname
+        # expand shapefile path or create one from polygon feature
+        if isinstance(feature, (str, os.PathLike)):
+            # try expanding absolute path to shapefile
+            shapefile_path = Path(feature).expanduser().absolute()
+            if not shapefile_path.is_file():
+                # look for shapefile in workspace
+                shapefile_path = self.model_ws / feature
+        elif isinstance(feature, list):
+            shapefile_path = self.model_ws / f"ad{len(self._addict)}.shp"
+            features_to_shapefile(feature, "polygon", shapefile_path)
         else:
-            shapefile = feature
+            raise ValueError(
+                f"Feature must be a pathlike (shapefile) or array-like of geometries"
+            )
 
-        self._addict[adname] = shapefile
-        sn = os.path.join(self.model_ws, f"{shapefile}.shp")
-        assert os.path.isfile(sn), f"Shapefile does not exist: {sn}"
+        # make sure shapefile exists
+        assert (
+            shapefile_path.is_file()
+        ), f"Shapefile does not exist: {shapefile_path}"
 
+        # store shapefile info
+        self._addict[shapefile_path.stem] = shapefile_path
         for k in layers:
-            self._active_domain[k] = adname
+            self._active_domain[k] = shapefile_path.stem
 
     def add_refinement_features(self, features, featuretype, level, layers):
         """
         Parameters
         ----------
-        features : str or array-like
+        features : str, path-like or array-like
             features can be
-                a str, the name of a shapefile in the workspace
+                a shapefile name (str) or Pathlike
                 a list of points, lines, or polygons
-                flopy.utils.geometry.Collection object
+                a flopy.utils.geometry.Collection object
                 a list of flopy.utils.geometry objects
-                shapely.geometry.Collection object
-                geojson.GeometryCollection object
+                a shapely.geometry.Collection object
+                a geojson.GeometryCollection object
                 a list of shapefile.Shape objects
-                shapefile.Shapes object
+                a shapefile.Shapes object
         featuretype : str
             Must be either 'point', 'line', or 'polygon'
         level : int
@@ -434,24 +444,33 @@ def add_refinement_features(self, features, featuretype, level, layers):
         self.nja = 0
 
         # Create shapefile or set shapefile to feature
-        rfname = f"rf{len(self._rfdict)}"
-        if isinstance(features, str):
-            shapefile = features
+        if isinstance(features, (str, os.PathLike)):
+            # try expanding absolute path to shapefile
+            shapefile_path = Path(features).expanduser().absolute()
+            if not shapefile_path.is_file():
+                # look for shapefile in workspace
+                shapefile_path = self.model_ws / features
         elif isinstance(features, (list, tuple, np.ndarray)):
-            rfname_w_path = os.path.join(self.model_ws, rfname)
-            features_to_shapefile(features, featuretype, rfname_w_path)
-            shapefile = f"{rfname}.shp"
+            shapefile_path = self.model_ws / f"rf{len(self._rfdict)}.shp"
+            features_to_shapefile(features, featuretype, shapefile_path)
         else:
             raise ValueError(
-                f"Features must be str (shapefile name) or array-like (of geometries)"
+                f"Features must be a pathlike (shapefile) or array-like of geometries"
             )
 
-        self._rfdict[rfname] = [shapefile, featuretype, level]
-        sn = os.path.join(self.model_ws, shapefile)
-        assert os.path.isfile(sn), f"Shapefile does not exist: {sn}"
+        # make sure shapefile exists
+        assert (
+            shapefile_path.is_file()
+        ), f"Shapefile does not exist: {shapefile_path}"
 
+        # store shapefile info
+        self._rfdict[shapefile_path.stem] = [
+            shapefile_path,
+            featuretype,
+            level,
+        ]
         for k in layers:
-            self._refinement_features[k].append(rfname)
+            self._refinement_features[k].append(shapefile_path.stem)
 
     def build(self, verbose=False):
         """