diff --git a/autotest/test_compare.py b/autotest/test_compare.py
index 18e64e207b..f96865cc3b 100644
--- a/autotest/test_compare.py
+++ b/autotest/test_compare.py
@@ -1,8 +1,25 @@
+import os
+
 import numpy as np
 import pytest
+from modflow_devtools.markers import requires_exe, requires_pkg
 
 from flopy.mf6.utils import MfGrdFile
-from flopy.utils.compare import _diffmax, _difftol
+from flopy.modflow import (
+    Modflow,
+    ModflowBas,
+    ModflowDis,
+    ModflowLpf,
+    ModflowOc,
+    ModflowPcg,
+    ModflowWel,
+)
+from flopy.utils.compare import (
+    _diffmax,
+    _difftol,
+    compare_budget,
+    compare_heads,
+)
 
 
 def test_diffmax():
@@ -39,9 +56,108 @@ def test_eval_bud_diff(example_data_path):
     # TODO: create/run minimal model, then compare budget files
 
 
-@pytest.mark.skip(reason="todo")
-def test_compare_budget():
-    pass
+@pytest.fixture
+def comparison_model_1(function_tmpdir):
+    nlay = 3
+    nrow = 3
+    ncol = 3
+    model_name = "t1"
+
+    ml = Modflow(
+        modelname=model_name,
+        model_ws=function_tmpdir,
+        verbose=True,
+        exe_name="mf2005",
+    )
+    dis = ModflowDis(
+        ml, nlay=nlay, nrow=nrow, ncol=ncol, top=0, botm=[-1.0, -2.0, -3.0]
+    )
+    ibound = np.ones((nlay, nrow, ncol), dtype=int)
+    ibound[0, 1, 1] = 0
+    ibound[0, 0, -1] = -1
+    bas = ModflowBas(ml, ibound=ibound)
+    lpf = ModflowLpf(ml, ipakcb=102)
+    wd = ModflowWel.get_empty(ncells=2, aux_names=["v1", "v2"])
+    wd["k"][0] = 2
+    wd["i"][0] = 2
+    wd["j"][0] = 2
+    wd["flux"][0] = -1000.0
+    wd["v1"][0] = 1.0
+    wd["v2"][0] = 2.0
+    wd["k"][1] = 2
+    wd["i"][1] = 1
+    wd["j"][1] = 1
+    wd["flux"][1] = -500.0
+    wd["v1"][1] = 200.0
+    wd["v2"][1] = 100.0
+    wel_data = {0: wd}
+    wel = ModflowWel(ml, stress_period_data=wel_data, dtype=wd.dtype)
+    oc = ModflowOc(ml)
+    pcg = ModflowPcg(ml)
+
+    ml.write_input()
+
+    # run the modflow-2005 model
+    success, buff = ml.run_model(silent=False)
+    assert success, "could not run MODFLOW-2005 model"
+
+    # load the model
+    m = Modflow.load(
+        f"{model_name}.nam",
+        model_ws=function_tmpdir,
+        verbose=True,
+        exe_name="mf2005",
+    )
+
+    wl = m.wel.stress_period_data[0]
+    assert np.array_equal(wel.stress_period_data[0], wl), (
+        "previous well package stress period data does not match "
+        "stress period data loaded."
+    )
+
+    # change model work space
+    pth = os.path.join(function_tmpdir, "flopy")
+    m.change_model_ws(new_pth=pth)
+
+    # remove the existing well package
+    m.remove_package("WEL")
+
+    # recreate well package with binary output
+    wel = ModflowWel(
+        m, stress_period_data=wel_data, binary=True, dtype=wd.dtype
+    )
+
+    m.write_input()
+
+    fn1 = function_tmpdir / "flopy" / f"{model_name}.nam"
+    fn0 = function_tmpdir / f"{model_name}.nam"
+    fhsum = function_tmpdir / f"{os.path.splitext(model_name)[0]}.head.out"
+    fbsum = function_tmpdir / f"{os.path.splitext(model_name)[0]}.budget.out"
+
+    return m, fn1, fn0, fhsum, fbsum
+
+
+@requires_exe("mf2005")
+def test_compare_budget_and_heads(comparison_model_1):
+    m, fn1, fn0, fhsum, fbsum = comparison_model_1
+    success, buff = m.run_model()
+    assert success, "could not run the new MODFLOW-2005 model"
+
+    # compare the files
+    assert compare_heads(
+        fn0, fn1, outfile=fhsum
+    ), "head comparison failure (pathlib.Path)"
+    assert compare_heads(
+        str(fn0), str(fn1), outfile=fhsum
+    ), "head comparison failure (str path)"
+    assert compare_budget(
+        fn0, fn1, max_incpd=0.1, max_cumpd=0.1, outfile=fbsum
+    ), "budget comparison failure (pathlib.Path)"
+    assert compare_budget(
+        str(fn0), str(fn1), max_incpd=0.1, max_cumpd=0.1, outfile=str(fbsum)
+    ), "budget comparison failure (str path)"
+
+    # todo test with files1 and files2 arguments
 
 
 @pytest.mark.skip(reason="todo")
diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py
index f1782d6ec7..e2b4dc61e3 100644
--- a/flopy/export/shapefile_utils.py
+++ b/flopy/export/shapefile_utils.py
@@ -81,9 +81,9 @@ def write_grid_shapefile(
         value to fill nans
     epsg : str, int
         epsg code
-    prj : str or PathLike
+    prj : str or PathLike, optional, default None
         projection file path
-    verbose : bool
+    verbose : bool, default False
         whether to print verbose output
 
     Returns
@@ -234,7 +234,7 @@ def model_attributes_to_shapefile(
     array_dict : dict of {name:2D array} pairs
        Additional 2D arrays to add as attributes to the shapefile.
        (default is None)
-    verbose : bool
+    verbose : bool, optional, default False
         whether to print verbose output
     **kwargs : keyword arguments
         modelgrid : fp.modflow.Grid object
diff --git a/flopy/export/utils.py b/flopy/export/utils.py
index 545e03c617..27670be8fa 100644
--- a/flopy/export/utils.py
+++ b/flopy/export/utils.py
@@ -1,5 +1,5 @@
 import os
-from typing import Union
+from typing import Optional, Union
 
 import numpy as np
 
@@ -1516,6 +1516,8 @@ def export_array(
     fieldname : str
         Attribute field name for array values (shapefile export only).
         (default 'values')
+    verbose : bool, optional, default False
+        whether to show verbose output
     kwargs:
         keyword arguments to np.savetxt (ascii)
         rasterio.open (GeoTIFF)
@@ -1662,7 +1664,8 @@ def export_contours(
     contours,
     fieldname="level",
     epsg=None,
-    prj=None,
+    prj: Optional[Union[str, os.PathLike]] = None,
+    verbose=False,
     **kwargs,
 ):
     """
@@ -1678,8 +1681,10 @@ def export_contours(
         gis attribute table field name
     epsg : int
         EPSG code. See https://www.epsg-registry.org/ or spatialreference.org
-    prj : str
+    prj : str or PathLike, optional, default None
         Existing projection file to be used with new shapefile.
+    verbose : bool, optional, default False
+        whether to show verbose output
     **kwargs : key-word arguments to flopy.export.shapefile_utils.recarray2shp
 
     Returns
@@ -1705,7 +1710,9 @@ def export_contours(
     # convert the dictionary to a recarray
     ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray)
 
-    recarray2shp(ra, geoms, filename, epsg=epsg, prj=prj, **kwargs)
+    recarray2shp(
+        ra, geoms, filename, epsg=epsg, prj=prj, verbose=verbose, **kwargs
+    )
 
 
 def export_contourf(
@@ -1713,7 +1720,7 @@ def export_contourf(
     contours,
     fieldname="level",
     epsg=None,
-    prj=None,
+    prj: Optional[Union[str, os.PathLike]] = None,
     verbose=False,
     **kwargs,
 ):
@@ -1732,8 +1739,10 @@ def export_contourf(
         the range represented by the polygon.  Default is 'level'.
     epsg : int
         EPSG code. See https://www.epsg-registry.org/ or spatialreference.org
-    prj : str
+    prj : str or PathLike, optional, default None
         Existing projection file to be used with new shapefile.
+    verbose : bool, optional, default False
+        whether to show verbose output
 
     **kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
 
@@ -1800,7 +1809,9 @@ def export_contourf(
     # Create recarray
     ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray)
 
-    recarray2shp(ra, geoms, filename, epsg=epsg, prj=prj, **kwargs)
+    recarray2shp(
+        ra, geoms, filename, epsg=epsg, prj=prj, verbose=verbose, **kwargs
+    )
 
 
 def export_array_contours(
@@ -1812,7 +1823,8 @@ def export_array_contours(
     levels=None,
     maxlevels=1000,
     epsg=None,
-    prj=None,
+    prj: Optional[Union[str, os.PathLike]] = None,
+    verbose=False,
     **kwargs,
 ):
     """
@@ -1836,8 +1848,10 @@ def export_array_contours(
         maximum number of contour levels
     epsg : int
         EPSG code. See https://www.epsg-registry.org/ or spatialreference.org
-    prj : str
+    prj : str or PathLike, optional, default None
         Existing projection file to be used with new shapefile.
+    verbose : bool, optional, default False
+        whether to show verbose output
     **kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
 
     """
@@ -1857,7 +1871,9 @@ def export_array_contours(
         levels = np.arange(imin, imax, interval)
     ax = plt.subplots()[-1]
     ctr = contour_array(modelgrid, ax, a, levels=levels)
-    export_contours(filename, ctr, fieldname, epsg, prj, **kwargs)
+    export_contours(
+        filename, ctr, fieldname, epsg, prj, verbose=verbose, **kwargs
+    )
     plt.close()
 
 
@@ -1871,7 +1887,6 @@ def contour_array(modelgrid, ax, a, **kwargs):
         modelgrid object
     ax : matplotlib.axes.Axes
         ax to add the contours
-
     a : np.ndarray
         array to contour
 
diff --git a/flopy/mbase.py b/flopy/mbase.py
index b4f5c21261..5a3fc15808 100644
--- a/flopy/mbase.py
+++ b/flopy/mbase.py
@@ -689,7 +689,7 @@ def __getattr__(self, item):
 
     def get_ext_dict_attr(
         self,
-        ext_unit_dict: Union[str, os.PathLike] = None,
+        ext_unit_dict=None,
         unit=None,
         filetype=None,
         pop_key=True,
@@ -725,7 +725,7 @@ def _output_msg(self, i, add=True):
     def add_output_file(
         self,
         unit,
-        fname: Union[str, os.PathLike] = None,
+        fname: Optional[Union[str, os.PathLike]] = None,
         extension="cbc",
         binflag=True,
         package=None,
@@ -808,6 +808,7 @@ def add_output(
             binary or not. (default is False)
 
         """
+        fname = str(fname)
         if fname in self.output_fnames:
             if self.verbose:
                 print(
@@ -834,7 +835,7 @@ def add_output(
             self._output_msg(-1, add=True)
 
     def remove_output(
-        self, fname: Union[str, os.PathLike, None] = None, unit=None
+        self, fname: Optional[Union[str, os.PathLike]] = None, unit=None
     ):
         """
         Remove an output file from the model by specifying either the
@@ -848,6 +849,7 @@ def remove_output(
             Unit number of output array
         """
         if fname is not None:
+            fname = str(fname)
             for i, e in enumerate(self.output_fnames):
                 if fname in e:
                     if self.verbose:
@@ -883,6 +885,7 @@ def get_output(
         unit : int, optional
             Unit number of output array
         """
+        fname = str(fname)
         if fname is not None:
             for i, e in enumerate(self.output_fnames):
                 if fname in e:
@@ -917,6 +920,7 @@ def set_output_attribute(
         """
         idx = None
         if fname is not None:
+            fname = str(fname)
             for i, e in enumerate(self.output_fnames):
                 if fname in e:
                     idx = i
@@ -1717,6 +1721,7 @@ def run_model(
         normal_msg[idx] = s.lower()
 
     # Check to make sure that program and namefile exist
+    exe_name = str(Path(exe_name).expanduser().absolute())
     exe = which(exe_name)
     if exe is None:
         if exe_name.lower().endswith(".exe"):
diff --git a/flopy/mf6/mfbase.py b/flopy/mf6/mfbase.py
index 1084420759..7706a1203d 100644
--- a/flopy/mf6/mfbase.py
+++ b/flopy/mf6/mfbase.py
@@ -191,7 +191,7 @@ class MFFileMgmt:
     Parameters
     ----------
 
-    path : str
+    path : str or PathLike
         Path on disk to the simulation
 
     Attributes
diff --git a/flopy/mf6/modflow/mfsimulation.py b/flopy/mf6/modflow/mfsimulation.py
index c96e2e5933..9f519212fe 100644
--- a/flopy/mf6/modflow/mfsimulation.py
+++ b/flopy/mf6/modflow/mfsimulation.py
@@ -951,7 +951,7 @@ def check(
 
         Parameters
         ----------
-        f : str or PathLike
+        f : str or PathLike, optional
             String defining file name or file handle for summary file
             of check method output. If str or pathlike, a file handle
             is created. If None, the method does not write results to
@@ -2455,7 +2455,7 @@ def _is_in_solution_group(self, item, index, any_idx_after=False):
 
     def plot(
         self,
-        model_list: Union[str, List[str]] = None,
+        model_list: Optional[Union[str, List[str]]] = None,
         SelPackList=None,
         **kwargs,
     ):
@@ -2467,9 +2467,9 @@ def plot(
 
         Parameters
         ----------
-            model_list: (list)
+            model_list: list, optional
                 List of model names to plot, if none all models will be plotted
-            SelPackList: (list)
+            SelPackList: list, optional
                 List of package names to plot, if none all packages will be
                 plotted
             kwargs:
diff --git a/flopy/mfusg/mfusg.py b/flopy/mfusg/mfusg.py
index 5408b0b175..ba26e67338 100644
--- a/flopy/mfusg/mfusg.py
+++ b/flopy/mfusg/mfusg.py
@@ -15,7 +15,7 @@ class MfUsg(Modflow):
 
     Parameters
     ----------
-    modelname : str, default "modflowusgtest".
+    modelname : str or PathLike, default "modflowusgtest".
         Name of model.  This string will be used to name the MODFLOW input
         that are created with write_model.
     namefile_ext : str, default "nam"
@@ -159,7 +159,7 @@ def load(
 
         Parameters
         ----------
-        f : str
+        f : str or PathLike
             Path to MODFLOW name file to load.
         version : str, default "mfusg"
             MODFLOW version. Must be "mfusg".
diff --git a/flopy/utils/check.py b/flopy/utils/check.py
index 3db4c85889..5b401375b3 100644
--- a/flopy/utils/check.py
+++ b/flopy/utils/check.py
@@ -18,6 +18,9 @@ class check:
     ----------
     package : object
         Instance of Package class.
+    f : str or PathLike, optional
+        Path to the summary file. If no path is provided, a summary
+        file is not created and results are only written to stdout.
     verbose : bool
         Boolean flag used to determine if check method results are
         written to the screen
diff --git a/flopy/utils/compare.py b/flopy/utils/compare.py
index 457c10200c..cd50bc6516 100644
--- a/flopy/utils/compare.py
+++ b/flopy/utils/compare.py
@@ -95,9 +95,9 @@ def compare_budget(
 
     Parameters
     ----------
-    namefile1 : str
+    namefile1 : str or PathLike, optional
         namefile path for base model
-    namefile2 : str
+    namefile2 : str or PathLike, optional
         namefile path for comparison model
     max_cumpd : float
         maximum percent discrepancy allowed for cumulative budget terms
@@ -108,11 +108,11 @@ def compare_budget(
     outfile : str
         budget comparison output file name. If outfile is None, no
         comparison output is saved. (default is None)
-    files1 : str
+    files1 : str, PathLike, or list, optional
         base model output file. If files1 is not None, results
         will be extracted from files1 and namefile1 will not be used.
         (default is None)
-    files2 : str
+    files2 : str, PathLike, or list, optional
         comparison model output file. If files2 is not None, results
         will be extracted from files2 and namefile2 will not be used.
         (default is None)
@@ -140,7 +140,7 @@ def compare_budget(
         lst_file = get_entries_from_namefile(namefile1, "list")
         lst_file1 = lst_file[0][0] if any(lst_file) else None
     else:
-        if isinstance(files1, str):
+        if isinstance(files1, (str, os.PathLike)):
             files1 = [files1]
         for file in files1:
             if (
@@ -154,7 +154,7 @@ def compare_budget(
         lst_file = get_entries_from_namefile(namefile2, "list")
         lst_file2 = lst_file[0][0] if any(lst_file) else None
     else:
-        if isinstance(files2, str):
+        if isinstance(files2, (str, os.PathLike)):
             files2 = [files2]
         for file in files2:
             if (
@@ -302,9 +302,9 @@ def compare_swrbudget(
 
     Parameters
     ----------
-    namefile1 : str
+    namefile1 : str or PathLike, optional
         namefile path for base model
-    namefile2 : str
+    namefile2 : str or PathLike, optional
         namefile path for comparison model
     max_cumpd : float
         maximum percent discrepancy allowed for cumulative budget terms
@@ -312,14 +312,14 @@ def compare_swrbudget(
     max_incpd : float
         maximum percent discrepancy allowed for incremental budget terms
         (default is 0.01)
-    outfile : str
+    outfile : str or PathLike, optional
         budget comparison output file name. If outfile is None, no
         comparison output is saved. (default is None)
-    files1 : str
+    files1 : str, PathLike, or list, optional
         base model output file. If files1 is not None, results
         will be extracted from files1 and namefile1 will not be used.
         (default is None)
-    files2 : str
+    files2 : str, PathLike, or list, optional
         comparison model output file. If files2 is not None, results
         will be extracted from files2 and namefile2 will not be used.
         (default is None)
@@ -534,7 +534,7 @@ def compare_heads(
     verbose : bool
         boolean indicating if verbose output should be written to the
         terminal (default is False)
-    exfile : str or PathLike
+    exfile : str or PathLike, optional
         path to a file with exclusion array data. Head differences will not
         be evaluated where exclusion array values are greater than zero.
         (default is None)
@@ -578,7 +578,7 @@ def compare_heads(
             status1 = entries[0][1] if any(entries) else None
 
     else:
-        if isinstance(files1, str):
+        if isinstance(files1, (str, os.PathLike)):
             files1 = [files1]
         for file in files1:
             if text.lower() == "head":
@@ -619,7 +619,7 @@ def compare_heads(
             hfpth2 = entries[0][0] if any(entries) else None
             status2 = entries[0][1] if any(entries) else None
     else:
-        if isinstance(files2, str):
+        if isinstance(files2, (str, os.PathLike)):
             files2 = [files2]
         for file in files2:
             if text2.lower() == "head":
@@ -683,7 +683,7 @@ def compare_heads(
     # get data from exclusion file
     if exfile is not None:
         e = None
-        if isinstance(exfile, str):
+        if isinstance(exfile, (str, os.PathLike)):
             try:
                 exd = np.genfromtxt(exfile).flatten()
             except:
@@ -899,23 +899,23 @@ def compare_concentrations(
 
     Parameters
     ----------
-    namefile1 : str
+    namefile1 : str or PathLike
         namefile path for base model
-    namefile2 : str
+    namefile2 : str or PathLike
         namefile path for comparison model
     precision : str
         precision for binary head file ("auto", "single", or "double")
         default is "auto"
     ctol : float
         maximum allowed concentration difference (default is 0.001)
-    outfile : str
+    outfile : str or PathLike, optional
         concentration comparison output file name. If outfile is None, no
         comparison output is saved. (default is None)
-    files1 : str
+    files1 : str, PathLike, or list, optional
         base model output file. If files1 is not None, results
         will be extracted from files1 and namefile1 will not be used.
         (default is None)
-    files2 : str
+    files2 : str, PathLike, or list, optional
         comparison model output file. If files2 is not None, results
         will be extracted from files2 and namefile2 will not be used.
         (default is None)
@@ -1131,23 +1131,23 @@ def compare_stages(
 
     Parameters
     ----------
-    namefile1 : str
+    namefile1 : str or PathLike
         namefile path for base model
-    namefile2 : str
+    namefile2 : str or PathLike
         namefile path for comparison model
     precision : str
         precision for binary head file ("auto", "single", or "double")
         default is "auto"
     htol : float
         maximum allowed stage difference (default is 0.001)
-    outfile : str
+    outfile : str or PathLike, optional
         head comparison output file name. If outfile is None, no
         comparison output is saved. (default is None)
-    files1 : str
+    files1 : str, PathLike, or list, optional
         base model output file. If files1 is not None, results
         will be extracted from files1 and namefile1 will not be used.
         (default is None)
-    files2 : str
+    files2 : str, PathLike, or list, optional
         comparison model output file. If files2 is not None, results
         will be extracted from files2 and namefile2 will not be used.
         (default is None)
@@ -1183,7 +1183,7 @@ def compare_stages(
                 sfpth1 = sfpth
                 break
     elif files1 is not None:
-        if isinstance(files1, str):
+        if isinstance(files1, (str, os.PathLike)):
             files1 = [files1]
         for file in files1:
             for ext in valid_ext:
@@ -1201,7 +1201,7 @@ def compare_stages(
                 sfpth2 = sfpth
                 break
     elif files2 is not None:
-        if isinstance(files2, str):
+        if isinstance(files2, (str, os.PathLike)):
             files2 = [files2]
         for file in files2:
             for ext in valid_ext:
@@ -1350,9 +1350,9 @@ def compare(
 
     Parameters
     ----------
-    namefile1 : str
+    namefile1 : str or PathLike, optional
         namefile path for base model
-    namefile2 : str
+    namefile2 : str or PathLike, optional
         namefile path for comparison model
     precision : str
         precision for binary head file ("auto", "single", or "double")
@@ -1365,17 +1365,17 @@ def compare(
         (default is 0.01)
     htol : float
         maximum allowed head difference (default is 0.001)
-    outfile1 : str
+    outfile1 : str or PathLike, optional
         budget comparison output file name. If outfile1 is None, no budget
         comparison output is saved. (default is None)
-    outfile2 : str
+    outfile2 : str or PathLike, optional
         head comparison output file name. If outfile2 is None, no head
         comparison output is saved. (default is None)
-    files1 : str
+    files1 : str, PathLike, or list, optional
         base model output file. If files1 is not None, results
         will be extracted from files1 and namefile1 will not be used.
         (default is None)
-    files2 : str
+    files2 : str, PathLike, or list, optional
         comparison model output file. If files2 is not None, results
         will be extracted from files2 and namefile2 will not be used.
         (default is None)