diff --git a/docs/examples/array_example.py b/docs/examples/array_example.py
index afc641f..e7d5e41 100644
--- a/docs/examples/array_example.py
+++ b/docs/examples/array_example.py
@@ -45,11 +45,12 @@
 constant = data_path / "constant.txt"
 external = data_path / "external.txt"
 shape = (1000, 100)
+type = "double"
 
 # Open and load a NumPy array representation
 
 fhandle = open(internal)
-imfa = MFArray.load(fhandle, data_path, shape, header=False)
+imfa = MFArray.load(fhandle, data_path, shape, type=type, header=False)
 
 # Get values
 
@@ -70,7 +71,7 @@
 plt.colorbar()
 
 fhandle = open(constant)
-cmfa = MFArray.load(fhandle, data_path, shape, header=False)
+cmfa = MFArray.load(fhandle, data_path, shape, type=type, header=False)
 cvals = cmfa.value
 plt.imshow(cvals[0:100])
 plt.colorbar()
@@ -93,7 +94,7 @@
 # External
 
 fhandle = open(external)
-emfa = MFArray.load(fhandle, data_path, shape, header=False)
+emfa = MFArray.load(fhandle, data_path, shape, type=type, header=False)
 evals = emfa.value
 evals
 
@@ -118,7 +119,9 @@
 
 fhandle = open(ilayered)
 shape = (3, 1000, 100)
-ilmfa = MFArray.load(fhandle, data_path, shape, header=False, layered=True)
+ilmfa = MFArray.load(
+    fhandle, data_path, shape, type=type, header=False, layered=True
+)
 vals = ilmfa.value
 
 ilmfa._value  # internal storage
@@ -165,7 +168,9 @@
 
 fhandle = open(clayered)
 shape = (3, 1000, 100)
-clmfa = MFArray.load(fhandle, data_path, shape, header=False, layered=True)
+clmfa = MFArray.load(
+    fhandle, data_path, shape, type=type, header=False, layered=True
+)
 
 clmfa._value
 
@@ -218,7 +223,9 @@
 
 fhandle = open(mlayered)
 shape = (3, 1000, 100)
-mlmfa = MFArray.load(fhandle, data_path, shape, header=False, layered=True)
+mlmfa = MFArray.load(
+    fhandle, data_path, shape, type=type, header=False, layered=True
+)
 
 mlmfa.how
 
diff --git a/flopy4/array.py b/flopy4/array.py
index 775ec16..4b80781 100644
--- a/flopy4/array.py
+++ b/flopy4/array.py
@@ -434,6 +434,16 @@ def load(cls, f, cwd, shape, header=True, **kwargs):
         model_shape = kwargs.pop("model_shape", None)
         params = kwargs.pop("blk_params", {})
         mempath = kwargs.pop("mempath", None)
+        atype = kwargs.get("type", None)
+
+        if atype is not None:
+            if atype == "integer":
+                dtype = np.int32
+            elif atype == "double":
+                dtype = np.float64
+        else:
+            raise ValueError("array spec type not defined")
+
         if model_shape and isinstance(shape, str):
             if shape == "(nodes)":
                 n = math.prod([x for x in model_shape])
@@ -459,7 +469,7 @@ def load(cls, f, cwd, shape, header=True, **kwargs):
             lshp = shape[1:]
             objs = []
             for _ in range(nlay):
-                mfa = cls._load(f, cwd, lshp, name)
+                mfa = cls._load(f, cwd, lshp, dtype=dtype, name=name)
                 objs.append(mfa)
 
             return MFArray(
@@ -474,11 +484,17 @@ def load(cls, f, cwd, shape, header=True, **kwargs):
         else:
             kwargs.pop("layered", None)
             return cls._load(
-                f, cwd, shape, layered=layered, name=name, **kwargs
+                f,
+                cwd,
+                shape,
+                layered=layered,
+                dtype=dtype,
+                name=name,
+                **kwargs,
             )
 
     @classmethod
-    def _load(cls, f, cwd, shape, layered=False, **kwargs):
+    def _load(cls, f, cwd, shape, layered=False, dtype=None, **kwargs):
         control_line = multi_line_strip(f).split()
 
         if CommonNames.iprn.lower() in control_line:
@@ -491,17 +507,20 @@ def _load(cls, f, cwd, shape, layered=False, **kwargs):
         clpos = 1
 
         if how == MFArrayType.internal:
-            array = cls.read_array(f)
+            array = cls.read_array(f, dtype)
 
         elif how == MFArrayType.constant:
-            array = float(control_line[clpos])
+            if dtype == np.float64:
+                array = float(control_line[clpos])
+            else:
+                array = int(control_line[clpos])
             clpos += 1
 
         elif how == MFArrayType.external:
             extpath = Path(control_line[clpos])
             fpath = cwd / extpath
             with open(fpath) as foo:
-                array = cls.read_array(foo)
+                array = cls.read_array(foo, dtype)
             clpos += 1
 
         else:
@@ -509,7 +528,10 @@ def _load(cls, f, cwd, shape, layered=False, **kwargs):
 
         factor = None
         if len(control_line) > 2:
-            factor = float(control_line[clpos + 1])
+            if dtype == np.float64:
+                factor = float(control_line[clpos + 1])
+            else:
+                factor = int(control_line[clpos + 1])
 
         return cls(
             shape,
@@ -521,7 +543,7 @@ def _load(cls, f, cwd, shape, layered=False, **kwargs):
         )
 
     @staticmethod
-    def read_array(f):
+    def read_array(f, dtype):
         """
         Read a MODFLOW 6 array from an open file
         into a flat NumPy array representation.
@@ -538,5 +560,5 @@ def read_array(f):
             astr.append(line)
 
         astr = StringIO(" ".join(astr))
-        array = np.genfromtxt(astr).ravel()
+        array = np.genfromtxt(astr, dtype=dtype).ravel()
         return array
diff --git a/flopy4/block.py b/flopy4/block.py
index e5722b2..7310768 100644
--- a/flopy4/block.py
+++ b/flopy4/block.py
@@ -225,6 +225,7 @@ def load(cls, f, **kwargs):
         name = None
         index = None
         found = False
+        period = False
         params = dict()
         members = cls.params
 
@@ -238,12 +239,17 @@ def load(cls, f, **kwargs):
             if line == "\n":
                 continue
             words = strip(line).lower().split()
-            key = words[0]
+            if period:
+                key = "stress_period_data"
+            else:
+                key = words[0]
             if key == "begin":
                 found = True
                 name = words[1]
                 if len(words) > 2 and str.isdigit(words[2]):
                     index = int(words[2])
+                if name == "period":
+                    period = True
             elif key == "end":
                 break
             elif found:
@@ -268,12 +274,14 @@ def load(cls, f, **kwargs):
                     # TODO: inject from model somehow?
                     # and remove special handling here
                     kwrgs["cwd"] = ""
+                    # kwrgs["type"] = param.type
                     kwrgs["mempath"] = f"{mempath}/{name}"
-                if ptype is not MFArray:
+                if ptype is not MFArray and ptype is not MFList:
                     kwrgs.pop("model_shape", None)
                     kwrgs.pop("blk_params", None)
 
                 params[param.name] = ptype.load(f, **kwrgs)
+                period = False
 
         return cls(name=name, index=index, params=params)
 
@@ -281,7 +289,7 @@ def write(self, f):
         """Write the block to file."""
         index = self.index if self.index is not None else ""
         begin = f"BEGIN {self.name.upper()} {index}\n"
-        end = f"END {self.name.upper()}\n"
+        end = f"END {self.name.upper()}\n\n"
 
         f.write(begin)
         super().write(f)
diff --git a/flopy4/compound.py b/flopy4/compound.py
index 3ea38c0..96e419c 100644
--- a/flopy4/compound.py
+++ b/flopy4/compound.py
@@ -7,7 +7,7 @@
 
 from flopy4.array import MFArray, MFArrayType
 from flopy4.param import MFParam, MFParams, MFReader
-from flopy4.scalar import MFDouble, MFInteger, MFScalar
+from flopy4.scalar import MFScalar
 from flopy4.utils import strip
 
 PAD = "  "
@@ -338,21 +338,39 @@ def load(cls, f, **kwargs) -> "MFList":
         """Load list input with the given component parameters from a file."""
 
         blk_params = kwargs.pop("blk_params", {})
+        model_shape = kwargs.pop("model_shape", None)
         params = kwargs.pop("params", None)
-        type = kwargs.pop("type", None)
         kwargs.pop("mname", None)
-        kwargs.pop("shape", None)
+        kwargs.pop("shape", None)  # e.g. maxbound
+
+        param_lists = []
+        param_cols = []
+        param_types = []
+        for k in list(params):
+            if params[k].name == "aux" or params[k].name == "boundname":
+                continue
+                # raise NotImplementedError(
+                #    "boundames and auxvars not yet supported in period blocks"
+                # )
+            pcols = 0
+            if params[k].shape is None or params[k].shape == "":
+                pcols = 1
+            elif params[k].shape == "(ncelldim)":
+                if model_shape:
+                    pcols = len(model_shape)
+                else:
+                    raise ValueError("model_shape not set")
+            else:
+                pcols = len(params[k].shape.split(","))
+            param_cols.append(pcols)
+            param_lists.append(list())
+            param_types.append(params[k].type)
 
         if list(params.items())[-1][1].shape == "(:)":
-            maxsplit = len(params) - 1
+            maxsplit = sum(param_cols) - 1
         else:
             maxsplit = -1
 
-        param_lists = []
-        # TODO: support multi-dimensional params
-        for i in range(len(params)):
-            param_lists.append(list())
-
         while True:
             pos = f.tell()
             line = f.readline()
@@ -361,9 +379,28 @@ def load(cls, f, **kwargs) -> "MFList":
                 break
             else:
                 tokens = strip(line).split(maxsplit=maxsplit)
-                assert len(tokens) == len(param_lists)
-                for i, token in enumerate(tokens):
-                    param_lists[i].append(token)
+                assert len(tokens) == sum(param_cols)
+                icol = 0
+                for i, p in enumerate(param_lists):
+                    if param_cols[i] == 1:
+                        if param_types[i] == "integer":
+                            param_lists[i].append(int(tokens[icol]))
+                        elif param_types[i] == "double":
+                            param_lists[i].append(float(tokens[icol]))
+                        else:
+                            param_lists[i].append(tokens[icol])
+                        icol += 1
+                    else:
+                        row_l = []
+                        for j in range(param_cols[i]):
+                            if param_types[i] == "integer":
+                                row_l.append(int(tokens[icol]))
+                            elif param_types[i] == "double":
+                                row_l.append(float(tokens[icol]))
+                            else:
+                                row_l.append(tokens[icol])
+                            icol += 1
+                        param_lists[i].append(row_l)
 
         if blk_params and "dimensions" in blk_params:
             nbound = blk_params.get("dimensions").get("nbound")
@@ -372,31 +409,41 @@ def load(cls, f, **kwargs) -> "MFList":
                     if len(param_list) > nbound:
                         raise ValueError("MFList nbound not satisfied")
 
-        list_params = MFList.create_list_params(params, param_lists, **kwargs)
-        return cls(list_params, type=type, **kwargs)
+        list_params = MFList.create_list_params(
+            params, param_lists, param_cols, **kwargs
+        )
+        return cls(list_params, **kwargs)
 
     @staticmethod
     def create_list_params(
         params: Dict[str, MFParam],
         param_lists: list,
+        param_cols: list,
         **kwargs,
     ) -> Dict[str, MFParam]:
         """Create the param dictionary"""
         idx = 0
         list_params = dict()
         for param_name, param in params.items():
-            if type(param) is MFDouble:
+            if param_name == "aux" or param_name == "boundname":
+                continue
+            shape = None
+            if param_cols[idx] == 1:
+                shape = len(param_lists[idx])
+            else:
+                shape = (len(param_lists[idx]), param_cols[idx])
+            if type(param) is MFArray and param.type == "double":
                 list_params[param_name] = MFArray(
-                    shape=len(param_lists[idx]),
+                    shape=shape,
                     array=np.array(param_lists[idx], dtype=np.float64),
                     how=MFArrayType.internal,
                     factor=1.0,
                     path=None,
                     **kwargs,
                 )
-            elif type(param) is MFInteger:
+            elif type(param) is MFArray and param.type == "integer":
                 list_params[param_name] = MFArray(
-                    shape=len(param_lists[idx]),
+                    shape=shape,
                     array=np.array(param_lists[idx], dtype=np.int32),
                     how=MFArrayType.internal,
                     factor=1,
@@ -406,7 +453,7 @@ def create_list_params(
             else:
                 list_params[param_name] = MFScalarList(
                     value=param_lists[idx],
-                    type=type(param),
+                    # type=type(param),
                     **kwargs,
                 )
 
@@ -427,5 +474,9 @@ def write(self, f, **kwargs):
         for i in range(count):
             line = f"{PAD}"
             for name, param in self.params.items():
-                line += f"{param.value[i]}\t"
+                if isinstance(param.value[i], np.ndarray):
+                    for v in param.value[i]:
+                        line += f"{v}\t"
+                else:
+                    line += f"{param.value[i]}\t"
             f.write(line + "\n")
diff --git a/flopy4/ispec/exg_gwfgwf.py b/flopy4/ispec/exg_gwfgwf.py
index 3444531..35cc9fe 100644
--- a/flopy4/ispec/exg_gwfgwf.py
+++ b/flopy4/ispec/exg_gwfgwf.py
@@ -463,7 +463,7 @@ class ExgGwfgwf(MFPackage):
     )
 
     aux = MFArray(
-        type = "array",
+        type = "double",
         block = "exchangedata",
         shape = "(naux)",
         reader = "urword",
diff --git a/flopy4/ispec/gwe_dis.py b/flopy4/ispec/gwe_dis.py
index 3530f28..712ccab 100644
--- a/flopy4/ispec/gwe_dis.py
+++ b/flopy4/ispec/gwe_dis.py
@@ -199,7 +199,7 @@ class GweDis(MFPackage):
     )
 
     delr = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol)",
         reader = "readarray",
@@ -211,7 +211,7 @@ class GweDis(MFPackage):
     )
 
     delc = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(nrow)",
         reader = "readarray",
@@ -223,7 +223,7 @@ class GweDis(MFPackage):
     )
 
     top = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow)",
         reader = "readarray",
@@ -235,7 +235,7 @@ class GweDis(MFPackage):
     )
 
     botm = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
@@ -247,7 +247,7 @@ class GweDis(MFPackage):
     )
 
     idomain = MFArray(
-        type = "array",
+        type = "integer",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
diff --git a/flopy4/ispec/gwf_chd.py b/flopy4/ispec/gwf_chd.py
new file mode 100644
index 0000000..1e0c838
--- /dev/null
+++ b/flopy4/ispec/gwf_chd.py
@@ -0,0 +1,270 @@
+# generated file
+from flopy4.array import MFArray
+from flopy4.compound import MFRecord, MFList
+from flopy4.package import MFPackage
+from flopy4.scalar import MFDouble, MFFilename, MFInteger, MFKeyword, MFString
+
+
+class GwfChd(MFPackage):
+    multipkg = False
+    stress = False
+    advanced = False
+
+    auxiliary = MFString(
+        type = "string",
+        block = "options",
+        shape = "(naux)",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to specify aux variables""",
+        description =
+"""REPLACE auxnames {'{#1}': 'Groundwater Flow'}""",
+    )
+
+    auxmultname = MFString(
+        type = "string",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""name of auxiliary variable for multiplier""",
+        description =
+"""REPLACE auxmultname {'{#1}': 'CHD head value'}""",
+    )
+
+    boundnames = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""""",
+        description =
+"""REPLACE boundnames {'{#1}': 'constant-head'}""",
+    )
+
+    print_input = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""print input to listing file""",
+        description =
+"""REPLACE print_input {'{#1}': 'constant-head'}""",
+    )
+
+    print_flows = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""print CHD flows to listing file""",
+        description =
+"""REPLACE print_flows {'{#1}': 'constant-head'}""",
+    )
+
+    save_flows = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""save CHD flows to budget file""",
+        description =
+"""REPLACE save_flows {'{#1}': 'constant-head'}""",
+    )
+
+    ts_filerecord = MFRecord(
+        type = "record",
+        params = {
+            "ts6": MFKeyword(),
+            "filein": MFKeyword(),
+            "ts6_filename": MFString(),
+        },
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""""",
+        description =
+"""""",
+    )
+
+    ts6 = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""head keyword""",
+        description =
+"""keyword to specify that record corresponds to a time-series file.""",
+    )
+
+    filein = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""file keyword""",
+        description =
+"""keyword to specify that an input filename is expected next.""",
+    )
+
+    ts6_filename = MFString(
+        type = "string",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""file name of time series information""",
+        description =
+"""REPLACE timeseriesfile {}""",
+    )
+
+    obs_filerecord = MFRecord(
+        type = "record",
+        params = {
+            "obs6": MFKeyword(),
+            "filein": MFKeyword(),
+            "obs6_filename": MFString(),
+        },
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""""",
+        description =
+"""""",
+    )
+
+    obs6 = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""obs keyword""",
+        description =
+"""keyword to specify that record corresponds to an observations file.""",
+    )
+
+    obs6_filename = MFString(
+        type = "string",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""obs6 input filename""",
+        description =
+"""REPLACE obs6_filename {'{#1}': 'constant-head'}""",
+    )
+
+    dev_no_newton = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""turn off Newton for unconfined cells""",
+        description =
+"""turn off Newton for unconfined cells""",
+    )
+
+    maxbound = MFInteger(
+        type = "integer",
+        block = "dimensions",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""maximum number of constant heads""",
+        description =
+"""REPLACE maxbound {'{#1}': 'constant-head'}""",
+    )
+
+    cellid = MFArray(
+        type = "integer",
+        block = "period",
+        shape = "(ncelldim)",
+        reader = "urword",
+        optional = False,
+        longname =
+"""cell identifier""",
+        description =
+"""REPLACE cellid {}""",
+    )
+
+    head = MFDouble(
+        type = "double",
+        block = "period",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""head value assigned to constant head""",
+        description =
+"""is the head at the boundary. If the Options block includes a
+TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values
+can be obtained from a time series by entering the time-series name in
+place of a numeric value.""",
+    )
+
+    aux = MFArray(
+        type = "double",
+        block = "period",
+        shape = "(naux)",
+        reader = "urword",
+        optional = True,
+        longname =
+"""auxiliary variables""",
+        description =
+"""REPLACE aux {'{#1}': 'constant head'}""",
+    )
+
+    boundname = MFString(
+        type = "string",
+        block = "period",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""constant head boundary name""",
+        description =
+"""REPLACE boundname {'{#1}': 'constant head boundary'}""",
+    )
+
+    stress_period_data = MFList(
+        type = "recarray",
+        params = {
+            "cellid": cellid,
+            "head": head,
+            "aux": aux,
+            "boundname": boundname,
+        },
+        block = "period",
+        shape = "(maxbound)",
+        reader = "urword",
+        optional = False,
+        longname =
+"""""",
+        description =
+"""""",
+    )
\ No newline at end of file
diff --git a/flopy4/ispec/gwf_dis.py b/flopy4/ispec/gwf_dis.py
index 2cd78a4..e21c00b 100644
--- a/flopy4/ispec/gwf_dis.py
+++ b/flopy4/ispec/gwf_dis.py
@@ -199,7 +199,7 @@ class GwfDis(MFPackage):
     )
 
     delr = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol)",
         reader = "readarray",
@@ -211,7 +211,7 @@ class GwfDis(MFPackage):
     )
 
     delc = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(nrow)",
         reader = "readarray",
@@ -223,7 +223,7 @@ class GwfDis(MFPackage):
     )
 
     top = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow)",
         reader = "readarray",
@@ -235,7 +235,7 @@ class GwfDis(MFPackage):
     )
 
     botm = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
@@ -247,7 +247,7 @@ class GwfDis(MFPackage):
     )
 
     idomain = MFArray(
-        type = "array",
+        type = "integer",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
diff --git a/flopy4/ispec/gwf_ic.py b/flopy4/ispec/gwf_ic.py
index 4b276fa..f1cf8d4 100644
--- a/flopy4/ispec/gwf_ic.py
+++ b/flopy4/ispec/gwf_ic.py
@@ -37,7 +37,7 @@ class GwfIc(MFPackage):
     )
 
     strt = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(nodes)",
         reader = "readarray",
diff --git a/flopy4/ispec/gwf_model.py b/flopy4/ispec/gwf_model.py
index b46f31c..2f13ab4 100644
--- a/flopy4/ispec/gwf_model.py
+++ b/flopy4/ispec/gwf_model.py
@@ -1,12 +1,16 @@
 # generated file
 from flopy4.model import MFModel
 from flopy4.resolver import Resolve
+from flopy4.ispec.gwf_chd import GwfChd
 from flopy4.ispec.gwf_dis import GwfDis
 from flopy4.ispec.gwf_ic import GwfIc
 from flopy4.ispec.gwf_nam import GwfNam
+from flopy4.ispec.gwf_npf import GwfNpf
 
 
 class GwfModel(MFModel, Resolve):
+    chd6 = GwfChd()
     dis6 = GwfDis()
     ic6 = GwfIc()
     nam6 = GwfNam()
+    npf6 = GwfNpf()
diff --git a/flopy4/ispec/gwf_npf.py b/flopy4/ispec/gwf_npf.py
new file mode 100644
index 0000000..818d523
--- /dev/null
+++ b/flopy4/ispec/gwf_npf.py
@@ -0,0 +1,615 @@
+# generated file
+from flopy4.array import MFArray
+from flopy4.compound import MFRecord, MFList
+from flopy4.package import MFPackage
+from flopy4.scalar import MFDouble, MFFilename, MFInteger, MFKeyword, MFString
+
+
+class GwfNpf(MFPackage):
+    multipkg = False
+    stress = False
+    advanced = False
+
+    save_flows = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to save NPF flows""",
+        description =
+"""keyword to indicate that budget flow terms will be written to the file
+specified with ``BUDGET SAVE FILE'' in Output Control.""",
+    )
+
+    print_flows = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to print NPF flows to listing file""",
+        description =
+"""keyword to indicate that calculated flows between cells will be
+printed to the listing file for every stress period time step in which
+``BUDGET PRINT'' is specified in Output Control. If there is no Output
+Control option and ``PRINT_FLOWS'' is specified, then flow rates are
+printed for the last time step of each stress period.  This option can
+produce extremely large list files because all cell-by-cell flows are
+printed.  It should only be used with the NPF Package for models that
+have a small number of cells.""",
+    )
+
+    alternative_cell_averaging = MFString(
+        type = "string",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""conductance weighting option""",
+        description =
+"""is a text keyword to indicate that an alternative method will be used
+for calculating the conductance for horizontal cell connections.  The
+text value for ALTERNATIVE_CELL_AVERAGING can be ``LOGARITHMIC'',
+``AMT-LMK'', or ``AMT-HMK''.  ``AMT-LMK'' signifies that the
+conductance will be calculated using arithmetic-mean thickness and
+logarithmic-mean hydraulic conductivity.  ``AMT-HMK'' signifies that
+the conductance will be calculated using arithmetic-mean thickness and
+harmonic-mean hydraulic conductivity. If the user does not specify a
+value for ALTERNATIVE_CELL_AVERAGING, then the harmonic-mean method
+will be used.  This option cannot be used if the XT3D option is
+invoked.""",
+    )
+
+    thickstrt = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to activate THICKSTRT option""",
+        description =
+"""indicates that cells having a negative ICELLTYPE are confined, and
+their cell thickness for conductance calculations will be computed as
+STRT-BOT rather than TOP-BOT.  This option should be used with caution
+as it only affects conductance calculations in the NPF Package.""",
+    )
+
+    cvoptions = MFRecord(
+        type = "record",
+        params = {
+            "variablecv": MFKeyword(),
+            "dewatered": MFKeyword(),
+        },
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""vertical conductance options""",
+        description =
+"""none""",
+    )
+
+    variablecv = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""keyword to activate VARIABLECV option""",
+        description =
+"""keyword to indicate that the vertical conductance will be calculated
+using the saturated thickness and properties of the overlying cell and
+the thickness and properties of the underlying cell.  If the DEWATERED
+keyword is also specified, then the vertical conductance is calculated
+using only the saturated thickness and properties of the overlying
+cell if the head in the underlying cell is below its top.  If these
+keywords are not specified, then the default condition is to calculate
+the vertical conductance at the start of the simulation using the
+initial head and the cell properties.  The vertical conductance
+remains constant for the entire simulation.""",
+    )
+
+    dewatered = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to activate DEWATERED option""",
+        description =
+"""If the DEWATERED keyword is specified, then the vertical conductance
+is calculated using only the saturated thickness and properties of the
+overlying cell if the head in the underlying cell is below its top.""",
+    )
+
+    perched = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to activate PERCHED option""",
+        description =
+"""keyword to indicate that when a cell is overlying a dewatered
+convertible cell, the head difference used in Darcy's Law is equal to
+the head in the overlying cell minus the bottom elevation of the
+overlying cell.  If not specified, then the default is to use the head
+difference between the two cells.""",
+    )
+
+    rewet_record = MFRecord(
+        type = "record",
+        params = {
+            "rewet": MFKeyword(),
+            "wetfct": MFDouble(),
+            "iwetit": MFInteger(),
+            "ihdwet": MFInteger(),
+        },
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""""",
+        description =
+"""""",
+    )
+
+    rewet = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""keyword to activate rewetting""",
+        description =
+"""activates model rewetting.  Rewetting is off by default.""",
+    )
+
+    wetfct = MFDouble(
+        type = "double",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""wetting factor to use for rewetting""",
+        description =
+"""is a keyword and factor that is included in the calculation of the
+head that is initially established at a cell when that cell is
+converted from dry to wet.""",
+    )
+
+    iwetit = MFInteger(
+        type = "integer",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""interval to use for rewetting""",
+        description =
+"""is a keyword and iteration interval for attempting to wet cells.
+Wetting is attempted every IWETIT iteration. This applies to outer
+iterations and not inner iterations. If IWETIT is specified as zero or
+less, then the value is changed to 1.""",
+    )
+
+    ihdwet = MFInteger(
+        type = "integer",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""flag to determine wetting equation""",
+        description =
+"""is a keyword and integer flag that determines which equation is used
+to define the initial head at cells that become wet.  If IHDWET is 0,
+h = BOT + WETFCT (hm - BOT). If IHDWET is not 0, h = BOT + WETFCT
+(THRESH).""",
+    )
+
+    xt3doptions = MFRecord(
+        type = "record",
+        params = {
+            "xt3d": MFKeyword(),
+            "rhs": MFKeyword(),
+        },
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to activate XT3D""",
+        description =
+"""none""",
+    )
+
+    xt3d = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""keyword to activate XT3D""",
+        description =
+"""keyword indicating that the XT3D formulation will be used.  If the RHS
+keyword is also included, then the XT3D additional terms will be added
+to the right-hand side.  If the RHS keyword is excluded, then the XT3D
+terms will be put into the coefficient matrix.  Use of XT3D will
+substantially increase the computational effort, but will result in
+improved accuracy for anisotropic conductivity fields and for
+unstructured grids in which the CVFD requirement is violated.  XT3D
+requires additional information about the shapes of grid cells.  If
+XT3D is active and the DISU Package is used, then the user will need
+to provide in the DISU Package the angldegx array in the
+CONNECTIONDATA block and the VERTICES and CELL2D blocks.""",
+    )
+
+    rhs = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to XT3D on right hand side""",
+        description =
+"""If the RHS keyword is also included, then the XT3D additional terms
+will be added to the right-hand side.  If the RHS keyword is excluded,
+then the XT3D terms will be put into the coefficient matrix.""",
+    )
+
+    save_specific_discharge = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to save specific discharge""",
+        description =
+"""keyword to indicate that x, y, and z components of specific discharge
+will be calculated at cell centers and written to the budget file,
+which is specified with ``BUDGET SAVE FILE'' in Output Control.  If
+this option is activated, then additional information may be required
+in the discretization packages and the GWF Exchange package (if GWF
+models are coupled).  Specifically, ANGLDEGX must be specified in the
+CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be
+specified for the GWF Exchange as an auxiliary variable.""",
+    )
+
+    save_saturation = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to save saturation""",
+        description =
+"""keyword to indicate that cell saturation will be written to the budget
+file, which is specified with ``BUDGET SAVE FILE'' in Output Control.
+Saturation will be saved to the budget file as an auxiliary variable
+saved with the DATA-SAT text label.  Saturation is a cell variable
+that ranges from zero to one and can be used by post processing
+programs to determine how much of a cell volume is saturated.  If
+ICELLTYPE is 0, then saturation is always one.""",
+    )
+
+    k22overk = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to indicate that specified K22 is a ratio""",
+        description =
+"""keyword to indicate that specified K22 is a ratio of K22 divided by K.
+If this option is specified, then the K22 array entered in the NPF
+Package will be multiplied by K after being read.""",
+    )
+
+    k33overk = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""keyword to indicate that specified K33 is a ratio""",
+        description =
+"""keyword to indicate that specified K33 is a ratio of K33 divided by K.
+If this option is specified, then the K33 array entered in the NPF
+Package will be multiplied by K after being read.""",
+    )
+
+    tvk_filerecord = MFRecord(
+        type = "record",
+        params = {
+            "tvk6": MFKeyword(),
+            "filein": MFKeyword(),
+            "tvk6_filename": MFString(),
+        },
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""""",
+        description =
+"""""",
+    )
+
+    tvk6 = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""tvk keyword""",
+        description =
+"""keyword to specify that record corresponds to a time-varying hydraulic
+conductivity (TVK) file.  The behavior of TVK and a description of the
+input file is provided separately.""",
+    )
+
+    filein = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""file keyword""",
+        description =
+"""keyword to specify that an input filename is expected next.""",
+    )
+
+    tvk6_filename = MFString(
+        type = "string",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""file name of TVK information""",
+        description =
+"""defines a time-varying hydraulic conductivity (TVK) input file.
+Records in the TVK file can be used to change hydraulic conductivity
+properties at specified times or stress periods.""",
+    )
+
+    export_array_ascii = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""export array variables to layered ascii files.""",
+        description =
+"""keyword that specifies input griddata arrays should be written to
+layered ascii output files.""",
+    )
+
+    export_array_netcdf = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""export array variables to netcdf output files.""",
+        description =
+"""keyword that specifies input griddata arrays should be written to the
+model output netcdf file.""",
+    )
+
+    dev_no_newton = MFKeyword(
+        type = "keyword",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""turn off Newton for unconfined cells""",
+        description =
+"""turn off Newton for unconfined cells""",
+    )
+
+    dev_omega = MFDouble(
+        type = "double",
+        block = "options",
+        shape = "",
+        reader = "urword",
+        optional = True,
+        longname =
+"""set saturation omega value""",
+        description =
+"""set saturation omega value""",
+    )
+
+    icelltype = MFArray(
+        type = "integer",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = False,
+        longname =
+"""confined or convertible indicator""",
+        description =
+"""flag for each cell that specifies how saturated thickness is treated.
+0 means saturated thickness is held constant;  $>$0 means saturated
+thickness varies with computed head when head is below the cell top;
+$<$0 means saturated thickness varies with computed head unless the
+THICKSTRT option is in effect.  When THICKSTRT is in effect, a
+negative value for ICELLTYPE indicates that the saturated thickness
+value used in conductance calculations in the NPF Package will be
+computed as STRT-BOT and held constant.  If the THICKSTRT option is
+not in effect, then negative values provided by the user for ICELLTYPE
+are automatically reassigned by the program to a value of one.""",
+    )
+
+    k = MFArray(
+        type = "double",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = False,
+        longname =
+"""hydraulic conductivity (L/T)""",
+        description =
+"""is the hydraulic conductivity.  For the common case in which the user
+would like to specify the horizontal hydraulic conductivity and the
+vertical hydraulic conductivity, then K should be assigned as the
+horizontal hydraulic conductivity, K33 should be assigned as the
+vertical hydraulic conductivity, and K22 and the three rotation angles
+should not be specified.  When more sophisticated anisotropy is
+required, then K corresponds to the K11 hydraulic conductivity axis.
+All included cells (IDOMAIN $>$ 0) must have a K value greater than
+zero.""",
+    )
+
+    k22 = MFArray(
+        type = "double",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = True,
+        longname =
+"""hydraulic conductivity of second ellipsoid axis""",
+        description =
+"""is the hydraulic conductivity of the second ellipsoid axis (or the
+ratio of K22/K if the K22OVERK option is specified); for an unrotated
+case this is the hydraulic conductivity in the y direction.  If K22 is
+not included in the GRIDDATA block, then K22 is set equal to K.  For a
+regular MODFLOW grid (DIS Package is used) in which no rotation angles
+are specified, K22 is the hydraulic conductivity along columns in the
+y direction. For an unstructured DISU grid, the user must assign
+principal x and y axes and provide the angle for each cell face
+relative to the assigned x direction.  All included cells (IDOMAIN $>$
+0) must have a K22 value greater than zero.""",
+    )
+
+    k33 = MFArray(
+        type = "double",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = True,
+        longname =
+"""hydraulic conductivity of third ellipsoid axis (L/T)""",
+        description =
+"""is the hydraulic conductivity of the third ellipsoid axis (or the
+ratio of K33/K if the K33OVERK option is specified); for an unrotated
+case, this is the vertical hydraulic conductivity.  When anisotropy is
+applied, K33 corresponds to the K33 tensor component.  All included
+cells (IDOMAIN $>$ 0) must have a K33 value greater than zero.""",
+    )
+
+    angle1 = MFArray(
+        type = "double",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = True,
+        longname =
+"""first anisotropy rotation angle (degrees)""",
+        description =
+"""is a rotation angle of the hydraulic conductivity tensor in degrees.
+The angle represents the first of three sequential rotations of the
+hydraulic conductivity ellipsoid. With the K11, K22, and K33 axes of
+the ellipsoid initially aligned with the x, y, and z coordinate axes,
+respectively, ANGLE1 rotates the ellipsoid about its K33 axis (within
+the x - y plane). A positive value represents counter-clockwise
+rotation when viewed from any point on the positive K33 axis, looking
+toward the center of the ellipsoid. A value of zero indicates that the
+K11 axis lies within the x - z plane. If ANGLE1 is not specified,
+default values of zero are assigned to ANGLE1, ANGLE2, and ANGLE3, in
+which case the K11, K22, and K33 axes are aligned with the x, y, and z
+axes, respectively.""",
+    )
+
+    angle2 = MFArray(
+        type = "double",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = True,
+        longname =
+"""second anisotropy rotation angle (degrees)""",
+        description =
+"""is a rotation angle of the hydraulic conductivity tensor in degrees.
+The angle represents the second of three sequential rotations of the
+hydraulic conductivity ellipsoid. Following the rotation by ANGLE1
+described above, ANGLE2 rotates the ellipsoid about its K22 axis (out
+of the x - y plane). An array can be specified for ANGLE2 only if
+ANGLE1 is also specified. A positive value of ANGLE2 represents
+clockwise rotation when viewed from any point on the positive K22
+axis, looking toward the center of the ellipsoid. A value of zero
+indicates that the K11 axis lies within the x - y plane. If ANGLE2 is
+not specified, default values of zero are assigned to ANGLE2 and
+ANGLE3; connections that are not user-designated as vertical are
+assumed to be strictly horizontal (that is, to have no z component to
+their orientation); and connection lengths are based on horizontal
+distances.""",
+    )
+
+    angle3 = MFArray(
+        type = "double",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = True,
+        longname =
+"""third anisotropy rotation angle (degrees)""",
+        description =
+"""is a rotation angle of the hydraulic conductivity tensor in degrees.
+The angle represents the third of three sequential rotations of the
+hydraulic conductivity ellipsoid. Following the rotations by ANGLE1
+and ANGLE2 described above, ANGLE3 rotates the ellipsoid about its K11
+axis. An array can be specified for ANGLE3 only if ANGLE1 and ANGLE2
+are also specified. An array must be specified for ANGLE3 if ANGLE2 is
+specified. A positive value of ANGLE3 represents clockwise rotation
+when viewed from any point on the positive K11 axis, looking toward
+the center of the ellipsoid. A value of zero indicates that the K22
+axis lies within the x - y plane.""",
+    )
+
+    wetdry = MFArray(
+        type = "double",
+        block = "griddata",
+        shape = "(nodes)",
+        reader = "readarray",
+        optional = True,
+        longname =
+"""wetdry threshold and factor""",
+        description =
+"""is a combination of the wetting threshold and a flag to indicate which
+neighboring cells can cause a cell to become wet. If WETDRY $<$ 0,
+only a cell below a dry cell can cause the cell to become wet. If
+WETDRY $>$ 0, the cell below a dry cell and horizontally adjacent
+cells can cause a cell to become wet. If WETDRY is 0, the cell cannot
+be wetted. The absolute value of WETDRY is the wetting threshold. When
+the sum of BOT and the absolute value of WETDRY at a dry cell is
+equaled or exceeded by the head at an adjacent cell, the cell is
+wetted.  WETDRY must be specified if ``REWET'' is specified in the
+OPTIONS block.  If ``REWET'' is not specified in the options block,
+then WETDRY can be entered, and memory will be allocated for it, even
+though it is not used.""",
+    )
\ No newline at end of file
diff --git a/flopy4/ispec/gwt_dis.py b/flopy4/ispec/gwt_dis.py
index b924243..80221fd 100644
--- a/flopy4/ispec/gwt_dis.py
+++ b/flopy4/ispec/gwt_dis.py
@@ -199,7 +199,7 @@ class GwtDis(MFPackage):
     )
 
     delr = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol)",
         reader = "readarray",
@@ -211,7 +211,7 @@ class GwtDis(MFPackage):
     )
 
     delc = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(nrow)",
         reader = "readarray",
@@ -223,7 +223,7 @@ class GwtDis(MFPackage):
     )
 
     top = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow)",
         reader = "readarray",
@@ -235,7 +235,7 @@ class GwtDis(MFPackage):
     )
 
     botm = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
@@ -247,7 +247,7 @@ class GwtDis(MFPackage):
     )
 
     idomain = MFArray(
-        type = "array",
+        type = "integer",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
diff --git a/flopy4/ispec/prt_dis.py b/flopy4/ispec/prt_dis.py
index c80c396..448a4e1 100644
--- a/flopy4/ispec/prt_dis.py
+++ b/flopy4/ispec/prt_dis.py
@@ -199,7 +199,7 @@ class PrtDis(MFPackage):
     )
 
     delr = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol)",
         reader = "readarray",
@@ -211,7 +211,7 @@ class PrtDis(MFPackage):
     )
 
     delc = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(nrow)",
         reader = "readarray",
@@ -223,7 +223,7 @@ class PrtDis(MFPackage):
     )
 
     top = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow)",
         reader = "readarray",
@@ -235,7 +235,7 @@ class PrtDis(MFPackage):
     )
 
     botm = MFArray(
-        type = "array",
+        type = "double",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
@@ -247,7 +247,7 @@ class PrtDis(MFPackage):
     )
 
     idomain = MFArray(
-        type = "array",
+        type = "integer",
         block = "griddata",
         shape = "(ncol, nrow, nlay)",
         reader = "readarray",
diff --git a/flopy4/ispec/prt_prp.py b/flopy4/ispec/prt_prp.py
index 6e38656..00f2486 100644
--- a/flopy4/ispec/prt_prp.py
+++ b/flopy4/ispec/prt_prp.py
@@ -283,111 +283,78 @@ class PrtPrp(MFPackage):
 cell is inactive at release time.""",
     )
 
