diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py
index ad27bf31c8..4076e6685f 100644
--- a/flopy/plot/plotutil.py
+++ b/flopy/plot/plotutil.py
@@ -1554,98 +1554,6 @@ def saturated_thickness(head, top, botm, laytyp, mask_values=None):
 
         return sat_thk
 
-    @staticmethod
-    def centered_specific_discharge(Qx, Qy, Qz, delr, delc, sat_thk):
-        """
-        DEPRECATED. Use postprocessing.get_specific_discharge() instead.
-
-        Using the MODFLOW discharge, calculate the cell centered specific
-        discharge by dividing by the flow width and then averaging
-        to the cell center.
-
-        Parameters
-        ----------
-        Qx : numpy.ndarray
-            MODFLOW 'flow right face'
-        Qy : numpy.ndarray
-            MODFLOW 'flow front face'.  The sign on this array will be flipped
-            by this function so that the y axis is positive to north.
-        Qz : numpy.ndarray
-            MODFLOW 'flow lower face'.  The sign on this array will be
-            flipped by this function so that the z axis is positive
-            in the upward direction.
-        delr : numpy.ndarray
-            MODFLOW delr array
-        delc : numpy.ndarray
-            MODFLOW delc array
-        sat_thk : numpy.ndarray
-            Saturated thickness for each cell
-
-        Returns
-        -------
-        (qx, qy, qz) : tuple of numpy.ndarrays
-            Specific discharge arrays that have been interpolated to cell centers.
-
-        """
-        warnings.warn(
-            "centered_specific_discharge() has been deprecated and will be "
-            "removed in version 3.3.5. Use "
-            "postprocessing.get_specific_discharge() instead.",
-            DeprecationWarning,
-        )
-
-        qx = None
-        qy = None
-        qz = None
-
-        if Qx is not None:
-            nlay, nrow, ncol = Qx.shape
-            qx = np.zeros(Qx.shape, dtype=Qx.dtype)
-
-            for k in range(nlay):
-                for j in range(ncol - 1):
-                    area = (
-                        delc[:]
-                        * 0.5
-                        * (sat_thk[k, :, j] + sat_thk[k, :, j + 1])
-                    )
-                    idx = area > 0.0
-                    qx[k, idx, j] = Qx[k, idx, j] / area[idx]
-
-            qx[:, :, 1:] = 0.5 * (qx[:, :, 0 : ncol - 1] + qx[:, :, 1:ncol])
-            qx[:, :, 0] = 0.5 * qx[:, :, 0]
-
-        if Qy is not None:
-            nlay, nrow, ncol = Qy.shape
-            qy = np.zeros(Qy.shape, dtype=Qy.dtype)
-
-            for k in range(nlay):
-                for i in range(nrow - 1):
-                    area = (
-                        delr[:]
-                        * 0.5
-                        * (sat_thk[k, i, :] + sat_thk[k, i + 1, :])
-                    )
-                    idx = area > 0.0
-                    qy[k, i, idx] = Qy[k, i, idx] / area[idx]
-
-            qy[:, 1:, :] = 0.5 * (qy[:, 0 : nrow - 1, :] + qy[:, 1:nrow, :])
-            qy[:, 0, :] = 0.5 * qy[:, 0, :]
-            qy = -qy
-
-        if Qz is not None:
-            qz = np.zeros(Qz.shape, dtype=Qz.dtype)
-            dr = delr.reshape((1, delr.shape[0]))
-            dc = delc.reshape((delc.shape[0], 1))
-            area = dr * dc
-            for k in range(nlay):
-                qz[k, :, :] = Qz[k, :, :] / area[:, :]
-            qz[1:, :, :] = 0.5 * (qz[0 : nlay - 1, :, :] + qz[1:nlay, :, :])
-            qz[0, :, :] = 0.5 * qz[0, :, :]
-            qz = -qz
-
-        return (qx, qy, qz)
-
 
 class UnstructuredPlotUtilities:
     """
@@ -2270,182 +2178,6 @@ def plot_shapefile(
     return pc
 
 
-def cvfd_to_patch_collection(verts, iverts):
-    """
-    Create a patch collection from control volume vertices and incidence list
-
-    Parameters
-    ----------
-    verts : ndarray
-        2d array of x and y points.
-    iverts : list of lists
-        should be of len(ncells) with a list of vertex numbers for each cell
-
-    """
-    warnings.warn(
-        "cvfd_to_patch_collection is deprecated and will be removed in "
-        "version 3.3.5. Use PlotMapView for plotting",
-        DeprecationWarning,
-    )
-
-    from matplotlib.collections import PatchCollection
-    from matplotlib.patches import Polygon
-
-    ptchs = []
-    for ivertlist in iverts:
-        points = []
-        for iv in ivertlist:
-            points.append((verts[iv, 0], verts[iv, 1]))
-        # close the polygon, if necessary
-        if ivertlist[0] != ivertlist[-1]:
-            iv = ivertlist[0]
-            points.append((verts[iv, 0], verts[iv, 1]))
-        ptchs.append(Polygon(points))
-    pc = PatchCollection(ptchs)
-    return pc
-
-
-def plot_cvfd(
-    verts,
-    iverts,
-    ax=None,
-    layer=0,
-    cmap="Dark2",
-    edgecolor="scaled",
-    facecolor="scaled",
-    a=None,
-    masked_values=None,
-    **kwargs,
-):
-    """
-    Generic function for plotting a control volume finite difference grid of
-    information.
-
-    Parameters
-    ----------
-    verts : ndarray
-        2d array of x and y points.
-    iverts : list of lists
-        should be of len(ncells) with a list of vertex number for each cell
-    ax : matplotlib.pylot axis
-        matplotlib.pyplot axis instance. Default is None
-    layer : int
-        layer to extract. Used in combination to the optional ncpl
-        parameter. Default is 0
-    cmap : string
-        Name of colormap to use for polygon shading (default is 'Dark2')
-    edgecolor : string
-        Color name.  (Default is 'scaled' to scale the edge colors.)
-    facecolor : string
-        Color name.  (Default is 'scaled' to scale the face colors.)
-    a : numpy.ndarray
-        Array to plot.
-    masked_values : iterable of floats, ints
-        Values to mask.
-    kwargs : dictionary
-        Keyword arguments that are passed to PatchCollection.set(``**kwargs``).
-        Some common kwargs would be 'linewidths', 'linestyles', 'alpha', etc.
-
-    Returns
-    -------
-    pc : matplotlib.collections.PatchCollection
-
-    Examples
-    --------
-
-    """
-    warnings.warn(
-        "plot_cvfd is deprecated and will be removed in version 3.3.5. "
-        "Use PlotMapView for plotting",
-        DeprecationWarning,
-    )
-
-    if "vmin" in kwargs:
-        vmin = kwargs.pop("vmin")
-    else:
-        vmin = None
-
-    if "vmax" in kwargs:
-        vmax = kwargs.pop("vmax")
-    else:
-        vmax = None
-
-    if "ncpl" in kwargs:
-        nlay = layer + 1
-        ncpl = kwargs.pop("ncpl")
-        if isinstance(ncpl, int):
-            i = int(ncpl)
-            ncpl = np.ones((nlay), dtype=int) * i
-        elif isinstance(ncpl, list) or isinstance(ncpl, tuple):
-            ncpl = np.array(ncpl)
-        i0 = 0
-        i1 = 0
-        for k in range(nlay):
-            i0 = i1
-            i1 = i0 + ncpl[k]
-        # retain iverts in selected layer
-        iverts = iverts[i0:i1]
-        # retain vertices in selected layer
-        tverts = []
-        for iv in iverts:
-            for iloc in iv:
-                tverts.append((verts[iloc, 0], verts[iloc, 1]))
-        verts = np.array(tverts)
-        # calculate offset for starting vertex in layer based on
-        # global vertex numbers
-        iadj = iverts[0][0]
-        # reset iverts to relative vertices in selected layer
-        tiverts = []
-        for iv in iverts:
-            i = []
-            for t in iv:
-                i.append(t - iadj)
-            tiverts.append(i)
-        iverts = tiverts
-    else:
-        i0 = 0
-        i1 = len(iverts)
-
-    # get current axis
-    if ax is None:
-        ax = plt.gca()
-    cm = plt.get_cmap(cmap)
-
-    pc = cvfd_to_patch_collection(verts, iverts)
-    pc.set(**kwargs)
-
-    # set colors
-    if a is None:
-        nshp = len(pc.get_paths())
-        cccol = cm(1.0 * np.arange(nshp) / nshp)
-        if facecolor == "scaled":
-            pc.set_facecolor(cccol)
-        else:
-            pc.set_facecolor(facecolor)
-        if edgecolor == "scaled":
-            pc.set_edgecolor(cccol)
-        else:
-            pc.set_edgecolor(edgecolor)
-    else:
-        pc.set_cmap(cm)
-        if masked_values is not None:
-            for mval in masked_values:
-                a = np.ma.masked_equal(a, mval)
-
-        # add NaN values to mask
-        a = np.ma.masked_where(np.isnan(a), a)
-
-        if edgecolor == "scaled":
-            pc.set_edgecolor("none")
-        else:
-            pc.set_edgecolor(edgecolor)
-        pc.set_array(a[i0:i1])
-        pc.set_clim(vmin=vmin, vmax=vmax)
-    # add the patch collection to the axis
-    ax.add_collection(pc)
-    return pc
-
-
 def _set_coord_info(mg, xul, yul, xll, yll, rotation):
     """
 
