diff --git a/autotest/test_gnc.py b/autotest/test_gnc.py
new file mode 100644
index 000000000..9b3b05157
--- /dev/null
+++ b/autotest/test_gnc.py
@@ -0,0 +1,339 @@
+"""
+Tests for ghost node correction (GNC) data computed from a grid.
+
+Cases:
+  - synthetic : a hand built grid with known contributing cells and factors.
+  - gridgen   : the computed data must reproduce gridgen's qtg.gnc.dat.
+"""
+
+import numpy as np
+import pytest
+from modflow_devtools.markers import requires_exe, requires_pkg
+
+import flopy
+from flopy.discretization import UnstructuredGrid, VertexGrid
+from flopy.utils.gnc import (
+    _check_gnc,
+    get_gnc,
+    get_gnc_dtype,
+    get_gridprops_gnc6,
+    get_numalphaj,
+)
+from flopy.utils.gridgen import Gridgen
+
+
+def synthetic_grid(connectivity=True):
+    """A coarse cell with four fine cells north and two fine cells east
+
+    Cell 0 is a size 4 cell centered on the origin.  Cells 1 and 2 are size 2
+    cells on its east face and cells 3 to 6 are size 1 cells on its north
+    face, so the north face has a two level refinement jump.
+    """
+    centers = [
+        (0.0, 0.0, 4.0),
+        (3.0, 1.0, 2.0),
+        (3.0, -1.0, 2.0),
+        (-1.5, 3.0, 1.0),
+        (-0.5, 3.0, 1.0),
+        (0.5, 3.0, 1.0),
+        (1.5, 3.0, 1.0),
+    ]
+    vertices, iverts = [], []
+    for x, y, size in centers:
+        half = size / 2.0
+        iv = []
+        for vx, vy in [
+            (x - half, y - half),
+            (x + half, y - half),
+            (x + half, y + half),
+            (x - half, y + half),
+        ]:
+            iv.append(len(vertices))
+            vertices.append([len(vertices), vx, vy])
+        iverts.append(iv)
+
+    # every cell is connected to cell 0 only
+    conn = {}
+    if connectivity:
+        conn["iac"] = np.array([7, 2, 2, 2, 2, 2, 2])
+        conn["ja"] = np.array([0, 1, 2, 3, 4, 5, 6, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0])
+    grid = UnstructuredGrid(
+        vertices=vertices,
+        iverts=iverts,
+        xcenters=np.array([c[0] for c in centers]),
+        ycenters=np.array([c[1] for c in centers]),
+        ncpl=np.array([len(centers)]),
+        **conn,
+    )
+    level = np.array([0, 1, 1, 2, 2, 2, 2])
+    return grid, level
+
+
+def gnc_key(gnc):
+    """Sorted view of ghost node records for order independent comparison"""
+    numalphaj = get_numalphaj(gnc)
+    nodes = np.sort(np.column_stack([gnc[f"j{i}"] for i in range(numalphaj)]), axis=1)
+    alpha = np.sort(
+        np.column_stack([gnc[f"alpha{i}"] for i in range(numalphaj)]), axis=1
+    )
+    key = np.column_stack([gnc["n"], gnc["m"], nodes, alpha])
+    return key[np.lexsort(key.T[::-1])]
+
+
+def build_gridgen(ws, nlay=3, level=3, layers=None, **kwargs):
+    """Build a gridgen grid with a refined block in the middle"""
+    from shapely.geometry import Polygon
+
+    botm = [1.0 - k * (1.0 / nlay) 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=10, ncol=10, delr=1.0, delc=1.0, top=1.0, botm=botm
+    )
+    g = Gridgen(gwf.modelgrid, model_ws=ws, **kwargs)
+    g.add_refinement_features(
+        [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])],
+        "polygon",
+        level,
+        range(nlay) if layers is None else layers,
+    )
+    g.build()
+    return g
+
+
+def test_get_gnc_dtype():
+    dtype = get_gnc_dtype(3)
+    assert dtype.names == (
+        "n",
+        "m",
+        "j0",
+        "j1",
+        "j2",
+        "alpha0",
+        "alpha1",
+        "alpha2",
+    )
+    assert get_numalphaj(np.recarray((0,), dtype=dtype)) == 3
+
+
+@pytest.mark.parametrize("use_level", [True, False])
+def test_get_gnc_synthetic(use_level):
+    grid, level = synthetic_grid()
+    gnc = get_gnc(grid, level=level if use_level else None)
+
+    # cell 0 has four contributing cells on its north face
+    assert get_numalphaj(gnc) == 4
+    assert len(gnc) == 3
+
+    # connections without a contributing cell on the offset side are dropped
+    assert sorted(zip(gnc["n"], gnc["m"])) == [(0, 1), (0, 5), (0, 6)]
+
+    records = {(rec["n"], rec["m"]): rec for rec in gnc}
+
+    # cell 1 is offset 1.0 from cell 0 and the north cells are 3.0 away, so
+    # the factors total 1/3 and are shared by four contributing cells
+    rec = records[(0, 1)]
+    assert [rec[f"j{i}"] for i in range(4)] == [3, 4, 5, 6]
+    assert np.allclose([rec[f"alpha{i}"] for i in range(4)], 1.0 / 12.0)
+
+    # cell 6 is offset 1.5 and the east cells are 3.0 away, so the factors
+    # total 0.5 shared by two contributing cells, the first of which is
+    # repeated three times to fill numalphaj
+    rec = records[(0, 6)]
+    assert [rec[f"j{i}"] for i in range(4)] == [1, 1, 1, 2]
+    assert np.allclose(
+        [rec[f"alpha{i}"] for i in range(4)], [1 / 12, 1 / 12, 1 / 12, 0.25]
+    )
+
+
+def factor_totals(gnc):
+    """Total contributing factor of each cell, keyed by ghost node"""
+    numalphaj = get_numalphaj(gnc)
+    totals = {}
+    for rec in gnc:
+        entry = totals.setdefault((rec["n"], rec["m"]), {})
+        for i in range(numalphaj):
+            node = rec[f"j{i}"]
+            entry[node] = entry.get(node, 0.0) + rec[f"alpha{i}"]
+    return totals
+
+
+def test_get_gnc_padding_preserves_factors():
+    grid, level = synthetic_grid()
+    unpadded = get_gnc(grid, level=level)
+    padded = get_gnc(grid, level=level, numalphaj=6)
+    assert get_numalphaj(unpadded) == 4
+    assert get_numalphaj(padded) == 6
+
+    # padding repeats a contributing cell and splits its factor, so the total
+    # factor of every cell must be unchanged
+    expected, actual = factor_totals(unpadded), factor_totals(padded)
+    assert expected.keys() == actual.keys()
+    for ghost, cells in expected.items():
+        assert cells.keys() == actual[ghost].keys()
+        for node, alpha in cells.items():
+            assert np.allclose(alpha, actual[ghost][node])
+        assert sum(cells.values()) < 1.0
+
+
+def test_get_gnc_numalphaj_too_small():
+    grid, level = synthetic_grid()
+    with pytest.raises(ValueError, match="more than numalphaj"):
+        get_gnc(grid, level=level, numalphaj=2)
+
+
+def test_get_gnc_supplied_connectivity():
+    grid, level = synthetic_grid()
+    bare, _ = synthetic_grid(connectivity=False)
+    assert bare.iac is None
+
+    # the synthetic cells do not share vertices, so connectivity cannot be
+    # built from shared edges and has to be supplied
+    assert len(get_gnc(bare, level=level)) == 0
+
+    gnc = get_gnc(bare, level=level, iac=grid.iac, ja=grid.ja)
+    assert np.array_equal(gnc_key(gnc), gnc_key(get_gnc(grid, level=level)))
+
+
+def test_check_gnc():
+    dtype = get_gnc_dtype(2)
+    gnc = np.recarray((1,), dtype=dtype)
+    gnc[0] = (0, 1, 2, 2, 0.6, 0.6)
+    with pytest.raises(ValueError, match="must be less than one"):
+        _check_gnc(gnc)
+
+    gnc[0] = (0, 1, 2, 2, 0.1, 0.1)
+    _check_gnc(gnc)
+
+    # cell 0 is connected to cell 2 but not to cell 1
+    ia = np.array([0, 2, 3, 5])
+    ja = np.array([0, 2, 1, 2, 0])
+    with pytest.raises(ValueError, match="is not connected to cell"):
+        _check_gnc(gnc, ia=ia, ja=ja)
+
+
+@pytest.mark.parametrize("n,m", [(5, 1), (1, 5), (-3, 1), (1, -3)])
+def test_check_gnc_node_out_of_range(n, m):
+    """A node outside the grid is reported rather than indexed"""
+    gnc = np.recarray((1,), dtype=get_gnc_dtype(1))
+    gnc[0] = (n, m, 0, 0.1)
+    # a negative node would otherwise wrap and check the wrong cell
+    with pytest.raises(ValueError, match="which is not a cell of a grid"):
+        _check_gnc(gnc, ia=np.array([0, 2, 3]), ja=np.array([0, 1, 1, 0]))
+
+
+def test_get_gridprops_gnc6_requires_ncpl():
+    grid, level = synthetic_grid()
+    gnc = get_gnc(grid, level=level)
+    with pytest.raises(ValueError, match="ncpl is required"):
+        get_gridprops_gnc6(gnc, dis_type="disv")
+    with pytest.raises(ValueError, match="Unknown dis_type"):
+        get_gridprops_gnc6(gnc, dis_type="dis")
+
+    gridprops = get_gridprops_gnc6(gnc, dis_type="disu")
+    assert gridprops["numalphaj"] == 4
+    assert gridprops["numgnc"] == len(gnc)
+    assert gridprops["gncdata"][0][0] == (gnc["n"][0],)
+
+
+@requires_exe("gridgen")
+@requires_pkg("shapely", "geopandas")
+@pytest.mark.parametrize(
+    "nlay,level,layers,smoothing",
+    [
+        (3, 3, None, 1),
+        (1, 2, None, 1),
+        (1, 4, None, 1),
+        (3, 3, [0], 1),
+        (1, 3, None, 2),
+    ],
+)
+def test_get_gnc_matches_gridgen(function_tmpdir, nlay, level, layers, smoothing):
+    g = build_gridgen(
+        function_tmpdir,
+        nlay=nlay,
+        level=level,
+        layers=layers,
+        smoothing_level_horizontal=smoothing,
+        smoothing_level_vertical=smoothing,
+    )
+    # gridgen writes the ghost node data whenever it exports a grid
+    expected = np.atleast_1d(
+        np.genfromtxt(function_tmpdir / "qtg.gnc.dat", dtype=get_gnc_dtype(2))
+    )
+    for name in ("n", "m", "j0", "j1"):
+        expected[name] -= 1
+    assert len(expected) > 0
+
+    grid = UnstructuredGrid(**g.get_gridprops_unstructuredgrid())
+    gnc = get_gnc(grid, numalphaj=2)
+
+    assert len(gnc) == len(expected)
+    # gridgen writes the factors with six significant digits
+    assert np.allclose(gnc_key(gnc), gnc_key(expected), atol=2.0e-6)
+
+
+def square_cells(centers):
+    """Build vertices and iverts for a list of (x, y, size) squares"""
+    vertices, iverts = [], []
+    for x, y, size in centers:
+        half = size / 2.0
+        iv = []
+        for vx, vy in [
+            (x - half, y - half),
+            (x + half, y - half),
+            (x + half, y + half),
+            (x - half, y + half),
+        ]:
+            iv.append(len(vertices))
+            vertices.append([len(vertices), vx, vy])
+        iverts.append(iv)
+    return vertices, iverts
+
+
+def test_get_gnc_ihc():
+    """A connection marked vertical is not corrected and does not contribute"""
+    grid, level = synthetic_grid()
+
+    # cells 0 and 1 are connected at positions 1 and 8 of ja
+    ihc = np.ones(len(grid.ja), dtype=int)
+    ihc[[1, 8]] = 0
+    gnc = get_gnc(grid, level=level, ihc=ihc)
+
+    # the ghost node on the 0 to 1 connection is gone, and cell 1 is no longer
+    # available as a contributing cell, which leaves cell 2 on its own
+    assert sorted((rec["n"], rec["m"]) for rec in gnc) == [(0, 5), (0, 6)]
+    assert get_numalphaj(gnc) == 1
+    assert np.allclose(sorted(gnc["alpha0"]), [1.0 / 6.0, 0.5])
+
+
+def test_get_gnc_aligned_finer_cell():
+    """A finer cell centered on the face of a coarse cell needs no ghost node"""
+    centers = [(0.0, 0.0, 4.0), (3.0, 0.0, 2.0)]
+    vertices, iverts = square_cells(centers)
+    grid = UnstructuredGrid(
+        vertices=vertices,
+        iverts=iverts,
+        xcenters=np.array([c[0] for c in centers]),
+        ycenters=np.array([c[1] for c in centers]),
+        ncpl=np.array([len(centers)]),
+        iac=np.array([2, 2]),
+        ja=np.array([0, 1, 1, 0]),
+    )
+    assert len(get_gnc(grid, level=np.array([0, 1]))) == 0
+
+
+def test_get_gnc_connectivity_needs_constant_ncpl():
+    """Shared edges cannot give the layer layout when ncpl varies by layer"""
+    centers = [(0.5, 0.5, 1.0), (1.5, 0.5, 1.0), (2.5, 0.5, 1.0), (0.5, 1.5, 1.0)]
+    vertices, iverts = square_cells(centers)
+    grid = UnstructuredGrid(
+        vertices=vertices,
+        iverts=iverts,
+        xcenters=np.array([c[0] for c in centers]),
+        ycenters=np.array([c[1] for c in centers]),
+        ncpl=np.array([3, 1]),
+    )
+    assert grid.iac is None
+    with pytest.raises(ValueError, match="different"):
+        get_gnc(grid)
diff --git a/flopy/mfusg/__init__.py b/flopy/mfusg/__init__.py
index 579db55f9..b5d4b621a 100644
--- a/flopy/mfusg/__init__.py
+++ b/flopy/mfusg/__init__.py
@@ -36,7 +36,6 @@
     "MfUsgDpt",
     "MfUsgEvt",
     "MfUsgGnc",
