diff --git a/autotest/framework.py b/autotest/framework.py
index 961a97c70ef..d429d7da121 100644
--- a/autotest/framework.py
+++ b/autotest/framework.py
@@ -1,6 +1,7 @@
 import os
 import shutil
 import time
+from itertools import repeat
 from pathlib import Path
 from subprocess import PIPE, STDOUT, Popen
 from traceback import format_exc
@@ -313,13 +314,7 @@ def __init__(
         workspace = Path(workspace).expanduser().absolute()
         assert workspace.is_dir(), f"{workspace} is not a valid directory"
         if verbose:
-            from pprint import pprint
-
             print("Initializing test", name, "in workspace", workspace)
-            contents = list(workspace.glob("*"))
-            if any(contents):
-                print(f"Workspace is not empty:")
-                pprint(contents)
 
         self.name = name
         self.workspace = workspace
@@ -339,7 +334,7 @@ def __init__(
         self.rclose = 0.001 if rclose is None else rclose
         self.overwrite = overwrite
         self.verbose = verbose
-        self.xfail = xfail
+        self.xfail = [xfail] if isinstance(xfail, bool) else xfail
 
     def __repr__(self):
         return self.name
@@ -632,7 +627,7 @@ def run_sim_or_model(
         workspace: Union[str, os.PathLike],
         target: Union[str, os.PathLike] = "mf6",
         xfail: bool = False,
-    ) -> bool:
+    ) -> Tuple[bool, List[str]]:
         """
         Run a simulation or model with FloPy.
 
@@ -660,14 +655,16 @@ def run_sim_or_model(
         self.cmp_namefile = (
             None
             if "mf6" in target.name or "libmf6" in target.name
-            else os.path.basename(nf) if nf else None
+            else os.path.basename(nf)
+            if nf
+            else None
         )
 
         # run the model
         try:
             # via MODFLOW API
             if "libmf6" in target.name and self.api_func:
-                success, _ = self.api_func(target, workspace)
+                success, buff = self.api_func(target, workspace)
             # via MF6 executable
             elif "mf6" in target.name:
                 # parallel test if configured
@@ -676,7 +673,7 @@ def run_sim_or_model(
                         f"Parallel test {self.name} on {self.ncpus} processes"
                     )
                     try:
-                        success, _ = run_parallel(
+                        success, buff = run_parallel(
                             workspace, target, self.ncpus
                         )
                     except Exception:
@@ -689,10 +686,11 @@ def run_sim_or_model(
                 else:
                     # otherwise serial run
                     try:
-                        success, _ = flopy.run_model(
+                        success, buff = flopy.run_model(
                             target,
                             self.workspace / "mfsim.nam",
                             model_ws=workspace,
+                            report=True,
                         )
                     except Exception:
                         warn(
@@ -704,8 +702,8 @@ def run_sim_or_model(
             else:
                 # non-MF6 model
                 try:
-                    success, _ = flopy.run_model(
-                        target, self.cmp_namefile, workspace
+                    success, buff = flopy.run_model(
+                        target, self.cmp_namefile, workspace, report=True
                     )
                 except Exception:
                     warn(f"{target} model failed:\n{format_exc()}")
@@ -734,7 +732,7 @@ def run_sim_or_model(
                 f"Unhandled error in comparison model {self.name}:\n{format_exc()}"
             )
 
-        return success
+        return success, buff
 
     def compare_output(self, compare):
         """
@@ -792,17 +790,33 @@ def run(self):
             sims = sims if isinstance(sims, Iterable) else [sims]
             sims = [sim for sim in sims if sim]  # filter Nones
             self.sims = sims
+            nsims = len(sims)
+            self.buffs = list(repeat(None, nsims))
+            assert len(self.xfail) in [
+                1,
+                nsims,
+            ], f"Invalid xfail: expected a single boolean or one for each model"
+            if len(self.xfail) == 1 and nsims:
+                self.xfail = list(repeat(self.xfail[0], nsims))
             write_input(*sims, overwrite=self.overwrite, verbose=self.verbose)
         else:
             self.sims = [MFSimulation.load(sim_ws=self.workspace)]
+            self.buffs = [None]
+            assert (
+                len(self.xfail) == 1
+            ), f"Invalid xfail: expected a single boolean or one for each model"
 
         # run models/simulations
-        for sim_or_model in self.sims:
+        for i, sim_or_model in enumerate(self.sims):
             workspace = get_workspace(sim_or_model)
             target = self._try_resolve(sim_or_model.exe_name, self.targets.mf6)
-            assert self.run_sim_or_model(
-                workspace, target, self.xfail
-            ), f"{'Simulation' if 'mf6' in str(target) else 'Model'} failed: {workspace}"
+            xfail = self.xfail[i]
+            success, buff = self.run_sim_or_model(workspace, target, xfail)
+            self.buffs[i] = buff  # store model output for assertions later
+            assert success, (
+                f"{'Simulation' if 'mf6' in str(target) else 'Model'} "
+                f"{'should have failed' if xfail else 'failed'}: {workspace}"
+            )
 
         # get expected output files from main simulation
         _, self.outp = get_mf6_files(
@@ -842,10 +856,11 @@ def run(self):
                     # todo: don't hardcode workspace or assume agreement with test case
                     # simulation workspace, set & access simulation workspaces directly
                     workspace = self.workspace / self.compare
-                    assert self.run_sim_or_model(
+                    success, _ = self.run_sim_or_model(
                         workspace,
                         self.targets.get(self.compare, self.targets.mf6),
-                    ), f"Comparison model failed: {workspace}"
+                    )
+                    assert success, f"Comparison model failed: {workspace}"
 
                 # compare model results, if enabled
                 if self.verbose:
diff --git a/autotest/prt/prt_test_utils.py b/autotest/prt/prt_test_utils.py
index b8743ccdc6e..9fa75bc9e4b 100644
--- a/autotest/prt/prt_test_utils.py
+++ b/autotest/prt/prt_test_utils.py
@@ -1,6 +1,6 @@
 import os
 from types import SimpleNamespace
-from typing import Optional, Tuple
+from typing import Tuple
 
 import flopy
 import matplotlib as mpl
@@ -12,107 +12,104 @@ def all_equal(series, val):
     return a[0] == val and (a[0] == a).all()
 
 
-def get_gwf_sim(
-    name, ws, mf6
-) -> Tuple[flopy.mf6.MFSimulation, SimpleNamespace]:
-    """
-    Simple GWF simulation for use/modification by PRT tests
-    """
-
-    # test case context
-    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,
-        # mp7 release points (cell-local coordinates)
-        releasepts_mp7=[
-            # node number, localx, localy, localz
-            (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5)
-            for i in range(9)
-        ],
-        # PRT release points (cell indices and global coordinates both required)
-        releasepts_prt=[
-            # particle index, k, i, j, x, y, z
-            [i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5]
-            for i in range(9)
-        ],
-    )
+class BasicDisCase:
+    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
+    releasepts_mp7 = [
+        # node number, localx, localy, localz
+        (0, float(f"0.{i + 1}"), float(f"0.{i + 1}"), 0.5)
+        for i in range(9)
+    ]
+    releasepts_prt = [
+        # particle index, k, i, j, x, y, z
+        [i, 0, 0, 0, float(f"0.{i + 1}"), float(f"9.{i + 1}"), 0.5]
+        for i in range(9)
+    ]
 
-    # create simulation
-    sim = flopy.mf6.MFSimulation(
-        sim_name=name,
-        exe_name=mf6,
-        version="mf6",
-        sim_ws=ws,
-    )
+    @staticmethod
+    def get_gwf_sim(name, ws, mf6) -> flopy.mf6.MFSimulation:
+        """
+        Simple GWF simulation for use/modification by PRT tests
+        """
+
+        # 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 tdis package
+        flopy.mf6.modflow.mftdis.ModflowTdis(
+            sim,
+            pname="tdis",
+            time_units="DAYS",
+            nper=BasicDisCase.nper,
+            perioddata=[
+                (BasicDisCase.perlen, BasicDisCase.nstp, BasicDisCase.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 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=BasicDisCase.nlay,
+            nrow=BasicDisCase.nrow,
+            ncol=BasicDisCase.ncol,
+        )
 
-    # create gwf initial conditions package
-    flopy.mf6.modflow.mfgwfic.ModflowGwfic(gwf, pname="ic")
+        # 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 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 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 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)
+        # create iterative model solution for gwf model
+        ims = flopy.mf6.ModflowIms(sim)
 
-    return sim, ctx
+        return sim
 
 
 def check_track_data(
diff --git a/autotest/prt/test_prt_exg01.py b/autotest/prt/test_prt_exg01.py
index 813ec90ee82..3ad60b2e3d4 100644
--- a/autotest/prt/test_prt_exg01.py
+++ b/autotest/prt/test_prt_exg01.py
@@ -25,23 +25,22 @@
 from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from prt_test_utils import check_budget_data, check_track_data, get_gwf_sim
+from prt_test_utils import BasicDisCase, check_budget_data, check_track_data
 
 from framework import TestFramework
 
 simname = "prtexg01"
-ex = [simname, f"{simname}bnms"]
+cases = [simname, f"{simname}bnms"]
 
 
-# model names
 def get_model_name(idx, mdl):
-    return f"{ex[idx]}_{mdl}"
+    return f"{cases[idx]}_{mdl}"
 
 
-def build_sim(idx, ws, mf6):
+def build_models(idx, test):
     # create simulation
-    name = ex[idx]
-    sim, ctx = get_gwf_sim(name, ws, mf6)
+    name = cases[idx]
+    sim = BasicDisCase.get_gwf_sim(name, test.workspace, test.targets.mf6)
 
     # create prt model
     prtname = get_model_name(idx, "prt")
@@ -51,19 +50,19 @@ def build_sim(idx, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=ctx.nlay,
-        nrow=ctx.nrow,
-        ncol=ctx.ncol,
+        nlay=BasicDisCase.nlay,
+        nrow=BasicDisCase.nrow,
+        ncol=BasicDisCase.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=BasicDisCase.porosity)
 
     # create prp package
     rpts = (
-        [r + [str(r[0] + 1)] for r in ctx.releasepts_prt]
+        [r + [str(r[0] + 1)] for r in BasicDisCase.releasepts_prt]
         if "bnms" in name
-        else ctx.releasepts_prt
+        else BasicDisCase.releasepts_prt
     )
     flopy.mf6.ModflowPrtprp(
         prt,
@@ -113,15 +112,15 @@ def build_sim(idx, ws, mf6):
     )
     sim.register_solution_package(ems, [prt.name])
 
-    return sim, ctx
+    return sim
 
 
-def build_mp7_sim(ctx, idx, ws, mp7, gwf):
+def build_mp7_sim(idx, ws, mp7, gwf):
     partdata = flopy.modpath.ParticleData(
-        partlocs=[p[0] for p in ctx.releasepts_mp7],
-        localx=[p[1] for p in ctx.releasepts_mp7],
-        localy=[p[2] for p in ctx.releasepts_mp7],
-        localz=[p[3] for p in ctx.releasepts_mp7],
+        partlocs=[p[0] for p in BasicDisCase.releasepts_mp7],
+        localx=[p[1] for p in BasicDisCase.releasepts_mp7],
+        localy=[p[2] for p in BasicDisCase.releasepts_mp7],
+        localz=[p[3] for p in BasicDisCase.releasepts_mp7],
         timeoffset=0,
         drape=0,
     )
@@ -139,7 +138,7 @@ def build_mp7_sim(ctx, idx, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=ctx.porosity,
+        porosity=BasicDisCase.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -153,41 +152,9 @@ def build_mp7_sim(ctx, idx, ws, mp7, gwf):
     return mp
 
 
-def eval_results(ctx, test):
-    print(f"Evaluating results for sim {test.name}")
-    simpath = Path(test.workspace)
-
-    # check budget data
-    check_budget_data(simpath / f"{test.name}_prt.lst", ctx.perlen, ctx.nper)
-
-    # check particle track data
-    prt_track_file = simpath / f"{test.name}_prt.trk"
-    prt_track_hdr_file = simpath / f"{test.name}_prt.trk.hdr"
-    prt_track_csv_file = simpath / f"{test.name}_prt.trk.csv"
-    assert prt_track_file.exists()
-    assert prt_track_hdr_file.exists()
-    assert prt_track_csv_file.exists()
-    check_track_data(
-        track_bin=prt_track_file,
-        track_hdr=prt_track_hdr_file,
-        track_csv=prt_track_csv_file,
-    )
-
-
-@pytest.mark.parametrize("idx, name", enumerate(ex))
-def test_mf6model(idx, name, function_tmpdir, targets):
-    ws = function_tmpdir
-    sim, ctx = build_sim(idx, str(ws), targets.mf6)
-    sim.write_simulation()
-
-    test = TestFramework(
-        name=name,
-        workspace=ws,
-        targets=targets,
-        check=lambda s: eval_results(ctx, s),
-        compare=None,
-    )
-    test.run()
+def check_output(idx, test):
+    name = test.name
+    ws = Path(test.workspace)
 
     # model names
     gwfname = get_model_name(idx, "gwf")
@@ -195,6 +162,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     mp7name = get_model_name(idx, "mp7")
 
     # extract model objects
+    sim = test.sims[0]
     gwf = sim.get_model(gwfname)
     prt = sim.get_model(prtname)
 
@@ -202,7 +170,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(ctx, idx, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(idx, ws, test.targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -252,7 +220,9 @@ 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", ctx.perlen, ctx.nper)
+    check_budget_data(
+        ws / f"{name}_prt.lst", BasicDisCase.perlen, BasicDisCase.nper
+    )
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
@@ -333,3 +303,16 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     # compare mf6 / mp7 pathline data
     assert mf6_pls.shape == mp7_pls.shape
     assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi01.py b/autotest/prt/test_prt_fmi01.py
index b97909914af..abb6ec84e63 100644
--- a/autotest/prt/test_prt_fmi01.py
+++ b/autotest/prt/test_prt_fmi01.py
@@ -12,7 +12,7 @@
 Two test cases are defined, one with the particle
 release (PRP) package option STOP_AT_WEAK_SINK on
 and one with the option off. No effect on results
-are expected because the model has no weak sinks.
+is expected, because the model has no weak sinks.
 (Motivated by an old bug in which particles were
 tracked improperly when this option was enabled,
 even with no weak sink cells in the vicinity.)
@@ -36,21 +36,29 @@
 from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from prt_test_utils import (all_equal, check_budget_data, check_track_data,
-                            get_gwf_sim, get_model_name, get_partdata,
-                            has_default_boundnames)
+from prt_test_utils import (
+    BasicDisCase,
+    all_equal,
+    check_budget_data,
+    check_track_data,
+    get_model_name,
+    get_partdata,
+    has_default_boundnames,
+)
+
+from framework import TestFramework
 
 simname = "prtfmi01"
-ex = [simname, f"{simname}saws"]
+cases = [simname, f"{simname}saws"]
 
 
-def build_prt_sim(ctx, name, ws, mf6):
+def build_prt_sim(name, gwf_ws, prt_ws, mf6):
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
@@ -58,8 +66,10 @@ def build_prt_sim(ctx, name, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=ctx.nper,
-        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
+        nper=BasicDisCase.nper,
+        perioddata=[
+            (BasicDisCase.perlen, BasicDisCase.nstp, BasicDisCase.tsmult)
+        ],
     )
 
     # create prt model
@@ -70,19 +80,19 @@ def build_prt_sim(ctx, name, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=ctx.nlay,
-        nrow=ctx.nrow,
-        ncol=ctx.ncol,
+        nlay=BasicDisCase.nlay,
+        nrow=BasicDisCase.nrow,
+        ncol=BasicDisCase.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=BasicDisCase.porosity)
 
     # convert mp7 to prt release points and check against expectation
-    partdata = get_partdata(prt.modelgrid, ctx.releasepts_mp7)
+    partdata = get_partdata(prt.modelgrid, BasicDisCase.releasepts_mp7)
     coords = partdata.to_coords(prt.modelgrid)
     releasepts = [(i, 0, 0, 0, c[0], c[1], c[2]) for i, c in enumerate(coords)]
-    assert np.allclose(ctx.releasepts_prt, releasepts)
+    assert np.allclose(BasicDisCase.releasepts_prt, releasepts)
 
     # create prp package
     prp_track_file = f"{prtname}.prp.trk"
@@ -112,8 +122,8 @@ def build_prt_sim(ctx, name, ws, mf6):
 
     # create the flow model interface
     gwfname = get_model_name(name, "gwf")
-    gwf_budget_file = f"{gwfname}.bud"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.bud"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -133,11 +143,8 @@ def build_prt_sim(ctx, name, ws, mf6):
     return sim
 
 
-def build_mp7_sim(ctx, name, ws, mp7, gwf):
-    # convert mp7 particledata to prt release points
-    partdata = get_partdata(gwf.modelgrid, ctx.releasepts_mp7)
-
-    # create modpath 7 simulation
+def build_mp7_sim(name, ws, mp7, gwf):
+    partdata = get_partdata(gwf.modelgrid, BasicDisCase.releasepts_mp7)
     mp7name = get_model_name(name, "mp7")
     pg = flopy.modpath.ParticleGroup(
         particlegroupname="G1",
@@ -152,7 +159,7 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=ctx.porosity,
+        porosity=BasicDisCase.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -166,33 +173,33 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     return mp
 
 
-@pytest.mark.parametrize("name", ex)
-def test_mf6model(name, function_tmpdir, targets):
-    # workspace
-    ws = function_tmpdir
+def build_models(idx, test):
+    gwfsim = BasicDisCase.get_gwf_sim(test.name, test.workspace, test.targets.mf6)
+    prtsim = build_prt_sim(
+        test.name, test.workspace, test.workspace / "prt", test.targets.mf6
+    )
+    return gwfsim, prtsim
+
 
-    # model names
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    mp7_ws = test.workspace / "mp7"
     gwfname = get_model_name(name, "gwf")
     prtname = get_model_name(name, "prt")
     mp7name = get_model_name(name, "mp7")
 
-    # build mf6 models
-    gwfsim, ctx = get_gwf_sim(name, ws, targets.mf6)
-    prtsim = build_prt_sim(ctx, name, ws, targets.mf6)
-
-    # run mf6 models
-    for sim in [gwfsim, prtsim]:
-        sim.write_simulation()
-        success, buff = sim.run_simulation(report=True)
-        assert success, pformat(buff)
-
-    # extract mf6 models and grid
+    # extract mf6 simulations/models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
     mg = gwf.modelgrid
 
-    # build mp7 model
-    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
+    # build/run mp7 model... can't run within framework as
+    # flopy needs gwf output files to write mp7 input files
+    mp7sim = build_mp7_sim(name, mp7_ws, test.targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -206,19 +213,19 @@ def test_mf6model(name, function_tmpdir, targets):
     prt_track_csv_file = f"{prtname}.trk.csv"
     prp_track_file = f"{prtname}.prp.trk"
     prp_track_csv_file = f"{prtname}.prp.trk.csv"
-    assert (ws / gwf_budget_file).is_file()
-    assert (ws / gwf_head_file).is_file()
-    assert (ws / prt_track_file).is_file()
-    assert (ws / prt_track_csv_file).is_file()
-    assert (ws / prp_track_file).is_file()
-    assert (ws / prp_track_csv_file).is_file()
+    assert (gwf_ws / gwf_budget_file).is_file()
+    assert (gwf_ws / gwf_head_file).is_file()
+    assert (prt_ws / prt_track_file).is_file()
+    assert (prt_ws / prt_track_csv_file).is_file()
+    assert (prt_ws / prp_track_file).is_file()
+    assert (prt_ws / prp_track_csv_file).is_file()
 
     # check mp7 output files exist
     mp7_pathline_file = f"{mp7name}.mppth"
-    assert (ws / mp7_pathline_file).is_file()
+    assert (mp7_ws / mp7_pathline_file).is_file()
 
     # load mp7 pathline results
-    plf = PathlineFile(ws / mp7_pathline_file)
+    plf = PathlineFile(mp7_ws / mp7_pathline_file)
     mp7_pls = pd.DataFrame(
         plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
     )
@@ -228,7 +235,7 @@ def test_mf6model(name, function_tmpdir, targets):
     mp7_pls["k"] = mp7_pls["k"] + 1
 
     # load mf6 pathline results
-    mf6_pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False)
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
 
     # make sure pathline df has "name" (boundname) column and default values
     assert "name" in mf6_pls
@@ -239,19 +246,25 @@ def test_mf6model(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", ctx.perlen, ctx.nper)
+    check_budget_data(
+        prt_ws / f"{name}_prt.lst", BasicDisCase.perlen, BasicDisCase.nper
+    )
 
     # check mf6 prt particle track data were written to binary/CSV files
     # and that different formats are equal
-    for track_csv in [ws / prt_track_csv_file, ws / prp_track_csv_file]:
+    for track_csv in [
+        prt_ws / prt_track_csv_file,
+        prt_ws / prp_track_csv_file,
+    ]:
         check_track_data(
-            track_bin=ws / prt_track_file,
-            track_hdr=ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
+            track_bin=prt_ws / prt_track_file,
+            track_hdr=prt_ws
+            / Path(prt_track_file.replace(".trk", ".trk.hdr")),
             track_csv=track_csv,
         )
 
     # extract head, budget, and specific discharge results from GWF model
-    hds = HeadFile(ws / gwf_head_file).get_data()
+    hds = HeadFile(gwf_ws / gwf_head_file).get_data()
     bud = gwf.output.budget()
     spdis = bud.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
@@ -297,7 +310,7 @@ def test_mf6model(name, function_tmpdir, targets):
 
     # view/save plot
     # plt.show()
-    plt.savefig(ws / f"test_{simname}.png")
+    plt.savefig(gwf_ws / f"test_{simname}.png")
 
     # convert mf6 pathlines to mp7 format
     mf6_pls = to_mp7_pathlines(mf6_pls)
@@ -321,3 +334,16 @@ def test_mf6model(name, function_tmpdir, targets):
     # compare mf6 / mp7 pathline data
     assert mf6_pls.shape == mp7_pls.shape
     assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi02.py b/autotest/prt/test_prt_fmi02.py
index 59174e755ba..d3169835522 100644
--- a/autotest/prt/test_prt_fmi02.py
+++ b/autotest/prt/test_prt_fmi02.py
@@ -36,11 +36,17 @@
 from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from prt_test_utils import (check_budget_data, check_track_data, get_gwf_sim,
-                            get_model_name)
+from prt_test_utils import (
+    BasicDisCase,
+    check_budget_data,
+    check_track_data,
+    get_model_name,
+)
+
+from framework import TestFramework
 
 simname = "prtfmi02"
-ex = [
+cases = [
     f"{simname}all",
     f"{simname}rel",
     f"{simname}trst",
@@ -99,13 +105,13 @@ def create_idomain(nlay, nrow, ncol):
     return idmn
 
 
-def build_prt_sim(ctx, name, ws, mf6):
+def build_prt_sim(name, gwf_ws, prt_ws, mf6):
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
@@ -113,8 +119,10 @@ def build_prt_sim(ctx, name, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=ctx.nper,
-        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
+        nper=BasicDisCase.nper,
+        perioddata=[
+            (BasicDisCase.perlen, BasicDisCase.nstp, BasicDisCase.tsmult)
+        ],
     )
 
     # create prt model
@@ -125,14 +133,16 @@ def build_prt_sim(ctx, name, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=ctx.nlay,
-        nrow=ctx.nrow,
-        ncol=ctx.ncol,
-        idomain=create_idomain(ctx.nlay, ctx.nrow, ctx.ncol),
+        nlay=BasicDisCase.nlay,
+        nrow=BasicDisCase.nrow,
+        ncol=BasicDisCase.ncol,
+        idomain=create_idomain(
+            BasicDisCase.nlay, BasicDisCase.nrow, BasicDisCase.ncol
+        ),
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=BasicDisCase.porosity)
 
     # create a prp package for groups a and b
     prps = [
@@ -161,8 +171,8 @@ def build_prt_sim(ctx, name, ws, mf6):
 
     # create the flow model interface
     gwfname = get_model_name(name, "gwf")
-    gwf_budget_file = f"{gwfname}.bud"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.bud"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -182,7 +192,7 @@ def build_prt_sim(ctx, name, ws, mf6):
     return sim
 
 
-def build_mp7_sim(ctx, name, ws, mp7, gwf):
+def build_mp7_sim(name, ws, mp7, gwf):
     mp7name = get_model_name(name, "mp7")
     mp7_pathline_file = f"{mp7name}.mppth"
     pgs = [
@@ -208,7 +218,7 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=ctx.porosity,
+        porosity=BasicDisCase.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -222,34 +232,39 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     return mp
 
 
-@pytest.mark.parametrize("name", ex)
-def test_mf6model(name, function_tmpdir, targets):
-    # workspace
-    ws = function_tmpdir
-
-    # model names
-    gwfname = get_model_name(name, "gwf")
-    prtname = get_model_name(name, "prt")
-    mp7name = get_model_name(name, "mp7")
-
-    # build mf6 model
-    gwfsim, ctx = get_gwf_sim(name, ws, targets.mf6)
-
+def build_models(idx, test):
+    # build gwf model
+    gwfsim = BasicDisCase.get_gwf_sim(
+        test.name, test.workspace, test.targets.mf6
+    )
     # add idomain
     gwf = gwfsim.get_model()
     dis = gwf.get_package("DIS")
-    dis.idomain = create_idomain(ctx.nlay, ctx.nrow, ctx.ncol)
+    dis.idomain = create_idomain(
+        BasicDisCase.nlay, BasicDisCase.nrow, BasicDisCase.ncol
+    )
 
     # build prt model
-    prtsim = build_prt_sim(ctx, name, ws, targets.mf6)
+    prtsim = build_prt_sim(
+        test.name, test.workspace, test.workspace / "prt", test.targets.mf6
+    )
+    return gwfsim, prtsim
+
 
-    # run mf6 models
-    for sim in [gwfsim, prtsim]:
-        sim.write_simulation()
-        success, buff = sim.run_simulation(report=True)
-        assert success, pformat(buff)
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    mp7_ws = test.workspace / "mp7"
+
+    # model names
+    gwfname = get_model_name(name, "gwf")
+    prtname = get_model_name(name, "prt")
+    mp7name = get_model_name(name, "mp7")
 
     # extract models
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
 
@@ -257,7 +272,7 @@ def test_mf6model(name, function_tmpdir, targets):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(name, mp7_ws, test.targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -270,16 +285,16 @@ def test_mf6model(name, function_tmpdir, targets):
     prt_track_file = f"{prtname}.trk"
     prt_track_csv_file = f"{prtname}.trk.csv"
     mp7_pathline_file = f"{mp7name}.mppth"
-    assert (ws / gwf_budget_file).is_file()
-    assert (ws / gwf_head_file).is_file()
-    assert (ws / prt_track_file).is_file()
-    assert (ws / prt_track_csv_file).is_file()
+    assert (gwf_ws / gwf_budget_file).is_file()
+    assert (gwf_ws / gwf_head_file).is_file()
+    assert (prt_ws / prt_track_file).is_file()
+    assert (prt_ws / prt_track_csv_file).is_file()
 
     # check mp7 output files exist
-    assert (ws / mp7_pathline_file).is_file()
+    assert (mp7_ws / mp7_pathline_file).is_file()
 
     # load mp7 pathline results
-    plf = PathlineFile(ws / mp7_pathline_file)
+    plf = PathlineFile(mp7_ws / mp7_pathline_file)
     mp7_pls = pd.DataFrame(
         plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
     )
@@ -289,7 +304,7 @@ def test_mf6model(name, function_tmpdir, targets):
     mp7_pls["k"] = mp7_pls["k"] + 1
 
     # load mf6 pathline results
-    mf6_pls = pd.read_csv(ws / prt_track_csv_file)
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file)
 
     # if event is ALL, output should be the same as MODPATH 7,
     # so continue with comparisons.
@@ -329,20 +344,22 @@ def all_equal(col, val):
     )
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
+    check_budget_data(
+        prt_ws / f"{name}_prt.lst", BasicDisCase.perlen, BasicDisCase.nper
+    )
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
-        track_bin=ws / prt_track_file,
-        track_hdr=ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
-        track_csv=ws / prt_track_csv_file,
+        track_bin=prt_ws / prt_track_file,
+        track_hdr=prt_ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
+        track_csv=prt_ws / prt_track_csv_file,
     )
 
     # check that particle names are particle indices
     # assert len(mf6_pldata) == len(mf6_pldata[mf6_pldata['irpt'].astype(str).eq(mf6_pldata['name'])])
 
     # get head, budget, and spdis results from GWF model
-    hds = HeadFile(ws / gwf_head_file).get_data()
+    hds = HeadFile(gwf_ws / gwf_head_file).get_data()
     bud = gwf.output.budget()
     spdis = bud.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
@@ -388,7 +405,7 @@ def all_equal(col, val):
 
     # view/save plot
     # plt.show()
-    plt.savefig(ws / f"test_{simname}.png")
+    plt.savefig(gwf_ws / f"test_{simname}.png")
 
     # check that cell numbers are correct
     for i, row in list(mf6_pls.iterrows()):
@@ -437,3 +454,16 @@ def all_equal(col, val):
     # compare mf6 / mp7 pathline data
     assert mf6_pls.shape == mp7_pls.shape
     assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi03.py b/autotest/prt/test_prt_fmi03.py
index b593c670b9b..fc155001427 100644
--- a/autotest/prt/test_prt_fmi03.py
+++ b/autotest/prt/test_prt_fmi03.py
@@ -37,11 +37,17 @@
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
 from matplotlib.collections import LineCollection
-from prt_test_utils import (check_budget_data, check_track_data, get_gwf_sim,
-                            get_model_name)
+from prt_test_utils import (
+    BasicDisCase,
+    check_budget_data,
+    check_track_data,
+    get_model_name,
+)
+
+from framework import TestFramework
 
 simname = "prtfmi03"
-ex = [f"{simname}_l1", f"{simname}_l2"]
+cases = [f"{simname}_l1", f"{simname}_l2"]
 stopzone_cells = [(0, 1, 8), (0, 8, 1)]
 
 
@@ -52,13 +58,13 @@ def create_izone(nlay, nrow, ncol):
     return izone
 
 
-def build_prt_sim(ctx, name, ws, mf6):
+def build_prt_sim(name, gwf_ws, prt_ws, mf6):
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
@@ -66,8 +72,10 @@ def build_prt_sim(ctx, name, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=ctx.nper,
-        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
+        nper=BasicDisCase.nper,
+        perioddata=[
+            (BasicDisCase.perlen, BasicDisCase.nstp, BasicDisCase.tsmult)
+        ],
     )
 
     # create prt model
@@ -75,25 +83,25 @@ def build_prt_sim(ctx, name, ws, mf6):
     prt = flopy.mf6.ModflowPrt(sim, modelname=prtname)
 
     # create prt discretization
-    ctx.nlay = int(name[-1])
-    botm = [ctx.top - (k + 1) for k in range(ctx.nlay)]
+    nlay = int(name[-1])
+    botm = [BasicDisCase.top - (k + 1) for k in range(nlay)]
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=ctx.nlay,
-        nrow=ctx.nrow,
-        ncol=ctx.ncol,
-        top=ctx.top,
+        nlay=nlay,
+        nrow=BasicDisCase.nrow,
+        ncol=BasicDisCase.ncol,
+        top=BasicDisCase.top,
         botm=botm,
     )
 
     # create mip package
-    izone = create_izone(ctx.nlay, ctx.nrow, ctx.ncol)
+    izone = create_izone(nlay, BasicDisCase.nrow, BasicDisCase.ncol)
     flopy.mf6.ModflowPrtmip(
         prt,
         pname="mip",
-        porosity=ctx.porosity,
-        izone=izone,  # if ctx.nlay == 1 else np.array([izone, izone]),
+        porosity=BasicDisCase.porosity,
+        izone=izone,
     )
 
     # create prp package
@@ -101,8 +109,8 @@ def build_prt_sim(ctx, name, ws, mf6):
         prt,
         pname="prp1",
         filename=f"{prtname}_1.prp",
-        nreleasepts=len(ctx.releasepts_prt),
-        packagedata=ctx.releasepts_prt,
+        nreleasepts=len(BasicDisCase.releasepts_prt),
+        packagedata=BasicDisCase.releasepts_prt,
         perioddata={0: ["FIRST"]},
         istopzone=1,
     )
@@ -119,8 +127,8 @@ def build_prt_sim(ctx, name, ws, mf6):
 
     # create the flow model interface
     gwfname = get_model_name(name, "gwf")
-    gwf_budget_file = f"{gwfname}.bud"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.bud"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -137,15 +145,15 @@ def build_prt_sim(ctx, name, ws, mf6):
     )
     sim.register_solution_package(ems, [prt.name])
 
-    return sim, ctx
+    return sim
 
 
-def build_mp7_sim(ctx, name, ws, mp7, gwf):
+def build_mp7_sim(name, ws, mp7, gwf):
     partdata = flopy.modpath.ParticleData(
-        partlocs=[p[0] for p in ctx.releasepts_mp7],
-        localx=[p[1] for p in ctx.releasepts_mp7],
-        localy=[p[2] for p in ctx.releasepts_mp7],
-        localz=[p[3] for p in ctx.releasepts_mp7],
+        partlocs=[p[0] for p in BasicDisCase.releasepts_mp7],
+        localx=[p[1] for p in BasicDisCase.releasepts_mp7],
+        localy=[p[2] for p in BasicDisCase.releasepts_mp7],
+        localz=[p[3] for p in BasicDisCase.releasepts_mp7],
         timeoffset=0,
         drape=0,
     )
@@ -163,10 +171,10 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=ctx.porosity,
+        porosity=BasicDisCase.porosity,
     )
-    ctx.nlay = int(name[-1])
-    izone = create_izone(ctx.nlay, ctx.nrow, ctx.ncol)
+    nlay = int(name[-1])
+    izone = create_izone(nlay, BasicDisCase.nrow, BasicDisCase.ncol)
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
         simulationtype="pathline",
@@ -174,7 +182,7 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
         budgetoutputoption="summary",
         stoptimeoption="extend",
         stopzone=1,
-        zones=izone,  # if ctx.nlay == 1 else np.array([izone, izone]),
+        zones=izone,
         zonedataoption="on",
         particlegroups=[pg],
     )
@@ -182,43 +190,45 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     return mp
 
 
-@pytest.mark.parametrize("name", ex)
-def test_mf6model(name, function_tmpdir, targets):
-    ws = function_tmpdir
-
-    # define model names
-    gwfname = get_model_name(name, "gwf")
-    prtname = get_model_name(name, "prt")
-    mp7name = get_model_name(name, "mp7")
-
-    # build mf6 simulations
-    gwfsim, ctx = get_gwf_sim(name, ws, targets.mf6)
+def build_models(idx, test):
+    gwfsim = BasicDisCase.get_gwf_sim(
+        test.name, test.workspace, test.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)]
+    nlay = int(test.name[-1])
+    botm = [BasicDisCase.top - (k + 1) for k in range(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
+        [list(repeat(b, BasicDisCase.nrow * BasicDisCase.ncol)) for b in botm]
+    ).reshape((nlay, BasicDisCase.nrow, BasicDisCase.ncol))
+    dis.nlay = nlay
     dis.botm.set_data(botm_data)
-    prtsim, ctx = build_prt_sim(ctx, name, ws, targets.mf6)
+    prtsim = build_prt_sim(
+        test.name, test.workspace, test.workspace / "prt", test.targets.mf6
+    )
+    return gwfsim, prtsim
 
-    # run mf6 simulations
-    for sim in [gwfsim, prtsim]:
-        sim.write_simulation()
-        success, buff = sim.run_simulation(report=True)
-        assert success, pformat(buff)
 
-    # extract mf6 models
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    mp7_ws = test.workspace / "mp7"
+
+    # model names
+    gwfname = get_model_name(name, "gwf")
+    prtname = get_model_name(name, "prt")
+    mp7name = get_model_name(name, "mp7")
+
+    # extract mf6 simulations/models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
-
-    # extract model grid
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(name, mp7_ws, test.targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -230,17 +240,17 @@ def test_mf6model(name, function_tmpdir, targets):
     gwf_head_file = f"{gwfname}.hds"
     prt_track_file = f"{prtname}.trk"
     prt_track_csv_file = f"{prtname}.trk.csv"
-    assert (ws / gwf_budget_file).is_file()
-    assert (ws / gwf_head_file).is_file()
-    assert (ws / prt_track_file).is_file()
-    assert (ws / prt_track_csv_file).is_file()
+    assert (gwf_ws / gwf_budget_file).is_file()
+    assert (gwf_ws / gwf_head_file).is_file()
+    assert (prt_ws / prt_track_file).is_file()
+    assert (prt_ws / prt_track_csv_file).is_file()
 
     # check mp7 output files exist
     mp7_pathline_file = f"{mp7name}.mppth"
-    assert (ws / mp7_pathline_file).is_file()
+    assert (mp7_ws / mp7_pathline_file).is_file()
 
     # load mp7 pathline results
-    plf = PathlineFile(ws / mp7_pathline_file)
+    plf = PathlineFile(mp7_ws / mp7_pathline_file)
     mp7_pls = pd.DataFrame(
         plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
     )
@@ -250,20 +260,22 @@ def test_mf6model(name, function_tmpdir, targets):
     mp7_pls["k"] = mp7_pls["k"] + 1
 
     # load mf6 pathline results
-    mf6_pls = pd.read_csv(ws / prt_track_csv_file)
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file)
 
     # check budget data were written to mf6 prt list file
-    check_budget_data(ws / f"{name}_prt.lst", ctx.perlen, ctx.nper)
+    check_budget_data(
+        prt_ws / f"{name}_prt.lst", BasicDisCase.perlen, BasicDisCase.nper
+    )
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
-        track_bin=ws / prt_track_file,
-        track_hdr=ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
-        track_csv=ws / prt_track_csv_file,
+        track_bin=prt_ws / prt_track_file,
+        track_hdr=prt_ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
+        track_csv=prt_ws / prt_track_csv_file,
     )
 
     # get head, budget, and spdis results from GWF model
-    hds = HeadFile(ws / gwf_head_file).get_data()
+    hds = HeadFile(gwf_ws / gwf_head_file).get_data()
     bud = gwf.output.budget()
     spdis = bud.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
@@ -351,7 +363,7 @@ def plot_stop_zone(nn, ax):
 
     # view/save plot
     # plt.show()
-    plt.savefig(ws / f"test_{name}_map.png")
+    plt.savefig(gwf_ws / f"test_{name}_map.png")
 
     # check that cell numbers are correct
     for i, row in list(mf6_pls.iterrows()):
@@ -400,3 +412,16 @@ def plot_stop_zone(nn, ax):
     # compare mf6 / mp7 pathline data
     assert mf6_pls.shape == mp7_pls.shape
     assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi04.py b/autotest/prt/test_prt_fmi04.py
index df5eab9ce36..b72573614f2 100644
--- a/autotest/prt/test_prt_fmi04.py
+++ b/autotest/prt/test_prt_fmi04.py
@@ -35,37 +35,44 @@
 from flopy.plot.plotutil import to_mp7_pathlines
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
-from prt_test_utils import (check_budget_data, check_track_data, get_gwf_sim,
-                            get_ireason_code, get_model_name)
+from prt_test_utils import (
+    BasicDisCase,
+    check_budget_data,
+    check_track_data,
+    get_ireason_code,
+    get_model_name,
+)
+
+from framework import TestFramework
 
 simname = "prtfmi04"
-ex = [simname, f"{simname}saws"]
+cases = [simname, f"{simname}saws"]
 
 
-def build_prt_sim(ctx, name, ws, mf6):
-    # output file names
+def build_prt_sim(name, gwf_ws, prt_ws, mf6):
+    # output files
     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"
+    gwf_budget_file = gwf_ws / f"{gwfname}.bud"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
+    prt_track_file = prt_ws / f"{prtname}.trk"
+    prt_track_csv_file = prt_ws / f"{prtname}.trk.csv"
 
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
-    pd = (ctx.perlen, ctx.nstp, ctx.tsmult)
+    pd = (BasicDisCase.perlen, BasicDisCase.nstp, BasicDisCase.tsmult)
     flopy.mf6.modflow.mftdis.ModflowTdis(
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=ctx.nper,
+        nper=BasicDisCase.nper,
         perioddata=[pd],
     )
 
@@ -76,21 +83,21 @@ def build_prt_sim(ctx, name, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=ctx.nlay,
-        nrow=ctx.nrow,
-        ncol=ctx.ncol,
+        nlay=BasicDisCase.nlay,
+        nrow=BasicDisCase.nrow,
+        ncol=BasicDisCase.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=BasicDisCase.porosity)
 
     # create prp package
     flopy.mf6.ModflowPrtprp(
         prt,
         pname="prp1",
         filename=f"{prtname}_1.prp",
-        nreleasepts=len(ctx.releasepts_prt),
-        packagedata=ctx.releasepts_prt,
+        nreleasepts=len(BasicDisCase.releasepts_prt),
+        packagedata=BasicDisCase.releasepts_prt,
         perioddata={0: ["FIRST"]},
         stop_at_weak_sink="saws" in name,
     )
@@ -123,15 +130,13 @@ def build_prt_sim(ctx, name, ws, mf6):
     return sim
 
 
-def build_mp7_sim(ctx, name, ws, mp7, gwf):
+def build_mp7_sim(name, ws, mp7, gwf):
     mp7name = f"{name}_mp7"
-    mp7_pathline_file = f"{mp7name}.mppth"
-
     partdata = flopy.modpath.ParticleData(
-        partlocs=[p[0] for p in ctx.releasepts_mp7],
-        localx=[p[1] for p in ctx.releasepts_mp7],
-        localy=[p[2] for p in ctx.releasepts_mp7],
-        localz=[p[3] for p in ctx.releasepts_mp7],
+        partlocs=[p[0] for p in BasicDisCase.releasepts_mp7],
+        localx=[p[1] for p in BasicDisCase.releasepts_mp7],
+        localy=[p[2] for p in BasicDisCase.releasepts_mp7],
+        localz=[p[3] for p in BasicDisCase.releasepts_mp7],
         timeoffset=0,
         drape=0,
     )
@@ -148,7 +153,7 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=ctx.porosity,
+        porosity=BasicDisCase.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -170,19 +175,11 @@ def get_different_rows(source_df, new_df):
     return changed_rows_df.drop("_merge", axis=1)
 
 
-@pytest.mark.parametrize("name", ex)
-def test_mf6model(name, function_tmpdir, targets):
-    # workspace
-    ws = function_tmpdir
-
-    # output file names
-    gwfname = get_model_name(name, "gwf")
-    prtname = get_model_name(name, "prt")
-    mp7name = get_model_name(name, "mp7")
-
+def build_models(idx, test):
     # build gwf model
-    gwfsim, ctx = get_gwf_sim(name, ws, targets.mf6)
-
+    gwfsim = BasicDisCase.get_gwf_sim(
+        test.name, test.workspace, test.targets.mf6
+    )
     # add wel package
     gwf = gwfsim.get_model()
     wells = [
@@ -197,21 +194,32 @@ def test_mf6model(name, function_tmpdir, targets):
     )
 
     # build prt model
-    prtsim = build_prt_sim(ctx, name, ws, targets.mf6)
+    prtsim = build_prt_sim(
+        test.name, test.workspace, test.workspace / "prt", test.targets.mf6
+    )
+    return gwfsim, prtsim
 
-    # run mf6 models
-    for sim in [gwfsim, prtsim]:
-        sim.write_simulation()
-        success, buff = sim.run_simulation(report=True)
-        assert success, pformat(buff)
+
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    mp7_ws = test.workspace / "mp7"
+
+    # model names
+    gwfname = get_model_name(name, "gwf")
+    prtname = get_model_name(name, "prt")
+    mp7name = get_model_name(name, "mp7")
 
     # extract mf6 models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(name, mp7_ws, test.targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -224,16 +232,16 @@ def test_mf6model(name, function_tmpdir, targets):
     prt_track_file = f"{prtname}.trk"
     prt_track_csv_file = f"{prtname}.trk.csv"
     mp7_pathline_file = f"{mp7name}.mppth"
-    assert (ws / gwf_budget_file).is_file()
-    assert (ws / gwf_head_file).is_file()
-    assert (ws / prt_track_file).is_file()
-    assert (ws / prt_track_csv_file).is_file()
+    assert (gwf_ws / gwf_budget_file).is_file()
+    assert (gwf_ws / gwf_head_file).is_file()
+    assert (prt_ws / prt_track_file).is_file()
+    assert (prt_ws / prt_track_csv_file).is_file()
 
     # check mp7 output files exist
-    assert (ws / mp7_pathline_file).is_file()
+    assert (mp7_ws / mp7_pathline_file).is_file()
 
     # load mp7 pathline results
-    plf = PathlineFile(ws / mp7_pathline_file)
+    plf = PathlineFile(mp7_ws / mp7_pathline_file)
     mp7_pls = pd.DataFrame(
         plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
     )
@@ -243,7 +251,7 @@ def test_mf6model(name, function_tmpdir, targets):
     mp7_pls["k"] = mp7_pls["k"] + 1
 
     # load mf6 pathline results
-    mf6_pls = pd.read_csv(ws / prt_track_csv_file)
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file)
 
     # if STOP_AT_WEAK_SINK disabled, check for an extra datum when particle exited weak sink
     wksk_irsn = get_ireason_code("WEAKSINK")
@@ -262,17 +270,19 @@ 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"{name}_prt.lst", ctx.perlen, ctx.nper)
+    check_budget_data(
+        prt_ws / f"{name}_prt.lst", BasicDisCase.perlen, BasicDisCase.nper
+    )
 
     # check mf6 prt particle track data were written to binary/CSV files
     check_track_data(
-        track_bin=ws / prt_track_file,
-        track_hdr=ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
-        track_csv=ws / prt_track_csv_file,
+        track_bin=prt_ws / prt_track_file,
+        track_hdr=prt_ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
+        track_csv=prt_ws / prt_track_csv_file,
     )
 
     # extract head, budget, and specific discharge results from GWF model
-    hds = HeadFile(ws / gwf_head_file).get_data()
+    hds = HeadFile(gwf_ws / gwf_head_file).get_data()
     bud = gwf.output.budget()
     spdis = bud.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
@@ -331,7 +341,7 @@ def all_equal(col, val):
 
     # view/save plot
     # plt.show()
-    plt.savefig(ws / f"test_{simname}.png")
+    plt.savefig(gwf_ws / f"test_{simname}.png")
 
     # convert mf6 pathlines to mp7 format
     mf6_pls = to_mp7_pathlines(mf6_pls)
@@ -359,3 +369,16 @@ def all_equal(col, val):
     # compare mf6 / mp7 pathline data
     assert mf6_pls.shape == mp7_pls.shape
     assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi05.py b/autotest/prt/test_prt_fmi05.py
index fa811407103..e7a58212791 100644
--- a/autotest/prt/test_prt_fmi05.py
+++ b/autotest/prt/test_prt_fmi05.py
@@ -38,14 +38,16 @@
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
 from prt_test_utils import (
+    BasicDisCase,
     all_equal,
     check_budget_data,
     check_track_data,
-    get_gwf_sim,
     get_model_name,
     get_partdata,
 )
 
+from framework import TestFramework
+
 simname = "prtfmi05"
 cases = [
     # options block options
@@ -76,13 +78,13 @@ def get_perioddata(name, periods=1, fraction=None) -> Optional[dict]:
     return {i: opt for i in range(periods)}
 
 
-def build_prt_sim(ctx, name, ws, mf6, fraction=None):
+def build_prt_sim(name, gwf_ws, prt_ws, mf6, fraction=None):
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
@@ -90,8 +92,10 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=ctx.nper,
-        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
+        nper=BasicDisCase.nper,
+        perioddata=[
+            (BasicDisCase.perlen, BasicDisCase.nstp, BasicDisCase.tsmult)
+        ],
     )
 
     # create prt model
@@ -102,20 +106,20 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=ctx.nlay,
-        nrow=ctx.nrow,
-        ncol=ctx.ncol,
+        nlay=BasicDisCase.nlay,
+        nrow=BasicDisCase.nrow,
+        ncol=BasicDisCase.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=BasicDisCase.porosity)
 
     # convert mp7 particledata to prt release points
-    partdata = get_partdata(prt.modelgrid, ctx.releasepts_mp7)
+    partdata = get_partdata(prt.modelgrid, BasicDisCase.releasepts_mp7)
     releasepts = list(partdata.to_prp(prt.modelgrid))
 
     # check release points match expectation
-    assert np.allclose(ctx.releasepts_prt, releasepts)
+    assert np.allclose(BasicDisCase.releasepts_prt, releasepts)
 
     # create prp package
     prp_track_file = f"{prtname}.prp.trk"
@@ -147,8 +151,8 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
 
     # create the flow model interface
     gwfname = get_model_name(name, "gwf")
-    gwf_budget_file = f"{gwfname}.bud"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.bud"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -168,11 +172,8 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
     return sim
 
 
-def build_mp7_sim(ctx, name, ws, mp7, gwf):
-    # convert mp7 particledata to prt release points
-    partdata = get_partdata(gwf.modelgrid, ctx.releasepts_mp7)
-
-    # create modpath 7 simulation
+def build_mp7_sim(name, ws, mp7, gwf):
+    partdata = get_partdata(gwf.modelgrid, BasicDisCase.releasepts_mp7)
     mp7name = get_model_name(name, "mp7")
     pg = flopy.modpath.ParticleGroup(
         particlegroupname="G1",
@@ -187,7 +188,7 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     )
     mpbas = flopy.modpath.Modpath7Bas(
         mp,
-        porosity=ctx.porosity,
+        porosity=BasicDisCase.porosity,
     )
     mpsim = flopy.modpath.Modpath7Sim(
         mp,
@@ -201,28 +202,33 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
     return mp
 
 
-@pytest.mark.parametrize("name", cases)
-@pytest.mark.parametrize("fraction", [0.5])
-def test_mf6model(name, function_tmpdir, targets, fraction):
-    # workspace
-    ws = function_tmpdir
+def build_models(idx, test, fraction):
+    # build mf6 models
+    gwfsim = BasicDisCase.get_gwf_sim(
+        test.name, test.workspace, test.targets.mf6
+    )
+    prtsim = build_prt_sim(
+        test.name,
+        test.workspace,
+        test.workspace / "prt",
+        test.targets.mf6,
+        fraction,
+    )
+    return gwfsim, prtsim
 
-    # model names
+
+def check_output(idx, test, fraction):
+    name = test.name
+    ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    mp7_ws = test.workspace / "mp7"
     gwfname = get_model_name(name, "gwf")
     prtname = get_model_name(name, "prt")
     mp7name = get_model_name(name, "mp7")
 
-    # build mf6 models
-    gwfsim, ctx = get_gwf_sim(name, ws, targets.mf6)
-    prtsim = build_prt_sim(ctx, name, ws, targets.mf6, fraction)
-
-    # run mf6 models
-    for sim in [gwfsim, prtsim]:
-        sim.write_simulation()
-        success, buff = sim.run_simulation(report=True)
-        assert success, pformat(buff)
-
-    # extract mf6 models
+    # extract mf6 simulations/models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
 
@@ -230,7 +236,7 @@ def test_mf6model(name, function_tmpdir, targets, fraction):
     mg = gwf.modelgrid
 
     # build mp7 model
-    mp7sim = build_mp7_sim(ctx, name, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(name, mp7_ws, test.targets.mp7, gwf)
 
     # run mp7 model
     mp7sim.write_input()
@@ -246,17 +252,17 @@ def test_mf6model(name, function_tmpdir, targets, fraction):
     prp_track_csv_file = f"{prtname}.prp.trk.csv"
     assert (ws / gwf_budget_file).is_file()
     assert (ws / gwf_head_file).is_file()
-    assert (ws / prt_track_file).is_file()
-    assert (ws / prt_track_csv_file).is_file()
-    assert (ws / prp_track_file).is_file()
-    assert (ws / prp_track_csv_file).is_file()
+    assert (prt_ws / prt_track_file).is_file()
+    assert (prt_ws / prt_track_csv_file).is_file()
+    assert (prt_ws / prp_track_file).is_file()
+    assert (prt_ws / prp_track_csv_file).is_file()
 
     # check mp7 output files exist
     mp7_pathline_file = f"{mp7name}.mppth"
-    assert (ws / mp7_pathline_file).is_file()
+    assert (mp7_ws / mp7_pathline_file).is_file()
 
     # load mp7 pathline results
-    plf = PathlineFile(ws / mp7_pathline_file)
+    plf = PathlineFile(mp7_ws / mp7_pathline_file)
     mp7_pls = pd.DataFrame(
         plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
     )
@@ -269,7 +275,7 @@ def test_mf6model(name, function_tmpdir, targets, fraction):
     mp7_pls["time"] = mp7_pls["time"] + fraction
 
     # load mf6 pathline results
-    mf6_pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False)
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
 
     # make sure pathline df has "name" (boundname) column and empty values
     assert "name" in mf6_pls
@@ -280,14 +286,20 @@ def test_mf6model(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", ctx.perlen, ctx.nper)
+    check_budget_data(
+        prt_ws / f"{name}_prt.lst", BasicDisCase.perlen, BasicDisCase.nper
+    )
 
     # check mf6 prt particle track data were written to binary/CSV files
     # and that different formats are equal
-    for track_csv in [ws / prt_track_csv_file, ws / prp_track_csv_file]:
+    for track_csv in [
+        prt_ws / prt_track_csv_file,
+        prt_ws / prp_track_csv_file,
+    ]:
         check_track_data(
-            track_bin=ws / prt_track_file,
-            track_hdr=ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
+            track_bin=prt_ws / prt_track_file,
+            track_hdr=prt_ws
+            / Path(prt_track_file.replace(".trk", ".trk.hdr")),
             track_csv=track_csv,
         )
 
@@ -362,3 +374,17 @@ def test_mf6model(name, function_tmpdir, targets, fraction):
     # compare mf6 / mp7 pathline data
     assert mf6_pls.shape == mp7_pls.shape
     assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+@pytest.mark.parametrize("fraction", [0.5])
+def test_mf6model(idx, name, function_tmpdir, targets, fraction):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t, fraction),
+        check=lambda t: check_output(idx, t, fraction),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi06.py b/autotest/prt/test_prt_fmi06.py
index 2bc117f1148..e7bb21c7385 100644
--- a/autotest/prt/test_prt_fmi06.py
+++ b/autotest/prt/test_prt_fmi06.py
@@ -28,8 +28,10 @@
     plot_nodes_and_vertices,
 )
 
+from framework import TestFramework
+
 simname = "prtfmi06"
-ex = [f"{simname}", f"{simname}bprp"]
+cases = [f"{simname}", f"{simname}bprp"]
 
 # model info
 nlay = 1
@@ -69,12 +71,8 @@
 ]
 
 
-def build_gwf_sim(idx, dir, mf6):
-    # model name
-    gwfname = f"{ex[idx]}_gwf"
-
-    # build MODFLOW 6 files
-    ws = dir
+def build_gwf_sim(idx, ws, mf6):
+    gwfname = f"{cases[idx]}_gwf"
     sim = flopy.mf6.MFSimulation(
         sim_name=gwfname, version="mf6", exe_name=mf6, sim_ws=ws
     )
@@ -157,14 +155,14 @@ def build_gwf_sim(idx, dir, mf6):
     return sim
 
 
-def build_prt_sim(idx, ws, mf6):
+def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
     # create simulation
-    name = ex[idx]
+    name = cases[idx]
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
@@ -173,7 +171,7 @@ def build_prt_sim(idx, ws, mf6):
     )
 
     # create prt model
