diff --git a/.docs/Notebooks/gnc_example.py b/.docs/Notebooks/gnc_example.py
new file mode 100644
index 000000000..206e08a36
--- /dev/null
+++ b/.docs/Notebooks/gnc_example.py
@@ -0,0 +1,328 @@
+# ---
+# jupyter:
+#   jupytext:
+#     notebook_metadata_filter: all
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.14.5
+#   kernelspec:
+#     display_name: Python 3 (ipykernel)
+#     language: python
+#     name: python3
+#   metadata:
+#     section: dis
+#     authors:
+#       - name: Joseph Hughes
+# ---
+
+# # Ghost Node Correction (GNC) Data for MODFLOW 6
+#
+# The control volume finite difference formulation used by MODFLOW assumes that the line connecting two cell centers crosses the shared face at a right angle through the middle of the face. A quadtree grid violates that assumption wherever a coarse cell connects to a finer cell, because the shared face is offset from the center of the coarse cell. The Ghost Node Correction (GNC) Package corrects the resulting error by interpolating the head at a ghost node, which is the point in the coarse cell that does lie on the perpendicular through the middle of the face.
+#
+# FloPy builds GNC Package input two ways, and we demonstrate both here. The first uses the ghost node data GRIDGEN writes when it exports a grid. The second computes the ghost node data from a model grid, a grid conforming array of refinement levels, and the grid connectivity, and does not require GRIDGEN.
+#
+# We also compare the ghost node correction against XT3D, which is the other MODFLOW 6 option for improving accuracy on a quadtree grid, in terms of both the answer and what the correction costs.
+
+# +
+import re
+import sys
+from pathlib import Path
+from tempfile import TemporaryDirectory
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+from shapely.geometry import Polygon
+
+import flopy
+from flopy.utils import flopy_io, get_gnc, get_gridprops_gnc6
+from flopy.utils.gridgen import Gridgen
+
+print(sys.version)
+print(f"numpy version: {np.__version__}")
+print(f"matplotlib version: {mpl.__version__}")
+print(f"flopy version: {flopy.__version__}")
+# -
+
+# The FloPy GRIDGEN module requires that the gridgen executable can be called using subprocess **(i.e., gridgen is in your path)**.
+
+gridgen_exe = flopy.which("gridgen")
+if gridgen_exe is None:
+    msg = (
+        "Warning, gridgen is not in your path. "
+        "When you create the gridgen object you will need to "
+        "provide a full path to the gridgen binary executable."
+    )
+    print(msg)
+else:
+    print(f"gridgen executable was found at: {flopy_io.relpath_safe(gridgen_exe)}")
+
+# +
+temp_dir = TemporaryDirectory()
+workspace = Path(temp_dir.name)
+gridgen_ws = workspace / "gridgen"
+gridgen_ws.mkdir(parents=True, exist_ok=True)
+
+print(f"Model workspace is : {flopy_io.scrub_login(str(workspace))}")
+print(f"Gridgen workspace is : {flopy_io.scrub_login(str(gridgen_ws))}")
+# -
+
+# ## Build the quadtree grid
+#
+# GRIDGEN works from a base MODFLOW grid. We use a 3 layer grid of 20 rows and 20 columns and refine a square in the middle of the grid by three levels, which produces cells one eighth the width of the base grid cells.
+
+# +
+nlay, nrow, ncol = 3, 20, 20
+delr = delc = 1.0
+top = 1.0
+botm = [top - (k + 1) * top / nlay for k in range(nlay)]
+
+base_grid = flopy.discretization.StructuredGrid(
+    delr=np.full(ncol, delr, dtype=float),
+    delc=np.full(nrow, delc, dtype=float),
+    top=np.full((nrow, ncol), top),
+    botm=np.array([np.full((nrow, ncol), b) for b in botm]),
+)
+# -
+
+# +
+center, half_width = ncol / 2.0, 3.0
+corners = [
+    (center - half_width, center - half_width),
+    (center + half_width, center - half_width),
+    (center + half_width, center + half_width),
+    (center - half_width, center + half_width),
+]
+
+g = Gridgen(base_grid, model_ws=str(gridgen_ws))
+g.add_refinement_features([Polygon(corners)], "polygon", 3, range(nlay))
+g.build(verbose=False)
+
+disv_gridprops = g.get_gridprops_disv()
+ncpl = disv_gridprops["ncpl"]
+print(f"Number of cells per layer: {ncpl}")
+# -
+
+# ## Ghost node data from GRIDGEN
+#
+# GRIDGEN computes the ghost node data whenever it exports a grid and writes it to the `qtg.gnc.dat` file. The `get_gnc()` method reads that file and returns a record array with zero-based node numbers, where cell `n` contains the ghost node, cell `m` is the connecting cell, and cells `j0` and `j1` are the contributing cells whose heads are interpolated.
+
+gnc = g.get_gnc()
+print(f"Number of ghost nodes: {len(gnc)}")
+print(gnc[:5])
+
+# GRIDGEN always writes two contributing cells. When a ghost node has only one contributing cell, that cell is repeated and its contributing factor is halved, which MODFLOW accumulates into the same matrix position. The contributing factors always sum to less than one, because one minus the sum is the factor applied to the head in cell `n`.
+
+alpha = gnc["alpha0"] + gnc["alpha1"]
+print(f"Contributing factors range from {alpha.min():.4f} to {alpha.max():.4f}")
+
+# A grid that is not refined has no ghost nodes and the record array is empty. The GNC Package should not be created in that case.
+
+# The `get_gridprops_gnc6()` method converts the node numbers to cellids and returns a dictionary that can be unpacked directly into the `ModflowGwfgnc` constructor. Cellids are built for a DISV grid here; pass `dis_type="disu"` for a DISU grid.
+
+gnc_gridprops = g.get_gridprops_gnc6(dis_type="disv")
+print(f"numgnc: {gnc_gridprops['numgnc']}")
+print(f"numalphaj: {gnc_gridprops['numalphaj']}")
+print(f"first record: {gnc_gridprops['gncdata'][0]}")
+
+# ## Ghost node data from the model grid
+#
+# The ghost node data can also be computed from the model grid, without running GRIDGEN. All that is needed is the grid, which provides the cell centers and the connectivity, and a grid conforming array of refinement levels. Cell areas are used when levels are not supplied.
+
+vgrid = flopy.discretization.VertexGrid(**g.get_gridprops_vertexgrid())
+
+# A vertex grid does not carry connectivity, so `get_gnc()` builds it from the cells that share an edge. Connectivity can also be passed with the `ia` or `iac` and `ja` arguments, which is what an unstructured grid already provides.
+#
+# The refinement level of each cell follows from the cell area, where level 0 is a base grid cell and each level halves the cell width. We compute the areas from the cell vertices so that nothing in this section depends on GRIDGEN. The levels are given for one layer, and `get_gnc()` applies them to every layer.
+
+
+# +
+def cell_area(icpl):
+    x, y = np.array(vgrid.get_cell_vertices(icpl)).T
+    return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
+
+
+area = np.array([cell_area(icpl) for icpl in range(ncpl)])
+level = np.round(np.log2(np.sqrt(area.max() / area))).astype(int)
+print(f"Refinement levels present: {np.unique(level)}")
+# -
+
+# The `get_gnc()` function returns the same record array that GRIDGEN wrote. We ask for two contributing cells so the result can be compared directly.
+
+gnc_grid = get_gnc(vgrid, level=level, numalphaj=2)
+print(f"Number of ghost nodes: {len(gnc_grid)}")
+print(gnc_grid[:5])
+
+
+# The two record arrays hold the same ghost nodes. We sort the records because the two routines visit the cells in a different order, and we sort the contributing cells within each record because the two cells are sometimes listed in the opposite order. That ordering does not matter, since MODFLOW accumulates the contribution of each cell.
+
+
+# +
+def sort_gnc(recarray):
+    nodes = np.sort(np.column_stack([recarray["j0"], recarray["j1"]]), axis=1)
+    alpha = np.sort(np.column_stack([recarray["alpha0"], recarray["alpha1"]]), axis=1)
+    key = np.column_stack([recarray["n"], recarray["m"], nodes, alpha])
+    return key[np.lexsort(key.T[::-1])]
+
+
+# GRIDGEN writes the contributing factors with six significant digits
+assert np.allclose(sort_gnc(gnc_grid), sort_gnc(gnc), atol=2e-6)
+print("The computed ghost node data matches the GRIDGEN ghost node data.")
+# -
+
+# The dictionary for the GNC Package is built with the `get_gridprops_gnc6()` function, which also verifies that each cell `n` is connected to cell `m` and that the contributing factors sum to less than one.
+
+gnc_gridprops = get_gridprops_gnc6(gnc_grid, dis_type="disv", ncpl=ncpl)
+print(f"numgnc: {gnc_gridprops['numgnc']}")
+
+# ## Where the ghost nodes are
+#
+# Every ghost node lies on a connection between a coarse cell and a finer cell, so the ghost nodes trace the boundary of the refined area. We plot the connections in the upper layer.
+
+# +
+fig, ax = plt.subplots(figsize=(7, 7))
+ax.set_aspect("equal")
+pmv = flopy.plot.PlotMapView(modelgrid=vgrid, ax=ax, layer=0)
+pmv.plot_grid(colors="0.5", lw=0.5)
+
+xc, yc = vgrid.xcellcenters, vgrid.ycellcenters
+for rec in gnc_grid[gnc_grid["n"] < ncpl]:
+    n, m = rec["n"], rec["m"]
+    ax.plot([xc[n], xc[m]], [yc[n], yc[m]], color="C3", lw=1.0, zorder=2)
+    ax.plot(xc[m], yc[m], "o", color="C3", ms=2.5, zorder=3)
+    for j in (rec["j0"], rec["j1"]):
+        ax.plot(xc[j], yc[j], "s", color="C0", ms=3.0, zorder=3)
+
+ax.plot([], [], color="C3", lw=1.0, label="ghost node connection")
+ax.plot([], [], "s", color="C0", ms=3.0, lw=0, label="contributing cell")
+ax.legend(loc="upper right", framealpha=1.0)
+ax.set_title("Ghost node connections in layer 1")
+# -
+
+# ## Effect of the correction
+#
+# We build the same model three ways and compare the results. The uncorrected model uses the standard formulation, the corrected model adds the GNC Package, and the third model uses XT3D. The correction is applied implicitly by default, so the BICGSTAB linear acceleration option is specified in the IMS Package.
+#
+# The model is confined and homogeneous, with constant heads on the left and right edges and no flow across the top and bottom edges. Head then varies linearly between the two constant head columns, which gives an exact solution to compare against.
+
+# +
+h_left, h_right = 1.0, 0.0
+xcenters = vgrid.xcellcenters
+left = [icpl for icpl in range(ncpl) if xcenters[icpl] < delr]
+right = [icpl for icpl in range(ncpl) if xcenters[icpl] > ncol - delr]
+
+chdspd = [[(k, icpl), h_left] for k in range(nlay) for icpl in left]
+chdspd += [[(k, icpl), h_right] for k in range(nlay) for icpl in right]
+
+x_left, x_right = xcenters[left].mean(), xcenters[right].mean()
+exact = h_left + (h_right - h_left) * (xcenters - x_left) / (x_right - x_left)
+exact = np.tile(exact, nlay)
+print(f"Number of constant head cells: {len(chdspd)}")
+# -
+
+
+# MODFLOW 6 reports the memory it allocates at the end of the simulation listing file, which we read back for each model along with the simulated heads.
+
+
+# +
+def run_model(name, gnc=False, xt3d=False):
+    ws = workspace / name
+    sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=str(ws), exe_name="mf6")
+    flopy.mf6.ModflowTdis(sim)
+    flopy.mf6.ModflowIms(
+        sim,
+        linear_acceleration="bicgstab",
+        inner_maximum=1000,
+        inner_dvclose=1e-10,
+        outer_dvclose=1e-10,
+    )
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name)
+    flopy.mf6.ModflowGwfdisv(gwf, **disv_gridprops)
+    flopy.mf6.ModflowGwfic(gwf, strt=0.5 * (h_left + h_right))
+    flopy.mf6.ModflowGwfnpf(gwf, xt3doptions=xt3d, icelltype=0, k=1.0)
+    flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdspd)
+    flopy.mf6.ModflowGwfoc(
+        gwf, head_filerecord=f"{name}.hds", saverecord=[("HEAD", "ALL")]
+    )
+    if gnc:
+        flopy.mf6.ModflowGwfgnc(gwf, **gnc_gridprops)
+    sim.write_simulation(silent=True)
+    success, buff = sim.run_simulation(silent=True)
+    assert success, f"{name} did not converge"
+
+    listing = (ws / "mfsim.lst").open().read()
+    memory = float(re.search(r"Total\s+([0-9.E+-]+)\s*\n\s*Virtual", listing).group(1))
+    return gwf.output.head().get_data().flatten(), memory
+
+
+# +
+heads, error, memory = {}, {}, {}
+for name, kwargs in [
+    ("uncorrected", {}),
+    ("gnc", {"gnc": True}),
+    ("xt3d", {"xt3d": True}),
+]:
+    heads[name], memory[name] = run_model(name, **kwargs)
+    error[name] = np.abs(heads[name] - exact)
+
+print(f"{'variant':14s}{'max error':>12s}{'rms error':>12s}{'memory, MB':>13s}")
+for name in ("uncorrected", "gnc", "xt3d"):
+    rms = np.sqrt((error[name] ** 2).mean())
+    print(f"{name:14s}{error[name].max():12.3e}{rms:12.3e}{memory[name]:13.1f}")
+# -
+
+# The ghost node correction removes about 18 times the head error introduced by the refinement. XT3D reproduces a linear head field exactly by construction, so it is exact on this problem; that is a property of this test rather than a general ranking of the two corrections.
+
+for name in ("gnc", "xt3d"):
+    removed = 1.0 - error[name].max() / error["uncorrected"].max()
+    print(f"{name:5s} removes {100 * removed:.1f} percent of the error")
+
+# ## Cost of the correction
+#
+# The two corrections reach a comparable answer by different means. XT3D replaces the flow calculation on every connection in the model, which extends the stencil of every cell. The ghost node correction only adds terms on the connections that have a ghost node, which are the connections between a coarse cell and a finer cell, and there are far fewer of those.
+
+print(f"Cells in the model: {ncpl * nlay}")
+print(f"Ghost nodes: {gnc_gridprops['numgnc']}")
+print(
+    f"Ghost nodes are on {100 * gnc_gridprops['numgnc'] / (ncpl * nlay):.1f} "
+    "percent of the cells"
+)
+
+# That shows up in the memory MODFLOW 6 allocates. XT3D nearly doubles it, because the extended stencil applies to every cell in the model. The ghost node correction adds a couple of percent. Both corrections remove nearly all of the error introduced by the refinement, and the ghost node correction does so in about half the memory.
+
+for name in ("gnc", "xt3d"):
+    print(
+        f"{name:5s} memory relative to the uncorrected model: "
+        f"{memory[name] / memory['uncorrected']:.3f}"
+    )
+print(f"gnc memory relative to xt3d: {memory['gnc'] / memory['xt3d']:.3f}")
+
+# Run times are not compared here. They depend on how many iterations the solver takes, and the ordering of the two corrections changes with the problem, so run time is not a reliable way to choose between them.
+#
+# The practical difference is in what each one asks of the user. XT3D is a single keyword in the NPF Package and needs no other input. The ghost node correction needs the ghost node data, which was the difficult part of using the GNC Package and is what the FloPy functionality shown in this notebook provides.
+
+# ## Where the error is
+#
+# The error in the uncorrected model is concentrated on the boundary of the refined area, which is where the ghost nodes are. The ghost node correction removes most of it.
+
+# +
+vmax = error["uncorrected"].max()
+
+fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)
+for ax, name in zip(axes, ("uncorrected", "gnc")):
+    ax.set_aspect("equal")
+    pmv = flopy.plot.PlotMapView(modelgrid=vgrid, ax=ax, layer=0)
+    cb = pmv.plot_array(error[name], cmap="magma_r", vmin=0.0, vmax=vmax)
+    pmv.plot_grid(colors="0.5", lw=0.3, alpha=0.5)
+    ax.set_title(f"{name}, layer 1")
+fig.colorbar(cb, ax=axes, shrink=0.7, label="absolute head error")
+# -
+
+# Clean up the temporary workspace.
+
+try:
+    temp_dir.cleanup()
+except (PermissionError, NotADirectoryError):
+    pass
diff --git a/.docs/Notebooks/mfusg_gnc_example.py b/.docs/Notebooks/mfusg_gnc_example.py
new file mode 100644
index 000000000..3f4c4d7aa
--- /dev/null
+++ b/.docs/Notebooks/mfusg_gnc_example.py
@@ -0,0 +1,202 @@
+# ---
+# jupyter:
+#   jupytext:
+#     notebook_metadata_filter: all
+#     text_representation:
+#       extension: .py
+#       format_name: light
+#       format_version: '1.5'
+#       jupytext_version: 1.14.5
+#   kernelspec:
+#     display_name: Python 3 (ipykernel)
+#     language: python
+#     name: python3
+#   metadata:
+#     section: mfusg
+#     authors:
+#       - name: Joseph Hughes
+# ---
+
+# # MODFLOW-USG: Ghost Node Correction (GNC) Data for a Quadtree Grid
+#
+# The control volume finite difference formulation used by MODFLOW-USG assumes that the line connecting two cell centers crosses the shared face at a right angle through the middle of the face. A quadtree grid violates that assumption wherever a coarse cell connects to a finer cell, because the shared face is offset from the center of the coarse cell. The Ghost Node Correction (GNC) Package corrects the resulting error by interpolating the head at a ghost node, which is the point in the coarse cell that does lie on the perpendicular through the middle of the face.
+#
+# GRIDGEN computes the ghost node data for a quadtree grid, and FloPy converts it to GNC Package input. We build a quadtree grid, create the GNC Package, and compare the corrected and uncorrected solutions.
+#
+# The same ghost node data can be computed from a model grid without running GRIDGEN, which is shown in the [MODFLOW 6 ghost node correction example](https://flopy.readthedocs.io/en/latest/Notebooks/gnc_example.html).
+
+# +
+import sys
+from pathlib import Path
+from tempfile import TemporaryDirectory
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+from shapely.geometry import Polygon
+
+import flopy
+from flopy.utils import flopy_io
+from flopy.utils.gridgen import Gridgen
+
+print(sys.version)
+print(f"numpy version: {np.__version__}")
+print(f"matplotlib version: {mpl.__version__}")
+print(f"flopy version: {flopy.__version__}")
+# -
+
+# The FloPy GRIDGEN module requires that the gridgen executable can be called using subprocess **(i.e., gridgen is in your path)**.
+
+gridgen_exe = flopy.which("gridgen")
+if gridgen_exe is None:
+    msg = (
+        "Warning, gridgen is not in your path. "
+        "When you create the gridgen object you will need to "
+        "provide a full path to the gridgen binary executable."
+    )
+    print(msg)
+else:
+    print(f"gridgen executable was found at: {flopy_io.relpath_safe(gridgen_exe)}")
+
+# +
+temp_dir = TemporaryDirectory()
+workspace = Path(temp_dir.name)
+gridgen_ws = workspace / "gridgen"
+gridgen_ws.mkdir(parents=True, exist_ok=True)
+
+print(f"Model workspace is : {flopy_io.scrub_login(str(workspace))}")
+print(f"Gridgen workspace is : {flopy_io.scrub_login(str(gridgen_ws))}")
+# -
+
+# ## Build the quadtree grid
+#
+# GRIDGEN works from a base MODFLOW grid. We use a single layer grid of 10 rows and 10 columns and refine a square in the middle of the grid by three levels, which produces cells one eighth the width of the base grid cells.
+
+# +
+nlay, nrow, ncol = 1, 10, 10
+delr = delc = 1.0
+
+base_grid = flopy.discretization.StructuredGrid(
+    delr=np.full(ncol, delr, dtype=float),
+    delc=np.full(nrow, delc, dtype=float),
+    top=np.full((nrow, ncol), 1.0),
+    botm=np.zeros((nlay, nrow, ncol)),
+)
+
+g = Gridgen(base_grid, model_ws=str(gridgen_ws))
+refinement = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
+g.add_refinement_features(refinement, "polygon", 3, layers=[0])
+g.build(verbose=False)
+
+disu_gridprops = g.get_gridprops_disu5()
+print(f"Number of cells: {g.get_nodes()}")
+# -
+
+# ## Ghost node data
+#
+# GRIDGEN computes the ghost node data whenever it exports a grid and writes it to the `qtg.gnc.dat` file. The `get_gnc()` method reads that file and returns a record array with zero-based node numbers, where cell `n` contains the ghost node, cell `m` is the connecting cell, and cells `j0` and `j1` are the contributing cells whose heads are interpolated.
+
+gnc = g.get_gnc()
+print(f"Number of ghost nodes: {len(gnc)}")
+print(gnc[:5])
+
+# Two contributing cells are always written. When a ghost node has only one contributing cell, that cell is repeated and its contributing factor is halved, which MODFLOW-USG accumulates into the same matrix position. MODFLOW-USG reads a fixed number of contributing cells per record and indexes `IBOUND` with each of them, so an unused slot cannot be filled with a dummy cell number of zero the way it can in MODFLOW 6. Repeating a cell keeps every slot valid.
+
+single = gnc["j0"] == gnc["j1"]
+print(f"{single.sum()} of {len(gnc)} ghost nodes have one contributing cell")
+
+# The contributing factors always sum to less than one, because one minus the sum is the factor applied to the head in cell `n`.
+
+alpha = gnc["alpha0"] + gnc["alpha1"]
+print(f"Contributing factors range from {alpha.min():.4f} to {alpha.max():.4f}")
+
+# The `get_gridprops_gnc5()` method returns a dictionary that can be unpacked directly into the `MfUsgGnc` constructor. GRIDGEN writes contributing factors rather than conductances, so `iflalphan` is always 0. The `i2kn` and `isymgncn` options can be set through the method.
+
+gnc_gridprops = g.get_gridprops_gnc5()
+for key in ("numgnc", "numalphaj", "i2kn", "isymgncn", "iflalphan"):
+    print(f"{key}: {gnc_gridprops[key]}")
+
+# ## Build and run the models
+#
+# We build the same model with and without the GNC Package. The default `isymgncn` of 0 updates the left-hand side matrix, which makes the matrix asymmetric, so the model is solved with the complex option of the SMS Package.
+
+# +
+chdspd = []
+for x, y, head in [(0.0, 10.0, 1.0), (10.0, 0.0, 0.0)]:
+    node = g.intersect([(x, y)], "point", 0)["nodenumber"][0]
+    chdspd.append([node, head, head])
+print(f"Constant head cells: {chdspd}")
+
+
+def build_model(name, gnc=False):
+    m = flopy.mfusg.MfUsg(
+        modelname=name,
+        model_ws=str(workspace / name),
+        exe_name="mfusg",
+        structured=False,
+    )
+    flopy.mfusg.MfUsgDisU(m, **disu_gridprops)
+    flopy.mfusg.MfUsgBas(m)
+    flopy.mfusg.MfUsgLpf(m)
+    flopy.modflow.ModflowChd(m, stress_period_data=chdspd)
+    flopy.mfusg.MfUsgSms(m, options="COMPLEX")
+    flopy.modflow.ModflowOc(m, stress_period_data={(0, 0): ["save head"]})
+    if gnc:
+        flopy.mfusg.MfUsgGnc(m, **gnc_gridprops)
+    return m
+
+
+# +
+heads = {}
+for name, gnc_flag in [("uncorrected", False), ("gnc", True)]:
+    m = build_model(name, gnc=gnc_flag)
+    m.write_input()
+    success, buff = m.run_model(silent=True)
+    assert success, f"{name} did not converge"
+    head_file = workspace / name / f"{name}.hds"
+    heads[name] = np.concatenate(flopy.utils.HeadUFile(head_file).get_data())
+    print(f"{name} converged")
+# -
+
+# The GNC Package file lists the cell containing the ghost node, the connecting cell, the two contributing cells, and the two contributing factors, using one-based node numbers.
+
+gnc_file = workspace / "gnc" / "gnc.gnc"
+print("".join(gnc_file.open().readlines()[:8]))
+
+# ## Effect of the correction
+#
+# The correction changes the simulated heads around the refined area, where the ghost nodes are.
+
+# +
+diff = heads["gnc"] - heads["uncorrected"]
+print(f"Maximum head difference: {np.abs(diff).max():.3e}")
+
+ugrid = flopy.discretization.UnstructuredGrid(**g.get_gridprops_unstructuredgrid())
+vmax = np.abs(diff).max()
+
+fig, axes = plt.subplots(1, 2, figsize=(11, 5), constrained_layout=True)
+
+ax = axes[0]
+ax.set_aspect("equal")
+pmv = flopy.plot.PlotMapView(modelgrid=ugrid, ax=ax, layer=0)
+cb = pmv.plot_array(heads["gnc"], cmap="jet")
+pmv.plot_grid(colors="0.5", lw=0.3, alpha=0.5)
+pmv.contour_array(heads["gnc"], levels=[0.2, 0.4, 0.6, 0.8], colors="white")
+ax.set_title("Corrected head")
+fig.colorbar(cb, ax=ax, shrink=0.7, label="head")
+
+ax = axes[1]
+ax.set_aspect("equal")
+pmv = flopy.plot.PlotMapView(modelgrid=ugrid, ax=ax, layer=0)
+cb = pmv.plot_array(diff, cmap="RdBu_r", vmin=-vmax, vmax=vmax)
+pmv.plot_grid(colors="0.5", lw=0.3, alpha=0.5)
+ax.set_title("Corrected minus uncorrected head")
+fig.colorbar(cb, ax=ax, shrink=0.7, label="head difference")
+# -
+
+# Clean up the temporary workspace.
+
+try:
+    temp_dir.cleanup()
+except (PermissionError, NotADirectoryError):
+    pass
diff --git a/.gitignore b/.gitignore
index ace40e661..5b252fbb2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -113,3 +113,7 @@ flopy/mf6/data/toml/
 
 # uv lockfile
 uv.lock
