From 2302b23b18be8ab580a318830b69cd72b99d247f Mon Sep 17 00:00:00 2001
From: Joseph Hughes <jdhughes99+1@gmail.com>
Date: Sat, 22 Aug 2026 23:22:53 -0500
Subject: [PATCH 01/10] feat(Mf6Splitter): split structured models into DISV
 models

split_model() and split_multi_model() take to_disv, which writes DISV models from a
structured model instead of the structured bounding box of each part. A DISV model
carries only the cells assigned to it, and stacks of cells that are inactive in every
layer are excluded, so a model no longer pads out to a bounding box full of inactive
cells. to_disv is ignored with a warning for a DISV model and is an error for a DISU
model.

The grid of a model is built one model at a time from the new
StructuredGrid.get_cell_iverts(), so peak memory scales with the largest part rather
than with the parent grid. save_node_mapping() records the grid type of the split
models so a structured to DISV mapping round trips.

Closes #2816
---
 autotest/test_grid.py                  |  18 ++
 autotest/test_model_splitter.py        | 277 +++++++++++++++++++++++++
 flopy/discretization/structuredgrid.py |  40 +++-
 flopy/mf6/utils/model_splitter.py      | 183 ++++++++++++++--
 4 files changed, 487 insertions(+), 31 deletions(-)

diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index 2cd4e9096..b3bcfee69 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -285,6 +285,24 @@ def test_unstructured_grid_get_cell_vertices():
     assert v1 == v4, "Positional and kwarg should match"
 
 
+def test_get_cell_iverts():
+    nrow, ncol = 4, 5
+    grid = StructuredGrid(delr=np.ones(ncol), delc=np.ones(nrow))
+
+    iverts = grid.get_cell_iverts()
+    assert iverts.shape == (grid.ncpl, 4)
+    assert np.array_equal(iverts.tolist(), grid.iverts)
+
+    nodes = [0, 7, grid.ncpl - 1]
+    assert np.array_equal(grid.get_cell_iverts(nodes), iverts[nodes])
+    assert np.array_equal(grid.get_cell_iverts(3), iverts[[3]])
+
+    # vertices are ordered clockwise from the upper left corner of a cell
+    for node, iv in enumerate(iverts):
+        i, j = divmod(node, ncol)
+        assert list(iv) == grid._build_structured_iverts(i, j)
+
+
 def test_get_lrc_get_node():
     nlay, nrow, ncol = 3, 4, 5
     nnodes = nlay * nrow * ncol
diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index 60bbb8330..f1d75b474 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -138,6 +138,283 @@ def test_unstructured_model_splitter(function_tmpdir):
     np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg)
 
 