-    prtname = f"{ex[idx]}_prt"
+    prtname = f"{cases[idx]}_prt"
     prt = flopy.mf6.ModflowPrt(sim, modelname=prtname)
 
     # create prt discretization
@@ -216,9 +214,9 @@ def build_prt_sim(idx, ws, mf6):
     )
 
     # create the flow model interface
-    gwfname = f"{ex[idx]}_gwf"
-    gwf_budget_file = f"{gwfname}.cbc"
-    gwf_head_file = f"{gwfname}.hds"
+    gwfname = f"{cases[idx]}_gwf"
+    gwf_budget_file = gwf_ws / f"{gwfname}.cbc"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -239,11 +237,8 @@ def build_prt_sim(idx, ws, mf6):
 
 
 def build_mp7_sim(idx, ws, mp7, gwf):
-    # convert mp7 particledata to prt release points
     partdata = get_partdata(gwf.modelgrid, releasepts_mp7)
-
-    # create modpath 7 simulation
-    mp7name = f"{ex[idx]}_mp7"
+    mp7name = f"{cases[idx]}_mp7"
     pg = flopy.modpath.ParticleGroup(
         particlegroupname="G1",
         particledata=partdata,
@@ -271,47 +266,38 @@ 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):
-    # workspace
-    ws = function_tmpdir
-
-    # test case name
-    name = ex[idx]
-
-    # model names
-    gwfname = f"{ex[idx]}_gwf"
-    prtname = f"{ex[idx]}_prt"
-    mp7name = f"{ex[idx]}_mp7"
+def build_models(idx, test):
+    gwfsim = build_gwf_sim(idx, test.workspace, test.targets.mf6)
+    prtsim = build_prt_sim(
+        idx, test.workspace, test.workspace / "prt", test.targets.mf6
+    )
+    return gwfsim, prtsim
 
