diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index 8b6b53a5b..3d6874a27 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -16,7 +16,13 @@
 from autotest.test_dis_cases import case_dis, case_disv
 from autotest.test_grid_cases import GridCases
 from flopy.discretization import StructuredGrid, UnstructuredGrid, VertexGrid
-from flopy.mf6 import MFSimulation
+from flopy.mf6 import (
+    MFSimulation,
+    ModflowGwf,
+    ModflowGwfdis,
+    ModflowGwfdisu,
+    ModflowGwfdisv,
+)
 from flopy.modflow import Modflow, ModflowDis
 from flopy.utils import import_optional_dependency
 from flopy.utils.crs import get_authority_crs
@@ -1903,3 +1909,90 @@ def test_unstructured_grid_get_node():
 
     with pytest.raises(IndexError, match=r"Node .* out of range"):
         ug.get_node(200)
+
+
+@pytest.mark.mf6
+def test_structured_mf6_gridprops(example_data_path):
+    sim = MFSimulation.load(sim_ws=example_data_path / "mf6-freyberg")
+    gwf = sim.get_model()
+    dis = gwf.dis
+    modelgrid = gwf.modelgrid
+
+    new_sim = MFSimulation()
+    new_gwf = ModflowGwf(new_sim)
+    new_dis = ModflowGwfdis(new_gwf, **modelgrid.dis_properties())
+    attrs = ("delc", "delr", "top", "botm", "idomain", "xorigin", "yorigin", "angrot")
+    for attr in attrs:
+        v0 = getattr(dis, attr).array
+        v1 = getattr(new_dis, attr).array
+        if attr in ("xorigin", "yorigin", "angrot") and v0 is None:
+            v0 = 0
+        np.testing.assert_allclose(
+            v0, v1, err_msg=f"{attr} not consistent with valid array data"
+        )
+
+
+@pytest.mark.mf6
+def test_vertex_mf6_gridprops(example_data_path):
+    sim = MFSimulation.load(sim_ws=example_data_path / "mf6" / "test003_gwftri_disv")
+    gwf = sim.get_model()
+    disv = gwf.disv
+    modelgrid = gwf.modelgrid
+
+    new_sim = MFSimulation()
+    new_gwf = ModflowGwf(new_sim)
+    new_disv = ModflowGwfdisv(new_gwf, **modelgrid.disv_properties())
+
+    attrs = (
+        "vertices",
+        "top",
+        "botm",
+        "idomain",
+        "xorigin",
+        "yorigin",
+        "angrot",
+        "cell2d",
+    )
+    for attr in attrs:
+        v0 = getattr(disv, attr).array
+        v1 = getattr(new_disv, attr).array
+        if attr in ("xorigin", "yorigin", "angrot") and v0 is None:
+            v0 = 0
+
+        if attr in ("cell2d", "vertices"):
+            for col in v0.dtype.names:
+                np.testing.assert_allclose(
+                    v0[col],
+                    v1[col],
+                    err_msg=f"{attr} not consistent with valid array data",
+                )
+        else:
+            np.testing.assert_allclose(
+                v0, v1, err_msg=f"{attr} not consistent with valid array data"
+            )
+
+
+@pytest.mark.mf6
+def test_unstructured_mf6_gridprops(example_data_path):
+    sim = MFSimulation.load(sim_ws=example_data_path / "mf6" / "test006_gwf3")
+    gwf = sim.get_model()
+    disu = gwf.disu
+    modelgrid = gwf.modelgrid
+
+    new_sim = MFSimulation()
+    new_gwf = ModflowGwf(new_sim)
+    dis_props = modelgrid.disu_properties()
+    dis_props["area"] = disu.area.array
+    dis_props["cl12"] = disu.cl12.array
+    new_disu = ModflowGwfdisu(new_gwf, **dis_props)
+
+    attrs = ("top", "bot", "iac", "ja", "nodes", "ihc", "xorigin", "yorigin", "angrot")
+    for attr in attrs:
+        v0 = getattr(disu, attr).array
+        v1 = getattr(new_disu, attr).array
+        if attr in ("xorigin", "yorigin", "angrot") and v0 is None:
+            v0 = 0
+
+        np.testing.assert_allclose(
+            v0, v1, err_msg=f"{attr} not consistent with valid array data"
+        )
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 4c0b74e50..fd74ade4e 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -535,6 +535,25 @@ def xyzextent(self):
             np.max(self.xyzvertices[2]),
         )
 
