diff --git a/autotest/test_mf6.py b/autotest/test_mf6.py
index a2e473167..6cc9e9c3e 100644
--- a/autotest/test_mf6.py
+++ b/autotest/test_mf6.py
@@ -651,6 +651,286 @@ def test_create_and_run_model(tmpdir):
     assert success
 
 
+@requires_exe("mf6")
+def test_get_set_data_record(tmpdir):
+    # names
+    sim_name = "testrecordsim"
+    model_name = "testrecordmodel"
+    exe_name = "mf6"
+
+    # set up simulation
+    tdis_name = f"{sim_name}.tdis"
+    sim = MFSimulation(
+        sim_name=sim_name, version="mf6", exe_name=exe_name, sim_ws=str(tmpdir)
+    )
+    tdis_rc = [(10.0, 4, 1.0), (6.0, 3, 1.0)]
+    tdis = mftdis.ModflowTdis(
+        sim, time_units="DAYS", nper=2, perioddata=tdis_rc
+    )
+
+    # create model instance
+    model = mfgwf.ModflowGwf(
+        sim, modelname=model_name, model_nam_file=f"{model_name}.nam"
+    )
+
+    # create solution and add the model
+    ims_package = mfims.ModflowIms(
+        sim,
+        print_option="ALL",
+        complexity="SIMPLE",
+        outer_dvclose=0.00001,
+        outer_maximum=50,
+        under_relaxation="NONE",
+        inner_maximum=30,
+        inner_dvclose=0.00001,
+        linear_acceleration="CG",
+        preconditioner_levels=7,
+        preconditioner_drop_tolerance=0.01,
+        number_orthogonalizations=2,
+    )
+    sim.register_ims_package(ims_package, [model_name])
+
+    # add packages to model
+    dis_package = mfgwfdis.ModflowGwfdis(
+        model,
+        length_units="FEET",
+        nlay=3,
+        nrow=10,
+        ncol=10,
+        delr=500.0,
+        delc=500.0,
+        top=100.0,
+        botm=[50.0, 10.0, -50.0],
+        filename=f"{model_name}.dis",
+    )
+    ic_package = mfgwfic.ModflowGwfic(
+        model,
+        strt=[100.0, 90.0, 80.0],
+        filename=f"{model_name}.ic",
+    )
+    npf_package = mfgwfnpf.ModflowGwfnpf(
+        model, save_flows=True, icelltype=1, k=50.0, k33=1.0
+    )
+
+    sto_package = mfgwfsto.ModflowGwfsto(
+        model, save_flows=True, iconvert=1, ss=0.000001, sy=0.15
+    )
+    # wel packages
+    period_one = ModflowGwfwel.stress_period_data.empty(
+        model,
+        maxbound=3,
+        aux_vars=["var1", "var2", "var3"],
+        boundnames=True,
+        timeseries=True,
+    )
+    period_one[0][0] = ((0, 9, 2), -50.0, -1, -2, -3, None)
+    period_one[0][1] = ((1, 4, 7), -100.0, 1, 2, 3, "well_1")
+    period_one[0][2] = ((1, 3, 2), -20.0, 4, 5, 6, "well_2")
+    period_two = ModflowGwfwel.stress_period_data.empty(
+        model,
+        maxbound=2,
+        aux_vars=["var1", "var2", "var3"],
+        boundnames=True,
+        timeseries=True,
+    )
+    period_two[0][0] = ((2, 3, 2), -80.0, 1, 2, 3, "well_2")
+    period_two[0][1] = ((2, 4, 7), -10.0, 4, 5, 6, "well_1")
+    stress_period_data = {}
+    stress_period_data[0] = period_one[0]
+    stress_period_data[1] = period_two[0]
+    wel_package = ModflowGwfwel(
+        model,
+        print_input=True,
+        print_flows=True,
+        auxiliary=[("var1", "var2", "var3")],
+        maxbound=5,
+        stress_period_data=stress_period_data,
+        boundnames=True,
+        save_flows=True,
+    )
+    # rch package
+    rch_period_list = []
+    for row in range(0, 10):
+        for col in range(0, 10):
+            rch_amt = (1 + row / 10) * (1 + col / 10)
+            rch_period_list.append(((0, row, col), rch_amt, 0.5))
+    rch_period = {}
+    rch_period[0] = rch_period_list
+    rch_package = ModflowGwfrch(
+        model,
+        fixed_cell=True,
+        auxiliary="MULTIPLIER",
+        auxmultname="MULTIPLIER",
+        print_input=True,
+        print_flows=True,
+        save_flows=True,
+        maxbound=54,
+        stress_period_data=rch_period,
+    )
+
+    # write simulation to new location
+    sim.set_all_data_external()
+    sim.write_simulation()
+
+    # test get_record, set_record for list data
+    wel = model.get_package("wel")
+    spd_record = wel.stress_period_data.get_record()
+    well_sp_1 = spd_record[0]
+    assert (
+        well_sp_1["filename"] == "testrecordmodel.wel_stress_period_data_1.txt"
+    )
+    assert well_sp_1["binary"] is False
+    assert well_sp_1["data"][0][0] == (0, 9, 2)
+    assert well_sp_1["data"][0][1] == -50.0
+    # modify
+    del well_sp_1["filename"]
+    well_sp_1["data"][0][0] = (1, 9, 2)
+    well_sp_2 = spd_record[1]
+    del well_sp_2["filename"]
+    well_sp_2["data"][0][0] = (1, 1, 1)
+    # save
+    spd_record[0] = well_sp_1
+    spd_record[1] = well_sp_2
+    wel.stress_period_data.set_record(spd_record)
+    # verify changes
+    spd_record = wel.stress_period_data.get_record()
+    well_sp_1 = spd_record[0]
+    assert "filename" not in well_sp_1
+    assert well_sp_1["data"][0][0] == (1, 9, 2)
+    assert well_sp_1["data"][0][1] == -50.0
+    well_sp_2 = spd_record[1]
+    assert "filename" not in well_sp_2
+    assert well_sp_2["data"][0][0] == (1, 1, 1)
+    spd = wel.stress_period_data.get_data()
+    assert spd[0][0][0] == (1, 9, 2)
+    # change well_sp_2 back to external
+    well_sp_2["filename"] = "wel_spd_data_2.txt"
+    spd_record[1] = well_sp_2
+    wel.stress_period_data.set_record(spd_record)
+    # change well_sp_2 data
+    spd[1][0][0] = (1, 2, 2)
+    wel.stress_period_data.set_data(spd)
+    # verify changes
+    spd_record = wel.stress_period_data.get_record()
+    well_sp_2 = spd_record[1]
+    assert well_sp_2["filename"] == "wel_spd_data_2.txt"
+    assert well_sp_2["data"][0][0] == (1, 2, 2)
+
+    # test get_data/set_data vs get_record/set_record
+    dis = model.get_package("dis")
+    botm = dis.botm.get_record()
+    assert len(botm) == 3
+    layer_2 = botm[1]
+    layer_3 = botm[2]
+    # verify layer 2
+    assert layer_2["filename"] == "testrecordmodel.dis_botm_layer2.txt"
+    assert layer_2["binary"] is False
+    assert layer_2["factor"] == 1.0
+    assert layer_2["iprn"] is None
+    assert layer_2["data"][0][0] == 10.0
+    # change and set layer 2
+    layer_2["filename"] = "botm_layer2.txt"
+    layer_2["binary"] = True
+    layer_2["iprn"] = 3
+    layer_2["factor"] = 2.0
+    layer_2["data"] = layer_2["data"] * 0.5
+    botm[1] = layer_2
+    # change and set layer 3
+    del layer_3["filename"]
+    layer_3["factor"] = 0.5
+    layer_3["data"] = layer_3["data"] * 2.0
+    botm[2] = layer_3
+    dis.botm.set_record(botm)
+
+    # get botm in two different ways, verifying changes made
+    botm_record = dis.botm.get_record()
+    layer_1 = botm_record[0]
+    assert layer_1["filename"] == "testrecordmodel.dis_botm_layer1.txt"
+    assert layer_1["binary"] is False
+    assert layer_1["iprn"] is None
+    assert layer_1["data"][0][0] == 50.0
+    layer_2 = botm_record[1]
+    assert layer_2["filename"] == "botm_layer2.txt"
+    assert layer_2["binary"] is True
+    assert layer_2["factor"] == 2.0
+    assert layer_2["iprn"] == 3
+    assert layer_2["data"][0][0] == 5.0
+    layer_3 = botm_record[2]
+    assert "filename" not in layer_3
+    assert layer_3["factor"] == 0.5
+    assert layer_3["data"][0][0] == -100.0
+    botm_data = dis.botm.get_data(apply_mult=True)
+    assert botm_data[0][0][0] == 50.0
+    assert botm_data[1][0][0] == 10.0
+    assert botm_data[2][0][0] == -50.0
+    botm_data = dis.botm.get_data()
+    assert botm_data[0][0][0] == 50.0
+    assert botm_data[1][0][0] == 5.0
+    assert botm_data[2][0][0] == -100.0
+    # modify and set botm data with set_data
+    botm_data[0][0][0] = 6.0
+    botm_data[1][0][0] = -8.0
+    botm_data[2][0][0] = -205.0
+    dis.botm.set_data(botm_data)
+    # verify that data changed and metadata did not change
+    botm_record = dis.botm.get_record()
+    layer_1 = botm_record[0]
+    assert layer_1["filename"] == "testrecordmodel.dis_botm_layer1.txt"
+    assert layer_1["binary"] is False
+    assert layer_1["iprn"] is None
+    assert layer_1["data"][0][0] == 6.0
+    assert layer_1["data"][0][1] == 50.0
+    layer_2 = botm_record[1]
+    assert layer_2["filename"] == "botm_layer2.txt"
+    assert layer_2["binary"] is True
+    assert layer_2["factor"] == 2.0
+    assert layer_2["iprn"] == 3
+    assert layer_2["data"][0][0] == -8.0
+    assert layer_2["data"][0][1] == 5.0
+    layer_3 = botm_record[2]
+    assert "filename" not in layer_3
+    assert layer_3["factor"] == 0.5
+    assert layer_3["data"][0][0] == -205.0
+    botm_data = dis.botm.get_data()
+    assert botm_data[0][0][0] == 6.0
+    assert botm_data[1][0][0] == -8.0
+    assert botm_data[2][0][0] == -205.0
+
+    spd_record = rch_package.stress_period_data.get_record()
+    assert 0 in spd_record
+    assert isinstance(spd_record[0], dict)
+    assert "filename" in spd_record[0]
+    assert (
+        spd_record[0]["filename"]
+        == "testrecordmodel.rch_stress_period_data_1.txt"
+    )
+    assert "binary" in spd_record[0]
+    assert spd_record[0]["binary"] is False
+    assert "data" in spd_record[0]
+    assert spd_record[0]["data"][0][0] == (0, 0, 0)
+    spd_record[0]["data"][0][0] = (0, 0, 8)
+    rch_package.stress_period_data.set_record(spd_record)
+
+    spd_data = rch_package.stress_period_data.get_data()
+    assert spd_data[0][0][0] == (0, 0, 8)
+    spd_data[0][0][0] = (0, 0, 7)
+    rch_package.stress_period_data.set_data(spd_data)
+
+    spd_record = rch_package.stress_period_data.get_record()
+    assert isinstance(spd_record[0], dict)
+    assert "filename" in spd_record[0]
+    assert (
+        spd_record[0]["filename"]
+        == "testrecordmodel.rch_stress_period_data_1.txt"
+    )
+    assert "binary" in spd_record[0]
+    assert spd_record[0]["binary"] is False
+    assert "data" in spd_record[0]
+    assert spd_record[0]["data"][0][0] == (0, 0, 7)
+
+    sim.write_simulation()
+
+
 @requires_exe("mf6")
 def test_output(tmpdir, example_data_path):
     ex_name = "test001e_UZF_3lay"