-    # build mf6 models
-    gwfsim = build_gwf_sim(idx, ws, targets.mf6)
-    prtsim = build_prt_sim(idx, ws, targets.mf6)
 
-    # run gwf model
-    gwfsim.write_simulation()
-    success, buff = gwfsim.run_simulation(report=True)
-    assert success, pformat(buff)
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    mp7_ws = test.workspace / "mp7"
+    gwfname = f"{name}_gwf"
+    prtname = f"{name}_prt"
+    mp7name = f"{name}_mp7"
 
-    # run prt model
-    prtsim.write_simulation()
-    success, buff = prtsim.run_simulation(report=True)
+    # if invalid release points, check for error message
     if "bprp" in name:
-        assert not success, pformat(buff)
+        buff = test.buffs[1]
         assert any("Error: release point" in l for l in buff)
         return
-    else:
-        assert success, pformat(buff)
 
-    # extract mf6 models
+    # extract mf6 simulations/models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
-
-    # extract model grid
     mg = gwf.modelgrid
 
     # todo build mp7 model
-    mp7sim = build_mp7_sim(idx, ws, targets.mp7, gwf)
+    mp7sim = build_mp7_sim(idx, mp7_ws, test.targets.mp7, gwf)
 
     # todo run mp7 model
     mp7sim.write_input()
