diff --git a/autotest/t027_test.py b/autotest/t027_test.py
index cd384044f0..33a125d162 100644
--- a/autotest/t027_test.py
+++ b/autotest/t027_test.py
@@ -378,8 +378,10 @@ def test_export():
 
     # netDF4 tests
     if netCDF4 is not None:
-        m.wel.export(os.path.join(cpth, "MNW2-Fig28_well.nc"))
-        m.mnw2.export(os.path.join(cpth, "MNW2-Fig28.nc"))
+        fcw = m.wel.export(os.path.join(cpth, "MNW2-Fig28_well.nc"))
+        fcw.write()
+        fcm = m.mnw2.export(os.path.join(cpth, "MNW2-Fig28.nc"))
+        fcm.write()
         fpth = os.path.join(cpth, "MNW2-Fig28.nc")
         nc = netCDF4.Dataset(fpth)
         assert np.array_equal(
@@ -509,7 +511,7 @@ def test_checks():
     test_make_well()
     test_blank_lines()
     test_make_package()
-    # test_export()
+    test_export()
     # test_checks()
     test_mnw1_load_write()
     test_mnw2_create_file()
diff --git a/autotest/t505_test.py b/autotest/t505_test.py
index 8eb9e14eca..9f548b6ef2 100644
--- a/autotest/t505_test.py
+++ b/autotest/t505_test.py
@@ -62,7 +62,7 @@
     os.makedirs(cpth)
 
 
-def np001():
+def test_np001():
     # init paths
     test_ex_name = "np001"
     model_name = "np001_mod"
@@ -405,7 +405,7 @@ def np001():
         budget_file = os.path.join(os.getcwd(), expected_cbc_file)
         budget_obj = bf.CellBudgetFile(budget_file, precision="double")
         budget_frf_valid = np.array(
-            budget_obj.get_data(text="FLOW-JA-FACE", full3D=True)
+            budget_obj.get_data(text="RIV", full3D=False)
         )
 
         # compare output to expected results
@@ -415,11 +415,8 @@ def np001():
         assert pymake.compare_heads(
             None, None, files1=head_file, files2=head_new, outfile=outfile
         )
-
-        budget_frf = sim.simulation_data.mfdata[
-            (model_name, "CBC", "FLOW-JA-FACE")
-        ]
-        assert array_util.array_comp(budget_frf_valid, budget_frf)
+        budget_frf = sim.simulation_data.mfdata[(model_name, "CBC", "RIV")]
+        assert array_util.riv_array_comp(budget_frf_valid, budget_frf)
 
         # clean up
         sim.delete_output_files()
@@ -443,7 +440,7 @@ def np001():
         budget_file = os.path.join(os.getcwd(), expected_cbc_file)
         budget_obj = bf.CellBudgetFile(budget_file, precision="double")
         budget_frf_valid = np.array(
-            budget_obj.get_data(text="FLOW-JA-FACE", full3D=True)
+            budget_obj.get_data(text="RIV", full3D=False)
         )
 
         # compare output to expected results
@@ -454,14 +451,21 @@ def np001():
             None, None, files1=head_file, files2=head_new, outfile=outfile
         )
 
-        budget_frf = sim.simulation_data.mfdata[
-            (model_name, "CBC", "FLOW-JA-FACE")
-        ]
-        assert array_util.array_comp(budget_frf_valid, budget_frf)
+        budget_frf = sim.simulation_data.mfdata[(model_name, "CBC", "RIV")]
+        assert array_util.riv_array_comp(budget_frf_valid, budget_frf)
 
         # clean up
         sim.delete_output_files()
 
+    # test rename all packages
+    rename_folder = os.path.join(run_folder, "rename")
+    sim.rename_all_packages("file_rename")
+    sim.set_sim_path(rename_folder)
+    sim.write_simulation()
+    if run:
+        sim.run_simulation()
+        sim.delete_output_files()
+
     try:
         error_occurred = False
         well_spd = {
@@ -543,6 +547,7 @@ def np001():
     )
     wel_package.write()
     mpath = sim.simulation_data.mfpath.get_model_path(model.name)
+    spath = sim.simulation_data.mfpath.get_sim_path()
     found_cellid = False
     with open(os.path.join(mpath, "np001_mod.wel"), "r") as fd:
         for line in fd:
@@ -560,6 +565,7 @@ def np001():
     well_spd = {0: [(-1, -1, -1, -2000.0), (0, 0, 7, -2.0)], 1: []}
     wel_package = ModflowGwfwel(
         model,
+        filename="file_rename.wel",
         print_input=True,
         print_flows=True,
         save_flows=True,
@@ -570,7 +576,7 @@ def np001():
     found_begin = False
     found_end = False
     text_between_begin_and_end = False
-    with open(os.path.join(mpath, "np001_mod.wel"), "r") as fd:
+    with open(os.path.join(mpath, "file_rename.wel"), "r") as fd:
         for line in fd:
             if line.strip().lower() == "begin period  2":
                 found_begin = True
@@ -587,7 +593,7 @@ def np001():
         test_ex_name,
         "mf6",
         exe_name,
-        run_folder,
+        spath,
         write_headers=False,
     )
     wel = test_sim.get_model().wel
@@ -611,7 +617,7 @@ def np001():
     return
 
 
-def np002():
+def test_np002():
     # init paths
     test_ex_name = "np002"
     model_name = "np002_mod"
@@ -780,13 +786,6 @@ def np002():
         npf_package = model_.get_package("npf")
         k = npf_package.k.array
 
-        # get expected results
-        budget_file = os.path.join(os.getcwd(), expected_cbc_file)
-        budget_obj = bf.CellBudgetFile(budget_file, precision="double")
-        budget_frf_valid = np.array(
-            budget_obj.get_data(text="FLOW JA FACE    ", full3D=True)
-        )
-
         # compare output to expected results
         head_file = os.path.join(os.getcwd(), expected_head_file)
         head_new = os.path.join(run_folder, "np002_mod.hds")
@@ -795,12 +794,6 @@ def np002():
             None, None, files1=head_file, files2=head_new, outfile=outfile
         )
 
-        array_util = PyListUtil()
-        budget_frf = sim.simulation_data.mfdata[
-            (model_name, "CBC", "FLOW-JA-FACE")
-        ]
-        assert array_util.array_comp(budget_frf_valid, budget_frf)
-
         # verify external text file was written correctly
         ext_file_path = os.path.join(run_folder, "initial_heads.txt")
         fd = open(ext_file_path, "r")
@@ -2505,28 +2498,6 @@ def test006_2models_gnc():
         stress_period_data=stress_period_data,
     )
 
-    gncrecarray = testutils.read_gncrecarray(os.path.join(pth, "gnc.txt"))
-    # test gnc delete
-    new_gncrecarray = gncrecarray[10:]
-    gnc_package = ModflowGwfgnc(
-        sim,
-        print_input=True,
-        print_flows=True,
-        numgnc=26,
-        numalphaj=1,
-        gncdata=new_gncrecarray,
-    )
-    sim.remove_package(gnc_package.package_type)
-
-    gnc_package = ModflowGwfgnc(
-        sim,
-        print_input=True,
-        print_flows=True,
-        numgnc=36,
-        numalphaj=1,
-        gncdata=gncrecarray,
-    )
-
     exgrecarray = testutils.read_exchangedata(os.path.join(pth, "exg.txt"))
 
     # build obs dictionary
@@ -2569,6 +2540,28 @@ def test006_2models_gnc():
         observations=gwf_obs,
     )
 
+    gncrecarray = testutils.read_gncrecarray(os.path.join(pth, "gnc.txt"))
+    # test gnc delete
+    new_gncrecarray = gncrecarray[10:]
+    gnc_package = ModflowGwfgnc(
+        sim,
+        print_input=True,
+        print_flows=True,
+        numgnc=26,
+        numalphaj=1,
+        gncdata=new_gncrecarray,
+    )
+    sim.remove_package(gnc_package.package_type)
+
+    gnc_package = ModflowGwfgnc(
+        sim,
+        print_input=True,
+        print_flows=True,
+        numgnc=36,
+        numalphaj=1,
+        gncdata=gncrecarray,
+    )
+
     # change folder to save simulation
     sim.simulation_data.mfpath.set_sim_path(run_folder)
 
@@ -2616,6 +2609,15 @@ def test006_2models_gnc():
         # clean up
         sim.delete_output_files()
 
+    # test rename all packages
+    rename_folder = os.path.join(run_folder, "rename")
+    sim.rename_all_packages("file_rename")
+    sim.set_sim_path(rename_folder)
+    sim.write_simulation()
+    if run:
+        sim.run_simulation()
+        sim.delete_output_files()
+
     return
 
 
@@ -3228,8 +3230,8 @@ def test_transport():
 
 
 if __name__ == "__main__":
-    np001()
-    np002()
+    test_np001()
+    test_np002()
     test004_bcfss()
     test005_advgw_tidal()
     test006_2models_gnc()
