diff --git a/flopy/modpath/mp7bas.py b/flopy/modpath/mp7bas.py
index e9f8344aa..2480b7985 100644
--- a/flopy/modpath/mp7bas.py
+++ b/flopy/modpath/mp7bas.py
@@ -55,24 +55,14 @@ def __init__(
 
         self._generate_heading()
 
-        if model.flowmodel.version == "mf6":
-            self.laytyp = Util2d(
-                self.parent,
-                (shape[0],),
-                np.int32,
-                model.laytyp,
-                name="bas - laytype",
-                locat=self.unit_number[0],
-            )
-        else:
-            self.laytyp = Util2d(
-                self.parent,
-                (shape[0],),
-                np.int32,
-                model.laytyp,
-                name="bas - laytype",
-                locat=self.unit_number[0],
-            )
+        self.laytyp = Util2d(
+            self.parent,
+            (shape[0],),
+            np.int32,
+            model.laytyp,
+            name="bas - laytype",
+            locat=self.unit_number[0],
+        )
         if model.flowmodel.version != "mf6":
             self.ibound = Util3d(
                 model,
diff --git a/flopy/modpath/mp7sim.py b/flopy/modpath/mp7sim.py
index 774bbe11a..96487c9bb 100644
--- a/flopy/modpath/mp7sim.py
+++ b/flopy/modpath/mp7sim.py
@@ -498,6 +498,14 @@ def __init__(
         self.timepointoption = timepointoption
         self.timepointdata = timepointdata
 
+        shape = self.parent.shape
+        if len(shape) == 3:
+            shape3d = shape
+        elif len(shape) == 2:
+            shape3d = (shape[0], 1, shape[1])
+        else:
+            shape3d = (1, 1, shape[0])
+
         # zonedataoption
         try:
             self.zonedataoption = onoffOpt[zonedataoption.lower()].value
@@ -517,7 +525,7 @@ def __init__(
                 )
             self.zones = Util3d(
                 model,
-                self.parent.shape,
+                shape3d,
                 np.int32,
                 zones,
                 name="zones",
@@ -541,7 +549,7 @@ def __init__(
                 )
             self.retardation = Util3d(
                 model,
-                self.parent.shape,
+                shape3d,
                 np.float32,
                 retardation,
                 name="retardation",