diff --git a/docs/examples/quickstart.py b/docs/examples/quickstart.py
index 3e64c1b4..d001872c 100644
--- a/docs/examples/quickstart.py
+++ b/docs/examples/quickstart.py
@@ -45,7 +45,7 @@
 
 assert chd.data["head"][0, 0] == 1.0
 assert chd.data.head.sel(per=0)[99] == 0.0
-assert np.allclose(chd.data.head[:, 1:99], np.full(98, 1e30))
+assert np.allclose(chd.data.head[:, 1:99], np.full(98, 3e30))
 
 assert gwf.dis.data.botm.sel(lay=0, col=0, row=0) == 0.0
 
diff --git a/docs/examples/twri.py b/docs/examples/twri.py
index 3f743e86..35d17d24 100644
--- a/docs/examples/twri.py
+++ b/docs/examples/twri.py
@@ -8,6 +8,21 @@
 
 import flopy4
 
+
+def plot_head(head, workspace):
+    import matplotlib.pyplot as plt
+
+    # Plot head results
+    plt.figure(figsize=(10, 6))
+    head.isel(layer=0, time=0).plot.contourf()
+    plt.title("Filled Contour Plot TWRI Head")
+    plt.xlabel("x")
+    plt.ylabel("y")
+    plt.grid(True)
+    plt.savefig(workspace / "head.png", dpi=300, bbox_inches="tight")
+    plt.close()
+
+
 # Timing
 time = flopy4.mf6.utils.time.Time.from_timestamps(
     ["2000-01-01", "2000-01-02", "2000-01-03", "2000-01-04"]
@@ -28,7 +43,7 @@
 grid = flopy4.mf6.utils.grid.StructuredGrid(
     nlay=nlay, nrow=nrow, ncol=ncol, top=top, botm=bottom, delr=delr, delc=delc, idomain=idomain
 )
-dims = {"nper": nper, **dict(grid.dataset.sizes)}  # TODO: temporary
+dims = {"nper": nper, "ncpl": nrow * ncol, **dict(grid.dataset.sizes)}  # TODO: temporary
 
 # Grid discretization
 # TODO: xorigin, yorigin
@@ -84,9 +99,9 @@
 )
 
 # Uniform recharge on the top layer
-rch_rate = np.stack(np.full((nlay, nrow, ncol), flopy4.mf6.constants.FILL_DNODATA))
+rch_rate = np.full((nlay, nrow, ncol), flopy4.mf6.constants.FILL_DNODATA)
 rate = np.repeat(np.expand_dims(rch_rate, axis=0), repeats=nper, axis=0)
-rate[0, 0, :] = 3.0e-8
+rate[0, 0, ...] = 3.0e-8
 rch = flopy4.mf6.gwf.Rch(recharge=rate.reshape(nper, -1), dims=dims)
 
 # Output control
@@ -103,24 +118,24 @@
 # Wells scattered throughout the model
 wel_q = -5.0
 wel_nodes = [
-    [2, 4, 10, -5.0],
-    [1, 3, 5, -5.0],
-    [1, 5, 11, -5.0],
-    [0, 8, 7, -5.0],
-    [0, 8, 9, -5.0],
-    [0, 8, 11, -5.0],
-    [0, 8, 13, -5.0],
-    [0, 10, 7, -5.0],
-    [0, 10, 9, -5.0],
-    [0, 10, 11, -5.0],
-    [0, 10, 13, -5.0],
-    [0, 12, 7, -5.0],
-    [0, 12, 9, -5.0],
-    [0, 12, 11, -5.0],
-    [0, 12, 13, -5.0],
+    [2, 4, 10],
+    [1, 3, 5],
+    [1, 5, 11],
+    [0, 8, 7],
+    [0, 8, 9],
+    [0, 8, 11],
+    [0, 8, 13],
+    [0, 10, 7],
+    [0, 10, 9],
+    [0, 10, 11],
+    [0, 10, 13],
+    [0, 12, 7],
+    [0, 12, 9],
+    [0, 12, 11],
+    [0, 12, 13],
 ]
 wel = flopy4.mf6.gwf.Wel(
-    q={"*": {(layer, row, col): wel_q for layer, row, col, wel_q in wel_nodes}},
+    q={"*": {(layer, row, col): wel_q for layer, row, col in wel_nodes}},
     dims=dims,
 )
 
@@ -181,4 +196,82 @@
 )
 
 # Plot head results
