diff --git a/autotest/t504_test.py b/autotest/t504_test.py
index d8411cead6..00e5617f00 100644
--- a/autotest/t504_test.py
+++ b/autotest/t504_test.py
@@ -622,6 +622,43 @@ def test006_2models_mvr():
         # clean up
         sim.delete_output_files()
 
+    # test load_only
+    model_package_check = ['ic', 'maw', 'npf', 'oc']
+    load_only_lists = [['ic6', 'npf6', 'oc', 'gwf6-gwf6', 'ims'],
+                       ['ic', 'maw', 'npf', 'gwf-gwf', 'ims'],
+                       ['ic', 'maw6', 'npf']]
+    for load_only in load_only_lists:
+        sim = MFSimulation.load(sim_name, 'mf6', exe_name, pth,
+                                load_only=load_only)
+        for model_name in model_names:
+            model = sim.get_model(model_name)
+            for package in model_package_check:
+                assert (package in model.package_type_dict or
+                        package in sim.package_type_dict) == \
+                       (package in load_only or '{}6'.format(package) in
+                        load_only)
+        assert (len(sim._exchange_files) > 0) == ('gwf6-gwf6' in load_only or
+                                                  'gwf-gwf' in load_only)
+        assert (len(sim._ims_files) > 0) == ('ims6' in load_only or
+                                             'ims' in load_only)
+
+    # load package by name
+    load_only_list = ['ic6', 'maw', 'npf_p1', 'oc_p2', 'ims']
+    sim = MFSimulation.load(sim_name, 'mf6', exe_name, pth,
+                            load_only=load_only_list)
+    model_parent = sim.get_model('parent')
+    model_child = sim.get_model('child')
+    assert 'oc' not in model_parent.package_type_dict
+    assert 'oc' in model_child.package_type_dict
+    assert 'npf' in model_parent.package_type_dict
+    assert 'npf' not in model_child.package_type_dict
+
+    if run:
+        # test running a runnable load_only case
+        sim = MFSimulation.load(sim_name, 'mf6', exe_name, pth,
+                                load_only=load_only_lists[0])
+        assert sim.run_simulation()[0]
+
     return
 
 
@@ -684,6 +721,26 @@ def test001e_uzf_3lay():
         outfile = os.path.join(save_folder, 'head_compare.dat')
         assert pymake.compare_heads(None, None, files1=head_file, files2=head_new, outfile=outfile)
 
+    # test load_only
+    model_package_check = ['chd', 'ic', 'npf', 'oc', 'sto', 'uzf']
+    load_only_lists = [['chd6', 'ic6', 'ims', 'npf6', 'obs', 'oc', 'sto'],
+                       ['chd6', 'ims', 'npf6', 'obs', 'oc', 'sto', 'uzf6'],
+                       ['chd', 'ic', 'npf', 'obs', 'sto'],
+                       ['ic6', 'ims', 'obs6', 'oc6']]
+    for load_only in load_only_lists:
+        sim = MFSimulation.load(model_name, 'mf6', exe_name, pth,
+                                load_only=load_only)
+        model = sim.get_model()
+        for package in model_package_check:
+            print(package)
+            assert (package in model.package_type_dict) == \
+                   (package in load_only or '{}6'.format(package) in load_only)
+    if run:
+        # test running a runnable load_only case
+        sim = MFSimulation.load(model_name, 'mf6', exe_name, pth,
+                                load_only=load_only_lists[0])
+        assert sim.run_simulation()[0]
+
 
 def test045_lake2tr():
     # init paths
@@ -699,8 +756,10 @@ def test045_lake2tr():
         os.makedirs(save_folder)
 
     expected_output_folder = os.path.join(pth, 'expected_output')
-    expected_head_file_a = os.path.join(expected_output_folder, 'lakeex2a_unch.hds')
-    expected_head_file_b = os.path.join(expected_output_folder, 'lakeex2a_adj.hds')
+    expected_head_file_a = os.path.join(expected_output_folder,
+                                        'lakeex2a_unch.hds')
+    expected_head_file_b = os.path.join(expected_output_folder,
+                                        'lakeex2a_adj.hds')
 
     # load simulation
     sim = MFSimulation.load(model_name, 'mf6', exe_name, pth)
@@ -868,12 +927,12 @@ def test027_timeseriestest():
 
 
 if __name__ == '__main__':
-    test006_gwf3()
     test001a_tharmonic()
     test001e_uzf_3lay()
     test003_gwfs_disv()
     test005_advgw_tidal()
     test006_2models_mvr()
+    test006_gwf3()
     test027_timeseriestest()
     test036_twrihfb()
     test045_lake1ss_table()
