diff --git a/autotest/test_mp7.py b/autotest/test_mp7.py
index f4604748c..edb82c48d 100644
--- a/autotest/test_mp7.py
+++ b/autotest/test_mp7.py
@@ -855,6 +855,54 @@ def test_mp7sim_replacement(function_tmpdir):
     assert success, buff
 
 
+@requires_exe("mf6", "mp7")
+@pytest.mark.parametrize(
+    "porosity_type", ("constant", "list", "array_1d", "array")
+)
+@pytest.mark.slow
+def test_mp7bas_porosity(ex01_mf6_model, porosity_type):
+    sim, function_tmpdir = ex01_mf6_model
+    sim.run_simulation()
+
+    mpnam = f"{ex01_mf6_model_name}_mp_forward"
+    gwf = sim.get_model(ex01_mf6_model_name)
+
+    if porosity_type == "constant":
+        porosity = 0.30
+    elif porosity_type == "list":
+        porosity = [0.35, 0.30, 0.25]  # constant per layer
+    elif porosity_type == "array_1d":
+        porosity = np.array([0.35, 0.30, 0.25])  # constant per layer
+    elif porosity_type == "array":
+        # Initialize porosity array based on dim of dis.botm
+        # lay_1=0.35, lay_2=0.30, lay_3=0.25
+        porosity = np.array([[[0.35]], [[0.30]], [[0.25]]]) * np.ones(
+            (gwf.dis.botm.array.shape)
+        )
+        # Diversify porosity values, adjust both halves of the model
+        porosity[:, :, : porosity.shape[2] // 2] -= 0.05
+        porosity[:, :, porosity.shape[2] // 2 :] += 0.05
+
+    mp = Modpath7.create_mp7(
+        modelname=mpnam,
+        trackdir="forward",
+        flowmodel=gwf,
+        exe_name="mp7",
+        model_ws=function_tmpdir,
+        rowcelldivisions=1,
+        columncelldivisions=1,
+        layercelldivisions=1,
+        porosity=porosity,
+    )
+
+    # Check mean of assigned porosity
+    assert np.isclose(np.mean(mp.get_package("MPBAS").porosity.array), 0.3)
+
+    mp.write_input()
+    success, _ = mp.run_model()
+    assert success, f"mp7 model ({mp.name}) did not run"
+
+
 @requires_exe("mf6", "mp7")
 def test_flopy_2223(function_tmpdir):
     mf6sim = Mp7Cases.mf6(function_tmpdir)
diff --git a/flopy/modpath/mp7.py b/flopy/modpath/mp7.py
index 75d4a0ce8..6eea1cf99 100644
--- a/flopy/modpath/mp7.py
+++ b/flopy/modpath/mp7.py
@@ -412,6 +412,7 @@ def create_mp7(
         rowcelldivisions=2,
         layercelldivisions=2,
         nodes=None,
+        porosity=0.30,
     ):
         """
         Create a default MODPATH 7 model using a passed flowmodel with
@@ -447,6 +448,8 @@ def create_mp7(
             direction (default is 2).
         nodes : int, list of ints, tuple of ints, or np.ndarray
             Nodes (zero-based) with particles. If  (default is node 0).
+        porosity: float or array of floats (nlay, nrow, ncol)
+            The porosity array (the default is 0.30).
 
         Returns
         -------
@@ -477,7 +480,7 @@ def create_mp7(
 
         # create MODPATH 7 basic file and add to the MODPATH 7
         # model instance (mp)
-        Modpath7Bas(mp, defaultiface=defaultiface)
+        Modpath7Bas(mp, porosity=porosity, defaultiface=defaultiface)
 
         # create particles
         if nodes is None: