diff --git a/autotest/t501_test.py b/autotest/t501_test.py
index 773fb6c6cc..aca0b65aac 100644
--- a/autotest/t501_test.py
+++ b/autotest/t501_test.py
@@ -165,6 +165,45 @@ def test_mf6_string_to_file_path():
             raise AssertionError("Relative path error")
 
 
+def test_mf6_subdir():
+    base_dir = base_test_dir(__file__, rel_path="temp", verbose=True)
+    test_setup = FlopyTestSetup(verbose=True, test_dirs=base_dir)
+
+    sim = flopy.mf6.MFSimulation(sim_ws=base_dir)
+    tdis = flopy.mf6.modflow.mftdis.ModflowTdis(sim)
+    gwf = flopy.mf6.ModflowGwf(sim, model_rel_path="level2")
+    ims = flopy.mf6.modflow.mfims.ModflowIms(sim)
+    sim.register_ims_package(ims, [])
+    dis = flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(gwf)
+    sim.set_all_data_external(external_data_folder="dat")
+    sim.write_simulation()
+
+    sim_r = flopy.mf6.MFSimulation.load(
+        "mfsim.nam",
+        sim_ws=sim.simulation_data.mfpath.get_sim_path(),
+    )
+    gwf_r = sim_r.get_model()
+    assert (
+        gwf.dis.delc.get_file_entry() == gwf_r.dis.delc.get_file_entry()
+    ), "Something wrong with model external paths"
+
+    sim_r.set_all_data_internal()
+    sim_r.set_all_data_external(
+        external_data_folder=os.path.join("dat", "dat_l2")
+    )
+    sim_r.write_simulation()
+
+    sim_r2 = flopy.mf6.MFSimulation.load(
+        "mfsim.nam",
+        sim_ws=sim_r.simulation_data.mfpath.get_sim_path(),
+    )
+    gwf_r2 = sim_r.get_model()
+    assert (
+        gwf_r.dis.delc.get_file_entry() == gwf_r2.dis.delc.get_file_entry()
+    ), "Something wrong with model external paths"
+
+
 if __name__ == "__main__":
-    test_mf6()
-    test_mf6_string_to_file_path()
+    # test_mf6()
+    # test_mf6_string_to_file_path()
+    test_mf6_subdir()
diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py
index c59cb3de5b..cd0627dfae 100644
--- a/flopy/mf6/data/mfdatastorage.py
+++ b/flopy/mf6/data/mfdatastorage.py
@@ -1511,7 +1511,13 @@ def store_external(
             ]
             if rel_path is not None and len(rel_path) > 0 and rel_path != ".":
                 # include model relative path in external file path
-                fp_relative = os.path.join(rel_path, file_path)
+                # only if model relative path is not already in external
+                #  file path i.e. when reading!
+                fp_rp_l = fp_relative.split(os.path.sep)
+                rp_l_r = rel_path.split(os.path.sep)[::-1]
+                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)
         if data is not None:
             if self.data_structure_type == DataStructureType.recarray: