diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py
index d7254151c8..6c56652243 100644
--- a/flopy/modflow/mf.py
+++ b/flopy/modflow/mf.py
@@ -747,7 +747,7 @@ def load(
             os.path.join(model_ws, f)
         )
 
-        ml = Modflow(
+        ml = cls(
             modelname,
             version=version,
             exe_name=exe_name,
diff --git a/flopy/modflow/mfag.py b/flopy/modflow/mfag.py
index ffe8e9ce17..00d812e788 100644
--- a/flopy/modflow/mfag.py
+++ b/flopy/modflow/mfag.py
@@ -891,7 +891,7 @@ def load(cls, f, model, nper=0, ext_unit_dict=None):
                         "Something went wrong at: {}".format(line)
                     )
 
-        return ModflowAg(
+        return cls(
             model,
             options=options,
             time_series=time_series,
diff --git a/flopy/modflow/mfbas.py b/flopy/modflow/mfbas.py
index d82e5a1d4e..30e502c722 100644
--- a/flopy/modflow/mfbas.py
+++ b/flopy/modflow/mfbas.py
@@ -398,7 +398,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True, **kwargs):
             )
 
         # create bas object and return
-        bas = ModflowBas(
+        bas = cls(
             model,
             ibound=ibound,
             strt=strt,
diff --git a/flopy/modflow/mfbcf.py b/flopy/modflow/mfbcf.py
index 0bc6d6ff99..baaa2c8cce 100644
--- a/flopy/modflow/mfbcf.py
+++ b/flopy/modflow/mfbcf.py
@@ -529,7 +529,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 model.add_pop_key_list(ipakcb)
 
         # create instance of bcf object
-        bcf = ModflowBcf(
+        bcf = cls(
             model,
             ipakcb=ipakcb,
             intercellt=intercellt,
diff --git a/flopy/modflow/mfchd.py b/flopy/modflow/mfchd.py
index 23e7e7e31c..5c8ff28e7f 100644
--- a/flopy/modflow/mfchd.py
+++ b/flopy/modflow/mfchd.py
@@ -267,7 +267,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         return Package.load(
             f,
             model,
-            ModflowChd,
+            cls,
             nper=nper,
             check=check,
             ext_unit_dict=ext_unit_dict,
diff --git a/flopy/modflow/mfde4.py b/flopy/modflow/mfde4.py
index 53658abf2c..09be4c0e1f 100644
--- a/flopy/modflow/mfde4.py
+++ b/flopy/modflow/mfde4.py
@@ -313,7 +313,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowDe4.ftype()
             )
 
-        de4 = ModflowDe4(
+        de4 = cls(
             model,
             itmx=itmx,
             mxup=mxup,
diff --git a/flopy/modflow/mfdis.py b/flopy/modflow/mfdis.py
index 231d1596d3..d9f9dc6ec3 100644
--- a/flopy/modflow/mfdis.py
+++ b/flopy/modflow/mfdis.py
@@ -809,45 +809,6 @@ def check(self, f=None, verbose=True, level=1, checktype=None):
         chk.summarize()
         return chk
 
-        # if f is not None:
-        #     if isinstance(f, str):
-        #         pth = os.path.join(self.parent.model_ws, f)
-        #         f = open(pth, 'w', 0)
-        #
-        # errors = False
-        # txt = '\n{} PACKAGE DATA VALIDATION:\n'.format(self.name[0])
-        # t = ''
-        # t1 = ''
-        # inactive = self.parent.bas6.ibound.array == 0
-        # # thickness errors
-        # d = self.thickness.array
-        # d[inactive] = 1.
-        # if d.min() <= 0:
-        #     errors = True
-        #     t = '{}  ERROR: Negative or zero cell thickness specified.\n'.format(t)
-        #     if level > 0:
-        #         idx = np.column_stack(np.where(d <= 0.))
-        #         t1 = self.level1_arraylist(idx, d, self.thickness.name, t1)
-        # else:
-        #     t = '{}  Specified cell thickness is OK.\n'.format(t)
-        #
-        # # add header to level 0 text
-        # txt += t
-        #
-        # if level > 0:
-        #     if errors:
-        #         txt += '\n  DETAILED SUMMARY OF {} ERRORS:\n'.format(self.name[0])
-        #         # add level 1 header to level 1 text
-        #         txt += t1
-        #
-        # # write errors to summary file
-        # if f is not None:
-        #     f.write('{}\n'.format(txt))
-        #
-        # # write errors to stdout
-        # if verbose:
-        #     print(txt)
-
     @classmethod
     def load(cls, f, model, ext_unit_dict=None, check=True):
         """
@@ -1067,7 +1028,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
             )
 
         # create dis object instance
-        dis = ModflowDis(
+        dis = cls(
             model,
             nlay=nlay,
             nrow=nrow,
diff --git a/flopy/modflow/mfdisu.py b/flopy/modflow/mfdisu.py
index 6e9bc4e31d..329767d0b5 100644
--- a/flopy/modflow/mfdisu.py
+++ b/flopy/modflow/mfdisu.py
@@ -798,7 +798,7 @@ def load(cls, f, model, ext_unit_dict=None, check=False):
             )
 
         # create dis object instance
-        disu = ModflowDisU(
+        disu = cls(
             model,
             nodes=nodes,
             nlay=nlay,
diff --git a/flopy/modflow/mfdrn.py b/flopy/modflow/mfdrn.py
index a3481e8d2d..aa3c527d03 100644
--- a/flopy/modflow/mfdrn.py
+++ b/flopy/modflow/mfdrn.py
@@ -328,7 +328,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         return Package.load(
             f,
             model,
-            ModflowDrn,
+            cls,
             nper=nper,
             check=check,
             ext_unit_dict=ext_unit_dict,
diff --git a/flopy/modflow/mfdrt.py b/flopy/modflow/mfdrt.py
index 97e2a75aee..af84946148 100644
--- a/flopy/modflow/mfdrt.py
+++ b/flopy/modflow/mfdrt.py
@@ -314,7 +314,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         return Package.load(
             f,
             model,
-            ModflowDrt,
+            cls,
             nper=nper,
             check=check,
             ext_unit_dict=ext_unit_dict,
diff --git a/flopy/modflow/mfevt.py b/flopy/modflow/mfevt.py
index 1e8f273cb0..c7842c5284 100644
--- a/flopy/modflow/mfevt.py
+++ b/flopy/modflow/mfevt.py
@@ -417,7 +417,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
         args["unitnumber"] = unitnumber
         args["filenames"] = filenames
 
-        evt = ModflowEvt(model, **args)
+        evt = cls(model, **args)
 
         # return evt object
         return evt
diff --git a/flopy/modflow/mffhb.py b/flopy/modflow/mffhb.py
index 5e77062d3e..3a7621e5ab 100644
--- a/flopy/modflow/mffhb.py
+++ b/flopy/modflow/mffhb.py
@@ -756,7 +756,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
         nfhbx1 = 0
         nfhbx2 = 0
 
-        fhb = ModflowFhb(
+        fhb = cls(
             model,
             nbdtim=nbdtim,
             nflw=nflw,
diff --git a/flopy/modflow/mfflwob.py b/flopy/modflow/mfflwob.py
index 96cc3945f9..f4642b9159 100755
--- a/flopy/modflow/mfflwob.py
+++ b/flopy/modflow/mfflwob.py
@@ -537,7 +537,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                 model.add_pop_key_list(iufbobsv)
 
         # create ModflowFlwob object instance
-        flwob = ModflowFlwob(
+        flwob = cls(
             model,
             iufbobsv=iufbobsv,
             tomultfb=tomultfb,
diff --git a/flopy/modflow/mfgage.py b/flopy/modflow/mfgage.py
index 826d0e11af..b5fb14ea9c 100644
--- a/flopy/modflow/mfgage.py
+++ b/flopy/modflow/mfgage.py
@@ -395,14 +395,13 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
         for file in files:
             filenames.append(os.path.basename(file))
 
-        gagepak = ModflowGage(
+        return cls(
             model,
             numgage=numgage,
             gage_data=gage_data,
             filenames=filenames,
             unitnumber=unitnumber,
         )
-        return gagepak
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfghb.py b/flopy/modflow/mfghb.py
index 65bc83ca69..c4b3ad535e 100644
--- a/flopy/modflow/mfghb.py
+++ b/flopy/modflow/mfghb.py
@@ -307,7 +307,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         return Package.load(
             f,
             model,
-            ModflowGhb,
+            cls,
             nper=nper,
             check=check,
             ext_unit_dict=ext_unit_dict,
diff --git a/flopy/modflow/mfgmg.py b/flopy/modflow/mfgmg.py
index d02e21243d..065ec8d8f9 100644
--- a/flopy/modflow/mfgmg.py
+++ b/flopy/modflow/mfgmg.py
@@ -410,8 +410,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 )
                 model.add_pop_key_list(iunitmhc)
 
-        # create the gmg object
-        gmg = ModflowGmg(
+        return cls(
             model,
             mxiter=mxiter,
             iiter=iiter,
@@ -429,7 +428,6 @@ def load(cls, f, model, ext_unit_dict=None):
             chglimit=chglimit,
             unitnumber=unitnumber,
         )
-        return gmg
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py
index 35d6358dbe..586a4e5b6b 100644
--- a/flopy/modflow/mfhfb.py
+++ b/flopy/modflow/mfhfb.py
@@ -412,7 +412,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowHfb.ftype()
             )
 
-        hfb = ModflowHfb(
+        return cls(
             model,
             nphfb=0,
             mxfb=0,
@@ -423,7 +423,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return hfb
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfhob.py b/flopy/modflow/mfhob.py
index e2518e9ed9..a035f26d2a 100755
--- a/flopy/modflow/mfhob.py
+++ b/flopy/modflow/mfhob.py
@@ -491,8 +491,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                     )
                     model.add_pop_key_list(iuhobsv)
 
-        # create hob object instance
-        hob = ModflowHob(
+        return cls(
             model,
             iuhobsv=iuhobsv,
             hobdry=hobdry,
@@ -502,8 +501,6 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
             filenames=filenames,
         )
 
-        return hob
-
     @staticmethod
     def ftype():
         return "HOB"
diff --git a/flopy/modflow/mfhyd.py b/flopy/modflow/mfhyd.py
index 34d928724b..784cabef5d 100644
--- a/flopy/modflow/mfhyd.py
+++ b/flopy/modflow/mfhyd.py
@@ -366,8 +366,8 @@ def load(cls, f, model, ext_unit_dict=None):
                 )
                 model.add_pop_key_list(ihydun)
 
-        # create hyd instance
-        hyd = ModflowHyd(
+        # return hyd instance
+        return cls(
             model,
             nhyd=nhyd,
             ihydun=ihydun,
@@ -377,9 +377,6 @@ def load(cls, f, model, ext_unit_dict=None):
             filenames=filenames,
         )
 
-        # return hyd instance
-        return hyd
-
     @staticmethod
     def ftype():
         return "HYD"
diff --git a/flopy/modflow/mflak.py b/flopy/modflow/mflak.py
index 8f8e02a3bb..d6ae3dcd15 100644
--- a/flopy/modflow/mflak.py
+++ b/flopy/modflow/mflak.py
@@ -908,7 +908,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
                     )
                     ipos += 1
 
-        lakpak = ModflowLak(
+        return cls(
             model,
             options=options,
             nlakes=nlakes,
@@ -928,7 +928,6 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return lakpak
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mflmt.py b/flopy/modflow/mflmt.py
index 4096b6506b..e3673cbb8b 100644
--- a/flopy/modflow/mflmt.py
+++ b/flopy/modflow/mflmt.py
@@ -257,7 +257,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowLmt.ftype()
             )
 
-        lmt = ModflowLmt(
+        return cls(
             model,
             output_file_name=output_file_name,
             output_file_unit=output_file_unit,
@@ -267,7 +267,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return lmt
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mflpf.py b/flopy/modflow/mflpf.py
index e24c8184c7..e2dfdb4314 100644
--- a/flopy/modflow/mflpf.py
+++ b/flopy/modflow/mflpf.py
@@ -723,7 +723,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                 model.add_pop_key_list(ipakcb)
 
         # create instance of lpf class
-        lpf = ModflowLpf(
+        lpf = cls(
             model,
             ipakcb=ipakcb,
             laytyp=laytyp,
diff --git a/flopy/modflow/mfmlt.py b/flopy/modflow/mfmlt.py
index ea1603f6db..ac0bf6a2bd 100644
--- a/flopy/modflow/mfmlt.py
+++ b/flopy/modflow/mfmlt.py
@@ -232,16 +232,13 @@ def load(cls, f, model, nrow=None, ncol=None, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowMlt.ftype()
             )
 
-        # create mlt dictionary
-        mlt = ModflowMlt(
+        return cls(
             model,
             mult_dict=mult_dict,
             unitnumber=unitnumber,
             filenames=filenames,
         )
 
-        return mlt
-
     @staticmethod
     def mult_function(mult_dict, line):
         """
diff --git a/flopy/modflow/mfmnw1.py b/flopy/modflow/mfmnw1.py
index 32a1e7c985..75b8ffabc7 100644
--- a/flopy/modflow/mfmnw1.py
+++ b/flopy/modflow/mfmnw1.py
@@ -278,7 +278,7 @@ def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
         if openfile:
             f.close()
 
-        return ModflowMnw1(
+        return cls(
             model,
             mxmnw=mxmnw,
             ipakcb=ipakcb,
diff --git a/flopy/modflow/mfmnw2.py b/flopy/modflow/mfmnw2.py
index 08eea66b55..d7f1603747 100644
--- a/flopy/modflow/mfmnw2.py
+++ b/flopy/modflow/mfmnw2.py
@@ -1464,7 +1464,7 @@ def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
                         filenames[1] = os.path.basename(value.filename)
                         model.add_pop_key_list(key)
 
-        return ModflowMnw2(
+        return cls(
             model,
             mnwmax=mnwmax,
             nodtot=nodtot,
diff --git a/flopy/modflow/mfmnwi.py b/flopy/modflow/mfmnwi.py
index 6b77fbb430..54e69090f5 100644
--- a/flopy/modflow/mfmnwi.py
+++ b/flopy/modflow/mfmnwi.py
@@ -281,7 +281,7 @@ def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
                 )
                 idx += 1
 
-        return ModflowMnwi(
+        return cls(
             model,
             wel1flag=wel1flag,
             qsumflag=qsumflag,
diff --git a/flopy/modflow/mfnwt.py b/flopy/modflow/mfnwt.py
index fbb1543b08..ee7a4e2ab2 100644
--- a/flopy/modflow/mfnwt.py
+++ b/flopy/modflow/mfnwt.py
@@ -556,7 +556,7 @@ def load(cls, f, model, ext_unit_dict=None):
         kwargs["filenames"] = filenames
 
         # create and return an instance of the nwt class
-        return ModflowNwt(model, **kwargs)
+        return cls(model, **kwargs)
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfoc.py b/flopy/modflow/mfoc.py
index adfea5467a..3d2bf63f47 100644
--- a/flopy/modflow/mfoc.py
+++ b/flopy/modflow/mfoc.py
@@ -1104,8 +1104,7 @@ def load(
         for u in unitnumber:
             model.add_pop_key_list(u)
 
-        # create instance of oc class
-        oc = ModflowOc(
+        return cls(
             model,
             ihedfm=ihedfm,
             iddnfm=iddnfm,
@@ -1118,8 +1117,6 @@ def load(
             filenames=filenames,
         )
 
-        return oc
-
     @staticmethod
     def ftype():
         return "OC"
diff --git a/flopy/modflow/mfpar.py b/flopy/modflow/mfpar.py
index de0b7abaea..cc5a1ce58d 100644
--- a/flopy/modflow/mfpar.py
+++ b/flopy/modflow/mfpar.py
@@ -193,8 +193,8 @@ def set_pval(self, model, ext_unit_dict):
 
         return
 
-    @classmethod
-    def load(cls, f, npar, verbose=False):
+    @staticmethod
+    def load(f, npar, verbose=False):
         """
         Load property parameters from an existing package.
 
diff --git a/flopy/modflow/mfparbc.py b/flopy/modflow/mfparbc.py
index 2793fef2e8..a12d9dd87a 100644
--- a/flopy/modflow/mfparbc.py
+++ b/flopy/modflow/mfparbc.py
@@ -118,9 +118,7 @@ def load(cls, f, npar, dt, model, ext_unit_dict=None, verbose=False):
                     pinst,
                 ]
 
-        # print bc_parms
-        bcpar = ModflowParBc(bc_parms)
-        return bcpar
+        return cls(bc_parms)
 
     @staticmethod
     def loadarray(f, npar, verbose=False):
diff --git a/flopy/modflow/mfpcg.py b/flopy/modflow/mfpcg.py
index 37ae54cad8..9b19e90a08 100644
--- a/flopy/modflow/mfpcg.py
+++ b/flopy/modflow/mfpcg.py
@@ -342,8 +342,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowPcg.ftype()
             )
 
-        # create instance of pcg class
-        pcg = ModflowPcg(
+        return cls(
             model,
             mxiter=mxiter,
             iter1=iter1,
@@ -360,7 +359,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return pcg
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfpcgn.py b/flopy/modflow/mfpcgn.py
index 0a19b484cc..26259ec594 100644
--- a/flopy/modflow/mfpcgn.py
+++ b/flopy/modflow/mfpcgn.py
@@ -558,7 +558,7 @@ def load(cls, f, model, ext_unit_dict=None):
                     ext_unit_dict, unit=ipunit
                 )
 
-        pcgn = ModflowPcgn(
+        return cls(
             model,
             iter_mo=iter_mo,
             iter_mi=iter_mi,
@@ -581,7 +581,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return pcgn
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfpks.py b/flopy/modflow/mfpks.py
index b1fd29e444..c49d31e221 100644
--- a/flopy/modflow/mfpks.py
+++ b/flopy/modflow/mfpks.py
@@ -287,8 +287,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowPks.ftype()
             )
 
-        pks = ModflowPks(model, unitnumber=unitnumber, filenames=filenames)
-        return pks
+        return cls(model, unitnumber=unitnumber, filenames=filenames)
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfpval.py b/flopy/modflow/mfpval.py
index 514c8555af..4493321007 100644
--- a/flopy/modflow/mfpval.py
+++ b/flopy/modflow/mfpval.py
@@ -218,13 +218,12 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowPval.ftype()
             )
 
-        pval = ModflowPval(
+        return cls(
             model,
             pval_dict=pval_dict,
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return pval
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfrch.py b/flopy/modflow/mfrch.py
index 022aaf64a8..49b7ce567e 100644
--- a/flopy/modflow/mfrch.py
+++ b/flopy/modflow/mfrch.py
@@ -527,7 +527,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
                 model.add_pop_key_list(ipakcb)
 
         # create recharge package instance
-        rch = ModflowRch(
+        rch = cls(
             model,
             nrchop=nrchop,
             ipakcb=ipakcb,
diff --git a/flopy/modflow/mfriv.py b/flopy/modflow/mfriv.py
index a238f94d23..e0ed25804d 100644
--- a/flopy/modflow/mfriv.py
+++ b/flopy/modflow/mfriv.py
@@ -389,7 +389,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         return Package.load(
             f,
             model,
-            ModflowRiv,
+            cls,
             nper=nper,
             check=check,
             ext_unit_dict=ext_unit_dict,
diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py
index 02aa59736a..23ededc759 100644
--- a/flopy/modflow/mfsfr2.py
+++ b/flopy/modflow/mfsfr2.py
@@ -1070,7 +1070,7 @@ def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
                         filenames[2] = os.path.basename(value.filename)
                         model.add_pop_key_list(key)
 
-        return ModflowSfr2(
+        return cls(
             model,
             nstrm=nstrm,
             nss=nss,
diff --git a/flopy/modflow/mfsip.py b/flopy/modflow/mfsip.py
index 6d8a4dc3ac..6db7ad019a 100644
--- a/flopy/modflow/mfsip.py
+++ b/flopy/modflow/mfsip.py
@@ -271,7 +271,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowSip.ftype()
             )
 
-        sip = ModflowSip(
+        return cls(
             model,
             mxiter=mxiter,
             nparm=nparm,
@@ -283,7 +283,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return sip
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfsms.py b/flopy/modflow/mfsms.py
index 5ec5cda8d8..738ab63446 100644
--- a/flopy/modflow/mfsms.py
+++ b/flopy/modflow/mfsms.py
@@ -629,7 +629,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowSms.ftype()
             )
 
-        sms = ModflowSms(
+        return cls(
             model,
             hclose=hclose,
             hiclose=hiclose,
@@ -664,7 +664,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return sms
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfsor.py b/flopy/modflow/mfsor.py
index 7c0e389679..fd91673f29 100644
--- a/flopy/modflow/mfsor.py
+++ b/flopy/modflow/mfsor.py
@@ -212,11 +212,8 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowSor.ftype()
             )
 
-        # create sor object
-        sor = ModflowSor(model, unitnumber=unitnumber, filenames=filenames)
-
         # return sor object
-        return sor
+        return cls(model, unitnumber=unitnumber, filenames=filenames)
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfstr.py b/flopy/modflow/mfstr.py
index b288445ef6..65c69e85ed 100644
--- a/flopy/modflow/mfstr.py
+++ b/flopy/modflow/mfstr.py
@@ -994,7 +994,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
                     ext_unit_dict, unit=abs(istcb2)
                 )
 
-        strpak = ModflowStr(
+        return cls(
             model,
             mxacts=mxacts,
             nss=nss,
@@ -1012,7 +1012,6 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return strpak
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfsub.py b/flopy/modflow/mfsub.py
index 12a8cdb1ed..10e9a13130 100644
--- a/flopy/modflow/mfsub.py
+++ b/flopy/modflow/mfsub.py
@@ -928,8 +928,8 @@ def load(cls, f, model, ext_unit_dict=None):
                         model.add_pop_key_list(unit)
                     ipos += 1
 
-        # create sub instance
-        sub = ModflowSub(
+        # return sub instance
+        return cls(
             model,
             ipakcb=ipakcb,
             isuboc=isuboc,
@@ -961,8 +961,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        # return sub instance
-        return sub
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfswi2.py b/flopy/modflow/mfswi2.py
index 44beea98e8..80b6a2937a 100644
--- a/flopy/modflow/mfswi2.py
+++ b/flopy/modflow/mfswi2.py
@@ -782,8 +782,8 @@ def load(cls, f, model, ext_unit_dict=None):
                     ext_unit_dict, unit=abs(iswiobs)
                 )
 
-        # create swi2 instance
-        swi2 = ModflowSwi2(
+        # return swi2 instance
+        return cls(
             model,
             nsrf=nsrf,
             istrat=istrat,
@@ -812,9 +812,6 @@ def load(cls, f, model, ext_unit_dict=None):
             filenames=filenames,
         )
 
-        # return swi2 instance
-        return swi2
-
     @staticmethod
     def ftype():
         return "SWI2"
diff --git a/flopy/modflow/mfswr1.py b/flopy/modflow/mfswr1.py
index 8c1f2379d9..a11cf5a7ad 100644
--- a/flopy/modflow/mfswr1.py
+++ b/flopy/modflow/mfswr1.py
@@ -184,11 +184,8 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowSwr1.ftype()
             )
 
-        # create swr1 object instance
-        swr1 = ModflowSwr1(model, unitnumber=unitnumber, filenames=filenames)
-
         # return swr object
-        return swr1
+        return cls(model, unitnumber=unitnumber, filenames=filenames)
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflow/mfswt.py b/flopy/modflow/mfswt.py
index e104aa55df..5984469863 100644
--- a/flopy/modflow/mfswt.py
+++ b/flopy/modflow/mfswt.py
@@ -931,8 +931,8 @@ def load(cls, f, model, ext_unit_dict=None):
                         model.add_pop_key_list(unit)
                     ipos += 1
 
-        # create sub-wt instance
-        swt = ModflowSwt(
+        # return sut-wt instance
+        return cls(
             model,
             ipakcb=ipakcb,
             iswtoc=iswtoc,
@@ -970,9 +970,6 @@ def load(cls, f, model, ext_unit_dict=None):
             filenames=filenames,
         )
 
-        # return sut-wt instance
-        return swt
-
     @staticmethod
     def ftype():
         return "SWT"
diff --git a/flopy/modflow/mfupw.py b/flopy/modflow/mfupw.py
index 2fe710e58f..0e8a071eae 100644
--- a/flopy/modflow/mfupw.py
+++ b/flopy/modflow/mfupw.py
@@ -609,7 +609,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                 model.add_pop_key_list(ipakcb)
 
         # create upw object
-        upw = ModflowUpw(
+        upw = cls(
             model,
             ipakcb=ipakcb,
             iphdry=iphdry,
diff --git a/flopy/modflow/mfuzf1.py b/flopy/modflow/mfuzf1.py
index fae888810b..49cf53d811 100644
--- a/flopy/modflow/mfuzf1.py
+++ b/flopy/modflow/mfuzf1.py
@@ -1131,7 +1131,7 @@ def load_util2d(name, dtype, per=None):
                     ipos += 1
 
         # create uzf object
-        return ModflowUzf1(
+        return cls(
             model,
             nuztop=nuztop,
             iuzfopt=iuzfopt,
diff --git a/flopy/modflow/mfwel.py b/flopy/modflow/mfwel.py
index c60b1d1cba..d40905d15c 100644
--- a/flopy/modflow/mfwel.py
+++ b/flopy/modflow/mfwel.py
@@ -427,7 +427,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         return Package.load(
             f,
             model,
-            ModflowWel,
+            cls,
             nper=nper,
             check=check,
             ext_unit_dict=ext_unit_dict,
diff --git a/flopy/modflow/mfzon.py b/flopy/modflow/mfzon.py
index 16db6a8e35..634e8ee94d 100644
--- a/flopy/modflow/mfzon.py
+++ b/flopy/modflow/mfzon.py
@@ -224,13 +224,12 @@ def load(cls, f, model, nrow=None, ncol=None, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowZon.ftype()
             )
 
-        zon = ModflowZon(
+        return cls(
             model,
             zone_dict=zone_dict,
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return zon
 
     @staticmethod
     def ftype():
diff --git a/flopy/modflowlgr/mflgr.py b/flopy/modflowlgr/mflgr.py
index a0bc937f0b..c48bae16b8 100644
--- a/flopy/modflowlgr/mflgr.py
+++ b/flopy/modflowlgr/mflgr.py
@@ -694,7 +694,8 @@ def load(
         if openfile:
             f.close()
 
-        lgr = ModflowLgr(
+        # return model object
+        return cls(
             version=version,
             exe_name=exe_name,
             modelname=modelname,
@@ -706,6 +707,3 @@ def load(
             children=children,
             children_data=children_data,
         )
-
-        # return model object
-        return lgr
diff --git a/flopy/mt3d/mt.py b/flopy/mt3d/mt.py
index 5c33a63147..8682bb3e0f 100644
--- a/flopy/mt3d/mt.py
+++ b/flopy/mt3d/mt.py
@@ -649,7 +649,7 @@ def load(
                     modelname, 50 * "-"
                 )
             )
-        mt = Mt3dms(
+        mt = cls(
             modelname=modelname,
             namefile_ext=modelname_extension,
             version=version,
diff --git a/flopy/mt3d/mtadv.py b/flopy/mt3d/mtadv.py
index efb61a35db..b97d0ac321 100644
--- a/flopy/mt3d/mtadv.py
+++ b/flopy/mt3d/mtadv.py
@@ -417,7 +417,7 @@ def load(cls, f, model, ext_unit_dict=None):
             )
 
         # Construct and return adv package
-        adv = Mt3dAdv(
+        return cls(
             model,
             mixelm=mixelm,
             percel=percel,
@@ -437,7 +437,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return adv
 
     @staticmethod
     def ftype():
diff --git a/flopy/mt3d/mtbtn.py b/flopy/mt3d/mtbtn.py
index 9b4ebf878f..e7c775fbb5 100644
--- a/flopy/mt3d/mtbtn.py
+++ b/flopy/mt3d/mtbtn.py
@@ -1242,7 +1242,7 @@ def load(cls, f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=Mt3dBtn.ftype()
             )
 
-        btn = Mt3dBtn(
+        return cls(
             model,
             MFStyleArr=MFStyleArr,
             DRYCell=DRYCell,
@@ -1293,7 +1293,6 @@ def load(cls, f, model, ext_unit_dict=None):
             filenames=filenames,
             **kwargs
         )
-        return btn
 
     @staticmethod
     def ftype():
diff --git a/flopy/mt3d/mtdsp.py b/flopy/mt3d/mtdsp.py
index 88def2e85d..a9437268f5 100644
--- a/flopy/mt3d/mtdsp.py
+++ b/flopy/mt3d/mtdsp.py
@@ -443,7 +443,7 @@ def load(
                 ext_unit_dict, filetype=Mt3dDsp.ftype()
             )
 
-        dsp = Mt3dDsp(
+        return cls(
             model,
             al=al,
             trpt=trpt,
@@ -454,7 +454,6 @@ def load(
             filenames=filenames,
             **kwargs
         )
-        return dsp
 
     @staticmethod
     def ftype():
diff --git a/flopy/mt3d/mtgcg.py b/flopy/mt3d/mtgcg.py
index 6b0a4fc793..5b75c6b621 100644
--- a/flopy/mt3d/mtgcg.py
+++ b/flopy/mt3d/mtgcg.py
@@ -240,7 +240,7 @@ def load(cls, f, model, ext_unit_dict=None):
             )
 
         # Construct and return gcg package
-        gcg = Mt3dGcg(
+        return cls(
             model,
             mxiter=mxiter,
             iter1=iter1,
@@ -252,7 +252,6 @@ def load(cls, f, model, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return gcg
 
     @staticmethod
     def ftype():
diff --git a/flopy/mt3d/mtlkt.py b/flopy/mt3d/mtlkt.py
index 489e1d6cd2..455d063275 100644
--- a/flopy/mt3d/mtlkt.py
+++ b/flopy/mt3d/mtlkt.py
@@ -512,7 +512,7 @@ def load(
                 model.add_pop_key_list(icbclk)
 
         # Construct and return LKT package
-        lkt = Mt3dLkt(
+        return cls(
             model,
             nlkinit=nlkinit,
             mxlkbc=mxlkbc,
@@ -524,7 +524,6 @@ def load(
             filenames=filenames,
             **kwargs
         )
-        return lkt
 
     @staticmethod
     def get_default_dtype(ncomp=1):
diff --git a/flopy/mt3d/mtrct.py b/flopy/mt3d/mtrct.py
index b603e5ac90..d8a043e07d 100644
--- a/flopy/mt3d/mtrct.py
+++ b/flopy/mt3d/mtrct.py
@@ -799,7 +799,7 @@ def load(
             )
 
         # Construct and return rct package
-        rct = Mt3dRct(
+        return cls(
             model,
             isothm=isothm,
             ireact=ireact,
@@ -815,7 +815,6 @@ def load(
             filenames=filenames,
             **kwargs
         )
-        return rct
 
     @staticmethod
     def ftype():
diff --git a/flopy/mt3d/mtsft.py b/flopy/mt3d/mtsft.py
index 2c4be25de8..7074a72199 100644
--- a/flopy/mt3d/mtsft.py
+++ b/flopy/mt3d/mtsft.py
@@ -802,7 +802,7 @@ def load(
                 model.add_pop_key_list(abs(ioutobs))
 
         # Construct and return SFT package
-        sft = Mt3dSft(
+        return cls(
             model,
             nsfinit=nsfinit,
             mxsfbc=mxsfbc,
@@ -824,7 +824,6 @@ def load(
             filenames=filenames,
             **kwargs
         )
-        return sft
 
     @staticmethod
     def ftype():
diff --git a/flopy/mt3d/mtssm.py b/flopy/mt3d/mtssm.py
index ec856bcf85..bc40a32bae 100644
--- a/flopy/mt3d/mtssm.py
+++ b/flopy/mt3d/mtssm.py
@@ -817,7 +817,7 @@ def load(
             )
 
         # Construct and return ssm package
-        ssm = Mt3dSsm(
+        return cls(
             model,
             crch=crch,
             cevt=cevt,
@@ -827,7 +827,6 @@ def load(
             filenames=filenames,
             **kwargs
         )
-        return ssm
 
     @staticmethod
     def ftype():
diff --git a/flopy/mt3d/mtuzt.py b/flopy/mt3d/mtuzt.py
index 4bda87a62d..79f674176d 100644
--- a/flopy/mt3d/mtuzt.py
+++ b/flopy/mt3d/mtuzt.py
@@ -779,7 +779,7 @@ def load(
                 model.add_pop_key_list(icbcuz)
 
         # Construct and return uzt package
-        uzt = Mt3dUzt(
+        return cls(
             model,
             icbcuz=icbcuz,
             iet=iet,
@@ -791,7 +791,6 @@ def load(
             filenames=filenames,
             **kwargs
         )
-        return uzt
 
     @staticmethod
     def ftype():
diff --git a/flopy/pakbase.py b/flopy/pakbase.py
index b491ccb7f7..519f3935be 100644
--- a/flopy/pakbase.py
+++ b/flopy/pakbase.py
@@ -849,8 +849,8 @@ def write_file(self, check=False):
         print("IMPLEMENTATION ERROR: write_file must be overloaded")
         return
 
-    @classmethod
-    def load(cls, f, model, pak_type, ext_unit_dict=None, **kwargs):
+    @staticmethod
+    def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
         """
         Default load method for standard boundary packages.
 
diff --git a/flopy/seawat/swt.py b/flopy/seawat/swt.py
index 9ab6cd632b..b13f9bb8ff 100644
--- a/flopy/seawat/swt.py
+++ b/flopy/seawat/swt.py
@@ -489,7 +489,7 @@ def load(
             modelname = f
 
         # create instance of a seawat model and load modflow and mt3dms models
-        ms = Seawat(
+        ms = cls(
             modelname=modelname,
             namefile_ext="nam",
             modflowmodel=None,
diff --git a/flopy/seawat/swtvdf.py b/flopy/seawat/swtvdf.py
index 40d714b303..c7021b427a 100644
--- a/flopy/seawat/swtvdf.py
+++ b/flopy/seawat/swtvdf.py
@@ -540,7 +540,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
             )
 
         # Construct and return vdf package
-        vdf = SeawatVdf(
+        return cls(
             model,
             mt3drhoflg=mt3drhoflg,
             mfnadvfd=mfnadvfd,
@@ -562,7 +562,6 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return vdf
 
     @staticmethod
     def ftype():
diff --git a/flopy/seawat/swtvsc.py b/flopy/seawat/swtvsc.py
index 40da05c6db..c4ece4595b 100644
--- a/flopy/seawat/swtvsc.py
+++ b/flopy/seawat/swtvsc.py
@@ -508,7 +508,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
             )
 
         # Construct and return vsc package
-        vsc = SeawatVsc(
+        return cls(
             model,
             mt3dmuflg=mt3dmuflg,
             viscmin=viscmin,
@@ -526,7 +526,6 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
             unitnumber=unitnumber,
             filenames=filenames,
         )
-        return vsc
 
     @staticmethod
     def ftype():
diff --git a/flopy/utils/optionblock.py b/flopy/utils/optionblock.py
index c0884e9af2..5926414f3d 100644
--- a/flopy/utils/optionblock.py
+++ b/flopy/utils/optionblock.py
@@ -410,7 +410,7 @@ def load_options(cls, options, package):
             else:
                 if openfile:
                     options.close()
-                return OptionBlock(options_line=option_line, package=package)
+                return cls(options_line=option_line, package=package)
 
 
 class OptionUtil(object):
diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py
index f3a1d02f10..7318398408 100644
--- a/flopy/utils/rasters.py
+++ b/flopy/utils/rasters.py
@@ -783,7 +783,7 @@ def load(cls, raster):
         bands = dataset.indexes
         meta = dataset.meta
 
-        return Raster(
+        return cls(
             array,
             bands,
             meta["crs"],
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 165ad24cfa..7d01c1f893 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -618,7 +618,7 @@ def from_namfile(cls, namefile):
         except:
             print("   could not remove start_datetime")
 
-        return SpatialReference(**attribs)
+        return cls(**attribs)
 
     @classmethod
     def from_gridspec(cls, gridspec_file, lenuni=0):
diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py
index 7d76ba7935..8be8be6a4b 100644
--- a/flopy/utils/util_array.py
+++ b/flopy/utils/util_array.py
@@ -968,8 +968,7 @@ def load(
                 array_format=array_format,
             )
             u2ds.append(u2d)
-        u3d = Util3d(model, shape, dtype, u2ds, name)
-        return u3d
+        return cls(model, shape, dtype, u2ds, name)
 
     def __mul__(self, other):
         if np.isscalar(other):
@@ -3028,7 +3027,7 @@ def load(
         )
 
         if cr_dict["type"] == "constant":
-            u2d = Util2d(
+            u2d = cls(
                 model,
                 shape,
                 dtype,
@@ -3066,7 +3065,7 @@ def load(
                     shape, f, dtype, bintype="Head"
                 )
             f.close()
-            u2d = Util2d(
+            u2d = cls(
                 model,
                 shape,
                 dtype,
@@ -3080,7 +3079,7 @@ def load(
 
         elif cr_dict["type"] == "internal":
             data = Util2d.load_txt(shape, f_handle, dtype, cr_dict["fmtin"])
-            u2d = Util2d(
+            u2d = cls(
                 model,
                 shape,
                 dtype,
@@ -3113,7 +3112,7 @@ def load(
                 header_data, data = Util2d.load_bin(
                     shape, ext_unit.filehandle, dtype, bintype="Head"
                 )
-            u2d = Util2d(
+            u2d = cls(
                 model,
                 shape,
                 dtype,
@@ -3129,7 +3128,7 @@ def load(
             model.pop_key_list.append(cr_dict["nunit"])
         elif cr_dict["type"] == "block":
             data = Util2d.load_block(shape, f_handle, dtype)
-            u2d = Util2d(
+            u2d = cls(
                 model,
                 shape,
                 dtype,