diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
index 4a45cb107..14e98076a 100644
--- a/.github/workflows/benchmark.yml
+++ b/.github/workflows/benchmark.yml
@@ -57,7 +57,6 @@ jobs:
         shell: bash -l {0}
         run: |
           pip install --upgrade pip
-          pip install https://github.com/modflowpy/pymake/zipball/master
           pip install xmipy
           pip install .
 
diff --git a/.github/workflows/commit.yml b/.github/workflows/commit.yml
index e9562c926..80f1ad968 100644
--- a/.github/workflows/commit.yml
+++ b/.github/workflows/commit.yml
@@ -194,7 +194,6 @@ jobs:
         if: runner.os == 'Windows'
         run: |
           pip install --upgrade pip
-          pip install https://github.com/modflowpy/pymake/zipball/master
           pip install xmipy
           pip install .
 
diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml
index b6dc20545..4cb4771eb 100644
--- a/.github/workflows/examples.yml
+++ b/.github/workflows/examples.yml
@@ -56,7 +56,6 @@ jobs:
         shell: bash -l {0}
         run: |
           pip install --upgrade pip
-          pip install https://github.com/modflowpy/pymake/zipball/master
           pip install xmipy
           pip install .
 
diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml
index cdbcdf4c8..349e8df22 100644
--- a/.github/workflows/regression.yml
+++ b/.github/workflows/regression.yml
@@ -56,7 +56,6 @@ jobs:
         shell: bash -l {0}
         run: |
           pip install --upgrade pip
-          pip install https://github.com/modflowpy/pymake/zipball/master
           pip install xmipy
           pip install .
 
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 6faebce7e..8885eb3ca 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -66,10 +66,12 @@ Before you submit your Pull Request (PR) consider the following guidelines:
 
     ```shell
     cd autotest
-    pytest -v ci_prepare.py
-    pytest -v
+    pytest -v -n auto
     ```