+@requires_exe("mf6")
+def test_structured_to_disv_model_splitter(function_tmpdir):
+    sim_path = get_example_data_path() / "mf6-freyberg"
+    split_path = function_tmpdir / "split_model"
+
+    sim = MFSimulation.load(sim_ws=sim_path)
+    sim.set_sim_path(function_tmpdir)
+    sim.write_simulation()
+    sim.run_simulation()
+
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    idomain = modelgrid.idomain.reshape((modelgrid.nlay, modelgrid.ncpl))
+
+    array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int)
+    array[modelgrid.nrow // 2 :, :] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_model(array, to_disv=True)
+
+    new_sim.set_sim_path(split_path)
+    new_sim.write_simulation()
+    new_sim.run_simulation()
+
+    ml0 = new_sim.get_model("freyberg_0")
+    ml1 = new_sim.get_model("freyberg_1")
+
+    for ml in (ml0, ml1):
+        if ml.modelgrid.grid_type != "vertex":
+            raise AssertionError("Split model is not a DISV model")
+
+    # every stack of cells that is active in a layer, and no other, is split
+    nactive = np.count_nonzero(np.any(idomain != 0, axis=0))
+    if ml0.modelgrid.ncpl + ml1.modelgrid.ncpl != nactive:
+        raise AssertionError("Split models do not include all active cells")
+
+    original_heads = gwf.output.head().get_alldata()[-1]
+    heads0 = ml0.output.head().get_alldata()[-1]
+    heads1 = ml1.output.head().get_alldata()[-1]
+
+    new_heads = mfsplit.reconstruct_array({0: heads0, 1: heads1})
+
+    idx = modelgrid.idomain != 0
+    err_msg = "Heads from original and split models do not match"
+    np.testing.assert_allclose(new_heads[idx], original_heads[idx], err_msg=err_msg)
+
+
+@requires_exe("mf6")
+def test_structured_to_disv_idomain(function_tmpdir):
+    nlay, nrow, ncol = 3, 10, 10
+    idomain = np.ones((nlay, nrow, ncol), dtype=int)
+    idomain[:, 0:2, 0:2] = 0  # inactive in every layer, excluded from the split
+    idomain[0, 5:, 5:] = 0  # inactive in one layer, carried into the split
+    idomain[2, 0:3, 7:] = 0
+
+    sim = flopy.mf6.MFSimulation(sim_name="ml", sim_ws=function_tmpdir, exe_name="mf6")
+    flopy.mf6.ModflowTdis(sim)
+    flopy.mf6.ModflowIms(sim, complexity="simple")
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="ml", save_flows=True)
+    flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=100.0,
+        delc=100.0,
+        top=30.0,
+        botm=[20.0, 10.0, 0.0],
+        idomain=idomain,
+    )
+    flopy.mf6.ModflowGwfnpf(gwf, k=1.0)
+    flopy.mf6.ModflowGwfic(gwf, strt=25.0)
+    chd = [[(0, i, 2), 25.0] for i in range(2, nrow)]
+    chd += [[(0, i, ncol - 1), 20.0] for i in range(0, 5)]
+    flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd)
+    flopy.mf6.ModflowGwfoc(
+        gwf,
+        head_filerecord="ml.hds",
+        saverecord=[("HEAD", "ALL")],
+    )
+    sim.write_simulation()
+    sim.run_simulation()
+
+    original_heads = gwf.output.head().get_alldata()[-1]
+
+    # split on a raveled array that covers every cell of the DIS grid
+    array = np.zeros((nrow, ncol), dtype=int)
+    array[:, ncol // 2 :] = 1
+    array = array.ravel()
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_model(array, to_disv=True)
+    new_sim.set_sim_path(function_tmpdir / "split_model")
+    new_sim.write_simulation()
+    new_sim.run_simulation()
+
+    heads = {}
+    ncpl = 0
+    oidomain = idomain.reshape((nlay, nrow * ncol))
+    for mkey in (0, 1):
+        ml = new_sim.get_model(f"ml_{mkey}")
+        grid = ml.modelgrid
+        nodes = mfsplit._grid_info[mkey][-1]
+        ncpl += grid.ncpl
+
+        # a stack that is inactive in one layer keeps its idomain
+        ninactive = np.count_nonzero(oidomain[:, nodes] == 0)
+        if np.count_nonzero(grid.idomain == 0) != ninactive:
+            raise AssertionError(f"Model {mkey} idomain was not remapped")
+
+        heads[mkey] = ml.output.head().get_alldata()[-1]
+
+    # a stack that is inactive in every layer is excluded
+    nactive = np.count_nonzero(np.any(oidomain != 0, axis=0))
+    if ncpl != nactive:
+        raise AssertionError("Split models do not include all active cells")
+
+    new_heads = mfsplit.reconstruct_array(heads)
+
+    idx = idomain != 0
+    err_msg = "Heads from original and split models do not match"
+    np.testing.assert_allclose(new_heads[idx], original_heads[idx], err_msg=err_msg)
+
+
+def test_structured_to_disv_reconstruct_recarray():
+    sim_path = get_example_data_path() / "mf6-freyberg"
+
+    sim = MFSimulation.load(sim_ws=sim_path)
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+
+    array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int)
+    array[modelgrid.nrow // 2 :, :] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_model(array, to_disv=True)
+
+    # riv and wel are split across both models, chd is in one model
+    for pkgtype in ("riv", "wel", "chd"):
+        original = getattr(gwf, pkgtype).stress_period_data.get_data(0)
+        recarrays = {}
+        for mkey in (0, 1):
+            ml = new_sim.get_model(f"freyberg_{mkey}")
+            if hasattr(ml, pkgtype):
+                recarrays[mkey] = getattr(ml, pkgtype).stress_period_data.get_data(0)
+
+        new_recarray = mfsplit.reconstruct_recarray(recarrays)
+
+        # cellids are returned as the layer, row, column of the original model
+        onodes = modelgrid.get_node([tuple(cid) for cid in original.cellid])
+        nnodes = modelgrid.get_node([tuple(cid) for cid in new_recarray.cellid])
+        oidx = np.argsort(onodes)
+        nidx = np.argsort(nnodes)
+
+        err_msg = f"Reconstructed {pkgtype} recarray does not match the original"
+        if not np.array_equal(np.array(onodes)[oidx], np.array(nnodes)[nidx]):
+            raise AssertionError(err_msg)
+
+        for name in original.dtype.names:
+            if name == "cellid":
+                continue
+            np.testing.assert_allclose(
+                new_recarray[name][nidx], original[name][oidx], err_msg=err_msg
+            )
+
+
+def test_structured_to_disv_model_geometry():
+    sim_path = get_example_data_path() / "mf6-freyberg"
+
+    sim = MFSimulation.load(sim_ws=sim_path)
+    gwf = sim.get_model()
+    gwf.modelgrid.set_coord_info(xoff=1000.0, yoff=2000.0, angrot=15.0)
+    modelgrid = gwf.modelgrid
+
+    array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int)
+    array[:, modelgrid.ncol // 2 :] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    mfsplit.split_model(array, to_disv=True)
+
+    xcenters = modelgrid.xcellcenters.ravel()
+    ycenters = modelgrid.ycellcenters.ravel()
+    areas = modelgrid.area.ravel()
+    for mkey, model in mfsplit._model_dict.items():
+        grid = model.modelgrid
+        nodes = mfsplit._grid_info[mkey][-1]
+
+        err_msg = f"Model {mkey} cells are not coincident with original cells"
+        np.testing.assert_allclose(grid.xcellcenters, xcenters[nodes], err_msg=err_msg)
+        np.testing.assert_allclose(grid.ycellcenters, ycenters[nodes], err_msg=err_msg)
+
+        # MODFLOW 6 calculates a CELL2D area that is positive only when the
+        # vertices are listed in clockwise order (Disv.f90 get_cell2d_area)
+        cell2d = model.disv.cell2d.array
+        vertices = model.disv.vertices.array
+        icvert = np.column_stack([cell2d[f"icvert_{iv}"] for iv in range(4)])
+        x = vertices["xv"][icvert]
+        y = vertices["yv"][icvert]
+        cell_area = -0.5 * np.sum(
+            x * np.roll(y, -1, axis=1) - np.roll(x, -1, axis=1) * y, axis=1
+        )
+        err_msg = f"Model {mkey} CELL2D vertices are not in clockwise order"
+        np.testing.assert_allclose(cell_area, areas[nodes], err_msg=err_msg)
+
+
+def test_to_disv_grid_types():
+    sim_path = get_example_data_path() / "mf6" / "test003_gwftri_disv"
+    sim = MFSimulation.load(sim_ws=sim_path)
+    modelgrid = sim.get_model().modelgrid
+
+    array = np.zeros((modelgrid.ncpl,), dtype=int)
+    array[0:85] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    with pytest.warns(UserWarning, match="already a DISV model"):
+        mfsplit.split_model(array, to_disv=True)
+
+    sim_path = get_example_data_path() / "mf6" / "test006_gwf3"
+    sim = MFSimulation.load(sim_ws=sim_path)
+    modelgrid = sim.get_model().modelgrid
+
+    array = np.zeros((modelgrid.nnodes,), dtype=int)
+    array[65:] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    with pytest.raises(ValueError, match="DISU"):
+        mfsplit.split_model(array, to_disv=True)
+
+
+@requires_exe("mf6")
+@requires_pkg("h5py")
+def test_save_load_node_mapping_structured_to_disv(function_tmpdir):
+    sim_path = get_example_data_path() / "mf6-freyberg"
+    new_sim_path = function_tmpdir / "split_model"
+    hdf_file = new_sim_path / "node_map.hdf5"
+
+    sim = MFSimulation.load(sim_ws=sim_path)
+    sim.set_sim_path(function_tmpdir)
+    sim.write_simulation()
+    sim.run_simulation()
+
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+    original_heads = gwf.output.head().get_alldata()[-1]
+
+    array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int)
+    array[modelgrid.nrow // 2 :, :] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_model(array, to_disv=True)
+    new_sim.set_sim_path(new_sim_path)
+    new_sim.write_simulation()
+    new_sim.run_simulation()
+    original_node_map = mfsplit._node_map
+
+    mfsplit.save_node_mapping(hdf_file)
+
+    new_sim2 = MFSimulation.load(sim_ws=new_sim_path)
+    mfsplit2 = Mf6Splitter.load_node_mapping(hdf_file)
+
+    if original_node_map != mfsplit2._node_map:
+        raise AssertionError("Node map read/write not returning proper values")
+
+    array_dict = {}
+    for mkey in (0, 1):
+        ml = new_sim2.get_model(f"freyberg_{mkey}")
+        if ml.modelgrid.grid_type != "vertex":
+            raise AssertionError("Loaded node mapping grid type is not vertex")
+        array_dict[mkey] = ml.output.head().get_alldata()[-1]
+
+    new_heads = mfsplit2.reconstruct_array(array_dict)
+
+    idx = modelgrid.idomain != 0
+    err_msg = "Heads from original and split models do not match"
+    np.testing.assert_allclose(new_heads[idx], original_heads[idx], err_msg=err_msg)
+
+
 @requires_exe("mf6")
 @pytest.mark.slow
 def test_model_with_lak_sfr_mvr(function_tmpdir):
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index 690cf306e..a6561f3a2 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -2001,21 +2001,43 @@ def get_plottable_layer_array(self, a, layer):
         assert plotarray.shape == required_shape, msg
         return plotarray
 
+    def get_cell_iverts(self, nodes=None):
+        """
+        Get the vertex numbers that define one or more model cells
+
+        Parameters
+        ----------
+        nodes : int, list, np.ndarray
+            optional two-dimensional node numbers. Vertex numbers are returned
+            for every cell in a layer when nodes is None.
+
+        Returns
+        -------
+        np.ndarray : array of vertex numbers with shape (nnodes, 4). The
+            vertices of a cell are ordered clockwise from its upper left
+            corner.
+
+        """
+        if nodes is None:
+            nodes = np.arange(self.ncpl, dtype=int)
+        else:
+            nodes = np.atleast_1d(nodes)
+
+        i, j = np.divmod(nodes, self.ncol)
+        iverts = np.empty((4, nodes.size), dtype=int)
+        iverts[0] = i * (self.ncol + 1) + j
+        iverts[1] = iverts[0] + 1
+        iverts[2] = iverts[1] + self.ncol + 1
+        iverts[3] = iverts[0] + self.ncol + 1
+        return iverts.T
+
     def _set_structured_iverts(self):
         """
         Build a list of the vertices that define each model cell and the x, y
         pair for each vertex
 
         """
-        rowarr = np.repeat(np.arange(self.nrow, dtype=int), self.ncol)
-        colarr = np.tile(np.arange(self.ncol, dtype=int), self.nrow)
-
-        iverts = np.empty((4, self.ncpl), dtype=int)
-        iverts[0] = rowarr * (self.ncol + 1) + colarr
-        iverts[1] = rowarr * (self.ncol + 1) + colarr + 1
-        iverts[2] = (rowarr + 1) * (self.ncol + 1) + colarr + 1
-        iverts[3] = (rowarr + 1) * (self.ncol + 1) + colarr
-        self._iverts = iverts.T.tolist()
+        self._iverts = self.get_cell_iverts().tolist()
         return
 
     def _build_structured_iverts(self, i, j):
diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py
index 8b333cb5e..afe98d51f 100644
--- a/flopy/mf6/utils/model_splitter.py
+++ b/flopy/mf6/utils/model_splitter.py
@@ -1,4 +1,5 @@
 import inspect
+import warnings
 
 import numpy as np
 
@@ -153,6 +154,8 @@ def __init__(self, sim, modelname=None):
             self._ncpl = self._modelgrid.nnodes
         self._shape = self._modelgrid.shape
         self._grid_type = self._modelgrid.grid_type
+        self._new_grid_type = self._grid_type
+        self._to_disv = False
         self._node_map = {}
         self._node_map_r = {}
         self._new_connections = None
@@ -223,6 +226,8 @@ def switch_models(self, modelname, remap_nodes=False):
 
         if remap_nodes:
             self._modelgrid = self._model.modelgrid
+            self._grid_type = self._modelgrid.grid_type
+            self._new_grid_type = self._grid_type
             self._node_map = {}
             self._node_map_r = {}
             self._new_connections = None
@@ -283,10 +288,12 @@ def _map_grids_to_hdf5(self, f):
         }
 
         mgg = f.create_group("modelgrids")
-        grid_type = f["grid_type"][0].decode("utf8")
         mkeys = f["mkey"][:]
         mnames = [self._modelname] + list(f["new_modelnames"][:])
         grids = [self._modelgrid, ] + [self._model_dict[mk].modelgrid for mk in mkeys]
+        grid_types = [f["grid_type"][0].decode("utf8")] + [
+            f["new_grid_type"][0].decode("utf8")
+        ] * len(mkeys)
 
         for ix, name in enumerate(mnames):
             if hasattr(name, "decode"):
@@ -295,6 +302,7 @@ def _map_grids_to_hdf5(self, f):
             # name = name.decode("utf8")
             grd_grp = mgg.create_group(name)
             grid = grids[ix]
+            grid_type = grid_types[ix]
             if grid_type == "structured":
                 delc = grid.delc
                 dc = grd_grp.create_dataset("delc", (len(delc),), dtype=float, **kwargs)
@@ -411,6 +419,9 @@ def save_node_mapping(self, filename):
         gt_ds = f.create_dataset("grid_type", (1,), dtype=string_dt)
         gt_ds[:] = [self._grid_type,]
 
+        ngt_ds = f.create_dataset("new_grid_type", (1,), dtype=string_dt)
+        ngt_ds[:] = [self._new_grid_type,]
+
         mname_ds = f.create_dataset("modelname", (1,), dtype=string_dt)
         mname_ds[:] = [self._modelname,]
 
@@ -561,10 +572,15 @@ def construct_modelgrid(f, name, grid_type):
         # construct representation of original model geometry
         modelname = f["modelname"][0].decode("utf8")
         grid_type = f["grid_type"][0].decode("utf8")
+        new_grid_type = grid_type
+        if "new_grid_type" in f:
+            new_grid_type = f["new_grid_type"][0].decode("utf8")
+
         grid = construct_modelgrid(f, modelname, grid_type)
         ml = FakeModel(modelname, grid)
         sim = FakeSim({modelname: ml})
         mfs = Mf6Splitter(sim)
+        mfs._new_grid_type = new_grid_type
         mfs._model_dict = {}
         mfs._new_ncpl = {}
 
@@ -573,10 +589,10 @@ def construct_modelgrid(f, name, grid_type):
         modelnames = [i.decode("utf8") for i in f["new_modelnames"][:]]
         for ix, mname in enumerate(modelnames):
             mkey = mkeys[ix]
-            sgrid = construct_modelgrid(f, mname, grid_type)
+            sgrid = construct_modelgrid(f, mname, new_grid_type)
             sml = FakeModel(mname, sgrid)
             mfs._model_dict[mkey] = sml
-            if grid_type in ("structured", "vertex"):
+            if new_grid_type in ("structured", "vertex"):
                 mfs._new_ncpl[mkey] = sgrid.ncpl
             else:
                 mfs._new_ncpl[mkey] = sgrid.nnodes
@@ -584,8 +600,9 @@ def construct_modelgrid(f, name, grid_type):
         f.close()
 
         # create additional splitting data efficient reconstruction
-        split_array = np.zeros((mfs._ncpl,), dtype=int)
-        model_array = np.zeros((mfs._ncpl,), dtype=int)
+        # cells that were not remapped keep a model number of -1
+        split_array = np.full((mfs._ncpl,), -1, dtype=int)
+        model_array = np.full((mfs._ncpl,), -1, dtype=int)
         for k, v in node_map.items():
             k = int(k)
             model_array[k] = v[0]
@@ -867,9 +884,9 @@ def reconstruct_recarray(self, recarrays):
             remapper = self.reversed_node_map[mkey]
             orec = recarray.copy()
             modelgrid = self._model_dict[mkey].modelgrid
-            if self._grid_type in ("structured", "vertex"):
+            if self._new_grid_type in ("structured", "vertex"):
                 layer = [i[0] for i in orec.cellid]
-                if self._grid_type == "structured":
+                if self._new_grid_type == "structured":
                     cellid = [(0, i[1], i[2]) for i in orec.cellid]
                     node = modelgrid.get_node(cellid)
                 else:
@@ -879,7 +896,7 @@ def reconstruct_recarray(self, recarrays):
 
             new_node = [remapper[i] for i in node if i in remapper]
 
-            if modelgrid.grid_type == "structured":
+            if self._grid_type == "structured":
                 if self._modelgrid is None:
                     new_cellid = list(
                         zip(*np.unravel_index(new_node, self._shape))
@@ -889,7 +906,7 @@ def reconstruct_recarray(self, recarrays):
                 new_cellid = [
                     (layer[ix], i[1], i[2]) for ix, i in enumerate(new_cellid)
                 ]
-            elif modelgrid.grid_type == "vertex":
+            elif self._grid_type == "vertex":
                 new_cellid = [(layer[ix], i) for ix, i in enumerate(new_node)]
             else:
                 new_cellid = [(i,) for i in new_node]
@@ -1013,7 +1030,7 @@ def _remap_nodes(self, array):
             )
 
         grid_info = {}
-        if self._modelgrid.grid_type == "structured":
+        if self._grid_type == "structured" and not self._to_disv:
             a = array.reshape(self._modelgrid.nrow, self._modelgrid.ncol)
             for m in np.unique(a):
                 cells = np.asarray(a == m).nonzero()
@@ -1038,8 +1055,16 @@ def _remap_nodes(self, array):
                     np.ravel(mapping),
                 ]
         else:
+            # a DISV model has no use for a stack of cells that is inactive
+            # in every layer, so those cells are not carried into it
+            mask = np.ones((self._ncpl,), dtype=bool)
+            if self._to_disv:
+                for arr in idomain:
+                    mask &= arr == 0
+                mask = ~mask
+
             for m in mkeys:
-                cells = np.asarray(array == m).nonzero()[0]
+                cells = np.asarray((array == m) & mask).nonzero()[0]
                 mapping = np.zeros((len(cells),), dtype=int)
                 mapping[:] = cells
                 grid_info[m] = [(len(cells),), None, None, mapping]
@@ -1075,8 +1100,14 @@ def _remap_nodes(self, array):
         exchange_meta = {i: {} for i in mkeys}
         usg_meta = {i: {} for i in mkeys}
         for node, conn in self._connection.items():
+            if node not in self._node_map:
+                continue
+
             mdl, nnode = self._node_map[node]
             for ix, cnode in enumerate(conn):
+                if cnode not in self._node_map:
+                    continue
+
                 cmdl, cnnode = self._node_map[cnode]
                 if cmdl == mdl:
                     if nnode in new_connections[mdl]["internal"]:
@@ -1163,7 +1194,7 @@ def _remap_nodes(self, array):
                                 ]
                             }
 
-        if self._modelgrid.grid_type in ("vertex", "unstructured"):
+        if self._grid_type in ("vertex", "unstructured"):
             self._map_verts_iverts(array)
 
         self._new_connections = new_connections
@@ -1327,6 +1358,58 @@ def _remap_cell2d(self, item, cell2d, mapped_data):
 
         return mapped_data
 
+    def _structured_to_disv(self, mapped_data):
+        """
+        Method to create the DISV grid of models split from a structured model
+
+        Parameters
+        ----------
+        mapped_data : dict
+            dictionary of remapped package data
+
+        Returns
+        -------
+            dict
+        """
+        modelgrid = self._modelgrid
+        ncol = modelgrid.ncol
+        xedge, yedge = modelgrid.xyedges
+        xcenter, ycenter = modelgrid.xycenters
+        vert_dtype = np.dtype([("iv", int), ("xv", float), ("yv", float)])
+        cell_dtype = np.dtype(
+            [("icell2d", int), ("xc", float), ("yc", float), ("ncvert", int)]
+            + [(f"icvert_{iv}", int) for iv in range(4)]
+        )
+
+        for mkey in self._model_dict.keys():
+            nodes = self._grid_info[mkey][-1]
+            # corners shared by cells in the model collapse to one vertex
+            iverts = modelgrid.get_cell_iverts(nodes)
+            overts, iverts = np.unique(iverts, return_inverse=True)
+            iverts = iverts.reshape((nodes.size, 4))
+
+            i, j = np.divmod(overts, ncol + 1)
+            vertices = np.recarray((overts.size,), dtype=vert_dtype)
+            vertices["iv"] = np.arange(overts.size, dtype=int)
+            vertices["xv"] = xedge[j]
+            vertices["yv"] = yedge[i]
+
+            i, j = np.divmod(nodes, ncol)
+            cell2d = np.recarray((nodes.size,), dtype=cell_dtype)
+            cell2d["icell2d"] = np.arange(nodes.size, dtype=int)
+            cell2d["xc"] = xcenter[j]
+            cell2d["yc"] = ycenter[i]
+            cell2d["ncvert"] = 4
+            for iv in range(4):
+                cell2d[f"icvert_{iv}"] = iverts[:, iv]
+
+            mapped_data[mkey]["ncpl"] = self._new_ncpl[mkey]
+            mapped_data[mkey]["nvert"] = overts.size
+            mapped_data[mkey]["vertices"] = vertices
+            mapped_data[mkey]["cell2d"] = cell2d
+
+        return mapped_data
+
     def _remap_filerecords(self, item, value, mapped_data, namfile=False):
         """
         Method to create new file record names and map them to their
@@ -3133,10 +3216,10 @@ def _new_node_to_cellid(self, model, new_node, layers, idx):
         """
 
         new_node = new_node[idx].astype(int)
-        if self._modelgrid.grid_type == "structured":
+        if self._new_grid_type == "structured":
             new_node += layers[idx] * model.modelgrid.ncpl
             new_cellids = model.modelgrid.get_lrc(new_node.astype(int))
-        elif self._modelgrid.grid_type == "vertex":
+        elif self._new_grid_type == "vertex":
             new_cellids = [tuple(cid) for cid in zip(layers[idx], new_node)]
 
         else:
@@ -3279,6 +3362,9 @@ def _remap_package(self, package, ismvr=False):
                     continue
 
                 if item in ("delr", "delc"):
+                    if self._to_disv:
+                        continue
+
                     for mkey, d in self._grid_info.items():
                         if item == "delr":
                             i0, i1 = d[2]
@@ -3288,6 +3374,9 @@ def _remap_package(self, package, ismvr=False):
                         mapped_data[mkey][item] = value.array[i0 : i1 + 1]
 
                 elif item in ("nrow", "ncol"):
+                    if self._to_disv:
+                        continue
+
                     for mkey, d in self._grid_info.items():
                         if item == "nrow":
                             i0, i1 = d[1]
@@ -3330,6 +3419,9 @@ def _remap_package(self, package, ismvr=False):
                 elif isinstance(value, mfdataarray.MFArray):
                     mapped_data = self._remap_array(item, value, mapped_data)
 
+            if self._to_disv:
+                mapped_data = self._structured_to_disv(mapped_data)
+
         elif isinstance(package, modflow.ModflowGwfhfb):
             mapped_data = self._remap_hfb(package, mapped_data)
 
@@ -3550,8 +3642,14 @@ def _remap_package(self, package, ismvr=False):
                         elif isinstance(value, mfdatalist.MFList):
                             mapped_data[mkey][item] = value.array
 
+        package_type = package.package_type
+        pname = package.name[0]
+        if self._to_disv and package_type == "dis":
+            package_type = "disv"
+            pname = None
+
         pak_cls = PackageContainer.package_factory(
-            package.package_type, self._model_type
+            package_type, self._model_type
         )
         paks = {}
         for mdl, data in mapped_data.items():
@@ -3570,7 +3668,7 @@ def _remap_package(self, package, ismvr=False):
                     self._new_sim.simulation_data.max_columns_of_data = max_cols
 
                 paks[mdl] = pak_cls(
-                    self._model_dict[mdl], pname=package.name[0], **data
+                    self._model_dict[mdl], pname=pname, **data
                 )
 
         if observations is not None:
@@ -3723,7 +3821,7 @@ def _create_exchanges(self):
 
                             node1 = exg[1]
                             for layer in range(self._modelgrid.nlay):
-                                if self._modelgrid.grid_type == "structured":
+                                if self._new_grid_type == "structured":
                                     tmpnode0 = node0 + (ncpl0 * layer)
                                     tmpnode1 = node1 + (ncpl1 * layer)
                                     cellidm0 = modelgrid0.get_lrc([tmpnode0])[
@@ -3732,7 +3830,7 @@ def _create_exchanges(self):
                                     cellidm1 = modelgrid1.get_lrc([tmpnode1])[
                                         0
                                     ]
-                                elif self._modelgrid.grid_type == "vertex":
+                                elif self._new_grid_type == "vertex":
                                     cellidm0 = (layer, node0)
                                     cellidm1 = (layer, node1)
                                 else:
@@ -3892,7 +3990,36 @@ def create_multi_model_exchanges(self, mname0, mname1):
             filename=filename,
         )
 
-    def split_model(self, array, sim_ws=None):
+    def _set_new_grid_type(self, to_disv):
+        """
+        Method to set the grid type of the models created by splitting
+
+        Parameters
+        ----------
+        to_disv : bool
+            flag to write DISV models when splitting a structured model
+
+        """
+        self._to_disv = False
+        self._new_grid_type = self._grid_type
+        if not to_disv:
+            return
+
+        if self._grid_type == "unstructured":
+            raise ValueError("DISV models cannot be split from a DISU model")
+
+        if self._grid_type == "vertex":
+            warnings.warn(
+                "the model is already a DISV model, to_disv is ignored",
+                UserWarning,
+                stacklevel=3,
+            )
+            return
+
+        self._to_disv = True
+        self._new_grid_type = "vertex"
+
+    def split_model(self, array, sim_ws=None, to_disv=False):
         """
         User method to split a model based on an array
 
@@ -3906,6 +4033,11 @@ def split_model(self, array, sim_ws=None):
             optional directory path for writing the new simulation to. This parameter
             is recommended when the model contains external files and the user would
             like to preserve external linkages while splitting.
+        to_disv : bool
+            optional flag to write DISV models when splitting a structured model.
+            A DISV model carries only the cells assigned to it, not the inactive
+            cells that pad a structured model to its bounding box. Stacks of cells
+            that are inactive in every layer are excluded. Default is False.
 
         Returns
         -------
@@ -3917,6 +4049,8 @@ def split_model(self, array, sim_ws=None):
                 "is part of a split simulation"
             )
 
+        self._set_new_grid_type(to_disv)
+
         if sim_ws is None:
             self._keep_external = False
             sim_ws = self._sim.sim_path
@@ -3970,7 +4104,7 @@ def split_model(self, array, sim_ws=None):
 
         return self._new_sim
 
-    def split_multi_model(self, array, sim_ws=None):
+    def split_multi_model(self, array, sim_ws=None, to_disv=False):
         """
         Method to split integrated models such as GWF-GWT or GWF-GWE models.
         Note: this method will not work to split multiple connected GWF models
@@ -3985,6 +4119,11 @@ def split_multi_model(self, array, sim_ws=None):
             optional directory path for writing the new simulation to. This parameter
             is recommended when the model contains external files and the user would
             like to preserve external linkages while splitting.
+        to_disv : bool
+            optional flag to write DISV models when splitting a structured model.
+            A DISV model carries only the cells assigned to it, not the inactive
+            cells that pad a structured model to its bounding box. Stacks of cells
+            that are inactive in every layer are excluded. Default is False.
 
         Returns
         -------
@@ -4048,10 +4187,10 @@ def split_multi_model(self, array, sim_ws=None):
             int(i): f"{gwf_base}_{i}" for i in model_labels
         }
 
-        new_sim = self.split_model(array)
+        new_sim = self.split_model(array, to_disv=to_disv)
         for mname in model_names[1:]:
             self.switch_models(modelname=mname, remap_nodes=False)
-            new_sim = self.split_model(array, sim_ws=sim_ws)
+            new_sim = self.split_model(array, sim_ws=sim_ws, to_disv=to_disv)
 
         for mbase in model_names[1:]:
             for label in model_labels:

From e98301ef6883977e128c081696257e29327dc8ef Mon Sep 17 00:00:00 2001
From: Joseph Hughes <jdhughes99+1@gmail.com>
Date: Sun, 23 Aug 2026 00:32:09 -0500
Subject: [PATCH 02/10] fix(structuredgrid): flatten and coerce nodes in
 get_cell_iverts()

np.atleast_1d() left a node array with a trailing dimension, for example shape
(n, 1), two dimensional, and the vertex arithmetic then failed to broadcast into
the one dimensional result rows. Ravel the node numbers and coerce them to int so
an array of any shape works as documented. The ravel is a view and asarray() is a
no-op for an int array, so the array is not copied.
---
 autotest/test_grid.py                  | 4 ++++
 flopy/discretization/structuredgrid.py | 2 +-
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index b3bcfee69..8861789fa 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -295,6 +295,10 @@ def test_get_cell_iverts():
 
     nodes = [0, 7, grid.ncpl - 1]
     assert np.array_equal(grid.get_cell_iverts(nodes), iverts[nodes])
+    assert np.array_equal(grid.get_cell_iverts(np.array(nodes)), iverts[nodes])
+    assert np.array_equal(
+        grid.get_cell_iverts(np.array(nodes).reshape((-1, 1))), iverts[nodes]
+    )
     assert np.array_equal(grid.get_cell_iverts(3), iverts[[3]])
 
     # vertices are ordered clockwise from the upper left corner of a cell
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index a6561f3a2..077017073 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -2021,7 +2021,7 @@ def get_cell_iverts(self, nodes=None):
         if nodes is None:
             nodes = np.arange(self.ncpl, dtype=int)
         else:
-            nodes = np.atleast_1d(nodes)
+            nodes = np.ravel(np.asarray(nodes, dtype=int))
 
         i, j = np.divmod(nodes, self.ncol)
         iverts = np.empty((4, nodes.size), dtype=int)

From 00e5aae347bf3cb014b5ebe234abdd746fab1511 Mon Sep 17 00:00:00 2001
From: Joseph Hughes <jdhughes99+1@gmail.com>
Date: Sun, 23 Aug 2026 01:37:43 -0500
Subject: [PATCH 03/10] refactor(Mf6Splitter): read the parent grid type from
 the modelgrid

