diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py
index 8f4d265a1..bc9ec5583 100644
--- a/autotest/test_binaryfile.py
+++ b/autotest/test_binaryfile.py
@@ -275,7 +275,7 @@ def test_load_binary_head_file(example_data_path):
 
 def test_plot_binary_head_file(example_data_path):
     hf = HeadFile(example_data_path / "freyberg" / "freyberg.githds")
-    hf.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
+    hf.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
 
     assert isinstance(hf.plot(), Axes)
     plt.close()
diff --git a/autotest/test_formattedfile.py b/autotest/test_formattedfile.py
index 0d31ed308..26ac20f45 100644
--- a/autotest/test_formattedfile.py
+++ b/autotest/test_formattedfile.py
@@ -69,7 +69,7 @@ def test_headfile_build_index(example_data_path):
 def test_formattedfile_reference(example_data_path):
     h = FormattedHeadFile(example_data_path / "mf2005_test" / "test1tr.githds")
     assert isinstance(h, FormattedHeadFile)
-    h.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
+    h.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
 
     assert isinstance(h.plot(masked_values=[6999.000]), Axes)
     plt.close()
diff --git a/autotest/test_modpathfile.py b/autotest/test_modpathfile.py
index 2bedeeb3f..9f667fac8 100644
--- a/autotest/test_modpathfile.py
+++ b/autotest/test_modpathfile.py
@@ -306,11 +306,10 @@ def test_get_destination_endpoint_data(
     )
 
 
-@pytest.mark.parametrize("longfieldname", [True, False])
 @requires_exe("mf6", "mp7")
-@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"})
-def test_write_shapefile(function_tmpdir, mp7_small, longfieldname):
-    from shapefile import Reader
+@requires_pkg("geopandas", "shapely")
+def test_write_shapefile(function_tmpdir, mp7_small):
+    import geopandas as gpd
 
     # setup and run model, then copy outputs to function_tmpdir
     sim, forward_model_name, _, _, _ = mp7_small
@@ -330,13 +329,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname):
     # define shapefile path
     shp_file = ws / "pathlines.shp"
 
-    # add a column to the pathline recarray
-    fieldname = "newfield" + ("longname" if longfieldname else "")
-    fieldval = "x"
-    pathlines = [
-        rfn.append_fields(pl, fieldname, list(repeat(fieldval, len(pl))), dtypes="|S1")
-        for pl in pathlines
-    ]
+    pline_names = [name[:10] for name in pathlines[0].dtype.names]
 
     # write the pathline recarray to shapefile
     pathline_file.write_shapefile(
@@ -350,8 +343,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname):
     assert shp_file.is_file()
 
     # load shapefile
-    with Reader(shp_file) as reader:
-        fieldnames = [f[0] for f in reader.fields[1:]]
-        fieldname = "newfiname_" if longfieldname else fieldname
-        assert fieldname in fieldnames
-        assert all(r[fieldname] == fieldval for r in reader.iterRecords())
+    gdf = gpd.read_file(shp_file)
+    fieldnames = list(gdf)
+    for fname in pline_names:
+        assert fname in fieldnames
diff --git a/flopy/mbase.py b/flopy/mbase.py
index 40cdbf700..ce7589f03 100644
--- a/flopy/mbase.py
+++ b/flopy/mbase.py
@@ -759,6 +759,39 @@ def _output_msg(self, i, add=True):
             f"{txt2} the output list."
         )
 