+
+# MODFLOW cell-by-cell budget output
+*.cbb
+*.CBB
diff --git a/autotest/test_gnc.py b/autotest/test_gnc.py
index 9b3b05157..db06980c4 100644
--- a/autotest/test_gnc.py
+++ b/autotest/test_gnc.py
@@ -6,6 +6,8 @@
   - gridgen   : the computed data must reproduce gridgen's qtg.gnc.dat.
 """
 
+import io
+
 import numpy as np
 import pytest
 from modflow_devtools.markers import requires_exe, requires_pkg
@@ -16,10 +18,11 @@
     _check_gnc,
     get_gnc,
     get_gnc_dtype,
+    get_gridprops_gnc5,
     get_gridprops_gnc6,
     get_numalphaj,
 )
-from flopy.utils.gridgen import Gridgen
+from flopy.utils.gridgen import Gridgen, get_ia_from_iac
 
 
 def synthetic_grid(connectivity=True):
@@ -337,3 +340,193 @@ def test_get_gnc_connectivity_needs_constant_ncpl():
     assert grid.iac is None
     with pytest.raises(ValueError, match="different"):
         get_gnc(grid)
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_get_gnc_inputs_agree(function_tmpdir):
+    g = build_gridgen(function_tmpdir)
+    expected = gnc_key(g.get_gnc())
+
+    grid = UnstructuredGrid(**g.get_gridprops_unstructuredgrid())
+    iac = g.get_iac()
+    ia, ja = get_ia_from_iac(iac), g.get_ja(iac.sum())
+
+    # refinement level of every cell, where level 0 is the base grid cell
+    area = g.get_area()
+    level = np.round(np.log2(np.sqrt(area.max() / area))).astype(int)
+    assert level.max() > 0
+
+    ncpl = g.get_gridprops_disv()["ncpl"]
+    vertex_grid = VertexGrid(**g.get_gridprops_vertexgrid())
+
+    for tag, gnc in [
+        ("level", get_gnc(grid, level=level, numalphaj=2)),
+        ("level per layer", get_gnc(grid, level=level[:ncpl], numalphaj=2)),
+        ("vertex grid", get_gnc(vertex_grid, ia=ia, ja=ja, numalphaj=2)),
+        ("iac", get_gnc(grid, iac=iac, ja=ja, numalphaj=2)),
+        # connectivity built from the grid, either the iac and ja an
+        # unstructured grid carries or the cells that share an edge
+        ("unstructured grid only", get_gnc(grid, numalphaj=2)),
+        ("vertex grid only", get_gnc(vertex_grid, numalphaj=2)),
+    ]:
+        assert np.allclose(gnc_key(gnc), expected, atol=2.0e-6), tag
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_get_gridprops_gnc_matches_gridgen(function_tmpdir):
+    g = build_gridgen(function_tmpdir, nlay=1)
+    iac = g.get_iac()
+    ia, ja = get_ia_from_iac(iac), g.get_ja(iac.sum())
+    ncpl = g.get_gridprops_disv()["ncpl"]
+
+    grid = UnstructuredGrid(**g.get_gridprops_unstructuredgrid())
+    gnc = get_gnc(grid, numalphaj=2)
+
+    gridprops = get_gridprops_gnc6(gnc, dis_type="disv", ncpl=ncpl, ia=ia, ja=ja)
+    expected = g.get_gridprops_gnc6(dis_type="disv")
+    assert gridprops["numgnc"] == expected["numgnc"]
+    assert gridprops["numalphaj"] == expected["numalphaj"]
+
+    gridprops = get_gridprops_gnc5(gnc, ia=ia, ja=ja)
+    expected = g.get_gridprops_gnc5()
+    assert gridprops["numgnc"] == expected["numgnc"]
+    assert gridprops["iflalphan"] == 0
+    assert gridprops["gncdata"].dtype == expected["gncdata"].dtype
+
+
+@pytest.mark.slow
+@requires_exe("mf6", "gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_mf6disv_gnc_padding(function_tmpdir):
+    """Repeating a contributing cell must not change the solution"""
+    g = build_gridgen(function_tmpdir, nlay=1)
+    disv_gridprops = g.get_gridprops_disv()
+    iac = g.get_iac()
+    ia, ja = get_ia_from_iac(iac), g.get_ja(iac.sum())
+    grid = UnstructuredGrid(**g.get_gridprops_unstructuredgrid())
+
+    chdspd = []
+    for x, y, head in [(0, 10, 1.0), (10, 0, 0.0)]:
+        ic = g.intersect([(x, y)], "point", 0)["nodenumber"][0]
+        chdspd.append([(0, ic), head])
+
+    def run(numalphaj):
+        gnc = get_gnc(grid, numalphaj=numalphaj)
+        assert get_numalphaj(gnc) == numalphaj
+        gridprops = get_gridprops_gnc6(
+            gnc, dis_type="disv", ncpl=disv_gridprops["ncpl"], ia=ia, ja=ja
+        )
+        ws = function_tmpdir / f"j{numalphaj}"
+        sim = flopy.mf6.MFSimulation(sim_name="m", sim_ws=ws, exe_name="mf6")
+        flopy.mf6.ModflowTdis(sim)
+        flopy.mf6.ModflowIms(
+            sim,
+            linear_acceleration="bicgstab",
+            inner_dvclose=1e-11,
+            outer_dvclose=1e-11,
+        )
+        gwf = flopy.mf6.ModflowGwf(sim, modelname="m")
+        flopy.mf6.ModflowGwfdisv(gwf, **disv_gridprops)
+        flopy.mf6.ModflowGwfic(gwf)
+        flopy.mf6.ModflowGwfnpf(gwf)
+        flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdspd)
+        flopy.mf6.ModflowGwfoc(
+            gwf, head_filerecord="m.hds", saverecord=[("HEAD", "ALL")]
+        )
+        flopy.mf6.ModflowGwfgnc(gwf, **gridprops)
+        sim.write_simulation()
+        success, buff = sim.run_simulation(silent=True)
+        assert success, "\n".join(buff[-25:])
+        return gwf.output.head().get_data().flatten()
+
+    assert np.allclose(run(2), run(4), atol=1e-8)
+
+
+@pytest.mark.slow
+@requires_exe("mfusg", "gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_mfusg_gnc_padding(function_tmpdir):
+    """Repeating a contributing cell must not change the solution"""
+    g = build_gridgen(function_tmpdir, nlay=1)
+    disu_gridprops = g.get_gridprops_disu5()
+    iac = g.get_iac()
+    ia, ja = get_ia_from_iac(iac), g.get_ja(iac.sum())
+    grid = UnstructuredGrid(**g.get_gridprops_unstructuredgrid())
+
+    chdspd = []
+    for x, y, head in [(0, 10, 1.0), (10, 0, 0.0)]:
+        ic = g.intersect([(x, y)], "point", 0)["nodenumber"][0]
+        chdspd.append([ic, head, head])
+
+    def run(numalphaj):
+        gridprops = get_gridprops_gnc5(get_gnc(grid, numalphaj=numalphaj), ia=ia, ja=ja)
+        ws = function_tmpdir / f"j{numalphaj}"
+        m = flopy.mfusg.MfUsg(
+            modelname="m", model_ws=ws, exe_name="mfusg", structured=False
+        )
+        flopy.mfusg.MfUsgDisU(m, **disu_gridprops)
+        flopy.mfusg.MfUsgBas(m)
+        flopy.mfusg.MfUsgLpf(m)
+        flopy.modflow.ModflowChd(m, stress_period_data=chdspd)
+        flopy.mfusg.MfUsgSms(m, options="COMPLEX")
+        flopy.modflow.ModflowOc(m, stress_period_data={(0, 0): ["save head"]})
+        flopy.mfusg.MfUsgGnc(m, **gridprops)
+        m.write_input()
+        success, buff = m.run_model(silent=True)
+        assert success, "\n".join(buff[-25:])
+        return np.concatenate(flopy.utils.HeadUFile(ws / "m.hds").get_data())
+
+    assert np.allclose(run(2), run(4), atol=1e-8)
+
+
+def test_mfusg_gnc_file_fields_stay_separated(function_tmpdir):
+    """A value that fills its format width must not run into the next field"""
+    model = flopy.mfusg.MfUsg(modelname="m", model_ws=function_tmpdir, structured=False)
+    dtype = flopy.mfusg.MfUsgGnc.get_default_dtype(2, 0)
+    gncdata = np.zeros(2, dtype=dtype)
+    gncdata[0] = (23, 33, 22, 22, 0.125, 0.166667)
+    # ten digit node numbers fill the %10d field width
+    gncdata[1] = (1234567889, 1234567889, 1234567889, 1234567889, 0.125, 0.125)
+    flopy.mfusg.MfUsgGnc(model, numgnc=2, numalphaj=2, gncdata=gncdata)
+    model.write_input()
+
+    # the list is read with URWORD, so every record must have one token per field
+    records = (function_tmpdir / "m.gnc").open().readlines()[2:]
+    for line in records:
+        if line.strip():
+            assert len(line.split()) == 6
+
+    # the contributing factors must not be truncated
+    assert np.allclose(float(records[0].split()[5]), 0.166667, atol=1e-6)
+
+
+def test_fmt_string_separates_free_format_fields():
+    """A free format list is read with URWORD, so its fields must be separated
+
+    A value that fills its format width runs into the next value when the
+    field formats are concatenated, which made a record unreadable.  Ten digit
+    node numbers fill the %10d field width.
+    """
+    from flopy.mfusg.cln_dtypes import MfUsgClnDtypes
+    from flopy.mfusg.mfusg import fmt_string
+
+    dtypes = {
+        "gnc": flopy.mfusg.MfUsgGnc.get_default_dtype(2, 0),
+        "cln node": MfUsgClnDtypes.get_clnnode_dtype(),
+    }
+    for name, dtype in dtypes.items():
+        record = np.zeros(1, dtype=dtype)
+        for field in dtype.names:
+            if np.issubdtype(dtype[field], np.integer):
+                record[0][field] = 1234567889
+
+        buff = io.StringIO()
+        np.savetxt(buff, record, fmt=fmt_string(record, free=True), delimiter="")
+        assert len(buff.getvalue().split()) == len(dtype.names), name
+
+        # a fixed format list is read by position, so it stays unseparated
+        buff = io.StringIO()
+        np.savetxt(buff, record, fmt=fmt_string(record, free=False), delimiter="")
+        assert len(buff.getvalue().split()) < len(dtype.names), name
diff --git a/autotest/test_gridgen.py b/autotest/test_gridgen.py
index d369a21cd..27c997aec 100644
--- a/autotest/test_gridgen.py
+++ b/autotest/test_gridgen.py
@@ -14,7 +14,7 @@
 from autotest.test_grid_cases import GridCases
 from flopy.discretization.unstructuredgrid import UnstructuredGrid
 from flopy.discretization.vertexgrid import VertexGrid
-from flopy.utils.gridgen import Gridgen
+from flopy.utils.gridgen import Gridgen, get_ia_from_iac
 
 
 @requires_exe("gridgen")
@@ -873,3 +873,290 @@ def test_flopy_issue_1492(function_tmpdir):
         pmv.contour_array(head, levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0)
         pmv.plot_vector(spdis["qx"], spdis["qy"], color="white")
         plt.show()
+
+
+def build_gnc_gridgen(ws, layers=None, nlay=3):
+    """Build a gridgen grid with a refined block in the middle"""
+    from shapely.geometry import Polygon
+
+    nrow = ncol = 10
+    top = 1.0
+    dz = top / nlay
+    botm = [top - k * dz for k in range(1, nlay + 1)]
+
+    sim = flopy.mf6.MFSimulation(sim_name="base", sim_ws=ws)
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="base")
+    flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=1.0,
+        delc=1.0,
+        top=top,
+        botm=botm,
+    )
+
+    g = Gridgen(gwf.modelgrid, model_ws=ws)
+    polys = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
+    g.add_refinement_features(
+        polys, "polygon", 3, range(nlay) if layers is None else layers
+    )
+    g.build()
+    return g
+
+
+@pytest.mark.parametrize("nrec", [0, 1, 3])
+def test_read_qtg_gnc_dat(function_tmpdir, nrec):
+    lines = [
+        "89\t125\t88\t88\t0.125\t0.125",
+        "129\t128\t174\t175\t0.166667\t0.166667",
+        "163\t164\t124\t124\t0.125\t0.125",
+    ][:nrec]
+    (function_tmpdir / "qtg.gnc.dat").write_text("\n".join(lines))
+
+    gnc = Gridgen.read_qtg_gnc_dat(function_tmpdir)
+
+    assert gnc.dtype.names == ("n", "m", "j0", "j1", "alpha0", "alpha1")
+    assert len(gnc) == nrec
+
+    if nrec > 0:
+        # node numbers are converted to zero-based, alphas are not modified
+        assert gnc["n"][0] == 88
+        assert gnc["m"][0] == 124
+        assert gnc["j0"][0] == gnc["j1"][0] == 87
+        assert gnc["alpha0"][0] == gnc["alpha1"][0] == 0.125
+    if nrec > 1:
+        assert gnc["j0"][1] == 173
+        assert gnc["j1"][1] == 174
+        assert np.allclose(gnc["alpha1"][1], 0.166667)
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_gnc_data(function_tmpdir):
+    g = build_gnc_gridgen(function_tmpdir)
+    gnc = g.get_gnc()
+
+    # one record per line of the file gridgen wrote
+    nlines = len(
+        [
+            line
+            for line in (function_tmpdir / "qtg.gnc.dat").read_text().splitlines()
+            if line.strip()
+        ]
+    )
+    assert len(gnc) == nlines > 0
+
+    nodes = g.get_nodes()
+    for name in ("n", "m", "j0", "j1"):
+        assert gnc[name].min() >= 0
+        assert gnc[name].max() < nodes
+
+    # the ghost node is always in the coarser of the two cells
+    area = g.get_area()
+    assert np.all(area[gnc["n"]] > area[gnc["m"]])
+
+    # contributing factors must sum to less than one
+    assert np.all(gnc["alpha0"] + gnc["alpha1"] < 1.0)
+
+    # n must be connected to m, and each j must be connected to n
+    iac = g.get_iac()
+    ia = get_ia_from_iac(iac)
+    ja = g.get_ja(iac.sum())
+    for rec in gnc:
+        neighbors = ja[ia[rec["n"]] : ia[rec["n"] + 1]]
+        assert rec["m"] in neighbors
+        assert rec["j0"] in neighbors
+        assert rec["j1"] in neighbors
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_gridprops_gnc6_disv(function_tmpdir):
+    g = build_gnc_gridgen(function_tmpdir)
+    gnc = g.get_gnc()
+    gridprops = g.get_gridprops_gnc6(dis_type="disv")
+
+    assert gridprops["numalphaj"] == 2
+    assert gridprops["numgnc"] == len(gnc) == len(gridprops["gncdata"])
+
+    ncpl = g.get_gridprops_disv()["ncpl"]
+    nlay = g.get_nlay()
+    for rec, (cellidn, cellidm, j0, j1, alpha0, alpha1) in zip(
+        gnc, gridprops["gncdata"]
+    ):
+        for node, cellid in zip(
+            (rec["n"], rec["m"], rec["j0"], rec["j1"]), (cellidn, cellidm, j0, j1)
+        ):
+            assert cellid == (node // ncpl, node % ncpl)
+            assert 0 <= cellid[0] < nlay
+            assert 0 <= cellid[1] < ncpl
+        # gridgen only computes horizontal corrections
+        assert cellidn[0] == cellidm[0] == j0[0] == j1[0]
+        assert (alpha0, alpha1) == (rec["alpha0"], rec["alpha1"])
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_gridprops_gnc6_disu(function_tmpdir):
+    # refining a single layer gives a different number of nodes per layer
+    g = build_gnc_gridgen(function_tmpdir, layers=[0])
+    gnc = g.get_gnc()
+    gridprops = g.get_gridprops_gnc6(dis_type="disu")
+
+    assert gridprops["numalphaj"] == 2
+    assert gridprops["numgnc"] == len(gnc)
+    for rec, (cellidn, cellidm, j0, j1, _, _) in zip(gnc, gridprops["gncdata"]):
+        assert (cellidn, cellidm, j0, j1) == (
+            (rec["n"],),
+            (rec["m"],),
+            (rec["j0"],),
+            (rec["j1"],),
+        )
+
+    # disv cellids cannot be built when nodes per layer are not constant
+    nodelay = g.get_nodelay()
+    assert nodelay.min() != nodelay.max()
+    with pytest.raises(ValueError, match="not the same for all layers"):
+        g.get_gridprops_gnc6(dis_type="disv")
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_gridprops_gnc6_invalid(function_tmpdir):
+    g = build_gnc_gridgen(function_tmpdir, nlay=1)
+
+    with pytest.raises(ValueError, match="Unknown dis_type"):
+        g.get_gridprops_gnc6(dis_type="dis")
+
+    # n and m must be connected
+    (function_tmpdir / "qtg.gnc.dat").write_text("1\t400\t2\t2\t0.125\t0.125\n")
+    with pytest.raises(ValueError, match="is not connected to cell"):
+        g.get_gridprops_gnc6(dis_type="disv")
+    assert g.get_gridprops_gnc6(dis_type="disv", check=False)["numgnc"] == 1
+
+    # contributing factors must sum to less than one
+    (function_tmpdir / "qtg.gnc.dat").write_text("24\t34\t23\t23\t0.6\t0.6\n")
+    with pytest.raises(ValueError, match="must be less than one"):
+        g.get_gridprops_gnc6(dis_type="disv")
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_gridprops_gnc5(function_tmpdir):
+    g = build_gnc_gridgen(function_tmpdir, nlay=1)
+    gnc = g.get_gnc()
+    gridprops = g.get_gridprops_gnc5()
+
+    assert gridprops["numalphaj"] == 2
+    assert gridprops["numgnc"] == len(gnc)
+    # gridgen writes contributing factors, not conductances
+    assert gridprops["iflalphan"] == 0
+    assert gridprops["i2kn"] == 0
+    assert gridprops["isymgncn"] == 0
+
+    gncdata = gridprops["gncdata"]
+    assert gncdata.dtype == flopy.mfusg.MfUsgGnc.get_default_dtype(2, 0)
+    assert np.array_equal(gncdata["NodeN"], gnc["n"])
+    assert np.array_equal(gncdata["NodeM"], gnc["m"])
+    assert np.array_equal(gncdata["Node0"], gnc["j0"])
+    assert np.array_equal(gncdata["Node1"], gnc["j1"])
+    assert np.allclose(gncdata["Alpha0"], gnc["alpha0"])
+    assert np.allclose(gncdata["Alpha1"], gnc["alpha1"])
+
+    gridprops = g.get_gridprops_gnc5(i2kn=1, isymgncn=1)
+    assert gridprops["i2kn"] == 1
+    assert gridprops["isymgncn"] == 1
+
+
+@pytest.mark.slow
+@requires_exe("mf6", "gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_mf6disv_gnc(function_tmpdir):
+    g = build_gnc_gridgen(function_tmpdir)
+    disv_gridprops = g.get_gridprops_disv()
+    gnc_gridprops = g.get_gridprops_gnc6(dis_type="disv")
+    assert gnc_gridprops["numgnc"] > 0
+
+    chdspd = []
+    for x, y, head in [(0, 10, 1.0), (10, 0, 0.0)]:
+        ra = g.intersect([(x, y)], "point", 0)
+        chdspd.append([(0, ra["nodenumber"][0]), head])
+
+    def run(tag, gnc=False, xt3d=False):
+        ws = function_tmpdir / tag
+        sim = flopy.mf6.MFSimulation(sim_name="m", sim_ws=ws, exe_name="mf6")
+        flopy.mf6.ModflowTdis(sim)
+        flopy.mf6.ModflowIms(
+            sim,
+            linear_acceleration="bicgstab",
+            inner_dvclose=1e-9,
+            outer_dvclose=1e-9,
+        )
+        gwf = flopy.mf6.ModflowGwf(sim, modelname="m", save_flows=True)
+        flopy.mf6.ModflowGwfdisv(gwf, **disv_gridprops)
+        flopy.mf6.ModflowGwfic(gwf)
+        flopy.mf6.ModflowGwfnpf(gwf, xt3doptions=xt3d)
+        flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdspd)
+        flopy.mf6.ModflowGwfoc(
+            gwf, head_filerecord="m.hds", saverecord=[("HEAD", "ALL")]
+        )
+        if gnc:
+            flopy.mf6.ModflowGwfgnc(gwf, **gnc_gridprops)
+        sim.write_simulation()
+        success, buff = sim.run_simulation(silent=True)
+        assert success, "\n".join(buff[-25:])
+        return gwf.output.head().get_data().flatten()
+
+    head_none = run("none")
+    head_gnc = run("gnc", gnc=True)
+    head_xt3d = run("xt3d", xt3d=True)
+
+    # the correction must move the solution toward the xt3d solution
+    err_none = np.abs(head_none - head_xt3d).max()
+    err_gnc = np.abs(head_gnc - head_xt3d).max()
+    assert err_gnc < err_none / 5.0, f"gnc {err_gnc} vs uncorrected {err_none}"
+
+
+@pytest.mark.slow
+@requires_exe("mfusg", "gridgen")
+@requires_pkg("shapely", "geopandas")
+def test_mfusg_gnc(function_tmpdir):
+    g = build_gnc_gridgen(function_tmpdir, nlay=1)
+    disu_gridprops = g.get_gridprops_disu5()
+    gnc_gridprops = g.get_gridprops_gnc5()
+    assert gnc_gridprops["numgnc"] > 0
+
+    chdspd = []
+    for x, y, head in [(0, 10, 1.0), (10, 0, 0.0)]:
+        ra = g.intersect([(x, y)], "point", 0)
+        chdspd.append([ra["nodenumber"][0], head, head])
+
+    def run(tag, gnc=False):
+        ws = function_tmpdir / tag
+        m = flopy.mfusg.MfUsg(
+            modelname="m", model_ws=ws, exe_name="mfusg", structured=False
+        )
+        flopy.mfusg.MfUsgDisU(m, **disu_gridprops)
+        flopy.mfusg.MfUsgBas(m)
+        flopy.mfusg.MfUsgLpf(m)
+        flopy.modflow.ModflowChd(m, stress_period_data=chdspd)
+        flopy.mfusg.MfUsgSms(m, options="COMPLEX")
+        flopy.modflow.ModflowOc(m, stress_period_data={(0, 0): ["save head"]})
+        if gnc:
+            flopy.mfusg.MfUsgGnc(m, **gnc_gridprops)
+        m.write_input()
+        success, buff = m.run_model(silent=True)
+        assert success, "\n".join(buff[-25:])
+        return np.concatenate(flopy.utils.HeadUFile(ws / "m.hds").get_data())
+
+    head_none = run("none")
+    head_gnc = run("gnc", gnc=True)
+    assert np.abs(head_none - head_gnc).max() > 0.0
+
+    # the written package must round trip gridgen's one-based node numbers
+    written = np.genfromtxt(function_tmpdir / "gnc" / "m.gnc", skip_header=2)
+    expected = np.genfromtxt(function_tmpdir / "qtg.gnc.dat")
+    assert np.array_equal(written[:, :4], expected[:, :4])
+    assert np.allclose(written[:, 4:], expected[:, 4:], atol=1e-6)
diff --git a/flopy/mfusg/mfusg.py b/flopy/mfusg/mfusg.py
index 91c136fd1..fdcf7012b 100644
--- a/flopy/mfusg/mfusg.py
+++ b/flopy/mfusg/mfusg.py
@@ -576,4 +576,7 @@ def fmt_string(array, free=False):
             raise TypeError(msg)
         else:
             raise TypeError(f"mfusg.fmt_string error: unknown vtype in field: {field}")
-    return "".join(fmts)
+    # a free format list is read with URWORD, so the fields are separated the
+    # way MfList.fmt_string separates them; a fixed format list is read by
+    # position and relies on the field widths
+    return (" " if free else "").join(fmts)
diff --git a/flopy/mfusg/mfusggnc.py b/flopy/mfusg/mfusggnc.py
index 2114e46db..800aaf454 100644
--- a/flopy/mfusg/mfusggnc.py
+++ b/flopy/mfusg/mfusggnc.py
@@ -174,10 +174,14 @@ def write_file(self, f=None, check=False):
 
         f_gnc.write(f"{self.heading}\n")
 
+        # options are keywords, so write them as words rather than as a list
+        options = (
+            self.options if isinstance(self.options, str) else " ".join(self.options)
+        )
         f_gnc.write(
             f" {0:9d} {0:9d} {self.numgnc:9d} {self.numalphaj:9d}"
             f" {self.i2kn:9d} {self.isymgncn:9d} {self.iflalphan:9d}"
-            f" {self.options}\n"
+            f" {options}\n"
         )
 
         gdata = self.gncdata.copy()
@@ -187,7 +191,9 @@ def write_file(self, f=None, check=False):
         for idx in range(self.numalphaj):
             gdata[f"Node{idx:d}"] += 1
 
-        np.savetxt(f_gnc, gdata, fmt=fmt_string(gdata), delimiter="")
+        # the gnc list is read with URWORD, so it is not fixed-width and the
+        # contributing factors do not need to be truncated to %10.2e
+        np.savetxt(f_gnc, gdata, fmt=fmt_string(gdata, free=True), delimiter="")
 
         f_gnc.write("\n")
         f_gnc.close()
diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py
index b4129b4a7..9cd1d97ff 100644
--- a/flopy/utils/__init__.py
+++ b/flopy/utils/__init__.py
@@ -29,7 +29,7 @@
 from .formattedfile import FormattedHeadFile
 
 get_modflow = get_modflow_module.run_main
-from .gnc import get_gnc, get_gridprops_gnc6
+from .gnc import get_gnc, get_gridprops_gnc5, get_gridprops_gnc6
 from .gridintersect import GridIntersect
 from .hfb_util import make_hfb_array
 from .mflistfile import (
diff --git a/flopy/utils/gnc.py b/flopy/utils/gnc.py
index 9d86b174d..099f59ef9 100644
--- a/flopy/utils/gnc.py
+++ b/flopy/utils/gnc.py
@@ -1,9 +1,12 @@
 """
 Ghost node correction (GNC) data for quadtree-like grids.
 