switch_models() already updates the modelgrid, so caching a second copy of its
grid type is redundant and goes stale when the model is switched. Read
self._modelgrid.grid_type where the grid type of the parent model is needed, and
leave switch_models() alone.

Also cover reconstruct_recarray() for a vertex model, which had no test.
---
 autotest/test_model_splitter.py   | 33 +++++++++++++++++++++++++++++++
 flopy/mf6/utils/model_splitter.py | 17 ++++++++--------
 2 files changed, 41 insertions(+), 9 deletions(-)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index f1d75b474..e41afc90e 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -304,6 +304,39 @@ def test_structured_to_disv_reconstruct_recarray():
             )
 
 
+def test_vertex_reconstruct_recarray():
+    sim_path = get_example_data_path() / "mf6" / "test003_gwftri_disv"
+
+    sim = MFSimulation.load(sim_ws=sim_path)
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+
+    array = np.zeros((modelgrid.ncpl,), dtype=int)
+    array[0:85] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_model(array)
+
+    for pname in ("chd_left", "chd_right"):
+        original = gwf.get_package(pname).stress_period_data.get_data(0)
+        recarrays = {}
+        for mkey in (0, 1):
+            ml = new_sim.get_model(f"gwf_1_{mkey}")
+            pkg = ml.get_package(pname)
+            if pkg is not None:
+                recarrays[mkey] = pkg.stress_period_data.get_data(0)
+
+        new_recarray = mfsplit.reconstruct_recarray(recarrays)
+
+        # cellids are returned as the layer, node of the original model
+        onodes = sorted(cid[-1] for cid in original.cellid)
+        nnodes = sorted(cid[-1] for cid in new_recarray.cellid)
+
+        err_msg = f"Reconstructed {pname} recarray does not match the original"
+        if onodes != nnodes:
+            raise AssertionError(err_msg)
+
+
 def test_structured_to_disv_model_geometry():
     sim_path = get_example_data_path() / "mf6-freyberg"
 
diff --git a/flopy/mf6/utils/model_splitter.py b/flopy/mf6/utils/model_splitter.py
index afe98d51f..137b041e2 100644
--- a/flopy/mf6/utils/model_splitter.py
+++ b/flopy/mf6/utils/model_splitter.py
@@ -226,8 +226,6 @@ def switch_models(self, modelname, remap_nodes=False):
 
         if remap_nodes:
             self._modelgrid = self._model.modelgrid
-            self._grid_type = self._modelgrid.grid_type
-            self._new_grid_type = self._grid_type
             self._node_map = {}
             self._node_map_r = {}
             self._new_connections = None
@@ -896,7 +894,7 @@ def reconstruct_recarray(self, recarrays):
 
             new_node = [remapper[i] for i in node if i in remapper]
 
-            if self._grid_type == "structured":
+            if self._modelgrid.grid_type == "structured":
                 if self._modelgrid is None:
                     new_cellid = list(
                         zip(*np.unravel_index(new_node, self._shape))
@@ -906,7 +904,7 @@ def reconstruct_recarray(self, recarrays):
                 new_cellid = [
                     (layer[ix], i[1], i[2]) for ix, i in enumerate(new_cellid)
                 ]
-            elif self._grid_type == "vertex":
+            elif self._modelgrid.grid_type == "vertex":
                 new_cellid = [(layer[ix], i) for ix, i in enumerate(new_node)]
             else:
                 new_cellid = [(i,) for i in new_node]
@@ -1030,7 +1028,7 @@ def _remap_nodes(self, array):
             )
 
         grid_info = {}
-        if self._grid_type == "structured" and not self._to_disv:
+        if self._modelgrid.grid_type == "structured" and not self._to_disv:
             a = array.reshape(self._modelgrid.nrow, self._modelgrid.ncol)
             for m in np.unique(a):
                 cells = np.asarray(a == m).nonzero()
@@ -1194,7 +1192,7 @@ def _remap_nodes(self, array):
                                 ]
                             }
 
-        if self._grid_type in ("vertex", "unstructured"):
+        if self._modelgrid.grid_type in ("vertex", "unstructured"):
             self._map_verts_iverts(array)
 
         self._new_connections = new_connections