+    def to_geo_dataframe(self, gdf=None, kper=0):
+        """
+        Method to build a Geodataframe from model inputs. Note: transient data
+        will only be exported for a single stress period.
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        if gdf is None:
+            modelgrid = self.modelgrid
+            if modelgrid is not None:
+                gdf = modelgrid.geo_dataframe
+            else:
+                raise AttributeError(
+                    "model does not have a grid instance, please supply a geodataframe"
+                )
+
+        for package in self.packagelist:
+            if package.package_type in ("hfb6",):
+                continue
+            if callable(getattr(package, "to_geo_dataframe", None)):
+                gdf = package.to_geo_dataframe(gdf, kper=kper, sparse=False)
+
+        return gdf
+
     def add_output_file(
         self,
         unit,
diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py
index e00c877b1..ef3e78508 100644
--- a/flopy/mf6/data/mfdataarray.py
+++ b/flopy/mf6/data/mfdataarray.py
@@ -323,6 +323,64 @@ def data(self):
         """Returns array data.  Calls get_data with default parameters."""
         return self._get_data()
 
+    def to_geo_dataframe(self, gdf=None, name=None, forgive=False, **kwargs):
+        """
+        Method to add an input array to a geopandas GeoDataFrame
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame object
+        name : str
+            optional attribute name, default uses util2d name
+        forgive : bool
+            optional flag to continue if data shape not compatible with GeoDataFrame
+
+        Returns
+        -------
+            geopandas GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if gdf is None:
+                if modelgrid is None:
+                    return gdf
+                gdf = modelgrid.geo_dataframe
+
+            if modelgrid is not None:
+                if modelgrid.grid_type != "unstructured":
+                    ncpl = modelgrid.ncpl
+                else:
+                    ncpl = modelgrid.nnodes
+            else:
+                ncpl = len(gdf)
+
+            if name is None:
+                name = self.name
+
+            data = self.array
+            if data is None:
+                return gdf
+
+            if data.size == ncpl:
+                gdf[name] = data.ravel()
+
+            elif data.size % ncpl == 0:
+                data = data.reshape((-1, ncpl))
+                for ix, arr in enumerate(data):
+                    aname = f"{name}_{ix}"
+                    gdf[aname] = arr
+            elif forgive:
+                return gdf
+            else:
+                raise ValueError(
+                    f"Data size {data.size} not compatible with dataframe length {ncpl}"
+                )
+
+            return gdf
+
     def new_simulation(self, sim_data):
         """Initialize MFArray object for a new simulation
 
@@ -1890,6 +1948,50 @@ def _build_period_data(
                     output[sp] = data
         return output
 
+    def to_geo_dataframe(self, gdf=None, kper=0, forgive=False, **kwargs):
+        """
+
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if gdf is None:
+                if modelgrid is None:
+                    return gdf
+                gdf = modelgrid.geo_dataframe
+
+            if modelgrid is not None:
+                if modelgrid.grid_type != "unstructured":
+                    ncpl = modelgrid.ncpl
+                else:
+                    ncpl = modelgrid.nnodes
+            else:
+                ncpl = len(gdf)
+
+            if self.array is None:
+                return gdf
+
+            name = f"{self.path[1]}_{self.name}"
+
+            data = self.get_data(key=kper, apply_mult=True)
+            if data.size == ncpl:
+                gdf[name] = data.ravel()
+
+            elif data.size % ncpl == 0:
+                data = data.reshape((-1, ncpl))
+                for ix, arr in enumerate(data):
+                    aname = f"{name}_{ix}"
+                    gdf[aname] = arr
+            elif forgive:
+                return gdf
+            else:
+                raise ValueError(
+                    f"Data size {data.size} not compatible with dataframe length {ncpl}"
+                )
+
+            return gdf
+
     def set_record(self, data_record):
         """Sets data and metadata at layer `layer` and time `key` to
         `data_record`.  For unlayered data do not pass in `layer`.
diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py
index 529177c41..3d05e79a2 100644
--- a/flopy/mf6/data/mfdatalist.py
+++ b/flopy/mf6/data/mfdatalist.py
@@ -145,6 +145,55 @@ def to_array(self, kper=0, mask=False):
         model_grid = self.data_dimensions.get_model_grid()
         return list_to_array(sarr, model_grid, kper, mask)
 
+    def to_geo_dataframe(self, gdf=None, sparse=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.geo_dataframe
+
+            data = self.to_array(mask=True)
+            if data is None:
+                return gdf
+
+            col_names = []
+            for name, array3d in data.items():
+                aname = f"{self.path[1].lower()}_{name}"
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        gdf[f"{aname}_{lay}"] = arr.ravel()
+                        col_names.append(f"{aname}_{lay}")
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def new_simulation(self, sim_data):
         """Initialize MFList object for a new simulation.
 
@@ -1596,6 +1645,57 @@ def to_array(self, kper=0, mask=False):
         """Returns list data as an array."""
         return super().to_array(kper, mask)
 
+    def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.geo_dataframe
+
+            data = self.to_array(kper=kper, mask=True)
+            if data is None:
+                return gdf
+
+            col_names = []
+            for name, array3d in data.items():
+                aname = f"{self.path[1].lower()}_{name}"
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        gdf[f"{aname}_{lay}"] = arr.ravel()
+                        col_names.append(f"{aname}_{lay}")
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def remove_transient_key(self, transient_key):
         """Remove transient stress period key.  Method is used
         internally by FloPy and is not intended to the end user.
diff --git a/flopy/mf6/data/mfdataplist.py b/flopy/mf6/data/mfdataplist.py
index 331177ef0..60c9162cf 100644
--- a/flopy/mf6/data/mfdataplist.py
+++ b/flopy/mf6/data/mfdataplist.py
@@ -839,6 +839,55 @@ def to_array(self, kper=0, mask=False):
         model_grid = self.data_dimensions.get_model_grid()
         return list_to_array(sarr, model_grid, kper, mask)
 
+    def to_geo_dataframe(self, gdf=None, sparse=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.geo_dataframe
+
+            data = self.to_array(mask=True)
+            if data is None:
+                return gdf
+
+            col_names = []
+            for name, array3d in data.items():
+                aname = f"{self.path[1].lower()}_{name}"
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        gdf[f"{aname}_{lay}"] = arr.ravel()
+                        col_names.append(f"{aname}_{lay}")
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def set_record(self, record, autofill=False, check_data=True):
         """Sets the contents of the data and metadata to "data_record".
         Data_record is a dictionary with has the following format:
@@ -1935,9 +1984,7 @@ def plot(
         )
 
 
-class MFPandasTransientList(
-    MFPandasList, mfdata.MFTransient, DataListInterface
-):
+class MFPandasTransientList(MFPandasList, mfdata.MFTransient, DataListInterface):
     """
     Provides an interface for the user to access and update MODFLOW transient
     pandas list data.
@@ -2019,6 +2066,55 @@ def to_array(self, kper=0, mask=False):
         """Returns list data as an array."""
         return super().to_array(kper, mask)
 
+    def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.geo_dataframe
+
+            data = self.to_array(kper=kper, mask=True)
+
+            col_names = []
+            for name, array3d in data.items():
+                aname = f"{self.path[1].lower()}_{name}"
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        gdf[f"{aname}_{lay}"] = arr.ravel()
+                        col_names.append(f"{aname}_{lay}")
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def remove_transient_key(self, transient_key):
         """Remove transient stress period key.  Method is used
         internally by FloPy and is not intended to the end user.
diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py
index 3260ca12c..bcbd600ea 100644
--- a/flopy/mf6/data/mfdatautil.py
+++ b/flopy/mf6/data/mfdatautil.py
@@ -202,6 +202,8 @@ def list_to_array(sarr, model_grid, kper=0, mask=False):
         cnt = np.zeros(shape, dtype=np.float64)
         for sp_rec in sarr:
             if sp_rec is not None:
+                if "cellid" not in sp_rec.dtype.names:
+                    return None
                 for rec in sp_rec:
                     arr[rec["cellid"]] += rec[name]
                     cnt[rec["cellid"]] += 1.0
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
index 0abc25cbd..0cb7e321a 100644
--- a/flopy/mf6/mfmodel.py
+++ b/flopy/mf6/mfmodel.py
@@ -789,6 +789,40 @@ def output(self):
         except AttributeError:
             return MF6Output(self, budgetkey=budgetkey)
 
+    def to_geo_dataframe(self, gdf=None, kper=0):
+        """
+        Method to build a Geodataframe from model inputs. Note: transient data
+        will only be exported for a single stress period.
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        if gdf is None:
+            modelgrid = self.modelgrid
+            if modelgrid is not None:
+                gdf = modelgrid.geo_dataframe
+            else:
+                raise AttributeError(
+                    "model does not have a grid instance, "
+                    "please supply a geodataframe"
+                )
+
+        for package in self.packagelist:
+            if package.package_type in ("hfb",):
+                continue
+            if callable(getattr(package, "to_geo_dataframe", None)):
+                gdf = package.to_geo_dataframe(gdf, kper=kper, sparse=False)
+
+        return gdf
+
     def export(self, f, **kwargs):
         """Method to export a model to a shapefile or netcdf file
 
diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py
index 4ee9d6175..d2b0af8cd 100644
--- a/flopy/mf6/mfpackage.py
+++ b/flopy/mf6/mfpackage.py
@@ -12,6 +12,7 @@
 from ..pakbase import PackageInterface
 from ..utils import datautil
 from ..utils.check import mf6check
+from ..utils.utl_import import import_optional_dependency
 from ..version import __version__
 from .coordinates import modeldimensions
 from .data import (
@@ -2086,6 +2087,76 @@ def package_filename_dict(self):
         )
         return self._package_container.package_filename_dict
 
+    def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs):
+        """
+        Method to create a GeoDataFrame from a modflow package
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        if gdf is None:
+            if isinstance(self.parent, ModelInterface):
+                modelgrid = self.parent.modelgrid
+                if modelgrid is not None:
+                    if self.package_type == "hfb":
+                        gpd = import_optional_dependency("geopandas")
+                        from ..plot.plotutil import hfb_data_to_linework
+
+                        recarray = self.stress_period_data.data[kper]
+                        lines = hfb_data_to_linework(recarray, modelgrid)
+                        geo_interface = {"type": "FeatureCollection"}
+                        features = [
+                            {
+                                "id": f"{ix}",
+                                "geometry": {"coordinates": line, "type": "LineString"},
+                                "properties": {}
+                            }
+                            for ix, line in enumerate(lines)
+                        ]
+                        geo_interface["features"] = features
+                        gdf = gpd.GeoDataFrame.from_features(geo_interface)
+
+                        for name in recarray.dtype.names:
+                            gdf[name] = recarray[name]
+
+                        return gdf
+
+                    else:
+                        gdf = modelgrid.geo_dataframe
+                else:
+                    raise AttributeError(
+                        "model does not have a grid instance, "
+                        "please supply a geodataframe"
+                    )
+            else:
+                raise AssertionError(
+                    "Package does not have a model instance, "
+                    "please supply a geodataframe"
+                )
+
+        for attr, value in self.__dict__.items():
+            if callable(getattr(value, "to_geo_dataframe", None)):
+                if isinstance(value, (ModelInterface, PackageInterface)):
+                    continue
+                # do not pass sparse in here, "sparsify" after all data has been
+                #  added to geodataframe
+                gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper, sparse=False)
+
+        if sparse:
+            col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")]
+            gdf = gdf.dropna(subset=col_names, how="all")
+            gdf = gdf.dropna(axis="columns", how="all")
+
+        return gdf
+
     def get_package(self, name=None, type_only=False, name_only=False):
         """
         Finds a package by package name, package key, package type, or partial
diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py
index e018f54d4..81df8a55c 100644
--- a/flopy/modflow/mfhfb.py
+++ b/flopy/modflow/mfhfb.py
@@ -14,6 +14,7 @@
 from ..pakbase import Package
 from ..utils.flopy_io import line_parse
 from ..utils.recarray_utils import create_empty_recarray
+from ..utils.utl_import import import_optional_dependency
 from .mfparbc import ModflowParBc as mfparbc
 
 
@@ -179,6 +180,47 @@ def _ncells(self):
         """
         return self.nhfbnp
 
+    def _get_hfb_lines(self):
+        """
+        Method to get hfb lines. Lines are ordered by recarray index
+
+        Returns
+        -------
+            list : list of hfb lines
+        """
+        from ..plot.plotutil import hfb_data_to_linework
+
+        return hfb_data_to_linework(self.hfb_data, self.parent.modelgrid)
+
+    def to_geo_dataframe(self, **kwargs):
+        """
+        Method to create a LineString based geodataframe of horizontal flow barriers
+
+        Returns
+        -------
+            GeoDataFrame
+
+        """
+        gpd = import_optional_dependency("geopandas")
+
+        lines = self._get_hfb_lines()
+        geo_interface = {"type": "FeatureCollection"}
+        features = [
+            {
+                "id": f"{ix}",
+                "geometry": {"coordinates": line, "type": "LineString"},
+                "properties": {},
+            }
+            for ix, line in enumerate(lines)
+        ]
+        geo_interface["features"] = features
+        gdf = gpd.GeoDataFrame.from_features(geo_interface)
+        hfbs = self.hfb_data
+        for name in hfbs.dtype.names:
+            gdf[name] = hfbs[name]
+
+        return gdf
+
     def write_file(self):
         """
         Write the package file.
@@ -214,6 +256,9 @@ def get_empty(ncells=0, aux_names=None, structured=True):
         been extended to include aux variables and associated
         aux names.
 