-head.isel(layer=0, time=0).plot.contourf()
+plot_head(head, workspace)
+
+# UPDATE SIM for array based inputs
+
+# update simulation with array based inputs
+LAYER_NODATA = np.full((nrow, ncol), flopy4.mf6.constants.FILL_DNODATA, dtype=float)
+GRID_NODATA = np.full((nlay, nrow, ncol), flopy4.mf6.constants.FILL_DNODATA, dtype=float)
+
+head = np.repeat(np.expand_dims(GRID_NODATA, axis=0), repeats=nper, axis=0)
+for i in range(nrow):
+    for k in range(nlay - 1):
+        head[0, k, i, 0] = 0.0
+chdg = flopy4.mf6.gwf.Chdg(
+    print_input=True,
+    print_flows=True,
+    save_flows=True,
+    head=head.reshape(nper, -1),
+    dims=dims,
+)
+
+# Drain in the center left of the model
+elev = np.repeat(np.expand_dims(GRID_NODATA, axis=0), repeats=nper, axis=0)
+cond = np.repeat(np.expand_dims(GRID_NODATA, axis=0), repeats=nper, axis=0)
+for j in range(9):
+    elev[0, 0, 7, j + 1] = elevation[j]
+    cond[0, 0, 7, j + 1] = conductance
+drng = flopy4.mf6.gwf.Drng(
+    print_input=True,
+    print_flows=True,
+    save_flows=True,
+    elev=elev.reshape(nper, -1),
+    cond=cond.reshape(nper, -1),
+    dims=dims,
+)
+
+# well
+q = np.repeat(np.expand_dims(GRID_NODATA, axis=0), repeats=nper, axis=0)
+for layer, row, col in wel_nodes:
+    q[0, layer, row, col] = wel_q
+welg = flopy4.mf6.gwf.Welg(
+    q=q.reshape(nper, -1),
+    dims=dims,
+)
+
+# recharge
+recharge = np.repeat(np.expand_dims(LAYER_NODATA, axis=0), repeats=nper, axis=0)
+recharge[0, ...] = 3.0e-8
+rcha = flopy4.mf6.gwf.Rcha(recharge=recharge.reshape(nper, -1), dims=dims)
+
+# remove list based inputs
+# TODO: show variations on removing packages
+gwf.chd.remove(chd)
+del gwf.drn[0]
+del gwf.wel[0]
+del gwf.rch[0]
+
+# add array based inputs
+# TODO: needs type checking and list consolidation support (see comments in gwf init)
+gwf.chd = [chdg]
+gwf.drng = [drng]
+gwf.welg = [welg]
+gwf.rcha = [rcha]
+
+# create new workspace
+workspace = Path(__file__).parent / "twri2"
+workspace.mkdir(parents=True, exist_ok=True)
+sim.workspace = workspace
+
+sim.write()
+sim.run()
+
+# Load head results
+head = flopy4.mf6.utils.open_hds(
+    workspace / f"{gwf.name}.hds",
+    workspace / f"{gwf.name}.dis.grb",
+)
+
+# Plot head results
+plot_head(head, workspace)
diff --git a/flopy4/mf6/binding.py b/flopy4/mf6/binding.py
index a3c05fa9..e84009df 100644
--- a/flopy4/mf6/binding.py
+++ b/flopy4/mf6/binding.py
@@ -33,6 +33,8 @@ def _get_binding_type(component: Component) -> str:
             elif isinstance(component, Solution):
                 return f"{component.slntype}6"
             else:
+                if len(cls_name) == 4 and (cls_name[3] == "g" or cls_name[3] == "a"):
+                    return f"{cls_name[0:3].upper()}6"
                 return f"{cls_name.upper()}6"
 
         def _get_binding_terms(component: Component) -> tuple[str, ...] | None:
diff --git a/flopy4/mf6/codec/writer/filters.py b/flopy4/mf6/codec/writer/filters.py
index 69010135..bf9fa0ac 100644
--- a/flopy4/mf6/codec/writer/filters.py
+++ b/flopy4/mf6/codec/writer/filters.py
@@ -10,7 +10,7 @@
 
 from flopy4.mf6.constants import FILL_DNODATA
 
-ArrayHow = Literal["constant", "internal", "external", "layered internal"]
+ArrayHow = Literal["constant", "internal", "external", "layered constant", "layered internal"]
 
 
 def array_how(value: xr.DataArray) -> ArrayHow:
@@ -29,6 +29,14 @@ def array_how(value: xr.DataArray) -> ArrayHow:
     if value.ndim <= 2:
         return "internal"
     if value.ndim == 3:
+        layer_const = True
+        for layer in range(value.shape[0]):
+            val_layer = value.isel(nlay=layer)
+            if val_layer.max() != val_layer.min():
+                layer_const = False
+                break
+        if layer_const:
+            return "layered constant"
         return "layered internal"
     raise ValueError(f"Arrays with ndim > 3 are not supported, got ndim={value.ndim}")
 
diff --git a/flopy4/mf6/codec/writer/templates/macros.jinja b/flopy4/mf6/codec/writer/templates/macros.jinja
index 6f6e1a7e..df572206 100644
--- a/flopy4/mf6/codec/writer/templates/macros.jinja
+++ b/flopy4/mf6/codec/writer/templates/macros.jinja
@@ -35,7 +35,7 @@
 {{ inset }}CONSTANT {{ value|array2const }}
 {% elif how == "layered constant" %}
 {% for layer in value -%}
-{{ inset }}CONSTANT {{ layer|array2const }}
+{{ "\n" ~ inset }}CONSTANT {{ layer|array2const }}
 {%- endfor %}
 {% elif how == "layered internal" %}
 {% for layer in value %}
diff --git a/flopy4/mf6/constants.py b/flopy4/mf6/constants.py
index a1b3b89e..e053bfed 100644
--- a/flopy4/mf6/constants.py
+++ b/flopy4/mf6/constants.py
@@ -2,5 +2,5 @@
 
 MF6 = "mf6"
 FILL_DEFAULT = np.nan
-FILL_DNODATA = 1e30
+FILL_DNODATA = 3e30
 LENBOUNDNAME = 40
diff --git a/flopy4/mf6/context.py b/flopy4/mf6/context.py
index 2f31e4a0..08c64b9f 100644
--- a/flopy4/mf6/context.py
+++ b/flopy4/mf6/context.py
@@ -10,9 +10,29 @@
 from flopy4.utils import to_path
 
 