+    @property
+    def area(self):
+        """
+        Returns a numpy array of cell areas calculated using the shoelace algorithm
+        """
+        # irregular_shape_patch
+        from ..plot.plotutil import UnstructuredPlotUtilities
+
+        # when looping through to create determinants, need to start at -1
+        xverts, yverts = self.cross_section_vertices
+        xverts, yverts = UnstructuredPlotUtilities.irregular_shape_patch(xverts, yverts)
+        area_x2 = np.zeros((1, len(xverts)))
+        for i in range(xverts.shape[-1]):
+            # calculate the determinant of each line in polygon
+            area_x2 += xverts[:, i - 1] * yverts[:, i] - yverts[:, i - 1] * xverts[:, i]
+
+        area = np.abs(area_x2 / 2.0)
+        return np.ravel(area)
+
     @property
     def grid_lines(self):
         raise NotImplementedError("must define grid_lines in child class")
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index 18bab69f8..f4efd000b 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -759,6 +759,44 @@ def map_polygons(self):
 
         return self._polygons
 
+    def dis_properties(self, mf2005=False):
+        """
+        Method to get DIS package properties
+
+        Parameters
+        ----------
+        mf2005 : bool
+            flag to get legacy mf2005/mfnwt discretization package properties
+            from the modelgrid object
+
+        Returns
+        -------
+        dict : dictionary of discretization properties that can be used to build a
+            DIS package
+        """
+        dis_props = {
+            "delc": self.__delc,
+            "delr": self.__delr,
+            "top": self.top,
+            "botm": self.botm,
+            "nlay": self.nlay,
+            "nrow": self.nrow,
+            "ncol": self.ncol,
+        }
+
+        if mf2005:
+            if self.is_valid:
+                dis_props["xul"] = self.xvertices[0, 0]
+                dis_props["yul"] = self.yvertices[0, 0]
+                dis_props["rotation"] = self.angrot
+        else:
+            dis_props["xorigin"] = self.xoffset
+            dis_props["yorigin"] = self.yoffset
+            dis_props["angrot"] = self.angrot
+            dis_props["idomain"] = self.idomain
+
+        return dis_props
+
     def to_geodataframe(self):
         """
         Returns a geopandas GeoDataFrame of the model grid
diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py
index 2d81e7b4a..48a8cc70f 100644
--- a/flopy/discretization/unstructuredgrid.py
+++ b/flopy/discretization/unstructuredgrid.py
@@ -591,6 +591,52 @@ def map_polygons(self):
 
         return copy.copy(self._polygons)
 
+    def disu_properties(self, mfusg=False):
+        """
+        Method that returns disu properties from the grid for constructing
+        DISU packages for MFUSG and MF6. Note: not all required information for
+        DISU construction is stored in the UnstructuredGrid class, CL12 and HWVA
+        is not available from the Grid. Please double-check the properties returned
+        on a case by case basis and fill in where necessary for individual applications
+
+        Returns
+        -------
+        dict : dictionary of unstructured discretization properties
+        """
+
+        dis_props = {
+            "top": self._top,
+            "bot": self._botm,
+            "iac": self._iac,
+            "ja": self._ja,
+        }
+        if mfusg:
+            dis_props["nodelay"] = self.ncpl
+            dis_props["ivc"] = np.where(self._ihc < 1, 1, 0)
+
+        else:
+            dis_props["nodes"] = self.nnodes
+            dis_props["ihc"] = self._ihc
+            dis_props["idomain"] = self.idomain
+            dis_props["xorigin"] = self.xoffset
+            dis_props["yorigin"] = self.yoffset
+            dis_props["angrot"] = self.angrot
+
+            if self.is_valid:
+                dis_props["vertices"] = self._vertices
+                cell2d = []
+                for ix, iv in enumerate(self._iverts):
+                    c2d = tuple(
+                        [ix + 1, self._xc[ix], self._yc[ix], len(iv)] + list(iv)
+                    )
+                    cell2d.append(c2d)
+                dis_props["cell2d"] = cell2d
+
+        if self.is_valid:
+            dis_props["area"] = self.area
+
+        return dis_props
+
     def to_geodataframe(self):
         """
         Returns a geopandas GeoDataFrame of the model grid
diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py
index ddb8da4b5..757917d21 100644
--- a/flopy/discretization/vertexgrid.py
+++ b/flopy/discretization/vertexgrid.py
@@ -300,6 +300,33 @@ def map_polygons(self):
 
         return copy.copy(self._polygons)
 
+    def disv_properties(self):
+        """
+        Method to get DISV package properties
+
+        Returns
+        -------
+        dict : dictionary of properties that can be used to build a DISV package
+        """
+        dis_props = {
+            "nlay": self.nlay,
+            "ncpl": self.ncpl,
+            "vertices": self._vertices,
+            "top": self.top,
+            "botm": self.botm,
+            "idomain": self.idomain,
+            "xorigin": self.xoffset,
+            "yorigin": self.yoffset,
+            "angrot": self.angrot,
+        }
+
+        if self._cell2d is None and self._cell1d is not None:
+            dis_props["cell1d"] = self._cell1d
+        else:
+            dis_props["cell2d"] = self._cell2d
+
+        return dis_props
+
     def to_geodataframe(self):
         """
         Returns a geopandas GeoDataFrame of the model grid