diff --git a/autotest/test_mf6.py b/autotest/test_mf6.py
index 5e848a2f68..080c939301 100644
--- a/autotest/test_mf6.py
+++ b/autotest/test_mf6.py
@@ -1685,12 +1685,41 @@ def test_sfr_connections(function_tmpdir, example_data_path):
 
         # reload simulation
         sim2 = MFSimulation.load(sim_ws=sim_ws)
-        sim.set_all_data_external()
-        sim.write_simulation()
-        success, buff = sim.run_simulation()
+        sim2.set_all_data_external()
+        sim2.write_simulation()
+        success, buff = sim2.run_simulation()
         assert (
             success
-        ), f"simulation {sim.name} did not run after being reloaded"
+        ), f"simulation {sim2.name} did not run after being reloaded"
+
+        # test sfr recarray data
+        model2 = sim2.get_model()
+        sfr2 = model2.get_package("sfr")
+        sfr_pd = sfr2.packagedata
+        rec_data = [
+            (0, 0, 0, 0, 1.0, 1.0, 0.01, 10.0, 1.0, 1.0, 1.0, 1, 1.0, 0),
+            (1, 0, 1, 0, 1.0, 1.0, 0.01, 10.0, 1.0, 1.0, 1.0, 2, 1.0, 0),
+        ]
+        rec_type = [
+            ("ifno", int),
+            ("layer", int),
+            ("row", int),
+            ("column", int),
+            ("rlen", float),
+            ("rwid", float),
+            ("rgrd", float),
+            ("rtp", float),
+            ("rbth", float),
+            ("rhk", float),
+            ("man", float),
+            ("nconn", int),
+            ("ustrf", float),
+            ("nvd", int),
+        ]
+        pkg_data = np.rec.array(rec_data, rec_type)
+        sfr_pd.set_record({"data": pkg_data})
+        data = sfr_pd.get_data()
+        assert data[0][1] == (0, 0, 0)
 
 
 @requires_exe("mf6")
@@ -2232,24 +2261,16 @@ def test_multi_model(function_tmpdir):
     assert fi_out[2][2] == "MIXED"
 
     spca1 = ModflowUtlspca(
-        gwt2,
-        filename="gwt_model_1.rch1.spc",
-        print_input=True
+        gwt2, filename="gwt_model_1.rch1.spc", print_input=True
     )
     spca2 = ModflowUtlspca(
-        gwt2,
-        filename="gwt_model_1.rch2.spc",
-        print_input=False
+        gwt2, filename="gwt_model_1.rch2.spc", print_input=False
     )
     spca3 = ModflowUtlspca(
-        gwt2,
-        filename="gwt_model_1.rch3.spc",
-        print_input=True
+        gwt2, filename="gwt_model_1.rch3.spc", print_input=True
     )
     spca4 = ModflowUtlspca(
-        gwt2,
-        filename="gwt_model_1.rch4.spc",
-        print_input=True
+        gwt2, filename="gwt_model_1.rch4.spc", print_input=True
     )
 
     # test writing and loading spca packages
diff --git a/flopy/mf6/coordinates/modeldimensions.py b/flopy/mf6/coordinates/modeldimensions.py
index c41003592b..5b9460ecf8 100644
--- a/flopy/mf6/coordinates/modeldimensions.py
+++ b/flopy/mf6/coordinates/modeldimensions.py
@@ -615,10 +615,10 @@ def _resolve_data_item_shape(
                                 if data is None:
                                     if (
                                         self.simulation_data.verbosity_level.value
-                                        >= VerbosityLevel.normal.value
+                                        >= VerbosityLevel.verbose.value
                                     ):
                                         print(
-                                            "WARNING: Unable to resolve "
+                                            "INFORMATION: Unable to resolve "
                                             "dimension of {} based on shape "
                                             '"{}".'.format(
                                                 data_item_struct.path, item[0]
@@ -651,10 +651,10 @@ def _resolve_data_item_shape(
                             else:
                                 if (
                                     self.simulation_data.verbosity_level.value
-                                    >= VerbosityLevel.normal.value
+                                    >= VerbosityLevel.verbose.value
                                 ):
                                     print(
-                                        "WARNING: Unable to resolve "
+                                        "INFORMATION: Unable to resolve "
                                         "dimension of {} based on shape "
                                         '"{}".'.format(
                                             data_item_struct.path, item[0]
diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py
index fa0836076e..8fef1c2b7e 100644
--- a/flopy/mf6/data/mfdatastorage.py
+++ b/flopy/mf6/data/mfdatastorage.py
@@ -1276,6 +1276,10 @@ def store_internal(
                 self.layer_storage.first_item().data_storage_type = (
                     DataStorageType.internal_array
                 )
+                if data is None or isinstance(data, np.recarray):
+                    if not self.tuple_cellids(data):
+                        # fix data so cellid is a single tuple
+                        data = self.make_tuple_cellids(data.tolist())
                 if data is None or isinstance(data, np.recarray):
                     if self._simulation_data.verify_data and check_data:
                         self._verify_list(data)
@@ -1602,9 +1606,15 @@ def make_tuple_cellids(self, data):
         return new_data
 
     def tuple_cellids(self, data):
+        if data is None or len(data) == 0:
+            return True
         for data_entry, cellid in zip(data[0], self.recarray_cellid_list):
             if cellid:
-                if isinstance(data_entry, int):
+                if (
+                    isinstance(data_entry, int)
+                    or isinstance(data_entry, np.int32)
+                    or isinstance(data_entry, np.int64)
+                ):
                     # cellid is stored in separate columns in the recarray
                     # (eg: one column for layer one column for row and
                     # one columne for column)