@@ -2484,68 +2216,6 @@ def _set_coord_info(mg, xul, yul, xll, yll, rotation):
     return mg
 
 
-def _depreciated_dis_handler(modelgrid, dis):
-    """
-    PlotMapView handler for the deprecated dis parameter
-    which adds top and botm information to the modelgrid
-
-    Parameter
-    ---------
-    modelgrid : fp.discretization.Grid object
-
-    dis : fp.modflow.ModflowDis object
-
-    Returns
-    -------
-    modelgrid : fp.discretization.Grid
-
-    """
-    # creates a new modelgrid instance with the dis information
-    from ..discretization import StructuredGrid, UnstructuredGrid, VertexGrid
-
-    warnings.warn(
-        "the dis parameter has been depreciated and will be removed in "
-        "version 3.3.5.",
-        PendingDeprecationWarning,
-    )
-    if modelgrid.grid_type == "vertex":
-        modelgrid = VertexGrid(
-            vertices=modelgrid.vertices,
-            cell2d=modelgrid.cell2d,
-            top=dis.top.array,
-            botm=dis.botm.array,
-            idomain=modelgrid.idomain,
-            xoff=modelgrid.xoffset,
-            yoff=modelgrid.yoffset,
-            angrot=modelgrid.angrot,
-        )
-    if modelgrid.grid_type == "unstructured":
-        modelgrid = UnstructuredGrid(
-            vertices=modelgrid._vertices,
-            iverts=modelgrid._iverts,
-            xcenters=modelgrid._xc,
-            ycenters=modelgrid._yc,
-            top=dis.top.array,
-            botm=dis.botm.array,
-            idomain=modelgrid.idomain,
-            xoff=modelgrid.xoffset,
-            yoff=modelgrid.yoffset,
-            angrot=modelgrid.angrot,
-        )
-    else:
-        modelgrid = StructuredGrid(
-            delc=dis.delc.array,
-            delr=dis.delr.array,
-            top=dis.top.array,
-            botm=dis.botm.array,
-            idomain=modelgrid.idomain,
-            xoff=modelgrid.xoffset,
-            yoff=modelgrid.yoffset,
-            angrot=modelgrid.angrot,
-        )
-    return modelgrid
-
-
 def advanced_package_bc_helper(pkg, modelgrid, kper):
     """
     Helper function for plotting boundary conditions from "advanced" packages