@@ -325,19 +311,19 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     prt_track_csv_file = f"{prtname}.trk.csv"
     prp_track_file = f"{prtname}.prp.trk"
     prp_track_csv_file = f"{prtname}.prp.trk.csv"
-    assert (ws / gwf_budget_file).is_file()
-    assert (ws / gwf_head_file).is_file()
-    assert (ws / prt_track_file).is_file()
-    assert (ws / prt_track_csv_file).is_file()
-    assert (ws / prp_track_file).is_file()
-    assert (ws / prp_track_csv_file).is_file()
+    assert (gwf_ws / gwf_budget_file).is_file()
+    assert (gwf_ws / gwf_head_file).is_file()
+    assert (prt_ws / prt_track_file).is_file()
+    assert (prt_ws / prt_track_csv_file).is_file()
+    assert (prt_ws / prp_track_file).is_file()
+    assert (prt_ws / prp_track_csv_file).is_file()
 
     # check mp7 output files exist
     mp7_pathline_file = f"{mp7name}.mppth"
-    assert (ws / mp7_pathline_file).is_file()
+    assert (mp7_ws / mp7_pathline_file).is_file()
 
     # load mp7 pathline results
-    plf = PathlineFile(ws / mp7_pathline_file)
+    plf = PathlineFile(mp7_ws / mp7_pathline_file)
     mp7_pls = pd.DataFrame(
         plf.get_destination_pathline_data(range(mg.nnodes), to_recarray=True)
     )