diff --git a/examples/data/mf6/test006_2models_mvr/model1.nam b/examples/data/mf6/test006_2models_mvr/model1.nam
index a2b9ae0a96..92ac07d816 100644
--- a/examples/data/mf6/test006_2models_mvr/model1.nam
+++ b/examples/data/mf6/test006_2models_mvr/model1.nam
@@ -5,9 +5,9 @@ END OPTIONS
 BEGIN PACKAGES
   DIS6           model1.dis
   IC6            model1.ic
-  NPF6           model1.npf
+  NPF6           model1.npf    npf_p1
   CHD6           model1.chd
   MAW6           model1.mawq
-  OC6            model1.oc
+  OC6            model1.oc    oc_p1
 END PACKAGES
 
diff --git a/examples/data/mf6/test006_2models_mvr/model2.nam b/examples/data/mf6/test006_2models_mvr/model2.nam
index b4982cf0b8..8783862ae2 100644
--- a/examples/data/mf6/test006_2models_mvr/model2.nam
+++ b/examples/data/mf6/test006_2models_mvr/model2.nam
@@ -7,6 +7,6 @@ BEGIN PACKAGES
   IC6            model2.ic
   NPF6           model2.npf
   MAW6           model2.mawq
-  OC6            model2.oc
+  OC6            model2.oc    oc_p2
 END PACKAGES
 
diff --git a/examples/data/mt3d_test/mf2kmt3d/HSSTest/hsstest.WEL b/examples/data/mt3d_test/mf2kmt3d/HSSTest/hsstest.WEL
index e9f3a5184d..d468925805 100644
--- a/examples/data/mt3d_test/mf2kmt3d/HSSTest/hsstest.WEL
+++ b/examples/data/mt3d_test/mf2kmt3d/HSSTest/hsstest.WEL
@@ -1,3 +1,4 @@
          2         0
+		 2
          1        16        16        1.
          1        16        31       -1.
diff --git a/flopy/mf6/mfbase.py b/flopy/mf6/mfbase.py
index a95eedc45b..74ff82ab2b 100644
--- a/flopy/mf6/mfbase.py
+++ b/flopy/mf6/mfbase.py
@@ -625,3 +625,46 @@ def get_package(self, name=None):
     def register_package(self, package):
         path = (package.package_name,)
         return (path, None)
+
+    @staticmethod
+    def _load_only_dict(load_only):
+        if load_only is None:
+            return None
+        if isinstance(load_only, dict):
+            return load_only
+        if not isinstance(load_only, collections.Iterable):
+            raise FlopyException('load_only must be iterable or None. '
+                                 'load_only value of "{}" is '
+                                 'invalid'.format(load_only))
+        load_only_dict = {}
+        for item in load_only:
+            load_only_dict[item.lower()] = True
+        return load_only_dict
+
+    @staticmethod
+    def _in_pkg_list(pkg_list, pkg_type, pkg_name):
+        if pkg_type is not None:
+            pkg_type = pkg_type.lower()
+        if pkg_name is not None:
+            pkg_name = pkg_name.lower()
+        if pkg_type in pkg_list or pkg_name in pkg_list:
+            return True
+
+        # split to make cases like "gwf6-gwf6" easier to process
+        pkg_type = pkg_type.split('-')
+        try:
+            # if there is a number on the end of the package try
+            # excluding it
+            int(pkg_type[0][-1])
+            for key in pkg_list.keys():
+                key = key.split('-')
+                if len(key) == len(pkg_type):
+                    matches = True
+                    for key_item, pkg_item in zip(key, pkg_type):
+                        if pkg_item[0:-1] != key_item and pkg_item != key_item:
+                            matches = False
+                    if matches:
+                        return True
+        except ValueError:
+            return False
+        return False
\ No newline at end of file
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
index 43fc099b57..49e0d1cbef 100644
--- a/flopy/mf6/mfmodel.py
+++ b/flopy/mf6/mfmodel.py
@@ -467,7 +467,8 @@ def verbose(self, verbose):
     @classmethod
     def load_base(cls, simulation, structure, modelname='NewModel',
                   model_nam_file='modflowtest.nam', mtype='gwf', version='mf6',
-                  exe_name='mf6.exe', strict=True, model_rel_path='.'):
+                  exe_name='mf6.exe', strict=True, model_rel_path='.',
+                  load_only=None):
         """
         Load an existing model.
 
@@ -493,6 +494,14 @@ def load_base(cls, simulation, structure, modelname='NewModel',
             strict mode when loading files
         model_rel_path : string
             relative path of model folder to simulation folder
+        load_only : list
+            list of package abbreviations or package names corresponding to
+            packages that flopy will load. default is None, which loads all
+            packages. the discretization packages will load regardless of this
+            setting. subpackages, like time series and observations, will also
+            load regardless of this setting.
+            example list: ['ic', 'maw', 'npf', 'oc', 'my_well_package_1']
+
         Returns
         -------
         model : MFModel
@@ -505,6 +514,10 @@ def load_base(cls, simulation, structure, modelname='NewModel',
                        version=version, exe_name=exe_name,
                        add_to_simulation=False, structure=structure,
                        model_rel_path=model_rel_path)
+
+        # build case consistent load_only dictionary for quick lookups
+        load_only = instance._load_only_dict(load_only)
+
         # load name file
         instance.name_file.load(strict)
 
@@ -527,9 +540,19 @@ def load_base(cls, simulation, structure, modelname='NewModel',
         sim_struct = mfstructure.MFStructure().sim_struct
         instance._ftype_num_dict = {}
         for ftype, fname, pname in packages_ordered:
+            ftype_orig = ftype
             ftype = ftype[0:-1].lower()
             if ftype in structure.package_struct_objs or ftype in \
               sim_struct.utl_struct_objs:
+                if load_only is not None and not \
+                        instance._in_pkg_list(priority_packages, ftype_orig,
+                                              pname) \
+                        and not instance._in_pkg_list(load_only, ftype_orig,
+                                                      pname):
+                    if simulation.simulation_data.verbosity_level.value >= \
+                            VerbosityLevel.normal.value:
+                        print('    skipping package {}...'.format(ftype))
+                    continue
                 if model_rel_path and model_rel_path != '.':
                     # strip off model relative path from the file path
                     filemgr = simulation.simulation_data.mfpath
@@ -596,7 +619,7 @@ def get_grid_type(self):
 
         Returns
         -------
-        grid type : DiscritizationType
+        grid type : DiscretizationType
         """
         package_recarray = self.name_file.packages
         structure = mfstructure.MFStructure()
