diff --git a/autotest/t076_modelgrid_thick.py b/autotest/t076_modelgrid_thick.py
index 44fd300cb..74d994b17 100644
--- a/autotest/t076_modelgrid_thick.py
+++ b/autotest/t076_modelgrid_thick.py
@@ -199,7 +199,53 @@ def test_unstructured_thick():
     return
 
 
+def test_ncb_thick():
+    nlay = 3
+    nrow = ncol = 15
+    laycbd = np.array([1, 2, 0], dtype=int)
+    ncb = np.count_nonzero(laycbd)
+    dx = dy = 150
+    delc = np.array([dy,] * nrow)
+    delr = np.array([dx, ] * ncol)
+    top = np.ones((15, 15))
+    botm = np.ones((nlay + ncb, nrow, ncol))
+    elevations = np.array([-10, -20, -40, -50, -70])[:, np.newaxis]
+    botm *= elevations[:, None]
+
+    modelgrid = StructuredGrid(
+        delc=delc,
+        delr=delr,
+        top=top,
+        botm=botm,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        laycbd=laycbd
+    )
+
+    thick = modelgrid.thick
+
+    if not thick.shape[0] == nlay + ncb:
+        raise AssertionError("grid thick attribute returns incorrect shape")
+
+    thick = modelgrid.remove_confining_beds(modelgrid.thick)
+    if not thick.shape == modelgrid.shape:
+        raise AssertionError("quasi3d confining beds not properly removed")
+
+    sat_thick = modelgrid.saturated_thick(modelgrid.thick)
+    if not sat_thick.shape == modelgrid.shape:
+        raise AssertionError(
+            "saturated_thickness confining beds not removed"
+        )
+
+    if sat_thick[1, 0, 0] != 20:
+        raise AssertionError(
+            "saturated_thickness is not properly indexing confining beds"
+        )
+
+
 if __name__ == "__main__":
     test_unstructured_thick()
     test_vertices_thick()
     test_structured_thick()
+    test_ncb_thick()
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 7b8e587a0..0fd809570 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -325,6 +325,11 @@ def saturated_thick(self, array, mask=None):
         thick = self.thick
         top = self.top_botm[:-1].reshape(thick.shape)
         bot = self.top_botm[1:].reshape(thick.shape)
+        thick = self.remove_confining_beds(thick)
+        top = self.remove_confining_beds(top)
+        bot = self.remove_confining_beds(bot)
+        array = self.remove_confining_beds(array)
+
         idx = np.where((array < top) & (array > bot))
         thick[idx] = array[idx] - bot[idx]
         idx = np.where(array <= bot)
@@ -451,6 +456,32 @@ def xyzvertices(self):
     def cross_section_vertices(self):
         return self.xyzvertices[0], self.xyzvertices[1]
 
+    def remove_confining_beds(self, array):
+        """
+        Method to remove confining bed layers from an array
+
+        Parameters
+        ----------
+        array : np.ndarray
+            array to remove quasi3d confining bed data from. Shape of axis 0
+            should be (self.lay + ncb) to remove beds
+        Returns
+        -------
+            np.ndarray
+        """
+        if self.laycbd is not None:
+            ncb = np.count_nonzero(self.laycbd)
+            if ncb > 0:
+                if array.shape[0] == self.shape[0] + ncb:
+                    cb = 0
+                    idx = []
+                    for ix, i in enumerate(self.laycbd):
+                        idx.append(ix + cb)
+                        if i > 0:
+                            cb += 1
+                    array = array[idx]
+        return array
+
     def cross_section_lay_ncpl_ncb(self, ncb):
         """
         Get PlotCrossSection compatible layers, ncpl, and ncb
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index 7aed76ee2..14ca47a9c 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -161,7 +161,7 @@ def __init__(
                 self.__nlay = nlay
             else:
                 if laycbd is not None:
-                    self.__nlay = len(botm) - np.sum(laycbd > 0)
+                    self.__nlay = len(botm) - np.count_nonzero(laycbd)
                 else:
                     self.__nlay = len(botm)
         else:
diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py
index 171aceb54..633cbe753 100644
--- a/flopy/modflow/mf.py
+++ b/flopy/modflow/mf.py
@@ -331,7 +331,7 @@ def modelgrid(self):
                 yoff=self._modelgrid.yoffset,
                 angrot=self._modelgrid.angrot,
                 nlay=self.dis.nlay,
-                laycbd=self.dis.laycbd,
+                laycbd=self.dis.laycbd.array,
             )
 
         # resolve offsets
diff --git a/flopy/utils/postprocessing.py b/flopy/utils/postprocessing.py
index e9a9b5c67..f97214f3f 100644
--- a/flopy/utils/postprocessing.py
+++ b/flopy/utils/postprocessing.py
@@ -707,12 +707,13 @@ def get_specific_discharge(
     if classical_budget:
         # get saturated thickness (head - bottom elev for unconfined layer)
         if head is None:
-            sat_thk = modelgrid.thick
+            sat_thk = modelgrid.remove_confining_beds(modelgrid.thick)
         else:
             sat_thk = modelgrid.saturated_thick(
                 head, mask=[model.hdry, model.hnoflo]
             )
-            sat_thk.shape = model.modelgrid.shape
+
+        sat_thk.shape = modelgrid.shape
 
         # inform modelgrid of no-flow and dry cells
         modelgrid = model.modelgrid