diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py
index 2a27183289..9f2ad3f1e0 100644
--- a/autotest/test_binaryfile.py
+++ b/autotest/test_binaryfile.py
@@ -1,4 +1,3 @@
-import os
 from itertools import repeat
 
 import numpy as np
@@ -14,7 +13,6 @@
     CellBudgetFile,
     HeadFile,
     HeadUFile,
-    UcnFile,
     Util2d,
 )
 from flopy.utils.binaryfile import (
diff --git a/autotest/test_mbase.py b/autotest/test_mbase.py
index 8a2eddcadc..8e3ee5d828 100644
--- a/autotest/test_mbase.py
+++ b/autotest/test_mbase.py
@@ -10,7 +10,6 @@
 from flopy.mbase import resolve_exe
 from flopy.utils.flopy_io import relpath_safe
 
-
 _system = system()
 
 
diff --git a/autotest/test_particledata.py b/autotest/test_particledata.py
index 1f5bdc85fd..8840ef42da 100644
--- a/autotest/test_particledata.py
+++ b/autotest/test_particledata.py
@@ -1,8 +1,51 @@
+from functools import reduce
+from itertools import chain
+
+import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
+import pytest
+from autotest.test_grid_cases import GridCases
+
+import flopy
+from flopy.discretization import StructuredGrid
+from flopy.mf6.modflow.mfsimulation import MFSimulation
+from flopy.modflow.mf import Modflow
+from flopy.modflow.mfdis import ModflowDis
+from flopy.modpath import (
+    CellDataType,
+    FaceDataType,
+    LRCParticleData,
+    Modpath7,
+    Modpath7Bas,
+    Modpath7Sim,
+    NodeParticleData,
+    ParticleData,
+    ParticleGroupNodeTemplate,
+)
+from flopy.utils.modpathfile import PathlineFile
+
+# utilities
+
+
+def get_nn(grid: StructuredGrid, k, i, j):
+    return k * grid.nrow * grid.ncol + i * grid.ncol + j
+
+
+def flatten(a):
+    return [
+        [
+            *chain.from_iterable(
+                xx if isinstance(xx, tuple) else [xx] for xx in x
+            )
+        ]
+        for x in a
+    ]
+
+
+# test constructors
 
-from flopy.modpath import ParticleData
 
-structured_plocs = [(1, 1, 1), (1, 1, 2)]
 structured_dtype = np.dtype(
     [
         ("k", "<i4"),
@@ -15,37 +58,696 @@
         ("drape", "<i4"),
     ]
 )
-structured_array = np.core.records.fromrecords(
-    [
-        (1, 1, 1, 0.5, 0.5, 0.5, 0.0, 0),
-        (1, 1, 2, 0.5, 0.5, 0.5, 0.0, 0),
-    ],
-    dtype=structured_dtype,
-)
 
 
-def test_particledata_structured_partlocs_as_list_of_tuples():
-    locs = structured_plocs
+def test_particledata_structured_ctor_with_partlocs_as_list_of_tuples():
+    locs = [(0, 1, 1), (0, 1, 2)]
     data = ParticleData(partlocs=locs, structured=True)
 
     assert data.particlecount == 2
     assert data.dtype == structured_dtype
-    assert np.array_equal(data.particledata, structured_array)
+    assert isinstance(data.particledata, pd.DataFrame)
+    assert np.array_equal(
+        data.particledata.to_records(index=False),
+        np.core.records.fromrecords(
+            [
+                (0, 1, 1, 0.5, 0.5, 0.5, 0.0, 0),
+                (0, 1, 2, 0.5, 0.5, 0.5, 0.0, 0),
+            ],
+            dtype=structured_dtype,
+        ),
+    )
 
 
-def test_particledata_structured_partlocs_as_ndarray():
-    locs = np.array(structured_plocs)
+def test_particledata_structured_ctor_with_partlocs_as_ndarray():
+    locs = np.array([(0, 1, 1), (0, 1, 2)])
     data = ParticleData(partlocs=locs, structured=True)
 
     assert data.particlecount == 2
     assert data.dtype == structured_dtype
-    assert np.array_equal(data.particledata, structured_array)
+    assert isinstance(data.particledata, pd.DataFrame)
+    assert np.array_equal(
+        data.particledata.to_records(index=False),
+        np.core.records.fromrecords(
+            [
+                (0, 1, 1, 0.5, 0.5, 0.5, 0.0, 0),
+                (0, 1, 2, 0.5, 0.5, 0.5, 0.0, 0),
+            ],
+            dtype=structured_dtype,
+        ),
+    )
 
 
-def test_particledata_structured_partlocs_as_list_of_lists():
-    locs = [list(p) for p in structured_plocs]
+def test_particledata_structured_ctor_with_partlocs_as_list_of_lists():
+    locs = [list(p) for p in [(0, 1, 1), (0, 1, 2)]]
     data = ParticleData(partlocs=locs, structured=True)
 
     assert data.particlecount == 2
     assert data.dtype == structured_dtype
-    assert np.array_equal(data.particledata, structured_array)
+    assert isinstance(data.particledata, pd.DataFrame)
+    assert np.array_equal(
+        data.particledata.to_records(index=False),
+        np.core.records.fromrecords(
+            [
+                (0, 1, 1, 0.5, 0.5, 0.5, 0.0, 0),
+                (0, 1, 2, 0.5, 0.5, 0.5, 0.0, 0),
+            ],
+            dtype=structured_dtype,
+        ),
+    )
+
+
+# test to_prp()
+
+
+def test_particledata_to_prp_dis_1():
+    # model grid
+    grid = GridCases().structured_small()
+
+    # particle data
+    cells = [(0, 1, 1), (0, 1, 2)]
+    part_data = ParticleData(partlocs=cells, structured=True)
+
+    # convert to global coordinates
+    rpts_prt = flatten(list(part_data.to_prp(grid)))
+
+    # check conversion
+    assert len(rpts_prt) == len(cells)
+    assert all(
+        len(c) == 7 for c in rpts_prt
+    )  # each coord should be a tuple (irpt, k, i, j, x, y, z)
+
+    # expected
+    exp = np.core.records.fromrecords(
+        [
+            (0, 1, 1, 0.5, 0.5, 0.5, 0.0, 0),
+            (0, 1, 2, 0.5, 0.5, 0.5, 0.0, 0),
+        ],
+        dtype=structured_dtype,
+    )
+
+    for ci, cell in enumerate(cells):
+        # check containing cell is correct
+        rpt = rpts_prt[ci]
+        assert cell == grid.intersect(*rpt[4:7])
+
+        # check global coords are equivalent to local
+        k, i, j = cell
+        verts = grid.get_cell_vertices(i, j)
+        xs, ys = list(zip(*verts))
+        minx, maxx = min(xs), max(xs)
+        miny, maxy = min(ys), max(ys)
+        minz, maxz = grid.botm[k, i, j], grid.top[i, j]
+        assert np.isclose(rpt[4], minx + (exp[ci][3] * (maxx - minx)))
+        assert np.isclose(rpt[5], miny + (exp[ci][4] * (maxy - miny)))
+        assert np.isclose(rpt[6], minz + (exp[ci][5] * (maxz - minz)))
+
+
+def test_particledata_to_prp_dis_9():
+    # minimal structured grid
+    grid = GridCases().structured_small()
+
+    # release points in mp7 format (using cell-local coordinates)
+    rpts_mp7 = [
+        # k, i, j, localx, localy, localz
+        # (0-based indexing converted to 1-based for mp7 by flopy)
+        (0, 0, 0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5)
+        for i in range(9)
+    ]
+
+    # create particle data
+    part_data = ParticleData(
+        partlocs=[p[:3] for p in rpts_mp7],
+        structured=True,
+        localx=[p[3] for p in rpts_mp7],
+        localy=[p[4] for p in rpts_mp7],
+        localz=[p[5] for p in rpts_mp7],
+        timeoffset=0,
+        drape=0,
+    )
+
+    # expected release points in PRT format
+    rpts_exp = [
+        # release point index, k, i, j, x, y,
+        [0, 0, 0, 0, 0.10000000149011612, 1.1000000014901161, 7.5],
+        [1, 0, 0, 0, 0.20000000298023224, 1.2000000029802322, 7.5],
+        [2, 0, 0, 0, 0.30000001192092896, 1.300000011920929, 7.5],
+        [3, 0, 0, 0, 0.4000000059604645, 1.4000000059604645, 7.5],
+        [4, 0, 0, 0, 0.5, 1.5, 7.5],
+        [5, 0, 0, 0, 0.6000000238418579, 1.600000023841858, 7.5],
+        [6, 0, 0, 0, 0.699999988079071, 1.699999988079071, 7.5],
+        [7, 0, 0, 0, 0.800000011920929, 1.800000011920929, 7.5],
+        [8, 0, 0, 0, 0.8999999761581421, 1.899999976158142, 7.5],
+    ]
+
+    # convert to prt format
+    rpts_prt = flatten(list(part_data.to_prp(grid)))
+    assert np.allclose(rpts_prt, rpts_exp, atol=1e-3)
+
+
+@pytest.mark.parametrize("localx", [None, 0.5, 0.25])
+@pytest.mark.parametrize("localy", [None, 0.5, 0.25])
+def test_particledata_to_prp_disv_1(localx, localy):
+    """
+    1 particle in bottom left cell, testing with default
+    location (middle), explicitly specifying middle, and
+    offset in x and y directions
+    """
+
+    # model grid
+    grid = GridCases().vertex_small()
+
+    # particle data
+    locs = [4]
+    localx = [localx] if localx else None
+    localy = [localy] if localy else None
+    part_data = ParticleData(
+        partlocs=locs,
+        structured=False,
+        particleids=range(len(locs)),
+        localx=localx,
+        localy=localy,
+    )
+
+    # convert to global coordinates
+    rpts_prt = flatten(list(part_data.to_prp(grid)))
+
+    # check conversion succeeded
+    assert len(rpts_prt) == len(locs)
+    assert all(
+        len(c) == 6 for c in rpts_prt
+    )  # each coord should be a tuple (irpt, k, j, x, y, z)
+    for ci, c in enumerate(rpts_prt):
+        assert np.isclose(c[3], localx[0] if localx else 0.5)  # check x
+        assert np.isclose(c[4], localy[0] if localy else 0.5)  # check y
+        assert np.isclose(c[5], 7.5)  # check z
+
+    # debugging: plot grid, cell centers, and particle location
+    # fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 10))
+    # ax.set_aspect("equal")
+    # grid.plot() # plot grid
+    # xc, yc = ( # plot cell centers
+    #     grid.get_xcellcenters_for_layer(0),
+    #     grid.get_ycellcenters_for_layer(0)
+    # )
+    # xc = xc.flatten()
+    # yc = yc.flatten()
+    # for i in range(grid.ncpl):
+    #     x, y = xc[i], yc[i]
+    #     nn = grid.intersect(x, y, 0)[1]
+    #     ax.plot(x, y, "ro")
+    #     ax.annotate(str(nn + 1), (x, y), color="r") # 1-based node numbering
+    # for c in coords:  # plot particle location(s)
+    #     ax.plot(c[0], c[1], "bo")
+    # plt.show()
+
+
+def test_particledata_to_prp_disv_9():
+    # minimal vertex grid
+    grid = GridCases().vertex_small()
+
+    # release points in mp7 format (using cell-local coordinates)
+    rpts_mp7 = [
+        # node number, localx, localy, localz
+        # (0-based indexing converted to 1-based for mp7 by flopy)
+        (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5)
+        for i in range(9)
+    ]
+
+    # create particle data
+    part_data = ParticleData(
+        partlocs=[p[0] for p in rpts_mp7],
+        structured=False,
+        localx=[p[1] for p in rpts_mp7],
+        localy=[p[2] for p in rpts_mp7],
+        localz=[p[3] for p in rpts_mp7],
+        timeoffset=0,
+        drape=0,
+    )
+
+    # expected release points in PRT format; below we will use flopy
+    # to convert from mp7 to prt format and make sure they are equal
+    rpts_exp = [
+        # particle index, k, j, x, y, z
+        # (0-based indexing converted to 1-based for mf6 by flopy)
+        (
+            i,
+            0,
+            0,
+            float(f"0.{i + 1}"),
+            float(f"2.{i + 1}"),
+            (grid.xyzextent[5] - grid.xyzextent[4]) / 2,
+        )
+        for i in range(9)
+    ]
+
+    # convert to prt format
+    rpts_prt = flatten(list(part_data.to_prp(grid)))
+    assert np.allclose(rpts_prt, rpts_exp, atol=1e-3)
+
+
+def test_lrcparticledata_to_prp_divisions_defaults():
+    sd_data = CellDataType()
+    regions = [[0, 0, 1, 0, 1, 1]]
+    part_data = LRCParticleData(
+        subdivisiondata=[sd_data], lrcregions=[regions]
+    )
+    grid = GridCases().structured_small()
+    rpts_prt = flatten(list(part_data.to_prp(grid)))
+    rpts_exp = [
+        [0, 0, 0, 1, 1.166666, 1.166666, 5.833333],
+        [1, 0, 0, 1, 1.166666, 1.166666, 7.5],
+        [2, 0, 0, 1, 1.166666, 1.166666, 9.166666],
+        [3, 0, 0, 1, 1.1666666, 1.5, 5.833333],
+        [4, 0, 0, 1, 1.1666666, 1.5, 7.5],
+        [5, 0, 0, 1, 1.1666666, 1.5, 9.166666],
+        [6, 0, 0, 1, 1.166666, 1.833333, 5.833333],
+        [7, 0, 0, 1, 1.166666, 1.833333, 7.5],
+        [8, 0, 0, 1, 1.166666, 1.833333, 9.166666],
+        [9, 0, 0, 1, 1.5, 1.166666, 5.833333],
+        [10, 0, 0, 1, 1.5, 1.166666, 7.5],
+        [11, 0, 0, 1, 1.5, 1.166666, 9.166666],
+        [12, 0, 0, 1, 1.5, 1.5, 5.833333],
+        [13, 0, 0, 1, 1.5, 1.5, 7.5],
+        [14, 0, 0, 1, 1.5, 1.5, 9.166666],
+        [15, 0, 0, 1, 1.5, 1.833333, 5.833333],
+        [16, 0, 0, 1, 1.5, 1.833333, 7.5],
+        [17, 0, 0, 1, 1.5, 1.833333, 9.166666],
+        [18, 0, 0, 1, 1.833333, 1.166666, 5.833333],
+        [19, 0, 0, 1, 1.833333, 1.166666, 7.5],
+        [20, 0, 0, 1, 1.833333, 1.166666, 9.166666],
+        [21, 0, 0, 1, 1.833333, 1.5, 5.833333],
+        [22, 0, 0, 1, 1.833333, 1.5, 7.5],
+        [23, 0, 0, 1, 1.833333, 1.5, 9.166666],
+        [24, 0, 0, 1, 1.833333, 1.833333, 5.833333],
+        [25, 0, 0, 1, 1.833333, 1.833333, 7.5],
+        [26, 0, 0, 1, 1.833333, 1.833333, 9.166666],
+        [27, 0, 1, 1, 1.166666, 0.166666, 5.833333],
+        [28, 0, 1, 1, 1.166666, 0.166666, 7.5],
+        [29, 0, 1, 1, 1.166666, 0.166666, 9.166666],
+        [30, 0, 1, 1, 1.166666, 0.5, 5.833333],
+        [31, 0, 1, 1, 1.166666, 0.5, 7.5],
+        [32, 0, 1, 1, 1.166666, 0.5, 9.166666],
+        [33, 0, 1, 1, 1.166666, 0.833333, 5.833333],
+        [34, 0, 1, 1, 1.166666, 0.833333, 7.5],
+        [35, 0, 1, 1, 1.166666, 0.833333, 9.166666],
+        [36, 0, 1, 1, 1.5, 0.166666, 5.833333],
+        [37, 0, 1, 1, 1.5, 0.166666, 7.5],
+        [38, 0, 1, 1, 1.5, 0.166666, 9.166666],
+        [39, 0, 1, 1, 1.5, 0.5, 5.833333],
+        [40, 0, 1, 1, 1.5, 0.5, 7.5],
+        [41, 0, 1, 1, 1.5, 0.5, 9.166666],
+        [42, 0, 1, 1, 1.5, 0.833333, 5.833333],
+        [43, 0, 1, 1, 1.5, 0.833333, 7.5],
+        [44, 0, 1, 1, 1.5, 0.833333, 9.166666],
+        [45, 0, 1, 1, 1.833333, 0.166666, 5.833333],
+        [46, 0, 1, 1, 1.833333, 0.166666, 7.5],
+        [47, 0, 1, 1, 1.833333, 0.166666, 9.166666],
+        [48, 0, 1, 1, 1.833333, 0.5, 5.833333],
+        [49, 0, 1, 1, 1.833333, 0.5, 7.5],
+        [50, 0, 1, 1, 1.833333, 0.5, 9.166666],
+        [51, 0, 1, 1, 1.833333, 0.833333, 5.833333],
+        [52, 0, 1, 1, 1.833333, 0.833333, 7.5],
+        [53, 0, 1, 1, 1.833333, 0.833333, 9.166666],
+    ]
+
+    num_cells = reduce(
+        sum,
+        [
+            (lrc[3] - lrc[0] + 1)
+            * (lrc[4] - lrc[1] + 1)
+            * (lrc[5] - lrc[2] + 1)
+            for lrc in regions
+        ],
+    )
+    act_len = len(rpts_prt)
+    exp_len = (
+        num_cells
+        * sd_data.rowcelldivisions
+        * sd_data.columncelldivisions
+        * sd_data.layercelldivisions
+    )
+    assert act_len == exp_len
+    assert np.allclose(rpts_prt, rpts_exp)
+
+
+def test_lrcparticledata_to_prp_divisions_custom():
+    # create particle data
+    rel_minl = 0
+    rel_maxl = 2
+    rel_minr = 2
+    rel_maxr = 3
+    rel_minc = 2
+    rel_maxc = 3
+    cell_data = flopy.modpath.CellDataType(
+        drape=0,
+        rowcelldivisions=5,
+        columncelldivisions=5,
+        layercelldivisions=1,
+    )
+    regions = [[rel_minl, rel_minr, rel_minc, rel_maxl, rel_maxr, rel_maxc]]
+    part_data = flopy.modpath.LRCParticleData(
+        subdivisiondata=[cell_data],
+        lrcregions=[regions],
+    )
+
+    # create grid
+    nlay, nrow, ncol = 3, 21, 20
+    delr = delc = 500.0
+    top = 350.0
+    botm = [220.0, 200.0, 0.0]
+    m = Modflow()
+    dis = ModflowDis(
+        m,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+    )
+    grid = m.modelgrid
+    # todo once StructuredGrid initializer fixed,
+    # remove dummy model/dis and uncomment below
+    # grid = StructuredGrid(
+    #     delc=delc,
+    #     delr=delr,
+    #     top=top,
+    #     botm=botm,
+    #     nlay=nlay,
+    #     nrow=nrow,
+    #     ncol=ncol
+    # )
+
+    # convert to PRP package data
+    rpts_prt = list(part_data.to_prp(grid))
+    rpts_prt = pd.DataFrame(rpts_prt)
+    assert len(rpts_prt) == 300
+
+    # plot release points (debugging)
+    # fig = plt.figure()
+    # ax = fig.add_subplot(projection='3d')
+    # ax.scatter(rpts_prt[4], rpts_prt[5], rpts_prt[6])
+    # plt.show()
+
+
+def test_lrcparticledata_to_prp_top_bottom():
+    rd = 1
+    cd = 1
+    sddata = FaceDataType(
+        horizontaldivisions1=0,
+        verticaldivisions1=0,
+        horizontaldivisions2=0,
+        verticaldivisions2=0,
+        horizontaldivisions3=0,
+        verticaldivisions3=0,
+        horizontaldivisions4=0,
+        verticaldivisions4=0,
+        rowdivisions5=rd,
+        columndivisions5=cd,
+        rowdivisions6=rd,
+        columndivisions6=cd,
+    )
+    lrcregions = [[0, 1, 1, 0, 1, 1]]
+    data = LRCParticleData(subdivisiondata=[sddata], lrcregions=[lrcregions])
+    grid = GridCases().structured_small()
+    rpts_prt = flatten(list(data.to_prp(grid)))
+
+    num_cells = len(
+        [
+            (lrc[3] - lrc[0]) * (lrc[4] - lrc[1]) * (lrc[5] - lrc[2])
+            for lrc in lrcregions
+        ]
+    )
+    assert (
+        len(rpts_prt) == num_cells * rd * cd * 2
+    )  # 1 particle each on top and bottom faces
+
+    # particle should be centered on each face
+    verts = grid.get_cell_vertices(1, 1)
+    xs, ys = list(zip(*verts))
+    for coord in rpts_prt:
+        assert np.isclose(coord[4], np.mean(xs))
+        assert np.isclose(coord[5], np.mean(ys))
+
+    # check elevation
+    assert rpts_prt[0][6] == grid.top_botm[1, 1, 1]
+    assert rpts_prt[1][6] == grid.top_botm[0, 1, 1]
+
+
+def test_lrcparticledata_to_prp_1_per_face():
+    sddata = FaceDataType(
+        horizontaldivisions1=1,
+        verticaldivisions1=1,
+        horizontaldivisions2=1,
+        verticaldivisions2=1,
+        horizontaldivisions3=1,
+        verticaldivisions3=1,
+        horizontaldivisions4=1,
+        verticaldivisions4=1,
+        rowdivisions5=1,
+        columndivisions5=1,
+        rowdivisions6=1,
+        columndivisions6=1,
+    )
+    lrcregions = [[0, 1, 1, 0, 1, 1]]
+    data = LRCParticleData(subdivisiondata=[sddata], lrcregions=[lrcregions])
+    grid = GridCases().structured_small()
+    rpts_prt = flatten(list(data.to_prp(grid)))
+    rpts_exp = [
+        # irpt, k, i, j, x, y, z
+        [0, 0, 1, 1, 1.0, 0.5, 7.5],
+        [1, 0, 1, 1, 2.0, 0.5, 7.5],
+        [2, 0, 1, 1, 1.5, 0.0, 7.5],
+        [3, 0, 1, 1, 1.5, 1.0, 7.5],
+        [4, 0, 1, 1, 1.5, 0.5, 5.0],
+        [5, 0, 1, 1, 1.5, 0.5, 10.0],
+    ]
+    num_cells = len(
+        [
+            (lrc[3] - lrc[0]) * (lrc[4] - lrc[1]) * (lrc[5] - lrc[2])
+            for lrc in lrcregions
+        ]
+    )
+    assert len(rpts_prt) == num_cells * 6  # 1 particle on each face
+    assert np.allclose(rpts_prt, rpts_exp)
+
+
+def test_nodeparticledata_to_prp_disv_defaults(
+    function_tmpdir, example_data_path
+):
+    """
+    This test loads a GWF simulation, runs it, and feeds it to an MP7 simulation
+    to get expected particle release locations. These could be hard-coded but it
+    is not expensive to run the test003 model, and this way we show how expected
+    values are obtained.
+    """
+
+    # create particle data
+    pdat = NodeParticleData()
+
+    # load gwf simulation, switch workspace, write input files, and run
+    sim = MFSimulation.load(
+        sim_ws=example_data_path / "mf6" / "test003_gwfs_disv"
+    )
+    gwf = sim.get_model("gwf_1")
+    grid = gwf.modelgrid
+    gwf_name = "gwf"
+    gwf_ws = function_tmpdir / gwf_name
+    gwf_ws.mkdir()
+    sim.set_sim_path(gwf_ws)
+    sim.write_simulation()
+    sim.run_simulation()
+
+    # create, write and run mp7 simulation to get expected start locations
+    mp7_name = "mp7"
+    mp7_ws = function_tmpdir / mp7_name
+    mp7_ws.mkdir()
+    pg = ParticleGroupNodeTemplate(
+        particlegroupname="G1",
+        particledata=pdat,
+        filename=f"{mp7_name}.sloc",
+    )
+    mp = Modpath7(
+        modelname=mp7_name,
+        flowmodel=gwf,
+        exe_name="mp7",
+        model_ws=mp7_ws,
+    )
+    mpbas = Modpath7Bas(
+        mp,
+        porosity=0.1,
+    )
+    mpsim = Modpath7Sim(
+        mp,
+        simulationtype="pathline",
+        trackingdirection="forward",
+        budgetoutputoption="summary",
+        stoptimeoption="total",
+        particlegroups=[pg],
+    )
+    mp.write_input()
+    mp.run_model()
+
+    # extract particle starting locations computed by mp7
+    mp_pl_file = mp7_ws / f"{mp7_name}.mppth"
+    assert mp_pl_file.exists()
+    plf = PathlineFile(mp_pl_file)
+    pldata = plf.get_alldata()
+    mp7_pls = pd.concat([pd.DataFrame(ra) for ra in pldata])
+    mp7_pls = mp7_pls.sort_values(by=["time", "particleid"]).head(27)
+    mp7_rpts = [
+        [0, r.k, r.x, r.y, r.z] for r in mp7_pls.itertuples()
+    ]  # omit rpt index
+    mp7_rpts.sort()
+
+    # convert particle data to prt format, flatten (remove cell ID tuples),
+    # remove irpt as it is not gauranteed to match, and sort
+    prt_rpts = flatten(list(pdat.to_prp(grid)))
+    prt_rpts = [r[1:] for r in prt_rpts]  #
+    prt_rpts.sort()
+    assert np.allclose(prt_rpts, mp7_rpts)
+
+
+def test_nodeparticledata_to_prp_dis_1_per_face():
+    sddata = FaceDataType(
+        horizontaldivisions1=1,
+        verticaldivisions1=1,
+        horizontaldivisions2=1,
+        verticaldivisions2=1,
+        horizontaldivisions3=1,
+        verticaldivisions3=1,
+        horizontaldivisions4=1,
+        verticaldivisions4=1,
+        rowdivisions5=1,
+        columndivisions5=1,
+        rowdivisions6=1,
+        columndivisions6=1,
+    )
+    grid = GridCases().structured_small()
+    nodes = [get_nn(grid, 0, 1, 1)]
+    data = NodeParticleData(subdivisiondata=sddata, nodes=[nodes])
+
+    rpts = flatten(list(data.to_prp(grid)))
+
+    num_cells = len(nodes)
+    assert len(rpts) == num_cells * 6
+
+
+def test_nodeparticledata_prp_disv_big(function_tmpdir):
+    Lx = 10000.0
+    Ly = 10500.0
+    nlay = 3
+    nrow = 21
+    ncol = 20
+    delr = Lx / ncol
+    delc = Ly / nrow
+    top = 400
+    botm = [220, 200, 0]
+
+    ms = flopy.modflow.Modflow()
+    dis = flopy.modflow.ModflowDis(
+        ms,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+    )
+
+    from flopy.utils.gridgen import Gridgen
+
+    # create Gridgen workspace
+    gridgen_ws = function_tmpdir
+    gridgen_ws.mkdir(parents=True, exist_ok=True)
+
+    # create Gridgen object
+    g = Gridgen(ms.modelgrid, model_ws=gridgen_ws)
+
+    # add polygon for each refinement level
+    outer_polygon = [
+        [(3500, 4000), (3500, 6500), (6000, 6500), (6000, 4000), (3500, 4000)]
+    ]
+    g.add_refinement_features([outer_polygon], "polygon", 1, range(nlay))
+    refshp0 = gridgen_ws / "rf0"
+
+    middle_polygon = [
+        [(4000, 4500), (4000, 6000), (5500, 6000), (5500, 4500), (4000, 4500)]
+    ]
+    g.add_refinement_features([middle_polygon], "polygon", 2, range(nlay))
+    refshp1 = gridgen_ws / "rf1"
+
+    inner_polygon = [
+        [(4500, 5000), (4500, 5500), (5000, 5500), (5000, 5000), (4500, 5000)]
+    ]
+    g.add_refinement_features([inner_polygon], "polygon", 3, range(nlay))
+    refshp2 = gridgen_ws / "rf2"
+
+    # build the grid
+    g.build(verbose=False)
+    grid_props = g.get_gridprops_vertexgrid()
+    disv_props = g.get_gridprops_disv()
+    grid = flopy.discretization.VertexGrid(**grid_props)
+
+    # define particle data
+    coords = (4718.45, 5281.25)
+    cells = grid.intersect(coords[0], coords[1], 0)
+    nodew = grid.ncpl * 2 + cells[0]
+    facedata = flopy.modpath.FaceDataType(
+        drape=0,
+        verticaldivisions1=10,
+        horizontaldivisions1=10,
+        verticaldivisions2=10,
+        horizontaldivisions2=10,
+        verticaldivisions3=10,
+        horizontaldivisions3=10,
+        verticaldivisions4=10,
+        horizontaldivisions4=10,
+        rowdivisions5=0,
+        columndivisions5=0,
+        rowdivisions6=4,
+        columndivisions6=4,
+    )
+    pgdata = flopy.modpath.NodeParticleData(
+        subdivisiondata=facedata, nodes=nodew
+    )
+
+    # convert to PRP package data
+    rpts_prt = flatten(list(pgdata.to_prp(grid)))
+    print(rpts_prt)
+
+
+# test write
+
+
+def test_lrcparticledata_write(function_tmpdir):
+    # create particle data
+    rel_minl = 0
+    rel_maxl = 2
+    rel_minr = 2
+    rel_maxr = 3
+    rel_minc = 2
+    rel_maxc = 3
+    cell_data = flopy.modpath.CellDataType(
+        drape=0,
+        rowcelldivisions=5,
+        columncelldivisions=5,
+        layercelldivisions=1,
+    )
+    regions = [[rel_minl, rel_minr, rel_minc, rel_maxl, rel_maxr, rel_maxc]]
+    part_data = flopy.modpath.LRCParticleData(
+        subdivisiondata=[cell_data],
+        lrcregions=[regions],
+    )
+
+    # write to a file
+    p = function_tmpdir / "f.txt"
+    with open(p, "w") as f:
+        part_data.write(f)
+
+    # check lines written
+    lines = open(p).readlines()
+    assert lines == ["1 1\n", "2 1 0\n", " 5 5 1\n", "1 3 3 3 4 4 \n"]
diff --git a/autotest/test_plot_particle_tracks.py b/autotest/test_plot_particle_tracks.py
index 6116d5fbf7..2a6ba088a0 100644
--- a/autotest/test_plot_particle_tracks.py
+++ b/autotest/test_plot_particle_tracks.py
@@ -172,3 +172,188 @@ def test_plot_map_view_mp6_endpoint(modpath_model):
     )
     # plt.show()
     assert isinstance(ep, PathCollection)
