diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
index 2dea681dfc..9c22c4e238 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 cbf1e3e5e9..781a3607a9 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 076ba8223b..b6772337d8 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 e1ec0a0abc..9934964278 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 6faebce7ee..8885eb3ca4 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 45d20fcad5..2373a31493 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 e7b65fbd43..fb7988b024 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 664ffb7046..9a14a365f0 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 7b2ac4ee25..d2b48bf53b 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 4daf494d75..2dc957a082 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 0ad095c28e..3758cf07d1 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 8a37df8af6..045cfc0603 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 0af16e34d9..4a9af3dad8 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 4cdf70fefa..27a3025b8e 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 c64b98e6d3..da0750a410 100644
--- a/flopy/utils/compare.py
+++ b/flopy/utils/compare.py
@@ -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)
@@ -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):
@@ -656,7 +654,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:
@@ -1006,7 +1003,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)
@@ -1225,7 +1221,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/crosssection.py b/flopy/utils/crosssection.py
new file mode 100644
index 0000000000..889363d3eb
--- /dev/null
+++ b/flopy/utils/crosssection.py
@@ -0,0 +1,280 @@
+import numpy as np
+
+# power for Manning's hydraulic radius term
+_mpow = 2.0 / 3.0
+
+
+def calculate_rectchan_mannings_discharge(
+    conversion_factor, roughness, slope, width, depth
+):
+    """
+    Calculate Manning's discharge for a rectangular channel.
+
+    """
+    area = width * depth
+    return conversion_factor * area * depth**_mpow * slope**0.5 / roughness
+
+
+# n-point cross-section functions
+def get_wetted_station(
+    x0,
+    x1,
+    h0,
+    h1,
+    depth,
+):
+    """Get the wetted length in the x-direction"""
+    # -- calculate the minimum and maximum depth
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    # -- if depth is less than or equal to the minimum value the
+    #    station length (xlen) is zero
+    if depth <= hmin:
+        x1 = x0
+    # -- if depth is between hmin and hmax, station length is less
+    #    than h1 - h0
+    elif depth < hmax:
+        xlen = x1 - x0
+        dlen = h1 - h0
+        if abs(dlen) > 0.0:
+            slope = xlen / dlen
+        else:
+            slope = 0.0
+        if h0 > h1:
+            dx = (depth - h1) * slope
+            xt = x1 + dx
+            xt0 = xt
+            xt1 = x1
+        else:
+            dx = (depth - h0) * slope
+            xt = x0 + dx
+            xt0 = x0
+            xt1 = xt
+        x0 = xt0
+        x1 = xt1
+    return x0, x1
+
+
+def get_wetted_perimeter(
+    x0,
+    x1,
+    h0,
+    h1,
+    depth,
+):
+    # -- calculate the minimum and maximum depth
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    # -- calculate the wetted perimeter for the segment
+    xlen = x1 - x0
+    if xlen > 0.0:
+        if depth > hmax:
+            dlen = hmax - hmin
+        else:
+            dlen = depth - hmin
+    else:
+        if depth > hmin:
+            dlen = min(depth, hmax) - hmin
+        else:
+            dlen = 0.0
+    return np.sqrt(xlen**2.0 + dlen**2.0)
+
+
+def get_wetted_area(x0, x1, h0, h1, depth):
+    # -- calculate the minimum and maximum depth
+    hmin = min(h0, h1)
+    hmax = max(h0, h1)
+
+    # -- calculate the wetted area for the segment
+    xlen = x1 - x0
+    area = 0.0
+    if xlen > 0.0:
+        # -- add the area above hmax
+        if depth > hmax:
+            area = xlen * (depth - hmax)
+        # -- add the area below zmax
+        if hmax != hmin and depth > hmin:
+            area += 0.5 * (depth - hmin)
+    return area
+
+
+def wetted_area(
+    x,
+    h,
+    depth,
+    verbose=False,
+):
+    area = 0.0
+    if x.shape[0] == 1:
+        area = x[0] * depth
+    else:
+        for idx in range(0, x.shape[0] - 1):
+            x0, x1 = x[idx], x[idx + 1]
+            h0, h1 = h[idx], h[idx + 1]
+
+            # get station data
+            x0, x1 = get_wetted_station(x0, x1, h0, h1, depth)
+
+            # get wetted area
+            a = get_wetted_area(x0, x1, h0, h1, depth)
+            area += a
+
+            # write to screen
+            if verbose:
+                print(
+                    f"{idx}->{idx + 1} ({x0},{x1}) - "
+                    f"perimeter={x1 - x0} - area={a}"
+                )
+
+    return area
+
+
+def wetted_perimeter(
+    x,
+    h,
+    depth,
+    verbose=False,
+):
+    perimeter = 0.0
+    if x.shape[0] == 1:
+        perimeter = x[0]
+    else:
+        for idx in range(0, x.shape[0] - 1):
+            x0, x1 = x[idx], x[idx + 1]
+            h0, h1 = h[idx], h[idx + 1]
+
+            # get station data
+            x0, x1 = get_wetted_station(x0, x1, h0, h1, depth)
+
+            # get wetted perimeter
+            perimeter += get_wetted_perimeter(x0, x1, h0, h1, depth)
+
+            # write to screen
+            if verbose:
+                print(f"{idx}->{idx + 1} ({x0},{x1}) - perimeter={x1 - x0}")
+
+    return perimeter
+
+
+def manningsq(
+    x,
+    h,
+    depth,
+    roughness=0.01,
+    slope=0.001,
+    conv=1.0,
+):
+    if isinstance(roughness, float):
+        roughness = np.ones(x.shape, dtype=float) * roughness
+    if x.shape[0] > 1:
+        q = 0.0
+        for i0 in range(x.shape[0] - 1):
+            i1 = i0 + 1
+            perimeter = get_wetted_perimeter(x[i0], x[i1], h[i0], h[i1], depth)
+            area = get_wetted_area(x[i0], x[i1], h[i0], h[i1], depth)
+            if perimeter > 0.0:
+                radius = area / perimeter
+                q += (
+                    conv
+                    * area
+                    * radius**_mpow
+                    * slope**0.5
+                    / roughness[i0]
+                )
+    else:
+        perimeter = wetted_perimeter(x, h, depth)
+        area = wetted_area(x, h, depth)
+        radius = 0.0
+        if perimeter > 0.0:
+            radius = area / perimeter
+        q = conv * area * radius**_mpow * slope**0.5 / roughness[0]
+    return q
+
+
+def get_depths(
+    flows,
+    x,
+    h,
+    roughness=0.01,
+    slope=0.001,
+    conv=1.0,
+    dd=1e-4,
+    verbose=False,
+):
+    if isinstance(flows, float):
+        flows = np.array([flows], dtype=float)
+    if isinstance(roughness, float):
+        roughness = np.ones(x.shape, dtype=float) * roughness
+    depths = np.zeros(flows.shape, dtype=float)
+    for idx, q in enumerate(flows):
+        depths[idx] = qtodepth(
+            x,
+            h,
+            q,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+            dd=dd,
+            verbose=False,
+        )
+
+    return depths
+
+
+def qtodepth(
+    x,
+    h,
+    q,
+    roughness=0.01,
+    slope=0.001,
+    conv=1.0,
+    dd=1e-4,
+    verbose=False,
+):
+    h0 = 0.0
+    q0 = manningsq(
+        x,
+        h,
+        h0,
+        roughness=roughness,
+        slope=slope,
+        conv=conv,
+    )
+    r = q0 - q
+
+    iter = 0
+    if verbose:
+        print(f"iteration {iter:>2d} - residual={r}")
+    while abs(r) > 1e-12:
+        q1 = manningsq(
+            x,
+            h,
+            h0 + dd,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+        )
+        dq = q1 - q0
+        if dq != 0.0:
+            derv = dd / (q1 - q0)
+        else:
+            derv = 0.0
+        h0 -= derv * r
+        q0 = manningsq(
+            x,
+            h,
+            h0,
+            roughness=roughness,
+            slope=slope,
+            conv=conv,
+        )
+        r = q0 - q
+
+        iter += 1
+        if verbose:
+            print(f"iteration {iter:>2d} - residual={r}")
+        if iter > 100:
+            break
+    return h0
diff --git a/setup.cfg b/setup.cfg
index 03d22cb432..9e4701f232 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -59,7 +59,6 @@ test =
     filelock
     jupyter
     jupytext
-    mfpymake
     modflow-devtools
     pytest
     pytest-benchmark