diff --git a/flopy/mf6/modflow/mfgwf.py b/flopy/mf6/modflow/mfgwf.py
index e6c2211dc0..07766a3895 100644
--- a/flopy/mf6/modflow/mfgwf.py
+++ b/flopy/mf6/modflow/mfgwf.py
@@ -96,7 +96,9 @@ def __init__(self, simulation, modelname='model', model_nam_file=None,
     @classmethod
     def load(cls, simulation, structure, modelname='NewModel',
              model_nam_file='modflowtest.nam', version='mf6',
-             exe_name='mf6.exe', strict=True, model_rel_path='.'):
+             exe_name='mf6.exe', strict=True, model_rel_path='.',
+             load_only=None):
         return mfmodel.MFModel.load_base(simulation, structure, modelname,
                                          model_nam_file, 'gwf', version,
-                                         exe_name, strict, model_rel_path)
+                                         exe_name, strict, model_rel_path,
+                                         load_only)
diff --git a/flopy/mf6/modflow/mfsimulation.py b/flopy/mf6/modflow/mfsimulation.py
index 71d6b2c6a3..a66d82a54a 100644
--- a/flopy/mf6/modflow/mfsimulation.py
+++ b/flopy/mf6/modflow/mfsimulation.py
@@ -441,7 +441,7 @@ def model_names(self):
 
     @classmethod
     def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe',