@@ -347,7 +333,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     mp7_pls["k"] = mp7_pls["k"] + 1
 
     # load mf6 pathline results
-    mf6_pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False)
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
 
     # make sure pathline df has "name" (boundname) column and default values
     assert "name" in mf6_pls
@@ -358,13 +344,13 @@ def test_mf6model(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, nstp)
+    check_budget_data(prt_ws / f"{name}_prt.lst", perlen, nper, nstp)
 
     # check mf6 prt particle track data were written to binary/CSV files
     # and that different formats are equal
     for track_bin, track_csv in zip(
-        [ws / prt_track_file, ws / prp_track_file],
-        [ws / prt_track_csv_file, ws / prp_track_csv_file],
+        [prt_ws / prt_track_file, prt_ws / prp_track_file],
+        [prt_ws / prt_track_csv_file, prt_ws / prp_track_csv_file],
     ):
         check_track_data(
             track_bin=track_bin,
@@ -373,7 +359,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
         )
 
     # extract head, budget, and specific discharge results from GWF model
-    hds = HeadFile(ws / gwf_head_file).get_data()
+    hds = HeadFile(gwf_ws / gwf_head_file).get_data()
     bud = gwf.output.budget()
     spdis = bud.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
@@ -423,7 +409,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
 
     # view/save plot
     # plt.show()