+
+
+# MP7
+
+
+simname = "test_plot"
+nlay = 1
+nrow = 10
+ncol = 10
+top = 1.0
+botm = [0.0]
+nper = 1
+perlen = 1.0
+nstp = 1
+tsmult = 1.0
+porosity = 0.1
+
+
+@pytest.fixture
+def mf6_gwf_sim(module_tmpdir):
+    gwfname = f"{simname}_gwf"
+
+    # create simulation
+    sim = flopy.mf6.MFSimulation(
+        sim_name=simname,
+        exe_name="mf6",
+        version="mf6",
+        sim_ws=module_tmpdir,
+    )
+
+    # create tdis package
+    flopy.mf6.modflow.mftdis.ModflowTdis(
+        sim,
+        pname="tdis",
+        time_units="DAYS",
+        nper=nper,
+        perioddata=[(perlen, nstp, tsmult)],
+    )
+
+    # create gwf model
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True)
+
+    # create gwf discretization
+    flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
+        gwf,
+        pname="dis",
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+    )
+
+    # create gwf initial conditions package
+    flopy.mf6.modflow.mfgwfic.ModflowGwfic(gwf, pname="ic")
+
+    # create gwf node property flow package
+    flopy.mf6.modflow.mfgwfnpf.ModflowGwfnpf(
+        gwf,
+        pname="npf",
+        save_saturation=True,
+        save_specific_discharge=True,
+    )
+
+    # create gwf chd package
+    spd = {
+        0: [[(0, 0, 0), 1.0, 1.0], [(0, 9, 9), 0.0, 0.0]],
+        1: [[(0, 0, 0), 0.0, 0.0], [(0, 9, 9), 1.0, 2.0]],
+    }
+    chd = flopy.mf6.ModflowGwfchd(
+        gwf,
+        pname="CHD-1",
+        stress_period_data=spd,
+        auxiliary=["concentration"],
+    )
+
+    # create gwf output control package
+    gwf_budget_file = f"{gwfname}.bud"
+    gwf_head_file = f"{gwfname}.hds"
+    oc = flopy.mf6.ModflowGwfoc(
+        gwf,
+        budget_filerecord=gwf_budget_file,
+        head_filerecord=gwf_head_file,
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    # create iterative model solution for gwf model
+    ims = flopy.mf6.ModflowIms(sim)
+
+    return sim
+
+
+@pytest.fixture
+def mp7_sim(function_tmpdir, mf6_gwf_sim):
+    pass
+
+
+@pytest.mark.skip(reason="todo")
+@requires_exe("mf6", "mp7")
+def test_plot_map_view_mp7_pathline(mp7_sim):
+    pass
+
+
+@pytest.mark.skip(reason="todo")
+@requires_exe("mf6", "mp7")
+def test_plot_cross_section_mp7_pathline(mp7_sim):
+    pass
+
+
+# MF6 PRT
+
+
+# @pytest.fixture
+# def mf6_prt_sim(function_tmpdir, mf6_gwf_sim):
+#     prtname = f"{simname}_prt"
+#
+#     # create prt model
+#     prt = flopy.mf6.ModflowPrt(mf6_gwf_sim, modelname=prtname)
+#
+#     # create prt discretization
+#     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
+#         prt,
+#         pname="dis",
+#         nlay=nlay,
+#         nrow=nrow,
+#         ncol=ncol,
+#     )
+#
+#     # create mip package
+#     flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+#
+#     # create prp package
+#     flopy.mf6.ModflowPrtprp(
+#         prt,
+#         pname="prp1",
+#         filename=f"{prtname}_1.prp",
+#         nreleasepts=len(releasepts),
+#         packagedata=releasepts,
+#         perioddata={0: ["FIRST"]},
+#     )
+#
+#     # create output control package
+#     flopy.mf6.ModflowPrtoc(
+#         prt,
+#         pname="oc",
+#         track_filerecord=[prt_track_file],
+#         trackcsv_filerecord=[prt_track_csv_file],
+#     )
+#
+#     # create a flow model interface
+#     # todo Fienen's report (crash when FMI created but not needed)
+#     # flopy.mf6.ModflowPrtfmi(
+#     #     prt,
+#     #     packagedata=[
+#     #         ("GWFHEAD", gwf_head_file),
+#     #         ("GWFBUDGET", gwf_budget_file),
+#     #     ],
+#     # )
+#
+#     # create exchange
+#     flopy.mf6.ModflowGwfprt(
+#         sim,
+#         exgtype="GWF6-PRT6",
+#         exgmnamea=gwfname,
+#         exgmnameb=prtname,
+#         filename=f"{gwfname}.gwfprt",
+#     )
+#
+#     # add explicit model solution
+#     ems = flopy.mf6.ModflowEms(
+#         sim,
+#         pname="ems",
+#         filename=f"{prtname}.ems",
+#     )
+#     sim.register_solution_package(ems, [prt.name])
+#
+#     return sim
+#
+#
+# @requires_exe("mf6")
+# def test_plot_map_view_prt_pathline(mf6_prt_sim):
+#     pass
+#
+#
+# @requires_exe("mf6")
+# def test_plot_cross_section_prt_pathline(mf6_prt_sim):
+#     pass
diff --git a/autotest/test_plotutil.py b/autotest/test_plotutil.py
new file mode 100644
index 0000000000..78b68cc710
--- /dev/null
+++ b/autotest/test_plotutil.py
@@ -0,0 +1,291 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from flopy.plot.plotutil import to_mp7_endpoints, to_mp7_pathlines
+
+# test PRT-MP7 pathline conversion functions
+# todo: define fields in a single location and reference from here
+# todo: support optional grid parameter to conversion functions
+# todo: test to_prt_pathlines() conversion and round-trip
+
+
+prt_pl_cols = []
+
+
+mp7_pl_cols = [
+    "particleid",
+    "particlegroup",
+    "sequencenumber",
+    "particleidloc",
+    "time",
+    "xloc",
+    "yloc",
+    "zloc",
+    "x",
+    "y",
+    "z",
+    "node",
+    "k",
+    "stressperiod",
+    "timestep",
+]
+
+
+mp7_ep_cols = [
+    "particleid",
+    "particlegroup",
+    "particleidloc",
+    "time",
+    "time0",
+    "xloc",
+    "xloc0",
+    "yloc",
+    "yloc0",
+    "zloc",
+    "zloc0",
+    "x0",
+    "y0",
+    "z0",
+    "x",
+    "y",
+    "z",
+    "node",
+    "node0",
+    "k",
+    "k0",
+    "zone",
+    "zone0",
+    "initialcellface",
+    "cellface",
+    "status",
+]
+
+
+pls = pd.DataFrame.from_records(
+    [
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            0,
+            0,
+            0,
+            0.0,
+            0.000000,
+            0.100000,
+            9.1,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            0,
+            1,
+            1,
+            0.0,
+            0.063460,
+            0.111111,
+            9.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            11,
+            0,
+            1,
+            1,
+            0.0,
+            0.830431,
+            0.184020,
+            8.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            21,
+            0,
+            1,
+            1,
+            0.0,
+            2.026390,
+            0.267596,
+            7.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            31,
+            0,
+            1,
+            1,
+            0.0,
+            3.704265,
+            0.360604,
+            6.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            60,
+            0,
+            1,
+            1,
+            0.0,
+            39.087992,
+            9.639587,
+            4.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            70,
+            0,
+            1,
+            1,
+            0.0,
+            40.765791,
+            9.732597,
+            3.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            80,
+            0,
+            1,
+            1,
+            0.0,
+            41.961755,
+            9.816110,
+            2.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            90,
+            0,
+            1,
+            1,
+            0.0,
+            42.728752,
+            9.888968,
+            1.0,
+            0.5,
+            "PRP000000001",
+        ],
+        [
+            1,
+            1,
+            1,
+            1,
+            1,
+            1,
+            100,
+            0,
+            5,
+            3,
+            0.0,
+            42.728752,
+            9.888968,
+            1.0,
+            0.5,
+            "PRP000000001",
+        ],
+    ],
+    columns=[
+        "kper",
+        "kstp",
+        "imdl",
+        "iprp",
+        "irpt",
+        "ilay",
+        "icell",
+        "izone",
+        "istatus",
+        "ireason",
+        "trelease",
+        "t",
+        "x",
+        "y",
+        "z",
+        "name",
+    ],
+)
+
+
+@pytest.mark.parametrize("dataframe", [True, False])
+def test_to_mp7_pathlines(dataframe):
+    inp_pls = pls if dataframe else pls.to_records(index=False)
+    mp7_pls = to_mp7_pathlines(inp_pls)
+    assert (
+        type(inp_pls)
+        == type(mp7_pls)
+        == (pd.DataFrame if dataframe else np.recarray)
+    )
+    assert len(mp7_pls) == 10
+    assert set(
+        dict(mp7_pls.dtypes).keys() if dataframe else mp7_pls.dtype.names
+    ) == set(mp7_pl_cols)
+
+
+@pytest.mark.parametrize("dataframe", [True, False])
+def test_to_mp7_endpoints(dataframe):
+    inp_pls = pls if dataframe else pls.to_records(index=False)
+    mp7_eps = to_mp7_endpoints(inp_pls)
+    assert len(mp7_eps) == 1
+    assert set(
+        dict(mp7_eps.dtypes).keys() if dataframe else mp7_eps.dtype.names
+    ) == set(mp7_ep_cols)
diff --git a/autotest/test_util_geometry.py b/autotest/test_util_geometry.py
index a727f7b5d9..2c70c730fe 100644
--- a/autotest/test_util_geometry.py
+++ b/autotest/test_util_geometry.py
@@ -1,7 +1,8 @@
 import numpy as np
 import pytest