+        Returns
+        -------
+            np.recarray
         """
         dtype = ModflowHfb.get_default_dtype(structured=structured)
         if aux_names is not None:
diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py
index 1abd2ce50..27e791b5e 100644
--- a/flopy/modflow/mfsfr2.py
+++ b/flopy/modflow/mfsfr2.py
@@ -623,6 +623,43 @@ def paths(self):
     def df(self):
         return pd.DataFrame(self.reach_data)
 
+    def to_geo_dataframe(self, gdf=None, sparse=True, **kwargs):
+        """
+        Method to export SFR reach data to a GeoDataFrame
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is True
+        """
+        modelgrid = self.parent.modelgrid
+        if modelgrid is None:
+            return gdf
+
+        if gdf is None:
+            gdf = modelgrid.geo_dataframe
+
+        df = self.df
+        if "k" in list(df):
+            df["node"] = df["node"] - (modelgrid.ncpl * df["k"])
+            df = df.drop(columns=["k", "i", "j"])
+
+        df["node"] += 1
+        gdf = gdf.merge(df, how="left", on="node")
+
+        if sparse:
+            col_names = [col for col in list(df) if col != "node"]
+            gdf = gdf.dropna(subset=col_names, how="all")
+            gdf = gdf.dropna(axis="columns", how="all")
+        else:
+            gdf = gdf.drop_duplicates(subset=["node"])
+            gdf = gdf.reset_index(drop=True)
+
+        return gdf
+
     def _make_graph(self):
         # get all segments and their outseg
         graph = {}
diff --git a/flopy/pakbase.py b/flopy/pakbase.py
index 8873254a6..8718f188d 100644
--- a/flopy/pakbase.py
+++ b/flopy/pakbase.py
@@ -673,6 +673,54 @@ def export(self, f, **kwargs):
 
         return export.utils.package_export(f, self, **kwargs)
 
+    def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs):
+        """
+        Method to create a GeoDataFrame from a modflow package
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        from .mbase import BaseModel
+
+        if gdf is None:
+            if isinstance(self.parent, BaseModel):
+                modelgrid = self.parent.modelgrid
+                if modelgrid is not None:
+                    gdf = modelgrid.geo_dataframe
+                else:
+                    raise AttributeError(
+                        "model does not have a grid instance, "
+                        "please supply a geodataframe"
+                    )
+            else:
+                raise AssertionError(
+                    "Package does not have a model instance, "
+                    "please supply a geodataframe"
+                )
+
+        for attr, value in self.__dict__.items():
+            if callable(getattr(value, "to_geo_dataframe", None)):
+                if isinstance(value, (BaseModel, PackageInterface)):
+                    continue
+                # do not pass sparse in here, make sparse after all data has been
+                #  added to geodataframe
+                gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper, sparse=False)
+
+        if sparse:
+            col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")]
+            gdf = gdf.dropna(subset=col_names, how="all")
+            gdf = gdf.dropna(axis="columns", how="all")
+
+        return gdf
+
     def _generate_heading(self):
         """Generate heading."""
         from . import __version__
diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py
index 8c1bfec41..17393f6f7 100644
--- a/flopy/plot/plotutil.py
+++ b/flopy/plot/plotutil.py
@@ -2954,3 +2954,71 @@ def to_prt_pathlines(
         return df
     else:
         return ret
+
+
+def hfb_data_to_linework(recarray, modelgrid):
+    """
+    Method to get lines representing horizontal flow barriers
+
+    Parameters
+    ----------
+    recarray : np.recarray
+        recarray of hfb input data
+    modelgrid : Grid
+        modelgrid instance
+
+    Returns
+    -------
+        list : list of line segments
+    """
+    iverts = modelgrid.iverts
+    verts = modelgrid.verts
+    nodes = []
+    if modelgrid.grid_type == "structured":
+        if "k" in recarray.dtype.names:
+            for rec in recarray:
+                node1 = modelgrid.get_node([(0, rec["irow1"], rec["icol1"])])[0]
+                node2 = modelgrid.get_node([(0, rec["irow2"], rec["icol2"])])[0]
+                nodes.append((node1, node2))
+        else:
+            for rec in recarray:
+                node1 = modelgrid.get_node([(0,) + rec["cellid1"][1:]])[0]
+                node2 = modelgrid.get_node([(0,) + rec["cellid2"][1:]])[0]
+                nodes.append((node1, node2))
+
+    elif modelgrid.grid_type == "vertex":
+        for rec in recarray:
+            nodes.append((rec["cellid1"][-1], rec["cellid2"][-1]))
+
+    else:
+        if "node1" in recarray.dtype.names:
+            nodes = list(zip(recarray["node1"], recarray["node2"]))
+        else:
+            for rec in recarray:
+                nodes.append((rec["cellid1"][0], rec["cellid2"][0]))
+
+    shared_edges = []
+    for node0, node1 in nodes:
+        iv0 = iverts[node0]
+        iv1 = iverts[node1]
+        edges = []
+        for ix in range(len(iv0)):
+            edges.append(tuple(sorted((iv0[ix - 1], iv0[ix]))))
+
+        for ix in range(len(iv1)):
+            edge = tuple(sorted((iv1[ix - 1], iv1[ix])))
+            if edge in edges:
+                shared_edges.append(edge)
+                break
+
+        if not shared_edges:
+            raise AssertionError(
+                f"No shared cell edges found. Cannot represent HFB "
+                f"for nodes {node0} and {node1}"
+            )
+
+    lines = []
+    for edge in shared_edges:
+        lines.append((tuple(verts[edge[0]]), tuple(verts[edge[1]])))
+
+    return lines
diff --git a/flopy/utils/binaryfile/__init__.py b/flopy/utils/binaryfile/__init__.py
index 5c9842b74..50ea1a278 100644
--- a/flopy/utils/binaryfile/__init__.py
+++ b/flopy/utils/binaryfile/__init__.py
@@ -1604,6 +1604,56 @@ def get_data(
             if t
         ]
 
+    def to_geo_dataframe(
+        self,
+        gdf=None,
+        modelgrid=None,
+        idx=None,
+        kstpkper=None,
+        totim=None,
+        text=None,
+    ):
+        if (idx, kstpkper, totim) == (None, None, None):
+            raise AssertionError(
+                "to_geo_dataframe() missing 1 required argument: "
+                "please provide 'idx', 'kstpkper', or 'totim'"
+            )
+
+        if gdf is None:
+            if modelgrid is None:
+                if self.modelgrid is None:
+                    raise AssertionError(
+                        "A geodataframe or modelgrid instance must be supplied"
+                    )
+                modelgrid = self.modelgrid
+
+            gdf = modelgrid.geo_dataframe
+
+        col_names = []
+        if text is None:
+            textlist = [i.decode() for i in self.textlist]
+            for text in textlist:
+                data = self.get_data(
+                    idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True
+                )
+
+                text = text.strip().lower().replace(" ", "_")
+                for ix, arr in enumerate(data[0]):
+                    name = f"{text}_{ix}"
+                    gdf[name] = arr.ravel()
+                    col_names.append(name)
+        else:
+            data = self.get_data(
+                idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True
+            )
+            text = text.strip().lower().replace(" ", "_")
+            for ix, arr in enumerate(data[0]):
+                name = f"{text}_{ix}"
+                gdf[name] = arr.ravel()
+                col_names.append(name)
+
+        return gdf
+
     def get_ts(self, idx, text=None, times=None):
         """
         Get a time series from the binary budget file.
diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py
index 183fb874e..78da71c2f 100644
--- a/flopy/utils/datafile.py
+++ b/flopy/utils/datafile.py
@@ -197,18 +197,18 @@ def __init__(self, filename: Union[str, os.PathLike], precision, verbose, **kwar
 
         self.model = None
         self.dis = None
-        self.mg = None
+        self.modelgrid = None
         if "model" in kwargs.keys():
             self.model = kwargs.pop("model")
-            self.mg = self.model.modelgrid
+            self.modelgrid = self.model.modelgrid
             self.dis = self.model.dis
         if "dis" in kwargs.keys():
             self.dis = kwargs.pop("dis")
-            self.mg = self.dis.parent.modelgrid
+            self.modelgrid = self.dis.parent.modelgrid
         if "tdis" in kwargs.keys():
             self.tdis = kwargs.pop("tdis")
         if "modelgrid" in kwargs.keys():
-            self.mg = kwargs.pop("modelgrid")
+            self.modelgrid = kwargs.pop("modelgrid")
         if len(kwargs.keys()) > 0:
             args = ",".join(kwargs.keys())
             raise Exception(f"LayerFile error: unrecognized kwargs: {args}")
@@ -218,8 +218,8 @@ def __init__(self, filename: Union[str, os.PathLike], precision, verbose, **kwar
 
         # now that we read the data and know nrow and ncol,
         # we can make a generic mg if needed
-        if self.mg is None:
-            self.mg = StructuredGrid(
+        if self.modelgrid is None:
+            self.modelgrid = StructuredGrid(
                 delc=np.ones((self.nrow,)),
                 delr=np.ones(self.ncol),
                 nlay=self.nlay,
@@ -240,6 +240,54 @@ def __enter__(self):
     def __exit__(self, *exc):
         self.close()
 
+    def to_geo_dataframe(
+        self, gdf=None, modelgrid=None, kstpkper=None, totim=None, attrib_name=None
+    ):
+        """
+        Generate a GeoDataFrame with data from a LayerFile instance
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional, existing geodataframe with NCPL geometries
+        modelgrid : Grid
+            optional modelgrid instance to generate a GeoDataFrame from
+        kstpkper : tuple of ints
+            A tuple containing the time step and stress period (kstp, kper).
+            These are zero-based kstp and kper values.
+        totim : float
+            The simulation time.
+        attrib_name : str
+            optional base name of attribute columns. (default is text attribute)
+
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if gdf is None:
+            if modelgrid is None:
+                if self.modelgrid is None:
+                    raise AssertionError(
+                        "A geodataframe or modelgrid instance must be supplied"
+                    )
+                modelgrid = self.modelgrid
+
+            gdf = modelgrid.geo_dataframe
+
+        array = np.atleast_3d(
+            self.get_data(kstpkper=kstpkper, totim=totim).transpose()
+        ).transpose()
+
+        if attrib_name is None:
+            attrib_name = self.text.decode()
+
+        for ix, arr in enumerate(array):
+            name = f"{attrib_name}_{ix}"
+            gdf[name] = np.ravel(arr)
+
+        return gdf
+
     def to_shapefile(
         self,
         filename: Union[str, os.PathLike],
@@ -287,7 +335,10 @@ def to_shapefile(
         >>> times = hdobj.get_times()
         >>> hdobj.to_shapefile('test_heads_sp6.shp', totim=times[-1])
         """
-
+        warnings.warn(
+            "to_shapefile() is deprecated and is being replaced by to_geo_dataframe()",
+            DeprecationWarning,
+        )
         plotarray = np.atleast_3d(
             self.get_data(kstpkper=kstpkper, totim=totim, mflay=mflay).transpose()
         ).transpose()
@@ -301,7 +352,7 @@ def to_shapefile(
 
         from ..export.shapefile_utils import write_grid_shapefile
 
-        write_grid_shapefile(filename, self.mg, attrib_dict, verbose)
+        write_grid_shapefile(filename, self.modelgrid, attrib_dict, verbose)
 
     def plot(
         self,
@@ -413,7 +464,7 @@ def plot(
             axes=axes,
             filenames=filenames,
             mflay=mflay,
-            modelgrid=self.mg,
+            modelgrid=self.modelgrid,
             **kwargs,
         )
 
diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py
index 6fbb930d2..938206de8 100644
--- a/flopy/utils/modpathfile.py
+++ b/flopy/utils/modpathfile.py
@@ -12,6 +12,7 @@
 from flopy.utils.particletrackfile import ParticleTrackFile
 
 from ..utils.flopy_io import loadtxt
+from ..utils.utl_import import import_optional_dependency
 
 
 class ModpathFile(ParticleTrackFile):
@@ -671,6 +672,61 @@ def get_destination_endpoint_data(self, dest_cells, source=False):
         inds = np.isin(raslice, dest_cells)
         return data[inds].copy().view(np.recarray)
 
+    def to_geo_dataframe(
+        self,
+        modelgrid,
+        data=None,
+        direction="ending",
+    ):
+        """
+        Create a geodataframe of particle starting / ending locations.
+
+        Parameters
+        ----------
+        modelgrid : flopy.discretization.grid instance
+            Used to scale and rotate Global x,y,z values in MODPATH Endpoint
+            file.
+        data : np.recarray
+            Record array of same form as that returned by EndpointFile.get_alldata.
+            (if none, EndpointFile.get_alldata() is exported).
+        direction : str
+            String defining if 'starting' or 'ending' particle locations should be
+            considered. (default is 'ending')
+        """
+        from ..utils import geometry
+
+        gpd = import_optional_dependency("geopandas")
+        shapely_geo = import_optional_dependency("shapely.geometry")
+        if data is None:
+            data = self.get_alldata()
+
+        if direction.lower() == "ending":
+            xcol, ycol, zcol = "x", "y", "z"
+        elif direction.lower() == "starting":
+            xcol, ycol, zcol = "x0", "y0", "z0"
+        else:
+            raise Exception(
+                'flopy.map.plot_endpoint direction must be "ending" or "starting".'
+            )
+        x, y = geometry.transform(
+            data[xcol],
+            data[ycol],
+            xoff=modelgrid.xoffset,
+            yoff=modelgrid.yoffset,
+            angrot_radians=modelgrid.angrot_radians,
+        )
+        z = data[zcol]
+
+        geoms = [shapely_geo.Point(p) for p in zip(x, y, z)]
+        gdf = gpd.GeoDataFrame(data, geometry=geoms, crs=modelgrid.crs)
+
+        # adjust to 1 based node numbers
+        for col in list(gdf):
+            if col in self.kijnames:
+                gdf[col] += 1
+
+        return gdf
+
     def write_shapefile(
         self,
         data=None,
@@ -715,44 +771,21 @@ def write_shapefile(
                 - ``epsg`` (int): use ``crs`` instead.
 
         """
-        from ..discretization import StructuredGrid
-        from ..export.shapefile_utils import recarray2shp
-        from ..utils import geometry
-        from ..utils.geometry import Point
+        import warnings
 
+        warnings.warn(
+            "write_shapefile is Deprecated, please use to_geo_dataframe() in the future"
+        )
         epd = (data if data is not None else endpoint_data).copy()
-        if epd is None:
-            epd = self.get_alldata()
+        gdf = self.to_geo_dataframe(modelgrid=mg, data=epd, direction=direction)
 
-        if direction.lower() == "ending":
-            xcol, ycol, zcol = "x", "y", "z"
-        elif direction.lower() == "starting":
-            xcol, ycol, zcol = "x0", "y0", "z0"
-        else:
-            raise Exception(
-                'flopy.map.plot_endpoint direction must be "ending" or "starting".'
-            )
-        if mg is None:
-            raise ValueError("A modelgrid object was not provided.")
-
-        if isinstance(mg, StructuredGrid):
-            x, y = geometry.transform(
-                epd[xcol],
-                epd[ycol],
-                xoff=mg.xoffset,
-                yoff=mg.yoffset,
-                angrot_radians=mg.angrot_radians,
-            )
-        else:
-            x, y = mg.get_coords(epd[xcol], epd[ycol])
-        z = epd[zcol]
+        if crs is not None:
+            if gdf.crs is None:
+                gdf = gdf.set_crs(crs)
+            else:
+                gdf = gdf.to_crs(crs)
 
-        geoms = [Point(x[i], y[i], z[i]) for i in range(len(epd))]
-        # convert back to one-based
-        for n in self.kijnames:
-            if n in epd.dtype.names:
-                epd[n] += 1
-        recarray2shp(epd, geoms, shpname=shpname, crs=crs, **kwargs)
+        gdf.to_file(shpname)
 
 
 class TimeseriesFile(ModpathFile):
diff --git a/flopy/utils/particletrackfile.py b/flopy/utils/particletrackfile.py
index c5df347c0..483ce1570 100644
--- a/flopy/utils/particletrackfile.py
+++ b/flopy/utils/particletrackfile.py
@@ -10,6 +10,8 @@
 import numpy as np
 from numpy.lib.recfunctions import stack_arrays
 
+from .utl_import import import_optional_dependency
+
 MIN_PARTICLE_TRACK_DTYPE = np.dtype(
     [
         ("x", np.float32),
@@ -183,6 +185,125 @@ def intersect(self, cells, to_recarray) -> np.recarray:
         """Find intersection of pathlines with cells."""
         pass
 
+    def to_geo_dataframe(
+        self,
+        modelgrid,
+        data=None,
+        one_per_particle=True,
+        direction="ending",
+    ):
+        """
+        Create a geodataframe of particle tracks.
+
+        Parameters
+        ----------
+        modelgrid : flopy.discretization.Grid instance
+            Used to scale and rotate Global x, y, z values.
+        data : np.recarray
+            Record array of same form as that returned by
+            get_alldata(). (if none, get_alldata() is exported).
+        one_per_particle : boolean (default True)
+            True writes a single LineString with a single set of attribute
+            data for each particle. False writes a record/geometry for each
+            pathline segment (each row in the PathLine file). This option can
+            be used to visualize attribute information (time, model layer,
+            etc.) across a pathline in a GIS.
+        direction : str
+            String defining if starting or ending particle locations should be
+            included in shapefile attribute information. Only used if
+            one_per_particle=False. (default is 'ending')
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from . import geometry
+
+        shapely_geo = import_optional_dependency("shapely.geometry")
+        gpd = import_optional_dependency("geopandas")
+
+        if data is None:
+            data = self._data.view(np.recarray)
+        else:
+            # convert pathline list to a single recarray
+            if isinstance(data, list):
+                s = data[0]
+                print(s.dtype)
+                for n in range(1, len(data)):
+                    s = stack_arrays((s, data[n]))
+                data = s.view(np.recarray)
+
+        data = data.copy()
+        data.sort(order=["particleid", "time"])
+
+        particles = np.unique(data.particleid)
+        geoms = []
+
+        # create a dict of attrs?
+        headings = ["particleid", "particlegroup", "time", "k", "i", "j", "node"]
+        attrs = []
+        for h in headings:
+            if h in data.dtype.names:
+                attrs.append(h)
+
+        if one_per_particle:
+            dfdata = {a: [] for a in attrs}
+            if direction == "ending":
+                idx = -1
+            else:
+                idx = 0
+
+            for p in particles:
+                ra = data[data.particleid == p]
+                for k in dfdata.keys():
+                    if k == "time":
+                        dfdata[k].append(np.max(ra[k]))
+                    else:
+                        dfdata[k].append(ra[k][idx])
+
+                x, y = geometry.transform(
+                    ra.x,
+                    ra.y,
+                    modelgrid.xoffset,
+                    modelgrid.yoffset,
+                    modelgrid.angrot_radians,
+                )
+                z = ra.z
+
+                line = list(zip(x, y, z))
+                geoms.append(shapely_geo.LineString(line))
+
+        else:
+            dfdata = {a: [] for a in data.dtype.names}
+            for pid in particles:
+                ra = data[data.particleid == pid]
+                x, y = geometry.transform(
+                    ra.x,
+                    ra.y,
+                    modelgrid.xoffset,
+                    modelgrid.yoffset,
+                    modelgrid.angrot_radians,
+                )
+                z = ra.z
+                geoms += [
+                    shapely_geo.LineString(
+                        [(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])]
+                    )
+                    for i in np.arange(1, (len(ra)))
+                ]
+                for k in dfdata.keys():
+                    dfdata[k].extend(ra[k][1:])
+
+        # now create a geodataframe
+        gdf = gpd.GeoDataFrame(dfdata, geometry=geoms, crs=modelgrid.crs)
+
+        # adjust to 1 based node numbers
+        for col in list(gdf):
+            if col in self.kijnames:
+                gdf[col] += 1
+
+        return gdf
+
     def write_shapefile(
         self,
         data=None,
@@ -227,102 +348,22 @@ def write_shapefile(
                - ``epsg`` (int): use ``crs`` instead.
 
         """
-        from ..export.shapefile_utils import recarray2shp
-        from . import geometry
-        from .geometry import LineString
-
-        series = data
-        if series is None:
-            series = self._data.view(np.recarray)
-        else:
-            # convert pathline list to a single recarray
-            if isinstance(series, list):
-                s = series[0]
-                print(s.dtype)
-                for n in range(1, len(series)):
-                    s = stack_arrays((s, series[n]))
-                series = s.view(np.recarray)
-
-        series = series.copy()
-        series.sort(order=["particleid", "time"])
+        import warnings
 
-        if mg is None:
-            raise ValueError("A modelgrid object was not provided.")
-
-        particles = np.unique(series.particleid)
-        geoms = []
-
-        # create dtype with select attributes in pth
-        names = series.dtype.names
-        dtype = []
-        atts = ["particleid", "particlegroup", "time", "k", "i", "j", "node"]
-        for att in atts:
-            if att in names:
-                t = np.int32
-                if att == "time":
-                    t = np.float32
-                dtype.append((att, t))
-        dtype = np.dtype(dtype)
-
-        # reset names to the selected names in the created dtype
-        names = dtype.names
-
-        # 1 geometry for each path
-        if one_per_particle:
-            loc_inds = 0
-            if direction == "ending":
-                loc_inds = -1
-
-            sdata = []
-            for pid in particles:
-                ra = series[series.particleid == pid]
-
-                x, y = geometry.transform(
-                    ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians
-                )
-                z = ra.z
-                geoms.append(LineString(list(zip(x, y, z))))
-
-                t = [pid]
-                if "particlegroup" in names:
-                    t.append(ra.particlegroup[0])
-                t.append(ra.time.max())
-                if "k" in names:
-                    t.append(ra.k[loc_inds])
-                if "node" in names:
-                    t.append(ra.node[loc_inds])
-                else:
-                    if "i" in names:
-                        t.append(ra.i[loc_inds])
-                    if "j" in names:
-                        t.append(ra.j[loc_inds])
-                sdata.append(tuple(t))
-
-            sdata = np.array(sdata, dtype=dtype).view(np.recarray)
-
-        # geometry for each row in PathLine file
-        else:
-            dtype = series.dtype
-            sdata = []
-            for pid in particles:
-                ra = series[series.particleid == pid]
-                if mg is not None:
-                    x, y = geometry.transform(
-                        ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians
-                    )
-                else:
-                    x, y = geometry.transform(ra.x, ra.y, 0, 0, 0)
-                z = ra.z
-                geoms += [
-                    LineString([(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])])
-                    for i in np.arange(1, (len(ra)))
-                ]
-                sdata += ra[1:].tolist()
-            sdata = np.array(sdata, dtype=dtype).view(np.recarray)
-
-        # convert back to one-based
-        for n in set(self.kijnames).intersection(set(sdata.dtype.names)):
-            sdata[n] += 1
+        warnings.warn(
+            "write_shapefile will be Deprecated, please use to_geo_dataframe()",
+            DeprecationWarning,
+        )
+        gdf = self.to_geo_dataframe(
+            modelgrid=mg,
+            data=data,
+            one_per_particle=one_per_particle,
+            direction=direction,
+        )
+        if crs is not None:
+            if gdf.crs is None:
+                gdf = gdf.set_crs(crs)
+            else:
+                gdf = gdf.to_crs(crs)
 
-        # write the final recarray to a shapefile
-        recarray2shp(sdata, geoms, shpname=shpname, crs=crs, **kwargs)
+        gdf.to_file(shpname)
diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py
index b83025a90..2e4175209 100644
--- a/flopy/utils/util_array.py
+++ b/flopy/utils/util_array.py
@@ -633,6 +633,45 @@ def data_type(self):
     def plottable(self):
         return True
 
+    def to_geo_dataframe(self, gdf=None, forgive=False, name=None, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        forgive : bool
+            optional flag to continue running and pass data that is not compatible
+            with the geodataframe shape
+        name : str
+            optional array name base. If None, method uses the .name attribute
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.geo_dataframe
+
+            if name is not None:
+                names = [name for _ in range(len(self.util_2ds))]
+            else:
+                names = self.name
+            for lay, u2d in enumerate(self.util_2ds):
+                name = f"{names[lay]}_{lay}"
+                gdf = u2d.to_geo_dataframe(gdf=gdf, name=name, forgive=forgive)
+
+            return gdf
+
     def export(self, f, **kwargs):
         from .. import export
 
@@ -1058,6 +1097,31 @@ def data_type(self):
     def plottable(self):
         return False
 
+    def to_geo_dataframe(self, gdf=None, kper=0, forgive=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        forgive : bool
+            boolean flag for sparse dataframe construction. Default is False
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        u3d = self.transient_3ds[kper]
+        # note: may need to provide a pass through name for u3d to avoid s.p.
+        # number being tacked on. Test this on a model with the LAK package...
+        name = self.name_base[:-1].lower()
+        gdf = u3d.to_geo_dataframe(gdf=gdf, forgive=forgive, name=name, **kwargs)
+        return gdf
+
     def get_zero_3d(self, kper):
         name = f"{self.name_base}{kper + 1}(filled zero)"
         return Util3d(
@@ -1398,6 +1462,29 @@ def from_4d(cls, model, pak_name, m4ds):
             name=name,
         )
 
+    def to_geo_dataframe(self, gdf=None, kper=0, forgive=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        forgive : bool
+            boolean flag for sparse dataframe construction. Default is False
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        u2d = self.transient_2ds[kper]
+        name = self.name_base[:-1]
+        gdf = u2d.to_geo_dataframe(gdf=gdf, name=name, forgive=forgive, **kwargs)
+        return gdf
+
     def __setattr__(self, key, value):
         if hasattr(self, "transient_2ds") and key == "cnstnt":
             # set cnstnt for each u2d
@@ -1875,6 +1962,56 @@ def _decide_how(self):
         else:
             self._how = "internal"
 
+    def to_geo_dataframe(self, gdf=None, name=None, forgive=False, **kwargs):
+        """
+        Method to add an input array to a geopandas GeoDataFrame
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame object
+        name : str
+            optional attribute name, default uses util2d name
+        forgive : bool
+            optional flag to continue if data shape not compatible with GeoDataFrame
+
+        Returns
+        -------
+            geopandas GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if gdf is None:
+                if modelgrid is None:
+                    return gdf
+                gdf = modelgrid.geo_dataframe
+
+            if modelgrid is not None:
+                if modelgrid.grid_type != "unstructured":
+                    ncpl = modelgrid.ncpl
+                else:
+                    ncpl = modelgrid.nnodes
+            else:
+                ncpl = len(gdf)
+
+            if name is None:
+                name = self.name
+
+            data = self.array
+
+            if data.size == ncpl:
+                gdf[name] = data.ravel()
+            elif forgive:
+                return gdf
+            else:
+                raise AssertionError(
+                    f"Data size {data.size} not compatible with dataframe length {ncpl}"
+                )
+
+            return gdf
+
     def plot(
         self,
         title=None,
diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py
index da88d8e1b..3b99bc86f 100644
--- a/flopy/utils/util_list.py
+++ b/flopy/utils/util_list.py
@@ -123,6 +123,56 @@ def get_empty(self, ncell=0):
         d = create_empty_recarray(ncell, self.dtype, default_value=-1.0e10)
         return d
 
+    def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.geo_dataframe
+
+            data = self.array
+
+            col_names = []
+            for name, array4d in data.items():
+                aname = f"{self.name[0].lower()}_{name}"
+                array = array4d[kper]
+                if modelgrid.grid_type == "unstructured":
+                    array = array.ravel()
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array[lay].ravel()
+                        gdf[f"{aname}_{lay}"] = arr.ravel()
+                        col_names.append(f"{aname}_{lay}")
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def export(self, f, **kwargs):
         from .. import export