-    plt.savefig(ws / f"test_{simname}.png")
+    plt.savefig(gwf_ws / f"test_{simname}.png")
 
     # convert mf6 pathlines to mp7 format
     mf6_pls = to_mp7_pathlines(mf6_pls)
@@ -447,7 +433,19 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     del mp7_pls["node"]
 
     # compare mf6 / mp7 pathline data
-    # import pdb
-    # pdb.set_trace()
     assert mf6_pls.shape == mp7_pls.shape
     assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+        xfail=[False, "bprp" in name],
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi07.py b/autotest/prt/test_prt_fmi07.py
index 3d2ac82dbc4..ab332e62475 100644
--- a/autotest/prt/test_prt_fmi07.py
+++ b/autotest/prt/test_prt_fmi07.py
@@ -15,19 +15,21 @@
 
 import flopy
 import pytest
-from prt_test_utils import get_gwf_sim, get_model_name, get_partdata
+from prt_test_utils import BasicDisCase, get_model_name, get_partdata
+
+from framework import TestFramework
 
 simname = "prtfmi07"
-ex = [simname]
+cases = [simname]
 
 
-def build_prt_sim(ctx, name, ws, mf6):
+def build_prt_sim(name, gwf_ws, prt_ws, mf6):
     # create simulation
     sim = flopy.mf6.MFSimulation(
         sim_name=name,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
@@ -35,8 +37,10 @@ def build_prt_sim(ctx, name, ws, mf6):
         sim,
         pname="tdis",
         time_units="DAYS",
-        nper=ctx.nper,
-        perioddata=[(ctx.perlen, ctx.nstp, ctx.tsmult)],
+        nper=BasicDisCase.nper,
+        perioddata=[
+            (BasicDisCase.perlen, BasicDisCase.nstp, BasicDisCase.tsmult)
+        ],
     )
 
     # create prt model