+from autotest.test_grid_cases import GridCases
 
-from flopy.utils.geometry import is_clockwise
+from flopy.utils.geometry import is_clockwise, point_in_polygon
 
 
 def test_does_isclockwise_work():
@@ -27,3 +28,62 @@ def test_does_isclockwise_work():
     rslt = is_clockwise(xv, yv)
 
     assert bool(rslt) is False, "is_clockwise() failed"
+
+
+def debug_plot(grid, cell, xpts, ypts, mask):
+    import matplotlib.pyplot as plt
+
+    grid.plot()
+    plt.plot(*list(zip(*cell)))
+    plt.scatter(xpts[mask], ypts[mask], c="red")
+    # plt.show()
+
+
+def test_point_in_polygon_interior():
+    grid = GridCases().structured_small()
+    cell = grid.verts[grid.iverts[0]].tolist()
+    xpts = grid.xcellcenters
+    ypts = grid.ycellcenters
+    mask = point_in_polygon(xpts, ypts, cell)
+    assert mask.sum() == 1
+    assert mask[0, 0] == True
+    debug_plot(grid, cell, xpts, ypts, mask)
+
+
+def test_point_in_polygon_vertices():
+    grid = GridCases().structured_small()
+    cell = grid.verts[grid.iverts[0]].tolist()
+    xpts, ypts = list(zip(*cell))
+    xpts = np.array([xpts])
+    ypts = np.array([ypts])
+    mask = point_in_polygon(xpts, ypts, cell)
+    assert mask.sum() == 1  # only bottom left (due to axis-alignment)
+    debug_plot(grid, cell, xpts, ypts, mask)
+
+    # move points inside boundary
+    xpts[0, 2] = xpts[0, 2] - 0.001
+    ypts[0, 2] = ypts[0, 2] + 0.001
+    xpts[0, 1] = xpts[0, 1] - 0.001
+    ypts[0, 1] = ypts[0, 1] - 0.001
+    xpts[0, 0] = xpts[0, 0] + 0.001
+    ypts[0, 0] = ypts[0, 0] - 0.001
+    mask = point_in_polygon(xpts, ypts, cell)
+    assert mask.sum() == 4  # expect all points
+    debug_plot(grid, cell, xpts, ypts, mask)
+
+
+def test_point_in_polygon_faces():
+    grid = GridCases().structured_small()
+    cell = grid.verts[grid.iverts[0]].tolist()
+    xpts_v, ypts_v = list(zip(*cell))
+    xpts_v = np.array([xpts_v])
+    ypts_v = np.array([ypts_v])
+    xpts = np.array(
+        [[xpts_v[0, 0], xpts_v[0, 2], np.mean(xpts_v), np.mean(xpts_v)]]
+    )
+    ypts = np.array(
+        [[np.mean(ypts_v), np.mean(ypts_v), ypts_v[0, 0], ypts_v[0, 2]]]
+    )
+    mask = point_in_polygon(xpts, ypts, cell)
+    assert mask.sum() == 2  # only inner faces
+    debug_plot(grid, cell, xpts, ypts, mask)
diff --git a/flopy/modpath/mp7particledata.py b/flopy/modpath/mp7particledata.py
index 9f8b8b393d..9d11a00f1c 100644
--- a/flopy/modpath/mp7particledata.py
+++ b/flopy/modpath/mp7particledata.py
@@ -4,8 +4,11 @@
 
 
 """
+from itertools import product
+from typing import Generator, Iterator, Tuple
 
 import numpy as np
+import pandas as pd
 from numpy.lib.recfunctions import unstructured_to_structured
 
 from ..utils.recarray_utils import create_empty_recarray
@@ -314,21 +317,17 @@ def __init__(
         self.particlecount = particledata.shape[0]
         self.particleidoption = particleidoption
         self.locationstyle = locationstyle
-        self.particledata = particledata
-
-        return
+        self.dtype = particledata.dtype
+        self.particledata = pd.DataFrame.from_records(particledata)
 
     def write(self, f=None):
         """
