diff --git a/autotest/conftest.py b/autotest/conftest.py
index c09eddcf199..8892bc2efc2 100644
--- a/autotest/conftest.py
+++ b/autotest/conftest.py
@@ -100,6 +100,11 @@ def original_regression(request) -> bool:
     return request.config.getoption("--original-regression")
 
 
+@pytest.fixture
+def plot(request) -> bool:
+    return request.config.getoption("--plot")
+
+
 @pytest.fixture(scope="session")
 def markers(pytestconfig) -> str:
     return pytestconfig.getoption("-m")
@@ -124,6 +129,12 @@ def pytest_addoption(parser):
         default=False,
         help="include netcdf test cases",
     )
+    parser.addoption(
+        "--plot",
+        action="store_true",
+        default=False,
+        help="make plots of model output",
+    )
 
 
 def pytest_collection_modifyitems(config, items):
diff --git a/autotest/framework.py b/autotest/framework.py
index 773394c7a31..161d2090f1b 100644
--- a/autotest/framework.py
+++ b/autotest/framework.py
@@ -224,6 +224,7 @@ def __init__(
         api_func: Optional[Callable] = None,
         build: Optional[Callable] = None,
         check: Optional[Callable] = None,
+        plot: Optional[Callable] = None,
         compare: Optional[str] = "auto",
         parallel=False,
         ncpus=1,
@@ -244,6 +245,7 @@ def __init__(
         self.targets = targets
         self.build = build
         self.check = check
+        self.plot = plot
         self.parallel = parallel
         self.ncpus = [ncpus] if isinstance(ncpus, int) else ncpus
         self.api_func = api_func
@@ -753,3 +755,9 @@ def run(self):
             if self.verbose:
                 print("Checking outputs")
             self.check(self)
+
+        # plot results, if enabled
+        if self.plot:
+            if self.verbose:
+                print("Plotting outputs")
+            self.plot(self)
diff --git a/autotest/test_chf_dfw_loop.py b/autotest/test_chf_dfw_loop.py
index c60e8e36f6a..77b227b43d9 100644
--- a/autotest/test_chf_dfw_loop.py
+++ b/autotest/test_chf_dfw_loop.py
@@ -322,8 +322,7 @@ def build_models(idx, test):
     return sim, None
 
 
-def make_plot(test):
-    print("making plots...")
+def plot_output(idx, test):
     import matplotlib.pyplot as plt
 
     name = test.name
@@ -436,49 +435,51 @@ def make_plot(test):
     fname = ws / "loop_network_flow.png"
     plt.savefig(fname)
 
+    # read grb and get locations of the cell center
+    fpth = test.workspace / f"{name}.disv1d.grb"
+    grb = flopy.mf6.utils.MfGrdFile(fpth)
+    cellx = grb._datadict["CELLX"]
+    celly = grb._datadict["CELLY"]
+
     fig = plt.figure(figsize=(6, 6))
     ax = fig.add_subplot(1, 1, 1)
     ax.set_xlim(-100, 6100)
     ax.set_ylim(-100, 6100)
     pmv = flopy.plot.PlotMapView(model=chf, ax=ax)
     pmv.plot_grid()
+
+    # plot and label vertices
     vertices = chf.disv1d.vertices.get_data()
     ax.plot(vertices["xv"], vertices["yv"], "bo")
     for iv, x, y in vertices:
-        ax.text(x, y, f"{iv + 1}")
-    cell1d = chf.disv1d.cell1d.get_data()
-    print(cell1d.dtype)
-    for row in cell1d:
-        ic = row["icell1d"]
-        ncvert = row["ncvert"]
-        if ncvert == 2:
-            n0 = row["icvert_0"]
-            n1 = row["icvert_1"]
-            x0 = vertices["xv"][n0]
-            x1 = vertices["xv"][n1]
-            y0 = vertices["yv"][n0]
-            y1 = vertices["yv"][n1]
-            xm = 0.5 * (x0 + x1)
-            ym = 0.5 * (y0 + y1)
-        elif ncvert == 3:
-            n0 = row["icvert_1"]
-            xm = vertices["xv"][n0] + 150
-            ym = vertices["yv"][n0] + 150
-        ax.text(xm, ym, f"{ic + 1}", color="red")
-    # raise Exception()
-    fname = ws / "grid.png"
+        ax.text(
+            x,
+            y,
+            f"{iv + 1}",
+            color="blue",
+            horizontalalignment="left",
+            verticalalignment="bottom",
+        )
+
+    # plot and label cell centers
+    ax.plot(cellx, celly, marker="o", mfc="none", mec="red", linewidth=0.0)
+    for icell, (x, y) in enumerate(zip(cellx, celly)):
+        ax.text(
+            x,
+            y,
+            f"{icell + 1}",
+            color="red",
+            horizontalalignment="right",
+            verticalalignment="top",
+        )
+
+    fname = ws / "loop_network_grid.png"
     plt.savefig(fname)
 
     return
 
 
 def check_output(idx, test):
-    print("evaluating model...")
-
-    makeplot = False
-    if makeplot:
-        make_plot(test)
-
     # read the observation output
     name = cases[idx]
     fpth = test.workspace / f"{name}.obs.csv"
@@ -532,12 +533,13 @@ def check_output(idx, test):
 
 @pytest.mark.developmode
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
         check=lambda t: check_output(idx, t),
+        plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,
     )
     test.run()