@@ -47,16 +51,16 @@ def build_prt_sim(ctx, name, ws, mf6):
     flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
         prt,
         pname="dis",
-        nlay=ctx.nlay,
-        nrow=ctx.nrow,
-        ncol=ctx.ncol,
+        nlay=BasicDisCase.nlay,
+        nrow=BasicDisCase.nrow,
+        ncol=BasicDisCase.ncol,
     )
 
     # create mip package
-    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=ctx.porosity)
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=BasicDisCase.porosity)
 
     # convert mp7 to prt release points and check against expectation
-    partdata = get_partdata(prt.modelgrid, ctx.releasepts_mp7)
+    partdata = get_partdata(prt.modelgrid, BasicDisCase.releasepts_mp7)
     coords = partdata.to_coords(prt.modelgrid)
     # bad cell indices!
     releasepts = [(i, 0, 1, 1, c[0], c[1], c[2]) for i, c in enumerate(coords)]
@@ -89,8 +93,8 @@ def build_prt_sim(ctx, name, ws, mf6):
 
     # create the flow model interface
     gwfname = get_model_name(name, "gwf")
-    gwf_budget_file = f"{gwfname}.bud"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.bud"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -110,22 +114,31 @@ def build_prt_sim(ctx, name, ws, mf6):
     return sim
 
 
-@pytest.mark.parametrize("name", ex)
-def test_mf6model(name, function_tmpdir, targets):
-    # workspace
-    ws = function_tmpdir
-
+def build_models(idx, test):
     # build mf6 models
-    gwfsim, ctx = get_gwf_sim(name, ws, targets.mf6)
-    prtsim = build_prt_sim(ctx, name, ws, targets.mf6)
-
-    # run gwf models
-    gwfsim.write_simulation()
-    success, buff = gwfsim.run_simulation(report=True)
-    assert success, pformat(buff)
-
-    # run prt model (expect failure)
-    prtsim.write_simulation()
-    success, buff = prtsim.run_simulation(report=True)
-    assert not success, pformat(buff)
+    gwfsim = BasicDisCase.get_gwf_sim(
+        test.name, test.workspace, test.targets.mf6
+    )
+    prtsim = build_prt_sim(
+        test.name, test.workspace, test.workspace / "prt", test.targets.mf6
+    )
+    return gwfsim, prtsim
+
+
+def check_output(idx, test):
+    buff = test.buffs[1]
     assert any("Error: release point" in l for l in buff)
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+        xfail=[False, True],
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi08.py b/autotest/prt/test_prt_fmi08.py
index 67efe573a53..ff12c7fefc6 100644
--- a/autotest/prt/test_prt_fmi08.py
+++ b/autotest/prt/test_prt_fmi08.py
@@ -19,7 +19,6 @@
 
 
 from pathlib import Path
-from pprint import pformat
 
 import flopy
 import matplotlib.cm as cm
@@ -34,8 +33,10 @@
 from prt_test_utils import get_model_name
 from shapely.geometry import LineString, Point
 
+from framework import TestFramework
+
 simname = "prtfmi08"
-ex = [simname, f"{simname}_wel"]
+cases = [simname, f"{simname}_wel"]
 xmin = 0.0
 xmax = 2000.0
 ymin = 0.0
@@ -62,7 +63,7 @@
 
 
 def get_grid(workspace, targets):
-    workspace.mkdir(exist_ok=True)
+    workspace.mkdir(exist_ok=True, parents=True)
     tri = Triangle(
         maximum_area=area_max,
         angle=angle_min,
@@ -165,13 +166,13 @@ def build_gwf_sim(name, ws, targets):
     return sim
 
 
-def build_prt_sim(name, ws, targets):
-    ws = Path(ws)
+def build_prt_sim(name, gwf_ws, prt_ws, targets):
+    prt_ws = Path(prt_ws)
     gwfname = get_model_name(name, "gwf")
     prtname = get_model_name(name, "prt")
 
     # create grid
-    grid = get_grid(ws / "grid", targets)
+    grid = get_grid(prt_ws / "grid", targets)
     gridprops = grid.get_gridprops_vertexgrid()
     vgrid = VertexGrid(**gridprops, nlay=1)
     ibd = np.zeros(vgrid.ncpl, dtype=int)
@@ -195,7 +196,7 @@ def build_prt_sim(name, ws, targets):
 
     # create simulation
     sim = flopy.mf6.MFSimulation(
-        sim_name=name, version="mf6", exe_name=targets.mf6, sim_ws=ws
+        sim_name=name, version="mf6", exe_name=targets.mf6, sim_ws=prt_ws
     )
     tdis = flopy.mf6.ModflowTdis(
         sim, time_units="DAYS", perioddata=[[1.0, 1, 1.0]]
@@ -233,8 +234,8 @@ def build_prt_sim(name, ws, targets):
         track_filerecord=[prt_track_file],
         trackcsv_filerecord=[prt_track_csv_file],
     )
-    gwf_budget_file = f"{gwfname}.bud"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.bud"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -251,19 +252,25 @@ def build_prt_sim(name, ws, targets):
     return sim
 
 
-@pytest.mark.parametrize("name", ex)
-def test_mf6model(name, function_tmpdir, targets):
-    # workspace
-    ws = function_tmpdir
+def build_models(idx, test):
+    gwfsim = build_gwf_sim(test.name, test.workspace, test.targets)
+    prtsim = build_prt_sim(
+        test.name, test.workspace, test.workspace / "prt", test.targets
+    )
+    return gwfsim, prtsim
+
 
-    # build mf6 models
-    gwfsim = build_gwf_sim(name, ws, targets)
-    prtsim = build_prt_sim(name, ws, targets)
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    gwfname = get_model_name(name, "gwf")
+    prtname = get_model_name(name, "prt")
+    mp7name = get_model_name(name, "mp7")
 
-    # run gwf
-    gwfsim.write_simulation()
-    success, buff = gwfsim.run_simulation(report=True)
-    assert success, pformat(buff)
+    # extract mf6 simulations/models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
 
     # get gwf output
     gwf = gwfsim.get_model()
@@ -272,15 +279,9 @@ def test_mf6model(name, function_tmpdir, targets):
     spdis = bdobj.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
 
-    # run prt
-    prtsim.write_simulation()
-    success, buff = prtsim.run_simulation(report=True)
-    assert success, pformat(buff)
-
     # get prt output
-    prtname = get_model_name(name, "prt")
     prt_track_csv_file = f"{prtname}.prp.trk.csv"
-    pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False)
+    pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
 
     plot_debug = False
     if plot_debug:
@@ -307,9 +308,9 @@ def test_mf6model(name, function_tmpdir, targets):
         plt.show()
 
         # plot in 3d with pyvista (via vtk)
+        import pyvista as pv
         from flopy.export.vtk import Vtk
         from flopy.plot.plotutil import to_mp7_pathlines
-        import pyvista as pv
 
         def get_meshes(model, pathlines):
             vtk = Vtk(model=model, binary=False, smooth=False)
@@ -340,3 +341,16 @@ def callback(mesh, value):
         p.camera.zoom(1)
         p.add_slider_widget(lambda v: callback(path_mesh, v), [0, 30202])
         p.show()
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/prt/test_prt_fmi09.py b/autotest/prt/test_prt_fmi09.py
index 551993c0b20..af79251e800 100644
--- a/autotest/prt/test_prt_fmi09.py
+++ b/autotest/prt/test_prt_fmi09.py
@@ -11,7 +11,6 @@
 
 
 from pathlib import Path
-from pprint import pformat
 
 import flopy
 import matplotlib.cm as cm
@@ -22,8 +21,10 @@
 from flopy.utils.binaryfile import HeadFile
 from prt_test_utils import all_equal, check_track_data, get_model_name
 
+from framework import TestFramework
+
 simname = "prtfmi09"
-ex = [simname, f"{simname}_drp"]
+cases = [simname, f"{simname}_drp"]
 nlay, nrow, ncol = 2, 1, 5
 chdheads = [25.0]
 nper = len(chdheads)
@@ -136,8 +137,8 @@ def build_gwf_sim(name, ws, mf6):
     return sim
 
 
-def build_prt_sim(name, ws, mf6):
-    ws = Path(ws)
+def build_prt_sim(name, gwf_ws, prt_ws, mf6):
+    prt_ws = Path(prt_ws)
     gwfname = get_model_name(name, "gwf")
     prtname = get_model_name(name, "prt")
 
