diff --git a/.docs/Notebooks/mf6_simple_model_example.py b/.docs/Notebooks/mf6_simple_model_example.py
index aecc57ed3..6e941a9b2 100644
--- a/.docs/Notebooks/mf6_simple_model_example.py
+++ b/.docs/Notebooks/mf6_simple_model_example.py
@@ -198,8 +198,6 @@
 
 # ### Post-Process Head Results
 #
-# Post-processing MODFLOW 6 results is still a work in progress.  There aren't any Flopy plotting functions built in yet, like they are for other MODFLOW versions.  So we need to plot the results using general Flopy capabilities.  We can also use some of the Flopy ModelMap capabilities for MODFLOW 6, but in order to do so, we need to manually create a SpatialReference object, that is needed for the plotting.  Examples of both approaches are shown below.
-#
 # First, a link to the heads file is created with `HeadFile`. The link can then be accessed with the `get_data` function, by specifying, in this case, the step number and period number for which we want to retrieve data. A three-dimensional array is returned of size `nlay, nrow, ncol`. Matplotlib contouring functions are used to make contours of the layers or a cross-section.
 
 # Read the binary head file and plot the results
diff --git a/flopy/mf6/utils/reference.py b/flopy/mf6/utils/reference.py
index 1f472b505..8db890fbd 100644
--- a/flopy/mf6/utils/reference.py
+++ b/flopy/mf6/utils/reference.py
@@ -1,6 +1,10 @@
 """
 Module spatial referencing for flopy model objects
 
+.. deprecated:: 3.9
+    This module will be removed in FloPy 3.10+. Use
+    the `flopy.discretization` module instead.
+
 """
 
 import numpy as np
@@ -10,6 +14,10 @@ class StructuredSpatialReference:
     """
     a simple class to locate the model grid in x-y space
 
+    .. deprecated:: 3.9
+        This class will be removed in FloPy 3.10+. Use
+        `StructuredGrid` instead.
+
     Parameters
     ----------
 
@@ -544,6 +552,10 @@ class VertexSpatialReference:
     """
     a simple class to locate the model grid in x-y space
 
+    .. deprecated:: 3.9
+        This class will be removed in FloPy 3.10+. Use
+        `VertexGrid` instead.
+
     Parameters
     ----------
     xvdict: dictionary
@@ -852,33 +864,36 @@ class SpatialReference:
     """
     A dynamic inheritance class that locates a gridded model in space
 
-        Parameters
-        ----------
-        delr : numpy ndarray
-            the model discretization delr vector
-        delc : numpy ndarray
-            the model discretization delc vector
-        lenuni : int
-            the length units flag from the discretization package
-        xul : float
-            the x coordinate of the upper left corner of the grid
-        yul : float
-           the y coordinate of the upper left corner of the grid
-        rotation : float
-            the counter-clockwise rotation (in degrees) of the grid
-        proj4_str: str
-            a PROJ4 string that identifies the grid in space. warning:
-            case sensitive!
-        xadj : float
-            vertex grid: x vertex adjustment factor
-        yadj : float
-            vertex grid: y vertex adjustment factor
-        xvdict: dict
-            dictionary of x-vertices by cellnum ex. {0: (0,1,1,0)}
-        yvdict: dict
-            dictionary of y-vertices by cellnum ex. {0: (1,1,0,0)}
-        distype: str
-            model grid discretization type
+    .. deprecated:: 3.9
+        This class will be removed in FloPy 3.10+.
+
+    Parameters
+    ----------
+    delr : numpy ndarray
+        the model discretization delr vector
+    delc : numpy ndarray
+        the model discretization delc vector
+    lenuni : int
+        the length units flag from the discretization package
+    xul : float
+        the x coordinate of the upper left corner of the grid
+    yul : float
+        the y coordinate of the upper left corner of the grid
+    rotation : float
+        the counter-clockwise rotation (in degrees) of the grid
+    proj4_str: str
+        a PROJ4 string that identifies the grid in space. warning:
+        case sensitive!
+    xadj : float
+        vertex grid: x vertex adjustment factor
+    yadj : float
+        vertex grid: y vertex adjustment factor
+    xvdict: dict
+        dictionary of x-vertices by cellnum ex. {0: (0,1,1,0)}
+    yvdict: dict
+        dictionary of y-vertices by cellnum ex. {0: (1,1,0,0)}
+    distype: str
+        model grid discretization type
     """
 
     def __new__(