-Ghost node data is computed from a grid, a grid conforming array of refinement
-levels, and the grid connectivity by :func:`get_gnc`, and is converted to
-MODFLOW 6 package input by :func:`get_gridprops_gnc6`.
+Ghost node data can be computed from a grid, a grid conforming array of
+refinement levels, and the grid connectivity, or read from gridgen output
+using :meth:`flopy.utils.gridgen.Gridgen.get_gnc`.  The record arrays produced
+by either route are converted to MODFLOW 6 and MODFLOW-USG package input by
+:func:`get_gridprops_gnc6` and :func:`get_gridprops_gnc5`.
+
 """
 
 import numpy as np
@@ -412,3 +415,69 @@ def cellid(node):
         "numalphaj": numalphaj,
         "gncdata": gncdata,
     }
+
+
+def get_gridprops_gnc5(gnc, i2kn=0, isymgncn=0, ia=None, ja=None, iac=None, check=True):
+    """
+    Get a dictionary of information needed to create a MODFLOW-USG GNC
+    Package.  The returned dictionary can be unpacked directly into the
+    MfUsgGnc constructor.
+
+    Parameters
+    ----------
+    gnc : np.recarray
+        Ghost node data with zero-based node numbers
+    i2kn : int
+        Apply the second-order correction to unconfined transmissivity
+        (default is 0).
+    isymgncn : int
+        Update the right-hand side vector for symmetric systems instead of
+        the left-hand side matrix (default is 0).
+    ia : array_like
+        Zero-based CRS row pointer, used to check connectivity
+    ja : array_like
+        Zero-based CRS column indices, used to check connectivity
+    iac : array_like
+        Number of connections per cell, used if ia is None
+    check : bool
+        Verify that each n-m pair is connected and that the contributing
+        factors sum to less than one (default is True).
+
+    Returns
+    -------
+    gridprops : dict
+
+    Notes
+    -----
+    Contributing factors are always written, so iflalphan is 0.  The default
+    asymmetric implementation requires an asymmetric solver.  numgnc is zero
+    for a grid without ghost nodes, in which case the package should not be
+    created.
+
+    """
+    # imported here because flopy.mfusg imports flopy.utils
+    from ..mfusg.mfusggnc import MfUsgGnc
+
+    if check:
+        _check_gnc(gnc, ia=ia, ja=ja, iac=iac)
+
+    numalphaj = get_numalphaj(gnc)
+    iflalphan = 0
+    gncdata = MfUsgGnc.get_empty(
+        numgnc=len(gnc), numalphaj=numalphaj, iflalphan=iflalphan
+    )
+    # MfUsgGnc.write_file converts to one-based node numbers
+    gncdata["NodeN"] = gnc["n"]
+    gncdata["NodeM"] = gnc["m"]
+    for i in range(numalphaj):
+        gncdata[f"Node{i}"] = gnc[f"j{i}"]
+        gncdata[f"Alpha{i}"] = gnc[f"alpha{i}"]
+
+    return {
+        "numgnc": len(gncdata),
+        "numalphaj": numalphaj,
+        "i2kn": i2kn,
+        "isymgncn": isymgncn,
+        "iflalphan": iflalphan,
+        "gncdata": gncdata,
+    }
diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py
index 1857e8457..2cd4205bb 100644
--- a/flopy/utils/gridgen.py
+++ b/flopy/utils/gridgen.py
@@ -15,6 +15,7 @@
 from ..modflow import ModflowDis
 from ..utils import import_optional_dependency
 from ..utils.flopy_io import relpath_safe
+from .gnc import get_gnc_dtype, get_gridprops_gnc5, get_gridprops_gnc6
 from .util_array import Util2d
 
 # todo
@@ -23,6 +24,12 @@
 # support an asciigrid option for top and bottom interpolation
 # add intersection capability
 
+# gridgen always writes two contributing cells per ghost node.  When only one
+# contributing cell exists it is repeated with alpha halved, which both
+# MODFLOW 6 and MODFLOW-USG accumulate into the same matrix position.
+GNC_NUMALPHAJ = 2
+GNC_DTYPE = get_gnc_dtype(GNC_NUMALPHAJ)
+
 
 def read1d(f, a):
     """