diff --git a/examples/Tutorials/modflow6data/tutorial08_mf6_data.py b/examples/Tutorials/modflow6data/tutorial08_mf6_data.py
index 791f1a454..c1e7b1876 100644
--- a/examples/Tutorials/modflow6data/tutorial08_mf6_data.py
+++ b/examples/Tutorials/modflow6data/tutorial08_mf6_data.py
@@ -145,25 +145,72 @@
 }
 
 # A list containing data for the three specified layers is then passed in to
-# the `botm` object's `set_data` method.
+# the `botm` object's `set_record` method.
 
-dis.botm.set_data([a0, a1, a2])
+dis.botm.set_record([a0, a1, a2])
 print(dis.botm.get_file_entry())
 
-# Note that we could have also specified `botm` this way as part of the
-# original `flopy.mf6.ModflowGwfdis` constructor:
+# The botm data and its attributes (filename, binary, factor, iprn) can be
+# retrieved as a dictionary of dictionaries using get_record.
 
-a0 = {"factor": 0.5, "iprn": 1, "data": np.ones((10, 10))}
-a1 = -100
-a2 = {
-    "filename": "dis.botm.3.bin",
-    "factor": 2.0,
-    "iprn": 1,
-    "data": -100 * np.ones((4, 5)),
-    "binary": True,
+botm_record = dis.botm.get_record()
+print("botm layer 1 record:")
+print(botm_record[0])
+print("\nbotm layer 2 record:")
+print(botm_record[1])
+print("\nbotm layer 3 record:")
+print(botm_record[2])
+
+# The botm record retrieved can be modified and then saved with set_record.
+# For example, the array data's "factor" can be modified and saved.
+
+botm_record[0]["factor"] = 0.6
+dis.botm.set_record(botm_record)
+
+# The updated value can then be retrieved.
+
+botm_record = dis.botm.get_record()
+print(f"botm layer 1 factor:  {botm_record[0]['factor']}")
+
+# The get_record and set_record methods can also be used with list data to get
+# and set the data and its "filename" and "binary" attributes. This is
+# demonstrated with the wel package.  First, a wel package is constructed.
+
+welspdict = {
+    0: {"filename": "well_sp1.txt", "data": [[(0, 0, 0), 0.25]]},
+    1: [[(0, 0, 0), 0.1]],
 }
-botm = [a0, a1, a2]
-flopy.mf6.ModflowGwfdis(gwf, nlay=3, nrow=10, ncol=10, botm=botm)
+wel = flopy.mf6.ModflowGwfwel(
+    gwf,
+    print_input=True,
+    print_flows=True,
+    stress_period_data=welspdict,
+    save_flows=False,
+)
+
+# The wel stress period data and associated "filename" and "binary" attributes
+# can be retrieved with get_record.
+
+spd_record = wel.stress_period_data.get_record()
+print("Stress period 1 record:")
+print(spd_record[0])
+print("\nStress period 2 record:")
+print(spd_record[1])
+
+# The wel data and associated attributes can be changed by modifying the
+# record and then saving it with the set_record method.
+
+spd_record[0]["filename"] = "well_package_sp1.txt"
+spd_record[0]["binary"] = True
+spd_record[1]["filename"] = "well_package_sp2.txt"
+wel.stress_period_data.set_record(spd_record)
+
+# The changes can be verified by calling get_record again.
+
+spd_record = wel.stress_period_data.get_record()
+print(f"New filename for stress period 1:  {spd_record[0]['filename']}")
+print(f"New binary flag for stress period 1:  {spd_record[0]['binary']}")
+print(f"New filename for stress period 2:  {spd_record[1]['filename']}")
 
 try:
     temp_dir.cleanup()
diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py
index 1fcb722c1..b57799ff0 100644
--- a/flopy/mf6/data/mfdataarray.py
+++ b/flopy/mf6/data/mfdataarray.py
@@ -2,6 +2,7 @@
 import inspect
 import os
 import sys
+import warnings
 
 import numpy as np
 
@@ -640,12 +641,12 @@ def store_internal(
                     )
                 factor = storage.layer_storage[current_layer].factor
                 internal_data = {
-                    "data": self._get_data(current_layer, True),
-                    "factor": factor,
+                    current_layer[0]: {
+                        "data": self._get_data(current_layer, True),
+                        "factor": factor,
+                    }
                 }
-                self._set_data(
-                    internal_data, layer=current_layer, check_data=False
-                )
+                self._set_record(internal_data)
             except Exception as ex:
                 type_, value_, traceback_ = sys.exc_info()
                 raise MFDataException(
@@ -703,12 +704,12 @@ def has_data(self, layer=None):
             )
 
     def get_data(self, layer=None, apply_mult=False, **kwargs):
-        """Returns the data associated with layer "layer_num".  If "layer_num"
+        """Returns the data associated with layer "layer".  If "layer"
         is None, returns all data.
 
         Parameters
         ----------
-            layer_num : int
+            layer : int
 
         Returns
         -------
@@ -753,6 +754,126 @@ def _get_data(self, layer=None, apply_mult=False, **kwargs):
                 )
         return None
 
