diff --git a/docs/source/examples/verify_mass_conservation.py b/docs/source/examples/verify_mass_conservation.py
index 78bf01d1..98223a49 100644
--- a/docs/source/examples/verify_mass_conservation.py
+++ b/docs/source/examples/verify_mass_conservation.py
@@ -72,7 +72,11 @@ def _verify_mass_conservation(self):
         total amount of landscape change (volume) is close to the total
         amount of sediment put into the domain.
         """
-        # calculate change in input sediment and deposit volume
+
+        # check the change in volume of the deposit for this timestep
+        act_deposit_volume = np.sum((self.eta - self.eta_init) * self.dx**2) * (
+            1 - self.porosity
+        )
         input_sed_volume = self.Qs0 * self.dt
         act_deposit_volume = np.sum((self.eta - self.eta_init) * self.dx**2)
 
diff --git a/pyDeltaRCM/_version.py b/pyDeltaRCM/_version.py
index d47030b1..93ae2154 100644
--- a/pyDeltaRCM/_version.py
+++ b/pyDeltaRCM/_version.py
@@ -3,4 +3,4 @@ def __version__() -> str:
     Private version declaration, gets assigned to pyDeltaRCM.__version__
     during import
     """
-    return "2.3.0"
+    return "2.3.1"
diff --git a/pyDeltaRCM/default.yml b/pyDeltaRCM/default.yml
index 942a940f..ff4b8595 100644
--- a/pyDeltaRCM/default.yml
+++ b/pyDeltaRCM/default.yml
@@ -178,6 +178,9 @@ stepmax:
 force_deposit:
   type: ['bool', int]
   default: False
+porosity:
+  type: ['float', 'int']
+  default: 0
 sand_frac_bc:
   type: ['float', 'int']
   default: 0
diff --git a/pyDeltaRCM/init_tools.py b/pyDeltaRCM/init_tools.py
index b577da5d..6f3370f9 100644
--- a/pyDeltaRCM/init_tools.py
+++ b/pyDeltaRCM/init_tools.py
@@ -671,6 +671,7 @@ def init_sediment_routers(self) -> None:
             self.stepmax,
             self.force_deposit,
             self.theta_mud,
+            self.porosity,
             self.mod_erosion,
         )
         # initialize the SandRouter object
@@ -693,6 +694,7 @@ def init_sediment_routers(self) -> None:
             self.stepmax,
             self.force_deposit,
             self.theta_sand,
+            self.porosity,
             self.mod_erosion,
         )
 
diff --git a/pyDeltaRCM/model.py b/pyDeltaRCM/model.py
index b67d6a80..db2e05ea 100644
--- a/pyDeltaRCM/model.py
+++ b/pyDeltaRCM/model.py
@@ -1351,6 +1351,21 @@ def force_deposit(self) -> bool:
     def force_deposit(self, force_deposit: bool) -> None:
         self._force_deposit = force_deposit
 
+    @property
+    def porosity(self) -> float:
+        """
+        `porosity` of deposited sediment.
+
+        Default is 0, no porosity. Must be between 0 and 1.
+        """
+        return self._porosity
+
+    @porosity.setter
+    def porosity(self, porosity: float) -> None:
+        if (porosity < 0) or (porosity > 1):
+            raise ValueError("Value for porosity must be between 0 and 1, inclusive.")
+        self._porosity = porosity
+
     @property
     def clobber_netcdf(self) -> bool:
         """
@@ -1454,7 +1469,9 @@ def time_step(self) -> float:
     @time_step.setter
     def time_step(self, new_time_step: float) -> None:
         if new_time_step * self.init_Np_sed < 100:
-            _msg = "Using a very small time step, so the delta might evolve very slowly."
+            _msg = (
+                "Using a very small time step, so the delta might evolve very slowly."
+            )
             self.log_warning(_msg)
             warnings.warn(UserWarning(_msg))
 
diff --git a/pyDeltaRCM/sed_tools.py b/pyDeltaRCM/sed_tools.py
index 0c348cb7..aba5aac8 100644
--- a/pyDeltaRCM/sed_tools.py
+++ b/pyDeltaRCM/sed_tools.py
@@ -341,6 +341,7 @@ def _get_weight_at_cell_sediment(
     ("u_max", float32),
     ("qs0", float32),
     ("_u0", float32),
+    ("porosity", float32),
     ("Vp_sed", float32),
     ("Vp_res", float32),
     ("Vp_dep_mud", float32[:, :]),
@@ -474,7 +475,7 @@ def _update_fields(self, Vp_change: float, px: int, py: int) -> None:
         #    determinations require several comparisons and repeated indexing,
         #    so we use a jitted "helper" function to do the operations.
         qw0 = self.qw[px, py]
-        eta_change = Vp_change / (self._dx * self._dx)
+        eta_change = Vp_change / (self._dx * self._dx) / (1 - self.porosity)
 
         eta = self.eta[px, py] + eta_change  # new bed
         depth = self.stage[px, py] - eta  # new depth
@@ -594,6 +595,7 @@ def __init__(
         stepmax,
         force_deposit,
         theta_sed: float,
+        porosity,
         mod_erosion,
     ) -> None:
         self._dt = _dt
@@ -621,6 +623,7 @@ def __init__(
         self.stepmax = stepmax
         self.force_deposit = force_deposit
         self.theta_sed = theta_sed
+        self.porosity = porosity
         self.mod_erosion = mod_erosion
 
     def run(
@@ -878,6 +881,7 @@ def __init__(
         stepmax,
         force_deposit,
         theta_sed: float,
+        porosity,
         mod_erosion,
     ) -> None:
         self._dt = _dt
@@ -900,6 +904,7 @@ def __init__(
         self.stepmax = stepmax
         self.force_deposit = force_deposit
         self.theta_sed = theta_sed
+        self.porosity = porosity
         self.mod_erosion = mod_erosion
 
     def run(