diff --git a/examples/data/mf6/create_tests/np002/expected_output/np002_mod.hds b/examples/data/mf6/create_tests/np002/expected_output/np002_mod.hds
index 3c4f05d62c..e874506985 100644
Binary files a/examples/data/mf6/create_tests/np002/expected_output/np002_mod.hds and b/examples/data/mf6/create_tests/np002/expected_output/np002_mod.hds differ
diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py
index a77cb3cf1e..1875d0c3f1 100644
--- a/flopy/mf6/data/mfdatalist.py
+++ b/flopy/mf6/data/mfdatalist.py
@@ -797,6 +797,8 @@ def _get_file_entry(
             self._crnt_line_num = 1
             for mflist_line in range(0, data_lines):
                 text_line = []
+
+                # data
                 index = 0
                 self._get_file_entry_record(
                     data_complete,
@@ -808,12 +810,12 @@ def _get_file_entry(
                     indent,
                 )
 
-                # include comments
+                # comments
                 if (
-                    mflist_line in storage.comments
-                    and storage.comments[mflist_line].text
+                    mflist_line + 1 in storage.comments
+                    and storage.comments[mflist_line + 1].text
                 ):
-                    text_line.append(storage.comments[mflist_line].text)
+                    text_line.append(storage.comments[mflist_line + 1].text)
 
                 file_entry.append(f"{indent}{indent.join(text_line)}\n")
                 self._crnt_line_num += 1
diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py
index 233ee92bfe..d3992c2057 100644
--- a/flopy/mf6/data/mfdatastorage.py
+++ b/flopy/mf6/data/mfdatastorage.py
@@ -907,10 +907,8 @@ def _set_list(
             elif "data" in data:
                 data = data["data"]
         if isinstance(data, list):
-            if (
-                len(data) > 0
-                and not isinstance(data[0], tuple)
-                and not isinstance(data[0], list)
+            if len(data) > 0 and (
+                not PyListUtil.is_iterable(data[0]) or isinstance(data[0], str)
             ):
                 # single line of data needs to be encapsulated in a tuple
                 data = [tuple(data)]
@@ -1738,6 +1736,10 @@ def external_to_external(
         )
 
     def external_to_internal(self, layer, store_internal=False):
+        # reset comments
+        self.pre_data_comments = None
+        self.comments = {}
+
         if layer is None:
             layer = 0
         # load data from external file
diff --git a/flopy/mf6/data/mffileaccess.py b/flopy/mf6/data/mffileaccess.py
index e7c1cb7626..e3fa5f9da0 100644
--- a/flopy/mf6/data/mffileaccess.py
+++ b/flopy/mf6/data/mffileaccess.py
@@ -87,8 +87,6 @@ def _read_pre_data_comments(
                     arr_line, self._path, self._simulation_data, line_num
                 )
 
-            storage.add_data_line_comment(arr_line, line_num)
-
             line = file_handle.readline()
             arr_line = PyListUtil.split_data_line(line)
         return line
@@ -1175,7 +1173,15 @@ def read_list_data_from_file(
         # read any pre-data commented lines
         while current_line and MFComment.is_comment(arr_line, True):
             arr_line.insert(0, "\n")
-            storage.add_data_line_comment(arr_line, line_num)
+            if storage.pre_data_comments is None:
+                storage.pre_data_comments = MFComment(
+                    " ".join(arr_line),
+                    self._path,
+                    self._simulation_data,
+                    line_num,
+                )
+            else:
+                storage.pre_data_comments.add_text(" ".join(arr_line))
             PyListUtil.reset_delimiter_used()
             current_line = file_handle.readline()
             arr_line = PyListUtil.split_data_line(current_line)
diff --git a/flopy/mf6/mfbase.py b/flopy/mf6/mfbase.py
index 9f87820f09..2364925053 100644
--- a/flopy/mf6/mfbase.py
+++ b/flopy/mf6/mfbase.py
@@ -647,6 +647,37 @@ def _remove_package(self, package):
         for key in items_to_remove:
             del self.simulation_data.mfdata[key]
 
+    def _rename_package(self, package, new_name):
+        # fix package_name_dict key
+        if (
+            package.package_name is not None
+            and package.package_name.lower() in self.package_name_dict
+        ):
+            del self.package_name_dict[package.package_name.lower()]
+        self.package_name_dict[new_name.lower()] = package
+        # fix package_key_dict key
+        new_package_path = package.path[:-1] + (new_name,)
+        del self.package_key_dict[package.path[-1].lower()]
+        self.package_key_dict[new_package_path.lower()] = package
+        # get keys to fix in main dictionary
+        main_dict = self.simulation_data.mfdata
+        items_to_fix = []
+        for key in main_dict:
+            is_subkey = True
+            for pitem, ditem in zip(package.path, key):
+                if pitem != ditem:
+                    is_subkey = False
+                    break
+            if is_subkey:
+                items_to_fix.append(key)
+
+        # fix keys in main dictionary
+        for key in items_to_fix:
+            new_key = (
+                package.path[:-1] + (new_name,) + key[len(package.path) - 1 :]
+            )
+            main_dict[new_key] = main_dict.pop(key)
+
     def get_package(self, name=None):
         """
         Finds a package by package name, package key, package type, or partial
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
index f0006fe560..17a14ec0ea 100644
--- a/flopy/mf6/mfmodel.py
+++ b/flopy/mf6/mfmodel.py
@@ -1071,14 +1071,15 @@ def remove_package(self, package_name):
             packages = [package_name]
         else:
             packages = self.get_package(package_name)
-            if not isinstance(packages, list):
+            if not isinstance(packages, list) and packages is not None:
                 packages = [packages]
+        if packages is None:
+            return
         for package in packages:
             if package.model_or_sim.name != self.name:
                 except_text = (
-                    "Package can not be removed from model {} "
-                    "since it is "
-                    "not part of "
+                    "Package can not be removed from model "
+                    "{self.model_name} since it is not part of it."
                 )
                 raise mfstructure.FlopyException(except_text)
 
@@ -1091,7 +1092,7 @@ def remove_package(self, package_name):
                 message = (
                     "Error occurred while reading package names "
                     "from name file in model "
-                    '"{}".'.format(self.name)
+                    f'"{self.name}"'
                 )
                 raise MFDataException(
                     mfdata_except=mfde,
@@ -1129,8 +1130,8 @@ def remove_package(self, package_name):
             except MFDataException as mfde:
                 message = (
                     "Error occurred while setting package names "
-                    'from name file in model "{}".  Package name '
-                    "data:\n{}".format(self.name, new_rec_array)
+                    f'from name file in model "{self.name}".  Package name '
+                    f"data:\n{new_rec_array}"
                 )
                 raise MFDataException(
                     mfdata_except=mfde,
@@ -1151,6 +1152,76 @@ def remove_package(self, package_name):
             for child_package in child_package_list:
                 self._remove_package_from_dictionaries(child_package)
 
+    def update_package_filename(self, package, new_name):
+        """
+        Updates the filename for a package.  For internal flopy use only.
+
+        Parameters
+        ----------
+        package : MFPackage
+            Package object
+        new_name : str
+            New package name
+        """
+        try:
+            # get namefile package data
+            package_data = self.name_file.packages.get_data()
+        except MFDataException as mfde:
+            message = (
+                "Error occurred while updating package names "
+                "from name file in model "
+                f'"{self.name}".'
+            )
+            raise MFDataException(
+                mfdata_except=mfde,
+                model=self.model_name,
+                package=self.name_file._get_pname(),
+                message=message,
+            )
+        try:
+            # update namefile package data with new name
+            new_rec_array = None
+            for item in package_data:
+                base, leaf = os.path.split(item[1])
+                if leaf == package.filename:
+                    item[1] = os.path.join(base, new_name)
+
+                if new_rec_array is None:
+                    new_rec_array = np.rec.array(
+                        [item.tolist()], package_data.dtype
+                    )
+                else:
+                    new_rec_array = np.hstack((item, new_rec_array))
+        except:
+            type_, value_, traceback_ = sys.exc_info()
+            raise MFDataException(
+                self.structure.get_model(),
+                self.structure.get_package(),
+                self._path,
+                "updating package filename",
+                self.structure.name,
+                inspect.stack()[0][3],
+                type_,
+                value_,
+                traceback_,
+                None,
+                self._simulation_data.debug,
+            )
+        try:
+            self.name_file.packages.set_data(new_rec_array)
+        except MFDataException as mfde:
+            message = (
+                "Error occurred while updating package names "
+                f'from name file in model "{self.name}".  Package name '
+                f"data:\n{new_rec_array}"
+            )
+            raise MFDataException(
+                mfdata_except=mfde,
+                model=self.model_name,
+                package=self.name_file._get_pname(),
+                message=message,
+            )
+
     def rename_all_packages(self, name):
         """Renames all package files in the model.
 
@@ -1161,8 +1232,11 @@ def rename_all_packages(self, name):
                 <name>.<package ext>.
 
         """
+        nam_filename = f"{name}.nam"
+        self.simulation.rename_model_namefile(self, nam_filename)
+        self.name_file.filename = nam_filename
+        self.model_nam_file = nam_filename
         package_type_count = {}
-        self.name_file.filename = f"{name}.nam"
         for package in self.packagelist:
             if package.package_type not in package_type_count:
                 package.filename = f"{name}.{package.package_type}"
diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py
index 6db2443e95..be54b007f8 100644
--- a/flopy/mf6/mfpackage.py
+++ b/flopy/mf6/mfpackage.py
@@ -1486,6 +1486,10 @@ def __init__(
         else:
             self.model_name = None
 
+        # a package must have a dfn_file_name
+        if not hasattr(self, "dfn_file_name"):
+            self.dfn_file_name = ""
+
         if model_or_sim.type != "Model" and model_or_sim.type != "Simulation":
             message = (
                 "Invalid model_or_sim parameter. Expecting either a "
@@ -1638,8 +1642,11 @@ def filename(self, fname):
             except Exception:
                 print(
                     "WARNING: Unable to update file name for parent"
-                    "package of {}.".format(self.name)
+                    f"package of {self.package_name}."
                 )
+        if self.model_or_sim is not None and fname is not None:
+            if self._package_type != "nam":
+                self.model_or_sim.update_package_filename(self, fname)
         self._filename = fname
 
     @property
@@ -2385,6 +2392,61 @@ def create_package_dimensions(self):
                             ),
                         ]
                         break
+            elif (
+                self.dfn_file_name[4:7] == "gnc"
+                and self.model_or_sim.type == "Simulation"
+            ):
+                # get exchange file name associated with gnc package
+                exg_file_name = None
+                for exg in self.model_or_sim.exchange_files:
+                    gnc_data = exg.gnc_filerecord.get_data()
+                    if (
+                        gnc_data is not None
+                        and gnc_data[0][0].lower() == self.filename.lower()
+                    ):
+                        exg_file_name = exg.filename
+                        self.parent = exg
+                if exg_file_name is None:
+                    raise Exception(
+                        "Can not create a simulation-level "
+                        "gnc file without a corresponding "
+                        "exchange file. Exchange file must be "
+                        "created first."
+                    )
+                # get models associated with exchange file from sim nam file
+                try:
+                    exchange_recarray_data = (
+                        self.model_or_sim.name_file.exchanges.get_data()
+                    )
+                except MFDataException as mfde:
+                    message = (
+                        "An error occurred while retrieving exchange "
+                        "data from the simulation name file.  The error "
+                        "occurred while processing gnc file "
+                        f'"{self.filename}".'
+                    )
+                    raise MFDataException(
+                        mfdata_except=mfde,
+                        package=self._get_pname(),
+                        message=message,
+                    )
+                assert exchange_recarray_data is not None
+                model_1 = None
+                model_2 = None
+                for exchange in exchange_recarray_data:
+                    if exchange[1] == exg_file_name:
+                        model_1 = exchange[2]
+                        model_2 = exchange[3]
+
+                # assign models to gnc package
+                model_dims = [
+                    modeldimensions.ModelDimensions(
+                        model_1, self._simulation_data
+                    ),
+                    modeldimensions.ModelDimensions(
+                        model_2, self._simulation_data
+                    ),
+                ]
             elif self.parent_file is not None:
                 model_dims = []
                 for md in self.parent_file.dimensions.model_dim:
@@ -2650,8 +2712,14 @@ def _update_filename(self, old_fname, new_fname):
         if file_record is not None:
             file_record_data = file_record[0]
             for item in file_record_data:
-                if item.lower() == old_fname.lower():
-                    new_file_record_data.append((new_fname,))
+                base, fname = os.path.split(item)
+                if fname.lower() == old_fname.lower():
+                    if base:
+                        new_file_record_data.append(
+                            (os.path.join(base, new_fname),)
+                        )
+                    else:
+                        new_file_record_data.append((new_fname,))
                 else:
                     new_file_record_data.append((item,))
         else:
diff --git a/flopy/mf6/modflow/mfsimulation.py b/flopy/mf6/modflow/mfsimulation.py
index c747e72baf..b19fa03bcd 100644
--- a/flopy/mf6/modflow/mfsimulation.py
+++ b/flopy/mf6/modflow/mfsimulation.py
@@ -1,8 +1,8 @@
 import errno
 import sys
 import inspect
-import collections
 import os.path
+import numpy as np
 from ...mbase import run_model
 from ..mfbase import (
     PackageContainer,
@@ -22,13 +22,13 @@
 from ..data.mfdatautil import MFComment
 
 
-class SimulationDict(collections.OrderedDict):
+class SimulationDict(dict):
     """
     Class containing custom dictionary for MODFLOW simulations.  Dictionary
     contains model data.  Dictionary keys are "paths" to the data that include
     the model and package containing the data.
 
-    Behaves as an OrderedDict with some additional features described below.
+    Behaves as an dict with some additional features described below.
 
     Parameters
     ----------
@@ -38,11 +38,12 @@ class SimulationDict(collections.OrderedDict):
     """
 
     def __init__(self, path=None):
-        collections.OrderedDict.__init__(self)
+        dict.__init__(self)
         self._path = path
 
     def __getitem__(self, key):
-        """Define the __getitem__ magic method.
+        """
+        Define the __getitem__ magic method.
 
         Parameters
         ----------
@@ -65,12 +66,13 @@ def __getitem__(self, key):
             return val.data
 
         if key in self:
-            val = collections.OrderedDict.__getitem__(self, key)
+            val = dict.__getitem__(self, key)
             return val
         return AttributeError(key)
 
     def __setitem__(self, key, val):
-        """Define the __setitem__ magic method.
+        """
+        Define the __setitem__ magic method.
 
         Parameters
         ----------
@@ -80,10 +82,11 @@ def __setitem__(self, key, val):
             MFData to store in dictionary
 
         """
-        collections.OrderedDict.__setitem__(self, key, val)
+        dict.__setitem__(self, key, val)
 
     def find_in_path(self, key_path, key_leaf):
-        """Attempt to find key_leaf in a partial key path key_path.
+        """
+        Attempt to find key_leaf in a partial key path key_path.
 
         Parameters
         ----------
@@ -117,7 +120,8 @@ def find_in_path(self, key_path, key_leaf):
         return None, None
 
     def output_keys(self, print_keys=True):
-        """Return a list of output data keys supported by the dictionary.
+        """
+        Return a list of output data keys supported by the dictionary.
 
         Parameters
         ----------
@@ -136,7 +140,8 @@ def output_keys(self, print_keys=True):
         return [key for key in x.dataDict]
 
     def input_keys(self):
-        """Return a list of input data keys.
+        """
+        Return a list of input data keys.
 
         Returns
         -------
@@ -148,7 +153,8 @@ def input_keys(self):
             print(key)
 
     def observation_keys(self):
-        """Return a list of observation keys.
+        """
+        Return a list of observation keys.
 
         Returns
         -------
@@ -159,7 +165,8 @@ def observation_keys(self):
         mfobservation.MFObservationRequester.getkeys(self, self._path)
 
     def keys(self):
-        """Return a list of all keys.
+        """
+        Return a list of all keys.
 
         Returns
         -------
@@ -276,7 +283,8 @@ def max_columns_of_data(self, val):
             self.max_columns_user_set = True
 
     def set_sci_note_upper_thres(self, value):
-        """Sets threshold number where any number larger than threshold
+        """
+        Sets threshold number where any number larger than threshold
         is represented in scientific notation.
 
         Parameters
@@ -289,7 +297,8 @@ def set_sci_note_upper_thres(self, value):
         self._update_str_format()
 
     def set_sci_note_lower_thres(self, value):
-        """Sets threshold number where any number smaller than threshold
+        """
+        Sets threshold number where any number smaller than threshold
         is represented in scientific notation.
 
         Parameters
@@ -302,7 +311,8 @@ def set_sci_note_lower_thres(self, value):
         self._update_str_format()
 
     def _update_str_format(self):
-        """Update floating point formatting strings."""
+        """
+        Update floating point formatting strings."""
         self.reg_format_str = f"{{:.{self.float_precision}E}}"
         self.sci_format_str = (
             f"{{:{self.float_characters}.{self.float_precision}f}}"
@@ -447,7 +457,8 @@ def __init__(
         self.valid = False
 
     def __getattr__(self, item):
-        """Override __getattr__ to allow retrieving models.
+        """
+        Override __getattr__ to allow retrieving models.
 
         __getattr__ is used to allow for getting models and packages as if
         they are attributes
@@ -489,7 +500,8 @@ def __getattr__(self, item):
             raise AttributeError(item)
 
     def __repr__(self):
-        """Override __repr__ to print custom string.
+        """
+        Override __repr__ to print custom string.
 
         Returns
         --------
@@ -500,7 +512,8 @@ def __repr__(self):
         return self._get_data_str(True)
 
     def __str__(self):
-        """Override __str__ to print custom string.
+        """
+        Override __str__ to print custom string.
 
         Returns
         --------
@@ -555,7 +568,8 @@ def _get_data_str(self, formal):
 
     @property
     def model_names(self):
-        """Return a list of model names associated with this simulation.
+        """
+        Return a list of model names associated with this simulation.
 
         Returns
         --------
@@ -564,6 +578,18 @@ def model_names(self):
         """
         return self._models.keys()
 
+    @property
+    def exchange_files(self):
+        """
+        Return list of exchange files associated with this simulation.
+
+        Returns
+        --------
+            list: list of exchange names
+
+        """
+        return self._exchange_files.values()
+
     @classmethod
     def load(
         cls,
@@ -577,7 +603,8 @@ def load(
         verify_data=False,
         write_headers=True,
     ):
-        """Load an existing model.
+        """
+        Load an existing model.
 
         Parameters
         ----------
@@ -910,7 +937,8 @@ def load_package(
         dict_package_name=None,
         parent_package=None,
     ):
-        """Load a package from a file.
+        """
+        Load a package from a file.
 
         Parameters
         ----------
@@ -1003,7 +1031,8 @@ def load_package(
             return package
 
     def register_ims_package(self, ims_file, model_list):
-        """Register an ims package with the simulation.
+        """
+        Register an ims package with the simulation.
 
         Parameters
             ims_file : MFPackage
@@ -1054,6 +1083,7 @@ def register_ims_package(self, ims_file, model_list):
                     )
                 self._remove_package(self._ims_files[file.filename])
                 del self._ims_files[file.filename]
+                break
         # register ims package
         if not in_simulation:
             self._add_package(ims_file, self._get_package_path(ims_file))
@@ -1123,20 +1153,126 @@ def register_ims_package(self, ims_file, model_list):
     @staticmethod
     def _rename_package_group(group_dict, name):
         package_type_count = {}
+        # first build an array to avoid key modification errors
+        package_array = []
         for package in group_dict.values():
+            package_array.append(package)
+        # update package file names and count
+        for package in package_array:
             if package.package_type not in package_type_count:
                 package.filename = f"{name}.{package.package_type}"
                 package_type_count[package.package_type] = 1
             else:
                 package_type_count[package.package_type] += 1
-                package.filename = "{}_{}.{}".format(
-                    name,
-                    package_type_count[package.package.package_type],
-                    package.package_type,
-                )
+                ptc = package_type_count[package.package_type]
+                package.filename = f"{name}_{ptc}.{package.package_type}"
+
+    def _rename_exchange_file(self, package, new_filename):
+        self._exchange_files[package.filename] = package
+        try:
+            exchange_recarray_data = self.name_file.exchanges.get_data()
+        except MFDataException as mfde:
+            message = (
+                "An error occurred while retrieving exchange "
+                "data from the simulation name file.  The error "
+                "occurred while registering exchange file "
+                f'"{package.filename}".'
+            )
+            raise MFDataException(
+                mfdata_except=mfde,
+                package=package._get_pname(),
+                message=message,
+            )
+        if exchange_recarray_data is not None:
+            for index, exchange in zip(
+                range(0, len(exchange_recarray_data)),
+                exchange_recarray_data,
+            ):
+                if exchange[1] == package.filename:
+                    # update existing exchange
+                    exchange_recarray_data[index][1] = new_filename
+                    ex_recarray = self.name_file.exchanges
+                    try:
+                        ex_recarray.set_data(exchange_recarray_data)
+                    except MFDataException as mfde:
+                        message = (
+                            "An error occurred while setting "
+                            "exchange data in the simulation name "
+                            "file.  The error occurred while "
+                            "registering the following "
+                            "values (exgtype, filename, "
+                            f'exgmnamea, exgmnameb): "{package.exgtype} '
+                            f"{package.filename} {package.exgmnamea}"
+                            f'{package.exgmnameb}".'
+                        )
+                        raise MFDataException(
+                            mfdata_except=mfde,
+                            package=package._get_pname(),
+                            message=message,
+                        )
+                    return
+
+    def _set_timing_block(self, file_name):
+        struct_root = mfstructure.MFStructure()
+        tdis_pkg = "tdis{}".format(struct_root.get_version_string())
+        tdis_attr = getattr(self.name_file, tdis_pkg)
+        try:
+            tdis_attr.set_data(file_name)
+        except MFDataException as mfde:
+            message = (
+                "An error occurred while setting the tdis package "
+                f'file name "{file_name}".  The error occurred while '
+                "registering the tdis package with the "
+                "simulation"
+            )
+            raise MFDataException(
+                mfdata_except=mfde,
+                package=file_name,
+                message=message,
+            )
+
+    def update_package_filename(self, package, new_name):
+        """
+        Updates internal arrays to be consistent with a new file name.
+        This is for internal flopy library use only.
+
+        Parameters
+        ----------
+            package: MFPackage
+                Package with new name
+            new_name: str
+                Package's new name
+
+        """
+        if (
+            self._tdis_file is not None
+            and package.filename == self._tdis_file.filename
+        ):
+            self._set_timing_block(new_name)
+        if package.filename in self._exchange_files:
+            self._exchange_files[new_name] = self._exchange_files.pop(
+                package.filename
+            )
+            self._rename_exchange_file(package, new_name)
+        if package.filename in self._ims_files:
+            self._ims_files[new_name] = self._ims_files.pop(package.filename)
+            self._update_ims_solution_group(package.filename, new_name)
+        if package.filename in self._ghost_node_files:
+            self._ghost_node_files[new_name] = self._ghost_node_files.pop(
+                package.filename
+            )
+        if package.filename in self._mover_files:
+            self._mover_files[new_name] = self._mover_files.pop(
+                package.filename
+            )
+        if package.filename in self._other_files:
+            self._other_files[new_name] = self._other_files.pop(
+                package.filename
+            )
 
     def rename_all_packages(self, name):
-        """Rename all packages with name as prefix.
+        """
+        Rename all packages with name as prefix.
 
         Parameters
         ----------
@@ -1206,7 +1342,8 @@ def set_all_data_internal(self, check_data=True):
     def write_simulation(
         self, ext_file_action=ExtFileAction.copy_relative_paths, silent=False
     ):
-        """Write the simulation to files.
+        """
+        Write the simulation to files.
 
         Parameters
             ext_file_action : ExtFileAction
@@ -1396,7 +1533,8 @@ def run_simulation(
         use_async=False,
         cargs=None,
     ):
-        """Run the simulation.
+        """
+        Run the simulation.
 
         Parameters
         ----------
@@ -1455,7 +1593,8 @@ def delete_output_files(self):
                 os.remove(path)
 
     def remove_package(self, package_name):
-        """Removes package from the simulation.  `package_name` can be the
+        """
+        Removes package from the simulation.  `package_name` can be the
         package's name, type, or package object to be removed from the model.
 
         Parameters
@@ -1480,7 +1619,7 @@ def remove_package(self, package_name):
                 del self._exchange_files[package.filename]
             if package.filename in self._ims_files:
                 del self._ims_files[package.filename]
-                self._remove_ims_soultion_group(package.filename)
+                self._update_ims_solution_group(package.filename)
             if package.filename in self._ghost_node_files:
                 del self._ghost_node_files[package.filename]
             if package.filename in self._mover_files:
@@ -1492,7 +1631,8 @@ def remove_package(self, package_name):
 
     @property
     def model_dict(self):
-        """Return a dictionary of models associated with this simulation.
+        """
+        Return a dictionary of models associated with this simulation.
 
         Returns
         --------
@@ -1503,7 +1643,8 @@ def model_dict(self):
         return self._models.copy()
 
     def get_model(self, model_name=None):
-        """Returns the models in the simulation with a given model name, name
+        """
+        Returns the models in the simulation with a given model name, name
         file name, or model type.
 
         Parameters
@@ -1532,7 +1673,8 @@ def get_model(self, model_name=None):
         return None
 
     def get_exchange_file(self, filename):
-        """Get a specified exchange file.
+        """
+        Get a specified exchange file.
 
         Parameters
         ----------
@@ -1551,7 +1693,8 @@ def get_exchange_file(self, filename):
             raise FlopyException(excpt_str)
 
     def get_mvr_file(self, filename):
-        """Get a specified mover file.
+        """
+        Get a specified mover file.
 
         Parameters
         ----------
@@ -1570,7 +1713,8 @@ def get_mvr_file(self, filename):
             raise FlopyException(excpt_str)
 
     def get_gnc_file(self, filename):
-        """Get a specified gnc file.
+        """
+        Get a specified gnc file.
 
         Parameters
         ----------
@@ -1588,8 +1732,51 @@ def get_gnc_file(self, filename):
             excpt_str = f'GNC file "{filename}" can not be found.'
             raise FlopyException(excpt_str)
 
+    def remove_exchange_file(self, package):
+        """
+        Removes the exchange file "package". This is for internal flopy
+        library use only.
+
+        Parameters
+        ----------
+            package: MFPackage
+                Exchange package to be removed
+
+        """
+        self._exchange_files[package.filename] = package
+        try:
+            exchange_recarray_data = self.name_file.exchanges.get_data()
+        except MFDataException as mfde:
+            message = (
+                "An error occurred while retrieving exchange "
+                "data from the simulation name file.  The error "
+                "occurred while registering exchange file "
+                f'"{package.filename}".'
+            )
+            raise MFDataException(
+                mfdata_except=mfde,
+                package=package._get_pname(),
+                message=message,
+            )
+        remove_indices = []
+        if exchange_recarray_data is not None:
+            for index, exchange in zip(
+                range(0, len(exchange_recarray_data)),
+                exchange_recarray_data,
+            ):
+                if (
+                    package.filename is not None
+                    and exchange[1] == package.filename
+                ):
+                    remove_indices.append(index)
+        if len(remove_indices) > 0:
+            self.name_file.exchanges.set_data(
+                np.delete(exchange_recarray_data, remove_indices)
+            )
+
     def register_exchange_file(self, package):
-        """Register an exchange package file with the simulation.  This is a
+        """
+        Register an exchange package file with the simulation.  This is a
         call-back method made from the package and should not be called
         directly.
 
@@ -1619,7 +1806,7 @@ def register_exchange_file(self, package):
                     "An error occurred while retrieving exchange "
                     "data from the simulation name file.  The error "
                     "occurred while registering exchange file "
-                    '"{}".'.format(package.filename)
+                    f'"{package.filename}".'
                 )
                 raise MFDataException(
                     mfdata_except=mfde,
@@ -1646,13 +1833,9 @@ def register_exchange_file(self, package):
                                 "file.  The error occurred while "
                                 "registering the following "
                                 "values (exgtype, filename, "
-                                'exgmnamea, exgmnameb): "{} {} {}'
-                                '{}".'.format(
-                                    exgtype,
-                                    package.filename,
-                                    exgmnamea,
-                                    exgmnameb,
-                                )
+                                f'exgmnamea, exgmnameb): "{exgtype} '
+                                f"{package.filename} {exgmnamea}"
+                                f'{exgmnameb}".'
                             )
                             raise MFDataException(
                                 mfdata_except=mfde,
@@ -1670,10 +1853,8 @@ def register_exchange_file(self, package):
                     "An error occurred while setting exchange data "
                     "in the simulation name file.  The error occurred "
                     "while registering the following values (exgtype, "
-                    'filename, exgmnamea, exgmnameb): "{} {} {}'
-                    '{}".'.format(
-                        exgtype, package.filename, exgmnamea, exgmnameb
-                    )
+                    f'filename, exgmnamea, exgmnameb): "{exgtype} '
+                    f'{package.filename} {exgmnamea} {exgmnameb}".'
                 )
                 raise MFDataException(
                     mfdata_except=mfde,
@@ -1687,6 +1868,75 @@ def register_exchange_file(self, package):
                 # resolve exchange package dimensions object
                 package.dimensions = package.create_package_dimensions()
 
+    def _remove_package_by_type(self, package):
+        pname = None
+        if package.package_name is not None:
+            pname = package.package_name.lower()
+        if (
+            package.package_type.lower() == "tdis"
+            and self._tdis_file is not None
+            and self._tdis_file in self._packagelist
+        ):
+            # tdis package already exists. there can be only one tdis
+            # package.  remove existing tdis package
+            if (
+                self.simulation_data.verbosity_level.value
+                >= VerbosityLevel.normal.value
+            ):
+                print(
+                    "WARNING: tdis package already exists. Replacing "
+                    "existing tdis package."
+                )
+            self._remove_package(self._tdis_file)
+        elif (
+            package.package_type.lower() == "gnc"
+            and package.filename in self._ghost_node_files
+            and self._ghost_node_files[package.filename] in self._packagelist
+        ):
+            # gnc package with same file name already exists.  remove old
+            # gnc package
+            if (
+                self.simulation_data.verbosity_level.value
+                >= VerbosityLevel.normal.value
+            ):
+                print(
+                    f"WARNING: gnc package with name {pname} already exists. "
+                    "Replacing existing gnc package."
+                )
+            self._remove_package(self._ghost_node_files[package.filename])
+            del self._ghost_node_files[package.filename]
+        elif (
+            package.package_type.lower() == "mvr"
+            and package.filename in self._mover_files
+            and self._mover_files[package.filename] in self._packagelist
+        ):
+            # mvr package with same file name already exists.  remove old
+            # mvr package
+            if (
+                self.simulation_data.verbosity_level.value
+                >= VerbosityLevel.normal.value
+            ):
+                print(
+                    f"WARNING: mvr package with name {pname} already exists. "
+                    "Replacing existing mvr package."
+                )
+            self._remove_package(self._mover_files[package.filename])
+            del self._mover_files[package.filename]
+        elif (
+            package.package_type.lower() != "ims"
+            and pname in self.package_name_dict
+        ):
+            if (
+                self.simulation_data.verbosity_level.value
+                >= VerbosityLevel.normal.value
+            ):
+                print(
+                    "WARNING: Package with name "
+                    f"{package.package_name.lower()} already exists.  "
+                    "Replacing existing package."
+                )
+            self._remove_package(self.package_name_dict[pname])
+
     def register_package(
         self,
         package,
@@ -1694,7 +1944,8 @@ def register_package(
         set_package_name=True,
         set_package_filename=True,
     ):
-        """Register a package file with the simulation.  This is a
+        """
+        Register a package file with the simulation.  This is a
         call-back method made from the package and should not be called
         directly.
 
@@ -1717,76 +1968,7 @@ def register_package(
         package.container_type = [PackageContainerType.simulation]
         path = self._get_package_path(package)
         if add_to_package_list and package.package_type.lower != "nam":
-            pname = None
-            if package.package_name is not None:
-                pname = package.package_name.lower()
-            if (
-                package.package_type.lower() == "tdis"
-                and self._tdis_file is not None
-                and self._tdis_file in self._packagelist
-            ):
-                # tdis package already exists. there can be only one tdis
-                # package.  remove existing tdis package
-                if (
-                    self.simulation_data.verbosity_level.value
-                    >= VerbosityLevel.normal.value
-                ):
-                    print(
-                        "WARNING: tdis package already exists. Replacing "
-                        "existing tdis package."
-                    )
-                self._remove_package(self._tdis_file)
-            elif (
-                package.package_type.lower() == "gnc"
-                and package.filename in self._ghost_node_files
-                and self._ghost_node_files[package.filename]
-                in self._packagelist
-            ):
-                # gnc package with same file name already exists.  remove old
-                # gnc package
-                if (
-                    self.simulation_data.verbosity_level.value
-                    >= VerbosityLevel.normal.value
-                ):
-                    print(
-                        "WARNING: gnc package with name {} already exists. "
-                        "Replacing existing gnc package"
-                        ".".format(pname)
-                    )
-                self._remove_package(self._ghost_node_files[package.filename])
-                del self._ghost_node_files[package.filename]
-            elif (
-                package.package_type.lower() == "mvr"
-                and package.filename in self._mover_files
-                and self._mover_files[package.filename] in self._packagelist
-            ):
-                # mvr package with same file name already exists.  remove old
-                # mvr package
-                if (
-                    self.simulation_data.verbosity_level.value
-                    >= VerbosityLevel.normal.value
-                ):
-                    print(
-                        "WARNING: mvr package with name {} already exists. "
-                        "Replacing existing mvr package"
-                        ".".format(pname)
-                    )
-                self._remove_package(self._mover_files[package.filename])
-                del self._mover_files[package.filename]
-            elif (
-                package.package_type.lower() != "ims"
-                and pname in self.package_name_dict
-            ):
-                if (
-                    self.simulation_data.verbosity_level.value
-                    >= VerbosityLevel.normal.value
-                ):
-                    print(
-                        "WARNING: Package with name {} already exists.  "
-                        "Replacing existing package"
-                        ".".format(package.package_name.lower())
-                    )
-                self._remove_package(self.package_name_dict[pname])
+            self._remove_package_by_type(package)
             if package.package_type.lower() != "ims":
                 # all but ims packages get added here.  ims packages are
                 # added during ims package registration
@@ -1795,23 +1977,7 @@ def register_package(
             return path, self.structure.name_file_struct_obj
         elif package.package_type.lower() == "tdis":
             self._tdis_file = package
-            struct_root = mfstructure.MFStructure()
-            tdis_pkg = f"tdis{struct_root.get_version_string()}"
-            tdis_attr = getattr(self.name_file, tdis_pkg)
-            try:
-                tdis_attr.set_data(package.filename)
-            except MFDataException as mfde:
-                message = (
-                    "An error occurred while setting the tdis package "
-                    'file name "{}".  The error occurred while '
-                    "registering the tdis package with the "
-                    "simulation".format(package.filename)
-                )
-                raise MFDataException(
-                    mfdata_except=mfde,
-                    package=package._get_pname(),
-                    message=message,
-                )
+            self._set_timing_block(package.filename)
             return (
                 path,
                 self.structure.package_struct_objs[
@@ -1880,8 +2046,29 @@ def register_package(
             print(excpt_str)
             raise FlopyException(excpt_str)
 
+    def rename_model_namefile(self, model, new_namefile):
+        """
+        Rename a model's namefile.  For internal flopy library use only.
+
+        Parameters
+        ----------
+           model : MFModel
+               Model object whose namefile to rename
+           new_namefile : str
+               Name of the new namefile
+
+        """
+        # update simulation name file
+        models = self.name_file.models.get_data()
+        for mdl in models:
+            path, name_file_name = os.path.split(mdl[1])
+            if name_file_name == model.name_file.filename:
+                mdl[1] = os.path.join(path, new_namefile)
+        self.name_file.models.set_data(models)
+
     def register_model(self, model, model_type, model_name, model_namefile):
-        """Add a model to the simulation. This is a call-back method made
+        """
+        Add a model to the simulation. This is a call-back method made
         from the package and should not be called directly.
 
         Parameters
@@ -1932,7 +2119,8 @@ def register_model(self, model, model_type, model_name, model_namefile):
         return self.structure.model_struct_objs[model_type]
 
     def get_ims_package(self, key):
-        """Get the ims package with the specified `key`.
+        """
+        Get the ims package with the specified `key`.
 
         Parameters
         ----------
@@ -1949,7 +2137,8 @@ def get_ims_package(self, key):
         return None
 
     def remove_model(self, model_name):
-        """Remove model with name `model_name` from the simulation
+        """
+        Remove model with name `model_name` from the simulation
 
         Parameters
         ----------
@@ -1964,7 +2153,8 @@ def remove_model(self, model_name):
         # Update simulation name file
 
     def is_valid(self):
-        """Checks the validity of the solution and all of its models and
+        """
+        Checks the validity of the solution and all of its models and
         packages.  Returns true if the solution is valid, false if it is not.
 
 
@@ -2023,7 +2213,7 @@ def _get_package_path(package):
         else:
             return (package.package_type,)
 
-    def _remove_ims_soultion_group(self, ims_file):
+    def _update_ims_solution_group(self, ims_file, new_name=None):
         solution_recarray = self.name_file.solutiongroup
         for solution_group_num in solution_recarray.get_active_key_list():
             try:
@@ -2041,7 +2231,11 @@ def _remove_ims_soultion_group(self, ims_file):
             new_array = []
             for record in rec_array:
                 if record.slnfname == ims_file:
-                    continue
+                    if new_name is not None:
+                        record.slnfname = new_name
+                        new_array.append(tuple(record))
+                    else:
+                        continue
                 else:
                     new_array.append(record)
 
@@ -2128,7 +2322,8 @@ def _is_in_solution_group(self, item, index):
         return False
 
     def plot(self, model_list=None, SelPackList=None, **kwargs):
-        """Plot simulation or models.
+        """
+        Plot simulation or models.
 
         Method to plot a whole simulation or a series of models
         that are part of a simulation.
diff --git a/flopy/utils/datautil.py b/flopy/utils/datautil.py
index 594b111e45..71bf50e4ca 100644
--- a/flopy/utils/datautil.py
+++ b/flopy/utils/datautil.py
@@ -260,6 +260,13 @@ def array_comp(self, first_array, second_array):
             return False
         return True
 
+    def riv_array_comp(self, first_array, second_array):
+        for line_first, line_second in zip(first_array, second_array):
+            diff = np.abs(line_first[0][2] - line_second[0][2])
+            if diff > self.max_error:
+                return False
+        return True
+
     @staticmethod
     def reset_delimiter_used():
         PyListUtil.delimiter_used = None