+    def get_record(self, layer=None):
+        """Returns the data record associated with layer "layer".  If "layer"
+        is None, returns all data.
+
+        Parameters
+        ----------
+            layer : int
+
+        Returns
+        -------
+            data_record : dict
+                Dictionary containing data record for specified layer or
+                dictionary containing layer numbers (keys) and data record
+                for that layer (values).
+
+        """
+        if self._get_storage_obj() is None:
+            self._data_storage = self._new_storage(False)
+        if isinstance(layer, int):
+            layer = (layer,)
+        storage = self._get_storage_obj()
+        if storage is not None:
+            try:
+                return storage.get_record(layer)
+            except Exception as ex:
+                type_, value_, traceback_ = sys.exc_info()
+                raise MFDataException(
+                    self.structure.get_model(),
+                    self.structure.get_package(),
+                    self._path,
+                    "getting record",
+                    self.structure.name,
+                    inspect.stack()[0][3],
+                    type_,
+                    value_,
+                    traceback_,
+                    None,
+                    self._simulation_data.debug,
+                    ex,
+                )
+        return None
+
+    def set_record(self, data_record):
+        """Sets the contents of the data and metadata to the contents of
+        'data_record`.  For unlayered data do not pass in `layer`.  For
+        unlayered data 'data_record' is a dictionary with either of the
+        following key/value pairs:
+        1) {'filename':filename, 'factor':fct, 'iprn':print, 'data':data} -
+        dictionary defining external file information
+        2) {'data':data, 'factor':fct, 'iprn':print)}- dictionary defining
+        internal information. Data that is layered can also be set by defining
+        For layered data data_record must be a dictionary of dictionaries
+        with zero-based layer numbers for the outer dictionary keys and the
+        inner dictionary as described above:
+           {0: {'data':data_lay_0, 'factor':fct_lay_0, 'iprn':prn_lay_0)},
+            1: {'data':data_lay_1, 'factor':fct_lay_1, 'iprn':prn_lay_1)},
+            2: {'data':data_lay_2, 'factor':fct_lay_2, 'iprn':prn_lay_2)}}
+
+        Parameters
+        ----------
+        data_record : dict
+            An dictionary of data record information or a dictionary of
+            layers (keys) with a dictionary of data record information for each
+            layer (values).
+        """
+        self._set_record(data_record)
+
+    def _set_record(self, data_record):
+        """Sets the contents of the data and metadata to the contents of
+        'data_record`.  For unlayered data do not pass in `layer`.  For
+        unlayered data 'data_record' is a dictionary with either of the
+        following key/value pairs:
+        1) {'filename':filename, 'factor':fct, 'iprn':print, 'data':data} -
+        dictionary defining external file information
+        2) {'data':data, 'factor':fct, 'iprn':print)}- dictionary defining
+        internal information. Data that is layered can also be set by defining
+        For layered data data_record must be a list of dictionaries or a
+        dictionary of dictionaries with zero-based layer numbers for the outer
+        dictionary keys and the inner dictionary as described above:
+           {0: {'data':data_lay_0, 'factor':fct_lay_0, 'iprn':prn_lay_0)},
+            1: {'data':data_lay_1, 'factor':fct_lay_1, 'iprn':prn_lay_1)},
+            2: {'data':data_lay_2, 'factor':fct_lay_2, 'iprn':prn_lay_2)}}
+
+        Parameters
+        ----------
+        data_record : dict
+            An dictionary of data record information or a dictionary of
+            layers (keys) with a dictionary of data record information for each
+            layer (values).
+        """
+        if isinstance(data_record, dict):
+            first_key = list(data_record.keys())[0]
+            if isinstance(first_key, int):
+                for layer, record in data_record.items():
+                    self._set_data(record, layer=layer, preserve_record=False)
+            else:
+                self._set_data(data_record, preserve_record=False)
+        elif type(data_record) == list:
+            for layer, record in enumerate(data_record):
+                self._set_data(record, layer=layer, preserve_record=False)
+        else:
+            message = (
+                "Unable to set record.  The contents of data_record must be"
+                "a dictionary or a list."
+            )
+            type_, value_, traceback_ = sys.exc_info()
+            raise MFDataException(
+                self._data_dimensions.structure.get_model(),
+                self._data_dimensions.structure.get_package(),
+                self._data_dimensions.structure.path,
+                "setting record",
+                self._data_dimensions.structure.name,
+                inspect.stack()[0][3],
+                type_,
+                value_,
+                traceback_,
+                message,
+                self._simulation_data.debug,
+            )
+
     def set_data(self, data, multiplier=None, layer=None):
         """Sets the contents of the data at layer `layer` to `data` with
         multiplier `multiplier`. For unlayered data do not pass in
@@ -760,15 +881,6 @@ def set_data(self, data, multiplier=None, layer=None):
         1) ndarray - numpy ndarray containing all of the data
         2) [data] - python list containing all of the data
         3) val - a single constant value to be used for all of the data
-        4) {'filename':filename, 'factor':fct, 'iprn':print, 'data':data} -
-        dictionary defining external file information
-        5) {'data':data, 'factor':fct, 'iprn':print) - dictionary defining
-        internal information. Data that is layered can also be set by defining
-        a list with a length equal to the number of layers in the model.
-        Each layer in the list contains the data as defined in the
-        formats above:
-            [layer_1_val, [layer_2_array_vals],
-            {'filename':file_with_layer_3_data, 'factor':fct, 'iprn':print}]
 
         Parameters
         ----------
@@ -782,12 +894,17 @@ def set_data(self, data, multiplier=None, layer=None):
         """
         self._set_data(data, multiplier, layer)
 