+def update_child_attr(instance, attribute, new_value):
+    """
+    Generalized function to update child attribute (e.g. workspace).
+
+    Args:
+        instance: The model instance
+        attribute: The attribute being set (from attrs on_setattr)
+        new_value: The new value being set
+
+    Returns:
+        The new_value (unchanged)
+    """
+
+    for child in instance.children.values():  # type: ignore
+        if hasattr(child, attribute.name):
+            setattr(child, attribute.name, new_value)
+
+    return new_value
+
+
 @xattree
 class Context(Component, ABC):
-    workspace: Path = field(default=None, converter=to_path)
+    workspace: Path = field(default=None, converter=to_path, on_setattr=update_child_attr)
 
     def __attrs_post_init__(self):
         super().__attrs_post_init__()
diff --git a/flopy4/mf6/converter/unstructure.py b/flopy4/mf6/converter/unstructure.py
index 76a3f95f..652d1090 100644
--- a/flopy4/mf6/converter/unstructure.py
+++ b/flopy4/mf6/converter/unstructure.py
@@ -7,9 +7,11 @@
 import xarray as xr
 import xattree
 from modflow_devtools.dfn.schema.block import block_sort_key
+from xattree import XatSpec
 
 from flopy4.mf6.binding import Binding
 from flopy4.mf6.component import Component
+from flopy4.mf6.constants import FILL_DNODATA
 from flopy4.mf6.context import Context
 from flopy4.mf6.spec import FileInOut
 
@@ -144,9 +146,82 @@ def oc_setting_data(rec):
     return data
 
 
+def _unstructure_block_param(
+    block_name: str,
+    field_name: str,
+    xatspec: XatSpec,
+    value: Component,
+    data: dict[str, Any],
+    blocks: dict,
+    period_data: dict,
+) -> None:
+    # Skip child components that have been processed as bindings
+    if isinstance(value, Context) and field_name in xatspec.children:
+        child_spec = xatspec.children[field_name]
+        if hasattr(child_spec, "metadata") and "block" in child_spec.metadata:  # type: ignore
+            if child_spec.metadata["block"] == block_name:  # type: ignore
+                return
+
+    # filter out empty values and false keywords, and convert:
+    #   - paths to records
+    #   - datetimes to ISO format
+    #   - filter out false keywords
+    #   - 'auxiliary' fields to tuples
+    #   - xarray DataArrays with 'nper' dim to dict of kper-sliced datasets
+    #   - other values to their original form
+    # TODO: use cattrs converters for field unstructuring?
+    match field_value := data[field_name]:
+        case None:
+            return
+        case bool():
+            if field_value:
+                blocks[block_name][field_name] = field_value
+        case Path():
+            field_spec = xatspec.attrs[field_name]
+            field_meta = getattr(field_spec, "metadata", {})
+            t = _path_to_tuple(field_name, field_value, inout=field_meta.get("inout", "fileout"))
+            # name may have changed e.g dropping '_file' suffix
+            blocks[block_name][t[0]] = t
+        case datetime():
+            blocks[block_name][field_name] = field_value.isoformat()
+        case t if (
+            field_name == "auxiliary" and hasattr(field_value, "values") and field_value is not None
+        ):
+            blocks[block_name][field_name] = tuple(field_value.values.tolist())
+        case xr.DataArray():
+            has_spatial_dims = any(
+                dim in field_value.dims for dim in ["nlay", "nrow", "ncol", "ncpl", "nodes"]
+            )
+            if has_spatial_dims:
+                field_value = _hack_structured_grid_dims(
+                    field_value,
+                    structured_grid_dims=value.data.dims,  # type: ignore
+                )
+            if "nper" in field_value.dims and block_name == "period":
+                if not np.issubdtype(field_value.dtype, np.number):
+                    dat = _hack_period_non_numeric(field_name, field_value)
+                    for n, v in dat.items():
+                        period_data[n] = v
+                else:
+                    period_data[field_name] = {
+                        kper: field_value.isel(nper=kper)  # type: ignore
+                        for kper in range(field_value.sizes["nper"])
+                    }
+            else:
+                blocks[block_name][field_name] = field_value
+        case _:
+            blocks[block_name][field_name] = field_value
+
+
 def unstructure_component(value: Component) -> dict[str, Any]:
-    from flopy4.mf6.constants import FILL_DNODATA
+    xatspec = xattree.get_xatspec(type(value))
+    if "readarraygrid" in xatspec.attrs or "readasarrays" in xatspec.attrs:
+        return _unstructure_array_component(value)
+    else:
+        return _unstructure_component(value)
 
+
+def _unstructure_array_component(value: Component) -> dict[str, Any]:
     blockspec = dict(sorted(value.dfn.blocks.items(), key=block_sort_key))  # type: ignore
     blocks: dict[str, dict[str, Any]] = {}
     xatspec = xattree.get_xatspec(type(value))
@@ -161,73 +236,57 @@ def unstructure_component(value: Component) -> dict[str, Any]:
     for block_name, block in blockspec.items():
         period_data = {}  # type: ignore
         period_blocks = {}  # type: ignore
-        period_block_name = None
 
         if block_name not in blocks:
             blocks[block_name] = {}
 
         for field_name in block.keys():