+        Write the particle data to disk.
 
         Parameters
         ----------
         f : fileobject
             Fileobject that is open with write access
-
-        Returns
-        -------
-
         """
         # validate that a valid file object was passed
         if not hasattr(f, "write"):
@@ -338,7 +337,7 @@ def write(self, f=None):
             )
 
         # particle data item 4 and 5
-        d = np.recarray.copy(self.particledata)
+        d = np.recarray.copy(self.particledata.to_records(index=False))
         lnames = [name.lower() for name in d.dtype.names]
         # Add one to the kij and node indices
         for idx in (
@@ -358,7 +357,108 @@ def write(self, f=None):
         for v in d:
             f.write(fmt.format(*v))
 
-        return
+    def to_coords(self, grid) -> Iterator[tuple]:
+        """
+        Compute global particle coordinates on the given grid.
+
+        Parameters
+        ----------
+        grid : The grid on which to locate particle release points.
+
+        Returns
+        -------
+            Generates coordinate tuples (x, y, z)
+        """
+
+        def cvt_xy(p, vs):
+            mn, mx = min(vs), max(vs)
+            span = mx - mn
+            return mn + span * p
+
+        if grid.grid_type == "structured":
+            if not hasattr(self.particledata, "k"):
+                raise ValueError(
+                    f"Particle representation is not structured but grid is"
+                )
+
+            def cvt_z(p, k, i, j):
+                mn, mx = (
+                    grid.botm[k, i, j],
+                    grid.top[i, j] if k == 0 else grid.botm[k - 1, i, j],
+                )
+                span = mx - mn
+                return mn + span * p
+
+            def convert(row) -> Tuple[float, float, float]:
+                verts = grid.get_cell_vertices(row.i, row.j)
+                xs, ys = list(zip(*verts))
+                return [
+                    cvt_xy(row.localx, xs),
+                    cvt_xy(row.localy, ys),
+                    cvt_z(row.localz, row.k, row.i, row.j),
+                ]
+
+        else:
+            if hasattr(self.particledata, "k"):
+                raise ValueError(
+                    f"Particle representation is structured but grid is not"
+                )
+
+            def cvt_z(p, nn):
+                k, j = grid.get_lni([nn])[0]
+                mn, mx = (
+                    grid.botm[k, j],
+                    grid.top[j] if k == 0 else grid.botm[k - 1, j],
+                )
+                span = mx - mn
+                return mn + span * p
+
+            def convert(row) -> Tuple[float, float, float]:
+                verts = grid.get_cell_vertices(row.node)
+                xs, ys = list(zip(*verts))
+                return [
+                    cvt_xy(row.localx, xs),
+                    cvt_xy(row.localy, ys),
+                    cvt_z(row.localz, row.node),
+                ]
+
+        for t in self.particledata.itertuples():
+            yield convert(t)
+
+    def to_prp(self, grid) -> Iterator[tuple]:
+        """
+        Convert particle data to PRT particle release point (PRP)
+        package data entries for the given grid. A model grid is
+        required because MODPATH supports several ways to specify
+        particle release locations by cell ID and subdivision info
+        or local coordinates, but PRT expects global coordinates.
+
+        Parameters
+        ----------
+        grid : The grid on which to locate particle release points.
+
+        Returns
+        -------
+            Generates PRT particle release point (PRP) package
+            data tuples: release point index, k, [i,] j, x, y, z.
+            If the grid is not structured, i is omitted and j is
+            the within-layer cell index for vertex grids.
+        """
+
+        for i, (t, c) in enumerate(
+            zip(
+                self.particledata.itertuples(index=False),
+                self.to_coords(grid),
+            )
+        ):
+            row = [i]  # release point index (irpt)
+            if "node" in self.particledata:
+                k, j = grid.get_lni([t.node])[0]
+                row.extend([(k, j)])
+            else:
+                row.extend([t.k, t.i, t.j])
+            row.extend(c)
+            yield tuple(row)
 
     def _get_dtype(self, structured, particleid):
         """