-    def _set_data(self, data, multiplier=None, layer=None, check_data=True):
+    def _set_data(
+        self,
+        data,
+        multiplier=None,
+        layer=None,
+        check_data=True,
+        preserve_record=True,
+    ):
         self._resync()
         if self._get_storage_obj() is None:
             self._data_storage = self._new_storage(False)
-        if multiplier is None:
-            multiplier = [self._get_storage_obj().get_default_mult()]
         if isinstance(layer, int):
             layer = (layer,)
         if isinstance(data, str):
@@ -833,7 +950,11 @@ def _set_data(self, data, multiplier=None, layer=None, check_data=True):
                         layer_data = aux_var_data
                     try:
                         storage.set_data(
-                            layer_data, [layer], multiplier, self._current_key
+                            layer_data,
+                            [layer],
+                            multiplier,
+                            self._current_key,
+                            preserve_record=preserve_record,
                         )
                     except Exception as ex:
                         type_, value_, traceback_ = sys.exc_info()
@@ -874,7 +995,11 @@ def _set_data(self, data, multiplier=None, layer=None, check_data=True):
         else:
             try:
                 storage.set_data(
-                    data, layer, multiplier, key=self._current_key
+                    data,
+                    layer,
+                    multiplier,
+                    key=self._current_key,
+                    preserve_record=preserve_record,
                 )
             except Exception as ex:
                 type_, value_, traceback_ = sys.exc_info()
@@ -1662,20 +1787,50 @@ def has_data(self, layer=None):
             self.get_data_prep(layer)
             return super().has_data()
 
-    def get_data(self, layer=None, apply_mult=True, **kwargs):
-        """Returns the data associated with stress period key `layer`.
-        If `layer` is None, returns all data for time `layer`.
+    def get_record(self, key=None):
+        """Returns the data record (data and metadata) associated with stress
+        period key `layer`.  If `layer` is None, returns all data for
+        time `layer`.
 
         Parameters
         ----------
-            layer : int
+            key : int
+                Zero-based stress period of data to return
+
+        """
+        if self._data_storage is not None and len(self._data_storage) > 0:
+            if key is None:
+                sim_time = self._data_dimensions.package_dim.model_dim[
+                    0
+                ].simulation_time
+                num_sp = sim_time.get_num_stress_periods()
+                return self._build_period_data(num_sp, get_record=True)
+            else:
+                self.get_data_prep(key)
+                return super().get_record()
+
+    def get_data(self, key=None, apply_mult=True, **kwargs):
+        """Returns the data associated with stress period key `key`.
+        If `layer` is None, returns all data for time `key`.
+
+        Parameters
+        ----------
+            key : int
                 Zero-based stress period of data to return
             apply_mult : bool
                 Whether to apply multiplier to data prior to returning it
 
         """
+        if "layer" in kwargs:
+            warnings.warn(
+                "The 'layer' parameter has been deprecated, use 'key' "
+                "instead.",
+                category=DeprecationWarning,
+            )
+            key = kwargs["layer"]
+
         if self._data_storage is not None and len(self._data_storage) > 0:
-            if layer is None:
+            if key is None:
                 sim_time = self._data_dimensions.package_dim.model_dim[
                     0
                 ].simulation_time
@@ -1683,31 +1838,54 @@ def get_data(self, layer=None, apply_mult=True, **kwargs):
                 if "array" in kwargs:
                     return self._get_array(num_sp, apply_mult, **kwargs)
                 else:
-                    output = None
-                    for sp in range(0, num_sp):
-                        data = None
-                        if sp in self._data_storage:
-                            self.get_data_prep(sp)
-                            data = super().get_data(
-                                apply_mult=apply_mult, **kwargs
-                            )
-                        if output is None:
-                            if "array" in kwargs:
-                                output = [data]
-                            else:
-                                output = {sp: data}
-                        else:
-                            if "array" in kwargs:
-                                output.append(data)
-                            else:
-                                output[sp] = data
-                    return output
+                    return self._build_period_data(
+                        num_sp, apply_mult, **kwargs
+                    )
             else:
-                self.get_data_prep(layer)
+                self.get_data_prep(key)
                 return super().get_data(apply_mult=apply_mult)
         else:
             return None
 
+    def _build_period_data(
+        self, num_sp, apply_mult=False, get_record=False, **kwargs
+    ):
+        output = None
+        for sp in range(0, num_sp):
+            data = None
+            if sp in self._data_storage:
+                self.get_data_prep(sp)
+                if get_record:
+                    data = super().get_record()
+                else:
+                    data = super().get_data(apply_mult=apply_mult, **kwargs)
+            if output is None:
+                if "array" in kwargs:
+                    output = [data]
+                else:
+                    output = {sp: data}
+            else:
+                if "array" in kwargs:
+                    output.append(data)
+                else:
+                    output[sp] = data
+        return output
+
+    def set_record(self, data_record):
+        """Sets data and metadata at layer `layer` and time `key` to
+        `data_record`.  For unlayered data do not pass in `layer`.
+
+        Parameters
+        ----------
+            data_record : dict
+                Data record being set.  Data must be a dictionary with keys as
+                zero-based stress periods and values as a dictionary of data
+                and metadata (factor, iprn, filename, binary, data) for a given
+                stress period.  How to define the dictionary of data and
+                metadata is described in the MFData class's set_record method.
+        """
+        self._set_data_record(data_record, is_record=True)
+
     def set_data(self, data, multiplier=None, layer=None, key=None):
         """Sets the contents of the data at layer `layer` and time `key` to
         `data` with multiplier `multiplier`. For unlayered data do not pass