-    "MfUsgGnc",
     "MfUsgLak",
     "MfUsgLpf",
     "MfUsgMdt",
diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py
index 0da0e3d10..b4129b4a7 100644
--- a/flopy/utils/__init__.py
+++ b/flopy/utils/__init__.py
@@ -29,6 +29,7 @@
 from .formattedfile import FormattedHeadFile
 
 get_modflow = get_modflow_module.run_main
+from .gnc import get_gnc, 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
new file mode 100644
index 000000000..9d86b174d
--- /dev/null
+++ b/flopy/utils/gnc.py
@@ -0,0 +1,414 @@
+"""
+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`.
+"""
+
+import numpy as np
+
+
+def get_gnc_dtype(numalphaj):
+    """
+    Get the record dtype for ghost node data with numalphaj contributing
+    cells
+
+    Parameters
+    ----------
+    numalphaj : int
+        Number of contributing cells per ghost node
+
+    Returns
+    -------
+    dtype : np.dtype
+
+    """
+    dtype = [("n", int), ("m", int)]
+    dtype += [(f"j{i}", int) for i in range(numalphaj)]
+    dtype += [(f"alpha{i}", float) for i in range(numalphaj)]
+    return np.dtype(dtype)
+
+
+def get_numalphaj(gnc):
+    """
+    Get the number of contributing cells in a ghost node record array
+
+    Parameters
+    ----------
+    gnc : np.recarray
+        Ghost node data
+
+    Returns
+    -------
+    numalphaj : int
+
+    """
+    return len([name for name in gnc.dtype.names if name.startswith("j")])
+
+
+def _gnc_nodes(gnc):
+    """Return the cellid column names of a ghost node record array"""
+    numalphaj = get_numalphaj(gnc)
+    return ["n", "m"] + [f"j{i}" for i in range(numalphaj)]
+
+
+def _get_ia(ia=None, iac=None):
+    """Build the zero-based ia array from ia or iac"""
+    if ia is not None:
+        return np.asarray(ia, dtype=int)
+    if iac is None:
+        return None
+    return np.concatenate(([0], np.cumsum(np.asarray(iac, dtype=int))))
+
+
+def _check_gnc(gnc, ia=None, ja=None, iac=None):
+    """
+    Raise if ghost node records are not valid MODFLOW input
+
+    Parameters
+    ----------
+    gnc : np.recarray
+        Ghost node data
+    ia : array_like
+        Zero-based CRS row pointer.  Connectivity is not checked if ja is
+        None or if both ia and iac are None.
+    ja : array_like
+        Zero-based CRS column indices
+    iac : array_like
+        Number of connections per cell, used if ia is None
+
+    """
+    ia = _get_ia(ia, iac)
+    numalphaj = get_numalphaj(gnc)
+    alpha = np.zeros(len(gnc))
+    for i in range(numalphaj):
+        alpha += gnc[f"alpha{i}"]
+    for irec, total in enumerate(alpha):
+        if total >= 1.0:
+            raise ValueError(
+                f"gnc record {irec}: contributing factors sum to {total}, "
+                "which must be less than one"
+            )
+
+    if ia is None or ja is None:
+        return
+    ia = np.asarray(ia, dtype=int)
+    ja = np.asarray(ja, dtype=int)
+    nodes = ia.shape[0] - 1
+    for irec, rec in enumerate(gnc):
+        n, m = rec["n"], rec["m"]
+        # a node number outside the grid would either index past the end of ia
+        # or, when it is negative, silently wrap and check the wrong cell
+        for name, node in (("n", n), ("m", m)):
+            if not 0 <= node < nodes:
+                raise ValueError(
+                    f"gnc record {irec}: cell {name} is {node}, which is not a "
+                    f"cell of a grid with {nodes} cells"
+                )
+        # MODFLOW 6 rejects a ghost node whose n-m connection is absent
+        if m not in ja[ia[n] : ia[n + 1]]:
+            raise ValueError(
+                f"gnc record {irec}: cell {n} is not connected to cell {m}"
+            )
+
+
+def _node_centers(modelgrid):
+    """Return cell center arrays with one value per node"""
+    xc = np.asarray(modelgrid.xcellcenters).ravel()
+    yc = np.asarray(modelgrid.ycellcenters).ravel()
+    nnodes = modelgrid.nnodes
+    if xc.shape[0] == nnodes:
+        return xc, yc
+    # a vertex grid stores one value per cell2d, repeated for every layer
+    nlay = nnodes // xc.shape[0]
+    return np.tile(xc, nlay), np.tile(yc, nlay)
+
+
+def _node_layers(modelgrid):
+    """Return the layer index of every node"""
+    ncpl = modelgrid.ncpl
+    if np.isscalar(ncpl):
+        ncpl = np.full(modelgrid.nlay, ncpl, dtype=int)
+    return np.repeat(np.arange(len(ncpl)), ncpl)
+
+
+def _shared_edge_connectivity(modelgrid):
+    """Build ia and ja from the cells that share an edge in every layer"""
+    neighbors = modelgrid.neighbors(method="rook")
+    ncpl = modelgrid.ncpl
+    if not np.isscalar(ncpl):
+        if ncpl.min() != ncpl.max():
+            raise ValueError(
+                "Connectivity cannot be built for a grid with a different "
+                "number of cells in each layer, supply ia or iac and ja"
+            )
+        ncpl = int(ncpl.min())
+    nlay = modelgrid.nnodes // ncpl
+
+    # ghost nodes are horizontal, so only the connections within a layer
+    # are needed and the same layout is repeated for every layer
+    iac, ja = [], []
+    for k in range(nlay):
+        for icpl in range(ncpl):
+            conn = sorted(neighbors.get(icpl, []))
+            iac.append(len(conn) + 1)
+            ja.append(k * ncpl + icpl)
+            ja.extend(k * ncpl + j for j in conn)
+    return _get_ia(iac=iac), np.array(ja, dtype=int)
+
+
+def _contributing_cells(dn, naxis, d_nm, rtol):
+    """
+    Select the cells that contribute to the ghost node on one connection
+
+    Parameters
+    ----------
+    dn : np.ndarray
+        Offset from cell n to each of its neighbors, one row per neighbor
+    naxis : np.ndarray
+        Connection axis of each neighbor of cell n
+    d_nm : np.ndarray
+        Offset from cell n to the connected cell m
+    rtol : float
+        Relative tolerance on the transverse offset
+
+    Returns
+    -------
+    sel : np.ndarray or None
+        Mask of the contributing neighbors, None when the connection needs no
+        ghost node
+    alpha : float
+        Total contributing factor
+
+    """
+    axis = int(np.argmax(np.abs(d_nm)))
+    trans = 1 - axis
+    offset = d_nm[trans]
+    # the cell centers line up across the face, so there is nothing to correct
+    if abs(offset) <= rtol * abs(d_nm[axis]):
+        return None, 0.0
+
+    # the contributing cells are the neighbors of n on the side of the offset
+    sel = (naxis == trans) & (np.sign(dn[:, trans]) == np.sign(offset))
+    if not sel.any():
+        return None, 0.0
+    return sel, abs(offset) / np.abs(dn[sel, trans]).mean()
+
+
+def _cell_areas(modelgrid):
+    """Return the area of every node from the cell vertices"""
+    verts = np.asarray(modelgrid.verts)
+    iverts = modelgrid.iverts
+    areas = np.empty(len(iverts))
+    for i, iv in enumerate(iverts):
+        iv = [j for j in iv if j is not None]
+        x, y = verts[iv, 0], verts[iv, 1]
+        areas[i] = 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
+    if areas.shape[0] != modelgrid.nnodes:
+        nlay = modelgrid.nnodes // areas.shape[0]
+        areas = np.tile(areas, nlay)
+    return areas
+
+
+def get_gnc(
+    modelgrid,
+    level=None,
+    ia=None,
+    ja=None,
+    iac=None,
+    ihc=None,
+    numalphaj=None,
+    rtol=1.0e-6,
+):
+    """
+    Compute ghost node correction data for a quadtree-like grid
+
+    A ghost node is added in cell n for every horizontal connection to a
+    finer cell m whose center is offset from the center of n transverse to
+    the connection.  The head at the ghost node is interpolated between cell
+    n and the neighbors of n on the side of the offset.
+
+    Parameters
+    ----------
+    modelgrid : flopy.discretization.UnstructuredGrid or VertexGrid
+        Grid the ghost nodes are computed for.  Connectivity is taken from
+        the grid when ia, ja, and iac are None, either from the iac and ja
+        the grid carries or from the cells that share an edge.
+    level : array_like
+        Grid conforming array of refinement levels, where a larger value is a
+        finer cell.  Cell areas are used if None.
+    ia : array_like
+        Zero-based CRS row pointer
+    ja : array_like
+        Zero-based CRS column indices, with the diagonal first in each row
+    iac : array_like
+        Number of connections per cell, used if ia is None
+    ihc : array_like
+        Connection type for each entry in ja, where 0 is a vertical
+        connection.  Connections between cells in different layers are
+        treated as vertical if None.
+    numalphaj : int
+        Number of contributing cells written per ghost node.  The largest
+        number found is used if None.  Records with fewer contributing cells
+        repeat a cell rather than pad with zeros, which MODFLOW-USG requires.
+    rtol : float
+        Relative tolerance used to decide whether the center of cell m is
+        offset from the center of cell n (default is 1.0e-6).
+
+    Returns
+    -------
+    gnc : np.recarray
+        Record array with fields n, m, j0 to j[numalphaj-1], and alpha0 to
+        alpha[numalphaj-1].  Node numbers are zero-based.
+
+    Notes
+    -----
+    Only horizontal corrections are computed, matching gridgen.
+
+    """
+    ia = _get_ia(ia, iac)
+    if ia is None:
+        ia = _get_ia(iac=getattr(modelgrid, "iac", None))
+    if ja is None:
+        ja = getattr(modelgrid, "ja", None)
+    if ia is None or ja is None:
+        # a vertex grid does not carry connectivity, so build it from the
+        # cells that share an edge
+        ia, ja = _shared_edge_connectivity(modelgrid)
+    ja = np.asarray(ja, dtype=int)
+
+    nnodes = modelgrid.nnodes
+    xc, yc = _node_centers(modelgrid)
+
+    # a larger size is a coarser cell
+    if level is None:
+        size = _cell_areas(modelgrid)
+    else:
+        size = -np.asarray(level, dtype=float).ravel()
+        if size.shape[0] != nnodes:
+            size = np.tile(size, nnodes // size.shape[0])
+
+    if ihc is None:
+        layer = _node_layers(modelgrid)
+        horizontal = layer[ja] == layer[np.repeat(np.arange(nnodes), np.diff(ia))]
+    else:
+        horizontal = np.asarray(ihc) != 0
+
+    records = []
+    for n in range(nnodes):
+        ipos = np.arange(ia[n] + 1, ia[n + 1])
+        ipos = ipos[horizontal[ipos]]
+        conn = ja[ipos]
+        if conn.size == 0:
+            continue
+
+        # column 0 is the x offset and column 1 the y offset to each neighbor
+        d = np.column_stack((xc[conn] - xc[n], yc[conn] - yc[n]))
+        axis = np.argmax(np.abs(d), axis=1)
+
+        for k, m in enumerate(conn):
+            if size[m] >= size[n]:
+                continue
+            sel, alpha = _contributing_cells(d, axis, d[k], rtol)
+            if sel is None:
+                continue
+            js = conn[sel]
+            records.append((n, m, js, alpha / js.size))
+
+    if numalphaj is None:
+        numalphaj = max((len(rec[2]) for rec in records), default=1)
+    dtype = get_gnc_dtype(numalphaj)
+
+    gnc = np.recarray((len(records),), dtype=dtype)
+    for irec, (n, m, js, alpha) in enumerate(records):
+        if js.size > numalphaj:
+            raise ValueError(
+                f"gnc record {irec}: cell {n} has {js.size} contributing cells, "
+                f"which is more than numalphaj of {numalphaj}"
+            )
+        # pad by repeating the first contributing cell and splitting its
+        # factor, which both MODFLOW 6 and MODFLOW-USG accumulate
+        nrepeat = numalphaj - js.size + 1
+        nodes = np.concatenate((np.repeat(js[:1], nrepeat), js[1:]))
+        alphas = np.concatenate(
+            (np.full(nrepeat, alpha / nrepeat), np.full(js.size - 1, alpha))
+        )
+        gnc["n"][irec] = n
+        gnc["m"][irec] = m
+        for i in range(numalphaj):
+            gnc[f"j{i}"][irec] = nodes[i]
+            gnc[f"alpha{i}"][irec] = alphas[i]
+
+    return gnc
+
+
+def get_gridprops_gnc6(
+    gnc, dis_type="disv", ncpl=None, ia=None, ja=None, iac=None, 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
+    ----------
+    gnc : np.recarray
+        Ghost node data with zero-based node numbers
+    dis_type : str
+        Discretization the cellids are built for.  Valid options are 'disv'
+        (default) and 'disu'.
+    ncpl : int
+        Number of cells per layer, required for 'disv'
+    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
+    -----
+    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.
+
+    """
+    if check:
+        _check_gnc(gnc, ia=ia, ja=ja, iac=iac)
+
+    dis_type = dis_type.lower()
+    if dis_type == "disv":
+        if ncpl is None:
+            raise ValueError("ncpl is required to build disv cellids")
+
+        def cellid(node):
+            return (node // ncpl, node % ncpl)
+    elif dis_type == "disu":
+
+        def cellid(node):
+            return (node,)
+    else:
+        raise ValueError(f"Unknown dis_type {dis_type}, expected 'disv' or 'disu'")
+
+    numalphaj = get_numalphaj(gnc)
+    names = _gnc_nodes(gnc)
+    gncdata = [
+        tuple(cellid(rec[name]) for name in names)
+        + tuple(rec[f"alpha{i}"] for i in range(numalphaj))
+        for rec in gnc
+    ]
+
+    return {
+        "numgnc": len(gncdata),
+        "numalphaj": numalphaj,
+        "gncdata": gncdata,
+    }