-   Note: the FloPy test suite requires the [pytest](https://pypi.org/project/pytest/) and [pymake](https://github.com/modflowpy/pymake) python packages. All the FloPy dependencies must also be installed for the tests to pass.
+   
+   **Note**: the FloPy test suite requires the [pytest](https://pypi.org/project/pytest/) and [modflow-devtools](https://github.com/MODFLOW-USGS/modflow-devtools) python packages. Be sure to install all optional dependencies as described in the [developer docs](DEVELOPER.md), as tests will be skipped if optional dependencies are not found.
+
+    **Note:** you will either need to exclude notebooks when checking in your changeset or restore them after running the full test suite, to avoid including large volumes of execution metadata with your PR.
 
 7. Commit your changes using a descriptive commit message that follows our
   [commit message conventions](#commit). Adherence to these conventions
diff --git a/autotest/regression/test_lgr.py b/autotest/regression/test_lgr.py
index 45d20fcad..2373a3149 100644
--- a/autotest/regression/test_lgr.py
+++ b/autotest/regression/test_lgr.py
@@ -7,16 +7,14 @@
 from modflow_devtools.markers import requires_exe, requires_pkg
 
 import flopy
+from flopy.utils.compare import compare_heads
 
 
 @flaky
 @requires_exe("mflgr")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test_simplelgr(function_tmpdir, example_data_path):
     """Test load and write of distributed MODFLOW-LGR example problem."""
-    import pymake
-
     mflgr_v2_ex3_path = example_data_path / "mflgr_v2" / "ex3"
 
     ws = function_tmpdir / mflgr_v2_ex3_path.stem
@@ -66,12 +64,12 @@ def test_simplelgr(function_tmpdir, example_data_path):
     print("compare parent results")
     pth0 = join(ws, "ex3_parent.nam")
     pth1 = join(model_ws2, "ex3_parent.nam")
-    success = pymake.compare_heads(pth0, pth1)
+    success = compare_heads(pth0, pth1)
     assert success, "parent heads do not match"
 
     # compare child results
     print("compare child results")
     pth0 = join(ws, "ex3_child.nam")
     pth1 = join(model_ws2, "ex3_child.nam")
-    success = pymake.compare_heads(pth0, pth1)
+    success = compare_heads(pth0, pth1)
     assert success, "child heads do not match"
diff --git a/autotest/regression/test_mf6.py b/autotest/regression/test_mf6.py
index e7b65fbd4..fb7988b02 100644
--- a/autotest/regression/test_mf6.py
+++ b/autotest/regression/test_mf6.py
@@ -48,17 +48,15 @@
 from flopy.mf6.mfbase import FlopyException, MFDataException
 from flopy.mf6.utils import testutils
 from flopy.utils import CellBudgetFile
+from flopy.utils.compare import compare_concentrations, compare_heads
 from flopy.utils.datautil import PyListUtil
 
 pytestmark = pytest.mark.mf6
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test_np001(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "np001"
     model_name = "np001_mod"
@@ -407,7 +405,7 @@ def test_np001(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = os.path.join(ws, "np001_mod 1.hds")
     outfile = os.path.join(ws, "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -451,7 +449,7 @@ def test_np001(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = os.path.join(run_folder_new, "np001_mod 1.hds")
     outfile = os.path.join(run_folder_new, "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -626,11 +624,8 @@ def test_np001(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test_np002(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "np002"
     model_name = "np002_mod"
@@ -811,7 +806,7 @@ def test_np002(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = os.path.join(ws, "np002_mod.hds")
     outfile = os.path.join(ws, "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -903,11 +898,8 @@ def test_np002(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test021_twri(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test021_twri"
     model_name = "twri"
@@ -1116,7 +1108,7 @@ def test021_twri(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = os.path.join(ws, "twri.hds")
     outfile = os.path.join(ws, "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -1129,12 +1121,9 @@ def test021_twri(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test005_create_tests_advgw_tidal(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test005_advgw_tidal"
     model_name = "AdvGW_tidal"
@@ -1674,7 +1663,7 @@ def test005_create_tests_advgw_tidal(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "AdvGW_tidal.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -1758,11 +1747,8 @@ def test005_create_tests_advgw_tidal(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test004_create_tests_bcfss(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test004_bcfss"
     model_name = "bcf2ss"
@@ -1943,7 +1929,7 @@ def test004_create_tests_bcfss(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = os.path.join(str(function_tmpdir), "bcf2ss.hds")
     outfile = os.path.join(str(function_tmpdir), "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -1956,11 +1942,8 @@ def test004_create_tests_bcfss(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test035_create_tests_fhb(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test035_fhb"
     model_name = "fhb2015"
@@ -2086,7 +2069,7 @@ def test035_create_tests_fhb(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "fhb2015_fhb.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -2099,11 +2082,9 @@ def test035_create_tests_fhb(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake", "shapefile")
+@requires_pkg("shapefile")
 @pytest.mark.regression
 def test006_create_tests_gwf3_disv(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test006_gwf3_disv"
     model_name = "flow"
@@ -2374,7 +2355,7 @@ def test006_create_tests_gwf3_disv(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "flow.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -2393,11 +2374,8 @@ def test006_create_tests_gwf3_disv(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test006_create_tests_2models_gnc(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test006_2models_gnc"
     model_name_1 = "model1"
@@ -2682,7 +2660,7 @@ def test006_create_tests_2models_gnc(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "model1.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_1,
@@ -2693,7 +2671,7 @@ def test006_create_tests_2models_gnc(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "model2.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_2,
@@ -2733,12 +2711,9 @@ def test006_create_tests_2models_gnc(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test050_create_tests_circle_island(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test050_circle_island"
     model_name = "ci"
@@ -2821,7 +2796,7 @@ def test050_create_tests_circle_island(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "ci.output.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -2834,15 +2809,12 @@ def test050_create_tests_circle_island(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.xfail(
     reason="possible python3.7/windows incompatibilities in testutils.read_std_array "
     "https://github.com/modflowpy/flopy/runs/7581629193?check_suite_focus=true#step:11:1753"
 )
 @pytest.mark.regression
 def test028_create_tests_sfr(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test028_sfr"
     model_name = "test1tr"
@@ -3095,7 +3067,7 @@ def test028_create_tests_sfr(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "test1tr.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -3109,11 +3081,8 @@ def test028_create_tests_sfr(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test_create_tests_transport(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test_transport"
     name = "mst03"
@@ -3328,7 +3297,7 @@ def test_create_tests_transport(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "gwf_mst03.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file,
@@ -3336,7 +3305,7 @@ def test_create_tests_transport(function_tmpdir, example_data_path):
         outfile=outfile,
     )
     conc_new = str(function_tmpdir / "gwt_mst03.ucn")
-    assert pymake.compare_concs(
+    assert compare_concentrations(
         None,
         None,
         files1=expected_conc_file,
@@ -3349,12 +3318,10 @@ def test_create_tests_transport(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake", "shapely")
+@requires_pkg("shapely")
 @pytest.mark.slow
 @pytest.mark.regression
 def test001a_tharmonic(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test001a_Tharmonic"
     model_name = "flow15"
@@ -3416,7 +3383,7 @@ def test001a_tharmonic(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = str(function_tmpdir / "flow15_flow.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None, None, files1=expected_head_file_a, files2=head_new
     )
 
@@ -3470,7 +3437,7 @@ def test001a_tharmonic(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = os.path.join(str(save_folder), "flow15_flow.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None, None, files1=expected_head_file_b, files2=head_new
     )
 
@@ -3481,11 +3448,8 @@ def test001a_tharmonic(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test003_gwfs_disv(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test003_gwfs_disv"
     model_name = "gwf_1"
@@ -3521,7 +3485,7 @@ def test003_gwfs_disv(function_tmpdir, example_data_path):
     )
 
     head_new = os.path.join(str(function_tmpdir), "model.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None, None, files1=expected_head_file_a, files2=head_new
     )
 
@@ -3563,7 +3527,7 @@ def test003_gwfs_disv(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = os.path.join(str(save_folder), "model.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None, None, files1=expected_head_file_b, files2=head_new
     )
 
@@ -3574,12 +3538,9 @@ def test003_gwfs_disv(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test005_advgw_tidal(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test005_advgw_tidal"
     model_name = "gwf_1"
@@ -3632,7 +3593,7 @@ def test005_advgw_tidal(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = os.path.join(str(function_tmpdir), "advgw_tidal.hds")
     outfile = os.path.join(str(function_tmpdir), "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_a,
@@ -3642,11 +3603,8 @@ def test005_advgw_tidal(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test006_gwf3(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test006_gwf3"
     model_name = "gwf_1"
@@ -3704,7 +3662,7 @@ def test006_gwf3(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = os.path.join(str(function_tmpdir), "flow.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_a,
@@ -3751,7 +3709,7 @@ def test006_gwf3(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = os.path.join(save_folder, "flow.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_b,
@@ -3799,7 +3757,7 @@ def test006_gwf3(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = os.path.join(save_folder, "flow.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_b,
@@ -3827,11 +3785,8 @@ def test006_gwf3(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test045_lake1ss_table(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test045_lake1ss_table"
     model_name = "lakeex1b"
@@ -3865,7 +3820,7 @@ def test045_lake1ss_table(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "lakeex1b.hds")
     outfile = str(function_tmpdir / "headcompare_a.txt")
-    success = pymake.compare_heads(
+    success = compare_heads(
         None,
         None,
         files1=expected_head_file_a,
@@ -3895,7 +3850,7 @@ def test045_lake1ss_table(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(save_folder / "lakeex1b.hds")
     outfile = str(function_tmpdir / "headcompare_b.txt")
-    success = pymake.compare_heads(
+    success = compare_heads(
         None,
         None,
         files1=expected_head_file_b,
@@ -3906,12 +3861,9 @@ def test045_lake1ss_table(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test006_2models_mvr(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test006_2models_mvr"
     sim_name = "test006_2models_mvr"
@@ -3946,7 +3898,7 @@ def test006_2models_mvr(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = str(ws / "model1.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_a,
@@ -3954,7 +3906,7 @@ def test006_2models_mvr(function_tmpdir, example_data_path):
     )
 
     head_new = str(ws / "model2.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_aa,
@@ -4033,7 +3985,7 @@ def test006_2models_mvr(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = os.path.join(save_folder, "model1.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_b,
@@ -4041,7 +3993,7 @@ def test006_2models_mvr(function_tmpdir, example_data_path):
     )
 
     head_new = os.path.join(save_folder, "model2.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_bb,
@@ -4095,7 +4047,6 @@ def test006_2models_mvr(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test001e_uzf_3lay(function_tmpdir, example_data_path):
@@ -4196,12 +4147,9 @@ def test001e_uzf_3lay(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test045_lake2tr(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test045_lake2tr"
     model_name = "lakeex2a"
@@ -4225,7 +4173,7 @@ def test045_lake2tr(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = str(function_tmpdir / "lakeex2a.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_a,
@@ -4261,7 +4209,7 @@ def test045_lake2tr(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = str(save_folder / "lakeex2a.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_b,
@@ -4271,11 +4219,8 @@ def test045_lake2tr(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test036_twrihfb(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test036_twrihfb"
     model_name = "twrihfb2015"
@@ -4304,7 +4249,7 @@ def test036_twrihfb(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = str(function_tmpdir / "twrihfb2015_output.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_a,
@@ -4346,7 +4291,7 @@ def test036_twrihfb(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = str(save_folder / "twrihfb2015_output.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_b,
@@ -4355,12 +4300,9 @@ def test036_twrihfb(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test027_timeseriestest(function_tmpdir, example_data_path):
-    import pymake
-
     # init paths
     test_ex_name = "test027_TimeseriesTest"
     model_name = "gwf_1"
@@ -4396,7 +4338,7 @@ def test027_timeseriestest(function_tmpdir, example_data_path):
     # compare output to expected results
     head_new = str(function_tmpdir / "timeseriestest.hds")
     outfile = str(function_tmpdir / "head_compare.dat")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_a,
@@ -4425,7 +4367,7 @@ def test027_timeseriestest(function_tmpdir, example_data_path):
 
     # compare output to expected results
     head_new = os.path.join(str(save_folder), "timeseriestest.hds")
-    assert pymake.compare_heads(
+    assert compare_heads(
         None,
         None,
         files1=expected_head_file_b,
diff --git a/autotest/regression/test_mf6_examples.py b/autotest/regression/test_mf6_examples.py
index 664ffb704..9a14a365f 100644
--- a/autotest/regression/test_mf6_examples.py
+++ b/autotest/regression/test_mf6_examples.py
@@ -6,12 +6,12 @@
 from modflow_devtools.markers import requires_exe, requires_pkg
 
 from flopy.mf6 import MFSimulation
+from flopy.utils.compare import compare_heads
 
 pytestmark = pytest.mark.mf6
 
 
 @requires_exe("mf6")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test_mf6_example_simulations(function_tmpdir, mf6_example_namfiles):
@@ -23,8 +23,6 @@ def test_mf6_example_simulations(function_tmpdir, mf6_example_namfiles):
     # function_tmpdir: function-scoped temporary directory fixture
     # mf6_example_namfiles: ordered list of namfiles for 1+ coupled models
 
-    import pymake
-
     # make sure we have at least 1 name file
     if len(mf6_example_namfiles) == 0:
         pytest.skip("No namfiles (expected ordered collection)")
@@ -82,7 +80,7 @@ def run_models():
             headfiles2 = [p for p in cmpdir.glob("*.hds")]
 
             # compare heads
-            assert pymake.compare_heads(
+            assert compare_heads(
                 None,
                 None,
                 precision="double",
diff --git a/autotest/regression/test_mfnwt.py b/autotest/regression/test_mfnwt.py
index 7b2ac4ee2..d2b48bf53 100644
--- a/autotest/regression/test_mfnwt.py
+++ b/autotest/regression/test_mfnwt.py
@@ -2,10 +2,11 @@
 
 import pytest
 from autotest.conftest import get_example_data_path
-from modflow_devtools.markers import requires_exe, requires_pkg
+from modflow_devtools.markers import requires_exe
 
 from flopy.modflow import Modflow, ModflowNwt, ModflowUpw
 from flopy.utils import parsenamefile
+from flopy.utils.compare import compare_budget, compare_heads
 
 
 def get_nfnwt_namfiles():
@@ -28,13 +29,10 @@ def get_nfnwt_namfiles():
 
 
 @requires_exe("mfnwt")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 @pytest.mark.parametrize("namfile", get_nfnwt_namfiles())
 def test_run_mfnwt_model(function_tmpdir, namfile):
-    import pymake
-
     # load a MODFLOW-2005 model, convert to a MFNWT model,
     # write it back out, run the MFNWT model, load the MFNWT model,
     # and compare the results.
@@ -136,11 +134,9 @@ def test_run_mfnwt_model(function_tmpdir, namfile):
     fn1 = os.path.join(pthf, namfile)
 
     fsum = str(function_tmpdir / f"{base_name}.head.out")
-    assert pymake.compare_heads(
-        fn0, fn1, outfile=fsum
-    ), "head comparison failure"
+    assert compare_heads(fn0, fn1, outfile=fsum), "head comparison failure"
 
     fsum = str(function_tmpdir / f"{base_name}.budget.out")
-    assert pymake.compare_budget(
+    assert compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     ), "budget comparison failure"
diff --git a/autotest/regression/test_modflow.py b/autotest/regression/test_modflow.py
index 4daf494d7..2dc957a08 100644
--- a/autotest/regression/test_modflow.py
+++ b/autotest/regression/test_modflow.py
@@ -8,6 +8,7 @@
 from modflow_devtools.markers import requires_exe, requires_pkg
 
 from flopy.modflow import Modflow, ModflowOc
+from flopy.utils.compare import compare_budget, compare_heads
 
 
 @pytest.fixture
@@ -21,12 +22,9 @@ def uzf_example_path(example_data_path):
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test_uzf_unit_numbers(function_tmpdir, uzf_example_path):
-    import pymake
-
     mfnam = "UZFtest2.nam"
     ws = str(function_tmpdir / "ws")
     copytree(uzf_example_path, ws)
@@ -76,19 +74,16 @@ def test_uzf_unit_numbers(function_tmpdir, uzf_example_path):
 
     # compare budget terms
     fsum = join(str(function_tmpdir), f"{splitext(mfnam)[0]}.budget.out")
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
     assert success, "budget comparison failure"
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test_unitnums(function_tmpdir, mf2005_test_path):
-    import pymake
-
     mfnam = "testsfr2_tab.nam"
     ws = str(function_tmpdir / "ws")
     copytree(mf2005_test_path, ws)
@@ -117,22 +112,19 @@ def test_unitnums(function_tmpdir, mf2005_test_path):
     fn1 = join(model_ws2, mfnam)
 
     fsum = join(ws, f"{splitext(mfnam)[0]}.budget.out")
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
     assert success, "budget comparison failure"
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 def test_gage(function_tmpdir, example_data_path):
     """
     test043 load and write of MODFLOW-2005 GAGE example problem
     """
-    import pymake
-
     pth = str(example_data_path / "mf2005_test")
     fpth = join(pth, "testsfr2_tab.nam")
     ws = str(function_tmpdir / "ws")
@@ -172,7 +164,6 @@ def test_gage(function_tmpdir, example_data_path):
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 @pytest.mark.parametrize(
@@ -183,8 +174,6 @@ def test_gage(function_tmpdir, example_data_path):
     ],
 )
 def test_mf2005pcgn(function_tmpdir, namfile):
-    import pymake
-
     ws = function_tmpdir / "ws"
     copytree(Path(namfile).parent, ws)
     nf = Path(namfile).name
@@ -215,26 +204,23 @@ def test_mf2005pcgn(function_tmpdir, namfile):
     fn1 = str(ws2 / nf)
 
     fsum = str(function_tmpdir / f"{Path(namfile).stem}.head.out")
-    success = pymake.compare_heads(fn0, fn1, outfile=fsum, htol=0.005)
+    success = compare_heads(fn0, fn1, outfile=fsum, htol=0.005)
     assert success, "head comparison failure"
 
     fsum = str(function_tmpdir / f"{Path(namfile).stem}.budget.out")
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
     assert success, "budget comparison failure"
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 @pytest.mark.parametrize(
     "namfile", [str(__example_data_path / "secp" / nf) for nf in ["secp.nam"]]
 )
 def test_mf2005gmg(function_tmpdir, namfile):
-    import pymake
-
     ws = function_tmpdir / "ws"
     copytree(Path(namfile).parent, ws)
     nf = Path(namfile).name
@@ -260,18 +246,17 @@ def test_mf2005gmg(function_tmpdir, namfile):
     fn1 = str(function_tmpdir / nf)
 
     fsum = str(function_tmpdir / f"{Path(namfile).stem}.head.out")
-    success = pymake.compare_heads(fn0, fn1, outfile=fsum)
+    success = compare_heads(fn0, fn1, outfile=fsum)
     assert success, "head comparison failure"
 
     fsum = str(function_tmpdir / f"{Path(namfile).stem}.budget.out")
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
     assert success, "budget comparison failure"
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.regression
 @pytest.mark.parametrize(
     "namfile",
@@ -281,8 +266,6 @@ def test_mf2005(function_tmpdir, namfile):
     """
     test045 load and write of MODFLOW-2005 GMG example problem
     """
-    import pymake
-
     compth = function_tmpdir / "flopy"
     ws = function_tmpdir / "ws"
     copytree(Path(namfile).parent, str(ws))
@@ -325,17 +308,17 @@ def test_mf2005(function_tmpdir, namfile):
 
     # compare heads
     fsum = str(ws / f"{Path(namfile).stem}.head.out")
-    success = pymake.compare_heads(fn0, fn1, outfile=fsum)
+    success = compare_heads(fn0, fn1, outfile=fsum)
     assert success, "head comparison failure"
 
     # compare heads
     fsum = str(ws / f"{Path(namfile).stem}.ddn.out")
-    success = pymake.compare_heads(fn0, fn1, outfile=fsum, text="drawdown")
+    success = compare_heads(fn0, fn1, outfile=fsum, text="drawdown")
     assert success, "head comparison failure"
 
     # compare budgets
     fsum = str(ws / f"{Path(namfile).stem}.budget.out")
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
     assert success, "budget comparison failure"
@@ -354,13 +337,10 @@ def test_mf2005(function_tmpdir, namfile):
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 @pytest.mark.parametrize("namfile", mf2005_namfiles)
 def test_mf2005fhb(function_tmpdir, namfile):
-    import pymake
-
     ws = str(function_tmpdir / "ws")
     copytree(Path(namfile).parent, ws)
 
@@ -382,24 +362,21 @@ def test_mf2005fhb(function_tmpdir, namfile):
     fn1 = join(str(function_tmpdir), Path(namfile).name)
 
     fsum = join(ws, f"{Path(namfile).stem}.head.out")
-    success = pymake.compare_heads(fn0, fn1, outfile=fsum)
+    success = compare_heads(fn0, fn1, outfile=fsum)
     assert success, "head comparison failure"
 
     fsum = join(ws, f"{Path(namfile).stem}.budget.out")
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
     assert success, "budget comparison failure"
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 @pytest.mark.parametrize("namfile", mf2005_namfiles)
 def test_mf2005_lake(function_tmpdir, namfile, mf2005_test_path):
-    import pymake
-
     ws = str(function_tmpdir / "ws")
 
     copytree(mf2005_test_path, ws)
@@ -433,7 +410,7 @@ def test_mf2005_lake(function_tmpdir, namfile, mf2005_test_path):
 
     fsum = join(ws, f"{Path(namfile).stem}.budget.out")
 
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
     assert success, "budget comparison failure"
diff --git a/autotest/regression/test_str.py b/autotest/regression/test_str.py
index 0ad095c28..3758cf07d 100644
--- a/autotest/regression/test_str.py
+++ b/autotest/regression/test_str.py
@@ -2,6 +2,7 @@
 from modflow_devtools.markers import requires_exe, requires_pkg
 
 from flopy.modflow import Modflow, ModflowOc, ModflowStr
+from flopy.utils.compare import compare_heads
 
 str_items = {
     0: {
@@ -13,11 +14,8 @@
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test_str_fixed_free(function_tmpdir, example_data_path):
-    import pymake
-
     mf2005_model_path = example_data_path / "mf2005_test"
 
     m = Modflow.load(
@@ -123,6 +121,6 @@ def test_str_fixed_free(function_tmpdir, example_data_path):
     # compare the fixed and free format head files
     fn1 = str(function_tmpdir / "str.nam")
     fn2 = str(function_tmpdir / "str.nam")
-    assert pymake.compare_heads(
+    assert compare_heads(
         fn1, fn2, verbose=True
     ), "fixed and free format input output head files are different"
diff --git a/autotest/regression/test_swi2.py b/autotest/regression/test_swi2.py
index 8a37df8af..045cfc060 100644
--- a/autotest/regression/test_swi2.py
+++ b/autotest/regression/test_swi2.py
@@ -2,9 +2,10 @@
 import shutil
 
 import pytest
-from modflow_devtools.markers import requires_exe, requires_pkg
+from modflow_devtools.markers import requires_exe
 
 from flopy.modflow import Modflow
+from flopy.utils.compare import compare_budget
 
 
 @pytest.fixture
@@ -13,15 +14,12 @@ def swi_path(example_data_path):
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.slow
 @pytest.mark.regression
 @pytest.mark.parametrize(
     "namfile", ["swiex1.nam", "swiex2_strat.nam", "swiex3.nam"]
 )
 def test_mf2005swi2(function_tmpdir, swi_path, namfile):
-    import pymake
-
     name = namfile.replace(".nam", "")
     ws = str(function_tmpdir / "ws")
     shutil.copytree(swi_path, ws)
@@ -49,7 +47,7 @@ def test_mf2005swi2(function_tmpdir, swi_path, namfile):
     fn1 = os.path.join(model_ws2, namfile)
 
     fsum = os.path.join(ws, f"{os.path.splitext(namfile)[0]}.budget.out")
-    success = pymake.compare_budget(
+    success = compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     )
 
diff --git a/autotest/regression/test_wel.py b/autotest/regression/test_wel.py
index 0af16e34d..4a9af3dad 100644
--- a/autotest/regression/test_wel.py
+++ b/autotest/regression/test_wel.py
@@ -13,14 +13,12 @@
     ModflowPcg,
     ModflowWel,
 )
+from flopy.utils.compare import compare_budget, compare_heads
 
 
 @requires_exe("mf2005")
-@requires_pkg("pymake")
 @pytest.mark.regression
 def test_binary_well(function_tmpdir):
-    import pymake
-
     nlay = 3
     nrow = 3
     ncol = 3
@@ -102,13 +100,11 @@ def test_binary_well(function_tmpdir):
     fsum = os.path.join(
         str(function_tmpdir), f"{os.path.splitext(mfnam)[0]}.head.out"
     )
-    assert pymake.compare_heads(
-        fn0, fn1, outfile=fsum
-    ), "head comparison failure"
+    assert compare_heads(fn0, fn1, outfile=fsum), "head comparison failure"
 
     fsum = os.path.join(
         str(function_tmpdir), f"{os.path.splitext(mfnam)[0]}.budget.out"
     )
-    assert pymake.compare_budget(
+    assert compare_budget(
         fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fsum
     ), "budget comparison failure"
diff --git a/docs/flopy_method_dependencies.md b/docs/flopy_method_dependencies.md
index 4cdf70fef..27a3025b8 100644
--- a/docs/flopy_method_dependencies.md
+++ b/docs/flopy_method_dependencies.md
@@ -1,35 +1,35 @@
 Additional dependencies to use optional FloPy helper methods are listed below.
 
-| Method                                                                               | Python Package                                                     |
-| ------------------------------------------------------------------------------------ | ------------------------------------------------------------------ |
-| `.plot_shapefile()`                                                                  | **Pyshp** >= 2.0.0                                                 |
-| `.to_shapefile()`                                                                    | **Pyshp** >= 2.0.0                                                 |
-| `.export(*.shp)`                                                                     | **Pyshp** >= 2.0.0                                                 |
-| `.export(*.nc)`                                                                      | **netcdf4** >= 1.1, and **python-dateutil** >= 2.4.0               |
-| `.export(*.tif)`                                                                     | **rasterio**                                                       |
-| `.export(*.asc)` in `flopy.utils.reference` `SpatialReference` class                 | **scipy.ndimage**                                                  |
-| `.interpolate()` in `flopy.utils.reference` `SpatialReference` class                 | **scipy.interpolate**                                              |
-| `.interpolate()` in `flopy.mf6.utils.reference` `StructuredSpatialReference` class   | **scipy.interpolate**                                              |
-| `._parse_units_from_proj4()` in `flopy.utils.reference` `SpatialReference` class     | **pyproj**                                                         |
-| `.get_dataframes()` in `flopy.utils.mflistfile` `ListBudget` class                   | **pandas** >= 0.15.0                                               |
-| `.get_dataframes()` in `flopy.utils.observationfile` `ObsFiles` class                | **pandas** >= 0.15.0                                               |
-| `.get_dataframes()` in `flopy.utils.sfroutputfile` `ModflowSfr2` class               | **pandas** >= 0.15.0                                               |
-| `.get_dataframes()` in `flopy.utils.util_list` `MfList` class                        | **pandas** >= 0.15.0                                               |
-| `.get_dataframes()` in `flopy.utils.zonebud` `ZoneBudget` class                      | **pandas** >= 0.15.0                                               |
-| `.pivot_keyarray()` in `flopy.mf6.utils.arrayutils` `AdvancedPackageUtil` class      | **pandas** >= 0.15.0                                               |
-| `._get_vertices()` in `flopy.mf6.utils.binaryfile_utils` `MFOutputRequester` class   | **pandas** >= 0.15.0                                               |
-| `.get_dataframe()` in `flopy.mf6.utils.mfobservation` `Observations` class           | **pandas** >= 0.15.0                                               |
-| `.df()` in `flopy.modflow.mfsfr2` `SfrFile` class                                    | **pandas** >= 0.15.0                                               |
-| `.time_coverage()` in `flopy.export.metadata` `acc` class - ***used if available***  | **pandas** >= 0.15.0                                               |
-| `.loadtxt()` in `flopy.utils.flopyio` - ***used if available***                      | **pandas** >= 0.15.0                                               |
-| `.generate_classes()` in `flopy.mf6.utils`                                           | [**pymake**](https://github.com/modflowpy/pymake)                  |
-| `GridIntersect()` in `flopy.utils.gridintersect`                                     | **shapely**                                                        |
-| `GridIntersect().plot_polygon()` in `flopy.utils.gridintersect`                      | **shapely** and **descartes**                                      |
-| `Raster()` in `flopy.utils.Raster`                                                   | **rasterio**, **rasterstats**, **affine**, and **scipy**                            |
-| `Raster().sample_polygon()` in `flopy.utils.Raster`                                  | **shapely**                                                        |
-| `Raster().crop()` in `flopy.utils.Raster`                                            | **shapely**                                                        |
-| `.array_at_verts()` in `flopy.discretization.structuredgrid` `StructuredGrid` class  | **scipy.interpolate**                                              |
-| `get_sciencebase_xml_metadata()` in `flopy.export.metadata` `acdd` class             | **defusedxml**                                                     |
-| `flopy.utils.geospatial_utils` `GeoSpatialUtil` class                                | **geojson**                                                        |
-| `flopy.utils.geospatial_utils` `GeoSpatialCollection` class                          | **geojson**                                                        |
-| `flopy.export.vtk` `Vtk` class                                                       | **vtk**                                                            |
+| Method                                                                               | Python Package                                                           |
+| ------------------------------------------------------------------------------------ |--------------------------------------------------------------------------|
+| `.plot_shapefile()`                                                                  | **Pyshp** >= 2.0.0                                                       |
+| `.to_shapefile()`                                                                    | **Pyshp** >= 2.0.0                                                       |
+| `.export(*.shp)`                                                                     | **Pyshp** >= 2.0.0                                                       |
+| `.export(*.nc)`                                                                      | **netcdf4** >= 1.1, and **python-dateutil** >= 2.4.0                     |
+| `.export(*.tif)`                                                                     | **rasterio**                                                             |
+| `.export(*.asc)` in `flopy.utils.reference` `SpatialReference` class                 | **scipy.ndimage**                                                        |
+| `.interpolate()` in `flopy.utils.reference` `SpatialReference` class                 | **scipy.interpolate**                                                    |
+| `.interpolate()` in `flopy.mf6.utils.reference` `StructuredSpatialReference` class   | **scipy.interpolate**                                                    |
+| `._parse_units_from_proj4()` in `flopy.utils.reference` `SpatialReference` class     | **pyproj**                                                               |
+| `.get_dataframes()` in `flopy.utils.mflistfile` `ListBudget` class                   | **pandas** >= 0.15.0                                                     |
+| `.get_dataframes()` in `flopy.utils.observationfile` `ObsFiles` class                | **pandas** >= 0.15.0                                                     |
+| `.get_dataframes()` in `flopy.utils.sfroutputfile` `ModflowSfr2` class               | **pandas** >= 0.15.0                                                     |
+| `.get_dataframes()` in `flopy.utils.util_list` `MfList` class                        | **pandas** >= 0.15.0                                                     |
+| `.get_dataframes()` in `flopy.utils.zonebud` `ZoneBudget` class                      | **pandas** >= 0.15.0                                                     |
+| `.pivot_keyarray()` in `flopy.mf6.utils.arrayutils` `AdvancedPackageUtil` class      | **pandas** >= 0.15.0                                                     |
+| `._get_vertices()` in `flopy.mf6.utils.binaryfile_utils` `MFOutputRequester` class   | **pandas** >= 0.15.0                                                     |
+| `.get_dataframe()` in `flopy.mf6.utils.mfobservation` `Observations` class           | **pandas** >= 0.15.0                                                     |
+| `.df()` in `flopy.modflow.mfsfr2` `SfrFile` class                                    | **pandas** >= 0.15.0                                                     |
+| `.time_coverage()` in `flopy.export.metadata` `acc` class - ***used if available***  | **pandas** >= 0.15.0                                                     |
+| `.loadtxt()` in `flopy.utils.flopyio` - ***used if available***                      | **pandas** >= 0.15.0                                                     |
+| `.generate_classes()` in `flopy.mf6.utils`                                           | [**modflow-devtools**](https://github.com/MODFLOW-USGS/modflow-devtools) |
+| `GridIntersect()` in `flopy.utils.gridintersect`                                     | **shapely**                                                              |
+| `GridIntersect().plot_polygon()` in `flopy.utils.gridintersect`                      | **shapely** and **descartes**                                            |
+| `Raster()` in `flopy.utils.Raster`                                                   | **rasterio**, **rasterstats**, **affine**, and **scipy**                 |
+| `Raster().sample_polygon()` in `flopy.utils.Raster`                                  | **shapely**                                                              |
+| `Raster().crop()` in `flopy.utils.Raster`                                            | **shapely**                                                              |
+| `.array_at_verts()` in `flopy.discretization.structuredgrid` `StructuredGrid` class  | **scipy.interpolate**                                                    |
+| `get_sciencebase_xml_metadata()` in `flopy.export.metadata` `acdd` class             | **defusedxml**                                                           |
+| `flopy.utils.geospatial_utils` `GeoSpatialUtil` class                                | **geojson**                                                              |
+| `flopy.utils.geospatial_utils` `GeoSpatialCollection` class                          | **geojson**                                                              |
+| `flopy.export.vtk` `Vtk` class                                                       | **vtk**                                                                  |
diff --git a/flopy/utils/compare.py b/flopy/utils/compare.py
index c64b98e6d..1cbf4a347 100644
--- a/flopy/utils/compare.py
+++ b/flopy/utils/compare.py
@@ -131,7 +131,7 @@ def compare_budget(
     lst_file1 = None
     if files1 is None:
         lst_file = get_entries_from_namefile(namefile1, "list")
-        lst_file1 = lst_file[0][0]
+        lst_file1 = lst_file[0][0] if any(lst_file) else None
     else:
         if isinstance(files1, str):
             files1 = [files1]
@@ -145,7 +145,7 @@ def compare_budget(
     lst_file2 = None
     if files2 is None:
         lst_file = get_entries_from_namefile(namefile2, "list")
-        lst_file2 = lst_file[0][0]
+        lst_file2 = lst_file[0][0] if any(lst_file) else None
     else:
         if isinstance(files2, str):
             files2 = [files2]
@@ -166,7 +166,6 @@ def compare_budget(
     # Open output file
     if outfile is not None:
         f = open(outfile, "w")
-        f.write("Created by pymake.autotest.compare\n")
 
     # Initialize SWR budget objects
     lst1obj = flopy.utils.MfusgListBudget(lst_file1)
@@ -337,7 +336,7 @@ def compare_swrbudget(
     list1 = None
     if files1 is None:
         lst = get_entries_from_namefile(namefile1, "list")
-        list1 = lst[0][0]
+        list1 = lst[0][0] if any(lst) else None
     else:
         for file in files1:
             if (
@@ -349,7 +348,7 @@ def compare_swrbudget(
     list2 = None
     if files2 is None:
         lst = get_entries_from_namefile(namefile2, "list")
-        list2 = lst[0][0]
+        list2 = lst[0][0] if any(lst) else None
     else:
         for file in files2:
             if (
@@ -388,7 +387,6 @@ def compare_swrbudget(
     # Open output file
     if outfile is not None:
         f = open(outfile, "w")
-        f.write("Created by pymake.autotest.compare\n")
 
     # Process cumulative and incremental
     for idx in range(2):
@@ -560,7 +558,7 @@ def compare_heads(
     if files1 is None:
         # Get oc info, and return if OC not included in models
         ocf1 = get_entries_from_namefile(namefile1, "OC")
-        if ocf1[0][0] is None:
+        if not any(ocf1) is None:
             return True
 
         hu1, hfpth1, du1, _ = flopy.modflow.ModflowOc.get_ocoutput_units(
@@ -572,7 +570,8 @@ def compare_heads(
             iut = du1
         if iut != 0:
             entries = get_entries_from_namefile(namefile1, unit=abs(iut))
-            hfpth1, status1 = entries[0][0], entries[0][1]
+            hfpth1 = entries[0][0] if any(entries) else None
+            status1 = entries[0][1] if any(entries) else None
 
     else:
         if isinstance(files1, str):
@@ -603,7 +602,7 @@ def compare_heads(
     if files2 is None:
         # Get oc info, and return if OC not included in models
         ocf2 = get_entries_from_namefile(namefile2, "OC")
-        if ocf2[0][0] is None:
+        if not any(ocf2):
             return True
 
         hu2, hfpth2, du2, dfpth2 = flopy.modflow.ModflowOc.get_ocoutput_units(
@@ -615,7 +614,8 @@ def compare_heads(
             iut = du2
         if iut != 0:
             entries = get_entries_from_namefile(namefile2, unit=abs(iut))
-            hfpth2, status2 = entries[0][0], entries[0][1]
+            hfpth2 = entries[0][0] if any(entries) else None
+            status2 = entries[0][1] if any(entries) else None
     else:
         if isinstance(files2, str):
             files2 = [files2]
@@ -656,7 +656,6 @@ def compare_heads(
     # Open output file
     if outfile is not None:
         f = open(outfile, "w")
-        f.write("Created by pymake.autotest.compare\n")
         f.write(f"Performing {text.upper()} to {text2.upper()} comparison\n")
 
         if exfile is not None:
@@ -953,7 +952,7 @@ def compare_concentrations(
     if files1 is None:
         for ext in valid_ext:
             ucn = get_entries_from_namefile(namefile1, extension=ext)
-            ufpth = ucn[0][0]
+            ufpth = ucn[0][0] if any(ucn) else None
             if ufpth is not None:
                 ufpth1 = ufpth
                 break
@@ -973,7 +972,7 @@ def compare_concentrations(
     if files2 is None:
         for ext in valid_ext:
             ucn = get_entries_from_namefile(namefile2, extension=ext)
-            ufpth = ucn[0][0]
+            ufpth = ucn[0][0] if any(ucn) else None
             if ufpth is not None:
                 ufpth2 = ufpth
                 break
@@ -1006,7 +1005,6 @@ def compare_concentrations(
     # Open output file
     if outfile is not None:
         f = open(outfile, "w")
-        f.write("Created by pymake.autotest.compare_concs\n")
 
     # Get stage objects
     uobj1 = flopy.utils.UcnFile(ufpth1, precision=precision, verbose=verbose)
@@ -1178,7 +1176,7 @@ def compare_stages(
     if namefile1 is not None:
         for ext in valid_ext:
             stg = get_entries_from_namefile(namefile1, extension=ext)
-            sfpth = stg[0][0]
+            sfpth = stg[0][0] if any(stg) else None
             if sfpth is not None:
                 sfpth1 = sfpth
                 break
@@ -1196,7 +1194,7 @@ def compare_stages(
     if namefile2 is not None:
         for ext in valid_ext:
             stg = get_entries_from_namefile(namefile2, extension=ext)
-            sfpth = stg[0][0]
+            sfpth = stg[0][0] if any(stg) else None
             if sfpth is not None:
                 sfpth2 = sfpth
                 break
@@ -1225,7 +1223,6 @@ def compare_stages(
     # Open output file
     if outfile is not None:
         f = open(outfile, "w")
-        f.write("Created by pymake.autotest.compare_stages\n")
 
     # Get stage objects
     sobj1 = flopy.utils.SwrStage(sfpth1, verbose=verbose)
diff --git a/flopy/utils/mfreadnam.py b/flopy/utils/mfreadnam.py
index c50a9497d..35fce16c9 100644
--- a/flopy/utils/mfreadnam.py
+++ b/flopy/utils/mfreadnam.py
@@ -292,9 +292,7 @@ def get_entries_from_namefile(
     -------
     entries : list of tuples
         list of tuples containing FTYPE, UNIT, FNAME, STATUS for each
-        namefile entry that meets a user-specified value. if no entries
-        are found, a single-element list containing a None-valued tuple is returned.
-
+        namefile entry that meets a user-specified value.
     """
     entries = []
     with open(path, "r") as f:
diff --git a/setup.cfg b/setup.cfg
index 5500d710c..744734b5a 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -59,7 +59,6 @@ test =
     filelock
     jupyter
     jupytext
-    mfpymake
     modflow-devtools
     pytest
     pytest-benchmark