-             sim_ws='.', strict=True, verbosity_level=1):
+             sim_ws='.', strict=True, verbosity_level=1, load_only=None):
         """
         Load an existing model.
 
@@ -465,6 +465,14 @@ def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe',
                     messages
                 2 : verbose mode with full error/warning/informational
                     messages.  this is ideal for debugging
+        load_only : list
+            list of package abbreviations or package names corresponding to
+            packages that flopy will load. default is None, which loads all
+            packages. the discretization packages will load regardless of this
+            setting. subpackages, like time series and observations, will also
+            load regardless of this setting.
+            example list: ['ic', 'maw', 'npf', 'oc', 'ims', 'gwf6-gwf6']
+
         Returns
         -------
         sim : MFSimulation object
@@ -480,6 +488,9 @@ def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe',
         if verbosity_level.value >= VerbosityLevel.normal.value:
             print('loading simulation...')
 
+        # build case consistent load_only dictionary for quick lookups
+        load_only = instance._load_only_dict(load_only)
+
         # load simulation name file
         if verbosity_level.value >= VerbosityLevel.normal.value:
             print('  loading simulation name file...')
@@ -520,7 +531,7 @@ def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe',
             instance._models[item[2]] = model_obj.load(
                 instance,
                 instance.structure.model_struct_objs[item[0].lower()], item[2],
-                name_file, version, exe_name, strict, path)
+                name_file, version, exe_name, strict, path, load_only)
 
         # load exchange packages and dependent packages
         try:
@@ -544,6 +555,14 @@ def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe',
                                       package='nam',
                                       message=message)
             for exgfile in exch_data:
+                if load_only is not None and not \
+                        instance._in_pkg_list(load_only, exgfile[0],
+                                              exgfile[2]):
+                    if instance.simulation_data.verbosity_level.value >= \
+                            VerbosityLevel.normal.value:
+                        print('    skipping package {}..'
+                              '.'.format(exgfile[0].lower()))
+                    continue
                 # get exchange type by removing numbers from exgtype
                 exchange_type = ''.join([char for char in exgfile[0] if
                                          not char.isdigit()]).upper()
@@ -604,6 +623,15 @@ def load(cls, sim_name='modflowsim', version='mf6', exe_name='mf6.exe',
                                   message=message)
         for solution_group in solution_group_dict.values():
             for solution_info in solution_group:
+                if load_only is not None and \
+                        not instance._in_pkg_list(load_only,
+                                                  solution_info[0],
+                                                  solution_info[2]):
+                    if instance.simulation_data.verbosity_level.value >= \
+                            VerbosityLevel.normal.value:
+                        print('    skipping package {}..'
+                              '.'.format(solution_info[0].lower()))
+                    continue
                 ims_file = mfims.ModflowIms(instance, filename=solution_info[1],
                                             pname=solution_info[2])
                 if verbosity_level.value >= VerbosityLevel.normal.value:
diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py
index 7450b7bb3a..4a7eb243e8 100644
--- a/flopy/modflow/mf.py
+++ b/flopy/modflow/mf.py
@@ -603,7 +603,7 @@ def load_results(self, **kwargs):
 
     @staticmethod
     def load(f, version='mf2005', exe_name='mf2005.exe', verbose=False,
-             model_ws='.', load_only=None, forgive=True, check=True):
+             model_ws='.', load_only=None, forgive=False, check=True):
         """
         Load an existing MODFLOW model.
 
diff --git a/flopy/modflow/mfbcf.py b/flopy/modflow/mfbcf.py
index cec2c20b76..b0660f6a62 100644
--- a/flopy/modflow/mfbcf.py
+++ b/flopy/modflow/mfbcf.py
@@ -321,6 +321,12 @@ def load(f, model, ext_unit_dict=None):
             istart = 0
             for k in range(nlay):
                 lcode = line[istart:istart + 2]
+                if lcode.strip() == '':
+                    # hit end of line before expected end of data
+                    # read next line
+                    line = f.readline()
+                    istart = 0
+                    lcode = line[istart:istart + 2]
                 lcode = lcode.replace(' ', '0')
                 t.append(lcode)
                 istart += 2
diff --git a/flopy/modflowlgr/mflgr.py b/flopy/modflowlgr/mflgr.py
index 022212c53f..fbbe4564b0 100644
--- a/flopy/modflowlgr/mflgr.py
+++ b/flopy/modflowlgr/mflgr.py
@@ -414,7 +414,7 @@ def change_model_ws(self, new_pth=None, reset_external=False):
 
     @staticmethod
     def load(f, version='mflgr', exe_name='mflgr.exe', verbose=False,
-             model_ws='.', load_only=None, forgive=True, check=True):
+             model_ws='.', load_only=None, forgive=False, check=True):
         """
         Load an existing model.
 
diff --git a/flopy/mt3d/mt.py b/flopy/mt3d/mt.py
index c19de1dcb8..0247759b64 100644
--- a/flopy/mt3d/mt.py
+++ b/flopy/mt3d/mt.py
@@ -554,6 +554,10 @@ def load(f, version='mt3dms', exe_name='mt3dms.exe', verbose=False,
             Filetype(s) to load (e.g. ['btn', 'adv'])
             (default is None, which means that all will be loaded)
 
+        forgive : bool, optional
+            Option to raise exceptions on package load failure, which can be
+            useful for debugging. Default False.
+
         modflowmodel : flopy.modflow.mf.Modflow
             This is a flopy Modflow model object upon which this Mt3dms
             model is based. (the default is None)
diff --git a/flopy/seawat/swt.py b/flopy/seawat/swt.py
index 34535c360c..d2b88dece9 100644
--- a/flopy/seawat/swt.py
+++ b/flopy/seawat/swt.py
@@ -428,11 +428,11 @@ def load(f, version='seawat', exe_name='swtv4', verbose=False,
                     verbose=verbose)
 
         mf = Modflow.load(f, version='mf2k', exe_name=None, verbose=verbose,
-                          model_ws=model_ws, load_only=load_only, forgive=True,
-                          check=False)
+                          model_ws=model_ws, load_only=load_only,
+                          forgive=False, check=False)
 
         mt = Mt3dms.load(f, version='mt3dms', exe_name=None, verbose=verbose,
-                         model_ws=model_ws, forgive=True)
+                         model_ws=model_ws, forgive=False)
 
         # set listing and global files using mf objects
         ms.lst = mf.lst