-            # Skip child components that have been processed as bindings
-            if isinstance(value, Context) and field_name in xatspec.children:
-                child_spec = xatspec.children[field_name]
-                if hasattr(child_spec, "metadata") and "block" in child_spec.metadata:  # type: ignore
-                    if child_spec.metadata["block"] == block_name:  # type: ignore
-                        continue
-
-            # filter out empty values and false keywords, and convert:
-            #   - paths to records
-            #   - datetimes to ISO format
-            #   - filter out false keywords
-            #   - 'auxiliary' fields to tuples
-            #   - xarray DataArrays with 'nper' dim to dict of kper-sliced datasets
-            #   - other values to their original form
-            # TODO: use cattrs converters for field unstructuring?
-            match field_value := data[field_name]:
-                case None:
-                    continue
-                case bool():
-                    if field_value:
-                        blocks[block_name][field_name] = field_value
-                case Path():
-                    field_spec = xatspec.attrs[field_name]
-                    field_meta = getattr(field_spec, "metadata", {})
-                    t = _path_to_tuple(
-                        field_name, field_value, inout=field_meta.get("inout", "fileout")
-                    )
-                    # name may have changed e.g dropping '_file' suffix
-                    blocks[block_name][t[0]] = t
-                case datetime():
-                    blocks[block_name][field_name] = field_value.isoformat()
-                case t if (
-                    field_name == "auxiliary"
-                    and hasattr(field_value, "values")
-                    and field_value is not None
-                ):
-                    blocks[block_name][field_name] = tuple(field_value.values.tolist())
-                case xr.DataArray() if "nper" in field_value.dims:
-                    has_spatial_dims = any(
-                        dim in field_value.dims for dim in ["nlay", "nrow", "ncol", "nodes"]
-                    )
-                    if has_spatial_dims:
-                        field_value = _hack_structured_grid_dims(
-                            field_value,
-                            structured_grid_dims=value.parent.data.dims,  # type: ignore
-                        )
-                    if block_name == "period":
-                        if not np.issubdtype(field_value.dtype, np.number):
-                            dat = _hack_period_non_numeric(field_name, field_value)
-                            for n, v in dat.items():
-                                period_data[n] = v
-                        else:
-                            period_data[field_name] = {
-                                kper: field_value.isel(nper=kper)  # type: ignore
-                                for kper in range(field_value.sizes["nper"])
-                            }
-                    else:
-                        blocks[block_name][field_name] = field_value
-
-                case _:
-                    blocks[block_name][field_name] = field_value
+            _unstructure_block_param(
+                block_name, field_name, xatspec, value, data, blocks, period_data
+            )
+
+        # invert key order, (arr_name, kper) -> (kper, arr_name)
+        for arr_name, periods in period_data.items():
+            for kper, arr in periods.items():
+                if kper not in period_blocks:
+                    period_blocks[kper] = {}
+                period_blocks[kper][arr_name] = arr
+
+        # setup indexed period blocks, combine arrays into datasets
+        for kper, block in period_blocks.items():
+            key = f"period {kper + 1}"
+            for arr_name, val in block.items():
+                if not np.all(val == FILL_DNODATA):
+                    if key not in blocks:
+                        blocks[key] = {}
+                    blocks[key][arr_name] = val
+
+    return {name: block for name, block in blocks.items() if name != "period"}
+
+
+def _unstructure_component(value: Component) -> dict[str, Any]:
+    blockspec = dict(sorted(value.dfn.blocks.items(), key=block_sort_key))  # type: ignore
+    blocks: dict[str, dict[str, Any]] = {}
+    xatspec = xattree.get_xatspec(type(value))
+    data = xattree.asdict(value)
+
+    # create child component binding blocks
+    blocks.update(_make_binding_blocks(value))
+
+    # process blocks in order, unstructuring fields as needed,
+    # then slice period data into separate kper-indexed blocks
+    # each of which contains a dataset indexed for that period.
+    for block_name, block in blockspec.items():
+        period_data = {}  # type: ignore
+        period_blocks = {}  # type: ignore
+
+        if block_name not in blocks:
+            blocks[block_name] = {}
+
+        for field_name in block.keys():
+            _unstructure_block_param(
+                block_name, field_name, xatspec, value, data, blocks, period_data
+            )
 
         # invert key order, (arr_name, kper) -> (kper, arr_name)
         for arr_name, periods in period_data.items():
@@ -267,6 +326,10 @@ def unstructure_component(value: Component) -> dict[str, Any]:
         if perioddata := blocks.get("perioddata", None):
             blocks["perioddata"] = {"perioddata": xr.Dataset(perioddata)}
 
+    # TODO: this fixes out of order blocks (e.g. model namefile) from
+    # blocks.update() child binding call above
+    blocks = dict(sorted(blocks.items(), key=block_sort_key))
+
     # total temporary hack! manually set solutiongroup 1.
     # TODO still need to support multiple..
     if "solutiongroup" in blocks:
diff --git a/flopy4/mf6/gwf/__init__.py b/flopy4/mf6/gwf/__init__.py
index 8c69b9a1..e1a59736 100644
--- a/flopy4/mf6/gwf/__init__.py
+++ b/flopy4/mf6/gwf/__init__.py
@@ -8,20 +8,39 @@
 from xattree import xattree
 
 from flopy4.mf6.gwf.chd import Chd
+from flopy4.mf6.gwf.chdg import Chdg
 from flopy4.mf6.gwf.dis import Dis
 from flopy4.mf6.gwf.drn import Drn
+from flopy4.mf6.gwf.drng import Drng
 from flopy4.mf6.gwf.ic import Ic
 from flopy4.mf6.gwf.npf import Npf
 from flopy4.mf6.gwf.oc import Oc
 from flopy4.mf6.gwf.rch import Rch
+from flopy4.mf6.gwf.rcha import Rcha
 from flopy4.mf6.gwf.sto import Sto
 from flopy4.mf6.gwf.wel import Wel