@@ -1256,6 +1263,31 @@ def get_anglex(self, fldr=None):
         anglex = np.where(fldr == 2, 4.712389, anglex)
         return anglex
 
+    def get_gnc(self):
+        """
+        Get the ghost node correction data computed by gridgen
+
+        Returns
+        -------
+        gnc : np.recarray
+            Record array with fields n, m, j0, j1, alpha0, and alpha1.  Node
+            numbers are zero-based.  Cell n is the cell containing the ghost
+            node, cell m is the connecting cell, and cells j0 and j1 are the
+            contributing cells.
+
+        Notes
+        -----
+        Gridgen only computes horizontal ghost node corrections and drops
+        records where a contributing cell is inactive.
+
+        """
+        return self.read_qtg_gnc_dat(model_ws=self.model_ws)
+
+    def _gnc_connectivity(self):
+        """Return the zero-based ia and ja arrays"""
+        iac = self.get_iac()
+        return get_ia_from_iac(iac), self.get_ja(iac.sum())
+
     def get_verts_iverts(self, ncells, verbose=False):
         """
         Return a 2d array of x and y vertices and a list of size ncells that
@@ -1511,6 +1543,96 @@ def get_gridprops_disv(self):
 
         return gridprops
 
+    def _gnc_ncpl(self, dis_type):
+        """Return the number of cells per layer for dis_type"""
+        if dis_type.lower() != "disv":
+            return None
+        nodelay = self.get_nodelay()
+        ncpl = nodelay.min()
+        if ncpl != nodelay.max():
+            raise ValueError(
+                "Cannot create DISV ghost node properties because the "
+                "number of cells is not the same for all layers"
+            )
+        return ncpl
+
+    def get_gridprops_gnc6(self, dis_type="disv", check=True):
+        """
+        Get a dictionary of information needed to create a MODFLOW 6 GNC
+        Package.  The returned dictionary can be unpacked directly into the
+        ModflowGwfgnc constructor.
+
+        Parameters
+        ----------
+        dis_type : str
+            Discretization the cellids are built for.  Valid options are
+            'disv' (default) and 'disu'.
+        check : bool
+            Verify that each n-m pair is connected and that the contributing
+            factors sum to less than one (default is True).
+
+        Returns
+        -------
+        gridprops : dict
+
+        Notes
+        -----
+        The correction is applied implicitly unless the explicit option is
+        set, so the BICGSTAB linear acceleration option should be specified
+        in the IMS Package.  numgnc is zero for a grid without ghost nodes,
+        in which case the package should not be created.
+
+        """
+        ia, ja = self._gnc_connectivity() if check else (None, None)
+        return get_gridprops_gnc6(
+            self.get_gnc(),
+            dis_type=dis_type,
+            ncpl=self._gnc_ncpl(dis_type),
+            ia=ia,
+            ja=ja,
+            check=check,
+        )
+
+    def get_gridprops_gnc5(self, i2kn=0, isymgncn=0, check=True):
+        """
+        Get a dictionary of information needed to create a MODFLOW-USG GNC
+        Package.  The returned dictionary can be unpacked directly into the
+        MfUsgGnc constructor.
+
+        Parameters
+        ----------
+        i2kn : int
+            Apply the second-order correction to unconfined transmissivity
+            (default is 0).
+        isymgncn : int
+            Update the right-hand side vector for symmetric systems instead
+            of the left-hand side matrix (default is 0).
+        check : bool
+            Verify that each n-m pair is connected and that the contributing
+            factors sum to less than one (default is True).
+
+        Returns
+        -------
+        gridprops : dict
+
+        Notes
+        -----
+        Gridgen writes contributing factors, so iflalphan is always 0.  The
+        default asymmetric implementation requires an asymmetric solver.
+        numgnc is zero for a grid without ghost nodes, in which case the
+        package should not be created.
+
+        """
+        ia, ja = self._gnc_connectivity() if check else (None, None)
+        return get_gridprops_gnc5(
+            self.get_gnc(),
+            i2kn=i2kn,
+            isymgncn=isymgncn,
+            ia=ia,
+            ja=ja,
+            check=check,
+        )
+
     def get_gridprops_vertexgrid(self):
         """
         Get a dictionary of information needed to create a flopy VertexGrid.
@@ -2103,3 +2225,27 @@ def read_qtg_fahl_dat(model_ws: Union[str, PathLike], nja: int):
         fname = os.path.join(model_ws, "qtg.fahl.dat")
         with open(fname, "r") as f:
             return read1d(f=f, a=np.empty((nja), dtype=np.float32))
+
+    @staticmethod
+    def read_qtg_gnc_dat(model_ws: Union[str, PathLike]):
+        """Read qtg.gnc.dat file
+
+        Parameters
+        ----------
+        model_ws : str or PathLike
+            Directory where file is stored
+
+        Returns
+        -------
+        np.recarray
+            Ghost node records with zero-based node numbers.  The record is
+            empty if gridgen did not find any ghost nodes.
+        """
+        fname = os.path.join(model_ws, "qtg.gnc.dat")
+        # gridgen writes an empty file when the grid has no ghost nodes
+        if os.path.getsize(fname) == 0:
+            return np.recarray((0,), dtype=GNC_DTYPE)
+        gnc = np.atleast_1d(np.genfromtxt(fname, dtype=GNC_DTYPE))
+        for name in ("n", "m", "j0", "j1"):
+            gnc[name] -= 1
+        return gnc.view(np.recarray)