diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py
index f8dea2e496..6c56652243 100644
--- a/flopy/modflow/mf.py
+++ b/flopy/modflow/mf.py
@@ -670,8 +670,9 @@ def load_results(self, **kwargs):
         else:
             return hdObj, ddObj, bdObj
 
-    @staticmethod
+    @classmethod
     def load(
+        cls,
         f,
         version="mf2005",
         exe_name="mf2005.exe",
@@ -746,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 e99f59bf16..00d812e788 100644
--- a/flopy/modflow/mfag.py
+++ b/flopy/modflow/mfag.py
@@ -685,8 +685,8 @@ def get_default_dtype(maxells=0, block="well"):
 
         return np.dtype(dtype)
 
-    @staticmethod
-    def load(f, model, nper=0, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=0, ext_unit_dict=None):
         """
         Method to load the AG package from file
 
@@ -891,7 +891,7 @@ def load(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 e22a1f95e7..30e502c722 100644
--- a/flopy/modflow/mfbas.py
+++ b/flopy/modflow/mfbas.py
@@ -286,8 +286,8 @@ def write_file(self, check=True):
         # Close file
         f_bas.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=True, **kwargs):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=True, **kwargs):
         """
         Load an existing package.
 
@@ -398,7 +398,7 @@ def load(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 a05dd1ece2..baaa2c8cce 100644
--- a/flopy/modflow/mfbcf.py
+++ b/flopy/modflow/mfbcf.py
@@ -320,8 +320,8 @@ def write_file(self, f=None):
                 f_bcf.write(self.wetdry[k].get_file_entry())
         f_bcf.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -529,7 +529,7 @@ def load(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/mfde4.py b/flopy/modflow/mfde4.py
index 95cce31b57..09be4c0e1f 100644
--- a/flopy/modflow/mfde4.py
+++ b/flopy/modflow/mfde4.py
@@ -227,8 +227,8 @@ def write_file(self):
             f.write("\n")
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -313,7 +313,7 @@ def load(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 210be95a23..cdfa7f86e8 100644
--- a/flopy/modflow/mfdis.py
+++ b/flopy/modflow/mfdis.py
@@ -848,8 +848,8 @@ def check(self, f=None, verbose=True, level=1, checktype=None):
         # if verbose:
         #     print(txt)
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=True):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=True):
         """
         Load an existing package.
 
@@ -1067,7 +1067,7 @@ def load(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 917fa70b07..329767d0b5 100644
--- a/flopy/modflow/mfdisu.py
+++ b/flopy/modflow/mfdisu.py
@@ -529,8 +529,8 @@ def zcentroids(self):
     def ncpl(self):
         return self.nodes / self.nlay
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=False):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=False):
         """
         Load an existing package.
 
@@ -798,7 +798,7 @@ def load(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/mfevt.py b/flopy/modflow/mfevt.py
index baee297ec5..c7842c5284 100644
--- a/flopy/modflow/mfevt.py
+++ b/flopy/modflow/mfevt.py
@@ -215,8 +215,8 @@ def write_file(self, f=None):
                 f_evt.write(ievt)
         f_evt.close()
 
-    @staticmethod
-    def load(f, model, nper=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -417,7 +417,7 @@ def load(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 adb7b5e5c8..3a7621e5ab 100644
--- a/flopy/modflow/mffhb.py
+++ b/flopy/modflow/mffhb.py
@@ -420,8 +420,8 @@ def write_file(self):
 
         f.close()
 
-    @staticmethod
-    def load(f, model, nper=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -756,7 +756,7 @@ def load(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 afb791df75..f4642b9159 100755
--- a/flopy/modflow/mfflwob.py
+++ b/flopy/modflow/mfflwob.py
@@ -380,8 +380,8 @@ def write_file(self):
 
         return
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=True):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=True):
         """
         Load an existing package.
 
@@ -537,7 +537,7 @@ def load(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 3afd728ff4..19e739ea10 100644
--- a/flopy/modflow/mfgage.py
+++ b/flopy/modflow/mfgage.py
@@ -286,8 +286,8 @@ def write_file(self):
         # close the gage file
         f.close()
 
-    @staticmethod
-    def load(f, model, nper=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -395,7 +395,7 @@ def load(f, model, nper=None, ext_unit_dict=None):
         for file in files:
             filenames.append(os.path.basename(file))
 
-        gagepak = ModflowGage(
+        gagepak = cls(
             model,
             numgage=numgage,
             gage_data=gage_data,
diff --git a/flopy/modflow/mfgmg.py b/flopy/modflow/mfgmg.py
index 187c8f635b..99355317c3 100644
--- a/flopy/modflow/mfgmg.py
+++ b/flopy/modflow/mfgmg.py
@@ -318,8 +318,8 @@ def write_file(self):
         f_gmg.write("{}\n".format(self.relax))
         f_gmg.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -411,7 +411,7 @@ def load(f, model, ext_unit_dict=None):
                 model.add_pop_key_list(iunitmhc)
 
         # create the gmg object
-        gmg = ModflowGmg(
+        gmg = cls(
             model,
             mxiter=mxiter,
             iiter=iiter,
diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py
index 0b88a384a5..3fdf14ace7 100644
--- a/flopy/modflow/mfhfb.py
+++ b/flopy/modflow/mfhfb.py
@@ -253,8 +253,8 @@ def get_default_dtype(structured=True):
     def get_sfac_columns():
         return ["hydchr"]
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -412,7 +412,7 @@ def load(f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowHfb.ftype()
             )
 
-        hfb = ModflowHfb(
+        hfb = cls(
             model,
             nphfb=0,
             mxfb=0,
diff --git a/flopy/modflow/mfhob.py b/flopy/modflow/mfhob.py
index 89722b41cc..769fb2f654 100755
--- a/flopy/modflow/mfhob.py
+++ b/flopy/modflow/mfhob.py
@@ -308,8 +308,8 @@ def write_file(self):
 
         return
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=True):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=True):
         """
         Load an existing package.
 
@@ -492,7 +492,7 @@ def load(f, model, ext_unit_dict=None, check=True):
                     model.add_pop_key_list(iuhobsv)
 
         # create hob object instance
-        hob = ModflowHob(
+        hob = cls(
             model,
             iuhobsv=iuhobsv,
             hobdry=hobdry,
diff --git a/flopy/modflow/mfhyd.py b/flopy/modflow/mfhyd.py
index 29c2a48273..a52d6d1a35 100644
--- a/flopy/modflow/mfhyd.py
+++ b/flopy/modflow/mfhyd.py
@@ -286,8 +286,8 @@ def get_default_dtype():
         )
         return dtype
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -367,7 +367,7 @@ def load(f, model, ext_unit_dict=None):
                 model.add_pop_key_list(ihydun)
 
         # create hyd instance
-        hyd = ModflowHyd(
+        hyd = cls(
             model,
             nhyd=nhyd,
             ihydun=ihydun,
diff --git a/flopy/modflow/mflak.py b/flopy/modflow/mflak.py
index 114d3a1494..268e159c79 100644
--- a/flopy/modflow/mflak.py
+++ b/flopy/modflow/mflak.py
@@ -639,8 +639,8 @@ def write_file(self):
         # close the lak file
         f.close()
 
-    @staticmethod
-    def load(f, model, nper=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -908,7 +908,7 @@ def load(f, model, nper=None, ext_unit_dict=None):
                     )
                     ipos += 1
 
-        lakpak = ModflowLak(
+        lakpak = cls(
             model,
             options=options,
             nlakes=nlakes,
diff --git a/flopy/modflow/mflmt.py b/flopy/modflow/mflmt.py
index ed905a6891..229363a629 100644
--- a/flopy/modflow/mflmt.py
+++ b/flopy/modflow/mflmt.py
@@ -170,8 +170,8 @@ def write_file(self):
 
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -257,7 +257,7 @@ def load(f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowLmt.ftype()
             )
 
-        lmt = ModflowLmt(
+        lmt = cls(
             model,
             output_file_name=output_file_name,
             output_file_unit=output_file_unit,
diff --git a/flopy/modflow/mflpf.py b/flopy/modflow/mflpf.py
index cd03e7507e..e2dfdb4314 100644
--- a/flopy/modflow/mflpf.py
+++ b/flopy/modflow/mflpf.py
@@ -437,8 +437,8 @@ def write_file(self, check=True, f=None):
         f.close()
         return
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=True):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=True):
         """
         Load an existing package.
 
@@ -723,7 +723,7 @@ def load(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 60392397fc..418c4967c2 100644
--- a/flopy/modflow/mfmlt.py
+++ b/flopy/modflow/mfmlt.py
@@ -128,8 +128,8 @@ def write_file(self):
         """
         pass
 
-    @staticmethod
-    def load(f, model, nrow=None, ncol=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nrow=None, ncol=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -233,7 +233,7 @@ def load(f, model, nrow=None, ncol=None, ext_unit_dict=None):
             )
 
         # create mlt dictionary
-        mlt = ModflowMlt(
+        mlt = cls(
             model,
             mult_dict=mult_dict,
             unitnumber=unitnumber,
diff --git a/flopy/modflow/mfmnw1.py b/flopy/modflow/mfmnw1.py
index 56bb27f0f5..75b8ffabc7 100644
--- a/flopy/modflow/mfmnw1.py
+++ b/flopy/modflow/mfmnw1.py
@@ -209,8 +209,8 @@ def get_default_dtype(structured=True):
         else:
             pass
 
-    @staticmethod
-    def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
 
         if model.verbose:
             sys.stdout.write("loading mnw1 package file...\n")
@@ -278,7 +278,7 @@ def load(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 1d18207a55..d7f1603747 100644
--- a/flopy/modflow/mfmnw2.py
+++ b/flopy/modflow/mfmnw2.py
@@ -1324,8 +1324,8 @@ def get_default_spd_dtype(structured=True):
             msg = "get_default_spd_dtype: unstructured model not supported"
             raise NotImplementedError(msg)
 
-    @staticmethod
-    def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
         """
 
         Parameters
@@ -1464,7 +1464,7 @@ def load(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 016f054951..54e69090f5 100644
--- a/flopy/modflow/mfmnwi.py
+++ b/flopy/modflow/mfmnwi.py
@@ -202,8 +202,8 @@ def __init__(
 
         self.parent.add_package(self)
 
-    @staticmethod
-    def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
 
         if model.verbose:
             sys.stdout.write("loading mnw2 package file...\n")
@@ -281,7 +281,7 @@ def load(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 c798a95599..ee7a4e2ab2 100644
--- a/flopy/modflow/mfnwt.py
+++ b/flopy/modflow/mfnwt.py
@@ -377,8 +377,8 @@ def write_file(self):
 
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -556,7 +556,7 @@ def load(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 db7e6a7e42..f943e2a850 100644
--- a/flopy/modflow/mfoc.py
+++ b/flopy/modflow/mfoc.py
@@ -733,8 +733,8 @@ def get_ocoutput_units(f, ext_unit_dict=None):
         # return
         return ihedun, fhead, iddnun, fddn
 
-    @staticmethod
-    def load(f, model, nper=None, nstp=None, nlay=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, nstp=None, nlay=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -1103,7 +1103,7 @@ def load(f, model, nper=None, nstp=None, nlay=None, ext_unit_dict=None):
             model.add_pop_key_list(u)
 
         # create instance of oc class
-        oc = ModflowOc(
+        oc = cls(
             model,
             ihedfm=ihedfm,
             iddnfm=iddnfm,
diff --git a/flopy/modflow/mfparbc.py b/flopy/modflow/mfparbc.py
index 4dd87f726a..f24d2f8065 100644
--- a/flopy/modflow/mfparbc.py
+++ b/flopy/modflow/mfparbc.py
@@ -41,8 +41,8 @@ def get(self, fkey):
                 return self.bc_parms[key]
         return None
 
-    @staticmethod
-    def load(f, npar, dt, model, ext_unit_dict=None, verbose=False):
+    @classmethod
+    def load(cls, f, npar, dt, model, ext_unit_dict=None, verbose=False):
         """
         Load bc property parameters from an existing bc package
         that uses list data (e.g. WEL, RIV, etc.).
@@ -119,7 +119,7 @@ def load(f, npar, dt, model, ext_unit_dict=None, verbose=False):
                 ]
 
         # print bc_parms
-        bcpar = ModflowParBc(bc_parms)
+        bcpar = cls(bc_parms)
         return bcpar
 
     @staticmethod
diff --git a/flopy/modflow/mfpcg.py b/flopy/modflow/mfpcg.py
index 19ebc76ebe..ead674e6b0 100644
--- a/flopy/modflow/mfpcg.py
+++ b/flopy/modflow/mfpcg.py
@@ -219,8 +219,8 @@ def write_file(self):
             f.write("\n")
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -343,7 +343,7 @@ def load(f, model, ext_unit_dict=None):
             )
 
         # create instance of pcg class
-        pcg = ModflowPcg(
+        pcg = cls(
             model,
             mxiter=mxiter,
             iter1=iter1,
diff --git a/flopy/modflow/mfpcgn.py b/flopy/modflow/mfpcgn.py
index c8289a372d..1fa090b01e 100644
--- a/flopy/modflow/mfpcgn.py
+++ b/flopy/modflow/mfpcgn.py
@@ -393,8 +393,8 @@ def write_file(self):
             f.write(line)
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -558,7 +558,7 @@ def load(f, model, ext_unit_dict=None):
                     ext_unit_dict, unit=ipunit
                 )
 
-        pcgn = ModflowPcgn(
+        pcgn = cls(
             model,
             iter_mo=iter_mo,
             iter_mi=iter_mi,
diff --git a/flopy/modflow/mfpks.py b/flopy/modflow/mfpks.py
index 38ce08fdae..437fcd44a2 100644
--- a/flopy/modflow/mfpks.py
+++ b/flopy/modflow/mfpks.py
@@ -228,8 +228,8 @@ def write_file(self):
         f.write("END\n")
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -287,7 +287,7 @@ def load(f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowPks.ftype()
             )
 
-        pks = ModflowPks(model, unitnumber=unitnumber, filenames=filenames)
+        pks = cls(model, unitnumber=unitnumber, filenames=filenames)
         return pks
 
     @staticmethod
diff --git a/flopy/modflow/mfpval.py b/flopy/modflow/mfpval.py
index 6a7c516b21..5893671f81 100644
--- a/flopy/modflow/mfpval.py
+++ b/flopy/modflow/mfpval.py
@@ -139,8 +139,8 @@ def __getitem__(self, item):
         else:
             return None
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -218,7 +218,7 @@ def load(f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowPval.ftype()
             )
 
-        pval = ModflowPval(
+        pval = cls(
             model,
             pval_dict=pval_dict,
             unitnumber=unitnumber,
diff --git a/flopy/modflow/mfrch.py b/flopy/modflow/mfrch.py
index 242b54bdf4..49b7ce567e 100644
--- a/flopy/modflow/mfrch.py
+++ b/flopy/modflow/mfrch.py
@@ -366,8 +366,8 @@ def write_file(self, check=True, f=None):
                     f_rch.write(file_entry_irch)
         f_rch.close()
 
-    @staticmethod
-    def load(f, model, nper=None, ext_unit_dict=None, check=True):
+    @classmethod
+    def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         """
         Load an existing package.
 
@@ -527,7 +527,7 @@ def load(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/mfsfr2.py b/flopy/modflow/mfsfr2.py
index 7bd24930d7..23ededc759 100644
--- a/flopy/modflow/mfsfr2.py
+++ b/flopy/modflow/mfsfr2.py
@@ -834,8 +834,8 @@ def get_default_segment_dtype():
             ]
         )
 
-    @staticmethod
-    def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
 
         if model.verbose:
             sys.stdout.write("loading sfr2 package file...\n")
@@ -1070,7 +1070,7 @@ def load(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 e851a72b69..c3a73b7e05 100644
--- a/flopy/modflow/mfsip.py
+++ b/flopy/modflow/mfsip.py
@@ -190,8 +190,8 @@ def write_file(self):
             )
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -271,7 +271,7 @@ def load(f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowSip.ftype()
             )
 
-        sip = ModflowSip(
+        sip = cls(
             model,
             mxiter=mxiter,
             nparm=nparm,
diff --git a/flopy/modflow/mfsms.py b/flopy/modflow/mfsms.py
index f411e01914..c9b097f4ab 100644
--- a/flopy/modflow/mfsms.py
+++ b/flopy/modflow/mfsms.py
@@ -412,8 +412,8 @@ def write_file(self):
         f.write("\n")
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -629,7 +629,7 @@ def load(f, model, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowSms.ftype()
             )
 
-        sms = ModflowSms(
+        sms = cls(
             model,
             hclose=hclose,
             hiclose=hiclose,
diff --git a/flopy/modflow/mfsor.py b/flopy/modflow/mfsor.py
index 7330d66bff..898b590be1 100644
--- a/flopy/modflow/mfsor.py
+++ b/flopy/modflow/mfsor.py
@@ -152,8 +152,8 @@ def write_file(self):
         f.write(line)
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -213,7 +213,7 @@ def load(f, model, ext_unit_dict=None):
             )
 
         # create sor object
-        sor = ModflowSor(model, unitnumber=unitnumber, filenames=filenames)
+        sor = cls(model, unitnumber=unitnumber, filenames=filenames)
 
         # return sor object
         return sor
diff --git a/flopy/modflow/mfstr.py b/flopy/modflow/mfstr.py
index 152f50dcb2..a7394b4663 100644
--- a/flopy/modflow/mfstr.py
+++ b/flopy/modflow/mfstr.py
@@ -681,8 +681,8 @@ def write_file(self):
         # close the str file
         f_str.close()
 
-    @staticmethod
-    def load(f, model, nper=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nper=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -994,7 +994,7 @@ def load(f, model, nper=None, ext_unit_dict=None):
                     ext_unit_dict, unit=abs(istcb2)
                 )
 
-        strpak = ModflowStr(
+        strpak = cls(
             model,
             mxacts=mxacts,
             nss=nss,
diff --git a/flopy/modflow/mfsub.py b/flopy/modflow/mfsub.py
index 9f386ad793..b1d5bf987d 100644
--- a/flopy/modflow/mfsub.py
+++ b/flopy/modflow/mfsub.py
@@ -586,8 +586,8 @@ def write_file(self, check=False, f=None):
         # close sub file
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -929,7 +929,7 @@ def load(f, model, ext_unit_dict=None):
                     ipos += 1
 
         # create sub instance
-        sub = ModflowSub(
+        sub = cls(
             model,
             ipakcb=ipakcb,
             isuboc=isuboc,
diff --git a/flopy/modflow/mfswi2.py b/flopy/modflow/mfswi2.py
index bc2e87b728..89998a27d2 100644
--- a/flopy/modflow/mfswi2.py
+++ b/flopy/modflow/mfswi2.py
@@ -531,8 +531,8 @@ def write_file(self, check=True, f=None):
         # close swi2 file
         f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """Load an existing package.
 
         Parameters
@@ -783,7 +783,7 @@ def load(f, model, ext_unit_dict=None):
                 )
 
         # create swi2 instance
-        swi2 = ModflowSwi2(
+        swi2 = cls(
             model,
             nsrf=nsrf,
             istrat=istrat,
diff --git a/flopy/modflow/mfswr1.py b/flopy/modflow/mfswr1.py
index d8f3288cfe..3dc85226ed 100644
--- a/flopy/modflow/mfswr1.py
+++ b/flopy/modflow/mfswr1.py
@@ -122,8 +122,8 @@ def write_file(self):
         # f.write('{0}\n'.format(self.heading))
         # f.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -185,7 +185,7 @@ def load(f, model, ext_unit_dict=None):
             )
 
         # create swr1 object instance
-        swr1 = ModflowSwr1(model, unitnumber=unitnumber, filenames=filenames)
+        swr1 = cls(model, unitnumber=unitnumber, filenames=filenames)
 
         # return swr object
         return swr1
diff --git a/flopy/modflow/mfswt.py b/flopy/modflow/mfswt.py
index e943591bb4..7d3113081b 100644
--- a/flopy/modflow/mfswt.py
+++ b/flopy/modflow/mfswt.py
@@ -599,8 +599,8 @@ def __init__(
         # add package to model
         self.parent.add_package(self)
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -932,7 +932,7 @@ def load(f, model, ext_unit_dict=None):
                     ipos += 1
 
         # create sub-wt instance
-        swt = ModflowSwt(
+        swt = cls(
             model,
             ipakcb=ipakcb,
             iswtoc=iswtoc,
diff --git a/flopy/modflow/mfupw.py b/flopy/modflow/mfupw.py
index 639395f2be..0e8a071eae 100644
--- a/flopy/modflow/mfupw.py
+++ b/flopy/modflow/mfupw.py
@@ -355,8 +355,8 @@ def write_file(self, check=True, f=None):
                 f_upw.write(self.laywet[k].get_file_entry())
         f_upw.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=True):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=True):
         """
         Load an existing package.
 
@@ -609,7 +609,7 @@ def load(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 48695ac0f3..49cf53d811 100644
--- a/flopy/modflow/mfuzf1.py
+++ b/flopy/modflow/mfuzf1.py
@@ -873,8 +873,8 @@ def write_transient(name):
 
         f_uzf.close()
 
-    @staticmethod
-    def load(f, model, ext_unit_dict=None, check=False):
+    @classmethod
+    def load(cls, f, model, ext_unit_dict=None, check=False):
         """
         Load an existing package.
 
@@ -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/mfzon.py b/flopy/modflow/mfzon.py
index 64147512c4..ed954958d0 100644
--- a/flopy/modflow/mfzon.py
+++ b/flopy/modflow/mfzon.py
@@ -131,8 +131,8 @@ def write_file(self):
         """
         return
 
-    @staticmethod
-    def load(f, model, nrow=None, ncol=None, ext_unit_dict=None):
+    @classmethod
+    def load(cls, f, model, nrow=None, ncol=None, ext_unit_dict=None):
         """
         Load an existing package.
 
@@ -224,7 +224,7 @@ def load(f, model, nrow=None, ncol=None, ext_unit_dict=None):
                 ext_unit_dict, filetype=ModflowZon.ftype()
             )
 
-        zon = ModflowZon(
+        zon = cls(
             model,
             zone_dict=zone_dict,
             unitnumber=unitnumber,