-    release_timesrecord = MFRecord(
-        type = "record",
-        params = {
-            "release_times": MFKeyword(),
-            "times": MFArray(shape="(unknown)"),
-        },
-        block = "options",
-        shape = "",
-        reader = "urword",
-        optional = True,
-        longname =
-"""""",
-        description =
-"""""",
-    )
-
-    release_times = MFKeyword(
+    dev_forceternary = MFKeyword(
         type = "keyword",
         block = "options",
         shape = "",
         reader = "urword",
         optional = False,
         longname =
-"""""",
-        description =
-"""keyword indicating release times will follow""",
-    )
-
-    times = MFArray(
-        type = "array",
-        block = "options",
-        shape = "(unknown)",
-        reader = "urword",
-        optional = False,
-        longname =
-"""release times""",
+"""force ternary tracking method""",
         description =
-"""times to release, relative to the beginning of the simulation.
-RELEASE_TIMES and RELEASE_TIMESFILE are mutually exclusive.""",
+"""force use of the ternary tracking method regardless of cell type in
+DISV grids.""",
     )
 
-    release_timesfilerecord = MFRecord(
-        type = "record",
-        params = {
-            "release_timesfile": MFKeyword(),
-            "timesfile": MFString(),
-        },
+    release_time_tolerance = MFDouble(
+        type = "double",
         block = "options",
         shape = "",
         reader = "urword",
         optional = True,
         longname =
-"""""",
+"""release time coincidence tolerance""",
         description =
-"""""",
+"""real number indicating the tolerance within which to consider adjacent
+release times coincident. Coincident release times will be merged into
+a single release time. The default is $epsilon times 10^11$, where
+$epsilon$ is machine precision.""",
     )
 
-    release_timesfile = MFKeyword(
-        type = "keyword",
-        block = "options",
-        shape = "",
-        reader = "urword",
-        optional = False,
-        longname =
-"""""",
-        description =
-"""keyword indicating release times file name will follow""",
-    )
-
-    timesfile = MFString(
-        type = "string",
+    release_time_frequency = MFDouble(
+        type = "double",
         block = "options",
         shape = "",
         reader = "urword",
-        optional = False,
+        optional = True,
         longname =
-"""file keyword""",
+"""release time frequency""",
         description =
-"""name of the release times file.  RELEASE_TIMES and RELEASE_TIMESFILE
-are mutually exclusive.""",
+"""real number indicating the time frequency at which to release
+particles. This option can be used to schedule releases at a regular
+interval for the duration of the simulation, starting at the
+simulation start time. The release schedule is the union of this
+option, the RELEASETIMES block, and PERIOD block RELEASESETTING
+selections. If none of these are provided, a single release time is
+configured at the beginning of the first time step of the simulation's
+first stress period.""",
     )
 
-    dev_forceternary = MFKeyword(
-        type = "keyword",
-        block = "options",
+    nreleasepts = MFInteger(
+        type = "integer",
+        block = "dimensions",
         shape = "",
         reader = "urword",
         optional = False,
         longname =
-"""force ternary tracking method""",
+"""number of particle release points""",
         description =
-"""force use of the ternary tracking method regardless of cell type in
-DISV grids.""",
+"""is the number of particle release points.""",
     )
 
-    nreleasepts = MFInteger(
+    nreleasetimes = MFInteger(
         type = "integer",
         block = "dimensions",
         shape = "",
         reader = "urword",
         optional = False,
         longname =
-"""number of particle release points""",
+"""number of particle release times""",
         description =
-"""is the number of particle release points.""",
+"""is the number of particle release times specified in the RELEASETIMES
+block. This is not necessarily the total number of release times;
+release times are the union of RELEASE_TIME_FREQUENCY, RELEASETIMES
+block, and PERIOD block RELEASESETTING selections.""",
     )
 
     irptno = MFInteger(
@@ -407,7 +374,7 @@ class PrtPrp(MFPackage):
     )
 
     cellid = MFArray(
-        type = "array",
+        type = "integer",
         block = "packagedata",
         shape = "(ncelldim)",
         reader = "urword",
@@ -496,6 +463,34 @@ class PrtPrp(MFPackage):
 """""",
     )
 
+    time = MFDouble(
+        type = "double",
+        block = "releasetimes",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""release time""",
+        description =
+"""real value that defines the release time with respect to the
+simulation start time.""",
+    )
+
+    releasetimes = MFList(
+        type = "recarray",
+        params = {
+            "time": time,
+        },
+        block = "releasetimes",
+        shape = "(nreleasetimes)",
+        reader = "urword",
+        optional = False,
+        longname =
+"""""",
+        description =
+"""""",
+    )
+
     all = MFKeyword(
         type = "keyword",
         block = "period",
@@ -505,8 +500,8 @@ class PrtPrp(MFPackage):
         longname =
 """""",
         description =
-"""keyword to indicate release of particles at the start of all time
-steps in the period.""",
+"""keyword to indicate release at the start of all time steps in the
+period.""",
     )
 
     first = MFKeyword(
@@ -518,10 +513,23 @@ class PrtPrp(MFPackage):
         longname =
 """""",
         description =
-"""keyword to indicate release of particles at the start of the first
-time step in the period. This keyword may be used in conjunction with
-other keywords to release particles at the start of multiple time
-steps.""",
+"""keyword to indicate release at the start of the first time step in the
+period. This keyword may be used in conjunction with other
+RELEASESETTING options.""",
+    )
+
+    last = MFKeyword(
+        type = "keyword",
+        block = "period",
+        shape = "",
+        reader = "urword",
+        optional = False,
+        longname =
+"""""",
+        description =
+"""keyword to indicate release at the start of the last time step in the
+period. This keyword may be used in conjunction with other
+RELEASESETTING options.""",
     )
 
     frequency = MFInteger(
@@ -533,13 +541,12 @@ class PrtPrp(MFPackage):
         longname =
 """""",
         description =
-"""release particles at the specified time step frequency. This keyword
-may be used in conjunction with other keywords to release particles at
-the start of multiple time steps.""",
+"""release at the specified time step frequency. This keyword may be used
+in conjunction with other RELEASESETTING options.""",
     )
 
     steps = MFArray(
-        type = "array",
+        type = "integer",
         block = "period",
         shape = "(<nstp)",
         reader = "urword",
@@ -547,13 +554,12 @@ class PrtPrp(MFPackage):
         longname =
 """""",
         description =
-"""release particles at the start of each step specified in STEPS. This
-keyword may be used in conjunction with other keywords to release
-particles at the start of multiple time steps.""",
+"""release at the start of each step specified in STEPS. This option may
+be used in conjunction with other RELEASESETTING options.""",
     )
 
     fraction = MFArray(
-        type = "array",
+        type = "double",
         block = "period",
         shape = "(<nstp)",
         reader = "urword",
@@ -567,5 +573,7 @@ class PrtPrp(MFPackage):
 used with ALL, FIRST, or FREQUENCY. When used with STEPS, FRACTION may
 be a single value or an array of the same length as STEPS. If a single
 FRACTION value is provided with STEPS, the fraction applies to all
-steps.""",
+steps. NOTE: The FRACTION option has been removed. For fine control
+over release timing, specify times explicitly using the RELEASETIMES
+block.""",
     )
\ No newline at end of file
diff --git a/flopy4/model.py b/flopy4/model.py
index 4c81899..97f170a 100644
--- a/flopy4/model.py
+++ b/flopy4/model.py
@@ -95,11 +95,15 @@ def load(cls, f, **kwargs):
 
         packages["nam6"] = GwfNam.load(f, **kwargs)
 
-        MFModel.load_packages(members, packages["nam6"], packages, **kwargs)
+        simpath = Path(f.name).parent
+        MFModel.load_packages(
+            simpath, members, packages["nam6"], packages, **kwargs
+        )
         return cls(name=mname, mempath=mempath, mtype=mtype, packages=packages)
 
     @staticmethod
     def load_packages(
+        simpath: Path,
         members,
         blocks: Dict[str, MFBlock],
         packages: Dict[str, MFPackage],
@@ -124,7 +128,7 @@ def load_packages(
                 fname = param.value["fname"][i]
                 pname = param.value["pname"][i]
                 package = members.get(ftype.lower(), None)
-                with open(fname, "r") as f:
+                with open(Path(simpath / fname), "r") as f:
                     kwargs["model_shape"] = model_shape
                     kwargs["mempath"] = f"{mempath}/{pname}"
                     kwargs["ftype"] = ftype.lower()
@@ -142,13 +146,17 @@ def load_packages(
                         nodes = packages[pname].params["nodes"]
                         model_shape = nodes
 
-    def write(self, basepath, **kwargs):
+    def write(self, simpath, **kwargs):
         """Write the model to files."""
-        path = Path(basepath)
-        for p in self._p:
-            if self._p[p].name.endswith(".nam"):
-                with open(path / self._p[p].name, "w") as f:
-                    self._p[p].write(f)
-            else:
-                with open(path / self._p[p].name, "w") as f:
-                    self._p[p].write(f)
+        fname_l = self._p["nam6"].params["packages"]["fname"]
+        # ftype_l = self._p["nam6"].params["packages"]["ftype"]
+        pname_l = self._p["nam6"].params["packages"]["pname"]
+
+        # write package files
+        for index, p in enumerate(pname_l):
+            with open(Path(simpath / fname_l[index]), "w") as f:
+                self._p[p].write(f)
+
+        # write name file
+        with open(Path(simpath / f"{self.name}.nam"), "w") as f:
+            self._p["nam6"].write(f)
diff --git a/flopy4/package.py b/flopy4/package.py
index a2c15c8..f81c2d8 100644
--- a/flopy4/package.py
+++ b/flopy4/package.py
@@ -243,6 +243,7 @@ def load(cls, f, **kwargs):
         pname = strip(mempath.split("/")[-1])
         ftype = kwargs.pop("ftype", None)
         ptype = ftype.replace("6", "")
+        kwargs.pop("mname", None)
         kwargs.pop("modeltype", None)
 
         while True:
diff --git a/flopy4/param.py b/flopy4/param.py
index 23bfba7..54d94e7 100644
--- a/flopy4/param.py
+++ b/flopy4/param.py
@@ -239,6 +239,7 @@ def write(self, f, **kwargs):
                 if len(self.params[param.name]):
                     param.write(f, **kwargs)
             elif param.type is None:
+                print(param)
                 raise TypeError(
                     f"Unknown specification type for param '{param.name}'"
                 )
diff --git a/flopy4/simulation.py b/flopy4/simulation.py
index 2157c96..980c7f8 100644
--- a/flopy4/simulation.py
+++ b/flopy4/simulation.py
@@ -47,16 +47,17 @@ def load(cls, f, **kwargs):
         solvers = dict()
         sim_name = "sim"
         mempath = sim_name
+        simpath = Path(f.name).parent
 
         kwargs["mempath"] = f"{mempath}"
         kwargs["ftype"] = "nam6"
 
         nam = SimNam.load(f, **kwargs)
 
-        tdis = MFSimulation.load_tdis(nam, **kwargs)
-        MFSimulation.load_models(nam, models, **kwargs)
-        MFSimulation.load_exchanges(nam, exchanges, **kwargs)
-        MFSimulation.load_solvers(nam, solvers, **kwargs)
+        tdis = MFSimulation.load_tdis(simpath, nam, **kwargs)
+        MFSimulation.load_models(simpath, nam, models, **kwargs)
+        MFSimulation.load_exchanges(simpath, nam, exchanges, **kwargs)
+        MFSimulation.load_solvers(simpath, nam, solvers, **kwargs)
 
         return cls(
             name=sim_name,
@@ -69,6 +70,7 @@ def load(cls, f, **kwargs):
 
     @staticmethod
     def load_tdis(
+        simpath: Path,
         blocks: Dict[str, MFBlock],
         mempath,
         **kwargs,
@@ -82,13 +84,14 @@ def load_tdis(
             tdis_fpth = param.value
             if tdis_fpth is None:
                 raise ValueError("Invalid mfsim.name TDIS6 specification")
-            with open(tdis_fpth, "r") as f:
+            with open(Path(simpath / tdis_fpth), "r") as f:
                 kwargs["mempath"] = f"{mempath}"
                 tdis = SimTdis.load(f, **kwargs)
         return tdis
 
     @staticmethod
     def load_models(
+        simpath: Path,
         blocks: Dict[str, MFBlock],
         models: Dict[str, MFModel],
         mempath,
@@ -116,13 +119,14 @@ def load_models(
                     model = PrtModel
                 else:
                     model = None
-                with open(mfname, "r") as f:
+                with open(Path(simpath / mfname), "r") as f:
                     kwargs["mtype"] = mtype.lower()
                     kwargs["mempath"] = f"{mempath}/{mname}"
                     models[mname] = model.load(f, **kwargs)
 
     @staticmethod
     def load_exchanges(
+        simpath: Path,
         blocks: Dict[str, MFBlock],
         exchanges: Dict[str, Dict],
         mempath,
@@ -157,6 +161,7 @@ def load_exchanges(
 
     @staticmethod
     def load_solvers(
+        simpath: Path,
         blocks: Dict[str, MFBlock],
         solvers: Dict[str, Dict],
         mempath,
@@ -181,18 +186,29 @@ def load_solvers(
                     sln = SlnIms
                 else:
                     sln = None
-                with open(slnfname, "r") as f:
+                with open(Path(simpath / slnfname), "r") as f:
                     slnname = slntype.replace("6", "")
                     slnname = f"{slnname}_{i}"
                     kwargs["mempath"] = f"{mempath}/{slnname}"
                     solvers[slnname] = sln.load(f, **kwargs)
 
-    def write(self, basepath, **kwargs):
+    def write(self, simpath, **kwargs):
         """Write the simulation to files."""
-        path = Path(basepath)
-        with open(path / "mfsim.nam", "w") as f:
-            self.nam.write(f, **kwargs)
-        with open(path / f"{self.name}.tdis", "w") as f:
+        path = Path(simpath)
+
+        fpath = Path(self.nam.params["tdis6"])
+        newpath = Path(path / fpath.name)
+        with open(newpath, "w") as f:
             self.tdis.write(f, **kwargs)
+
+        # slntypes = self.nam.params["solutiongroup"]["slntype"]
+        slnfnames = self.nam.params["solutiongroup"]["slnfname"]
+        for index, sln in enumerate(self.solvers):
+            with open(path / slnfnames[index], "w") as f:
+                self.solvers[sln].write(f, **kwargs)
+
         for model in self.models:
-            self.models[model].write(basepath, **kwargs)
+            self.models[model].write(simpath, **kwargs)
+
+        with open(path / "mfsim.nam", "w") as f:
+            self.nam.write(f, **kwargs)
diff --git a/spec/dfn/gwf-chd.dfn b/spec/dfn/gwf-chd.dfn
new file mode 100644
index 0000000..56e3ca3
--- /dev/null
+++ b/spec/dfn/gwf-chd.dfn
@@ -0,0 +1,222 @@
+# --------------------- gwf chd options ---------------------
+# flopy multi-package
+# package-type stress-package
+
+block options
+name auxiliary
+type string
+shape (naux)
+reader urword
+optional true
+longname keyword to specify aux variables
+description REPLACE auxnames {'{#1}': 'Groundwater Flow'}
+
+block options
+name auxmultname
+type string
+shape
+reader urword
+optional true
+longname name of auxiliary variable for multiplier
+description REPLACE auxmultname {'{#1}': 'CHD head value'}
+
+block options
+name boundnames
+type keyword
+shape
+reader urword
+optional true
+longname
+description REPLACE boundnames {'{#1}': 'constant-head'}
+
+block options
+name print_input
+type keyword
+reader urword
+optional true
+longname print input to listing file
+description REPLACE print_input {'{#1}': 'constant-head'}
+mf6internal iprpak
+
+block options
+name print_flows
+type keyword
+reader urword
+optional true
+longname print CHD flows to listing file
+description REPLACE print_flows {'{#1}': 'constant-head'}
+mf6internal iprflow
+
+block options
+name save_flows
+type keyword
+reader urword
+optional true
+longname save CHD flows to budget file
+description REPLACE save_flows {'{#1}': 'constant-head'}
+mf6internal ipakcb
+
+block options
+name ts_filerecord
+type record ts6 filein ts6_filename
+shape
+reader urword
+tagged true
+optional true
+longname
+description
+
+block options
+name ts6
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname head keyword
+description keyword to specify that record corresponds to a time-series file.
+
+block options
+name filein
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname file keyword
+description keyword to specify that an input filename is expected next.
+
+block options
+name ts6_filename
+type string
+preserve_case true
+in_record true
+reader urword
+optional false
+tagged false
+longname file name of time series information
+description REPLACE timeseriesfile {}
+
+block options
+name obs_filerecord
+type record obs6 filein obs6_filename
+shape
+reader urword
+tagged true
+optional true
+longname
+description
+
+block options
+name obs6
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname obs keyword
+description keyword to specify that record corresponds to an observations file.
+
+block options
+name obs6_filename
+type string
+preserve_case true
+in_record true
+tagged false
+reader urword
+optional false
+longname obs6 input filename
+description REPLACE obs6_filename {'{#1}': 'constant-head'}
+
+# dev options
+block options
+name dev_no_newton
+type keyword
+reader urword
+optional true
+longname turn off Newton for unconfined cells
+description turn off Newton for unconfined cells
+mf6internal inewton
+
+# --------------------- gwf chd dimensions ---------------------
+
+block dimensions
+name maxbound
+type integer
+reader urword
+optional false
+longname maximum number of constant heads
+description REPLACE maxbound {'{#1}': 'constant-head'}
+
+
+# --------------------- gwf chd period ---------------------
+
+block period
+name iper
+type integer
+block_variable True
+in_record true
+tagged false
+shape
+valid
+reader urword
+optional false
+longname stress period number
+description REPLACE iper {}
+
+block period
+name stress_period_data
+type recarray cellid head aux boundname
+shape (maxbound)
+reader urword
+longname
+description
+mf6internal spd
+
+block period
+name cellid
+type integer
+shape (ncelldim)
+tagged false
+in_record true
+reader urword
+longname cell identifier
+description REPLACE cellid {}
+
+block period
+name head
+type double precision
+shape
+tagged false
+in_record true
+reader urword
+time_series true
+longname head value assigned to constant head
+description is the head at the boundary. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value.
+
+block period
+name aux
+type double precision
+in_record true
+tagged false
+shape (naux)
+reader urword
+optional true
+time_series true
+longname auxiliary variables
+description REPLACE aux {'{#1}': 'constant head'}
+mf6internal auxvar
+
+block period
+name boundname
+type string
+shape
+tagged false
+in_record true
+reader urword
+optional true
+longname constant head boundary name
+description REPLACE boundname {'{#1}': 'constant head boundary'}
diff --git a/spec/dfn/gwf-npf.dfn b/spec/dfn/gwf-npf.dfn
new file mode 100644
index 0000000..3937ea2
--- /dev/null
+++ b/spec/dfn/gwf-npf.dfn
@@ -0,0 +1,355 @@
+# --------------------- gwf npf options ---------------------
+
+block options
+name save_flows
+type keyword
+reader urword
+optional true
+longname keyword to save NPF flows
+description keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control.
+mf6internal ipakcb
+
+block options
+name print_flows
+type keyword
+reader urword
+optional true
+longname keyword to print NPF flows to listing file
+description keyword to indicate that calculated flows between cells will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period.  This option can produce extremely large list files because all cell-by-cell flows are printed.  It should only be used with the NPF Package for models that have a small number of cells.
+mf6internal iprflow
+
+block options
+name alternative_cell_averaging
+type string
+valid logarithmic amt-lmk amt-hmk
+reader urword
+optional true
+longname conductance weighting option
+description is a text keyword to indicate that an alternative method will be used for calculating the conductance for horizontal cell connections.  The text value for ALTERNATIVE\_CELL\_AVERAGING can be ``LOGARITHMIC'', ``AMT-LMK'', or ``AMT-HMK''.  ``AMT-LMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and logarithmic-mean hydraulic conductivity.  ``AMT-HMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and harmonic-mean hydraulic conductivity. If the user does not specify a value for ALTERNATIVE\_CELL\_AVERAGING, then the harmonic-mean method will be used.  This option cannot be used if the XT3D option is invoked.
+mf6internal cellavg
+
+block options
+name thickstrt
+type keyword
+reader urword
+optional true
+longname keyword to activate THICKSTRT option
+description indicates that cells having a negative ICELLTYPE are confined, and their cell thickness for conductance calculations will be computed as STRT-BOT rather than TOP-BOT.  This option should be used with caution as it only affects conductance calculations in the NPF Package.
+mf6internal ithickstrt
+
+block options
+name cvoptions
+type record variablecv dewatered
+reader urword
+optional true
+longname vertical conductance options
+description none
+
+block options
+name variablecv
+in_record true
+type keyword
+reader urword
+longname keyword to activate VARIABLECV option
+description keyword to indicate that the vertical conductance will be calculated using the saturated thickness and properties of the overlying cell and the thickness and properties of the underlying cell.  If the DEWATERED keyword is also specified, then the vertical conductance is calculated using only the saturated thickness and properties of the overlying cell if the head in the underlying cell is below its top.  If these keywords are not specified, then the default condition is to calculate the vertical conductance at the start of the simulation using the initial head and the cell properties.  The vertical conductance remains constant for the entire simulation.
+mf6internal ivarcv
+
+block options
+name dewatered
+in_record true
+type keyword
+reader urword
+optional true
+longname keyword to activate DEWATERED option
+description If the DEWATERED keyword is specified, then the vertical conductance is calculated using only the saturated thickness and properties of the overlying cell if the head in the underlying cell is below its top.
+mf6internal idewatcv
+
+block options
+name perched
+type keyword
+reader urword
+optional true
+longname keyword to activate PERCHED option
+description keyword to indicate that when a cell is overlying a dewatered convertible cell, the head difference used in Darcy's Law is equal to the head in the overlying cell minus the bottom elevation of the overlying cell.  If not specified, then the default is to use the head difference between the two cells.
+mf6internal iperched
+
+block options
+name rewet_record
+type record rewet wetfct iwetit ihdwet
+reader urword
+optional true
+longname
+description
+
+block options
+name rewet
+type keyword
+in_record true
+reader urword
+optional false
+longname keyword to activate rewetting
+description activates model rewetting.  Rewetting is off by default.
+mf6internal irewet
+
+block options
+name wetfct
+type double precision
+in_record true
+reader urword
+optional false
+longname wetting factor to use for rewetting
+description is a keyword and factor that is included in the calculation of the head that is initially established at a cell when that cell is converted from dry to wet.
+
+block options
+name iwetit
+type integer
+in_record true
+reader urword
+optional false
+longname interval to use for rewetting
+description is a keyword and iteration interval for attempting to wet cells. Wetting is attempted every IWETIT iteration. This applies to outer iterations and not inner iterations. If IWETIT is specified as zero or less, then the value is changed to 1.
+
+block options
+name ihdwet
+type integer
+in_record true
+reader urword
+optional false
+longname flag to determine wetting equation
+description is a keyword and integer flag that determines which equation is used to define the initial head at cells that become wet.  If IHDWET is 0, h = BOT + WETFCT (hm - BOT). If IHDWET is not 0, h = BOT + WETFCT (THRESH).
+
+block options
+name xt3doptions
+type record xt3d rhs
+reader urword
+optional true
+longname keyword to activate XT3D
+description none
+
+block options
+name xt3d
+in_record true
+type keyword
+reader urword
+longname keyword to activate XT3D
+description keyword indicating that the XT3D formulation will be used.  If the RHS keyword is also included, then the XT3D additional terms will be added to the right-hand side.  If the RHS keyword is excluded, then the XT3D terms will be put into the coefficient matrix.  Use of XT3D will substantially increase the computational effort, but will result in improved accuracy for anisotropic conductivity fields and for unstructured grids in which the CVFD requirement is violated.  XT3D requires additional information about the shapes of grid cells.  If XT3D is active and the DISU Package is used, then the user will need to provide in the DISU Package the angldegx array in the CONNECTIONDATA block and the VERTICES and CELL2D blocks.
+mf6internal ixt3d
+
+block options
+name rhs
+in_record true
+type keyword
+reader urword
+optional true
+longname keyword to XT3D on right hand side
+description If the RHS keyword is also included, then the XT3D additional terms will be added to the right-hand side.  If the RHS keyword is excluded, then the XT3D terms will be put into the coefficient matrix.
+mf6internal ixt3drhs
+
+block options
+name save_specific_discharge
+type keyword
+reader urword
+optional true
+longname keyword to save specific discharge
+description keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable.
+mf6internal isavspdis
+
+block options
+name save_saturation
+type keyword
+reader urword
+optional true
+longname keyword to save saturation
+description keyword to indicate that cell saturation will be written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  Saturation will be saved to the budget file as an auxiliary variable saved with the DATA-SAT text label.  Saturation is a cell variable that ranges from zero to one and can be used by post processing programs to determine how much of a cell volume is saturated.  If ICELLTYPE is 0, then saturation is always one.
+mf6internal isavsat
+
+block options
+name k22overk
+type keyword
+reader urword
+optional true
+longname keyword to indicate that specified K22 is a ratio
+description keyword to indicate that specified K22 is a ratio of K22 divided by K.  If this option is specified, then the K22 array entered in the NPF Package will be multiplied by K after being read.
+mf6internal ik22overk
+
+block options
+name k33overk
+type keyword
+reader urword
+optional true
+longname keyword to indicate that specified K33 is a ratio
+description keyword to indicate that specified K33 is a ratio of K33 divided by K.  If this option is specified, then the K33 array entered in the NPF Package will be multiplied by K after being read.
+mf6internal ik33overk
+
+block options
+name tvk_filerecord
+type record tvk6 filein tvk6_filename
+shape
+reader urword
+tagged true
+optional true
+longname
+description
+
+block options
+name tvk6
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname tvk keyword
+description keyword to specify that record corresponds to a time-varying hydraulic conductivity (TVK) file.  The behavior of TVK and a description of the input file is provided separately.
+
+block options
+name filein
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname file keyword
+description keyword to specify that an input filename is expected next.
+
+block options
+name tvk6_filename
+type string
+preserve_case true
+in_record true
+reader urword
+optional false
+tagged false
+longname file name of TVK information
+description defines a time-varying hydraulic conductivity (TVK) input file.  Records in the TVK file can be used to change hydraulic conductivity properties at specified times or stress periods.
+
+block options
+name export_array_ascii
+type keyword
+reader urword
+optional true
+mf6internal export_ascii
+longname export array variables to layered ascii files.
+description keyword that specifies input griddata arrays should be written to layered ascii output files.
+
+block options
+name export_array_netcdf
+type keyword
+reader urword
+optional true
+mf6internal export_nc
+longname export array variables to netcdf output files.
+description keyword that specifies input griddata arrays should be written to the model output netcdf file.
+
+# dev options
+
+block options
+name dev_no_newton
+type keyword
+reader urword
+optional true
+longname turn off Newton for unconfined cells
+description turn off Newton for unconfined cells
+mf6internal inewton
+
+block options
+name dev_omega
+type double precision
+reader urword
+optional true
+longname set saturation omega value
+description set saturation omega value
+mf6internal satomega
+
+# --------------------- gwf npf griddata ---------------------
+
+block griddata
+name icelltype
+type integer
+shape (nodes)
+valid
+reader readarray
+layered true
+optional
+longname confined or convertible indicator
+description flag for each cell that specifies how saturated thickness is treated.  0 means saturated thickness is held constant;  $>$0 means saturated thickness varies with computed head when head is below the cell top; $<$0 means saturated thickness varies with computed head unless the THICKSTRT option is in effect.  When THICKSTRT is in effect, a negative value for ICELLTYPE indicates that the saturated thickness value used in conductance calculations in the NPF Package will be computed as STRT-BOT and held constant.  If the THICKSTRT option is not in effect, then negative values provided by the user for ICELLTYPE are automatically reassigned by the program to a value of one.
+default_value 0
+
+block griddata
+name k
+type double precision
+shape (nodes)
+valid
+reader readarray
+layered true
+optional
+longname hydraulic conductivity (L/T)
+description is the hydraulic conductivity.  For the common case in which the user would like to specify the horizontal hydraulic conductivity and the vertical hydraulic conductivity, then K should be assigned as the horizontal hydraulic conductivity, K33 should be assigned as the vertical hydraulic conductivity, and K22 and the three rotation angles should not be specified.  When more sophisticated anisotropy is required, then K corresponds to the K11 hydraulic conductivity axis.  All included cells (IDOMAIN $>$ 0) must have a K value greater than zero.
+default_value 1.0
+
+block griddata
+name k22
+type double precision
+shape (nodes)
+valid
+reader readarray
+layered true
+optional true
+longname hydraulic conductivity of second ellipsoid axis
+description is the hydraulic conductivity of the second ellipsoid axis (or the ratio of K22/K if the K22OVERK option is specified); for an unrotated case this is the hydraulic conductivity in the y direction.  If K22 is not included in the GRIDDATA block, then K22 is set equal to K.  For a regular MODFLOW grid (DIS Package is used) in which no rotation angles are specified, K22 is the hydraulic conductivity along columns in the y direction. For an unstructured DISU grid, the user must assign principal x and y axes and provide the angle for each cell face relative to the assigned x direction.  All included cells (IDOMAIN $>$ 0) must have a K22 value greater than zero.
+
+block griddata
+name k33
+type double precision
+shape (nodes)
+valid
+reader readarray
+layered true
+optional true
+longname hydraulic conductivity of third ellipsoid axis (L/T)
+description is the hydraulic conductivity of the third ellipsoid axis (or the ratio of K33/K if the K33OVERK option is specified); for an unrotated case, this is the vertical hydraulic conductivity.  When anisotropy is applied, K33 corresponds to the K33 tensor component.  All included cells (IDOMAIN $>$ 0) must have a K33 value greater than zero.
+
+block griddata
+name angle1
+type double precision
+shape (nodes)
+valid
+reader readarray
+layered true
+optional true
+longname first anisotropy rotation angle (degrees)
+description is a rotation angle of the hydraulic conductivity tensor in degrees. The angle represents the first of three sequential rotations of the hydraulic conductivity ellipsoid. With the K11, K22, and K33 axes of the ellipsoid initially aligned with the x, y, and z coordinate axes, respectively, ANGLE1 rotates the ellipsoid about its K33 axis (within the x - y plane). A positive value represents counter-clockwise rotation when viewed from any point on the positive K33 axis, looking toward the center of the ellipsoid. A value of zero indicates that the K11 axis lies within the x - z plane. If ANGLE1 is not specified, default values of zero are assigned to ANGLE1, ANGLE2, and ANGLE3, in which case the K11, K22, and K33 axes are aligned with the x, y, and z axes, respectively.
+
+block griddata
+name angle2
+type double precision
+shape (nodes)
+valid
+reader readarray
+layered true
+optional true
+longname second anisotropy rotation angle (degrees)
+description is a rotation angle of the hydraulic conductivity tensor in degrees. The angle represents the second of three sequential rotations of the hydraulic conductivity ellipsoid. Following the rotation by ANGLE1 described above, ANGLE2 rotates the ellipsoid about its K22 axis (out of the x - y plane). An array can be specified for ANGLE2 only if ANGLE1 is also specified. A positive value of ANGLE2 represents clockwise rotation when viewed from any point on the positive K22 axis, looking toward the center of the ellipsoid. A value of zero indicates that the K11 axis lies within the x - y plane. If ANGLE2 is not specified, default values of zero are assigned to ANGLE2 and ANGLE3; connections that are not user-designated as vertical are assumed to be strictly horizontal (that is, to have no z component to their orientation); and connection lengths are based on horizontal distances.
+
+block griddata
+name angle3
+type double precision
+shape (nodes)
+valid
+reader readarray
+layered true
+optional true
+longname third anisotropy rotation angle (degrees)
+description is a rotation angle of the hydraulic conductivity tensor in degrees. The angle represents the third of three sequential rotations of the hydraulic conductivity ellipsoid. Following the rotations by ANGLE1 and ANGLE2 described above, ANGLE3 rotates the ellipsoid about its K11 axis. An array can be specified for ANGLE3 only if ANGLE1 and ANGLE2 are also specified. An array must be specified for ANGLE3 if ANGLE2 is specified. A positive value of ANGLE3 represents clockwise rotation when viewed from any point on the positive K11 axis, looking toward the center of the ellipsoid. A value of zero indicates that the K22 axis lies within the x - y plane.
+
+block griddata
+name wetdry
+type double precision
+shape (nodes)
+valid
+reader readarray
+layered true
+optional true
+longname wetdry threshold and factor
+description is a combination of the wetting threshold and a flag to indicate which neighboring cells can cause a cell to become wet. If WETDRY $<$ 0, only a cell below a dry cell can cause the cell to become wet. If WETDRY $>$ 0, the cell below a dry cell and horizontally adjacent cells can cause a cell to become wet. If WETDRY is 0, the cell cannot be wetted. The absolute value of WETDRY is the wetting threshold. When the sum of BOT and the absolute value of WETDRY at a dry cell is equaled or exceeded by the head at an adjacent cell, the cell is wetted.  WETDRY must be specified if ``REWET'' is specified in the OPTIONS block.  If ``REWET'' is not specified in the options block, then WETDRY can be entered, and memory will be allocated for it, even though it is not used.
diff --git a/spec/dfn/prt-prp.dfn b/spec/dfn/prt-prp.dfn
index 6affe2c..bb5edd4 100644
--- a/spec/dfn/prt-prp.dfn
+++ b/spec/dfn/prt-prp.dfn
@@ -168,77 +168,29 @@ longname drape
 description is a text keyword to indicate that if a particle's release point is in a cell that happens to be inactive at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is inactive at release time.
 
 block options
-name release_timesrecord
-type record release_times times
-shape
-reader urword
-tagged true
-optional true
-longname
-description
-
-block options
-name release_times
+name dev_forceternary
 type keyword
 reader urword
-in_record true
-tagged true
-shape
-longname
-description keyword indicating release times will follow
+optional false
+longname force ternary tracking method
+description force use of the ternary tracking method regardless of cell type in DISV grids.
+mf6internal ifrctrn
 
 block options
-name times
+name release_time_tolerance
 type double precision
-shape (unknown)
-reader urword
-in_record true
-tagged false
-repeating true
-longname release times
-description times to release, relative to the beginning of the simulation.  RELEASE\_TIMES and RELEASE\_TIMESFILE are mutually exclusive.
-
-block options
-name release_timesfilerecord
-type record release_timesfile timesfile
-shape
 reader urword
-tagged true
 optional true
-longname
-description
-
-block options
-name release_timesfile
-type keyword
-reader urword
-in_record true
-tagged true
-shape
-longname
-description keyword indicating release times file name will follow
-
-block options
-name timesfile
-type string
-preserve_case true
-shape
-in_record true
-reader urword
-tagged false
-optional false
-longname file keyword
-description name of the release times file.  RELEASE\_TIMES and RELEASE\_TIMESFILE are mutually exclusive.
+longname release time coincidence tolerance
+description real number indicating the tolerance within which to consider adjacent release times coincident. Coincident release times will be merged into a single release time. The default is $\epsilon \times 10^11$, where $\epsilon$ is machine precision.
 
 block options
-name dev_forceternary
-type keyword
+name release_time_frequency
+type double precision
 reader urword
-optional false
-longname force ternary tracking method
-description force use of the ternary tracking method regardless of cell type in DISV grids.
-mf6internal ifrctrn
-
+optional true
+longname release time frequency
+description real number indicating the time frequency at which to release particles. This option can be used to schedule releases at a regular interval for the duration of the simulation, starting at the simulation start time. The release schedule is the union of this option, the RELEASETIMES block, and PERIOD block RELEASESETTING selections. If none of these are provided, a single release time is configured at the beginning of the first time step of the simulation's first stress period.
 
 # --------------------- prt prp dimensions ---------------------
 
@@ -250,6 +202,14 @@ optional false
 longname number of particle release points
 description is the number of particle release points.
 
+block dimensions
+name nreleasetimes
+type integer
+reader urword
+optional false
+longname number of particle release times
+description is the number of particle release times specified in the RELEASETIMES block. This is not necessarily the total number of release times; release times are the union of RELEASE\_TIME\_FREQUENCY, RELEASETIMES block, and PERIOD block RELEASESETTING selections.
+
 # --------------------- prt prp packagedata ---------------------
 
 block packagedata
@@ -322,6 +282,26 @@ optional true
 longname release point name
 description name of the particle release point. BOUNDNAME is an ASCII character variable that can contain as many as 40 characters. If BOUNDNAME contains spaces in it, then the entire name must be enclosed within single quotes.
 
+# --------------------- prt prp releasetimes ---------------
+
+block releasetimes
+name releasetimes
+type recarray time
+shape (nreleasetimes)
+reader urword
+longname
+description
+
+block releasetimes
+name time
+type double precision
+shape
+tagged false
+in_record true
+reader urword
+longname release time
+description real value that defines the release time with respect to the simulation start time.
+
 # --------------------- prt prp period ---------------------
 
 block period
@@ -353,7 +333,7 @@ tagged false
 in_record true
 reader urword
 longname
-description specifies when to release particles within the stress period.  Overrides package-level RELEASETIME option, which applies to all stress periods. By default, RELEASESETTING configures particles for release at the beginning of the specified time steps. For time-offset releases, provide a FRACTION value.
+description specifies time steps at which to release a particle. A particle is released at the beginning of each specified time step. For fine control over release timing, specify times explicitly using the RELEASETIMES block. If the beginning of a specified time step coincides with a release time specified in the RELEASETIMES block or configured via RELEASE\_TIME\_FREQUENCY, only one particle is released at that time. Coincidence is evaluated up to the tolerance specified in RELEASE\_TIME\_TOLERANCE, or $\epsilon \times 10^11$ by default, where $\epsilon$ is machine precision. If no release times are configured via this setting, the RELEASETIMES block, or the RELEASE\_TIME\_FREQUENCY option, a single release time is configured at the beginning of the first time step of the simulation's first stress period.
 
 block period
 name all
@@ -362,7 +342,7 @@ shape
 in_record true
 reader urword
 longname
-description keyword to indicate release of particles at the start of all time steps in the period.
+description keyword to indicate release at the start of all time steps in the period.
 
 block period
 name first
@@ -371,7 +351,16 @@ shape
 in_record true
 reader urword
 longname
-description keyword to indicate release of particles at the start of the first time step in the period. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps.
+description keyword to indicate release at the start of the first time step in the period. This keyword may be used in conjunction with other RELEASESETTING options.
+
+block period
+name last
+type keyword
+shape
+in_record true
+reader urword
+longname
+description keyword to indicate release at the start of the last time step in the period. This keyword may be used in conjunction with other RELEASESETTING options.
 
 block period
 name frequency
@@ -381,7 +370,7 @@ tagged true
 in_record true
 reader urword
 longname
-description release particles at the specified time step frequency. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps.
+description release at the specified time step frequency. This keyword may be used in conjunction with other RELEASESETTING options.
 
 block period
 name steps
@@ -391,7 +380,7 @@ tagged true
 in_record true
 reader urword
 longname
-description release particles at the start of each step specified in STEPS. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps.
+description release at the start of each step specified in STEPS. This option may be used in conjunction with other RELEASESETTING options.
 
 block period
 name fraction
@@ -401,5 +390,6 @@ tagged true
 in_record true
 reader urword
 optional true
+removed 6.5.1
 longname
-description release particles after the specified fraction of the time step has elapsed. If FRACTION is not set, particles are released at the start of the specified time step(s). FRACTION must be a single value when used with ALL, FIRST, or FREQUENCY. When used with STEPS, FRACTION may be a single value or an array of the same length as STEPS. If a single FRACTION value is provided with STEPS, the fraction applies to all steps.
\ No newline at end of file
+description release particles after the specified fraction of the time step has elapsed. If FRACTION is not set, particles are released at the start of the specified time step(s). FRACTION must be a single value when used with ALL, FIRST, or FREQUENCY. When used with STEPS, FRACTION may be a single value or an array of the same length as STEPS. If a single FRACTION value is provided with STEPS, the fraction applies to all steps. NOTE: The FRACTION option has been removed. For fine control over release timing, specify times explicitly using the RELEASETIMES block.
diff --git a/spec/mf6pkg.template b/spec/mf6pkg.template
index 2df2378..69eac35 100644
--- a/spec/mf6pkg.template
+++ b/spec/mf6pkg.template
@@ -29,7 +29,7 @@ class {{c}}{{s}}(MFPackage):
         {%- else %}
 
     {{pname}} = MFArray(
-        type = "array",
+        type = "double",
         {%- endif %}
       {%- elif 'integer' == params[pname].type %}
         {%- if '' == params[pname].shape %}
@@ -39,7 +39,7 @@ class {{c}}{{s}}(MFPackage):
         {%- else %}
 
     {{pname}} = MFArray(
-        type = "array",
+        type = "integer",
         {%- endif %}
       {%- else %}
         {%- set tokens = params[pname].type.split(' ') %}
diff --git a/spec/toml/gwf-chd.toml b/spec/toml/gwf-chd.toml
new file mode 100644
index 0000000..d0a3d31
--- /dev/null
+++ b/spec/toml/gwf-chd.toml
@@ -0,0 +1,347 @@
+component = "GWF"
+subcomponent = "CHD"
+blocknames = [ "options", "dimensions", "period",]
+multipkg = false
+stress = true
+advanced = false
+multi = true
+
+[block.options.auxiliary]
+type = "string"
+block_variable = false
+valid = []
+shape = "(naux)"
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to specify aux variables"
+description = "REPLACE auxnames {'{#1}': 'Groundwater Flow'}"
+deprecated = ""
+
+[block.options.auxmultname]
+type = "string"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "name of auxiliary variable for multiplier"
+description = "REPLACE auxmultname {'{#1}': 'CHD head value'}"
+deprecated = ""
+
+[block.options.boundnames]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = ""
+description = "REPLACE boundnames {'{#1}': 'constant-head'}"
+deprecated = ""
+
+[block.options.print_input]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "print input to listing file"
+description = "REPLACE print_input {'{#1}': 'constant-head'}"
+deprecated = ""
+
+[block.options.print_flows]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "print CHD flows to listing file"
+description = "REPLACE print_flows {'{#1}': 'constant-head'}"
+deprecated = ""
+
+[block.options.save_flows]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "save CHD flows to budget file"
+description = "REPLACE save_flows {'{#1}': 'constant-head'}"
+deprecated = ""
+
+[block.options.ts_filerecord]
+type = "record ts6 filein ts6_filename"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = ""
+description = ""
+deprecated = ""
+
+[block.options.ts6]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "head keyword"
+description = "keyword to specify that record corresponds to a time-series file."
+deprecated = ""
+
+[block.options.filein]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "file keyword"
+description = "keyword to specify that an input filename is expected next."
+deprecated = ""
+
+[block.options.ts6_filename]
+type = "string"
+block_variable = false
+valid = []
+shape = ""
+tagged = false
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = true
+numeric_index = false
+longname = "file name of time series information"
+description = "REPLACE timeseriesfile {}"
+deprecated = ""
+
+[block.options.obs_filerecord]
+type = "record obs6 filein obs6_filename"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = ""
+description = ""
+deprecated = ""
+
+[block.options.obs6]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "obs keyword"
+description = "keyword to specify that record corresponds to an observations file."
+deprecated = ""
+
+[block.options.obs6_filename]
+type = "string"
+block_variable = false
+valid = []
+shape = ""
+tagged = false
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = true
+numeric_index = false
+longname = "obs6 input filename"
+description = "REPLACE obs6_filename {'{#1}': 'constant-head'}"
+deprecated = ""
+
+[block.options.dev_no_newton]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "turn off Newton for unconfined cells"
+description = "turn off Newton for unconfined cells"
+deprecated = ""
+
+[block.dimensions.maxbound]
+type = "integer"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "maximum number of constant heads"
+description = "REPLACE maxbound {'{#1}': 'constant-head'}"
+deprecated = ""
+
+[block.period.cellid]
+type = "integer"
+block_variable = false
+valid = []
+shape = "(ncelldim)"
+tagged = false
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "cell identifier"
+description = "REPLACE cellid {}"
+deprecated = ""
+
+[block.period.head]
+type = "double"
+block_variable = false
+valid = []
+shape = ""
+tagged = false
+in_record = true
+layered = false
+time_series = true
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "head value assigned to constant head"
+description = "is the head at the boundary. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value."
+deprecated = ""
+
+[block.period.aux]
+type = "double"
+block_variable = false
+valid = []
+shape = "(naux)"
+tagged = false
+in_record = true
+layered = false
+time_series = true
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "auxiliary variables"
+description = "REPLACE aux {'{#1}': 'constant head'}"
+deprecated = ""
+
+[block.period.boundname]
+type = "string"
+block_variable = false
+valid = []
+shape = ""
+tagged = false
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "constant head boundary name"
+description = "REPLACE boundname {'{#1}': 'constant head boundary'}"
+deprecated = ""
+
+[block.period.stress_period_data]
+type = "recarray cellid head aux boundname"
+block_variable = false
+valid = []
+shape = "(maxbound)"
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = ""
+description = ""
+deprecated = ""
diff --git a/spec/toml/gwf-npf.toml b/spec/toml/gwf-npf.toml
new file mode 100644
index 0000000..8e97362
--- /dev/null
+++ b/spec/toml/gwf-npf.toml
@@ -0,0 +1,620 @@
+component = "GWF"
+subcomponent = "NPF"
+blocknames = [ "options", "griddata",]
+multipkg = false
+stress = false
+advanced = false
+
+[block.options.save_flows]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to save NPF flows"
+description = "keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control."
+deprecated = ""
+
+[block.options.print_flows]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to print NPF flows to listing file"
+description = "keyword to indicate that calculated flows between cells will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period.  This option can produce extremely large list files because all cell-by-cell flows are printed.  It should only be used with the NPF Package for models that have a small number of cells."
+deprecated = ""
+
+[block.options.alternative_cell_averaging]
+type = "string"
+block_variable = false
+valid = [ "logarithmic", "amt-lmk", "amt-hmk",]
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "conductance weighting option"
+description = "is a text keyword to indicate that an alternative method will be used for calculating the conductance for horizontal cell connections.  The text value for ALTERNATIVE_CELL_AVERAGING can be ``LOGARITHMIC'', ``AMT-LMK'', or ``AMT-HMK''.  ``AMT-LMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and logarithmic-mean hydraulic conductivity.  ``AMT-HMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and harmonic-mean hydraulic conductivity. If the user does not specify a value for ALTERNATIVE_CELL_AVERAGING, then the harmonic-mean method will be used.  This option cannot be used if the XT3D option is invoked."
+deprecated = ""
+
+[block.options.thickstrt]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to activate THICKSTRT option"
+description = "indicates that cells having a negative ICELLTYPE are confined, and their cell thickness for conductance calculations will be computed as STRT-BOT rather than TOP-BOT.  This option should be used with caution as it only affects conductance calculations in the NPF Package."
+deprecated = ""
+
+[block.options.cvoptions]
+type = "record variablecv dewatered"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "vertical conductance options"
+description = "none"
+deprecated = ""
+
+[block.options.variablecv]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "keyword to activate VARIABLECV option"
+description = "keyword to indicate that the vertical conductance will be calculated using the saturated thickness and properties of the overlying cell and the thickness and properties of the underlying cell.  If the DEWATERED keyword is also specified, then the vertical conductance is calculated using only the saturated thickness and properties of the overlying cell if the head in the underlying cell is below its top.  If these keywords are not specified, then the default condition is to calculate the vertical conductance at the start of the simulation using the initial head and the cell properties.  The vertical conductance remains constant for the entire simulation."
+deprecated = ""
+
+[block.options.dewatered]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to activate DEWATERED option"
+description = "If the DEWATERED keyword is specified, then the vertical conductance is calculated using only the saturated thickness and properties of the overlying cell if the head in the underlying cell is below its top."
+deprecated = ""
+
+[block.options.perched]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to activate PERCHED option"
+description = "keyword to indicate that when a cell is overlying a dewatered convertible cell, the head difference used in Darcy's Law is equal to the head in the overlying cell minus the bottom elevation of the overlying cell.  If not specified, then the default is to use the head difference between the two cells."
+deprecated = ""
+
+[block.options.rewet_record]
+type = "record rewet wetfct iwetit ihdwet"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = ""
+description = ""
+deprecated = ""
+
+[block.options.rewet]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "keyword to activate rewetting"
+description = "activates model rewetting.  Rewetting is off by default."
+deprecated = ""
+
+[block.options.wetfct]
+type = "double"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "wetting factor to use for rewetting"
+description = "is a keyword and factor that is included in the calculation of the head that is initially established at a cell when that cell is converted from dry to wet."
+deprecated = ""
+
+[block.options.iwetit]
+type = "integer"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "interval to use for rewetting"
+description = "is a keyword and iteration interval for attempting to wet cells. Wetting is attempted every IWETIT iteration. This applies to outer iterations and not inner iterations. If IWETIT is specified as zero or less, then the value is changed to 1."
+deprecated = ""
+
+[block.options.ihdwet]
+type = "integer"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "flag to determine wetting equation"
+description = "is a keyword and integer flag that determines which equation is used to define the initial head at cells that become wet.  If IHDWET is 0, h = BOT + WETFCT (hm - BOT). If IHDWET is not 0, h = BOT + WETFCT (THRESH)."
+deprecated = ""
+
+[block.options.xt3doptions]
+type = "record xt3d rhs"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to activate XT3D"
+description = "none"
+deprecated = ""
+
+[block.options.xt3d]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "keyword to activate XT3D"
+description = "keyword indicating that the XT3D formulation will be used.  If the RHS keyword is also included, then the XT3D additional terms will be added to the right-hand side.  If the RHS keyword is excluded, then the XT3D terms will be put into the coefficient matrix.  Use of XT3D will substantially increase the computational effort, but will result in improved accuracy for anisotropic conductivity fields and for unstructured grids in which the CVFD requirement is violated.  XT3D requires additional information about the shapes of grid cells.  If XT3D is active and the DISU Package is used, then the user will need to provide in the DISU Package the angldegx array in the CONNECTIONDATA block and the VERTICES and CELL2D blocks."
+deprecated = ""
+
+[block.options.rhs]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to XT3D on right hand side"
+description = "If the RHS keyword is also included, then the XT3D additional terms will be added to the right-hand side.  If the RHS keyword is excluded, then the XT3D terms will be put into the coefficient matrix."
+deprecated = ""
+
+[block.options.save_specific_discharge]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to save specific discharge"
+description = "keyword to indicate that x, y, and z components of specific discharge will be calculated at cell centers and written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  If this option is activated, then additional information may be required in the discretization packages and the GWF Exchange package (if GWF models are coupled).  Specifically, ANGLDEGX must be specified in the CONNECTIONDATA block of the DISU Package; ANGLDEGX must also be specified for the GWF Exchange as an auxiliary variable."
+deprecated = ""
+
+[block.options.save_saturation]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to save saturation"
+description = "keyword to indicate that cell saturation will be written to the budget file, which is specified with ``BUDGET SAVE FILE'' in Output Control.  Saturation will be saved to the budget file as an auxiliary variable saved with the DATA-SAT text label.  Saturation is a cell variable that ranges from zero to one and can be used by post processing programs to determine how much of a cell volume is saturated.  If ICELLTYPE is 0, then saturation is always one."
+deprecated = ""
+
+[block.options.k22overk]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to indicate that specified K22 is a ratio"
+description = "keyword to indicate that specified K22 is a ratio of K22 divided by K.  If this option is specified, then the K22 array entered in the NPF Package will be multiplied by K after being read."
+deprecated = ""
+
+[block.options.k33overk]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "keyword to indicate that specified K33 is a ratio"
+description = "keyword to indicate that specified K33 is a ratio of K33 divided by K.  If this option is specified, then the K33 array entered in the NPF Package will be multiplied by K after being read."
+deprecated = ""
+
+[block.options.tvk_filerecord]
+type = "record tvk6 filein tvk6_filename"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = ""
+description = ""
+deprecated = ""
+
+[block.options.tvk6]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "tvk keyword"
+description = "keyword to specify that record corresponds to a time-varying hydraulic conductivity (TVK) file.  The behavior of TVK and a description of the input file is provided separately."
+deprecated = ""
+
+[block.options.filein]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "file keyword"
+description = "keyword to specify that an input filename is expected next."
+deprecated = ""
+
+[block.options.tvk6_filename]
+type = "string"
+block_variable = false
+valid = []
+shape = ""
+tagged = false
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = true
+numeric_index = false
+longname = "file name of TVK information"
+description = "defines a time-varying hydraulic conductivity (TVK) input file.  Records in the TVK file can be used to change hydraulic conductivity properties at specified times or stress periods."
+deprecated = ""
+
+[block.options.export_array_ascii]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "export array variables to layered ascii files."
+description = "keyword that specifies input griddata arrays should be written to layered ascii output files."
+deprecated = ""
+
+[block.options.export_array_netcdf]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "export array variables to netcdf output files."
+description = "keyword that specifies input griddata arrays should be written to the model output netcdf file."
+deprecated = ""
+
+[block.options.dev_no_newton]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "turn off Newton for unconfined cells"
+description = "turn off Newton for unconfined cells"
+deprecated = ""
+
+[block.options.dev_omega]
+type = "double"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "set saturation omega value"
+description = "set saturation omega value"
+deprecated = ""
+
+[block.griddata.icelltype]
+type = "integer"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = false
+preserve_case = false
+default_value = "0"
+numeric_index = false
+longname = "confined or convertible indicator"
+description = "flag for each cell that specifies how saturated thickness is treated.  0 means saturated thickness is held constant;  $>$0 means saturated thickness varies with computed head when head is below the cell top; $<$0 means saturated thickness varies with computed head unless the THICKSTRT option is in effect.  When THICKSTRT is in effect, a negative value for ICELLTYPE indicates that the saturated thickness value used in conductance calculations in the NPF Package will be computed as STRT-BOT and held constant.  If the THICKSTRT option is not in effect, then negative values provided by the user for ICELLTYPE are automatically reassigned by the program to a value of one."
+deprecated = ""
+
+[block.griddata.k]
+type = "double"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = false
+preserve_case = false
+default_value = "1.0"
+numeric_index = false
+longname = "hydraulic conductivity (L/T)"
+description = "is the hydraulic conductivity.  For the common case in which the user would like to specify the horizontal hydraulic conductivity and the vertical hydraulic conductivity, then K should be assigned as the horizontal hydraulic conductivity, K33 should be assigned as the vertical hydraulic conductivity, and K22 and the three rotation angles should not be specified.  When more sophisticated anisotropy is required, then K corresponds to the K11 hydraulic conductivity axis.  All included cells (IDOMAIN $>$ 0) must have a K value greater than zero."
+deprecated = ""
+
+[block.griddata.k22]
+type = "double"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "hydraulic conductivity of second ellipsoid axis"
+description = "is the hydraulic conductivity of the second ellipsoid axis (or the ratio of K22/K if the K22OVERK option is specified); for an unrotated case this is the hydraulic conductivity in the y direction.  If K22 is not included in the GRIDDATA block, then K22 is set equal to K.  For a regular MODFLOW grid (DIS Package is used) in which no rotation angles are specified, K22 is the hydraulic conductivity along columns in the y direction. For an unstructured DISU grid, the user must assign principal x and y axes and provide the angle for each cell face relative to the assigned x direction.  All included cells (IDOMAIN $>$ 0) must have a K22 value greater than zero."
+deprecated = ""
+
+[block.griddata.k33]
+type = "double"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "hydraulic conductivity of third ellipsoid axis (L/T)"
+description = "is the hydraulic conductivity of the third ellipsoid axis (or the ratio of K33/K if the K33OVERK option is specified); for an unrotated case, this is the vertical hydraulic conductivity.  When anisotropy is applied, K33 corresponds to the K33 tensor component.  All included cells (IDOMAIN $>$ 0) must have a K33 value greater than zero."
+deprecated = ""
+
+[block.griddata.angle1]
+type = "double"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "first anisotropy rotation angle (degrees)"
+description = "is a rotation angle of the hydraulic conductivity tensor in degrees. The angle represents the first of three sequential rotations of the hydraulic conductivity ellipsoid. With the K11, K22, and K33 axes of the ellipsoid initially aligned with the x, y, and z coordinate axes, respectively, ANGLE1 rotates the ellipsoid about its K33 axis (within the x - y plane). A positive value represents counter-clockwise rotation when viewed from any point on the positive K33 axis, looking toward the center of the ellipsoid. A value of zero indicates that the K11 axis lies within the x - z plane. If ANGLE1 is not specified, default values of zero are assigned to ANGLE1, ANGLE2, and ANGLE3, in which case the K11, K22, and K33 axes are aligned with the x, y, and z axes, respectively."
+deprecated = ""
+
+[block.griddata.angle2]
+type = "double"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "second anisotropy rotation angle (degrees)"
+description = "is a rotation angle of the hydraulic conductivity tensor in degrees. The angle represents the second of three sequential rotations of the hydraulic conductivity ellipsoid. Following the rotation by ANGLE1 described above, ANGLE2 rotates the ellipsoid about its K22 axis (out of the x - y plane). An array can be specified for ANGLE2 only if ANGLE1 is also specified. A positive value of ANGLE2 represents clockwise rotation when viewed from any point on the positive K22 axis, looking toward the center of the ellipsoid. A value of zero indicates that the K11 axis lies within the x - y plane. If ANGLE2 is not specified, default values of zero are assigned to ANGLE2 and ANGLE3; connections that are not user-designated as vertical are assumed to be strictly horizontal (that is, to have no z component to their orientation); and connection lengths are based on horizontal distances."
+deprecated = ""
+
+[block.griddata.angle3]
+type = "double"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "third anisotropy rotation angle (degrees)"
+description = "is a rotation angle of the hydraulic conductivity tensor in degrees. The angle represents the third of three sequential rotations of the hydraulic conductivity ellipsoid. Following the rotations by ANGLE1 and ANGLE2 described above, ANGLE3 rotates the ellipsoid about its K11 axis. An array can be specified for ANGLE3 only if ANGLE1 and ANGLE2 are also specified. An array must be specified for ANGLE3 if ANGLE2 is specified. A positive value of ANGLE3 represents clockwise rotation when viewed from any point on the positive K11 axis, looking toward the center of the ellipsoid. A value of zero indicates that the K22 axis lies within the x - y plane."
+deprecated = ""
+
+[block.griddata.wetdry]
+type = "double"
+block_variable = false
+valid = []
+shape = "(nodes)"
+tagged = true
+in_record = false
+layered = true
+time_series = false
+reader = "readarray"
+optional = true
+preserve_case = false
+numeric_index = false
+longname = "wetdry threshold and factor"
+description = "is a combination of the wetting threshold and a flag to indicate which neighboring cells can cause a cell to become wet. If WETDRY $<$ 0, only a cell below a dry cell can cause the cell to become wet. If WETDRY $>$ 0, the cell below a dry cell and horizontally adjacent cells can cause a cell to become wet. If WETDRY is 0, the cell cannot be wetted. The absolute value of WETDRY is the wetting threshold. When the sum of BOT and the absolute value of WETDRY at a dry cell is equaled or exceeded by the head at an adjacent cell, the cell is wetted.  WETDRY must be specified if ``REWET'' is specified in the OPTIONS block.  If ``REWET'' is not specified in the options block, then WETDRY can be entered, and memory will be allocated for it, even though it is not used."
+deprecated = ""
diff --git a/spec/toml/prt-prp.toml b/spec/toml/prt-prp.toml
index 504b69b..f48011f 100644
--- a/spec/toml/prt-prp.toml
+++ b/spec/toml/prt-prp.toml
@@ -1,6 +1,6 @@
 component = "PRT"
 subcomponent = "PRP"
-blocknames = [ "options", "dimensions", "packagedata", "period",]
+blocknames = [ "options", "dimensions", "packagedata", "releasetimes", "period",]
 multipkg = false
 stress = false
 advanced = false
@@ -312,61 +312,27 @@ longname = "drape"
 description = "is a text keyword to indicate that if a particle's release point is in a cell that happens to be inactive at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is inactive at release time."
 deprecated = ""
 
-[block.options.release_timesrecord]
-type = "record release_times times"
-block_variable = false
-valid = []
-shape = ""
-tagged = true
-in_record = false
-layered = false
-time_series = false
-reader = "urword"
-optional = true
-preserve_case = false
-numeric_index = false
-longname = ""
-description = ""
-deprecated = ""
-
-[block.options.release_times]
+[block.options.dev_forceternary]
 type = "keyword"
 block_variable = false
 valid = []
 shape = ""
 tagged = true
-in_record = true
+in_record = false
 layered = false
 time_series = false
 reader = "urword"
 optional = false
 preserve_case = false
 numeric_index = false
-longname = ""
-description = "keyword indicating release times will follow"
+longname = "force ternary tracking method"
+description = "force use of the ternary tracking method regardless of cell type in DISV grids."
 deprecated = ""
 
-[block.options.times]
+[block.options.release_time_tolerance]
 type = "double"
 block_variable = false
 valid = []
-shape = "(unknown)"
-tagged = false
-in_record = true
-layered = false
-time_series = false
-reader = "urword"
-optional = false
-preserve_case = false
-numeric_index = false
-longname = "release times"
-description = "times to release, relative to the beginning of the simulation.  RELEASE_TIMES and RELEASE_TIMESFILE are mutually exclusive."
-deprecated = ""
-
-[block.options.release_timesfilerecord]
-type = "record release_timesfile timesfile"
-block_variable = false
-valid = []
 shape = ""
 tagged = true
 in_record = false
@@ -376,46 +342,29 @@ reader = "urword"
 optional = true
 preserve_case = false
 numeric_index = false
-longname = ""
-description = ""
+longname = "release time coincidence tolerance"
+description = "real number indicating the tolerance within which to consider adjacent release times coincident. Coincident release times will be merged into a single release time. The default is $epsilon times 10^11$, where $epsilon$ is machine precision."
 deprecated = ""
 
-[block.options.release_timesfile]
-type = "keyword"
+[block.options.release_time_frequency]
+type = "double"
 block_variable = false
 valid = []
 shape = ""
 tagged = true
-in_record = true
+in_record = false
 layered = false
 time_series = false
 reader = "urword"
-optional = false
+optional = true
 preserve_case = false
 numeric_index = false
-longname = ""
-description = "keyword indicating release times file name will follow"
+longname = "release time frequency"
+description = "real number indicating the time frequency at which to release particles. This option can be used to schedule releases at a regular interval for the duration of the simulation, starting at the simulation start time. The release schedule is the union of this option, the RELEASETIMES block, and PERIOD block RELEASESETTING selections. If none of these are provided, a single release time is configured at the beginning of the first time step of the simulation's first stress period."
 deprecated = ""
 
-[block.options.timesfile]
-type = "string"
-block_variable = false
-valid = []
-shape = ""
-tagged = false
-in_record = true
-layered = false
-time_series = false
-reader = "urword"
-optional = false
-preserve_case = true
-numeric_index = false
-longname = "file keyword"
-description = "name of the release times file.  RELEASE_TIMES and RELEASE_TIMESFILE are mutually exclusive."
-deprecated = ""
-
-[block.options.dev_forceternary]
-type = "keyword"
+[block.dimensions.nreleasepts]
+type = "integer"
 block_variable = false
 valid = []
 shape = ""
@@ -427,11 +376,11 @@ reader = "urword"
 optional = false
 preserve_case = false
 numeric_index = false
-longname = "force ternary tracking method"
-description = "force use of the ternary tracking method regardless of cell type in DISV grids."
+longname = "number of particle release points"
+description = "is the number of particle release points."
 deprecated = ""
 
-[block.dimensions.nreleasepts]
+[block.dimensions.nreleasetimes]
 type = "integer"
 block_variable = false
 valid = []
@@ -444,8 +393,8 @@ reader = "urword"
 optional = false
 preserve_case = false
 numeric_index = false
-longname = "number of particle release points"
-description = "is the number of particle release points."
+longname = "number of particle release times"
+description = "is the number of particle release times specified in the RELEASETIMES block. This is not necessarily the total number of release times; release times are the union of RELEASE_TIME_FREQUENCY, RELEASETIMES block, and PERIOD block RELEASESETTING selections."
 deprecated = ""
 
 [block.packagedata.irptno]
@@ -567,6 +516,40 @@ longname = ""
 description = ""
 deprecated = ""
 
+[block.releasetimes.time]
+type = "double"
+block_variable = false
+valid = []
+shape = ""
+tagged = false
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = "release time"
+description = "real value that defines the release time with respect to the simulation start time."
+deprecated = ""
+
+[block.releasetimes.releasetimes]
+type = "recarray time"
+block_variable = false
+valid = []
+shape = "(nreleasetimes)"
+tagged = true
+in_record = false
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = ""
+description = ""
+deprecated = ""
+
 [block.period.releasesetting]
 type = "keystring all first frequency steps fraction"
 block_variable = false
@@ -581,7 +564,7 @@ optional = false
 preserve_case = false
 numeric_index = false
 longname = ""
-description = "specifies when to release particles within the stress period.  Overrides package-level RELEASETIME option, which applies to all stress periods. By default, RELEASESETTING configures particles for release at the beginning of the specified time steps. For time-offset releases, provide a FRACTION value."
+description = "specifies time steps at which to release a particle. A particle is released at the beginning of each specified time step. For fine control over release timing, specify times explicitly using the RELEASETIMES block. If the beginning of a specified time step coincides with a release time specified in the RELEASETIMES block or configured via RELEASE_TIME_FREQUENCY, only one particle is released at that time. Coincidence is evaluated up to the tolerance specified in RELEASE_TIME_TOLERANCE, or $epsilon times 10^11$ by default, where $epsilon$ is machine precision. If no release times are configured via this setting, the RELEASETIMES block, or the RELEASE_TIME_FREQUENCY option, a single release time is configured at the beginning of the first time step of the simulation's first stress period."
 deprecated = ""
 
 [block.period.all]
@@ -598,7 +581,7 @@ optional = false
 preserve_case = false
 numeric_index = false
 longname = ""
-description = "keyword to indicate release of particles at the start of all time steps in the period."
+description = "keyword to indicate release at the start of all time steps in the period."
 deprecated = ""
 
 [block.period.first]
@@ -615,7 +598,24 @@ optional = false
 preserve_case = false
 numeric_index = false
 longname = ""
-description = "keyword to indicate release of particles at the start of the first time step in the period. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps."
+description = "keyword to indicate release at the start of the first time step in the period. This keyword may be used in conjunction with other RELEASESETTING options."
+deprecated = ""
+
+[block.period.last]
+type = "keyword"
+block_variable = false
+valid = []
+shape = ""
+tagged = true
+in_record = true
+layered = false
+time_series = false
+reader = "urword"
+optional = false
+preserve_case = false
+numeric_index = false
+longname = ""
+description = "keyword to indicate release at the start of the last time step in the period. This keyword may be used in conjunction with other RELEASESETTING options."
 deprecated = ""
 
 [block.period.frequency]
@@ -632,7 +632,7 @@ optional = false
 preserve_case = false
 numeric_index = false
 longname = ""
-description = "release particles at the specified time step frequency. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps."
+description = "release at the specified time step frequency. This keyword may be used in conjunction with other RELEASESETTING options."
 deprecated = ""
 
 [block.period.steps]
@@ -649,7 +649,7 @@ optional = false
 preserve_case = false
 numeric_index = false
 longname = ""
-description = "release particles at the start of each step specified in STEPS. This keyword may be used in conjunction with other keywords to release particles at the start of multiple time steps."
+description = "release at the start of each step specified in STEPS. This option may be used in conjunction with other RELEASESETTING options."
 deprecated = ""
 
 [block.period.fraction]
@@ -666,7 +666,7 @@ optional = true
 preserve_case = false
 numeric_index = false
 longname = ""
-description = "release particles after the specified fraction of the time step has elapsed. If FRACTION is not set, particles are released at the start of the specified time step(s). FRACTION must be a single value when used with ALL, FIRST, or FREQUENCY. When used with STEPS, FRACTION may be a single value or an array of the same length as STEPS. If a single FRACTION value is provided with STEPS, the fraction applies to all steps."
+description = "release particles after the specified fraction of the time step has elapsed. If FRACTION is not set, particles are released at the start of the specified time step(s). FRACTION must be a single value when used with ALL, FIRST, or FREQUENCY. When used with STEPS, FRACTION may be a single value or an array of the same length as STEPS. If a single FRACTION value is provided with STEPS, the fraction applies to all steps. NOTE: The FRACTION option has been removed. For fine control over release timing, specify times explicitly using the RELEASETIMES block."
 deprecated = ""
 
 [block.period.perioddata]
diff --git a/test/test_array.py b/test/test_array.py
index cb4ac14..a366964 100644
--- a/test/test_array.py
+++ b/test/test_array.py
@@ -13,7 +13,7 @@ def test_array_load_1d(tmp_path):
     with open(fpth, "w") as f:
         f.write(f"{name.upper()}\n{how}\n{value}\n")
     with open(fpth, "r") as f:
-        array = MFArray.load(f, cwd=tmp_path, shape=(3))
+        array = MFArray.load(f, cwd=tmp_path, shape=(3), type="double")
         assert array.name == name
         assert np.allclose(array.value, np.array(v))
 
@@ -24,6 +24,7 @@ def test_array_load_3d(tmp_path):
     how = "INTERNAL"
     v = [[[1.0, 2.0, 3.0]], [[4.0, 5.0, 6.0]], [[7.0, 8.0, 9.0]]]
     value = ""
+
     for a in v:
         for b in a:
             value += " ".join(str(x) for x in b)
@@ -32,7 +33,7 @@ def test_array_load_3d(tmp_path):
     with open(fpth, "w") as f:
         f.write(f"{name.upper()}\n{how}\n{value}\n")
     with open(fpth, "r") as f:
-        array = MFArray.load(f, cwd=tmp_path, shape=(3, 1, 3))
+        array = MFArray.load(f, cwd=tmp_path, shape=(3, 1, 3), type="double")
         assert array.name == name
         assert np.allclose(array.value, np.array(v))
 
@@ -44,6 +45,7 @@ def test_array_load_3d_external(tmp_path):
     how = "OPEN/CLOSE"
     v = [[[1.0, 2.0, 3.0]], [[4.0, 5.0, 6.0]], [[7.0, 8.0, 9.0]]]
     value = ""
+
     for a in v:
         for b in a:
             value += " ".join(str(x) for x in b)
@@ -54,7 +56,7 @@ def test_array_load_3d_external(tmp_path):
     with open(fpth, "w") as f:
         f.write(f"{name.upper()}\n{how} {extfpth}\n")
     with open(fpth, "r") as f:
-        array = MFArray.load(f, cwd=tmp_path, shape=(3, 1, 3))
+        array = MFArray.load(f, cwd=tmp_path, shape=(3, 1, 3), type="double")
         assert array.name == name
         assert np.allclose(array.value, np.array(v))
 
@@ -65,6 +67,7 @@ def test_array_load_layered(tmp_path):
     how = "INTERNAL"
     v = [[[1.0, 2.0, 3.0]], [[4.0, 5.0, 6.0]], [[7.0, 8.0, 9.0]]]
     value = ""
+
     for a in v:
         for b in a:
             value += "  " + f"{how}\n"
@@ -74,6 +77,6 @@ def test_array_load_layered(tmp_path):
     with open(fpth, "w") as f:
         f.write(f"{name.upper()} LAYERED\n{value}")
     with open(fpth, "r") as f:
-        array = MFArray.load(f, cwd=tmp_path, shape=(3, 1, 3))
+        array = MFArray.load(f, cwd=tmp_path, shape=(3, 1, 3), type="double")
         assert array.name == name
         assert np.allclose(array.value, np.array(v))
diff --git a/test/test_dfn.py b/test/test_dfn.py
index bb83e39..03c5feb 100644
--- a/test/test_dfn.py
+++ b/test/test_dfn.py
@@ -18,8 +18,14 @@ def test_dfn_load(tmp_path):
     assert dfn.component == "prt"
     assert dfn.subcomponent == "prp"
     assert type(dfn.dfn) is dict
-    assert len(dfn) == 4
-    assert dfn.blocknames == ["options", "dimensions", "packagedata", "period"]
+    assert len(list(dfn.dfn["block"])) == 5
+    assert dfn.blocknames == [
+        "options",
+        "dimensions",
+        "packagedata",
+        "releasetimes",
+        "period",
+    ]
 
     for b in dfn.blocknames:
         block_d = dfn[b]
@@ -44,13 +50,9 @@ def test_dfn_load(tmp_path):
         "stop_at_weak_sink",
         "istopzone",
         "drape",
-        "release_timesrecord",
-        "release_times",
-        "times",
-        "release_timesfilerecord",
-        "release_timesfile",
-        "timesfile",
         "dev_forceternary",
+        "release_time_tolerance",
+        "release_time_frequency",
     ]
 
     assert dfn.param("options", "drape") == {
diff --git a/test/test_list.py b/test/test_list.py
index 5e17851..a44f324 100644
--- a/test/test_list.py
+++ b/test/test_list.py
@@ -11,8 +11,8 @@ class TestBlock(MFBlock):
     testblock = MFList(
         params={
             "s": MFString(),
-            "i": MFInteger(),
-            "d": MFDouble(),
+            "i": MFInteger(type="integer"),
+            "d": MFDouble(type="double"),
         },
         description="recarray",
         optional=False,
@@ -67,12 +67,13 @@ def test_list_load1(tmp_path):
     assert isinstance(TestBlock.testblock.params["i"], MFInteger)
     assert isinstance(TestBlock.testblock.params["d"], MFDouble)
     assert isinstance(in_list.params["testblock"]["s"], list)
-    assert isinstance(in_list.params["testblock"]["i"], np.ndarray)
-    assert isinstance(in_list.params["testblock"]["d"], np.ndarray)
+    assert isinstance(in_list.params["testblock"]["i"], list)
+    assert isinstance(in_list.params["testblock"]["d"], list)
     assert isinstance(in_list.params["testblock"]["s"][0], str)
-    assert isinstance(in_list.params["testblock"]["i"][0], np.int32)
-    assert isinstance(in_list.params["testblock"]["d"][0], np.float64)
+    assert isinstance(in_list.params["testblock"]["i"][0], int)
+    assert isinstance(in_list.params["testblock"]["d"][0], float)
     assert in_list.params["testblock"]["s"] == ["model", "exch", "sim"]
+    print(in_list.params["testblock"]["i"])
     assert np.allclose(in_list.params["testblock"]["i"], np.array([1, 2, 2]))
     assert np.allclose(
         in_list.params["testblock"]["d"], np.array([2.0, 3.0, 3.0])
diff --git a/test/test_package.py b/test/test_package.py
index 9e74d4e..9d186b0 100644
--- a/test/test_package.py
+++ b/test/test_package.py
@@ -34,7 +34,7 @@ class TestPackage(MFPackage):
         optional=False,
     )
     a = MFArray(
-        block="packagedata", type="array", description="array", shape=(3)
+        block="packagedata", type="double", description="array", shape=(3)
     )
 
 
@@ -73,6 +73,7 @@ class TestGwfIc(MFPackage):
         optional=False,
         # shape="(nodes)",
         shape=(10),
+        type="double",
     )
 
 
@@ -112,6 +113,7 @@ class TestGwfDis(MFPackage):
         optional=False,
         # shape="(ncol)",
         shape=(5),
+        type="double",
     )
     delc = MFArray(
         block="griddata",
@@ -120,6 +122,7 @@ class TestGwfDis(MFPackage):
         optional=False,
         # shape="(nrow)",
         shape=(5),
+        type="double",
     )
     top = MFArray(
         block="griddata",
@@ -129,6 +132,7 @@ class TestGwfDis(MFPackage):
         optional=False,
         # shape="(ncol, nrow)",
         shape=(5, 5),
+        type="double",
     )
     botm = MFArray(
         block="griddata",
@@ -137,6 +141,7 @@ class TestGwfDis(MFPackage):
         optional=False,
         # shape="(ncol, nrow, nlay)",
         shape=(3, 5, 5),
+        type="double",
     )
     idomain = MFArray(
         block="griddata",
@@ -155,6 +160,7 @@ class TestGwfDis(MFPackage):
         optional=False,
         # shape="(ncol, nrow, nlay)",
         shape=(3, 5, 5),
+        type="integer",
     )
 
 
diff --git a/test/test_sim.py b/test/test_sim.py
index e97f6b9..382285c 100644
--- a/test/test_sim.py
+++ b/test/test_sim.py
@@ -1,3 +1,6 @@
+import os
+import subprocess
+
 import numpy as np
 
 from flopy4.simulation import MFSimulation
@@ -26,7 +29,7 @@ def test_load_sim(tmp_path):
         f.write("    CONSTANT  -100.00000000\n")
         f.write("    CONSTANT  -150.00000000\n")
         f.write("    CONSTANT  -350.00000000\n")
-        f.write("END GRIDDATA\nn\\n")
+        f.write("END GRIDDATA\n")
 
     ic_fpth = tmp_path / f"{name}.ic"
     strt = np.linspace(0.0, 30.0, num=300)
@@ -107,7 +110,7 @@ def test_load_sim(tmp_path):
         f.write("  INNER_MAXIMUM  300\n")
         f.write("  INNER_DVCLOSE  1.00000000E-09\n")
         # TODO: fails
-        # f.write("  inner_rclose  1.00000000E-06\n")
+        # f.write("  INNER_RCLOSE  1.00000000E-06\n")
         f.write("  LINEAR_ACCELERATION  bicgstab\n")
         f.write("  RELAXATION_FACTOR       1.00000000\n")
         f.write("  SCALING_METHOD  none\n")
@@ -141,15 +144,15 @@ def test_load_sim(tmp_path):
     assert "exchanges" in s.nam
     assert "solutiongroup" in s.nam
     assert "tdis6" in s.nam["timing"]
-    # assert s.nam["timing"]["tdis6"].value == f"{tmp_path}/sim.tdis"
+    assert s.nam["timing"]["tdis6"].value == f"{tmp_path}/sim.tdis"
     assert "mtype" in s.nam["models"].params["models"]
     assert "mfname" in s.nam["models"].params["models"]
     assert "mname" in s.nam["models"].params["models"]
     assert s.nam["models"].params["models"]["mtype"][0] == "GWF6"
-    # assert (
-    #    s.nam["models"].params["models"]["mfname"][0]
-    #    == f"{tmp_path}/{name}.nam"
-    # )
+    assert (
+        s.nam["models"].params["models"]["mfname"][0]
+        == f"{tmp_path}/{name}.nam"
+    )
     assert s.nam["models"].params["models"]["mname"][0] == f"{name}"
     assert "slntype" in s.nam["solutiongroup"].params["solutiongroup"]
     assert "slnfname" in s.nam["solutiongroup"].params["solutiongroup"]
@@ -157,10 +160,10 @@ def test_load_sim(tmp_path):
     assert (
         s.nam["solutiongroup"].params["solutiongroup"]["slntype"][0] == "ims6"
     )
-    # assert (
-    #    s.nam["solutiongroup"].params["solutiongroup"]["slnfname"][0]
-    #    == f"{tmp_path}/{name}.ims"
-    # )
+    assert (
+        s.nam["solutiongroup"].params["solutiongroup"]["slnfname"][0]
+        == f"{tmp_path}/{name}.ims"
+    )
     assert (
         s.nam["solutiongroup"].params["solutiongroup"]["slnmnames"][0]
         == f"{name}"
@@ -338,8 +341,146 @@ def test_load_sim(tmp_path):
         f"sim/{name}/dis/dimensions/ncol"
     )
 
-    import os
+    write_dir = tmp_path / "write"
+    os.makedirs(write_dir)
+    s.write(write_dir)
+
+
+def test_load_chd01(tmp_path):
+    name = "gwf_chd01"
+
+    nlay = 1
+    nrow = 1
+    ncol = 100
+    dis_fpth = tmp_path / f"{name}.dis"
+    with open(dis_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("END OPTIONS\n\n")
+        f.write("BEGIN DIMENSIONS\n")
+        f.write(f"  NLAY  {nlay}\n")
+        f.write(f"  NROW  {nrow}\n")
+        f.write(f"  NCOL  {ncol}\n")
+        f.write("END DIMENSIONS\n\n")
+        f.write("BEGIN GRIDDATA\n")
+        f.write("  DELR\n    CONSTANT  1.00000000\n")
+        f.write("  DELC\n    CONSTANT  1.00000000\n")
+        f.write("  TOP\n    CONSTANT  1.00000000\n")
+        f.write("  BOTM\n    CONSTANT  0.00000000\n")
+        f.write("  IDOMAIN\n  INTERNAL FACTOR 1\n")
+        f.write(
+            "    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1"
+            " 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1"
+            " 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1"
+            " 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1"
+            " 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1\n"
+        )
+        f.write("END GRIDDATA\n")
+
+    ic_fpth = tmp_path / f"{name}.ic"
+    with open(ic_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("END OPTIONS\n\n")
+        f.write("BEGIN GRIDDATA\n")
+        f.write("  STRT\n    CONSTANT  1.00000000\n")
+        f.write("END GRIDDATA\n")
+
+    npf_fpth = tmp_path / f"{name}.npf"
+    with open(npf_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("  SAVE_SPECIFIC_DISCHARGE\n")
+        f.write("END OPTIONS\n\n")
+        f.write("BEGIN GRIDDATA\n")
+        f.write("  ICELLTYPE\n    CONSTANT  0\n")
+        f.write("  K\n    CONSTANT  1.00000000\n")
+        f.write("  K33\n    CONSTANT  1.00000000\n")
+        f.write("END GRIDDATA\n")
+
+    chd_fpth = tmp_path / f"{name}.chd"
+    with open(chd_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("  PRINT_FLOWS\n")
+        f.write("END OPTIONS\n\n")
+        f.write("BEGIN DIMENSIONS\n")
+        f.write("  MAXBOUND  2\n")
+        f.write("END DIMENSIONS\n\n")
+        f.write("BEGIN PERIOD 1\n")
+        f.write("  1 1 1 1.00000000E+00\n")
+        f.write("  1 1 100 0.00000000E+00\n")
+        f.write("END PERIOD 1\n")
+
+    nam_fpth = tmp_path / f"{name}.nam"
+    with open(nam_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("  SAVE_FLOWS\n")
+        f.write("END OPTIONS\n")
+        f.write("\n")
+        f.write("BEGIN PACKAGES\n")
+        f.write(f"  DIS6  {name}.dis  dis\n")
+        f.write(f"  IC6  {name}.ic  ic\n")
+        f.write(f"  NPF6  {name}.npf  npf\n")
+        f.write(f"  CHD6  {name}.chd  chd-1\n")
+        # f.write(f"  OC6  {name}.oc  oc\n")
+        f.write("END PACKAGES\n")
+
+    tdis_fpth = tmp_path / "chd01.tdis"
+    with open(tdis_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("  TIME_UNITS  days\n")
+        f.write("END OPTIONS\n\n")
+        f.write("BEGIN DIMENSIONS\n")
+        f.write("  NPER  1\n")
+        f.write("END DIMENSIONS\n\n")
+        f.write("BEGIN PERIODDATA\n")
+        f.write("  5.00000000  1       1.00000000\n")
+        f.write("END PERIODDATA\n\n")
+
+    ims_fpth = tmp_path / f"{name}.ims"
+    with open(ims_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("  PRINT_OPTION  summary\n")
+        f.write("END OPTIONS\n\n")
+        f.write("BEGIN NONLINEAR\n")
+        f.write("  OUTER_DVCLOSE  1.00000000E-06\n")
+        f.write("  OUTER_MAXIMUM  100\n")
+        f.write("  UNDER_RELAXATION  none\n")
+        f.write("END NONLINEAR\n\n")
+        f.write("BEGIN LINEAR\n")
+        f.write("  INNER_MAXIMUM  300\n")
+        f.write("  INNER_DVCLOSE  1.00000000E-06\n")
+        # TODO: fails
+        # f.write("  inner_rclose  1.00000000E-06\n")
+        f.write("  LINEAR_ACCELERATION  cg\n")
+        f.write("  RELAXATION_FACTOR       1.00000000\n")
+        f.write("  SCALING_METHOD  none\n")
+        f.write("  REORDERING_METHOD  none\n")
+        f.write("END LINEAR\n\n")
+
+    sim_fpth = tmp_path / "mfsim.nam"
+    with open(sim_fpth, "w") as f:
+        f.write("BEGIN OPTIONS\n")
+        f.write("END OPTIONS\n\n")
+        f.write("BEGIN TIMING\n")
+        f.write("  TDIS6  chd01.tdis\n")
+        f.write("END TIMING\n\n")
+        f.write("BEGIN MODELS\n")
+        f.write(f"  GWF6  {name}.nam  {name}\n")
+        f.write("END MODELS\n\n")
+        f.write("BEGIN EXCHANGES\n")
+        f.write("END EXCHANGES\n\n")
+        f.write("BEGIN SOLUTIONGROUP 1\n")
+        f.write(f"  ims6  {name}.ims  {name}\n")
+        f.write("END SOLUTIONGROUP 1\n\n")
+
+    s = None
+    with open(sim_fpth, "r") as f:
+        s = MFSimulation.load(f)
 
     write_dir = tmp_path / "write"
     os.makedirs(write_dir)
     s.write(write_dir)
+
+    os.chdir(write_dir)
+    s = subprocess.run(["which", "mf6"])
+    if s.returncode == 0:
+        subprocess.run(["mf6"])
+        subprocess.run(["diff", f"./{name}.lst", f"../{name}.lst"])