@@ -4000,15 +3998,16 @@ def _set_new_grid_type(self, to_disv):
             flag to write DISV models when splitting a structured model
 
         """
+        grid_type = self._modelgrid.grid_type
         self._to_disv = False
-        self._new_grid_type = self._grid_type
+        self._new_grid_type = grid_type
         if not to_disv:
             return
 
-        if self._grid_type == "unstructured":
+        if grid_type == "unstructured":
             raise ValueError("DISV models cannot be split from a DISU model")
 
-        if self._grid_type == "vertex":
+        if grid_type == "vertex":
             warnings.warn(
                 "the model is already a DISV model, to_disv is ignored",
                 UserWarning,

From 479baf5c5f722018a2d20079f311cea400920f78 Mon Sep 17 00:00:00 2001
From: Joseph Hughes <jdhughes99+1@gmail.com>
Date: Sun, 23 Aug 2026 09:30:59 -0500
Subject: [PATCH 04/10] test(model_splitter): cover advanced packages and
 transient data with to_disv

The advanced packages remap cellids through _new_node_to_cellid(), which returns
a layer and node pair for a DISV model, and no test exercised that path. Split the
lake2tr model, which has LAK, SFR, MVR, observation, and array based recharge and
evapotranspiration packages, to DISV models. Run the transient array test for both
grid types so the transient arrays and stress period data are covered as well.
---
 autotest/test_model_splitter.py | 53 +++++++++++++++++++++++++++++++--
 1 file changed, 51 insertions(+), 2 deletions(-)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index e41afc90e..8ebfffeef 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -484,6 +484,54 @@ def test_model_with_lak_sfr_mvr(function_tmpdir):
     np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg)
 
 
+@requires_exe("mf6")
+@pytest.mark.slow
+def test_structured_to_disv_with_lak_sfr_mvr(function_tmpdir):
+    sim_path = get_example_data_path() / "mf6" / "test045_lake2tr"
+
+    sim = MFSimulation.load(sim_ws=sim_path)
+    sim.set_sim_path(function_tmpdir)
+    sim.write_simulation()
+    sim.run_simulation()
+
+    gwf = sim.get_model()
+    modelgrid = gwf.modelgrid
+
+    array = np.zeros((modelgrid.nrow, modelgrid.ncol), dtype=int)
+    array[0:14, :] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_model(array, to_disv=True)
+
+    new_sim.set_sim_path(function_tmpdir / "split_model")
+    new_sim.write_simulation()
+    new_sim.run_simulation()
+
+    original_heads = gwf.output.head().get_alldata()[-1]
+
+    ml0 = new_sim.get_model("lakeex2a_0")
+    ml1 = new_sim.get_model("lakeex2a_1")
+    for ml in (ml0, ml1):
+        if ml.modelgrid.grid_type != "vertex":
+            raise AssertionError("Split model is not a DISV model")
+
+        # the advanced packages remap cellids through the new grid type
+        for pkgtype in ("lak", "sfr", "mvr", "evta", "rcha"):
+            if ml.get_package(pkgtype) is None:
+                raise AssertionError(f"{pkgtype} package was not split")
+
+    heads0 = ml0.output.head().get_alldata()[-1]
+    heads1 = ml1.output.head().get_alldata()[-1]
+
+    new_heads = mfsplit.reconstruct_array({0: heads0, 1: heads1})
+
+    idx = modelgrid.idomain != 0
+    err_msg = "Heads from original and split models do not match"
+    np.testing.assert_allclose(
+        new_heads[idx], original_heads[idx], atol=1e-4, err_msg=err_msg
+    )
+
+
 @requires_exe("mf6")
 @requires_pkg("pymetis")
 @pytest.mark.slow
@@ -969,7 +1017,8 @@ def test_empty_ssm(function_tmpdir):
 
 
 @requires_exe("mf6")
-def test_transient_array(function_tmpdir):
+@pytest.mark.parametrize("to_disv", [False, True])
+def test_transient_array(function_tmpdir, to_disv):
     name = "tarr"
     new_sim_path = function_tmpdir / f"{name}_split_model"
     nper = 3
@@ -1047,7 +1096,7 @@ def test_transient_array(function_tmpdir):
     sarr = np.ones((nrow, ncol), dtype=int)
     sarr[:, int(ncol / 2) :] = 2
     mfsplit = Mf6Splitter(sim)
-    new_sim = mfsplit.split_model(sarr)
+    new_sim = mfsplit.split_model(sarr, to_disv=to_disv)
 
     for name in new_sim.model_names:
         g = new_sim.get_model(name)

From 9e55a6ff717aab0cde9f0586332139d10845a56d Mon Sep 17 00:00:00 2001
From: Joseph Hughes <jdhughes99+1@gmail.com>
Date: Sun, 23 Aug 2026 09:34:41 -0500
Subject: [PATCH 05/10] test(model_splitter): cover split_multi_model() with
 to_disv

Split a GWF-GWT simulation into DISV models and check that both models of each
pair are DISV and that the concentrations reconstruct to the original model.
---
 autotest/test_model_splitter.py | 90 +++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index 8ebfffeef..7d446b5cf 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -484,6 +484,96 @@ def test_model_with_lak_sfr_mvr(function_tmpdir):
     np.testing.assert_allclose(new_heads, original_heads, err_msg=err_msg)
 
 
+@requires_exe("mf6")
+def test_structured_to_disv_multi_model(function_tmpdir):
+    nlay, nrow, ncol = 1, 10, 10
+    idomain = np.ones((nlay, nrow, ncol), dtype=int)
+    idomain[:, 0:2, 0:2] = 0
+
+    sim = flopy.mf6.MFSimulation(sim_name="mm", sim_ws=function_tmpdir, exe_name="mf6")
+    flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(100.0, 10, 1.0)])
+    ims_gwf = flopy.mf6.ModflowIms(sim, complexity="simple", filename="gwf.ims")
+    ims_gwt = flopy.mf6.ModflowIms(
+        sim,
+        complexity="simple",
+        linear_acceleration="bicgstab",
+        filename="gwt.ims",
+    )
+
+    dis_kwargs = {
+        "nlay": nlay,
+        "nrow": nrow,
+        "ncol": ncol,
+        "delr": 100.0,
+        "delc": 100.0,
+        "top": 10.0,
+        "botm": [0.0],
+        "idomain": idomain,
+    }
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="gwf", save_flows=True)
+    flopy.mf6.ModflowGwfdis(gwf, **dis_kwargs)
+    flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True, k=1.0)
+    flopy.mf6.ModflowGwfic(gwf, strt=10.0)
+    flopy.mf6.ModflowGwfsto(gwf, ss=1e-5, iconvert=0)
+    chd = [[(0, i, 2), 10.0, 1.0] for i in range(2, nrow)]
+    chd += [[(0, i, ncol - 1), 9.0, 0.0] for i in range(nrow)]
+    flopy.mf6.ModflowGwfchd(
+        gwf, pname="chd-1", auxiliary=["conc"], stress_period_data=chd
+    )
+    flopy.mf6.ModflowGwfoc(
+        gwf,
+        head_filerecord="gwf.hds",
+        budget_filerecord="gwf.cbc",
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    gwt = flopy.mf6.ModflowGwt(sim, modelname="gwt")
+    flopy.mf6.ModflowGwtdis(gwt, **dis_kwargs)
+    flopy.mf6.ModflowGwtic(gwt, strt=0.0)
+    flopy.mf6.ModflowGwtmst(gwt, porosity=0.2)
+    flopy.mf6.ModflowGwtadv(gwt, scheme="upstream")
+    flopy.mf6.ModflowGwtssm(gwt, sources=[["chd-1", "AUX", "conc"]])
+    flopy.mf6.ModflowGwtoc(
+        gwt,
+        concentration_filerecord="gwt.ucn",
+        saverecord=[("CONCENTRATION", "ALL")],
+    )
+    flopy.mf6.ModflowGwfgwt(sim, exgmnamea="gwf", exgmnameb="gwt")
+    sim.register_ims_package(ims_gwf, ["gwf"])
+    sim.register_ims_package(ims_gwt, ["gwt"])
+    sim.write_simulation()
+    sim.run_simulation()
+
+    original_conc = gwt.output.concentration().get_alldata()[-1]
+
+    array = np.zeros((nrow, ncol), dtype=int)
+    array[:, ncol // 2 :] = 1
+
+    mfsplit = Mf6Splitter(sim)
+    new_sim = mfsplit.split_multi_model(array, to_disv=True)
+    new_sim.set_sim_path(function_tmpdir / "split_model")
+    new_sim.write_simulation()
+    new_sim.run_simulation()
+
+    conc = {}
+    for mkey in (0, 1):
+        for mname in (f"gwf_{mkey}", f"gwt_{mkey}"):
+            if new_sim.get_model(mname).modelgrid.grid_type != "vertex":
+                raise AssertionError(f"Model {mname} is not a DISV model")
+
+        conc[mkey] = new_sim.get_model(f"gwt_{mkey}").output.concentration()
+        conc[mkey] = conc[mkey].get_alldata()[-1]
+
+    new_conc = mfsplit.reconstruct_array(conc)
+
+    idx = idomain != 0
+    err_msg = "Concentrations from original and split models do not match"
+    np.testing.assert_allclose(
+        new_conc[idx], original_conc[idx], atol=1e-6, err_msg=err_msg
+    )
+
+
 @requires_exe("mf6")
 @pytest.mark.slow
 def test_structured_to_disv_with_lak_sfr_mvr(function_tmpdir):

From b046319f414e4d4e204c621677c349d8fa8ab19e Mon Sep 17 00:00:00 2001
From: Joseph Hughes <jdhughes99+1@gmail.com>
Date: Sun, 23 Aug 2026 09:42:13 -0500
Subject: [PATCH 06/10] test(model_splitter): cover idomain vertical
 passthrough cells with to_disv

A cell with an idomain of -1 passes flow vertically and is not inactive, so its
stack is carried into the DISV model and its idomain is kept. Add passthrough
cells to the idomain test and compare heads only where cells are active.
---
 autotest/test_model_splitter.py | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index 7d446b5cf..14556d44a 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -192,6 +192,7 @@ def test_structured_to_disv_idomain(function_tmpdir):
     idomain[:, 0:2, 0:2] = 0  # inactive in every layer, excluded from the split
     idomain[0, 5:, 5:] = 0  # inactive in one layer, carried into the split
     idomain[2, 0:3, 7:] = 0
+    idomain[1, 3:7, 3:7] = -1  # vertical passthrough, carried into the split
 
     sim = flopy.mf6.MFSimulation(sim_name="ml", sim_ws=function_tmpdir, exe_name="mf6")
     flopy.mf6.ModflowTdis(sim)
@@ -243,10 +244,14 @@ def test_structured_to_disv_idomain(function_tmpdir):
         nodes = mfsplit._grid_info[mkey][-1]
         ncpl += grid.ncpl
 
-        # a stack that is inactive in one layer keeps its idomain
-        ninactive = np.count_nonzero(oidomain[:, nodes] == 0)
-        if np.count_nonzero(grid.idomain == 0) != ninactive:
-            raise AssertionError(f"Model {mkey} idomain was not remapped")
+        # a stack that is inactive or passes flow vertically in one layer
+        # keeps its idomain
+        for value in (0, -1):
+            nexpected = np.count_nonzero(oidomain[:, nodes] == value)
+            if np.count_nonzero(grid.idomain == value) != nexpected:
+                raise AssertionError(
+                    f"Model {mkey} idomain {value} cells were not remapped"
+                )
 
         heads[mkey] = ml.output.head().get_alldata()[-1]
 
@@ -257,7 +262,7 @@ def test_structured_to_disv_idomain(function_tmpdir):
 
     new_heads = mfsplit.reconstruct_array(heads)
 
-    idx = idomain != 0
+    idx = idomain > 0
     err_msg = "Heads from original and split models do not match"
     np.testing.assert_allclose(new_heads[idx], original_heads[idx], err_msg=err_msg)
 

From f8533a3d8f8878e445f138da679660ad1896eb9d Mon Sep 17 00:00:00 2001
From: Joshua Larsen <jlarsen@usgs.gov>
Date: Thu, 3 Sep 2026 11:25:10 -0700
Subject: [PATCH 07/10] Update test_model_splitter.py

update tests to call `_node_map_arr`  to follow new node mapping data type convention
---
 autotest/test_model_splitter.py | 25 ++++++++++++++-----------
 1 file changed, 14 insertions(+), 11 deletions(-)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index 7f2d15ac5..8ab3f7ca8 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -429,16 +429,18 @@ def test_save_load_node_mapping_structured_to_disv(function_tmpdir):
     new_sim.set_sim_path(new_sim_path)
     new_sim.write_simulation()
     new_sim.run_simulation()
-    original_node_map = mfsplit._node_map
+    original_node_map = mfsplit._node_map_arr
 
     mfsplit.save_node_mapping(hdf_file)
 
     new_sim2 = MFSimulation.load(sim_ws=new_sim_path)
     mfsplit2 = Mf6Splitter.load_node_mapping(hdf_file)
 
-    if original_node_map != mfsplit2._node_map:
-        raise AssertionError("Node map read/write not returning proper values")
-
+    np.testing.assert_allclose(
+        original_node_map,
+        mfsplit2._node_map_arr,
+        err_msg="Node map read/write not returning proper values",
+    )
     array_dict = {}
     for mkey in (0, 1):
         ml = new_sim2.get_model(f"freyberg_{mkey}")
@@ -698,20 +700,21 @@ def test_save_load_node_mapping_structured(function_tmpdir):
     new_sim.set_sim_path(new_sim_path)
     new_sim.write_simulation()
     new_sim.run_simulation()
-    original_node_map = mfsplit._node_map
+    original_node_map = mfsplit._node_map_arr
 
     mfsplit.save_node_mapping(hdf_file)
 
     new_sim2 = MFSimulation.load(sim_ws=new_sim_path)
 
     mfsplit2 = Mf6Splitter.load_node_mapping(hdf_file)
-    saved_node_map = mfsplit2._node_map
-
-    for k, v1 in original_node_map.items():
-        v2 = saved_node_map[k]
-        if not v1 == v2:
-            raise AssertionError("Node map read/write not returning proper values")
+    saved_node_map = mfsplit2._node_map_arr
 
+    np.testing.assert_allclose(
+        original_node_map,
+        saved_node_map,
+        err_msg="Node map read/write not returning proper values",
+    )
+    
     array_dict = {}
     for model in range(nparts):
         ml = new_sim2.get_model(f"freyberg_{model}")

From becb77db701d4b6555ad006cf73a69df9fdcd9f2 Mon Sep 17 00:00:00 2001
From: Joshua Larsen <jlarsen@usgs.gov>
Date: Thu, 3 Sep 2026 11:28:53 -0700
Subject: [PATCH 08/10] linting update

---
 autotest/test_model_splitter.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index e0fd20efc..ff58248f4 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -441,6 +441,7 @@ def test_save_load_node_mapping_structured_to_disv(function_tmpdir):
         mfsplit2._node_map_arr,
         err_msg="Node map read/write not returning proper values",
     )
+    
     array_dict = {}
     for mkey in (0, 1):
         ml = new_sim2.get_model(f"freyberg_{mkey}")

From 72a64e68e0c679333e8f4fda6f8dff0d8397d34f Mon Sep 17 00:00:00 2001
From: Joshua Larsen <jlarsen@usgs.gov>
Date: Thu, 3 Sep 2026 11:42:44 -0700
Subject: [PATCH 09/10] linting update

---
 autotest/test_model_splitter.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index ff58248f4..0b0e37460 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -714,7 +714,7 @@ def test_save_load_node_mapping_structured(function_tmpdir):
         saved_node_map,
         err_msg="Node map read/write not returning proper values",
     )
-    
+
     array_dict = {}
     for model in range(nparts):
         ml = new_sim2.get_model(f"freyberg_{model}")

From c08d8beed630a0b547a190eea88e135f9c3077e4 Mon Sep 17 00:00:00 2001
From: Joshua Larsen <jlarsen@usgs.gov>
Date: Thu, 3 Sep 2026 14:31:09 -0700
Subject: [PATCH 10/10] linting updates

---
 autotest/test_model_splitter.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index 202bdfbb9..763946e4c 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -580,7 +580,6 @@ def test_structured_to_disv_multi_model(function_tmpdir):
     np.testing.assert_allclose(
         new_conc[idx], original_conc[idx], atol=1e-6, err_msg=err_msg
     )
-    
 
 
 def test_hfb_model_splitter(function_tmpdir):
@@ -636,7 +635,7 @@ def test_hfb_model_splitter(function_tmpdir):
     new_sim.set_sim_path(function_tmpdir / "split_model")
     new_sim.write_simulation()
     new_sim.run_simulation()
-    
+
     heads = {}
     nbarrier = 0
     for mkey in (0, 1):