@@ -146,7 +147,7 @@ def build_prt_sim(name, ws, mf6):
         sim_name=prtname,
         exe_name=mf6,
         version="mf6",
-        sim_ws=ws,
+        sim_ws=prt_ws,
     )
 
     # create tdis package
@@ -198,8 +199,8 @@ def build_prt_sim(name, ws, mf6):
     )
 
     # create the flow model interface
-    gwf_budget_file = f"{gwfname}.cbc"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.cbc"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -219,29 +220,23 @@ def build_prt_sim(name, ws, mf6):
     return sim
 
 
-@pytest.mark.parametrize("name", ex)
-def test_mf6model(name, function_tmpdir, targets):
-    # workspace
-    ws = function_tmpdir
+def build_models(idx, test):
+    gwfsim = build_gwf_sim(test.name, test.workspace, test.targets.mf6)
+    prtsim = build_prt_sim(test.name, test.workspace, test.workspace / "prt", test.targets.mf6)
+    return gwfsim, prtsim
 
-    # determine if drape is enabled
-    drape = "drp" in name
 
-    # model names
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
     gwfname = get_model_name(name, "gwf")
     prtname = get_model_name(name, "prt")
+    drape = "drp" in name
 
-    # build mf6 models
-    gwfsim = build_gwf_sim(name, ws, targets.mf6)
-    prtsim = build_prt_sim(name, ws, targets.mf6)
-
-    # run mf6 models
-    for sim in [gwfsim, prtsim]:
-        sim.write_simulation()
-        success, buff = sim.run_simulation(report=True)
-        assert success, pformat(buff)
-
-    # extract mf6 models and grid
+    # extract mf6 simulations/models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
     gwf = gwfsim.get_model(gwfname)
     prt = prtsim.get_model(prtname)
     mg = gwf.modelgrid
@@ -253,15 +248,15 @@ def test_mf6model(name, function_tmpdir, targets):
     prt_track_csv_file = f"{prtname}.trk.csv"
     prp_track_file = f"{prtname}.prp.trk"
     prp_track_csv_file = f"{prtname}.prp.trk.csv"
-    assert (ws / gwf_budget_file).is_file()
-    assert (ws / gwf_head_file).is_file()
-    assert (ws / prt_track_file).is_file()
-    assert (ws / prt_track_csv_file).is_file()
-    assert (ws / prp_track_file).is_file()
-    assert (ws / prp_track_csv_file).is_file()
+    assert (gwf_ws / gwf_budget_file).is_file()
+    assert (gwf_ws / gwf_head_file).is_file()
+    assert (prt_ws / prt_track_file).is_file()
+    assert (prt_ws / prt_track_csv_file).is_file()
+    assert (prt_ws / prp_track_file).is_file()
+    assert (prt_ws / prp_track_csv_file).is_file()
 
     # load mf6 pathline results
-    mf6_pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False)
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
 
     # make sure all mf6 pathline data have correct model and PRP index (1)
     assert all_equal(mf6_pls["imdl"], 1)
@@ -272,15 +267,15 @@ def test_mf6model(name, function_tmpdir, targets):
 
     # check mf6 prt particle track data were written to binary/CSV files
     # and that different formats are equal
-    for track_csv in [ws / prt_track_csv_file, ws / prp_track_csv_file]:
+    for track_csv in [prt_ws / prt_track_csv_file, prt_ws / prp_track_csv_file]:
         check_track_data(
-            track_bin=ws / prt_track_file,
-            track_hdr=ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
+            track_bin=prt_ws / prt_track_file,
+            track_hdr=prt_ws / Path(prt_track_file.replace(".trk", ".trk.hdr")),
             track_csv=track_csv,
         )
 
     # extract head, budget, and specific discharge results from GWF model
-    hds = HeadFile(ws / gwf_head_file).get_data()
+    hds = HeadFile(gwf_ws / gwf_head_file).get_data()
     bud = gwf.output.budget()
     spdis = bud.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
@@ -308,7 +303,7 @@ def test_mf6model(name, function_tmpdir, targets):
 
     # view/save plot
     # plt.show()
-    plt.savefig(ws / f"test_{simname}.png")
+    plt.savefig(gwf_ws / f"test_{simname}.png")
 
     if drape:
         assert mf6_pls.shape[0] == 36
@@ -317,3 +312,16 @@ def test_mf6model(name, function_tmpdir, targets):
         assert mf6_pls.shape[0] == 9
         # istatus=8 permanently unreleased
         assert mf6_pls.istatus.eq(8).all()
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
\ No newline at end of file
diff --git a/autotest/prt/test_prt_fmi10.py b/autotest/prt/test_prt_fmi10.py
index df3bd50ff5a..57800f99eb7 100644
--- a/autotest/prt/test_prt_fmi10.py
+++ b/autotest/prt/test_prt_fmi10.py
@@ -13,7 +13,6 @@
 
 from math import isclose
 from pathlib import Path
-from pprint import pformat
 
 import flopy
 import matplotlib.pyplot as plt
@@ -26,6 +25,8 @@
 from prt_test_utils import get_model_name
 from shapely.geometry import LineString
 
+from framework import TestFramework
+
 simname = "prtfmi10"
 cases = [f"{simname}l2r", f"{simname}diag"]
 angle = 30
@@ -44,7 +45,7 @@ def chdhead(x):
 
 
 def get_tri(workspace, targets) -> Triangle:
-    workspace.mkdir(exist_ok=True)
+    workspace.mkdir(exist_ok=True, parents=True)
     tri = Triangle(
         angle=angle,
         maximum_area=max_area,
@@ -118,13 +119,13 @@ def build_gwf_sim(name, ws, targets):
     return sim
 
 
-def build_prt_sim(name, ws, targets):
-    ws = Path(ws)
+def build_prt_sim(name, gwf_ws, prt_ws, targets):
+    prt_ws = Path(prt_ws)
     gwfname = get_model_name(name, "gwf")
     prtname = get_model_name(name, "prt")
 
     # create grid
-    tri = get_tri(ws / "grid", targets)
+    tri = get_tri(prt_ws / "grid", targets)
     grid = VertexGrid(tri)
     gi = GridIntersect(grid)
 
@@ -135,7 +136,7 @@ def build_prt_sim(name, ws, targets):
 
     # create simulation
     sim = flopy.mf6.MFSimulation(
-        sim_name=name, version="mf6", exe_name=targets.mf6, sim_ws=ws
+        sim_name=name, version="mf6", exe_name=targets.mf6, sim_ws=prt_ws
     )
     tdis = flopy.mf6.ModflowTdis(
         sim, time_units="DAYS", perioddata=[[1.0, 1, 1.0]]
@@ -184,8 +185,8 @@ def build_prt_sim(name, ws, targets):
         track_filerecord=[prt_track_file],
         trackcsv_filerecord=[prt_track_csv_file],
     )
-    gwf_budget_file = f"{gwfname}.cbc"
-    gwf_head_file = f"{gwfname}.hds"
+    gwf_budget_file = gwf_ws / f"{gwfname}.cbc"
+    gwf_head_file = gwf_ws / f"{gwfname}.hds"
     flopy.mf6.ModflowPrtfmi(
         prt,
         packagedata=[
@@ -202,19 +203,28 @@ def build_prt_sim(name, ws, targets):
     return sim
 
 
-@pytest.mark.parametrize("name", cases)
-def test_mf6model(name, function_tmpdir, targets):
-    # workspace
-    ws = function_tmpdir
+def build_models(idx, test):
+    gwfsim = build_gwf_sim(test.name, test.workspace, test.targets)
+    prtsim = build_prt_sim(
+        test.name, test.workspace, test.workspace / "prt", test.targets
+    )
+    return gwfsim, prtsim
+
 
-    # build mf6 models
-    gwfsim = build_gwf_sim(name, ws, targets)
-    prtsim = build_prt_sim(name, ws, targets)
+def check_output(idx, test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    gwfname = get_model_name(name, "gwf")
+    prtname = get_model_name(name, "prt")
+    drape = "drp" in name
 
-    # run gwf
-    gwfsim.write_simulation()
-    success, buff = gwfsim.run_simulation(report=True)
-    assert success, pformat(buff)
+    # extract mf6 simulations/models and grid
+    gwfsim = test.sims[0]
+    prtsim = test.sims[1]
+    gwf = gwfsim.get_model(gwfname)
+    prt = prtsim.get_model(prtname)
+    mg = gwf.modelgrid
 
     # get gwf output
     gwf = gwfsim.get_model()
@@ -223,15 +233,10 @@ def test_mf6model(name, function_tmpdir, targets):
     spdis = bdobj.get_data(text="DATA-SPDIS")[0]
     qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
 
-    # run prt
-    prtsim.write_simulation()
-    success, buff = prtsim.run_simulation(report=True)
-    assert success, pformat(buff)
-
     # get prt output
     prtname = get_model_name(name, "prt")
     prt_track_csv_file = f"{prtname}.prp.trk.csv"
-    pls = pd.read_csv(ws / prt_track_csv_file, na_filter=False)
+    pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
     endpts = (
         pls.sort_values("t")
         .groupby(["imdl", "iprp", "irpt", "trelease"])
@@ -269,3 +274,16 @@ def test_mf6model(name, function_tmpdir, targets):
         assert pls.shape == (112, 16)
         assert endpts.shape == (2, 16)
         assert set(endpts.icell) == {111, 144}
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+        compare=None,
+    )
+    test.run()
diff --git a/autotest/test_gwf_utl03_obs01.py b/autotest/test_gwf_utl03_obs01.py
index e2badc32ec9..8141508cf7e 100644
--- a/autotest/test_gwf_utl03_obs01.py
+++ b/autotest/test_gwf_utl03_obs01.py
@@ -153,8 +153,8 @@ def build_model(idx, dir, exe):
     return sim, mc
 
 
-def build_models(idx, test, exe):
-    sim, mc = build_model(idx, test.workspace, exe)
+def build_models(idx, test):
+    sim, mc = build_model(idx, test.workspace, test.targets.mf6)
     sim.write_simulation()
     mc.write_simulation()
     hack_binary_obs(idx, test.workspace)
@@ -221,7 +221,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
-        build=lambda t: build_models(idx, t, targets.mf6),
+        build=lambda t: build_models(idx, t),
         check=lambda t: check_output(idx, t),
         targets=targets,
         overwrite=False,