+from flopy4.mf6.gwf.welg import Welg
 from flopy4.mf6.model import Model
 from flopy4.mf6.spec import field, path
 from flopy4.mf6.utils import open_cbc, open_hds
 from flopy4.utils import to_path
 
-__all__ = ["Gwf", "Chd", "Dis", "Drn", "Ic", "Npf", "Oc", "Sto", "Wel", "Rch"]
+__all__ = [
+    "Gwf",
+    "Chd",
+    "Chdg",
+    "Dis",
+    "Drn",
+    "Drng",
+    "Ic",
+    "Npf",
+    "Oc",
+    "Rch",
+    "Rcha",
+    "Sto",
+    "Wel",
+    "Welg",
+]
 
 
 def convert_grid(value):
@@ -78,10 +97,18 @@ def budget(self):
     oc: Oc | None = field(block="packages", default=None)
     npf: Npf | None = field(block="packages", default=None)
     sto: Sto | None = field(block="packages", default=None)
+    # TODO: implement type check for all lists (and singletons?)
     chd: list[Chd] = field(block="packages")
-    wel: list[Wel] = field(block="packages")
+    chdg: list[Chdg] = field(block="packages")
+    # TODO: consolidate all package flavors to single list
+    #       based on hydrologic feature
+    # chd: List[Union[Chd, Chdg]] = field(block="packages")
     drn: list[Drn] = field(block="packages")
+    drng: list[Drng] = field(block="packages")
     rch: list[Rch] = field(block="packages")
+    rcha: list[Rcha] = field(block="packages")
+    wel: list[Wel] = field(block="packages")
+    welg: list[Welg] = field(block="packages")
     output: Output = attrs.field(
         default=attrs.Factory(lambda self: Gwf.Output(self), takes_self=True)
     )
diff --git a/flopy4/mf6/gwf/chdg.py b/flopy4/mf6/gwf/chdg.py
new file mode 100644
index 00000000..943bbfc0
--- /dev/null
+++ b/flopy4/mf6/gwf/chdg.py
@@ -0,0 +1,46 @@
+from pathlib import Path
+from typing import ClassVar, Optional
+
+import numpy as np
+from numpy.typing import NDArray
+from xattree import xattree
+
+from flopy4.mf6.package import Package
+from flopy4.mf6.spec import array, field, path
+from flopy4.mf6.utils.grid import update_maxbound
+from flopy4.utils import to_path
+
+
+@xattree
+class Chdg(Package):
+    multi_package: ClassVar[bool] = True
+    auxiliary: Optional[list[str]] = array(block="options", default=None)
+    auxmultname: Optional[str] = field(block="options", default=None)
+    print_input: bool = field(block="options", default=False)
+    print_flows: bool = field(block="options", default=False)
+    readarraygrid: bool = field(block="options", default=True)
+    save_flows: bool = field(block="options", default=False)
+    obs_filerecord: Optional[Path] = path(
+        block="options", default=None, converter=to_path, inout="fileout"
+    )
+    export_array_netcdf: bool = field(block="options", default=False)
+    dev_no_newton: bool = field(default=False, block="options")
+    maxbound: Optional[int] = field(block="dimensions", default=None, init=False)
+    head: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "nodes",
+        ),
+        default=None,
+        on_setattr=update_maxbound,
+    )
+    aux: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "nodes",
+        ),
+        default=None,
+        on_setattr=update_maxbound,
+    )
diff --git a/flopy4/mf6/gwf/dis.py b/flopy4/mf6/gwf/dis.py
index cf8c4416..7b42bbb4 100644
--- a/flopy4/mf6/gwf/dis.py
+++ b/flopy4/mf6/gwf/dis.py
@@ -75,9 +75,15 @@ class Dis(Package):
         scope="gwf",
         init=False,
     )
+    ncpl: int = dim(
+        coord="lnode",
+        scope="gwf",
+        init=False,
+    )
 
     def __attrs_post_init__(self):
         self.nodes = self.ncol * self.nrow * self.nlay
+        self.ncpl = self.ncol * self.nrow
         super().__attrs_post_init__()
 
     def to_grid(self) -> StructuredGrid:
diff --git a/flopy4/mf6/gwf/drng.py b/flopy4/mf6/gwf/drng.py
new file mode 100644
index 00000000..68b579f5
--- /dev/null
+++ b/flopy4/mf6/gwf/drng.py
@@ -0,0 +1,56 @@
+from pathlib import Path
+from typing import ClassVar, Optional
+
+import numpy as np
+from numpy.typing import NDArray
+from xattree import xattree
+
+from flopy4.mf6.package import Package
+from flopy4.mf6.spec import array, field, path
+from flopy4.mf6.utils.grid import update_maxbound
+from flopy4.utils import to_path
+
+
+@xattree
+class Drng(Package):
+    multi_package: ClassVar[bool] = True
+    auxiliary: Optional[list[str]] = array(block="options", default=None)
+    auxmultname: Optional[str] = field(block="options", default=None)
+    print_input: bool = field(block="options", default=False)
+    print_flows: bool = field(block="options", default=False)
+    readarraygrid: bool = field(block="options", default=True)
+    save_flows: bool = field(block="options", default=False)
+    obs_filerecord: Optional[Path] = path(
+        block="options", default=None, converter=to_path, inout="fileout"
+    )
+    export_array_netcdf: bool = field(block="options", default=False)
+    mover: bool = field(block="options", default=False)
+    dev_cubic_scaling: bool = field(default=False, block="options")
+    maxbound: Optional[int] = field(block="dimensions", default=None, init=False)
+    elev: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "nodes",
+        ),
+        default=None,
+        on_setattr=update_maxbound,
+    )
+    cond: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "nodes",
+        ),
+        default=None,
+        on_setattr=update_maxbound,
+    )
+    aux: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "nodes",
+        ),
+        default=None,
+        on_setattr=update_maxbound,
+    )
diff --git a/flopy4/mf6/gwf/rcha.py b/flopy4/mf6/gwf/rcha.py
new file mode 100644
index 00000000..2e8642e1
--- /dev/null
+++ b/flopy4/mf6/gwf/rcha.py
@@ -0,0 +1,52 @@
+from pathlib import Path
+from typing import ClassVar, Optional
+
+import numpy as np
+from numpy.typing import NDArray
+from xattree import xattree
+
+from flopy4.mf6.package import Package
+from flopy4.mf6.spec import array, field, path
+from flopy4.utils import to_path
+
+
+@xattree
+class Rcha(Package):
+    multi_package: ClassVar[bool] = True
+    fixed_cell: bool = field(block="options", default=False)
+    auxiliary: Optional[list[str]] = array(block="options", default=None)
+    auxmultname: Optional[str] = field(block="options", default=None)
+    print_input: bool = field(block="options", default=False)
+    print_flows: bool = field(block="options", default=False)
+    readasarrays: bool = field(block="options", default=True)
+    save_flows: bool = field(block="options", default=False)
+    tas_filerecord: Optional[Path] = path(
+        block="options", default=None, converter=to_path, inout="filein"
+    )
+    obs_filerecord: Optional[Path] = path(
+        block="options", default=None, converter=to_path, inout="fileout"
+    )
+    irch: Optional[NDArray[np.int64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "ncpl",
+        ),
+        default=None,
+    )
+    recharge: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "ncpl",
+        ),
+        default=None,
+    )
+    aux: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "ncpl",
+        ),
+        default=None,
+    )
diff --git a/flopy4/mf6/gwf/welg.py b/flopy4/mf6/gwf/welg.py
new file mode 100644
index 00000000..438b03b3
--- /dev/null
+++ b/flopy4/mf6/gwf/welg.py
@@ -0,0 +1,52 @@
+from pathlib import Path
+from typing import ClassVar, Optional
+
+import numpy as np
+from numpy.typing import NDArray
+from xattree import xattree
+
+from flopy4.mf6.package import Package
+from flopy4.mf6.spec import array, field, path
+from flopy4.mf6.utils.grid import update_maxbound
+from flopy4.utils import to_path
+
+
+@xattree
+class Welg(Package):
+    multi_package: ClassVar[bool] = True
+    auxiliary: Optional[list[str]] = array(block="options", default=None)
+    auxmultname: Optional[str] = field(block="options", default=None)
+    print_input: bool = field(block="options", default=False)
+    print_flows: bool = field(block="options", default=False)
+    readarraygrid: bool = field(block="options", default=True)
+    save_flows: bool = field(block="options", default=False)
+    auto_flow_reduce: float = field(block="options", default=None)
+    afrcsv_filerecord: Optional[Path] = path(
+        block="options", default=None, converter=to_path, inout="fileout"
+    )
+    ts_filerecord: Optional[Path] = path(
+        block="options", default=None, converter=to_path, inout="filein"
+    )
+    obs_filerecord: Optional[Path] = path(
+        block="options", default=None, converter=to_path, inout="fileout"
+    )
+    mover: bool = field(block="options", default=False)
+    maxbound: Optional[int] = field(block="dimensions", default=None, init=False)
+    q: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "nodes",
+        ),
+        default=None,
+        on_setattr=update_maxbound,
+    )
+    aux: Optional[NDArray[np.float64]] = array(
+        block="period",
+        dims=(
+            "nper",
+            "nodes",
+        ),
+        default=None,
+        on_setattr=update_maxbound,
+    )
diff --git a/flopy4/mf6/simulation.py b/flopy4/mf6/simulation.py
index 2ddebde5..aa8a8a16 100644
--- a/flopy4/mf6/simulation.py
+++ b/flopy4/mf6/simulation.py
@@ -4,7 +4,7 @@
 from modflow_devtools.misc import cd, run_cmd
 from xattree import xattree
 
-from flopy4.mf6.context import Context
+from flopy4.mf6.context import Context, update_child_attr
 from flopy4.mf6.exchange import Exchange
 from flopy4.mf6.model import Model
 from flopy4.mf6.solution import Solution
@@ -32,6 +32,8 @@ def default_filename(self) -> str:
         return "mfsim.nam"
 
     def __attrs_post_init__(self):
+        from attrs import fields_dict
+
         super().__attrs_post_init__()
         if self.filename != "mfsim.nam":
             if self.filename is not None:
@@ -40,8 +42,9 @@ def __attrs_post_init__(self):
                     UserWarning,
                 )
             self.filename = "mfsim.nam"
-        for model in self.models.values():
-            model.workspace = self.workspace
+        fields = fields_dict(type(self))
+        field = fields["workspace"]
+        update_child_attr(self, field, self.workspace)
 
     @property
     def time(self) -> Time:
diff --git a/pixi.lock b/pixi.lock
index 70be5004..721b6978 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -104,7 +104,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/a7/c2/fe1e52489ae3122415c51f387e221dd0773709bad6c6cdaa599e8a2c5185/urllib3-2.5.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/50/bd/c336448be43d40be28e71f2e0f3caf7ccb28e2755c58f4c02c065bfe3e8e/WebOb-1.8.9-py2.py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -199,7 +199,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/a7/c2/fe1e52489ae3122415c51f387e221dd0773709bad6c6cdaa599e8a2c5185/urllib3-2.5.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/50/bd/c336448be43d40be28e71f2e0f3caf7ccb28e2755c58f4c02c065bfe3e8e/WebOb-1.8.9-py2.py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -297,7 +297,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/a7/c2/fe1e52489ae3122415c51f387e221dd0773709bad6c6cdaa599e8a2c5185/urllib3-2.5.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/50/bd/c336448be43d40be28e71f2e0f3caf7ccb28e2755c58f4c02c065bfe3e8e/WebOb-1.8.9-py2.py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -535,7 +535,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -755,7 +755,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -977,7 +977,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -1186,7 +1186,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -1379,7 +1379,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -1573,7 +1573,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -1780,7 +1780,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -1970,7 +1970,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -2161,7 +2161,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -2367,7 +2367,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -2560,7 +2560,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -2754,7 +2754,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/34/db/b10e48aa8fff7407e67470363eac595018441cf32d5e1001567a7aeba5d2/websocket_client-1.9.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/3f/0e/fa3b193432cfc60c93b42f3be03365f5f909d2b3ea410295cf36df739e31/widgetsnbextension-4.0.15-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c3/78/4d6d68555a92cb97b4c192759c4ab585c5cb23490f64d4ddf12c66a3b051/xarray-2025.10.1-py3-none-any.whl
-      - pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+      - pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
       - pypi: https://files.pythonhosted.org/packages/5b/8f/447cc9cb57456d786204af0f450ffb920039104c5eff6626337c9f403bd1/xyzservices-2025.10.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/1a/71/9de7229515a53d1cc5705ca9c411530f711a2242f962214d9dbfe2741aa4/zarr-3.1.3-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/80/ab/11a76c1e2126084fde2639514f24e6111b789b0bfa4fc6264a8975c7e1f1/zict-3.0.0-py2.py3-none-any.whl
@@ -10094,7 +10094,7 @@ packages:
   - types-requests ; extra == 'types'
   - types-setuptools ; extra == 'types'
   requires_python: '>=3.11'
-- pypi: git+https://github.com/wpbonelli/xattree.git#05778bdd17bec769c5bbd5672ae06bcd242f2c6a
+- pypi: git+https://github.com/wpbonelli/xattree.git#d77984d86795bf382c6ce7adc7ffa0904d61e79f
   name: xattree
   version: 0.1.0.dev0
   requires_dist:
diff --git a/test/test_mf6_codec.py b/test/test_mf6_codec.py
index 84db7e65..f0cde5c9 100644
--- a/test/test_mf6_codec.py
+++ b/test/test_mf6_codec.py
@@ -2,9 +2,11 @@
 
 from pprint import pprint
 
+import numpy as np
 import pytest
 
 from flopy4.mf6.codec import dumps, loads
+from flopy4.mf6.constants import FILL_DNODATA
 from flopy4.mf6.converter import COMPONENT_CONVERTER
 
 
@@ -249,14 +251,57 @@ def test_dumps_chd():
     assert len(lines) == 2
     assert "1 1 1 10.0" in dumped  # First CHD cell - node 1
     assert "1 10 10 20.0" in dumped  # Second CHD cell - node 100
-    assert "1e+30" not in dumped
-    assert "1.0e+30" not in dumped
+    assert "3e+30" not in dumped
+    assert "3.0e+30" not in dumped
 
     loaded = loads(dumped)
     print("CHD load:")
     pprint(loaded)
 
 
+def test_dumps_chdg():
+    from flopy4.mf6.gwf import Chdg, Dis, Gwf
+
+    nlay = 1
+    nrow = 10
+    ncol = 10
+
+    dis = Dis(nlay=nlay, nrow=nrow, ncol=ncol)
+    gwf = Gwf(dis=dis)
+
+    head = np.full((nlay, nrow, ncol), FILL_DNODATA, dtype=float)
+    head[0, 0, 0] = 1.0
+    head[0, 9, 9] = 0.0
+    chd = Chdg(
+        parent=gwf,
+        head=np.expand_dims(head.ravel(), axis=0),
+        save_flows=True,
+        print_input=True,
+        dims={"nper": 1},
+    )
+
+    dumped = dumps(COMPONENT_CONVERTER.unstructure(chd))
+    print("CHD dump:")
+    print(dumped)
+
+    assert "BEGIN PERIOD 1" in dumped
+    assert "END PERIOD 1" in dumped
+
+    period_section = dumped.split("BEGIN PERIOD 1")[1].split("END PERIOD 1")[0].strip()
+    lines = [line.strip() for line in period_section.split("\n") if line.strip()]
+
+    assert len(lines) == 12
+    assert "READARRAYGRID" in dumped
+    assert "MAXBOUND 2" in dumped
+    dump_data = [[float(x) for x in line.split()] for line in lines[2:12]]
+    dump_head = np.array(dump_data)
+    assert np.allclose(head, dump_head)
+
+    loaded = loads(dumped)
+    print("CHDG load:")
+    pprint(loaded)
+
+
 def test_dumps_wel():
     from flopy4.mf6.gwf import Dis, Gwf, Wel
 
@@ -291,8 +336,8 @@ def test_dumps_wel():
     assert "1 3 4 -100.0" in dumped  # (0,2,3) -> node 24
     assert "2 6 8 -50.0" in dumped  # (1,5,7) -> node 158
     assert "3 9 2 25.0" in dumped  # (2,8,1) -> node 282