@@ -421,7 +521,7 @@ def _fmt_string(self):
 
         """
         fmts = []
-        for field in self.particledata.dtype.descr:
+        for field in self.dtype.descr:
             vtype = field[1][1].lower()
             if vtype == "i" or vtype == "b":
                 fmts.append("{:9d}")
@@ -542,7 +642,6 @@ def __init__(
         self.columndivisions5 = columndivisions5
         self.rowdivisions6 = rowdivisions6
         self.columndivisions6 = columndivisions6
-        return
 
     def write(self, f=None):
         """
@@ -637,7 +736,6 @@ def __init__(
         self.columncelldivisions = columncelldivisions
         self.rowcelldivisions = rowcelldivisions
         self.layercelldivisions = layercelldivisions
-        return
 
     def write(self, f=None):
         """
@@ -668,45 +766,206 @@ def write(self, f=None):
         )
         f.write(line)
 
-        return
+
+def get_release_points(subdivisiondata, grid, k=None, i=None, j=None, nn=None):
+    if nn is None and (k is None or i is None or j is None):
+        raise ValueError(
+            f"A cell (node) must be specified by indices (for structured grids) or node number (for vertex/unstructured)"
+        )
+
+    rpts = []
+    template = [k, i, j] if nn is None else [nn]
+
+    # get cell coords and span in each dimension
+    if not (k is None or i is None or j is None):
+        verts = grid.get_cell_vertices(i, j)
+        minz, maxz = grid.botm[k, i, j], grid.top[i, j]
+    elif nn is not None:
+        verts = grid.get_cell_vertices(nn)
+        if grid.grid_type == "structured":
+            k, i, j = grid.get_lrc([nn])[0]
+            minz, maxz = (
+                grid.botm[k, i, j],
+                grid.top[i, j] if k == 0 else grid.botm[k - 1, i, j],
+            )
+        else:
+            k, j = grid.get_lni([nn])[0]
+            minz, maxz = (
+                grid.botm[k, j],
+                grid.top[j] if k == 0 else grid.botm[k - 1, j],
+            )
+    else:
+        raise ValueError(
+            f"A cell (node) must be specified by indices (for structured grids) or node number (for vertex/unstructured)"
+        )
+    xs, ys = list(zip(*verts))
+    minx, maxx = min(xs), max(xs)
+    miny, maxy = min(ys), max(ys)
+    xspan = maxx - minx
+    yspan = maxy - miny
+    zspan = maxz - minz
+
+    if isinstance(subdivisiondata, FaceDataType):
+        # x1 (west)
+        if (
+            subdivisiondata.verticaldivisions1 > 0
+            and subdivisiondata.horizontaldivisions1 > 0
+        ):
+            yincr = yspan / subdivisiondata.horizontaldivisions1
+            ylocs = [
+                (miny + (yincr * 0.5) + (yincr * d))
+                for d in range(subdivisiondata.horizontaldivisions1)
+            ]
+            zincr = zspan / subdivisiondata.verticaldivisions1
+            zlocs = [
+                (minz + (zincr * 0.5) + (zincr * d))
+                for d in range(subdivisiondata.verticaldivisions1)
+            ]
+            prod = list(product(*[ylocs, zlocs]))
+            rpts = rpts + [template + [minx, p[0], p[1]] for p in prod]
+
+        # x2 (east)
+        if (
+            subdivisiondata.verticaldivisions2 > 0
+            and subdivisiondata.horizontaldivisions2 > 0
+        ):
+            yincr = yspan / subdivisiondata.horizontaldivisions2
+            ylocs = [
+                (miny + (yincr * 0.5) + (yincr * d))
+                for d in range(subdivisiondata.horizontaldivisions2)
+            ]
+            zincr = zspan / subdivisiondata.verticaldivisions2
+            zlocs = [
+                (minz + (zincr * 0.5) + (zincr * d))
+                for d in range(subdivisiondata.verticaldivisions2)
+            ]
+            prod = list(product(*[ylocs, zlocs]))
+            rpts = rpts + [template + [maxx, p[0], p[1]] for p in prod]
+
+        # y1 (south)
+        if (
+            subdivisiondata.verticaldivisions3 > 0
+            and subdivisiondata.horizontaldivisions3 > 0
+        ):
+            xincr = xspan / subdivisiondata.horizontaldivisions3
+            xlocs = [
+                (minx + (xincr * 0.5) + (xincr * rd))
+                for rd in range(subdivisiondata.horizontaldivisions3)
+            ]
+            zincr = zspan / subdivisiondata.verticaldivisions3
+            zlocs = [
+                (minz + (zincr * 0.5) + (zincr * d))
+                for d in range(subdivisiondata.verticaldivisions3)
+            ]
+            prod = list(product(*[xlocs, zlocs]))
+            rpts = rpts + [template + [p[0], miny, p[1]] for p in prod]
+
+        # y2 (north)
+        if (
+            subdivisiondata.verticaldivisions4 > 0
+            and subdivisiondata.horizontaldivisions4 > 0
+        ):
+            xincr = xspan / subdivisiondata.horizontaldivisions4
+            xlocs = [
+                (minx + (xincr * 0.5) + (xincr * rd))
+                for rd in range(subdivisiondata.horizontaldivisions4)
+            ]
+            zincr = zspan / subdivisiondata.verticaldivisions4
+            zlocs = [
+                (minz + (zincr * 0.5) + (zincr * d))
+                for d in range(subdivisiondata.verticaldivisions4)
+            ]
+            prod = list(product(*[xlocs, zlocs]))
+            rpts = rpts + [template + [p[0], maxy, p[1]] for p in prod]
+
+        # z1 (bottom)
+        if (
+            subdivisiondata.rowdivisions5 > 0
+            and subdivisiondata.columndivisions5 > 0
+        ):
+            xincr = xspan / subdivisiondata.columndivisions5
+            xlocs = [
+                (minx + (xincr * 0.5) + (xincr * rd))
+                for rd in range(subdivisiondata.columndivisions5)
+            ]
+            yincr = yspan / subdivisiondata.rowdivisions5
+            ylocs = [
+                (miny + (yincr * 0.5) + (yincr * rd))
+                for rd in range(subdivisiondata.rowdivisions5)
+            ]
+            prod = list(product(*[xlocs, ylocs]))
+            rpts = rpts + [template + [p[0], p[1], minz] for p in prod]
+
+        # z2 (top)
+        if (
+            subdivisiondata.rowdivisions6 > 0
+            and subdivisiondata.columndivisions6 > 0
+        ):
+            xincr = xspan / subdivisiondata.columndivisions6
+            xlocs = [
+                (minx + (xincr * 0.5) + (xincr * rd))
+                for rd in range(subdivisiondata.columndivisions6)
+            ]
+            yincr = yspan / subdivisiondata.rowdivisions6
+            ylocs = [
+                (miny + (yincr * 0.5) + (yincr * rd))
+                for rd in range(subdivisiondata.rowdivisions6)
+            ]
+            prod = list(product(*[xlocs, ylocs]))
+            rpts = rpts + [template + [p[0], p[1], maxz] for p in prod]
+    elif isinstance(subdivisiondata, CellDataType):
+        xincr = xspan / subdivisiondata.columncelldivisions
+        xlocs = [
+            (minx + (xincr * 0.5) + (xincr * rd))
+            for rd in range(subdivisiondata.columncelldivisions)
+        ]
+        yincr = yspan / subdivisiondata.rowcelldivisions
+        ylocs = [
+            (miny + (yincr * 0.5) + (yincr * d))
+            for d in range(subdivisiondata.rowcelldivisions)
+        ]
+        zincr = zspan / subdivisiondata.layercelldivisions
+        zlocs = [
+            (minz + (zincr * 0.5) + (zincr * d))
+            for d in range(subdivisiondata.layercelldivisions)
+        ]
+        prod = list(product(*[xlocs, ylocs, zlocs]))
+        rpts = rpts + [template + [p[0], p[1], p[2]] for p in prod]
+    else:
+        raise ValueError(
+            f"Unsupported subdivision data type: {type(subdivisiondata)}"
+        )
+
+    return rpts
 
 
 class LRCParticleData:
     """
-    Layer, row, column particle data template class to create MODPATH 7
-    particle location input style 2 on cell faces (templatesubdivisiontype = 1)
-    and/or in cells (templatesubdivisiontype = 2). Particle locations for this
-    template are specified by layer, row, column regions.
+    MODPATH 7 particle release location template class for particle input style 2.
+    Assigns particles to locations on cell faces (templatesubdivisiontype=1) and/or
+    in cells (templatesubdivisiontype=2) for cells specified by (layer, row, column).
 
     Parameters
     ----------
-    subdivisiondata : FaceDataType, CellDataType or list of FaceDataType
-                      and/or CellDataType types
-        FaceDataType, CellDataType, or a list of FaceDataType and/or
-        CellDataTypes that are used to create one or more particle templates
-        in a particle group. If subdivisiondata is None, a default CellDataType
-        with 27 particles per cell will be created (default is None).
-    lrcregions : list of lists tuples or np.ndarrays
-        Layer, row, column (zero-based) regions with particles created using
-        the specified template parameters. A region is defined as a list/tuple
-        of minlayer, minrow, mincolumn, maxlayer, maxrow, maxcolumn values.
-        If subdivisiondata is a list, a list/tuple or array of layer, row,
-        column regions with the same length as subdivision data must be
-        provided. If lrcregions is None, particles will be placed in
-        the first model cell (default is None).
+    subdivisiondata : FaceDataType, CellDataType or array-like of such, optional
+        Particle template(s) defining how particles are arranged within each cell.
+        If None, defaults to CellDataType with 27 particles per cell (default is None).
+    lrcregions : array-like of array-like, optional
+        0-based regions (minlayer, minrow, mincolumn, maxlayer, maxrow, maxcolumn).
+        If subdivisiondata is array-like, regions must be the same length.
+        If None, particles are placed in the first model cell (default is None).
 
     Examples
     --------
 
     >>> import flopy
-    >>> pg = flopy.modpath.LRCParticleData(lrcregions=[0, 0, 0, 3, 10, 10])
+    >>> pg = flopy.modpath.LRCParticleData(lrcregions=[[0, 0, 0, 3, 10, 10]])
 
     """
 
     def __init__(self, subdivisiondata=None, lrcregions=None):
         """
         Class constructor
-
         """
         self.name = "LRCParticleData"
 
@@ -729,7 +988,7 @@ def __init__(self, subdivisiondata=None, lrcregions=None):
                 )
 
         # validate lrcregions data
-        if isinstance(lrcregions, (list, tuple)):
+        if isinstance(lrcregions, (list, tuple, np.ndarray)):
             # determine if the list or tuple contains lists or tuples
             alllsttup = all(
                 isinstance(el, (list, tuple, np.ndarray)) for el in lrcregions
@@ -780,7 +1039,6 @@ def __init__(self, subdivisiondata=None, lrcregions=None):
         self.totalcellregioncount = totalcellregioncount
         self.subdivisiondata = subdivisiondata
         self.lrcregions = lrcregions
-        return
 
     def write(self, f=None):
         """
@@ -805,46 +1063,112 @@ def write(self, f=None):
         # item 2
         f.write(f"{self.particletemplatecount} {self.totalcellregioncount}\n")
 
-        for sd, lrcregion in zip(self.subdivisiondata, self.lrcregions):
+        for sd, region in zip(self.subdivisiondata, self.lrcregions):
             # item 3
             f.write(
-                f"{sd.templatesubdivisiontype} {lrcregion.shape[0]} {sd.drape}\n"
+                f"{sd.templatesubdivisiontype} {region.shape[0]} {sd.drape}\n"
             )
 
             # item 4 or 5
             sd.write(f)
 
             # item 6
-            for row in lrcregion:
+            for row in region:
                 line = ""
                 for lrc in row:
                     line += f"{lrc + 1} "
                 line += "\n"
                 f.write(line)
 
-        return
+    def to_coords(self, grid) -> Iterator[tuple]:
+        """
+        Compute global particle coordinates on the given grid.
+
+        Parameters
+        ----------
+        grid : The grid on which to locate particle release points.
+
+        Returns
+        -------
+            Generator of coordinate tuples (x, y, z)
+        """
+
+        for region in self.lrcregions:
+            for row in region:
+                mink, mini, minj, maxk, maxi, maxj = row
+                for k in range(mink, maxk + 1):
+                    for i in range(mini, maxi + 1):
+                        for j in range(minj, maxj + 1):
+                            for sd in self.subdivisiondata:
+                                for rpt in get_release_points(
+                                    sd, grid, k, i, j
+                                ):
+                                    yield (rpt[3], rpt[4], rpt[5])
+
+    def to_prp(self, grid) -> Iterator[tuple]:
+        """
+        Convert particle data to PRT particle release point (PRP)
+        package data entries for the given grid. A model grid is
+        required because MODPATH supports several ways to specify
+        particle release locations by cell ID and subdivision info
+        or local coordinates, but PRT expects global coordinates.
+
+        Parameters
+        ----------
+        grid : The grid on which to locate particle release points.
+
+        Returns
+        -------
+            Generates PRT particle release point (PRP) package
+            data tuples: release point index, k, i, j, x, y, z
+        """
+
+        if grid.grid_type != "structured":
+            raise ValueError(
+                f"Particle representation is structured but grid is not"
+            )
+
+        irpt_offset = 0
+        for region in self.lrcregions:
+            for row in region:
+                mink, mini, minj, maxk, maxi, maxj = row
+                for k in range(mink, maxk + 1):
+                    for i in range(mini, maxi + 1):
+                        for j in range(minj, maxj + 1):
+                            for sd in self.subdivisiondata:
+                                for irpt, rpt in enumerate(
+                                    get_release_points(sd, grid, k, i, j)
+                                ):
+                                    assert rpt[0] == k
+                                    assert rpt[1] == i
+                                    assert rpt[2] == j
+                                    yield (
+                                        irpt_offset + irpt,
+                                        k,
+                                        i,
+                                        j,
+                                        rpt[3],
+                                        rpt[4],
+                                        rpt[5],
+                                    )
+                                irpt_offset += irpt + 1
 
 
 class NodeParticleData:
     """
-    Node particle data template class to create MODPATH 7 particle location
-    input style 3 on cell faces (templatesubdivisiontype = 1) and/or in cells
-    (templatesubdivisiontype = 2). Particle locations for this template are
-    specified by nodes.
+    MODPATH 7 particle release location template class for particle input style 3.
+    Assigns particles to locations on cell faces (templatesubdivisiontype=1) and/or
+    in cells (templatesubdivisiontype=2) for cells specified by node number
 
     Parameters
     ----------
-    subdivisiondata : FaceDataType, CellDataType or list of FaceDataType
-                      and/or CellDataType types
-        FaceDataType, CellDataType, or a list of FaceDataType and/or
-        CellDataTypes that are used to create one or more particle templates
-        in a particle group. If subdivisiondata is None, a default CellDataType
-        with 27 particles per cell will be created (default is None).
-    nodes : int, list of ints, tuple of ints, or np.ndarray
-        Nodes (zero-based) with particles created using the specified template
-        parameters. If subdivisiondata is a list, a list of nodes with the same
-        length as subdivision data must be provided. If nodes is None,
-        particles will be placed in the first model cell (default is None).
+    subdivisiondata : FaceDataType, CellDataType array-like of such, optional
+        Particle template(s) defining how particles are arranged within each cell.
+        If None, defaults to CellDataType with 27 particles per cell (default is None).
+    nodes : int or array-like of ints, optional
+        0-based node numbers. If subdivisiondata is array-like, nodes must be array-
+        like of the same length. If None, particles are placed in the first model cell
+        (default is None).
 
     Examples
     --------
@@ -892,10 +1216,10 @@ def __init__(self, subdivisiondata=None, nodes=None):
             if len(nodes.shape) == 1:
                 nodes = nodes.reshape(1, nodes.shape[0])
             # convert to a list of numpy arrays
-            t = []
-            for idx in range(nodes.shape[0]):
-                t.append(np.array(nodes[idx, :], dtype=np.int32))
-            nodes = t
+            nodes = [
+                np.array(nodes[i, :], dtype=np.int32)
+                for i in range(nodes.shape[0])
+            ]
         elif isinstance(nodes, (list, tuple)):
             # convert a single list/tuple to a list of tuples if only one
             # entry in subdivisiondata
@@ -939,7 +1263,6 @@ def __init__(self, subdivisiondata=None, nodes=None):
         self.totalcellcount = totalcellcount
         self.subdivisiondata = subdivisiondata
         self.nodedata = nodes
-        return
 
     def write(self, f=None):
         """
@@ -985,4 +1308,53 @@ def write(self, f=None):
                     line += "\n"
             f.write(line)
 
-        return
+    def to_coords(self, grid) -> Iterator[tuple]:
+        """
+        Compute global particle coordinates on the given grid.
+
+        Parameters
+        ----------
+        grid : The grid on which to locate particle release points.
+
+        Returns
+        -------
+            Generator of coordinate tuples (x, y, z)
+        """
+
+        for sd in self.subdivisiondata:
+            for nd in self.nodedata:
+                rpts = get_release_points(sd, grid, nn=int(nd[0]))
+                for i, rpt in enumerate(rpts):
+                    yield (rpt[1], rpt[2], rpt[3])
+
+    def to_prp(self, grid) -> Iterator[tuple]:
+        """
+        Convert particle data to PRT particle release point (PRP)
+        package data entries for the given grid. A model grid is
+        required because MODPATH supports several ways to specify
+        particle release locations by cell ID and subdivision info
+        or local coordinates, but PRT expects global coordinates.
+
+        Parameters
+        ----------
+        grid : The grid on which to locate particle release points.
+
+        Returns
+        -------
+            Generator of PRT particle release point (PRP) package
+            data tuples: release point index, k, j, x, y, z
+        """
+
+        for sd in self.subdivisiondata:
+            for nd in self.nodedata:
+                rpts = get_release_points(sd, grid, nn=int(nd[0]))
+                for irpt, rpt in enumerate(rpts):
+                    row = [irpt]
+                    if grid.grid_type == "structured":
+                        k, i, j = grid.get_lrc([rpt[0]])[0]
+                        row.extend([k, i, j])
+                    else:
+                        k, j = grid.get_lni([rpt[0]])[0]
+                        row.extend([(k, j)])
+                    row.extend([rpt[1], rpt[2], rpt[3]])
+                    yield tuple(row)
diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py
index 1e50ee69e8..e510f33629 100644
--- a/flopy/plot/crosssection.py
+++ b/flopy/plot/crosssection.py
@@ -6,10 +6,12 @@
 import numpy as np
 import pandas as pd
 from matplotlib.patches import Polygon
+from numpy.lib.recfunctions import stack_arrays
 
 from ..utils import geometry, import_optional_dependency
 from ..utils.geospatial_utils import GeoSpatialUtil
 from . import plotutil
+from .plotutil import to_mp7_endpoints, to_mp7_pathlines
 
 warnings.simplefilter("always", PendingDeprecationWarning)
 
@@ -1054,15 +1056,27 @@ def plot_pathline(
         self, pl, travel_time=None, method="cell", head=None, **kwargs
     ):
         """
-        Plot the MODPATH pathlines
+        Plot particle pathlines. Compatible with MODFLOW 6 PRT particle track
+        data format, or MODPATH 6 or 7 pathline data format.
 
         Parameters
         ----------
-        pl : list of rec arrays or a single rec array
-            rec array or list of rec arrays is data returned from
-            modpathfile PathlineFile get_data() or get_alldata()
-            methods. Data in rec array is 'x', 'y', 'z', 'time',
-            'k', and 'particleid'.
+        pl : list of recarrays or dataframes, or a single recarray or dataframe
+            Particle pathline data. If a list of recarrays or dataframes,
+            each must contain the path of only a single particle. If just
+            one recarray or dataframe, it should contain the paths of all
+            particles. The flopy.utils.modpathfile.PathlineFile.get_data()
+            or get_alldata() return value may be passed directly as this
+            argument.
+
+            For MODPATH 6 or 7 pathlines, columns must include 'x', 'y', 'z',
+            'time', 'k', and 'particleid'. Additional columns are ignored.
+
+            For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z',
+            't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional
+            columns are ignored. Note that MODFLOW 6 PRT does not assign to
+            particles a unique ID, but infers particle identity from 'imdl',
+            'iprp', 'irpt', and 'trelease' combos (i.e. via composite key).
         travel_time : float or str
             travel_time is a travel time selection for the displayed
             pathlines. If a float is passed then pathlines with times
@@ -1086,32 +1100,40 @@ def plot_pathline(
         Returns
         -------
         lc : matplotlib.collections.LineCollection
-
+            The pathlines added to the plot.
         """
 
         from matplotlib.collections import LineCollection
 
-        # make sure pathlines is a list
+        # make sure pl is a list
         if not isinstance(pl, list):
-            pids = np.unique(pl["particleid"])
-            if len(pids) > 1:
-                pl = [pl[pl["particleid"] == pid] for pid in pids]
-            else:
-                pl = [pl]
+            if not isinstance(pl, (np.ndarray, pd.DataFrame)):
+                raise TypeError(
+                    "Pathline data must be a list of recarrays or dataframes, "
+                    f"or a single recarray or dataframe, got {type(pl)}"
+                )
+            pl = [pl]
 
-        # make sure each element in pl is a recarray
+        # convert prt to mp7 format
         pl = [
-            p.to_records(index=False) if isinstance(p, pd.DataFrame) else p
+            to_mp7_pathlines(
+                p.to_records(index=False) if isinstance(p, pd.DataFrame) else p
+            )
             for p in pl
         ]
 
+        # merge pathlines then split on particleid
+        pls = stack_arrays(pl, asrecarray=True, usemask=False)
+        pids = np.unique(pls["particleid"])
+        pl = [pls[pls["particleid"] == pid] for pid in pids]
+
+        # configure plot settings
         marker = kwargs.pop("marker", None)
         markersize = kwargs.pop("markersize", None)
         markersize = kwargs.pop("ms", markersize)
         markercolor = kwargs.pop("markercolor", None)
         markerevery = kwargs.pop("markerevery", 1)
         ax = kwargs.pop("ax", self.ax)
-
         if "colors" not in kwargs:
             kwargs["colors"] = "0.5"
 
@@ -1176,7 +1198,7 @@ def plot_timeseries(
         self, ts, travel_time=None, method="cell", head=None, **kwargs
     ):
         """
-        Plot the MODPATH timeseries.
+        Plot the MODPATH timeseries. Not compatible with MODFLOW 6 PRT.
 
         Parameters
         ----------
@@ -1221,22 +1243,64 @@ def plot_endpoint(
         **kwargs,
     ):
         """
+        Plot particle endpoints. Compatible with MODFLOW 6 PRT particle
+        track data format, or MODPATH 6 or 7 endpoint data format.
 
         Parameters
         ----------
-
+        ep : recarray or dataframe
+            A numpy recarray with the endpoint particle data from the
+            MODPATH endpoint file.
+
+            For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z',
+            't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional
+            columns are ignored. Note that MODFLOW 6 PRT does not assign to
+            particles a unique ID, but infers particle identity from 'imdl',
+            'iprp', 'irpt', and 'trelease' combos (i.e. via composite key).
+        direction : str
+            String defining if starting or ending particle locations should be
+            considered. (default is 'ending')
+        selection : tuple
+            tuple that defines the zero-base layer, row, column location
+            (l, r, c) to use to make a selection of particle endpoints.
+            The selection could be a well location to determine capture zone
+            for the well. If selection is None, all particle endpoints for
+            the user-sepcified direction will be plotted. (default is None)
+        selection_direction : str
+            String defining is a selection should be made on starting or
+            ending particle locations. If selection is not None and
+            selection_direction is None, the selection direction will be set
+            to the opposite of direction. (default is None)
+
+        kwargs : ax, c, s or size, colorbar, colorbar_label, shrink. The
+            remaining kwargs are passed into the matplotlib scatter
+            method. If colorbar is True a colorbar will be added to the plot.
+            If colorbar_label is passed in and colorbar is True then
+            colorbar_label will be passed to the colorbar set_label()
+            method. If shrink is passed in and colorbar is True then
+            the colorbar size will be set using shrink.
 
         Returns
         -------
+        sp : matplotlib.collections.PathCollection
+            The PathCollection added to the plot.
         """
-        ax = kwargs.pop("ax", self.ax)
 
-        # colorbar kwargs
+        # convert ep to recarray if needed
+        if isinstance(ep, pd.DataFrame):
+            ep = ep.to_records(index=False)
+
+        # convert ep from prt to mp7 format if needed
+        if "t" in ep.dtype.names:
+            from .plotutil import to_mp7_endpoints
+
+            ep = to_mp7_endpoints(ep)
+
+        # configure plot settings
+        ax = kwargs.pop("ax", self.ax)
         createcb = kwargs.pop("colorbar", False)
         colorbar_label = kwargs.pop("colorbar_label", "Endpoint Time")
         shrink = float(kwargs.pop("shrink", 1.0))
-
-        # marker kwargs
         s = kwargs.pop("s", np.sqrt(50))
         s = float(kwargs.pop("size", s)) ** 2.0
 
diff --git a/flopy/plot/map.py b/flopy/plot/map.py
index 27a7ab76db..30d278c10d 100644
--- a/flopy/plot/map.py
+++ b/flopy/plot/map.py
@@ -6,9 +6,11 @@
 import pandas as pd
 from matplotlib.collections import LineCollection, PathCollection
 from matplotlib.path import Path
+from numpy.lib.recfunctions import stack_arrays
 
 from ..utils import geometry
 from . import plotutil
+from .plotutil import to_mp7_endpoints, to_mp7_pathlines
 
 warnings.simplefilter("always", PendingDeprecationWarning)
 
@@ -700,7 +702,8 @@ def plot_vector(
 
     def plot_pathline(self, pl, travel_time=None, **kwargs):
         """
-        Plot MODPATH pathlines.
+        Plot particle pathlines. Compatible with MODFLOW 6 PRT particle track
+        data format, or MODPATH 6 or 7 pathline data format.
 
         Parameters
         ----------
@@ -708,11 +711,18 @@ def plot_pathline(self, pl, travel_time=None, **kwargs):
             Particle pathline data. If a list of recarrays or dataframes,
             each must contain the path of only a single particle. If just
             one recarray or dataframe, it should contain the paths of all
-            particles. Pathline data returned from PathlineFile.get_data()
-            or get_alldata() can be passed directly as this argument. Data
-            columns should be 'x', 'y', 'z', 'time', 'k', and 'particleid'
-            at minimum. Additional columns are ignored. The 'particleid'
-            column must be unique to each particle path.
+            particles. The flopy.utils.modpathfile.PathlineFile.get_data()
+            or get_alldata() return value may be passed directly as this
+            argument.
+
+            For MODPATH 6 or 7 pathlines, columns must include 'x', 'y', 'z',
+            'time', 'k', and 'particleid'. Additional columns are ignored.
+
+            For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z',
+            't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional
+            columns are ignored. Note that MODFLOW 6 PRT does not assign to
+            particles a unique ID, but infers particle identity from 'imdl',
+            'iprp', 'irpt', and 'trelease' combos (i.e. via composite key).
         travel_time : float or str
             Travel time selection. If a float, then pathlines with total
             time less than or equal to the given value are plotted. If a
@@ -733,19 +743,29 @@ def plot_pathline(self, pl, travel_time=None, **kwargs):
 
         from matplotlib.collections import LineCollection
 
-        # make sure pathlines is a list
+        # make sure pl is a list
         if not isinstance(pl, list):
-            pids = np.unique(pl["particleid"])
-            if len(pids) > 1:
-                pl = [pl[pl["particleid"] == pid] for pid in pids]
-            else:
-                pl = [pl]
+            if not isinstance(pl, (np.ndarray, pd.DataFrame)):
+                raise TypeError(
+                    "Pathline data must be a list of recarrays or dataframes, "
+                    f"or a single recarray or dataframe, got {type(pl)}"
+                )
+            pl = [pl]
 
+        # convert prt to mp7 format
         pl = [
-            p.to_records(index=False) if isinstance(p, pd.DataFrame) else p
+            to_mp7_pathlines(
+                p.to_records(index=False) if isinstance(p, pd.DataFrame) else p
+            )
             for p in pl
         ]
 
+        # merge pathlines then split on particleid
+        pls = stack_arrays(pl, asrecarray=True, usemask=False)
+        pids = np.unique(pls["particleid"])
+        pl = [pls[pls["particleid"] == pid] for pid in pids]
+
+        # configure layer
         if "layer" in kwargs:
             kon = kwargs.pop("layer")
             if isinstance(kon, bytes):
@@ -758,19 +778,21 @@ def plot_pathline(self, pl, travel_time=None, **kwargs):
         else:
             kon = self.layer
 
+        # configure plot settings
         marker = kwargs.pop("marker", None)
         markersize = kwargs.pop("markersize", None)
         markersize = kwargs.pop("ms", markersize)
         markercolor = kwargs.pop("markercolor", None)
         markerevery = kwargs.pop("markerevery", 1)
         ax = kwargs.pop("ax", self.ax)
-
         if "colors" not in kwargs:
             kwargs["colors"] = "0.5"
 
+        # compose pathlines
         linecol = []
         markers = []
         for p in pl:
+            # filter by travel time
             tp = plotutil.filter_modpath_by_travel_time(p, travel_time)
 
             # transform data!
@@ -781,8 +803,10 @@ def plot_pathline(self, pl, travel_time=None, **kwargs):
                 self.mg.yoffset,
                 self.mg.angrot_radians,
             )
+
             # build polyline array
             arr = np.vstack((x0r, y0r)).T
+
             # select based on layer
             if kon >= 0:
                 kk = p["k"].copy().reshape(p.shape[0], 1)
@@ -790,7 +814,8 @@ def plot_pathline(self, pl, travel_time=None, **kwargs):
                 arr = np.ma.masked_where((kk != kon), arr)
             else:
                 arr = np.ma.asarray(arr)
-            # append line to linecol if there is some unmasked segment
+
+            # append pathline if there are any unmasked segments
             if not arr.mask.all():
                 linecol.append(arr)
                 if not arr.mask.all():
@@ -799,6 +824,7 @@ def plot_pathline(self, pl, travel_time=None, **kwargs):
                         for xy in arr[::markerevery]:
                             if not np.all(xy.mask):
                                 markers.append(xy)
+
         # create line collection
         lc = None
         if len(linecol) > 0:
@@ -815,13 +841,15 @@ def plot_pathline(self, pl, travel_time=None, **kwargs):
                     ms=markersize,
                 )
 
+        # set axis limits
         ax.set_xlim(self.extent[0], self.extent[1])
         ax.set_ylim(self.extent[2], self.extent[3])
+
         return lc
 
     def plot_timeseries(self, ts, travel_time=None, **kwargs):
         """
-        Plot MODPATH timeseries.
+        Plot MODPATH 6 or 7 timeseries. Incompatible with MODFLOW 6 PRT.
 
         Parameters
         ----------
@@ -865,13 +893,20 @@ def plot_endpoint(
         **kwargs,
     ):
         """