@@ -1730,7 +1908,11 @@ def set_data(self, data, multiplier=None, layer=None, key=None):
                 Zero based stress period to assign data too.  Does not apply
                 if `data` is a dictionary.
         """
+        self._set_data_record(data, multiplier, layer, key)
 
+    def _set_data_record(
+        self, data, multiplier=None, layer=None, key=None, is_record=False
+    ):
         if isinstance(data, dict):
             # each item in the dictionary is a list for one stress period
             # the dictionary key is the stress period the list is for
@@ -1741,7 +1923,10 @@ def set_data(self, data, multiplier=None, layer=None, key=None):
                     del_keys.append(key)
                 else:
                     self._set_data_prep(list_item, key)
-                    super().set_data(list_item, multiplier, layer)
+                    if is_record:
+                        super().set_record(list_item)
+                    else:
+                        super().set_data(list_item, multiplier, layer)
             for key in del_keys:
                 del data[key]
         else:
diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py
index e6b4abe42..fe19e1dd0 100644
--- a/flopy/mf6/data/mfdatalist.py
+++ b/flopy/mf6/data/mfdatalist.py
@@ -318,7 +318,7 @@ def store_internal(
                 internal_data = {
                     "data": data,
                 }
-                self._set_data(internal_data, check_data=check_data)
+                self._set_record(internal_data, check_data=check_data)
 
     def has_data(self, key=None):
         """Returns whether this MFList has any data associated with it."""
@@ -383,6 +383,36 @@ def get_data(self, apply_mult=False, **kwargs):
         """
         return self._get_data(apply_mult, **kwargs)
 
+    def get_record(self):
+        """Returns the list's data and metadata in a dictionary.  Data is in
+        key "data" and metadata in keys "filename" and "binary".
+
+        Returns
+        -------
+            data_record : dict
+
+        """
+        try:
+            if self._get_storage_obj() is None:
+                return None
+            return self._get_storage_obj().get_record()
+        except Exception as ex:
+            type_, value_, traceback_ = sys.exc_info()
+            raise MFDataException(
+                self.structure.get_model(),
+                self.structure.get_package(),
+                self._path,
+                "getting record",
+                self.structure.name,
+                inspect.stack()[0][3],
+                type_,
+                value_,
+                traceback_,
+                None,
+                self._simulation_data.debug,
+                ex,
+            )
+
     def _get_min_record_entries(self, data=None):
         try:
             if isinstance(data, dict) and "data" in data:
@@ -410,7 +440,9 @@ def _get_min_record_entries(self, data=None):
             )
         return len(type_list)
 
-    def _set_data(self, data, autofill=False, check_data=True):
+    def _set_data(
+        self, data, autofill=False, check_data=True, preserve_record=True
+    ):
         # set data
         self._resync()
         try:
@@ -418,7 +450,10 @@ def _set_data(self, data, autofill=False, check_data=True):
                 self._data_storage = self._new_storage()
             # store data
             self._get_storage_obj().set_data(
-                data, autofill=autofill, check_data=check_data
+                data,
+                autofill=autofill,
+                check_data=check_data,
+                preserve_record=preserve_record,
             )
         except Exception as ex:
             type_, value_, traceback_ = sys.exc_info()
@@ -539,13 +574,11 @@ def _check_line_size(self, data_line, min_line_size):
             )
 
     def set_data(self, data, autofill=False, check_data=True):
-        """Sets the contents of the data to "data" with.  Data can have the
+        """Sets the contents of the data to "data".  Data can have the
         following formats:
             1) recarray - recarray containing the datalist
             2) [(line_one), (line_two), ...] - list where each line of the
                datalist is a tuple within the list
-            3) {'filename':filename, factor=fct, iprn=print_code, data=data}
-               - dictionary defining the external file containing the datalist.
         If the data is transient, a dictionary can be used to specify each
         stress period where the dictionary key is <stress period> - 1 and
         the dictionary value is the datalist data defined above:
@@ -563,6 +596,44 @@ def set_data(self, data, autofill=False, check_data=True):
         """
         self._set_data(data, autofill, check_data=check_data)
 
+    def set_record(self, data_record, autofill=False, check_data=True):
+        """Sets the contents of the data and metadata to "data_record".
+        Data_record is a dictionary with has the following format:
+            {'filename':filename, 'binary':True/False, 'data'=data}
+        To store to file include 'filename' in the dictionary.
+
+        Parameters
+        ----------
+            data_record : ndarray/list/dict
+                Data and metadata to set
+            autofill : bool
+                Automatically correct data
+            check_data : bool
+                Whether to verify the data
+
+        """
+        self._set_record(data_record, autofill, check_data)
+
+    def _set_record(self, data_record, autofill=False, check_data=True):
+        """Sets the contents of the data and metadata to "data_record".
+        Data_record is a dictionary with has the following format:
+            {'filename':filename, 'data'=data}
+        To store to file include 'filename' in the dictionary.
+
+        Parameters
+        ----------
+            data_record : ndarray/list/dict
+                Data and metadata to set
+            autofill : bool
+                Automatically correct data
+            check_data : bool
+                Whether to verify the data
+
+        """
+        self._set_data(
+            data_record, autofill, check_data=check_data, preserve_record=False
+        )
+
     def append_data(self, data):
         """Appends "data" to the end of this list.  Assumes data is in a format
         that can be appended directly to a numpy recarray.