-    assert "1e+30" not in dumped
-    assert "1.0e+30" not in dumped
+    assert "3e+30" not in dumped
+    assert "3.0e+30" not in dumped
 
     loaded = loads(dumped)
     print("WEL load:")
@@ -356,8 +401,8 @@ def test_dumps_drn():
     assert "1 2 2 12.0 1.5" in dumped  # Period 2: (0,1,1)
     assert "1 3 4 9.0 0.8" in dumped  # Period 2: (0,2,3)
     assert "2 4 3 7.0 2.2" in dumped  # Period 2: (1,3,2)
-    assert "1e+30" not in dumped
-    assert "1.0e+30" not in dumped
+    assert "3e+30" not in dumped
+    assert "3.0e+30" not in dumped
 
     loaded = loads(dumped)
     print("DRN load:")
@@ -369,9 +414,9 @@ def test_dumps_npf():
 
     dis = Dis(nlay=2, nrow=5, ncol=5)
     gwf = Gwf(dis=dis)
-    drn = Npf(parent=gwf, cvoptions=Npf.CvOptions(dewatered=True), k=1.0)
+    npf = Npf(parent=gwf, cvoptions=Npf.CvOptions(dewatered=True), k=1.0)
 
-    dumped = dumps(COMPONENT_CONVERTER.unstructure(drn))
+    dumped = dumps(COMPONENT_CONVERTER.unstructure(npf))
     print("NPF dump:")
     print(dumped)
 
@@ -407,8 +452,8 @@ def test_dumps_chd_2():
     assert "100.0" in dumped  # Left boundary
     assert "95.0" in dumped  # Right boundary
     assert "98.0" in dumped  # Bottom boundary
-    assert "1e+30" not in dumped
-    assert "1.0e+30" not in dumped
+    assert "3e+30" not in dumped
+    assert "3.0e+30" not in dumped
 
     loaded = loads(dumped)
     print("CHD load:")
@@ -450,8 +495,8 @@ def test_dumps_wel_with_aux():
     # node q aux_value
     assert "1 2 3 -75.0 1.0" in dumped  # (0,1,2) -> node 8, q=-75.0, aux=1.0
     assert "2 4 5 -25.0 2.0" in dumped  # (1,3,4) -> node 45, q=-25.0, aux=2.0
-    assert "1e+30" not in dumped
-    assert "1.0e+30" not in dumped
+    assert "3e+30" not in dumped
+    assert "3.0e+30" not in dumped
 
     loaded = loads(dumped)
     print("WEL+aux load:")
diff --git a/test/test_mf6_component.py b/test/test_mf6_component.py
index 03cd9ff1..28c62e5c 100644
--- a/test/test_mf6_component.py
+++ b/test/test_mf6_component.py
@@ -11,7 +11,7 @@
 
 from flopy4.mf6.component import COMPONENTS
 from flopy4.mf6.constants import FILL_DNODATA, LENBOUNDNAME
-from flopy4.mf6.gwf import Chd, Dis, Gwf, Ic, Npf, Oc
+from flopy4.mf6.gwf import Chd, Chdg, Dis, Gwf, Ic, Npf, Oc
 from flopy4.mf6.ims import Ims
 from flopy4.mf6.simulation import Simulation
 from flopy4.mf6.tdis import Tdis
@@ -398,6 +398,57 @@ def test_quickstart(function_tmpdir):
     sim.run()
 
 
+def test_quickstart_grid(function_tmpdir):
+    sim_name = "quickstart"
+    gwf_name = "mymodel"
+
+    # dimensions
+    nlay = 1
+    nrow = 10
+    ncol = 10
+    nstp = 1
+
+    time = Time(perlen=[1.0], nstp=[1], tsmult=[1.0])
+    # time = Time(perlen=[1.0, 1.0], nstp=[1, 1], tsmult=[1.0, 1.0])
+    ims = Ims(models=[gwf_name])
+    dis = Dis(
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        top=1.0,
+        botm=0.0,
+    )
+    sim = Simulation(
+        tdis=time,
+        workspace=function_tmpdir,
+        name=sim_name,
+        solutions={"ims": ims},
+    )
+    gwf = Gwf(parent=sim, dis=dis, name=gwf_name)
+    ic = Ic(parent=gwf)
+    oc = Oc(
+        parent=gwf,
+        budget_file=f"{gwf_name}.bud",
+        head_file=f"{gwf_name}.hds",
+        save_head=["all"],
+        save_budget=["all"],
+    )
+    npf = Npf(parent=gwf, icelltype=0, k=1.0)
+
+    # chd grid based input, step 1 data head array
+    head = np.full((nlay, nrow, ncol), FILL_DNODATA, dtype=float)
+    head[0, 0, 0] = 1.0
+    head[0, 9, 9] = 0.0
+    # TODO: support dict style input keyed on SP with separate grid arrays
+    chd = Chdg(parent=gwf, head=np.expand_dims(head.ravel(), axis=0))
+    # headnone = np.full((nlay, nrow, ncol), FILL_DNODATA, dtype=float)
+    # ts_head = np.stack((head.ravel(), headnone.ravel()), axis=0)
+    # chd = Chdg(parent=gwf, head=ts_head)
+
+    sim.write()
+    sim.run()
+
+
 def test_write_ascii(function_tmpdir):
     sim_name = "sim"
     gwf_name = "gwf"