diff --git a/autotest/prt_test_utils.py b/autotest/prt/prt_test_utils.py
similarity index 61%
rename from autotest/prt_test_utils.py
rename to autotest/prt/prt_test_utils.py
index 92506baa570..cc7d2957d05 100644
--- a/autotest/prt_test_utils.py
+++ b/autotest/prt/prt_test_utils.py
@@ -1,5 +1,6 @@
 import os
-from typing import Union
+from types import SimpleNamespace
+from typing import Tuple, Union
 
 import flopy
 import numpy as np
@@ -39,7 +40,7 @@ def check_track_data(
     # check each column separately to avoid:
     # TypeError: The DType <class 'numpy._FloatAbstractDType'> could not be promoted by <class 'numpy.dtype[void]'>
     for k in data_bin.dtype.names:
-        if k == 'name':
+        if k == "name":
             continue
         assert np.allclose(data_bin[k], data_csv[k], equal_nan=True)
 
@@ -93,18 +94,23 @@ def get_output_event(case_name):
         if "wksk" in case_name
         else "TERMINATE"
         if "terminate" in case_name
-        else "ALL" # default
+        else "ALL"  # default
     )
 
 
 def get_ireason_code(output_event):
     return (
-        0 if output_event == "RELEASE"
-        else 1 if output_event == "TRANSIT"
-        else 2 if output_event == "TIMESTEP"
-        else 3 if output_event == "TERMINATE"
-        else 4 if output_event == "WEAKSINK"
-        else -1 # default
+        0
+        if output_event == "RELEASE"
+        else 1
+        if output_event == "TRANSIT"
+        else 2
+        if output_event == "TIMESTEP"
+        else 3
+        if output_event == "TERMINATE"
+        else 4
+        if output_event == "WEAKSINK"
+        else -1  # default
     )
 
 
@@ -140,3 +146,89 @@ def get_partdata(grid, rpts):
             timeoffset=0,
             drape=0,
         )
+
+
+def build_gwf_sim(
+    name, ws, mf6
+) -> Tuple[flopy.mf6.MFSimulation, SimpleNamespace]:
+    ctx = SimpleNamespace(
+        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,
+    )
+
+    # create simulation
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        exe_name=mf6,
+        version="mf6",
+        sim_ws=ws,
+    )
+
+    # create tdis package
+    flopy.mf6.modflow.mftdis.ModflowTdis(
+        sim,
+        pname="tdis",
+        time_units="DAYS",
+        nper=ctx.nper,
+        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
+    )
+
+    # create gwf model
+    gwfname = f"{name}_gwf"
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True)
+
+    # create gwf discretization
+    flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
+        gwf,
+        pname="dis",
+        nlay=ctx.nlay,
+        nrow=ctx.nrow,
+        ncol=ctx.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
+    # output file names
+    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, ctx
diff --git a/autotest/test_prt_exg01.py b/autotest/prt/test_prt_exg01.py
similarity index 90%
rename from autotest/test_prt_exg01.py
rename to autotest/prt/test_prt_exg01.py
index ad601a36e34..db6ecd1bae4 100644
--- a/autotest/test_prt_exg01.py
+++ b/autotest/prt/test_prt_exg01.py
@@ -2,8 +2,8 @@
 Test GWF and PRT models in the same simulation
 with an exchange.
 
-The grid is a 10x10 square with a single layer.
-The same flow system shown on the FloPy readme.
+The grid is a 10x10 square with a single layer,
+the same flow system shown on the FloPy readme.
 Particles are released from the top left cell.
 
 Results are compared against a MODPATH 7 model.
@@ -21,39 +21,21 @@
 import numpy as np
 import pandas as pd
 import pytest
+from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from flopy.plot.plotutil import to_mp7_pathlines
+from prt_test_utils import (all_equal, build_gwf_sim, check_budget_data,
+                            check_track_data, has_default_boundnames)
 
 from framework import TestFramework
-from prt_test_utils import (
-    all_equal,
-    check_budget_data,
-    check_track_data,
-    has_default_boundnames,
-)
 from simulation import TestSimulation
 
-
 # simulation/model names
-simname = "prtexg1"
+simname = "prtexg01"
 
 # test cases
 ex = [simname, f"{simname}bnms"]
 
-
-# model info
-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
-
 # release points
 releasepts = [
     # index, k, i, j, x, y, z
@@ -75,11 +57,9 @@ def get_model_name(idx, mdl):
 
 
 def build_sim(idx, ws, mf6):
-    from test_prt_fmi01 import build_gwf_sim
-
     # create simulation
     name = ex[idx]
-    sim = build_gwf_sim(name, ws, mf6)
+    sim, ctx = build_gwf_sim(name, ws, mf6)
 
     # create prt model
     prtname = get_model_name(idx, "prt")
@@ -89,13 +69,13 @@ def build_sim(idx, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=nlay,
-        nrow=nrow,
-        ncol=ncol,
+        nlay=ctx.nlay,
+        nrow=ctx.nrow,
+        ncol=ctx.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
 
     # create prp package
     rpts = (
@@ -124,7 +104,7 @@ def build_sim(idx, ws, mf6):
     )
 
     # create a flow model interface
-    # todo Fienen's report (crash when FMI created but not needed)
+    # todo Mike Fienen's report (crash when FMI created but not needed)
     # flopy.mf6.ModflowPrtfmi(
     #     prt,
     #     packagedata=[
@@ -151,10 +131,10 @@ def build_sim(idx, ws, mf6):
     )
     sim.register_solution_package(ems, [prt.name])
 
-    return sim
+    return sim, ctx
 
 
-def build_mp7_sim(idx, ws, mp7, gwf):
+def build_mp7_sim(ctx, idx, ws, mp7, gwf):
     partdata = flopy.modpath.ParticleData(
         partlocs=[p[0] for p in releasepts_mp7],
         localx=[p[1] for p in releasepts_mp7],
@@ -177,7 +157,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=porosity,
+        porosity=ctx.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -191,12 +171,12 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     return mp
 
 
-def eval_results(sim):
+def eval_results(ctx, sim):
     print(f"Evaluating results for sim {sim.name}")
     simpath = Path(sim.simpath)
 
     # check budget data
-    check_budget_data(simpath / f"{sim.name}_prt.lst", perlen, nper)
+    check_budget_data(simpath / f"{sim.name}_prt.lst", ctx.perlen, ctx.nper)
 
     # check particle track data
     prt_track_file = simpath / f"{sim.name}_prt.trk"
@@ -215,7 +195,7 @@ def eval_results(sim):
 @pytest.mark.parametrize("idx, name", enumerate(ex))
 def test_mf6model(idx, name, function_tmpdir, targets):
     ws = function_tmpdir
-    sim = build_sim(idx, str(ws), targets.mf6)
+    sim, ctx = build_sim(idx, str(ws), targets.mf6)
     sim.write_simulation()
 
     test = TestFramework()
@@ -223,7 +203,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
         TestSimulation(
             name=name,
             exe_dict=targets,
-            exfunc=eval_results,
+            exfunc=lambda s: eval_results(ctx, s),
             idxsim=0,
             make_comparison=False,
         ),
@@ -243,7 +223,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(idx, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(ctx, idx, ws, targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -293,7 +273,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
         assert pd.isna(mf6_pls["name"]).all()
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{name}_prt.lst", perlen, nper)
+    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
diff --git a/autotest/test_prt_fmi01.py b/autotest/prt/test_prt_fmi01.py
similarity index 79%
rename from autotest/test_prt_fmi01.py
rename to autotest/prt/test_prt_fmi01.py
index 8107247548a..9b28b8b5595 100644
--- a/autotest/test_prt_fmi01.py
+++ b/autotest/prt/test_prt_fmi01.py
@@ -32,17 +32,12 @@
 import numpy as np
 import pandas as pd
 import pytest
+from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from flopy.plot.plotutil import to_mp7_pathlines
-from prt_test_utils import (
-    all_equal,
-    check_budget_data,
-    check_track_data,
-    get_partdata,
-    has_default_boundnames,
-)
-
+from prt_test_utils import (all_equal, build_gwf_sim, check_budget_data,
+                            check_track_data, get_partdata,
+                            has_default_boundnames)
 
 # simulation name
 simname = "prtfmi01"
@@ -56,18 +51,6 @@ def get_model_name(idx, mdl):
     return f"{ex[idx]}_{mdl}"
 
 
-# model info
-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
-
 # release points in mp7 format (using local coordinates)
 releasepts_mp7 = [
     # node number, localx, localy, localz
@@ -86,78 +69,7 @@ def get_model_name(idx, mdl):
 ]
 
 
-def build_gwf_sim(name, ws, mf6):
-    # create simulation
-    sim = flopy.mf6.MFSimulation(
-        sim_name=name,
-        exe_name=mf6,
-        version="mf6",
-        sim_ws=ws,
-    )
-
-    # create tdis package
-    flopy.mf6.modflow.mftdis.ModflowTdis(
-        sim,
-        pname="tdis",
-        time_units="DAYS",
-        nper=nper,
-        perioddata=[(perlen, nstp, tsmult)],
-    )
-
-    # create gwf model
-    gwfname = f"{name}_gwf"
-    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
-    # output file names
-    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
-
-
-def build_prt_sim(idx, ws, mf6):
+def build_prt_sim(ctx, idx, ws, mf6):
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=ex[idx],
@@ -171,8 +83,8 @@ def build_prt_sim(idx, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=nper,
-        perioddata=[(perlen, nstp, tsmult)],
+        nper=ctx.nper,
+        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
     )
 
     # create prt model
@@ -183,13 +95,13 @@ def build_prt_sim(idx, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=nlay,
-        nrow=nrow,
-        ncol=ncol,
+        nlay=ctx.nlay,
+        nrow=ctx.nrow,
+        ncol=ctx.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
 
     # convert mp7 particledata to prt release points
     partdata = get_partdata(prt.modelgrid, releasepts_mp7)
@@ -248,7 +160,7 @@ def build_prt_sim(idx, ws, mf6):
     return sim
 
 
-def build_mp7_sim(idx, ws, mp7, gwf):
+def build_mp7_sim(ctx, idx, ws, mp7, gwf):
     # convert mp7 particledata to prt release points
     partdata = get_partdata(gwf.modelgrid, releasepts_mp7)
 
@@ -267,7 +179,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=porosity,
+        porosity=ctx.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -282,7 +194,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
 
 
 @pytest.mark.parametrize("idx, name", list(enumerate(ex)))
-def test_prt_fmi01(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets):
     # workspace
     ws = function_tmpdir
 
@@ -295,8 +207,8 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets):
     mp7name = get_model_name(idx, "mp7")
 
     # build mf6 models
-    gwfsim = build_gwf_sim(name, ws, targets.mf6)
-    prtsim = build_prt_sim(idx, ws, targets.mf6)
+    gwfsim, ctx = build_gwf_sim(name, ws, targets.mf6)
+    prtsim = build_prt_sim(ctx, idx, ws, targets.mf6)
 
     # run mf6 models
     for sim in [gwfsim, prtsim]:
@@ -312,7 +224,7 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(idx, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(ctx, idx, ws, targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -359,7 +271,7 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets):
     assert all_equal(mf6_pls["iprp"], 1)
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{name}_prt.lst", perlen, nper)
+    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
 
     # check mf6 prt particle track data were written to binary/CSV files
     # and that different formats are equal
diff --git a/autotest/test_prt_fmi02.py b/autotest/prt/test_prt_fmi02.py
similarity index 76%
rename from autotest/test_prt_fmi02.py
rename to autotest/prt/test_prt_fmi02.py
index 260c48a7904..f24d31b2b33 100644
--- a/autotest/test_prt_fmi02.py
+++ b/autotest/prt/test_prt_fmi02.py
@@ -32,22 +32,14 @@
 import numpy as np
 import pandas as pd
 import pytest
+from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from flopy.plot.plotutil import to_mp7_pathlines
-
-from prt_test_utils import (
-    check_budget_data,
-    check_track_data,
-    get_output_event,
-)
-
+from prt_test_utils import (build_gwf_sim, check_budget_data, check_track_data,
+                            get_output_event)
 
 # simulation/model names
 simname = "prtfmi02"
-gwfname = f"{simname}_gwf"
-prtname = f"{simname}_prt"
-mp7name = f"{simname}_mp7"
 
 # test cases
 ex = [
@@ -58,24 +50,14 @@
     f"{simname}wksk",
 ]
 
-# output file names
-gwf_budget_file = f"{gwfname}.bud"
-gwf_head_file = f"{gwfname}.hds"
-prt_track_file = f"{prtname}.trk"
-prt_track_csv_file = f"{prtname}.trk.csv"
-mp7_pathline_file = f"{mp7name}.mppth"
-
-# model info
-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
+
+# function to create idomain from grid dimensions
+def create_idomain(nlay, nrow, ncol):
+    idmn = np.ones((nlay, nrow, ncol), dtype=int)
+    idmn[0, 0, 9] = 0
+    idmn[0, 9, 0] = 0
+    return idmn
+
 
 # release points
 # todo: define for mp7 first, then use flopy utils to convert to global coords for mf6 prt
@@ -104,84 +86,16 @@
     for i in range(5)
 ]
 
-# idomain
-idomain = np.ones((nlay, nrow, ncol), dtype=int)
-idomain[0, 0, 9] = 0
-idomain[0, 9, 0] = 0
-# idomain = idomain.ravel()
 
+def build_prt_sim(ctx, name, ws, mf6):
+    gwfname = f"{name}_gwf"
+    prtname = f"{name}_prt"
+    gwf_budget_file = f"{gwfname}.bud"
+    gwf_head_file = f"{gwfname}.hds"
+    prt_track_file = f"{prtname}.trk"
+    prt_track_csv_file = f"{prtname}.trk.csv"
 
-def build_gwf_sim(idx, ws, mf6):
     # create simulation
-    sim = flopy.mf6.MFSimulation(
-        sim_name=ex[idx],
-        exe_name=mf6,
-        version="mf6",
-        sim_ws=ws,
-    )
-
-    # 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,
-        idomain=idomain,
-    )
-
-    # 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
-    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
-
-
-def build_prt_sim(idx, ws, mf6):
-    # create simulation
-    name = ex[idx]
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
@@ -194,25 +108,27 @@ def build_prt_sim(idx, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=nper,
-        perioddata=[(perlen, nstp, tsmult)],
+        nper=ctx.nper,
+        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
     )
 
     # create prt model
     prt = flopy.mf6.ModflowPrt(sim, modelname=prtname)
 
     # create prt discretization
+    # idomain
+
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=nlay,
-        nrow=nrow,
-        ncol=ncol,
-        idomain=idomain,
+        nlay=ctx.nlay,
+        nrow=ctx.nrow,
+        ncol=ctx.ncol,
+        idomain=create_idomain(ctx.nlay, ctx.nrow, ctx.ncol),
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
 
     # create prp packages
     flopy.mf6.ModflowPrtprp(
@@ -259,10 +175,13 @@ def build_prt_sim(idx, ws, mf6):
     )
     sim.register_solution_package(ems, [prt.name])
 
-    return sim
+    return sim, ctx
 
 
-def build_mp7_sim(ws, mp7, gwf):
+def build_mp7_sim(ctx, name, ws, mp7, gwf):
+    mp7name = f"{name}_mp7"
+    mp7_pathline_file = f"{mp7name}.mppth"
+
     pd_a = flopy.modpath.ParticleData(
         partlocs=[p[0] for p in releasepts_mp7_a],
         localx=[p[1] for p in releasepts_mp7_a],
@@ -297,7 +216,7 @@ def build_mp7_sim(ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=porosity,
+        porosity=ctx.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -311,26 +230,31 @@ def build_mp7_sim(ws, mp7, gwf):
     return mp
 
 
-@pytest.mark.parametrize("idx, name", enumerate(ex))
-def test_prt_fmi02(idx, name, function_tmpdir, targets):
+@pytest.mark.parametrize("name", ex)
+def test_mf6model(name, function_tmpdir, targets):
     ws = function_tmpdir
+    gwfname = f"{name}_gwf"
+    prtname = f"{name}_prt"
+    mp7name = f"{name}_mp7"
+    gwf_budget_file = f"{gwfname}.bud"
+    gwf_head_file = f"{gwfname}.hds"
+    prt_track_file = f"{prtname}.trk"
+    prt_track_csv_file = f"{prtname}.trk.csv"
+    mp7_pathline_file = f"{mp7name}.mppth"
 
     # build mf6 simulations
-    gwfsim = build_gwf_sim(idx, ws, targets.mf6)
-    prtsim = build_prt_sim(idx, ws, targets.mf6)
+    gwfsim, ctx = build_gwf_sim(name, ws, targets.mf6)
+    gwf = gwfsim.get_model()
+    dis = gwf.get_package("DIS")
+    dis.idomain = create_idomain(ctx.nlay, ctx.nrow, ctx.ncol)
 
-    # write mf6 simulation input files
-    gwfsim.write_simulation()
-    # run mf6 gwf simulation
-    success, _ = gwfsim.run_simulation()
-    assert success
+    prtsim, _ = build_prt_sim(ctx, name, ws, targets.mf6)
 
-    # write mf6 prt simulation input files
-    prtsim.write_simulation()
-
-    # run mf6 prt simulation
-    success, _ = prtsim.run_simulation()
-    assert success
+    # run mf6 models
+    for sim in [gwfsim, prtsim]:
+        sim.write_simulation()
+        success, _ = sim.run_simulation()
+        assert success
 
     # extract models
     gwf = gwfsim.get_model(gwfname)
@@ -340,7 +264,7 @@ def test_prt_fmi02(idx, name, function_tmpdir, targets):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -375,7 +299,7 @@ def test_prt_fmi02(idx, name, function_tmpdir, targets):
     # if event is TRANSIT, expect full results minus start loc.
     # if event is TIMESTEP or WEAKSINK, output should be empty.
     # in either case, return early and skip MP7 comparison.
-    event = get_output_event(ex[idx])
+    event = get_output_event(name)
     if event == "RELEASE" or event == "TERMINATE":
         assert len(mf6_pls) == len(releasepts_a) + len(releasepts_b)
         return
@@ -404,7 +328,7 @@ def all_equal(col, val):
     )
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{simname}_prt.lst", perlen, nper)
+    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
diff --git a/autotest/test_prt_fmi03.py b/autotest/prt/test_prt_fmi03.py
similarity index 75%
rename from autotest/test_prt_fmi03.py
rename to autotest/prt/test_prt_fmi03.py
index ec8c297f63e..188ae92fa25 100644
--- a/autotest/test_prt_fmi03.py
+++ b/autotest/prt/test_prt_fmi03.py
@@ -23,6 +23,7 @@
 """
 
 
+from itertools import repeat
 from pathlib import Path
 
 import flopy
@@ -31,35 +32,34 @@
 import numpy as np
 import pandas as pd
 import pytest
+from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from flopy.plot.plotutil import to_mp7_pathlines
 from matplotlib.collections import LineCollection
-from prt_test_utils import check_budget_data, check_track_data
-
+from prt_test_utils import build_gwf_sim, check_budget_data, check_track_data
 
 # simulation name
 simname = "prtfmi03"
 
 # test cases
-ex = [f"{simname}_1l", f"{simname}_2l"]
+ex = [f"{simname}_l1", f"{simname}_l2"]
 
 
 # model names
-def get_model_name(idx, mdl):
-    return f"{ex[idx]}_{mdl}"
+def get_model_name(name, mdl):
+    return f"{name}_{mdl}"
 
 
-# model info
-nlay = 1
-nrow = 10
-ncol = 10
-top = 1.0
-nper = 1
-perlen = 1.0
-nstp = 1
-tsmult = 1.0
-porosity = 0.1
+# izone
+stopzone_cells = [(0, 1, 8), (0, 8, 1)]
+
+
+def create_izone(nlay, nrow, ncol):
+    izone = np.zeros((nlay, nrow, ncol), dtype=int)
+    for iz in stopzone_cells:
+        izone[iz] = 1
+    return izone
+
 
 # release points
 # todo: define for mp7 first, then use flopy utils to convert to global coords for mf6 prt
@@ -76,90 +76,11 @@ def get_model_name(idx, mdl):
     for i in range(9)
 ]
 
-# zone array and stop zone cell ids
-stopzone_cells = [(0, 1, 8), (0, 8, 1)]
-izone = np.zeros((nlay, nrow, ncol), dtype=int)
-for iz in stopzone_cells:
-    izone[iz] = 1
-
-
-def build_gwf_sim(idx, ws, mf6):
-    # create simulation
-    sim = flopy.mf6.MFSimulation(
-        sim_name=ex[idx],
-        exe_name=mf6,
-        version="mf6",
-        sim_ws=ws,
-    )
-
-    # create tdis package
-    flopy.mf6.modflow.mftdis.ModflowTdis(
-        sim,
-        pname="tdis",
-        time_units="DAYS",
-        nper=nper,
-        perioddata=[(perlen, nstp, tsmult)],
-    )
-
-    # create gwf model
-    gwfname = get_model_name(idx, "gwf")
-    gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True)
-
-    # create gwf discretization
-    botm = [top - (k + 1) for k in range(nlay + idx)]
-    flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
-        gwf,
-        pname="dis",
-        nlay=nlay + idx,
-        nrow=nrow,
-        ncol=ncol,
-        top=top,
-        botm=botm,
-    )
-
-    # 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
-
 
-def build_prt_sim(idx, ws, mf6):
+def build_prt_sim(ctx, name, ws, mf6):
     # create simulation
     sim = flopy.mf6.MFSimulation(
-        sim_name=ex[idx],
+        sim_name=name,
         exe_name=mf6,
         version="mf6",
         sim_ws=ws,
@@ -170,32 +91,34 @@ def build_prt_sim(idx, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=nper,
-        perioddata=[(perlen, nstp, tsmult)],
+        nper=ctx.nper,
+        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
     )
 
     # create prt model
-    prtname = get_model_name(idx, "prt")
+    prtname = get_model_name(name, "prt")
     prt = flopy.mf6.ModflowPrt(sim, modelname=prtname)
 
     # create prt discretization
-    botm = [top - (k + 1) for k in range(nlay + idx)]
+    ctx.nlay = int(name[-1])
+    botm = [ctx.top - (k + 1) for k in range(ctx.nlay)]
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=nlay + idx,
-        nrow=nrow,
-        ncol=ncol,
-        top=top,
+        nlay=ctx.nlay,
+        nrow=ctx.nrow,
+        ncol=ctx.ncol,
+        top=ctx.top,
         botm=botm,
     )
 
     # create mip package
+    izone = create_izone(ctx.nlay, ctx.nrow, ctx.ncol)
     flopy.mf6.ModflowPrtmip(
         prt,
         pname="mip",
-        porosity=porosity,
-        izone=izone if idx == 0 else np.array([izone, izone]),
+        porosity=ctx.porosity,
+        izone=izone,  # if ctx.nlay == 1 else np.array([izone, izone]),
     )
 
     # create prp package
@@ -220,7 +143,7 @@ def build_prt_sim(idx, ws, mf6):
     )
 
     # create the flow model interface
-    gwfname = get_model_name(idx, "gwf")
+    gwfname = get_model_name(name, "gwf")
     gwf_budget_file = f"{gwfname}.bud"
     gwf_head_file = f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
@@ -239,10 +162,10 @@ def build_prt_sim(idx, ws, mf6):
     )
     sim.register_solution_package(ems, [prt.name])
 
-    return sim
+    return sim, ctx
 
 
-def build_mp7_sim(idx, ws, mp7, gwf):
+def build_mp7_sim(ctx, name, ws, mp7, gwf):
     partdata = flopy.modpath.ParticleData(
         partlocs=[p[0] for p in releasepts_mp7],
         localx=[p[1] for p in releasepts_mp7],
@@ -251,7 +174,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
         timeoffset=0,
         drape=0,
     )
-    mp7name = get_model_name(idx, "mp7")
+    mp7name = get_model_name(name, "mp7")
     pg = flopy.modpath.ParticleGroup(
         particlegroupname="G1",
         particledata=partdata,
@@ -265,8 +188,10 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=porosity,
+        porosity=ctx.porosity,
     )
+    ctx.nlay = int(name[-1])
+    izone = create_izone(ctx.nlay, ctx.nrow, ctx.ncol)
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
         simulationtype="pathline",
@@ -274,7 +199,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
         budgetoutputoption="summary",
         stoptimeoption="extend",
         stopzone=1,
-        zones=izone if idx == 0 else np.array([izone, izone]),
+        zones=izone,  # if ctx.nlay == 1 else np.array([izone, izone]),
         zonedataoption="on",
         particlegroups=[pg],
     )
@@ -282,18 +207,27 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     return mp
 
 
-@pytest.mark.parametrize("idx, name", list(enumerate(ex)))
-def test_mf6model(idx, name, function_tmpdir, targets):
+@pytest.mark.parametrize("name", ex)
+def test_mf6model(name, function_tmpdir, targets):
     ws = function_tmpdir
 
     # define model names
-    gwfname = get_model_name(idx, "gwf")
-    prtname = get_model_name(idx, "prt")
-    mp7name = get_model_name(idx, "mp7")
+    gwfname = get_model_name(name, "gwf")
+    prtname = get_model_name(name, "prt")
+    mp7name = get_model_name(name, "mp7")
 
     # build mf6 simulations
-    gwfsim = build_gwf_sim(idx, ws, targets.mf6)
-    prtsim = build_prt_sim(idx, ws, targets.mf6)
+    gwfsim, ctx = build_gwf_sim(name, ws, targets.mf6)
+    gwf = gwfsim.get_model()
+    dis = gwf.get_package("DIS")
+    ctx.nlay = int(name[-1])
+    botm = [ctx.top - (k + 1) for k in range(ctx.nlay)]
+    botm_data = np.array(
+        [list(repeat(b, ctx.nrow * ctx.ncol)) for b in botm]
+    ).reshape((ctx.nlay, ctx.nrow, ctx.ncol))
+    dis.nlay = ctx.nlay
+    dis.botm.set_data(botm_data)
+    prtsim, ctx = build_prt_sim(ctx, name, ws, targets.mf6)
 
     # run mf6 simulations
     for sim in [gwfsim, prtsim]:
@@ -309,7 +243,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(idx, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -344,7 +278,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     mf6_pls = pd.read_csv(ws / prt_track_csv_file)
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{name}_prt.lst", perlen, nper)
+    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
diff --git a/autotest/test_prt_fmi04.py b/autotest/prt/test_prt_fmi04.py
similarity index 77%
rename from autotest/test_prt_fmi04.py
rename to autotest/prt/test_prt_fmi04.py
index c0fd22701d6..d827bbc5719 100644
--- a/autotest/test_prt_fmi04.py
+++ b/autotest/prt/test_prt_fmi04.py
@@ -31,45 +31,18 @@
 import numpy as np
 import pandas as pd
 import pytest
+from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from flopy.plot.plotutil import to_mp7_pathlines
-
-from prt_test_utils import (
-    check_budget_data,
-    check_track_data,
-    get_ireason_code,
-)
-
+from prt_test_utils import (build_gwf_sim, check_budget_data, check_track_data,
+                            get_ireason_code)
 
 # simulation/model names
 simname = "prtfmi04"
-gwfname = f"{simname}_gwf"
-prtname = f"{simname}_prt"
-mp7name = f"{simname}_mp7"
 
 # test cases
 ex = [simname, f"{simname}saws"]
 
-# output file names
-gwf_budget_file = f"{gwfname}.bud"
-gwf_head_file = f"{gwfname}.hds"
-prt_track_file = f"{prtname}.trk"
-prt_track_csv_file = f"{prtname}.trk.csv"
-mp7_pathline_file = f"{mp7name}.mppth"
-
-# model info
-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
-
 # release points
 # todo: define for mp7 first, then use flopy utils to convert to global coords for mf6 prt
 releasepts = [
@@ -86,103 +59,30 @@
 ]
 
 
-def build_gwf_sim(idx, ws, mf6):
-    # create simulation
-    sim = flopy.mf6.MFSimulation(
-        sim_name=ex[idx],
-        exe_name=mf6,
-        version="mf6",
-        sim_ws=ws,
-    )
-
-    # create tdis package
-    pd = (perlen, nstp, tsmult)
-    flopy.mf6.modflow.mftdis.ModflowTdis(
-        sim,
-        pname="tdis",
-        time_units="DAYS",
-        nper=nper,
-        perioddata=[pd],
-    )
-
-    # 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"],
-    )
+def build_prt_sim(ctx, name, ws, mf6):
+    # output file names
+    gwfname = f"{name}_gwf"
+    prtname = f"{name}_prt"
+    gwf_budget_file = f"{gwfname}.bud"
+    gwf_head_file = f"{gwfname}.hds"
+    prt_track_file = f"{prtname}.trk"
+    prt_track_csv_file = f"{prtname}.trk.csv"
 
-    # create gwf wel package
-    wells = [
-        # k, i, j, q
-        (0, 4, 4, -0.1),
-    ]
-    wel = flopy.mf6.ModflowGwfwel(
-        gwf,
-        maxbound=len(wells),
-        save_flows=True,
-        stress_period_data={0: wells, 1: wells},
-    )
-
-    # create gwf output control package
-    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
-
-
-def build_prt_sim(idx, ws, mf6):
     # create simulation
-    name = ex[idx]
     sim = flopy.mf6.MFSimulation(
-        sim_name=ex[idx],
+        sim_name=name,
         exe_name=mf6,
         version="mf6",
         sim_ws=ws,
     )
 
     # create tdis package
-    pd = (perlen, nstp, tsmult)
+    pd = (ctx.perlen, ctx.nstp, ctx.tsmult)
     flopy.mf6.modflow.mftdis.ModflowTdis(
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=nper,
+        nper=ctx.nper,
         perioddata=[pd],
     )
 
@@ -193,13 +93,13 @@ def build_prt_sim(idx, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=nlay,
-        nrow=nrow,
-        ncol=ncol,
+        nlay=ctx.nlay,
+        nrow=ctx.nrow,
+        ncol=ctx.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
 
     # create prp package
     flopy.mf6.ModflowPrtprp(
@@ -237,11 +137,13 @@ def build_prt_sim(idx, ws, mf6):
     )
     sim.register_solution_package(ems, [prt.name])
 
-    return sim
+    return sim, ctx
 
 
-def build_mp7_sim(idx, ws, mp7, gwf):
-    name = ex[idx]
+def build_mp7_sim(ctx, name, ws, mp7, gwf):
+    mp7name = f"{name}_mp7"
+    mp7_pathline_file = f"{mp7name}.mppth"
+
     partdata = flopy.modpath.ParticleData(
         partlocs=[p[0] for p in releasepts_mp7],
         localx=[p[1] for p in releasepts_mp7],
@@ -263,7 +165,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=porosity,
+        porosity=ctx.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -285,13 +187,25 @@ def get_different_rows(source_df, new_df):
     return changed_rows_df.drop("_merge", axis=1)
 
 
-@pytest.mark.parametrize("idx, name", enumerate(ex))
-def test_prt_fmi04(idx, name, function_tmpdir, targets):
+@pytest.mark.parametrize("name", ex)
+def test_mf6model(name, function_tmpdir, targets):
     ws = function_tmpdir
 
     # build mf6 models
-    gwfsim = build_gwf_sim(idx, ws, targets.mf6)
-    prtsim = build_prt_sim(idx, ws, targets.mf6)
+    gwfsim, ctx = build_gwf_sim(name, ws, targets.mf6)
+    # create gwf wel package
+    gwf = gwfsim.get_model()
+    wells = [
+        # k, i, j, q
+        (0, 4, 4, -0.1),
+    ]
+    wel = flopy.mf6.ModflowGwfwel(
+        gwf,
+        maxbound=len(wells),
+        save_flows=True,
+        stress_period_data={0: wells, 1: wells},
+    )
+    prtsim, ctx = build_prt_sim(ctx, name, ws, targets.mf6)
 
     # run mf6 models
     for sim in [gwfsim, prtsim]:
@@ -299,6 +213,16 @@ def test_prt_fmi04(idx, name, function_tmpdir, targets):
         success, _ = sim.run_simulation()
         assert success
 
+    # output file names
+    gwfname = f"{name}_gwf"
+    prtname = f"{name}_prt"
+    mp7name = f"{name}_mp7"
+    gwf_budget_file = f"{gwfname}.bud"
+    gwf_head_file = f"{gwfname}.hds"
+    prt_track_file = f"{prtname}.trk"
+    prt_track_csv_file = f"{prtname}.trk.csv"
+    mp7_pathline_file = f"{mp7name}.mppth"
+
     # extract model objects
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
@@ -307,7 +231,7 @@ def test_prt_fmi04(idx, name, function_tmpdir, targets):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(idx, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -353,7 +277,7 @@ def all_equal(col, val):
     assert all_equal(mf6_pls["iprp"], 1)
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{simname}_prt.lst", perlen, nper)
+    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
diff --git a/autotest/test_prt_fmi05.py b/autotest/prt/test_prt_fmi05.py
similarity index 78%
rename from autotest/test_prt_fmi05.py
rename to autotest/prt/test_prt_fmi05.py
index fcbdbd80cca..d835e1030df 100644
--- a/autotest/test_prt_fmi05.py
+++ b/autotest/prt/test_prt_fmi05.py
@@ -33,16 +33,11 @@
 import numpy as np
 import pandas as pd
 import pytest
+from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from flopy.plot.plotutil import to_mp7_pathlines
-from prt_test_utils import (
-    all_equal,
-    check_budget_data,
-    check_track_data,
-    has_default_boundnames,
-)
-
+from prt_test_utils import (all_equal, build_gwf_sim, check_budget_data,
+                            check_track_data, has_default_boundnames)
 
 # simulation name
 simname = "prtfmi05"
@@ -59,21 +54,9 @@
 
 
 # model names
-def get_model_name(idx, mdl):
-    return f"{ex[idx]}_{mdl}"
-
-
-# model info
-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
+def get_model_name(name, mdl):
+    return f"{name}_{mdl}"
+
 
 # release points in mp7 format (using local coordinates)
 releasepts_mp7 = [
@@ -105,7 +88,7 @@ def get_partdata(grid):
     )
 
 
-def build_gwf_sim(name, ws, mf6):
+def build_prt_sim(ctx, name, ws, mf6, fraction=None):
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
@@ -119,96 +102,25 @@ def build_gwf_sim(name, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=nper,
-        perioddata=[(perlen, nstp, tsmult)],
-    )
-
-    # create gwf model
-    gwfname = f"{name}_gwf"
-    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
-    # output file names
-    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
-
-
-def build_prt_sim(idx, ws, mf6, fraction=None):
-    # create simulation
-    sim = flopy.mf6.MFSimulation(
-        sim_name=ex[idx],
-        exe_name=mf6,
-        version="mf6",
-        sim_ws=ws,
-    )
-
-    # create tdis package
-    flopy.mf6.modflow.mftdis.ModflowTdis(
-        sim,
-        pname="tdis",
-        time_units="DAYS",
-        nper=nper,
-        perioddata=[(perlen, nstp, tsmult)],
+        nper=ctx.nper,
+        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
     )
 
     # create prt model
-    prtname = get_model_name(idx, "prt")
+    prtname = get_model_name(name, "prt")
     prt = flopy.mf6.ModflowPrt(sim, modelname=prtname)
 
     # create prt discretization
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=nlay,
-        nrow=nrow,
-        ncol=ncol,
+        nlay=ctx.nlay,
+        nrow=ctx.nrow,
+        ncol=ctx.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
 
     # convert mp7 particledata to prt release points
     partdata = get_partdata(prt.modelgrid)
@@ -264,7 +176,7 @@ def get_perioddata(name, periods=1) -> Optional[dict]:
     )
 
     # create the flow model interface
-    gwfname = get_model_name(idx, "gwf")
+    gwfname = get_model_name(name, "gwf")
     gwf_budget_file = f"{gwfname}.bud"
     gwf_head_file = f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
@@ -283,15 +195,15 @@ def get_perioddata(name, periods=1) -> Optional[dict]:
     )
     sim.register_solution_package(ems, [prt.name])
 
-    return sim
+    return sim, ctx
 
 
-def build_mp7_sim(idx, ws, mp7, gwf):
+def build_mp7_sim(ctx, name, ws, mp7, gwf):
     # convert mp7 particledata to prt release points
     partdata = get_partdata(gwf.modelgrid)
 
     # create modpath 7 simulation
-    mp7name = get_model_name(idx, "mp7")
+    mp7name = get_model_name(name, "mp7")
     pg = flopy.modpath.ParticleGroup(
         particlegroupname="G1",
         particledata=partdata,
@@ -305,7 +217,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=porosity,
+        porosity=ctx.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -319,23 +231,20 @@ def build_mp7_sim(idx, ws, mp7, gwf):
     return mp
 
 
-@pytest.mark.parametrize("idx, name", list(enumerate(ex)))
+@pytest.mark.parametrize("name", ex)
 @pytest.mark.parametrize("fraction", [0.5])
-def test_prt_fmi01(idx, name, function_tmpdir, targets, fraction):
+def test_mf6model(name, function_tmpdir, targets, fraction):
     # workspace
     ws = function_tmpdir
 
-    # test case name
-    name = ex[idx]
-
     # model names
-    gwfname = get_model_name(idx, "gwf")
-    prtname = get_model_name(idx, "prt")
-    mp7name = get_model_name(idx, "mp7")
+    gwfname = get_model_name(name, "gwf")
+    prtname = get_model_name(name, "prt")
+    mp7name = get_model_name(name, "mp7")
 
     # build mf6 models
-    gwfsim = build_gwf_sim(name, ws, targets.mf6)
-    prtsim = build_prt_sim(idx, ws, targets.mf6, fraction)
+    gwfsim, ctx = build_gwf_sim(name, ws, targets.mf6)
+    prtsim, ctx = build_prt_sim(ctx, name, ws, targets.mf6, fraction)
 
     # run mf6 models
     for sim in [gwfsim, prtsim]:
@@ -351,7 +260,7 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets, fraction):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(idx, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -401,7 +310,7 @@ def test_prt_fmi01(idx, name, function_tmpdir, targets, fraction):
     assert all_equal(mf6_pls["iprp"], 1)
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{name}_prt.lst", perlen, nper)
+    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
 
     # check mf6 prt particle track data were written to binary/CSV files
     # and that different formats are equal
diff --git a/autotest/test_prt_fmi06.py b/autotest/prt/test_prt_fmi06.py
similarity index 99%
rename from autotest/test_prt_fmi06.py
rename to autotest/prt/test_prt_fmi06.py
index fb269f8db57..ae82cd7185f 100644
--- a/autotest/test_prt_fmi06.py
+++ b/autotest/prt/test_prt_fmi06.py
@@ -264,7 +264,7 @@ def build_mp7_sim(idx, ws, mp7, gwf):
 
 
 @pytest.mark.parametrize("idx, name", list(enumerate(ex)))
-def test_prt_fmi06(idx, name, function_tmpdir, targets):
+def test_mf6model(idx, name, function_tmpdir, targets):
     # workspace
     ws = function_tmpdir
 
diff --git a/autotest/test_prt_notebooks.py b/autotest/prt/test_prt_notebooks.py
similarity index 100%
rename from autotest/test_prt_notebooks.py
rename to autotest/prt/test_prt_notebooks.py
index 1c78ce02171..610e9a7fb2d 100644
--- a/autotest/test_prt_notebooks.py
+++ b/autotest/prt/test_prt_notebooks.py
@@ -5,16 +5,16 @@
 from pprint import pprint
 from warnings import warn
 
-import pytest
 import numpy as np
 import pandas as pd
-from numpy.lib.recfunctions import stack_arrays
+import pytest
 from flaky import flaky
-from flopy.utils import PathlineFile
 from flopy.mf6 import MFSimulation
 from flopy.plot.plotutil import to_mp7_pathlines
+from flopy.utils import PathlineFile
 from modflow_devtools.markers import excludes_platform
 from modflow_devtools.misc import run_cmd, set_env
+from numpy.lib.recfunctions import stack_arrays
 
 from conftest import project_root_path