diff --git a/autotest/regression/test_mf6.py b/autotest/regression/test_mf6.py
index f46e47193..3ffc52f9a 100644
--- a/autotest/regression/test_mf6.py
+++ b/autotest/regression/test_mf6.py
@@ -3600,6 +3600,187 @@ def test005_advgw_tidal(function_tmpdir, example_data_path):
     )
 
 
+@requires_exe("mf6")
+@pytest.mark.regression
+def test006_2models_different_dis(function_tmpdir, example_data_path):
+    # init paths
+    test_ex_name = "test006_2models_diff_dis"
+    model_name_1 = "model1"
+    model_name_2 = "model2"
+    pth = example_data_path / "mf6" / "create_tests" / test_ex_name
+
+    expected_output_folder = os.path.join(pth, "expected_output")
+    expected_head_file_1 = os.path.join(expected_output_folder, "model1.hds")
+    expected_head_file_2 = os.path.join(expected_output_folder, "model2.hds")
+
+    # create simulation
+    sim = MFSimulation(
+        sim_name=test_ex_name, version="mf6", exe_name="mf6", sim_ws=pth
+    )
+    tdis_rc = [(1.0, 1, 1.0)]
+    tdis_package = ModflowTdis(
+        sim, time_units="DAYS", nper=1, perioddata=tdis_rc
+    )
+    model_1 = ModflowGwf(
+        sim,
+        modelname=model_name_1,
+        model_nam_file=f"{model_name_1}.nam",
+    )
+    model_2 = ModflowGwf(
+        sim,
+        modelname=model_name_2,
+        model_nam_file=f"{model_name_2}.nam",
+    )
+    ims_package = ModflowIms(
+        sim,
+        print_option="SUMMARY",
+        outer_dvclose=0.00000001,
+        outer_maximum=1000,
+        under_relaxation="NONE",
+        inner_maximum=1000,
+        inner_dvclose=0.00000001,
+        rcloserecord=0.01,
+        linear_acceleration="BICGSTAB",
+        scaling_method="NONE",
+        reordering_method="NONE",
+        relaxation_factor=0.97,
+    )
+    sim.register_ims_package(ims_package, [model_1.name, model_2.name])
+    dis_package = ModflowGwfdis(
+        model_1,
+        length_units="METERS",
+        nlay=1,
+        nrow=7,
+        ncol=7,
+        idomain=1,
+        delr=100.0,
+        delc=100.0,
+        top=0.0,
+        botm=-100.0,
+        filename=f"{model_name_1}.dis",
+    )
+
+    vertices = testutils.read_vertices(os.path.join(pth, "vertices.txt"))
+    c2drecarray = testutils.read_cell2d(os.path.join(pth, "cell2d.txt"))
+    disv_package = ModflowGwfdisv(
+        model_2,
+        ncpl=121,
+        nlay=1,
+        nvert=148,
+        top=0.0,
+        botm=-40.0,
+        idomain=1,
+        vertices=vertices,
+        cell2d=c2drecarray,
+        filename=f"{model_name_2}.disv",
+    )
+    ic_package_1 = ModflowGwfic(
+        model_1, strt=1.0, filename=f"{model_name_1}.ic"
+    )
+    ic_package_2 = ModflowGwfic(
+        model_2, strt=1.0, filename=f"{model_name_2}.ic"
+    )
+    npf_package_1 = ModflowGwfnpf(
+        model_1, save_flows=True, perched=True, icelltype=0, k=1.0, k33=1.0
+    )
+    npf_package_2 = ModflowGwfnpf(
+        model_2, save_flows=True, perched=True, icelltype=0, k=1.0, k33=1.0
+    )
+    oc_package_1 = ModflowGwfoc(
+        model_1,
+        budget_filerecord="model1.cbc",
+        head_filerecord="model1.hds",
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+        printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+    oc_package_2 = ModflowGwfoc(
+        model_2,
+        budget_filerecord="model2.cbc",
+        head_filerecord="model2.hds",
+        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+        printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    # build periodrecarray for chd package
+    set_1 = [0, 7, 14, 18, 22, 26, 33]
+    set_2 = [6, 13, 17, 21, 25, 32, 39]
+    stress_period_data = []
+    for value in range(0, 7):
+        stress_period_data.append(((0, value, 0), 1.0))
+    for value in range(0, 7):
+        stress_period_data.append(((0, value, 6), 0.0))
+    chd_package = ModflowGwfchd(
+        model_1,
+        print_input=True,
+        print_flows=True,
+        save_flows=True,
+        maxbound=30,
+        stress_period_data=stress_period_data,
+    )
+    exgrecarray = testutils.read_exchangedata(
+        os.path.join(pth, "exg.txt"), 3, 2
+    )
+
+    # build obs dictionary
+    gwf_obs = {
+        ("gwfgwf_obs.csv"): [
+            ("gwf-1-3-2_1-1-1", "flow-ja-face", (0, 2, 1), (0, 0, 0)),
+            ("gwf-1-3-2_1-2-1", "flow-ja-face", (0, 2, 1), (0, 1, 0)),
+        ]
+    }
+
+    exg_package = ModflowGwfgwf(
+        sim,
+        print_input=True,
+        print_flows=True,
+        save_flows=True,
+        auxiliary="testaux",
+        nexg=9,
+        exchangedata=exgrecarray,
+        exgtype="gwf6-gwf6",
+        exgmnamea=model_name_1,
+        exgmnameb=model_name_2,
+        observations=gwf_obs,
+    )
+
+    gnc_path = os.path.join("gnc", "test006_2models_gnc.gnc")
+    gncrecarray = testutils.read_gncrecarray(
+        os.path.join(pth, "gnc.txt"), 3, 2
+    )
+    gnc_package = exg_package.gnc.initialize(
+        filename=gnc_path,
+        print_input=True,
+        print_flows=True,
+        numgnc=9,
+        numalphaj=1,
+        gncdata=gncrecarray,
+    )
+
+    # change folder to save simulation
+    sim.set_sim_path(function_tmpdir)
+
+    # write simulation to new location
+    sim.write_simulation()
+    # run simulation
+    success, buff = sim.run_simulation()
+    assert success
+
+    sim2 = MFSimulation.load(sim_ws=sim.sim_path)
+    exh = sim2.get_package("gwfgwf")
+    exh_data = exh.exchangedata.get_data()
+    assert exh_data[0][0] == (0, 2, 1)
+    assert exh_data[0][1] == (0, 0)
+    assert exh_data[3][0] == (0, 3, 1)
+    assert exh_data[3][1] == (0, 3)
+    gnc = sim2.get_package("gnc")
+    gnc_data = gnc.gncdata.get_data()
+    assert gnc_data[0][0] == (0, 2, 1)
+    assert gnc_data[0][1] == (0, 0)
+    assert gnc_data[0][2] == (0, 1, 1)
+
+    sim.delete_output_files()
+
+
 @requires_exe("mf6")
 @pytest.mark.regression
 def test006_gwf3(function_tmpdir, example_data_path):
diff --git a/examples/data/mf6/create_tests/test006_2models_diff_dis/cell2d.txt b/examples/data/mf6/create_tests/test006_2models_diff_dis/cell2d.txt
new file mode 100644
index 000000000..7e65843df
--- /dev/null
+++ b/examples/data/mf6/create_tests/test006_2models_diff_dis/cell2d.txt
@@ -0,0 +1,121 @@
+1 50.0 650.0 5 1 2 10 9 1 
+2 150.0 650.0 5 2 3 11 10 2 
+3 250.0 650.0 5 3 4 12 11 3 
+4 350.0 650.0 5 4 5 13 12 4 
+5 450.0 650.0 5 5 6 14 13 5 
+6 550.0 650.0 5 6 7 15 14 6 
+7 650.0 650.0 5 7 8 16 15 7 
+8 50.0 550.0 5 9 10 18 17 9 
+9 150.0 550.0 5 10 11 19 18 10 
+10 250.0 550.0 7 11 12 22 21 20 19 11 
+11 350.0 550.0 7 12 13 25 24 23 22 12 
+12 450.0 550.0 7 13 14 28 27 26 25 13 
+13 550.0 550.0 5 14 15 29 28 14 
+14 650.0 550.0 5 15 16 30 29 15 
+15 50.0 450.0 5 17 18 52 51 17 
+16 150.0 450.0 7 18 19 31 41 53 52 18 
+17 550.0 450.0 7 28 29 63 62 50 40 28 
+18 650.0 450.0 5 29 30 64 63 29 
+19 50.0 350.0 5 51 52 86 85 51 
+20 150.0 350.0 7 52 53 65 75 87 86 52 
+21 550.0 350.0 7 62 63 97 96 84 74 62 
+22 650.0 350.0 5 63 64 98 97 63 
+23 50.0 250.0 5 85 86 120 119 85 
+24 150.0 250.0 7 86 87 99 109 121 120 86 
+25 550.0 250.0 7 96 97 131 130 118 108 96 
+26 650.0 250.0 5 97 98 132 131 97 
+27 50.0 150.0 5 119 120 134 133 119 
+28 150.0 150.0 5 120 121 135 134 120 
+29 250.0 150.0 7 121 122 123 124 136 135 121 
+30 350.0 150.0 7 124 125 126 127 137 136 124 
+31 450.0 150.0 7 127 128 129 130 138 137 127 
+32 550.0 150.0 5 130 131 139 138 130 
+33 650.0 150.0 5 131 132 140 139 131 
+34 50.0 50.0 5 133 134 142 141 133 
+35 150.0 50.0 5 134 135 143 142 134 
+36 250.0 50.0 5 135 136 144 143 135 
+37 350.0 50.0 5 136 137 145 144 136 
+38 450.0 50.0 5 137 138 146 145 137 
+39 550.0 50.0 5 138 139 147 146 138 
+40 650.0 50.0 5 139 140 148 147 139 
+41 216.666666667 483.333333333 5 19 20 32 31 19 
+42 250.0 483.333333333 5 20 21 33 32 20 
+43 283.333333333 483.333333333 5 21 22 34 33 21 
+44 316.666666667 483.333333333 5 22 23 35 34 22 
+45 350.0 483.333333333 5 23 24 36 35 23 
+46 383.333333333 483.333333333 5 24 25 37 36 24 
+47 416.666666667 483.333333333 5 25 26 38 37 25 
+48 450.0 483.333333333 5 26 27 39 38 26 
+49 483.333333333 483.333333333 5 27 28 40 39 27 
+50 216.666666667 450.0 5 31 32 42 41 31 
+51 250.0 450.0 5 32 33 43 42 32 
+52 283.333333333 450.0 5 33 34 44 43 33 
+53 316.666666667 450.0 5 34 35 45 44 34 
+54 350.0 450.0 5 35 36 46 45 35 
+55 383.333333333 450.0 5 36 37 47 46 36 
+56 416.666666667 450.0 5 37 38 48 47 37 
+57 450.0 450.0 5 38 39 49 48 38 
+58 483.333333333 450.0 5 39 40 50 49 39 
+59 216.666666667 416.666666667 5 41 42 54 53 41 
+60 250.0 416.666666667 5 42 43 55 54 42 
+61 283.333333333 416.666666667 5 43 44 56 55 43 
+62 316.666666667 416.666666667 5 44 45 57 56 44 
+63 350.0 416.666666667 5 45 46 58 57 45 
+64 383.333333333 416.666666667 5 46 47 59 58 46 
+65 416.666666667 416.666666667 5 47 48 60 59 47 
+66 450.0 416.666666667 5 48 49 61 60 48 
+67 483.333333333 416.666666667 5 49 50 62 61 49 
+68 216.666666667 383.333333333 5 53 54 66 65 53 
+69 250.0 383.333333333 5 54 55 67 66 54 
+70 283.333333333 383.333333333 5 55 56 68 67 55 
+71 316.666666667 383.333333333 5 56 57 69 68 56 
+72 350.0 383.333333333 5 57 58 70 69 57 
+73 383.333333333 383.333333333 5 58 59 71 70 58 
+74 416.666666667 383.333333333 5 59 60 72 71 59 
+75 450.0 383.333333333 5 60 61 73 72 60 
+76 483.333333333 383.333333333 5 61 62 74 73 61 
+77 216.666666667 350.0 5 65 66 76 75 65 
+78 250.0 350.0 5 66 67 77 76 66 
+79 283.333333333 350.0 5 67 68 78 77 67 
+80 316.666666667 350.0 5 68 69 79 78 68 
+81 350.0 350.0 5 69 70 80 79 69 
+82 383.333333333 350.0 5 70 71 81 80 70 
+83 416.666666667 350.0 5 71 72 82 81 71 
+84 450.0 350.0 5 72 73 83 82 72 
+85 483.333333333 350.0 5 73 74 84 83 73 
+86 216.666666667 316.666666667 5 75 76 88 87 75 
+87 250.0 316.666666667 5 76 77 89 88 76 
+88 283.333333333 316.666666667 5 77 78 90 89 77 
+89 316.666666667 316.666666667 5 78 79 91 90 78 
+90 350.0 316.666666667 5 79 80 92 91 79 
+91 383.333333333 316.666666667 5 80 81 93 92 80 
+92 416.666666667 316.666666667 5 81 82 94 93 81 
+93 450.0 316.666666667 5 82 83 95 94 82 
+94 483.333333333 316.666666667 5 83 84 96 95 83 
+95 216.666666667 283.333333333 5 87 88 100 99 87 
+96 250.0 283.333333333 5 88 89 101 100 88 
+97 283.333333333 283.333333333 5 89 90 102 101 89 
+98 316.666666667 283.333333333 5 90 91 103 102 90 
+99 350.0 283.333333333 5 91 92 104 103 91 
+100 383.333333333 283.333333333 5 92 93 105 104 92 
+101 416.666666667 283.333333333 5 93 94 106 105 93 
+102 450.0 283.333333333 5 94 95 107 106 94 
+103 483.333333333 283.333333333 5 95 96 108 107 95 
+104 216.666666667 250.0 5 99 100 110 109 99 
+105 250.0 250.0 5 100 101 111 110 100 
+106 283.333333333 250.0 5 101 102 112 111 101 
+107 316.666666667 250.0 5 102 103 113 112 102 
+108 350.0 250.0 5 103 104 114 113 103 
+109 383.333333333 250.0 5 104 105 115 114 104 
+110 416.666666667 250.0 5 105 106 116 115 105 
+111 450.0 250.0 5 106 107 117 116 106 
+112 483.333333333 250.0 5 107 108 118 117 107 
+113 216.666666667 216.666666667 5 109 110 122 121 109 
+114 250.0 216.666666667 5 110 111 123 122 110 
+115 283.333333333 216.666666667 5 111 112 124 123 111 
+116 316.666666667 216.666666667 5 112 113 125 124 112 
+117 350.0 216.666666667 5 113 114 126 125 113 
+118 383.333333333 216.666666667 5 114 115 127 126 114 
+119 416.666666667 216.666666667 5 115 116 128 127 115 
+120 450.0 216.666666667 5 116 117 129 128 116 
+121 483.333333333 216.666666667 5 117 118 130 129 117 
\ No newline at end of file
diff --git a/examples/data/mf6/create_tests/test006_2models_diff_dis/exg.txt b/examples/data/mf6/create_tests/test006_2models_diff_dis/exg.txt
new file mode 100644
index 000000000..29122b239
--- /dev/null
+++ b/examples/data/mf6/create_tests/test006_2models_diff_dis/exg.txt
@@ -0,0 +1,9 @@
+     1 3 2  1 1    1  50. 16.67 33.33   100.99 
+     1 3 2  1 2   1  50. 16.67 33.33   100.99
+     1 3 2  1 3    1  50. 16.67 33.33   100.99
+     1 4 2  1 4    1  50. 16.67 33.33   100.99
+     1 4 2  1 5   1  50. 16.67 33.33   100.99
+     1 4 2  1 6    1  50. 16.67 33.33   100.99
+     1 5 2  1 7   1  50. 16.67 33.33   100.99
+     1 5 2  1 8    1  50. 16.67 33.33   100.99
+     1 5 2  1 9    1  50. 16.67 33.33   100.99                             
diff --git a/examples/data/mf6/create_tests/test006_2models_diff_dis/gnc.txt b/examples/data/mf6/create_tests/test006_2models_diff_dis/gnc.txt
new file mode 100644
index 000000000..1e139d0ac
--- /dev/null
+++ b/examples/data/mf6/create_tests/test006_2models_diff_dis/gnc.txt
@@ -0,0 +1,9 @@
+  1 3 2  1 1  1 2 2    0.333333333333
+  1 3 2  1 2  0 0 0    0.333333333333
+  1 3 2  1 3  1 4 2    0.333333333333
+  1 4 2  1 4  1 3 2    0.333333333333
+  1 4 2  1 5  0 0 0    0.333333333333
+  1 4 2  1 6  1 5 2    0.333333333333
+  1 5 2  1 7  1 4 2    0.333333333333
+  1 5 2  1 8  0 0 0    0.333333333333
+  1 5 2  1 9  1 6 2    0.333333333333
diff --git a/examples/data/mf6/create_tests/test006_2models_diff_dis/vertices.txt b/examples/data/mf6/create_tests/test006_2models_diff_dis/vertices.txt
new file mode 100644
index 000000000..7afb23a51
--- /dev/null
+++ b/examples/data/mf6/create_tests/test006_2models_diff_dis/vertices.txt
@@ -0,0 +1,148 @@
+1 0.0 700.0
+2 100.0 700.0
+3 200.0 700.0
+4 300.0 700.0
+5 400.0 700.0
+6 500.0 700.0
+7 600.0 700.0
+8 700.0 700.0
+9 0.0 600.0
+10 100.0 600.0
+11 200.0 600.0
+12 300.0 600.0
+13 400.0 600.0
+14 500.0 600.0
+15 600.0 600.0
+16 700.0 600.0
+17 0.0 500.0
+18 100.0 500.0
+19 200.0 500.0
+20 233.333333333 500.0
+21 266.666666667 500.0
+22 300.0 500.0
+23 333.333333333 500.0
+24 366.666666667 500.0
+25 400.0 500.0
+26 433.333333333 500.0
+27 466.666666667 500.0
+28 500.0 500.0
+29 600.0 500.0
+30 700.0 500.0
+31 200.0 466.666666667
+32 233.333333333 466.666666667
+33 266.666666667 466.666666667
+34 300.0 466.666666667
+35 333.333333333 466.666666667
+36 366.666666667 466.666666667
+37 400.0 466.666666667
+38 433.333333333 466.666666667
+39 466.666666667 466.666666667
+40 500.0 466.666666667
+41 200.0 433.333333333
+42 233.333333333 433.333333333
+43 266.666666667 433.333333333
+44 300.0 433.333333333
+45 333.333333333 433.333333333
+46 366.666666667 433.333333333
+47 400.0 433.333333333
+48 433.333333333 433.333333333
+49 466.666666667 433.333333333
+50 500.0 433.333333333
+51 0.0 400.0
+52 100.0 400.0
+53 200.0 400.0
+54 233.333333333 400.0
+55 266.666666667 400.0
+56 300.0 400.0
+57 333.333333333 400.0
+58 366.666666667 400.0
+59 400.0 400.0
+60 433.333333333 400.0
+61 466.666666667 400.0
+62 500.0 400.0
+63 600.0 400.0
+64 700.0 400.0
+65 200.0 366.666666667
+66 233.333333333 366.666666667
+67 266.666666667 366.666666667
+68 300.0 366.666666667
+69 333.333333333 366.666666667
+70 366.666666667 366.666666667
+71 400.0 366.666666667
+72 433.333333333 366.666666667
+73 466.666666667 366.666666667
+74 500.0 366.666666667
+75 200.0 333.333333333
+76 233.333333333 333.333333333
+77 266.666666667 333.333333333
+78 300.0 333.333333333
+79 333.333333333 333.333333333
+80 366.666666667 333.333333333
+81 400.0 333.333333333
+82 433.333333333 333.333333333
+83 466.666666667 333.333333333
+84 500.0 333.333333333
+85 0.0 300.0
+86 100.0 300.0
+87 200.0 300.0
+88 233.333333333 300.0
+89 266.666666667 300.0
+90 300.0 300.0
+91 333.333333333 300.0
+92 366.666666667 300.0
+93 400.0 300.0
+94 433.333333333 300.0
+95 466.666666667 300.0
+96 500.0 300.0
+97 600.0 300.0
+98 700.0 300.0
+99 200.0 266.666666667
+100 233.333333333 266.666666667
+101 266.666666667 266.666666667
+102 300.0 266.666666667
+103 333.333333333 266.666666667
+104 366.666666667 266.666666667
+105 400.0 266.666666667
+106 433.333333333 266.666666667
+107 466.666666667 266.666666667
+108 500.0 266.666666667
+109 200.0 233.333333333
+110 233.333333333 233.333333333
+111 266.666666667 233.333333333
+112 300.0 233.333333333
+113 333.333333333 233.333333333
+114 366.666666667 233.333333333
+115 400.0 233.333333333
+116 433.333333333 233.333333333
+117 466.666666667 233.333333333
+118 500.0 233.333333333
+119 0.0 200.0
+120 100.0 200.0
+121 200.0 200.0
+122 233.333333333 200.0
+123 266.666666667 200.0
+124 300.0 200.0
+125 333.333333333 200.0
+126 366.666666667 200.0
+127 400.0 200.0
+128 433.333333333 200.0
+129 466.666666667 200.0
+130 500.0 200.0
+131 600.0 200.0
+132 700.0 200.0
+133 0.0 100.0
+134 100.0 100.0
+135 200.0 100.0
+136 300.0 100.0
+137 400.0 100.0
+138 500.0 100.0
+139 600.0 100.0
+140 700.0 100.0
+141 0.0 0.0
+142 100.0 0.0
+143 200.0 0.0
+144 300.0 0.0
+145 400.0 0.0
+146 500.0 0.0
+147 600.0 0.0
+148 700.0 0.0
\ No newline at end of file
diff --git a/flopy/mf6/coordinates/modeldimensions.py b/flopy/mf6/coordinates/modeldimensions.py
index fadb11a12..c41003592 100644
--- a/flopy/mf6/coordinates/modeldimensions.py
+++ b/flopy/mf6/coordinates/modeldimensions.py
@@ -66,15 +66,17 @@ def unlock(self):
         self.locked = False
         self.package_dim.unlock()
 
-    def get_model_grid(self, data_item_num=None):
+    def get_model_grid(self, data_item_num=None, model_num=None):
         if self.locked:
-            if self.model_grid is None:
+            if self.model_grid is None or model_num is not None:
                 self.model_grid = self.get_model_dim(
-                    data_item_num
+                    data_item_num, model_num
                 ).get_model_grid()
             return self.model_grid
         else:
-            return self.get_model_dim(data_item_num).get_model_grid()
+            return self.get_model_dim(
+                data_item_num, model_num
+            ).get_model_grid()
 
     def get_data_shape(
         self,
@@ -100,24 +102,37 @@ def model_subspace_size(self, subspace_string="", data_item_num=None):
             subspace_string
         )
 
-    def get_model_dim(self, data_item_num):
+    def get_model_dim(self, data_item_num, model_num=None):
         if (
             self.package_dim.model_dim is None
-            or data_item_num is None
+            or (data_item_num is None and model_num is None)
             or len(self.package_dim.model_dim) == 1
         ):
             return self.package_dim.model_dim[0]
         else:
-            if not (len(self.structure.data_item_structures) > data_item_num):
-                raise FlopyException(
-                    'Data item index "{}" requested which '
-                    "is greater than the maximum index of"
-                    "{}.".format(
-                        data_item_num,
-                        len(self.structure.data_item_structures) - 1,
+            if model_num is None:
+                model_num = self.structure.data_item_structures[data_item_num][
+                    -1
+                ]
+                if not (
+                    len(self.structure.data_item_structures) > data_item_num
+                ):
+                    raise FlopyException(
+                        'Data item index "{}" requested which '
+                        "is greater than the maximum index of"
+                        "{}.".format(
+                            data_item_num,
+                            len(self.structure.data_item_structures) - 1,
+                        )
                     )
-                )
-            model_num = self.structure.data_item_structures[data_item_num][-1]
+            else:
+                if not len(self.package_dim.model_dim) > model_num:
+                    raise FlopyException(
+                        f'Model item index "{model_num}" requested which '
+                        "is greater than the maximum index of"
+                        f"{len(self.package_dim.model_dim)}."
+                    )
+
             if DatumUtil.is_int(model_num):
                 return self.package_dim.model_dim[int(model_num)]
 
diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py
index c9635b458..7cffa30bc 100644
--- a/flopy/mf6/data/mfdatalist.py
+++ b/flopy/mf6/data/mfdatalist.py
@@ -7,7 +7,7 @@
 
 from ...datbase import DataListInterface, DataType
 from ...mbase import ModelInterface
-from ...utils import datautil
+from ...utils.datautil import DatumUtil
 from ..data import mfdata, mfstructure
 from ..mfbase import ExtFileAction, MFDataException, VerbosityLevel
 from ..utils.mfenums import DiscretizationType
@@ -1043,14 +1043,21 @@ def _get_file_entry_record(
                         data_val = data_line[index]
                         if data_item.is_cellid or (
                             data_item.possible_cellid
-                            and storage._validate_cellid([data_val], 0)
+                            and storage._validate_cellid(
+                                [data_val], 0, data_item
+                            )
                         ):
                             if (
                                 data_item.shape is not None
                                 and len(data_item.shape) > 0
                                 and data_item.shape[0] == "ncelldim"
                             ):
-                                model_grid = data_dim.get_model_grid()
+                                model_num = DatumUtil.cellid_model_num(
+                                    data_item,
+                                    self.structure.model_data,
+                                    self._data_dimensions.package_dim.model_dim,
+                                )
+                                model_grid = data_dim.get_model_grid(model_num)
                                 cellid_size = (
                                     model_grid.get_num_spatial_coordinates()
                                 )
@@ -1058,9 +1065,7 @@ def _get_file_entry_record(
                                     resolved_shape, cellid_size
                                 )
                         data_size = 1
-                        if len(
-                            resolved_shape
-                        ) == 1 and datautil.DatumUtil.is_int(
+                        if len(resolved_shape) == 1 and DatumUtil.is_int(
                             resolved_shape[0]
                         ):
                             data_size = int(resolved_shape[0])
@@ -1722,7 +1727,7 @@ def store_as_external_file(
                 or replace_existing_external
             ):
                 fname, ext = os.path.splitext(external_file_path)
-                if datautil.DatumUtil.is_int(sp):
+                if DatumUtil.is_int(sp):
                     full_name = f"{fname}_{sp + 1}{ext}"
                 else:
                     full_name = f"{fname}_{sp}{ext}"
diff --git a/flopy/mf6/data/mfdatastorage.py b/flopy/mf6/data/mfdatastorage.py
index b235e38ad..491c2e20c 100644
--- a/flopy/mf6/data/mfdatastorage.py
+++ b/flopy/mf6/data/mfdatastorage.py
@@ -1569,25 +1569,35 @@ def _build_recarray(self, data, key, autofill):
             self._verify_list(new_data)
         return new_data
 
+    def _get_cellid_size(self, data_item_name):
+        model_num = DatumUtil.cellid_model_num(
+            data_item_name,
+            self.data_dimensions.structure.model_data,
+            self.data_dimensions.package_dim.model_dim,
+        )
+        model_grid = self.data_dimensions.get_model_grid(model_num=model_num)
+        return model_grid.get_num_spatial_coordinates()
+
     def make_tuple_cellids(self, data):
         # convert cellids from individual layer, row, column fields into
         # tuples (layer, row, column)
-        data_dim = self.data_dimensions
-        model_grid = data_dim.get_model_grid()
-        cellid_size = model_grid.get_num_spatial_coordinates()
-
         new_data = []
         current_cellid = ()
         for line in data:
+            data_idx = 0
             new_line = []
             for item, is_cellid in zip(line, self.recarray_cellid_list_ex):
                 if is_cellid:
+                    cellid_size = self._get_cellid_size(
+                        self._recarray_type_list[data_idx][0],
+                    )
                     current_cellid += (item,)
                     if len(current_cellid) == cellid_size:
                         new_line.append(current_cellid)
                         current_cellid = ()
                 else:
                     new_line.append(item)
+                    data_idx += 1
             new_data.append(tuple(new_line))
         return new_data
 
@@ -2092,14 +2102,14 @@ def resolve_shape_list(
                 return False, True
         return False, False
 
-    def _validate_cellid(self, arr_line, data_index):
+    def _validate_cellid(self, arr_line, data_index, data_item):
         if not self.data_dimensions.structure.model_data:
             # not model data so this is not a cell id
             return False
         if arr_line is None:
             return False
+        cellid_size = self._get_cellid_size(data_item.name)
         model_grid = self.data_dimensions.get_model_grid()
-        cellid_size = model_grid.get_num_spatial_coordinates()
         if cellid_size + data_index > len(arr_line):
             return False
         for index, dim_size in zip(
@@ -2353,9 +2363,8 @@ def _verify_list(self, data):
                         # this is a cell id.  verify that it contains the
                         # correct number of integers
                         if cellid_size is None:
-                            model_grid = datadim.get_model_grid()
-                            cellid_size = (
-                                model_grid.get_num_spatial_coordinates()
+                            cellid_size = self._get_cellid_size(
+                                self._recarray_type_list[index][0]
                             )
                         if (
                             cellid_size != 1
@@ -2934,9 +2943,7 @@ def build_type_list(
                             ):
                                 # A cellid is a single entry (tuple) in the
                                 # recarray.  Adjust dimensions accordingly.
-                                data_dim = self.data_dimensions
-                                grid = data_dim.get_model_grid()
-                                size = grid.get_num_spatial_coordinates()
+                                size = self._get_cellid_size(data_item.name)
                                 data_item.remove_cellid(resolved_shape, size)
                         if not data_item.optional or not min_size:
                             for index in range(0, resolved_shape[0]):
@@ -2973,8 +2980,7 @@ def _append_type_lists(self, name, data_type, iscellid):
         if iscellid and self._model_or_sim.model_type is not None:
             # write each part of the cellid out as a separate entry
             # to _recarray_list_list_ex
-            model_grid = self.data_dimensions.get_model_grid()
-            cellid_size = model_grid.get_num_spatial_coordinates()
+            cellid_size = self._get_cellid_size(name)
             # determine header for different grid types
             if cellid_size == 1:
                 self._do_ex_list_append(name, int, iscellid)
diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py
index dfe847060..1fa241ea0 100644
--- a/flopy/mf6/data/mfdatautil.py
+++ b/flopy/mf6/data/mfdatautil.py
@@ -186,7 +186,12 @@ def to_string(
             and is_cellid
             and data_dim.get_model_dim(None).model_name is not None
         ):
-            model_grid = data_dim.get_model_grid()
+            model_num = DatumUtil.cellid_model_num(
+                data_item.name,
+                data_dim.structure.model_data,
+                data_dim.package_dim.model_dim,
+            )
+            model_grid = data_dim.get_model_grid(model_num=model_num)
             cellid_size = model_grid.get_num_spatial_coordinates()
             if len(val) != cellid_size:
                 message = (
diff --git a/flopy/mf6/data/mffileaccess.py b/flopy/mf6/data/mffileaccess.py
index 5a2238aca..94efd5953 100644
--- a/flopy/mf6/data/mffileaccess.py
+++ b/flopy/mf6/data/mffileaccess.py
@@ -2071,7 +2071,7 @@ def _append_data_list(
             self._last_line_info.append([])
         if data_item.is_cellid or (
             data_item.possible_cellid
-            and storage._validate_cellid(arr_line, data_index)
+            and storage._validate_cellid(arr_line, data_index, data_item)
         ):
             if self._data_dimensions is None:
                 comment = (
@@ -2096,8 +2096,16 @@ def _append_data_list(
                     comment,
                     self._simulation_data.debug,
                 )
+            # in case of multiple model grids, determine which one to use
+            model_num = DatumUtil.cellid_model_num(
+                data_item.name,
+                struct.model_data,
+                self._data_dimensions.package_dim.model_dim,
+            )
+            model_grid = self._data_dimensions.get_model_grid(
+                model_num=model_num
+            )
             # read in the entire cellid
-            model_grid = self._data_dimensions.get_model_grid()
             cellid_size = model_grid.get_num_spatial_coordinates()
             cellid_tuple = ()
             if (
diff --git a/flopy/mf6/modflow/mfsimulation.py b/flopy/mf6/modflow/mfsimulation.py
index 3e9304ffc..1a292b7d6 100644
--- a/flopy/mf6/modflow/mfsimulation.py
+++ b/flopy/mf6/modflow/mfsimulation.py
@@ -885,6 +885,7 @@ def load(
                         f"  loading exchange package {exchange_file._get_pname()}..."
                     )
                 exchange_file.load(strict)
+                instance._add_package(exchange_file, exchange_file.path)
                 instance._exchange_files[exgfile[1]] = exchange_file
 
         # load simulation packages
diff --git a/flopy/mf6/utils/testutils.py b/flopy/mf6/utils/testutils.py
index 9e1e8676c..eeb66bf17 100644
--- a/flopy/mf6/utils/testutils.py
+++ b/flopy/mf6/utils/testutils.py
@@ -34,26 +34,27 @@ def read_cell2d(cell2d_file):
     return c2drecarray
 
 
-def read_exchangedata(gwf_file, cellid_size=3):
+def read_exchangedata(gwf_file, cellid_size=3, cellid_size_2=3):
     exgrecarray = []
     fd = open(gwf_file, "r")
     for line in fd:
         linesp = line.strip().split()
+        cellid_tot = cellid_size + cellid_size_2
         exgrecarray.append(
             (
                 make_int_tuple(linesp[0:cellid_size]),
-                make_int_tuple(linesp[cellid_size : cellid_size * 2]),
-                int(linesp[cellid_size * 2]),
-                float(linesp[cellid_size * 2 + 1]),
-                float(linesp[cellid_size * 2 + 2]),
-                float(linesp[cellid_size * 2 + 3]),
-                float(linesp[cellid_size * 2 + 4]),
+                make_int_tuple(linesp[cellid_size:cellid_tot]),
+                int(linesp[cellid_tot]),
+                float(linesp[cellid_tot + 1]),
+                float(linesp[cellid_tot + 2]),
+                float(linesp[cellid_tot + 3]),
+                float(linesp[cellid_tot + 4]),
             )
         )
     return exgrecarray
 
 
-def read_gncrecarray(gnc_file, cellid_size=3):
+def read_gncrecarray(gnc_file, cellid_size=3, cellid_size_2=3):
     gncrecarray = []
     fd = open(gnc_file, "r")
     for line in fd:
@@ -61,9 +62,17 @@ def read_gncrecarray(gnc_file, cellid_size=3):
         gncrecarray.append(
             (
                 make_int_tuple(linesp[0:cellid_size]),
-                make_int_tuple(linesp[cellid_size : cellid_size * 2]),
-                make_int_tuple(linesp[cellid_size * 2 : cellid_size * 3]),
-                float(linesp[cellid_size * 3]),
+                make_int_tuple(
+                    linesp[cellid_size : cellid_size + cellid_size_2]
+                ),
+                make_int_tuple(
+                    linesp[
+                        cellid_size
+                        + cellid_size_2 : cellid_size * 2
+                        + cellid_size_2
+                    ]
+                ),
+                float(linesp[cellid_size * 2 + cellid_size_2]),
             )
         )
     return gncrecarray
diff --git a/flopy/utils/datautil.py b/flopy/utils/datautil.py
index f0df79f1a..53844ae46 100644
--- a/flopy/utils/datautil.py
+++ b/flopy/utils/datautil.py
@@ -83,6 +83,26 @@ def is_basic_type(obj):
             return True
         return False
 
+    @staticmethod
+    def cellid_model_num(data_item_name, model_data, model_dim):
+        # determine which model to use based on cellid name
+        # contains hard coded relationship between data item names and
+        # model number
+        # TODO: Incorporate this into the DFNs
+        if model_data:
+            return None
+        if data_item_name.startswith("cellidm") and len(data_item_name) > 7:
+            model_num = data_item_name[7:]
+            if DatumUtil.is_int(model_num):
+                return int(model_num) - 1
+        if (
+            data_item_name == "cellidn" or data_item_name == "cellidsj"
+        ) and len(model_dim) > 0:
+            return 0
+        elif data_item_name == "cellidm" and len(model_dim) > 1:
+            return 1
+        return None
+
 
 class PyListUtil:
     """