-        Plot MODPATH endpoints.
+        Plot particle endpoints. Compatible with MODFLOW 6 PRT particle
+        track data format, or MODPATH 6 or 7 endpoint data format.
 
         Parameters
         ----------
         ep : recarray or dataframe
             A numpy recarray with the endpoint particle data from the
-            MODPATH endpoint file
+            MODPATH endpoint file.
+
+            For MODFLOW 6 PRT pathlines, columns must include 'x', 'y', 'z',
+            't', 'trelease', 'imdl', 'iprp', 'irpt', and 'ilay'. Additional
+            columns are ignored. Note that MODFLOW 6 PRT does not assign to
+            particles a unique ID, but infers particle identity from 'imdl',
+            'iprp', 'irpt', and 'trelease' combos (i.e. via composite key).
         direction : str
             String defining if starting or ending particle locations should be
             considered. (default is 'ending')
@@ -902,6 +937,17 @@ def plot_endpoint(
 
         """
 
+        # convert ep to recarray if needed
+        if isinstance(ep, pd.DataFrame):
+            ep = ep.to_records(index=False)
+
+        # convert ep from prt to mp7 format if needed
+        if "t" in ep.dtype.names:
+            from .plotutil import to_mp7_endpoints
+
+            ep = to_mp7_endpoints(ep)
+
+        # parse selection options
         ax = kwargs.pop("ax", self.ax)
         tep, _, xp, yp = plotutil.parse_modpath_selection_options(
             ep, direction, selection, selection_direction
diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py
index 52001cbd87..1776916779 100644
--- a/flopy/plot/plotutil.py
+++ b/flopy/plot/plotutil.py
@@ -10,6 +10,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 
 from ..datbase import DataInterface, DataType
 from ..utils import Util3d, import_optional_dependency
@@ -2591,3 +2592,347 @@ def parse_modpath_selection_options(
         tep = ep.copy()
 
     return tep, istart, xp, yp
+
+
+# define minimum expected fields
+_mp7_pathline_fields = ["x", "y", "z", "time", "k", "particleid"]
+_prt_pathline_fields = [
+    "x",
+    "y",
+    "z",
+    "t",
+    "trelease",
+    "imdl",
+    "iprp",
+    "irpt",
+    "ilay",
+]
+
+
+def to_mp7_pathlines(
+    data: Union[np.recarray, pd.DataFrame]
+) -> Union[np.recarray, pd.DataFrame]:
+    """
+    Convert MODFLOW 6 PRT pathline data to MODPATH 7 pathline format.
+
+    Parameters
+    ----------
+    data : np.recarray or pd.DataFrame
+        MODFLOW 6 PRT pathline data
+
+    Returns
+    -------
+    np.recarray or pd.DataFrame (consistent with input type)
+    """
+
+    # determine return type
+    ret_type = type(data)
+
+    # convert to dataframe if needed
+    if not isinstance(data, pd.DataFrame):
+        data = pd.DataFrame(data)
+
+    # check format
+    dt = data.dtypes
+    if not (
+        all(n in dt for n in _mp7_pathline_fields)
+        or all(n in dt for n in _prt_pathline_fields)
+    ):
+        raise ValueError(
+            "Pathline data must contain all of the following columns: "
+            f"{_mp7_pathline_fields} for MODPATH 7, or "
+            f"{_prt_pathline_fields} for MODFLOW 6 PRT"
+        )
+
+    # return early if already in MP7 format
+    if "t" not in dt:
+        return (
+            data if ret_type == pd.DataFrame else data.to_records(index=False)
+        )
+
+    # assign a unique particle index column incrementing an integer
+    # for each unique combination of irpt, iprp, imdl, and trelease
+    data = data.sort_values(["imdl", "iprp", "irpt", "trelease"])
+    particles = data.groupby(["imdl", "iprp", "irpt", "trelease"])
+    seqn_key = "sequencenumber"
+    data[seqn_key] = particles.ngroup()
+
+    # convert to recarray
+    data = data.to_records(index=False)
+
+    # define mp7 dtype
+    mp7_dtypes = np.dtype(
+        [
+            ("particleid", np.int32),  # same as sequencenumber
+            ("particlegroup", np.int32),
+            (
+                seqn_key,
+                np.int32,
+            ),  # mp7 sequencenumber (globally unique auto-generated ID)
+            (
+                "particleidloc",
+                np.int32,
+            ),  # mp7 particle ID (unique within a group, user-assigned or autogenerated)
+            ("time", np.float32),
+            ("x", np.float32),
+            ("y", np.float32),
+            ("z", np.float32),
+            ("k", np.int32),
+            ("node", np.int32),
+            ("xloc", np.float32),
+            ("yloc", np.float32),
+            ("zloc", np.float32),
+            ("stressperiod", np.int32),
+            ("timestep", np.int32),
+        ]
+    )
+
+    # build mp7 format recarray
+    ret = np.core.records.fromarrays(
+        [
+            data[seqn_key],
+            data["iprp"],
+            data[seqn_key],
+            data["irpt"],
+            data["t"],
+            data["x"],
+            data["y"],
+            data["z"],
+            data["ilay"],
+            data["icell"],
+            # todo local coords (xloc, yloc, zloc)
+            np.zeros(data.shape[0]),
+            np.zeros(data.shape[0]),
+            np.zeros(data.shape[0]),
+            data["kper"],
+            data["kstp"],
+        ],
+        dtype=mp7_dtypes,
+    )
+
+    return pd.DataFrame(ret) if ret_type == pd.DataFrame else ret
+
+
+def to_mp7_endpoints(
+    data: Union[np.recarray, pd.DataFrame]
+) -> Union[np.recarray, pd.DataFrame]:
+    """
+    Convert MODFLOW 6 PRT pathline data to MODPATH 7 endpoint format.
+
+    Parameters
+    ----------
+    data : np.recarray or pd.DataFrame
+        MODFLOW 6 PRT pathline data
+
+    Returns
+    -------
+    np.recarray or pd.DataFrame (consistent with input type)
+    """
+
+    # determine return type
+    ret_type = type(data)
+
+    # convert to dataframe if needed
+    if isinstance(data, np.recarray):
+        data = pd.DataFrame(data)
+
+    # assign a unique particle index column incrementing an integer
+    # for each unique combination of irpt, iprp, imdl, and trelease
+    data = data.sort_values(["imdl", "iprp", "irpt", "trelease"])
+    particles = data.groupby(["imdl", "iprp", "irpt", "trelease"])
+    seqn_key = "sequencenumber"
+    data[seqn_key] = particles.ngroup()
+
+    # select startpoints and endpoints, sorting by sequencenumber
+    startpts = (
+        data.sort_values("t").groupby(seqn_key).head(1).sort_values(seqn_key)
+    )
+    endpts = (
+        data.sort_values("t").groupby(seqn_key).tail(1).sort_values(seqn_key)
+    )
+
+    # add columns for
+    pairings = [
+        # initial coordinates
+        ("x0", "x"),
+        ("y0", "y"),
+        ("z0", "z"),
+        # initial zone
+        ("zone0", "izone"),
+        # initial node number
+        ("node0", "icell"),
+        # initial layer
+        ("k0", "ilay"),
+    ]
+    conditions = [
+        startpts[seqn_key].eq(row[seqn_key]) for _, row in startpts.iterrows()
+    ]
+    for fl, fr in pairings:
+        endpts[fl] = np.select(conditions, startpts[fr].to_numpy())
+
+    # convert to recarray
+    endpts = endpts.to_records(index=False)
+
+    # define mp7 dtype
+    mp7_dtype = np.dtype(
+        [
+            (
+                "particleid",
+                np.int32,
+            ),  # mp7 sequencenumber (globally unique auto-generated ID)
+            ("particlegroup", np.int32),
+            (
+                "particleidloc",
+                np.int32,
+            ),  # mp7 particle ID (unique within a group, user-assigned or autogenerated)
+            ("status", np.int32),
+            ("time0", np.float32),
+            ("time", np.float32),
+            ("node0", np.int32),
+            ("k0", np.int32),
+            ("xloc0", np.float32),
+            ("yloc0", np.float32),
+            ("zloc0", np.float32),
+            ("x0", np.float32),
+            ("y0", np.float32),
+            ("z0", np.float32),
+            ("zone0", np.int32),
+            ("initialcellface", np.int32),
+            ("node", np.int32),
+            ("k", np.int32),
+            ("xloc", np.float32),
+            ("yloc", np.float32),
+            ("zloc", np.float32),
+            ("x", np.float32),
+            ("y", np.float32),
+            ("z", np.float32),
+            ("zone", np.int32),
+            ("cellface", np.int32),
+        ]
+    )
+
+    # build mp7 format recarray
+    ret = np.core.records.fromarrays(
+        [
+            endpts["sequencenumber"],
+            endpts["iprp"],
+            endpts["irpt"],
+            endpts["istatus"],
+            endpts["trelease"],
+            endpts["t"],
+            endpts["node0"],
+            endpts["k0"],
+            # todo initial local coords (xloc0, yloc0, zloc0)
+            np.zeros(endpts.shape[0]),
+            np.zeros(endpts.shape[0]),
+            np.zeros(endpts.shape[0]),
+            endpts["x0"],
+            endpts["y0"],
+            endpts["z0"],
+            endpts["zone0"],
+            np.zeros(endpts.shape[0]),  # todo initial cell face?
+            endpts["icell"],
+            endpts["ilay"],
+            # todo local coords (xloc, yloc, zloc)
+            np.zeros(endpts.shape[0]),
+            np.zeros(endpts.shape[0]),
+            np.zeros(endpts.shape[0]),
+            endpts["x"],
+            endpts["y"],
+            endpts["z"],
+            endpts["izone"],
+            np.zeros(endpts.shape[0]),  # todo cell face?
+        ],
+        dtype=mp7_dtype,
+    )
+
+    return pd.DataFrame(ret) if ret_type == pd.DataFrame else ret
+
+
+def to_prt_pathlines(
+    data: Union[np.recarray, pd.DataFrame]
+) -> Union[np.recarray, pd.DataFrame]:
+    """
+    Convert MODPATH 7 pathline or endpoint data to MODFLOW 6 PRT pathline format.
+
+    Parameters
+    ----------
+    data : np.recarray or pd.DataFrame
+        MODPATH 7 pathline or endpoint data
+
+    Returns
+    -------
+    np.recarray or pd.DataFrame (consistent with input type)
+    """
+
+    # determine return type
+    ret_type = type(data)
+
+    # convert to dataframe if needed
+    if isinstance(data, np.recarray):
+        data = pd.DataFrame(data)
+
+    # check format
+    dt = data.dtypes
+    if not (
+        all(n in dt for n in _mp7_pathline_fields)
+        or all(n in dt for n in _prt_pathline_fields)
+    ):
+        raise ValueError(
+            "Pathline data must contain all of the following columns: "
+            f"{_mp7_pathline_fields} for MODPATH 7, or "
+            f"{_prt_pathline_fields} for MODFLOW 6 PRT"
+        )
+
+    # return early if already in PRT format
+    if "t" in dt:
+        return (
+            data if ret_type == pd.DataFrame else data.to_records(index=False)
+        )
+
+    # convert to recarray
+    data = data.to_records(index=False)
+
+    # define prt dtype
+    prt_dtypes = np.dtype(
+        [
+            ("kper", np.int32),
+            ("kstp", np.int32),
+            ("imdl", np.int32),
+            ("iprp", np.int32),
+            ("irpt", np.int32),
+            ("ilay", np.int32),
+            ("icell", np.int32),
+            ("izone", np.int32),
+            ("istatus", np.int32),
+            ("ireason", np.int32),
+            ("trelease", np.float32),
+            ("t", np.float32),
+            ("x", np.float32),
+            ("y", np.float32),
+            ("z", np.float32),
+        ]
+    )
+
+    # build prt format recarray
+    ret = np.core.records.fromarrays(
+        [
+            data["stressperiod"],
+            data["timestep"],
+            np.zeros(data.shape[0]),
+            np.zeros(data.shape[0]),
+            data["particlegroup"],
+            data["particleid"],
+            data["k"],
+            data["node"],
+            np.zeros(data.shape[0]),  # todo k
+            np.zeros(data.shape[0]),  # todo trelease?
+            data["t"],
+            data["x"],
+            data["y"],
+            data["z"],
+        ],
+        dtype=prt_dtypes,
+    )
+
+    return pd.DataFrame(ret) if ret_type == pd.DataFrame else ret
diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py
index cf384fa359..875cda8944 100644
--- a/flopy/utils/binaryfile.py
+++ b/flopy/utils/binaryfile.py
@@ -11,7 +11,7 @@
 import os
 import warnings
 from pathlib import Path
-from typing import Optional, Union
+from typing import List, Optional, Union
 
 import numpy as np
 
@@ -1580,7 +1580,7 @@ def get_data(
         text=None,
         paknam=None,
         full3D=False,
-    ):
+    ) -> Union[List, np.ndarray]:
         """
         Get data from the binary budget file.