@@ -1587,10 +1658,10 @@ def store_as_external_file(
         for sp in self._data_storage.keys():
             self._current_key = sp
             layer_storage = self._get_storage_obj().layer_storage
-            if (
-                layer_storage.get_total_size() > 0
-                and self._get_storage_obj().layer_storage[0].data_storage_type
+            if layer_storage.get_total_size() > 0 and (
+                self._get_storage_obj().layer_storage[0].data_storage_type
                 != DataStorageType.external_file
+                or replace_existing_external
             ):
                 fname, ext = os.path.splitext(external_file_path)
                 if datautil.DatumUtil.is_int(sp):
@@ -1643,6 +1714,34 @@ def has_data(self, key=None):
             self.get_data_prep(key)
             return super().has_data()
 
+    def get_record(self, key=None):
+        """Returns the data for stress period `key`.  If no key is specified
+        returns all records in a dictionary with zero-based stress period
+        numbers as keys.  See MFList's get_record documentation for more
+        information on the format of each record returned.
+
+        Parameters
+        ----------
+            key : int
+                Zero-based stress period to return data from.
+
+        Returns
+        -------
+            data_record : dict
+
+        """
+        if self._data_storage is not None and len(self._data_storage) > 0:
+            if key is None:
+                output = {}
+                for key in self._data_storage.keys():
+                    self.get_data_prep(key)
+                    output[key] = super().get_record()
+                return output
+            self.get_data_prep(key)
+            return super().get_record()
+        else:
+            return None
+
     def get_data(self, key=None, apply_mult=False, **kwargs):
         """Returns the data for stress period `key`.
 
@@ -1686,6 +1785,29 @@ def get_data(self, key=None, apply_mult=False, **kwargs):
         else:
             return None
 
+    def set_record(self, data_record, autofill=False, check_data=True):
+        """Sets the contents of the data based on the contents of
+        'data_record`.
+
+        Parameters
+        ----------
+        data_record : dict
+            Data_record being set.  Data_record must be a dictionary with
+            keys as zero-based stress periods and values as dictionaries
+            containing the data and metadata.  See MFList's set_record
+            documentation for more information on the format of the values.
+        autofill : bool
+            Automatically correct data
+        check_data : bool
+            Whether to verify the data
+        """
+        self._set_data_record(
+            data_record,
+            autofill=autofill,
+            check_data=check_data,
+            is_record=True,
+        )
+
     def set_data(self, data, key=None, autofill=False):
         """Sets the contents of the data at time `key` to `data`.
 
@@ -1694,7 +1816,7 @@ def set_data(self, data, key=None, autofill=False):
         data : dict, recarray, list
             Data being set.  Data can be a dictionary with keys as
             zero-based stress periods and values as the data.  If data is
-            an recarray or list of tuples, it will be assigned to the the
+            a recarray or list of tuples, it will be assigned to the
             stress period specified in `key`.  If any is set to None, that
             stress period of data will be removed.
         key : int
@@ -1703,9 +1825,14 @@ def set_data(self, data, key=None, autofill=False):
         autofill : bool
             Automatically correct data.
         """
+        self._set_data_record(data, key, autofill)
+
+    def _set_data_record(
+        self, data, key=None, autofill=False, check_data=False, is_record=False
+    ):
         self._cache_model_grid = True
         if isinstance(data, dict):
-            if "filename" not in data:
+            if "filename" not in data and "data" not in data:
                 # each item in the dictionary is a list for one stress period
                 # the dictionary key is the stress period the list is for
                 del_keys = []
@@ -1726,9 +1853,12 @@ def set_data(self, data, key=None, autofill=False):
                         else:
                             check = True
                         self._set_data_prep(list_item, key)
-                        super().set_data(
-                            list_item, autofill=autofill, check_data=check
-                        )
+                        if is_record:
+                            super().set_record(list_item, autofill, check_data)
+                        else:
+                            super().set_data(
+                                list_item, autofill=autofill, check_data=check
+                            )
                 for key in del_keys:
                     del data[key]
             else:
@@ -1736,6 +1866,25 @@ def set_data(self, data, key=None, autofill=False):
                 self._set_data_prep(data["data"], key)
                 super().set_data(data, autofill)
         else:
+            if is_record:
+                comment = (
+                    "Set record method requires that data_record is a "
+                    "dictionary."
+                )
+                type_, value_, traceback_ = sys.exc_info()
+                raise MFDataException(
+                    self.structure.get_model(),
+                    self.structure.get_package(),
+                    self._path,
+                    "setting data record",
+                    self.structure.name,
+                    inspect.stack()[0][3],
+                    type_,
+                    value_,
+                    traceback_,
+                    comment,
+                    self._simulation_data.debug,
+                )
             if key is None:
                 # search for a key
                 new_key_index = self.structure.first_non_keyword_index()
diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py
index 33c495263..20dec5f53 100644
--- a/flopy/mf6/data/mfdatastorage.py
+++ b/flopy/mf6/data/mfdatastorage.py
@@ -313,8 +313,6 @@ def __init__(
         self.layer_storage = MultiList(
             shape=layer_shape, callback=self._create_layer
         )
-        # self.layer_storage = [LayerStorage(self, x, data_storage_type)
-        #                      for x in range(layer_shape)]
         self.data_structure_type = data_structure_type
         package_dim = self.data_dimensions.package_dim
         self.in_model = (
@@ -633,6 +631,35 @@ def get_const_val(self, layer=None):
                 )
             return self.layer_storage[layer].get_data_const_val()
 
+    def get_record(self, layer=None):
+        if layer is None:
+            if self.layered:
+                record_dict = {}
+                for lay_num in self.layer_storage.indexes():
+                    record_dict[lay_num[0]] = self._get_record_layer(lay_num)
+                return record_dict
+            else:
+                return self._get_record_layer(0)
+        else:
+            return self._get_record_layer(layer)
+
+    def _get_record_layer(self, lay_num):
+        if (
+            self.layer_storage[lay_num].data_storage_type
+            == DataStorageType.external_file
+        ):
+            layer_dict = {
+                "filename": self.layer_storage[lay_num].fname,
+                "binary": self.layer_storage[lay_num].binary,
+            }
+        else:
+            layer_dict = {}
+        if self.data_structure_type == DataStructureType.ndarray:
+            layer_dict["factor"] = self.layer_storage[lay_num].factor
+            layer_dict["iprn"] = self.layer_storage[lay_num].iprn
+        layer_dict["data"] = self.get_data(lay_num, False, False)
+        return layer_dict
+
     def has_data(self, layer=None):
         ret_val = self._access_data(layer, False)
         return ret_val is not None and ret_val is not False
@@ -680,7 +707,7 @@ def _access_data(self, layer, return_data=False, apply_mult=True):
             == DataStorageType.external_file
         ):
             if return_data:
-                return self.external_to_internal(layer)
+                return self.external_to_internal(layer, apply_mult=apply_mult)
             else:
                 return True
         else:
@@ -871,19 +898,35 @@ def set_data(
         key=None,
         autofill=False,
         check_data=False,
+        preserve_record=False,
     ):
-        if multiplier is None:
-            multiplier = [1.0]
         if (
             self.data_structure_type == DataStructureType.recarray
             or self.data_structure_type == DataStructureType.scalar
         ):
-            self._set_list(data, layer, multiplier, key, autofill, check_data)
+            self._set_list(
+                data,
+                layer,
+                multiplier,
+                key,
+                autofill,
+                check_data,
+                preserve_record,
+            )
         else:
-            self._set_array(data, layer, multiplier, key, autofill)
+            self._set_array(
+                data, layer, multiplier, key, autofill, preserve_record
+            )
 
     def _set_list(
-        self, data, layer, multiplier, key, autofill, check_data=False
+        self,
+        data,
+        layer,
+        multiplier,
+        key,
+        autofill,
+        check_data=False,
+        preserve_record=False,
     ):
         if isinstance(data, dict):
             if "filename" in data:
@@ -918,6 +961,21 @@ def _set_list(
             ):
                 # single line of data needs to be encapsulated in a tuple
                 data = [tuple(data)]
+        if preserve_record:
+            if (
+                self.layer_storage[0].data_storage_type
+                == DataStorageType.external_file
+                and self.layer_storage[0].fname
+            ):
+                # build dictionary with current record data and
+                # store externally
+                data_record = {
+                    "filename": self.layer_storage[0].fname,
+                    "binary": self.layer_storage[0].binary,
+                    "data": data,
+                }
+                self.process_open_close_line(data_record, layer)
+                return
         self.store_internal(
             data,
             layer,
@@ -926,9 +984,12 @@ def _set_list(
             key=key,
             autofill=autofill,
             check_data=check_data,
+            preserve_record=preserve_record,
         )
 
-    def _set_array(self, data, layer, multiplier, key, autofill):
+    def _set_array(
+        self, data, layer, multiplier, key, autofill, preserve_record=False
+    ):
         # make a list out of a single item
         if (
             isinstance(data, int)
@@ -937,8 +998,40 @@ def _set_array(self, data, layer, multiplier, key, autofill):
         ):
             data = [data]
 
-        # check for possibility of multi-layered data
         success = False
+        if preserve_record:
+            if isinstance(data, np.ndarray):
+                # try to store while preserving the structure of the
+                # existing record
+                if self.layer_storage.get_total_size() > 1:
+                    if len(data) == self.layer_storage.get_total_size():
+                        # break ndarray into layers and store
+                        success = self._set_array_by_layer(
+                            data, multiplier, key, preserve_record
+                        )
+                else:
+                    # try to store as a single layer
+                    success = self._set_array_layer(
+                        data, layer, multiplier, key, preserve_record
+                    )
+            elif isinstance(data, dict):
+                first_key = list(data.keys())[0]
+                if isinstance(first_key, int):
+                    for layer_num, data_layer in data.items():
+                        success = self._set_array_layer(
+                            data_layer,
+                            layer_num,
+                            multiplier,
+                            key,
+                            preserve_record,
+                        )
+
+        if not success:
+            # storing while preserving the record failed, try storing as a
+            # new record
+            preserve_record = False
+
+        # check for possibility of multi-layered data
         layer_num = 0
         if (
             layer is None
@@ -946,22 +1039,14 @@ def _set_array(self, data, layer, multiplier, key, autofill):
             and len(data) == self.layer_storage.get_total_size()
             and not isinstance(data, dict)
         ):
-            # loop through list and try to store each list entry as a layer
-            success = True
-            for layer_num, layer_data in enumerate(data):
-                if (
-                    not isinstance(layer_data, list)
-                    and not isinstance(layer_data, dict)
-                    and not isinstance(layer_data, np.ndarray)
-                ):
-                    layer_data = [layer_data]
-                layer_index = self.layer_storage.nth_index(layer_num)
-                success = success and self._set_array_layer(
-                    layer_data, layer_index, multiplier, key
-                )
+            success = self._set_array_by_layer(
+                data, multiplier, key, preserve_record
+            )
         if not success:
             # try to store as a single layer
-            success = self._set_array_layer(data, layer, multiplier, key)
+            success = self._set_array_layer(
+                data, layer, multiplier, key, preserve_record
+            )
         self.layered = bool(self.layer_storage.get_total_size() > 1)
         if not success:
             message = (
@@ -984,11 +1069,58 @@ def _set_array(self, data, layer, multiplier, key, autofill):
                 self._simulation_data.debug,
             )
 
-    def _set_array_layer(self, data, layer, multiplier, key):
+    def _set_array_by_layer(self, data, multiplier, key, preserve_record):
+        # loop through list and try to store each list entry as a layer
+        success = True
+        for layer_num, layer_data in enumerate(data):
+            if (
+                not isinstance(layer_data, list)
+                and not isinstance(layer_data, dict)
+                and not isinstance(layer_data, np.ndarray)
+            ):
+                layer_data = [layer_data]
+            layer_index = self.layer_storage.nth_index(layer_num)
+            success = success and self._set_array_layer(
+                layer_data, layer_index, multiplier, key, preserve_record
+            )
+        return success
+
+    def _set_array_layer(self, data, layer, multiplier, key, preserve_record):
         # look for a single constant value
         data_type = self.data_dimensions.structure.get_datum_type(
             return_enum_type=True
         )
+
+        if isinstance(data, np.ndarray) and preserve_record:
+            # store data and preserve record
+            layer = self._layer_prep(layer)
+            if not self.layer_storage.in_shape(layer):
+                return False
+            if (
+                self.layer_storage[layer].data_storage_type
+                == DataStorageType.external_file
+            ):
+                self.store_external(
+                    self.layer_storage[layer].fname,
+                    layer,
+                    self.layer_storage[layer].factor,
+                    self.layer_storage[layer].iprn,
+                    data=data,
+                    do_not_verify=True,
+                    binary=self.layer_storage[layer].binary,
+                    preserve_record=preserve_record,
+                )
+            else:
+                self.store_internal(
+                    data,
+                    layer,
+                    False,
+                    multiplier,
+                    key=key,
+                    preserve_record=preserve_record,
+                )
+            return True
+
         if not isinstance(data, dict) and not isinstance(data, str):
             if self._calc_data_size(data, 2) == 1 and self._is_type(
                 data[0], data_type
@@ -999,6 +1131,7 @@ def _set_array_layer(self, data, layer, multiplier, key):
 
         # look for internal and open/close data
         if isinstance(data, dict):
+            data_la = None
             if "data" in data:
                 if (
                     isinstance(data["data"], int)
@@ -1007,7 +1140,7 @@ def _set_array_layer(self, data, layer, multiplier, key):
                 ):
                     # data should always in in a list/array
                     data["data"] = [data["data"]]
-
+                data_la = data["data"]
             if "filename" in data:
                 multiplier, iprn, binary = self.process_open_close_line(
                     data, layer
@@ -1018,6 +1151,7 @@ def _set_array_layer(self, data, layer, multiplier, key):
                     layer,
                     [multiplier],
                     print_format=iprn,
+                    data=data_la,
                     binary=binary,
                     do_not_verify=True,
                 )
@@ -1124,9 +1258,14 @@ def store_internal(
         autofill=False,
         print_format=None,
         check_data=False,
+        preserve_record=False,
     ):
-        if multiplier is None:
-            multiplier = [self.get_default_mult()]
+        if multiplier is None and layer is not None:
+            layer = self._layer_prep(layer)
+            if self.layer_storage.in_shape(layer) and preserve_record:
+                multiplier = [self.layer_storage[layer].factor]
+            else:
+                multiplier = [self.get_default_mult()]
         if self.data_structure_type == DataStructureType.recarray:
             if (
                 self.layer_storage.first_item().data_storage_type
@@ -1222,6 +1361,12 @@ def store_internal(
         else:
             layer, multiplier = self._store_prep(layer, multiplier)
             dimensions = self.get_data_dimensions(layer)
+            if preserve_record:
+                factor = self.layer_storage[layer].factor
+                adjustment = multiplier / factor
+                if adjustment != 1.0:
+                    # convert numbers to be multiplied by the original factor
+                    data = data * adjustment
             if const:
                 self.layer_storage[
                     layer
@@ -1258,8 +1403,9 @@ def store_internal(
                         message,
                         self._simulation_data.debug,
                     )
-            self.layer_storage[layer].factor = multiplier
-            self.layer_storage[layer].iprn = print_format
+            if not preserve_record:
+                self.layer_storage[layer].factor = multiplier
+                self.layer_storage[layer].iprn = print_format
 
     def _resolve_data_line(self, data, key):
         if len(self._recarray_type_list) > 1:
@@ -1502,16 +1648,21 @@ def store_external(
         data=None,
         do_not_verify=False,
         binary=False,
+        preserve_record=False,
     ):
-        if multiplier is None:
-            multiplier = [self.get_default_mult()]
+        if multiplier is None and layer is not None:
+            layer = self._layer_prep(layer)
+            if self.layer_storage.in_shape(layer) and preserve_record:
+                multiplier = [self.layer_storage[layer].factor]
+            else:
+                multiplier = [self.get_default_mult()]
         layer_new, multiplier = self._store_prep(layer, multiplier)
 
         # pathing to external file
         data_dim = self.data_dimensions
         model_name = data_dim.package_dim.model_dim[0].model_name
         fp_relative = file_path
-        if model_name is not None:
+        if model_name is not None and fp_relative is not None:
             rel_path = self._simulation_data.mfpath.model_relative_path[
                 model_name
             ]
@@ -1524,7 +1675,20 @@ def store_external(
                 for i, rp in enumerate(rp_l_r):
                     if rp != fp_rp_l[len(rp_l_r) - i - 1]:
                         fp_relative = os.path.join(rp, fp_relative)
-        fp = self._simulation_data.mfpath.resolve_path(fp_relative, model_name)
+            fp = self._simulation_data.mfpath.resolve_path(
+                fp_relative, model_name
+            )
+        else:
+            fp = os.path.join(
+                self._simulation_data.mfpath.get_sim_path(), fp_relative
+            )
+        if layer_new in self.layer_storage:
+            old_ext_file = self.layer_storage[layer_new].external
+            old_binary = self.layer_storage[layer_new].binary
+            if preserve_record and old_ext_file:
+                # use old file settings
+                fp = old_ext_file
+                binary = old_binary
         if data is not None:
             if self.data_structure_type == DataStructureType.recarray:
                 # create external file and write file entry to the file
@@ -1538,6 +1702,7 @@ def store_external(
                     None,
                     False,
                     print_format,
+                    preserve_record,
                 )
                 if binary:
                     file_access = MFFileAccessList(
@@ -1589,6 +1754,15 @@ def store_external(
                 # set as external data
                 self.layer_storage.first_item().internal_data = None
             else:
+                # if self.layer_storage.in_shape(layer_new):
+                #    factor = self.layer_storage[layer_new].factor
+                # if preserve_record:
+                #    adjustment = multiplier / factor
+                #    if adjustment != 1.0:
+                # convert numbers to be multiplied by the
+                # original factor
+                #        data = data * adjustment
+
                 # store data externally in file
                 data_size = self.get_data_size(layer_new)
                 data_type = data_dim.structure.data_item_structures[0].type
@@ -1636,10 +1810,10 @@ def store_external(
                         data_type,
                         data_size,
                     )
-                self.layer_storage[layer_new].factor = multiplier
+                if not preserve_record:
+                    self.layer_storage[layer_new].factor = multiplier
                 self.layer_storage[layer_new].internal_data = None
                 self.layer_storage[layer_new].data_const_value = None
-
         else:
             if self.data_structure_type == DataStructureType.recarray:
                 self.layer_storage.first_item().internal_data = None
@@ -1753,7 +1927,9 @@ def external_to_external(
             binary=binary,
         )
 
-    def external_to_internal(self, layer, store_internal=False):
+    def external_to_internal(
+        self, layer, store_internal=False, apply_mult=True
+    ):
         # reset comments
         self.pre_data_comments = None
         self.comments = {}
@@ -1790,7 +1966,7 @@ def external_to_internal(self, layer, store_internal=False):
                     layer,
                     read_file,
                 )[0]
-            if self.layer_storage[layer].factor is not None:
+            if apply_mult and self.layer_storage[layer].factor is not None:
                 data_out = data_out * self.layer_storage[layer].factor
 
             if store_internal:
@@ -2865,6 +3041,19 @@ def _has_layer_dim(self):
             or "nodes" in self.data_dimensions.structure.shape
         )
 
+    def _layer_prep(self, layer):
+        if layer is None:
+            # layer is none means the data provided is for all layers or this
+            # is not layered data
+            layer = (0,)
+            self.layer_storage.list_shape = (1,)
+            self.layer_storage.multi_dim_list = [
+                self.layer_storage.first_item()
+            ]
+        if isinstance(layer, int):
+            layer = (layer,)
+        return layer
+
     def _store_prep(self, layer, multiplier):
         if not (layer is None or self.layer_storage.in_shape(layer)):
             message = f"Layer {layer} is not a valid layer."
@@ -2882,25 +3071,23 @@ def _store_prep(self, layer, multiplier):
                 message,
                 self._simulation_data.debug,
             )
-        if layer is None:
-            # layer is none means the data provided is for all layers or this
-            # is not layered data
-            layer = (0,)
-            self.layer_storage.list_shape = (1,)
-            self.layer_storage.multi_dim_list = [
-                self.layer_storage.first_item()
-            ]
-        mult_ml = MultiList(multiplier)
-        if not mult_ml.in_shape(layer):
-            if multiplier[0] is None:
-                multiplier = self.get_default_mult()
-            else:
-                multiplier = multiplier[0]
+        layer = self._layer_prep(layer)
+        if multiplier is None:
+            multiplier = self.get_default_mult()
         else:
-            if mult_ml.first_item() is None:
-                multiplier = self.get_default_mult()
+            if isinstance(multiplier, float):
+                multiplier = [multiplier]
+            mult_ml = MultiList(multiplier)
+            if not mult_ml.in_shape(layer):
+                if multiplier[0] is None:
+                    multiplier = self.get_default_mult()
+                else:
+                    multiplier = multiplier[0]
             else:
-                multiplier = mult_ml.first_item()
+                if mult_ml.first_item() is None:
+                    multiplier = self.get_default_mult()
+                else:
+                    multiplier = mult_ml.first_item()
 
         return layer, multiplier
 
diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py
index cce00c691..b00c30ef5 100644
--- a/flopy/mf6/mfpackage.py
+++ b/flopy/mf6/mfpackage.py
@@ -1299,6 +1299,7 @@ def set_all_data_external(
                 and dataset.enabled
             ):
                 file_path = f"{base_name}_{dataset.structure.name}.txt"
+                replace_existing_external = False
                 if external_data_folder is not None:
                     # get simulation root path
                     root_path = self._simulation_data.mfpath.get_sim_path()
@@ -1317,9 +1318,10 @@ def set_all_data_external(
                         # create new external data folder
                         os.makedirs(full_path)
                     file_path = os.path.join(external_data_folder, file_path)
+                    replace_existing_external = True
                 dataset.store_as_external_file(
                     file_path,
-                    replace_existing_external=False,
+                    replace_existing_external=replace_existing_external,
                     check_data=check_data,
                 )
 
diff --git a/flopy/utils/datautil.py b/flopy/utils/datautil.py
index c8bada4bb..e57a5161d 100644
--- a/flopy/utils/datautil.py
+++ b/flopy/utils/datautil.py
@@ -582,6 +582,8 @@ def get_total_size(self):
         return shape_size
 
     def in_shape(self, indexes):
+        if isinstance(indexes, int):
+            indexes = [indexes]
         for index, item in zip(indexes, self.list_shape):
             if index > item:
                 return False