diff --git a/flopy4/mf6/netcdf.py b/flopy4/mf6/netcdf.py
index 8cd61ee3..65ef12d3 100644
--- a/flopy4/mf6/netcdf.py
+++ b/flopy4/mf6/netcdf.py
@@ -24,22 +24,39 @@
 from flopy4.version import __version__
 
 
-def _cf_var_attrs(dims: list[str], mesh: str | None, grid) -> dict:
+def _cf_var_attrs(dims: list[str], mesh: str | None, grid, layer: int | None = None) -> dict:
     """Return {"attrs": {...}, "encoding": {...}} for a NetCDF data variable.
 
     In mesh context x/y are abstract row/col indices, not geographic — only nmesh_face
     vars get coordinates linking. In structured context x/y ARE geographic.
+
+    layer is the 1-based layer number for a per-layer variable, or None
+    otherwise. z_l{layer} is appended to coordinates only when set -- the
+    only case it shares nmesh_face with the variable (CF-1.13 5.2).
     """
     attrs: dict[str, str] = {}
     encoding: dict[str, object] = {}
     has_crs = grid is not None and getattr(grid, "crs", None) is not None
+    # grid_mapping is CF-required only on variables spanning the full 2D
+    # horizontal spatial extent (CF-1.13 5.6); 1D dimension-definition arrays
+    # such as dis_delr/dis_delc are not georeferenced fields and must not
+    # carry it, matching MF6's own ncvar_gridmap (DisNCStructured.f90).
+    has_full_extent = ("x" in dims and "y" in dims) or "nmesh_face" in dims
 
-    if has_crs:
+    if has_crs and has_full_extent:
         attrs["grid_mapping"] = "projection"
     if "nmesh_face" in dims:
-        attrs["coordinates"] = "mesh_face_x mesh_face_y"
+        if layer:
+            attrs["coordinates"] = f"mesh_face_x mesh_face_y z_l{layer}"
+        else:
+            attrs["coordinates"] = "mesh_face_x mesh_face_y"
         attrs["mesh"] = "mesh"
         attrs["location"] = "face"
+    elif "layer" in dims:
+        # Structured per-layer fields link to the z (cell center elevation)
+        # coordinate, matching MF6's own DisNCStructured.f90 -- non-layered
+        # fields (e.g. dis_top, a y/x-only 2D array) do not.
+        attrs["coordinates"] = "z"
 
     return {"attrs": attrs, "encoding": encoding}
 
@@ -209,12 +226,22 @@ def from_model(
             raise ValueError("model must have a 'name' attribute")
 
         modeltype = model.__class__.__name__.lower()
-        attrs = {"title": f"{model.name.upper()} model input"}
+        # Matches MF6's own NCModel.f90 title convention (model-type-specific
+        # fragment + " array input" for the pre-solve/validate-mode INPUT
+        # file this class always produces); models MF6 itself doesn't
+        # support for NetCDF export (anything but GWF/GWT/GWE) fall back to
+        # a generic title.
+        _title_fragment = {"gwf": "hydraulic head", "gwt": "concentration", "gwe": "temperature"}
+        if modeltype in _title_fragment:
+            title = f"{model.name.upper()} {_title_fragment[modeltype]} array input"
+        else:
+            title = f"{model.name.upper()} model input"
+        attrs = {"title": title}
         packages = []
         distype = None
 
         if netcdf_format == NetCDFFormat.LAYERED_MESH:
-            attrs["mesh"] = NetCDFFormat.LAYERED_MESH.value
+            attrs["modflow_mesh"] = NetCDFFormat.LAYERED_MESH.value
 
         # Resolve nper from time arg, simulation tdis, or model's data dims.
         if time is not None:
@@ -260,11 +287,12 @@ def from_model(
                     val = getattr(package, f.name)
                     if val is None:
                         continue
-                    arr = np.asarray(val, dtype=np.float64)
+                    _dtype = _PKG_DTYPE_MAP.get(to_field_type(f.type), np.float64)
+                    arr = np.asarray(val, dtype=_dtype)
                     # Only broadcast scalars to full grid for nodes-shaped fields
                     shape_meta = f.metadata.get("shape", ())
                     if "nodes" in shape_meta and arr.size < _nodes:
-                        arr = np.full(_nodes, float(arr.ravel()[0]))
+                        arr = np.full(_nodes, arr.ravel()[0], dtype=_dtype)
                     p["params"].append({"name": f.name, "data": arr})
                 elif f.metadata.get("fk"):
                     # dynamically named arrays (RCHA's aux): one param per auxiliary name
@@ -285,6 +313,14 @@ def from_model(
                                 "data": dense(named_periods, _nper, carry_forward=False),
                             }
                         )
+                elif (
+                    f.metadata.get("block") == "period" and f.metadata.get("reader") == "readarray"
+                ):
+                    val = getattr(package, f.name)
+                    if val is None:
+                        continue
+                    _dtype = _PKG_DTYPE_MAP.get(to_field_type(f.type), np.float64)
+                    p["params"].append({"name": f.name, "data": np.asarray(val, dtype=_dtype)})
                 else:
                     # period arrays: {kper: array}, periods not given filled.
                     # Time-array series references stay in the package file.
@@ -345,12 +381,12 @@ def to_xarray(self) -> xr.Dataset:
         meta = self.model_dump(by_alias=True)
 
         if self._grid is not None and self._time is not None:  # type: ignore
-            conventions = "CF-1.11"  # type: ignore
-            if meta["attrs"]["mesh"] is not None:
+            conventions = "CF-1.13"  # type: ignore
+            if meta["attrs"]["modflow_mesh"] is not None:
                 conventions = f"{conventions} UGRID-1.0"
             _fmt = (
                 NetCDFFormat.LAYERED_MESH
-                if meta["attrs"]["mesh"] is not None
+                if meta["attrs"]["modflow_mesh"] is not None
                 else NetCDFFormat.STRUCTURED
             )
             dss.append(self._grid.to_xarray(netcdf_format=_fmt, modeltime=self._time))
@@ -413,7 +449,13 @@ def validate_attrs(cls, v: dict[str, str]) -> dict[str, str]:
         """
         validate model (dataset) scoped attributes dictionary
         """
+        # title is free-text (e.g. "GWFMODEL hydraulic head array input") and must
+        # keep its original casing, matching MF6 -- unlike the other keys/values
+        # here, which are lowercase identifiers by convention.
+        title = v.get("title")
         v = lower(v)
+        if title is not None:
+            v["title"] = title
         return v
 
     @staticmethod
@@ -431,7 +473,9 @@ def _backfill_meta(meta: dict, context: dict, verbose: bool = True) -> dict:
 
         _packages = []
         for pkg in _meta["packages"]:
-            pkgctx = {"mesh": _meta["attrs"]["mesh"]} if "mesh" in _meta["attrs"] else {}
+            pkgctx = (
+                {"mesh": _meta["attrs"]["modflow_mesh"]} if "modflow_mesh" in _meta["attrs"] else {}
+            )
             pkgctx["modelname"] = _meta["modelname"]
             pkgctx["gridtype"] = _meta["gridtype"]
             pkgctx |= context
@@ -443,17 +487,17 @@ def _backfill_meta(meta: dict, context: dict, verbose: bool = True) -> dict:
 
 class NetCDFModelAttrs(BaseModel):
     # order of params dictates when data added to info dict
-    mesh: str | None = Field(default=None)
+    modflow_mesh: str | None = Field(default=None)
     modflow_grid: str = Field()
     modflow_model: str = Field()
 
     model_config = ConfigDict(extra="allow")
 
-    @field_validator("mesh", mode="before")
+    @field_validator("modflow_mesh", mode="before")
     @classmethod
-    def validate_mesh(cls, v: str | None, info: ValidationInfo) -> str | None:
+    def validate_modflow_mesh(cls, v: str | None, info: ValidationInfo) -> str | None:
         """
-        validate model mesh attribute
+        validate model modflow_mesh attribute
         """
         if v is not None:
             if v.lower() != "layered":
@@ -728,6 +772,7 @@ def to_xarray(self) -> xr.Dataset:
             [str(d) for d in ds[varname].dims],
             mesh,
             self._context.get("grid"),
+            layer=meta["attrs"].get("layer"),
         )
         ds[varname].attrs.update(cf["attrs"])
         ds[varname].encoding.update(cf["encoding"])
diff --git a/flopy4/mf6/utils/crs.py b/flopy4/mf6/utils/crs.py
new file mode 100644
index 00000000..fb0f9aca
--- /dev/null
+++ b/flopy4/mf6/utils/crs.py
@@ -0,0 +1,149 @@
+from typing import Any
+
+import numpy as np
+from pyproj import CRS
+
+# CF-standard numeric grid_mapping parameters shared across all three
+# supported projected grid_mapping_name values.
+_SHARED_CF_PARAMS = (
+    "longitude_of_central_meridian",
+    "latitude_of_projection_origin",
+    "false_easting",
+    "false_northing",
+    "semi_major_axis",
+    "inverse_flattening",
+)
+
+
+def cf_grid_mapping_params(crs: CRS) -> tuple[dict[str, Any], str | None]:
+    """
+    Return CF-standard numeric grid_mapping parameters for *crs*, in
+    addition to wkt/crs_wkt: some netCDF readers -- notably ArcGIS's
+    classic netCDF connector -- position a projected grid using these
+    individual attributes rather than parsing wkt/crs_wkt text, and
+    default every parameter to 0 if absent.
+
+    Covers transverse_mercator, lambert_conformal_conic (2SP only -- CF
+    has no scale_factor attribute for this method, so a 1SP CRS cannot
+    be exactly represented), and albers_conical_equal_area. Always
+    derive these from the original, unwrapped CRS: pyproj.CRS.to_cf()
+    does not "see through" a wrapped DerivedProjectedCRS.
+
+    Parameters
+    ----------
+    crs : pyproj.CRS
+        The (unwrapped) CRS to extract parameters from.
+
+    Returns
+    -------
+    tuple[dict, str | None]
+        The parameters (empty if the grid_mapping_name is unsupported or
+        unrecognized), and a warning message if the grid_mapping_name is
+        recognized but unsupported (currently: 1SP lambert_conformal_conic
+        variants), else None.
+    """
+    cf = crs.to_cf()
+    gmn = cf.get("grid_mapping_name")
+
+    if gmn == "transverse_mercator":
+        keys = ("scale_factor_at_central_meridian", *_SHARED_CF_PARAMS)
+    elif gmn in ("lambert_conformal_conic", "albers_conical_equal_area"):
+        method = crs.coordinate_operation.method_name if crs.coordinate_operation else ""
+        if "1SP" in method:
+            return {}, (
+                "1SP lambert_conformal_conic variants are not supported for "
+                "CF grid_mapping parameter extraction (CF has no "
+                "scale_factor attribute for this method); some netCDF "
+                "readers may not correctly position this grid."
+            )
+        keys = ("standard_parallel", *_SHARED_CF_PARAMS)
+    else:
+        return {}, None
+
+    params: dict[str, Any] = {}
+    for k in keys:
+        if k in cf:
+            params[k] = np.asarray(cf[k]) if k == "standard_parallel" else cf[k]
+    return params, None
+
+
+def wrap_rotated_crs(
+    crs: CRS,
+    xorigin: float,
+    yorigin: float,
+    angrot: float,
+) -> CRS | None:
+    """
+    Wrap a projected CRS in a derived CRS encoding MF6 grid rotation.
+
+    Builds a DerivedProjectedCRS over *crs* with an EPSG:9624 (affine
+    parametric transformation) deriving conversion parameterized from
+    xorigin/yorigin/angrot (degrees), matching MF6's DisNCStructured.f90
+    CRS_WKT rotation-wrapping exactly. Per ISO 19111, a deriving
+    conversion is directed base->derived, so the parameters encode the
+    world->local (inverse) rotation; a CRS-aware consumer applies the
+    inverse to resolve true position from local coordinates.
+
+    Parameters
+    ----------
+    crs : pyproj.CRS
+        The base CRS to wrap. Must be a projected CRS.
+    xorigin, yorigin : float
+        Grid origin in *crs*.
+    angrot : float
+        Grid rotation angle in degrees.
+
+    Returns
+    -------
+    pyproj.CRS | None
+        The wrapped CRS, or None if *crs* is not a projected CRS.
+    """
+    if not crs.is_projected:
+        return None
+
+    ang = np.radians(angrot)
+    a0 = -(xorigin * np.cos(ang) + yorigin * np.sin(ang))
+    a1 = np.cos(ang)
+    a2 = np.sin(ang)
+    b0 = xorigin * np.sin(ang) - yorigin * np.cos(ang)
+    b1 = -np.sin(ang)
+    b2 = np.cos(ang)
+
+    derived = {
+        "type": "DerivedProjectedCRS",
+        "name": "MODFLOW 6 rotated grid CRS",
+        "base_crs": crs.to_json_dict(),
+        "conversion": {
+            "name": "MODFLOW 6 grid rotation",
+            "method": {
+                "name": "Affine parametric transformation",
+                "id": {"authority": "EPSG", "code": 9624},
+            },
+            "parameters": [
+                {"name": "A0", "value": a0, "unit": "metre"},
+                {"name": "A1", "value": a1, "unit": "unity"},
+                {"name": "A2", "value": a2, "unit": "unity"},
+                {"name": "B0", "value": b0, "unit": "metre"},
+                {"name": "B1", "value": b1, "unit": "unity"},
+                {"name": "B2", "value": b2, "unit": "unity"},
+            ],
+        },
+        "coordinate_system": {
+            "subtype": "Cartesian",
+            "axis": [
+                {
+                    "name": "Easting",
+                    "abbreviation": "X",
+                    "direction": "east",
+                    "unit": "metre",
+                },
+                {
+                    "name": "Northing",
+                    "abbreviation": "Y",
+                    "direction": "north",
+                    "unit": "metre",
+                },
+            ],
+        },
+    }
+    return CRS.from_json_dict(derived)
diff --git a/flopy4/mf6/utils/grid.py b/flopy4/mf6/utils/grid.py
index 24f124fd..b57a2737 100644
--- a/flopy4/mf6/utils/grid.py
+++ b/flopy4/mf6/utils/grid.py
@@ -563,6 +563,21 @@ def _layered_mesh_dataset(self, ds, modeltime=None, configuration=None):
         ds["layer"].attrs["axis"] = "Z"
         ds["layer"].encoding["_FillValue"] = None
 
+        # z_l1, z_l2, ...: per-layer cell center elevation, always written.
+        # Split per layer because face-indexed variables have no layer
+        # dimension to reference a combined z(layer, nmesh_face) (CF-1.13 5.2).
+        _z = self._coords["z"].values
+        for k in range(self.nlay):
+            _varname = f"z_l{k + 1}"
+            ds = ds.assign({_varname: (["nmesh_face"], _z[k].flatten())})
+            if _units is not None:
+                ds[_varname].attrs["units"] = _units
+            ds[_varname].attrs["standard_name"] = "altitude"
+            ds[_varname].attrs["positive"] = "up"
+            ds[_varname].attrs["long_name"] = f"cell center elevation (layer {k + 1})"
+            ds[_varname].attrs["layer"] = k + 1
+            ds[_varname].encoding["_FillValue"] = None
+
         # mesh container variable
         ds = ds.assign({"mesh": ([], np.int64(1))})
         ds["mesh"].attrs["cf_role"] = "mesh_topology"
@@ -610,8 +625,11 @@ def _layered_mesh_dataset(self, ds, modeltime=None, configuration=None):
         ds["mesh_face_y"].attrs["long_name"] = "Northing"
         ds["mesh_face_y"].attrs["bounds"] = "mesh_face_ybnds"
         ds["mesh_face_y"].encoding["_FillValue"] = None
-        ds["mesh_face_xbnds"].encoding["_FillValue"] = FILL_FLOAT64
-        ds["mesh_face_ybnds"].encoding["_FillValue"] = FILL_FLOAT64
+        # No explicit _FillValue attribute: matches MF6 (MeshNCModel.f90), which
+        # pads with NF90_FILL_DOUBLE in the data itself but never declares the
+        # attribute on mesh_face_xbnds/ybnds.
+        ds["mesh_face_xbnds"].encoding["_FillValue"] = None
+        ds["mesh_face_ybnds"].encoding["_FillValue"] = None
 
         # mesh face nodes
         var_d = {
@@ -654,13 +672,21 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
         # no UGRID mesh variable).
         _units = _coord_units(self.lenuni, self.crs)
 
-        xc = self.xoffset + self.xycenters[0]
-        yc = self.yoffset + self.xycenters[1]
-        # z = [float(x) for x in range(1, self.nlay + 1)]
+        # For a rotated grid, x/y (and their bounds) are grid-local -- zero
+        # offset, not xoffset/yoffset -- matching MF6's DisNCStructured.f90
+        # convention exactly. wrap_rotated_crs's derived-CRS parameters are
+        # anchored to this same zero-based local coordinate system, so x/y
+        # must stay zero-based here for that wrap to resolve the correct
+        # real-world position; offsetting them would silently break it.
+        _xoff = 0.0 if self.angrot else self.xoffset
+        _yoff = 0.0 if self.angrot else self.yoffset
+
+        xc = _xoff + self.xycenters[0]
+        yc = _yoff + self.xycenters[1]
 
         # set coordinate var bounds
         x_bnds = []
-        xv = self.xoffset + self.xyedges[0]
+        xv = _xoff + self.xyedges[0]
         for idx, val in enumerate(xv):
             if idx + 1 < len(xv):
                 bnd = []
@@ -669,7 +695,7 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
                 x_bnds.append(bnd)
 
         y_bnds = []
-        yv = self.yoffset + self.xyedges[1]
+        yv = _yoff + self.xyedges[1]
         for idx, val in enumerate(yv):
             if idx + 1 < len(yv):
                 bnd = []
@@ -698,16 +724,62 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
         ds["time"].attrs["standard_name"] = "time"
         ds["time"].attrs["long_name"] = "time"
         ds["time"].encoding["_FillValue"] = None
+        # Resolve the CRS (if any) once, up front, and attempt the
+        # rotation wrap once -- both the x/y axis/standard_name decision
+        # below and the projection variable's crs_wkt (further down) need
+        # the result, and pyproj.CRS.to_cf()/is_projected do not "see
+        # through" a wrapped DerivedProjectedCRS, so grid_mapping_name and
+        # any CF numeric grid_mapping parameters must always come from the
+        # original, unwrapped CRS.
+        _crs = None
+        if (
+            configuration is not None
+            and "wkt" in configuration
+            and configuration["wkt"] is not None
+        ):
+            from pyproj import CRS as ProjCRS
+
+            _crs = ProjCRS.from_wkt(configuration["wkt"])
+        elif self.crs is not None:
+            _crs = self.crs
+
+        _wrapped_crs = None
+        if self.angrot:
+            import warnings
+
+            from flopy4.mf6.utils.crs import wrap_rotated_crs
+
+            if _crs is not None:
+                _wrapped_crs = wrap_rotated_crs(_crs, self.xoffset, self.yoffset, self.angrot)
+            if _wrapped_crs is None:
+                warnings.warn(
+                    "rotated grid: no rotation-positioning information will "
+                    "be written to the output CRS; x/y coordinate values "
+                    "will not reflect the true rotated position "
+                    + (
+                        "(CRS is not a projected CRS)"
+                        if _crs is not None
+                        else "(no CRS configured)"
+                    )
+                )
+
+        # axis/standard_name assert that x/y hold the true projected
+        # position, which only holds once crs_wkt encodes the rotation
+        # (see wrap_rotated_crs above); otherwise x/y are grid-local.
+        _axis_ok = not self.angrot or _wrapped_crs is not None
+
         if _units is not None:
             ds["y"].attrs["units"] = _units
-        ds["y"].attrs["axis"] = "Y"
-        ds["y"].attrs["standard_name"] = "projection_y_coordinate"
+        if _axis_ok:
+            ds["y"].attrs["axis"] = "Y"
+            ds["y"].attrs["standard_name"] = "projection_y_coordinate"
         ds["y"].attrs["long_name"] = "Northing"
         ds["y"].attrs["bounds"] = "y_bnds"
         if _units is not None:
             ds["x"].attrs["units"] = _units
-        ds["x"].attrs["axis"] = "X"
-        ds["x"].attrs["standard_name"] = "projection_x_coordinate"
+        if _axis_ok:
+            ds["x"].attrs["axis"] = "X"
+            ds["x"].attrs["standard_name"] = "projection_x_coordinate"
         ds["x"].attrs["long_name"] = "Easting"
         ds["x"].attrs["bounds"] = "x_bnds"
         ds["x"].encoding["_FillValue"] = None
@@ -720,29 +792,30 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
         ds["x_bnds"].encoding["_FillValue"] = None
         ds["y_bnds"].encoding["_FillValue"] = None
 
-        # Write projection variable whenever CRS is available.
+        # z: cell center elevation, distinct from the discrete layer index.
+        # Always written, independent of CRS/NCF configuration.
+        ds = ds.assign({"z": (["layer", "y", "x"], self._coords["z"].values)})
+        if _units is not None:
+            ds["z"].attrs["units"] = _units
+        ds["z"].attrs["standard_name"] = "altitude"
+        ds["z"].attrs["positive"] = "up"
+        ds["z"].attrs["long_name"] = "cell center elevation"
+        ds["z"].encoding["_FillValue"] = None
+
+        # Write projection variable whenever CRS is available (_crs and
+        # _wrapped_crs were already resolved above, for the x/y
+        # axis/standard_name decision).
         # Lat/lon auxiliary coordinates are intentionally omitted: GDAL-based
         # tools (QGIS, ArcGIS) misplace projected rasters when 2D lat/lon arrays
         # coexist with projected x/y dimension coordinates. The NCF subpackage
         # (separate text file) still carries lat/lon for MODFLOW 6 output.
-        _crs = None
-        if (
-            configuration is not None
-            and "wkt" in configuration
-            and configuration["wkt"] is not None
-        ):
-            from pyproj import CRS as ProjCRS
-
-            _crs = ProjCRS.from_wkt(configuration["wkt"])
-        elif self.crs is not None:
-            _crs = self.crs
-
         if _crs is not None:
-            from pyproj import CRS as ProjCRS  # noqa: F811 (already imported above if configured)
             from pyproj.enums import WktVersion
 
+            # wkt (WKT1) is always the original, unwrapped CRS -- WKT1 has
+            # no derived-CRS syntax, so it cannot encode grid rotation.
             _wkt1 = _crs.to_wkt(WktVersion.WKT1_GDAL)
-            _wkt2 = _crs.to_wkt(WktVersion.WKT2_2019)
+            _wkt2 = (_wrapped_crs or _crs).to_wkt(WktVersion.WKT2_2019)
             ds["x"].attrs["grid_mapping"] = "projection"
             ds["y"].attrs["grid_mapping"] = "projection"
             ds = ds.assign({"projection": ([], np.int64(1))})
@@ -752,32 +825,15 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
             if _gmn:
                 ds["projection"].attrs["grid_mapping_name"] = _gmn
 
-            # GDAL-compatible georeferencing: GDAL's CF driver does not build a
-            # geotransform from projected 1D coordinate variables. Writing these
-            # two attributes on the projection variable ensures correct placement
-            # in GDAL-based tools (QGIS, ArcGIS Pro via GDAL). MF6 ignores them.
-            # Note: GDAL reads GeoTransform from the grid_mapping variable, not
-            # from global attrs — NC_GLOBAL#GeoTransform is ignored for extent.
-            # Derive effective pixel sizes from the actual grid bounds so that
-            # variable-spacing grids get the correct bounding box in GDAL.
-            # Using x[1]-x[0] only works for uniform grids; outer cells in
-            # Frenchman-Flat-style grids are much coarser than interior cells.
-            _x_left = float(x_bnds[0][0])
-            _x_right = float(x_bnds[-1][1])
-            _y_top = float(y_bnds[0][1])
-            _y_bot = float(y_bnds[-1][0])
-            _dx_eff = (_x_right - _x_left) / len(xc)
-            _dy_eff = (_y_bot - _y_top) / len(yc)  # negative for north-up
-            _gt = [
-                _x_left,  # upper-left x
-                _dx_eff,  # effective x pixel size
-                0.0,
-                _y_top,  # upper-left y
-                0.0,
-                _dy_eff,  # effective y pixel size (negative for north-up)
-            ]
-            ds["projection"].attrs["GeoTransform"] = " ".join(str(v) for v in _gt)
-            ds["projection"].attrs["spatial_ref"] = _wkt1
+            from flopy4.mf6.utils.crs import cf_grid_mapping_params
+
+            _params, _warning = cf_grid_mapping_params(_crs)
+            if _warning:
+                import warnings
+
+                warnings.warn(_warning)
+            for k, v in _params.items():
+                ds["projection"].attrs[k] = v
 
         return ds
 
@@ -1210,6 +1266,21 @@ def to_xarray(
             ds["layer"].attrs["axis"] = "Z"
             ds["layer"].encoding["_FillValue"] = None
 
+            # z_l1, z_l2, ...: per-layer cell center elevation, always written.
+            # Split per layer because face-indexed variables have no layer
+            # dimension to reference a combined z(layer, nmesh_face) (CF-1.13 5.2).
+            _z = self._coords["z"].values
+            for k in range(self.nlay):
+                _varname = f"z_l{k + 1}"
+                ds = ds.assign({_varname: (["nmesh_face"], _z[k])})
+                if _units is not None:
+                    ds[_varname].attrs["units"] = _units
+                ds[_varname].attrs["standard_name"] = "altitude"
+                ds[_varname].attrs["positive"] = "up"
+                ds[_varname].attrs["long_name"] = f"cell center elevation (layer {k + 1})"
+                ds[_varname].attrs["layer"] = k + 1
+                ds[_varname].encoding["_FillValue"] = None
+
             # mesh container variable
             ds = ds.assign({"mesh": ([], np.int64(1))})
             ds["mesh"].attrs["cf_role"] = "mesh_topology"
@@ -1257,8 +1328,11 @@ def to_xarray(
             ds["mesh_face_y"].attrs["long_name"] = "Northing"
             ds["mesh_face_y"].attrs["bounds"] = "mesh_face_ybnds"
             ds["mesh_face_y"].encoding["_FillValue"] = None
-            ds["mesh_face_xbnds"].encoding["_FillValue"] = FILL_FLOAT64
-            ds["mesh_face_ybnds"].encoding["_FillValue"] = FILL_FLOAT64
+            # No explicit _FillValue attribute: matches MF6 (MeshNCModel.f90), which
+            # pads with NF90_FILL_DOUBLE in the data itself but never declares the
+            # attribute on mesh_face_xbnds/ybnds.
+            ds["mesh_face_xbnds"].encoding["_FillValue"] = None
+            ds["mesh_face_ybnds"].encoding["_FillValue"] = None
 
             # mesh face nodes
             var_d = {
diff --git a/flopy4/mf6/utils/heads_reader.py b/flopy4/mf6/utils/heads_reader.py
index 9c658fb2..41400ec7 100644
--- a/flopy4/mf6/utils/heads_reader.py
+++ b/flopy4/mf6/utils/heads_reader.py
@@ -194,9 +194,13 @@ def _open_hds_netcdf(
     # Time is already CF-encoded datetime64; load eagerly (small coordinate).
     time_values = ds["time"].values
 
-    # --- Conventional CF structured: no "mesh" global attribute ---
+    # --- Conventional CF structured: no "modflow_mesh" global attribute ---
     # Single head(time, z, y, x) variable; x/y are dimension coordinates.
-    if "mesh" not in ds.attrs:
+    # "modflow_mesh" is the current attribute name (renamed from the bare
+    # "mesh" global attr, which collided with UGRID's reserved per-variable
+    # "mesh" attribute in some CF readers); "mesh" is still checked as a
+    # fallback for files written before the rename.
+    if "modflow_mesh" not in ds.attrs and "mesh" not in ds.attrs:
         if grid_type != "STRUCTURED":
             raise ValueError(
                 f"Conventional CF format (no 'mesh' attribute) is only supported "
diff --git a/flopy4/mf6/utils/netcdf_postprocess.py b/flopy4/mf6/utils/netcdf_postprocess.py
index 303963b0..41fbf56e 100644
--- a/flopy4/mf6/utils/netcdf_postprocess.py
+++ b/flopy4/mf6/utils/netcdf_postprocess.py
@@ -1,14 +1,26 @@
 """Post-processing utilities to bring MF6 output NC files into CF/CRS parity
 with flopy4 input NC files.
 
-MF6 writes a ``projection`` variable but omits some attributes required for
-full CF-1.11 compliance and GDAL-based tool support.
-
-- **Mesh output**: missing ``crs_wkt`` and ``grid_mapping_name``.
-- **Structured output**: missing ``wkt``, ``grid_mapping_name``, and the GDAL
-  georeferencing attributes (``GeoTransform`` / ``spatial_ref``) needed for
-  correct placement in QGIS and other GDAL-based tools.  ArcGIS Pro does not
-  require these attributes — it reads ``crs_wkt`` directly from raw MF6 output.
+MF6 6.8.0+ already writes ``crs_wkt``/``grid_mapping_name`` (mesh) and
+``wkt``/``crs_wkt``/``grid_mapping``/``grid_mapping_name`` (structured) --
+this module is then a no-op bridge for those.  Two things it still does:
+
+- **Legacy pre-6.8.0 files**: backfill missing ``wkt``/``crs_wkt``/
+  ``grid_mapping_name``, and fix ``crs_wkt`` holding WKT1 instead of WKT2
+  (a real bug present through 6.7.0, fixed in 6.8.0).
+- **Structured output, any release before CF grid_mapping numeric
+  parameter support**: add the CF-standard numeric grid_mapping
+  parameters (e.g. ``false_easting``, ``scale_factor_at_central_meridian``)
+  that ArcGIS's classic netCDF connector reads directly, for
+  ``transverse_mercator``, ``lambert_conformal_conic`` (2SP only), and
+  ``albers_conical_equal_area``.
+
+This module cannot backfill rotation-positioning for a rotated
+(ANGROT != 0) structured grid on a legacy file: MF6 never writes
+xorigin/yorigin/angrot as retrievable output attributes, and for a
+rotated grid the file's own x/y coordinate values are grid-local with no
+recoverable real-world position. Rotation-positioning is only available
+from a release where MF6 wraps crs_wkt in a derived CRS at write time.
 
 Usage::
 
@@ -33,8 +45,8 @@
 def _apply_mesh_crs_attrs(ds: xr.Dataset) -> xr.Dataset:
     """Add ``crs_wkt`` and ``grid_mapping_name`` to the ``projection`` variable.
 
-    MF6 mesh output already writes ``wkt``; this brings the variable into full
-    CF-1.11 parity with flopy4 input files.
+    No-op on MF6 6.8.0+, which already writes both. Bridges legacy
+    pre-6.8.0 files that only have ``wkt``.
     """
     if "projection" not in ds:
         return ds
@@ -53,7 +65,7 @@ def _apply_mesh_crs_attrs(ds: xr.Dataset) -> xr.Dataset:
     crs = ProjCRS.from_wkt(wkt)
     cf = crs.to_cf()
 
-    # wkt on mesh output is WKT1; crs_wkt must be WKT2 per CF-1.11
+    # wkt on mesh output is WKT1; crs_wkt must be WKT2 per CF-1.13
     ds["projection"].attrs.setdefault("crs_wkt", crs.to_wkt(WktVersion.WKT2_2019))
     gmn = cf.get("grid_mapping_name")
     if gmn:
@@ -63,11 +75,9 @@ def _apply_mesh_crs_attrs(ds: xr.Dataset) -> xr.Dataset:
 
 
 def _apply_structured_crs_attrs(ds: xr.Dataset) -> xr.Dataset:
-    """Add ``wkt``, ``grid_mapping_name``, ``GeoTransform``, and ``spatial_ref``
-    to the ``projection`` variable of an MF6 structured output NC file.
-
-    MF6 structured output already writes ``crs_wkt`` and ``grid_mapping`` on
-    x/y/head; this adds the remaining attrs needed for GDAL-based tools.
+    """Add ``wkt``/``grid_mapping_name`` (legacy bridge), fix ``crs_wkt``
+    (legacy bug), and add the CF-standard numeric grid_mapping parameters
+    (for a release before MF6 wrote them) to the ``projection`` variable.
     """
     if "projection" not in ds:
         return ds
@@ -83,12 +93,26 @@ def _apply_structured_crs_attrs(ds: xr.Dataset) -> xr.Dataset:
 
     from pyproj.enums import WktVersion
 
-    crs = ProjCRS.from_wkt(wkt)
-    cf = crs.to_cf()
+    from flopy4.mf6.utils.crs import cf_grid_mapping_params
 
-    # MF6 structured output writes WKT1 to crs_wkt; overwrite with WKT2 per CF-1.11.
-    # wkt and spatial_ref remain WKT1 for GDAL/legacy-tool compatibility.
-    _wkt1 = crs.to_wkt(WktVersion.WKT1_GDAL)
+    crs = ProjCRS.from_wkt(wkt)
+    # If crs_wkt is already a wrapped DerivedProjectedCRS (a rotated grid
+    # from a release that wraps it at write time), grid_mapping_name/CF
+    # numeric parameters must come from its base CRS -- to_cf() does not
+    # "see through" the wrapper. An ordinary ProjectedCRS also has a
+    # source_crs (its underlying geographic CRS), so only unwrap when
+    # crs itself is actually a DerivedProjectedCRS.
+    base_crs = crs
+    if crs.type_name == "Derived Projected CRS" and crs.source_crs is not None:
+        base_crs = crs.source_crs
+    cf = base_crs.to_cf()
+
+    # Pre-6.8.0 MF6 wrote WKT1 to crs_wkt (bug, fixed in 6.8.0); overwrite
+    # unconditionally -- idempotent on already-correct 6.8.0+ files.
+    # wkt remains WKT1 for legacy-tool compatibility. Derived from
+    # base_crs, not crs: WKT1 has no derived-CRS syntax, so converting an
+    # already-wrapped DerivedProjectedCRS to WKT1 raises a CRSError.
+    _wkt1 = base_crs.to_wkt(WktVersion.WKT1_GDAL)
     _wkt2 = crs.to_wkt(WktVersion.WKT2_2019)
     ds["projection"].attrs["crs_wkt"] = _wkt2
     ds["projection"].attrs.setdefault("wkt", _wkt1)
@@ -96,35 +120,13 @@ def _apply_structured_crs_attrs(ds: xr.Dataset) -> xr.Dataset:
     if gmn:
         ds["projection"].attrs.setdefault("grid_mapping_name", gmn)
 
-    # Derive GeoTransform from x_bnds/y_bnds if available, otherwise from
-    # cell-centre spacing. GDAL reads GeoTransform from the grid_mapping
-    # variable (not global attrs) to set the raster extent.
-    if "x_bnds" in ds and "y_bnds" in ds:
-        xb = ds["x_bnds"].values
-        yb = ds["y_bnds"].values
-        x_left = float(xb[0, 0])
-        x_right = float(xb[-1, 1])
-        y_top = float(yb[0, 1])  # y_bnds[row, 1] = top of row
-        y_bot = float(yb[-1, 0])  # y_bnds[row, 0] = bottom of row
-        ncol = ds.sizes.get("x", xb.shape[0])
-        nrow = ds.sizes.get("y", yb.shape[0])
-        dx_eff = (x_right - x_left) / ncol
-        dy_eff = (y_bot - y_top) / nrow  # negative for north-up
-    elif "x" in ds and "y" in ds:
-        x = ds["x"].values
-        y = ds["y"].values
-        dx = float(x[1] - x[0]) if len(x) > 1 else 1.0
-        dy = float(y[1] - y[0]) if len(y) > 1 else 1.0
-        x_left = float(x[0]) - 0.5 * dx
-        y_top = float(y[0]) - 0.5 * dy
-        dx_eff = dx
-        dy_eff = dy
-    else:
-        return ds
+    params, warning = cf_grid_mapping_params(base_crs)
+    if warning:
+        import warnings
 
-    gt = [x_left, dx_eff, 0.0, y_top, 0.0, dy_eff]
-    ds["projection"].attrs.setdefault("GeoTransform", " ".join(str(v) for v in gt))
-    ds["projection"].attrs.setdefault("spatial_ref", _wkt1)
+        warnings.warn(warning)
+    for k, v in params.items():
+        ds["projection"].attrs.setdefault(k, v)
 
     return ds
 
@@ -133,11 +135,10 @@ def postprocess_mesh_nc(
     path: Union[str, Path],
     out: Union[str, Path, None] = None,
 ) -> Path:
-    """Post-process an MF6 UGRID/mesh output NC file for CF-1.11 compliance.
+    """Post-process an MF6 UGRID/mesh output NC file for CF-1.13 compliance.
 
-    Adds the missing ``crs_wkt`` and ``grid_mapping_name`` attributes to the
-    ``projection`` variable so the file matches the conventions written by
-    flopy4 for input files.
+    Backfills ``crs_wkt`` and ``grid_mapping_name`` on the ``projection``
+    variable for legacy pre-6.8.0 files; a no-op on 6.8.0+.
 
     Parameters
     ----------
@@ -168,14 +169,13 @@ def postprocess_structured_nc(
     path: Union[str, Path],
     out: Union[str, Path, None] = None,
 ) -> Path:
-    """Post-process an MF6 CF-structured output NC file for full CF-1.11 and
-    GDAL compliance.
+    """Post-process an MF6 CF-structured output NC file for CF-1.13/GDAL compliance.
 
-    Adds the missing ``wkt``, ``grid_mapping_name``, ``GeoTransform``, and
-    ``spatial_ref`` attributes to the ``projection`` variable so the file
-    matches the conventions written by flopy4 for input files and is correctly
-    placed by QGIS and other GDAL-based tools.  ArcGIS Pro reads ``crs_wkt``
-    directly from raw MF6 output and does not require this post-processing.
+    Backfills ``wkt``/``grid_mapping_name`` (legacy pre-6.8.0 bridge, no-op
+    on 6.8.0+) and the CF-standard numeric grid_mapping parameters that
+    ArcGIS's classic netCDF connector reads directly (for a release
+    before MF6 wrote them). Does not backfill rotation-positioning for a
+    rotated grid -- see the module docstring.
 
     Parameters
     ----------
@@ -213,7 +213,7 @@ def main():
     parser = argparse.ArgumentParser(
         prog="ncfix",
         description=(
-            "Post-process an MF6 output NetCDF file to add missing CF-1.11 / GDAL "
+            "Post-process an MF6 output NetCDF file to add missing CF-1.13 / GDAL "
             "attributes to the projection variable."
         ),
     )
diff --git a/pyproject.toml b/pyproject.toml
index 9fabf80d..11bb5c8d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -117,6 +117,7 @@ ignore = [
 
 [project.scripts]
 flopy4 = "flopy4.cli:main"
+ncfix = "flopy4.mf6.utils.netcdf_postprocess:main"
 
 [project.entry-points.flopy4]
 plot = "flopy4.singledispatch.plot_int"
diff --git a/test/mf6/test_mf6_netcdf.py b/test/mf6/test_mf6_netcdf.py
index 8ead6e46..fe9e7143 100644
--- a/test/mf6/test_mf6_netcdf.py
+++ b/test/mf6/test_mf6_netcdf.py
@@ -64,7 +64,7 @@ def test_model_nomesh():
 
     assert ds.attrs["modflow_grid"] == "structured"
     assert ds.attrs["modflow_model"] == "gwf6: gwfmodel"
-    assert "mesh" not in ds.attrs
+    assert "modflow_mesh" not in ds.attrs
     assert "welg_0_q" in ds
     assert "welg_0_concentration" in ds
     assert np.allclose(ds["welg_0_q"].values, welg_0_q)
@@ -121,7 +121,7 @@ def test_package_nomesh():
 
     # assert ds.attrs["modflow_grid"] == "structured"
     # assert ds.attrs["modflow_model"] == "gwf6: gwfmodel"
-    assert "mesh" not in ds.attrs
+    assert "modflow_mesh" not in ds.attrs
     assert "welg_0_q" in ds
     assert "welg_0_concentration" in ds
     assert "welg_0_temperature" in ds
@@ -208,7 +208,7 @@ def test_model_mesh():
         "modeltype": "gwf6",
         "modelname": "gwfmodel",
         "gridtype": "structured",
-        "attrs": {"mesh": "layered"},
+        "attrs": {"modflow_mesh": "layered"},
         "packages": packages,
     }
 
@@ -220,7 +220,7 @@ def test_model_mesh():
 
     assert ds.attrs["modflow_grid"] == "structured"
     assert ds.attrs["modflow_model"] == "gwf6: gwfmodel"
-    assert ds.attrs["mesh"] == "layered"
+    assert ds.attrs["modflow_mesh"] == "layered"
     assert "dis_delr" in ds
     assert "dis_delc" in ds
     assert "rcha_0_recharge" in ds
diff --git a/test/mf6/test_mf6_netcdf_cf_conventions.py b/test/mf6/test_mf6_netcdf_cf_conventions.py
index 2cb88fbc..1cdbba8a 100644
--- a/test/mf6/test_mf6_netcdf_cf_conventions.py
+++ b/test/mf6/test_mf6_netcdf_cf_conventions.py
@@ -1,10 +1,10 @@
 """
-Verify CF-1.11 and UGRID-1.0 attribute compliance on flopy4-generated
+Verify CF-1.13 and UGRID-1.0 attribute compliance on flopy4-generated
 NetCDF input datasets for each supported grid/format combination:
 
-  - Structured DIS  (CF-1.11, NetCDFFormat.STRUCTURED)
-  - Layered Mesh DIS  (CF-1.11 + UGRID-1.0, NetCDFFormat.LAYERED_MESH via StructuredGrid)
-  - Layered Mesh DISV (CF-1.11 + UGRID-1.0, NetCDFFormat.LAYERED_MESH via VertexGrid)
+  - Structured DIS  (CF-1.13, NetCDFFormat.STRUCTURED)
+  - Layered Mesh DIS  (CF-1.13 + UGRID-1.0, NetCDFFormat.LAYERED_MESH via StructuredGrid)
+  - Layered Mesh DISV (CF-1.13 + UGRID-1.0, NetCDFFormat.LAYERED_MESH via VertexGrid)
 
 These are unit-level tests that exercise grid.to_xarray() and
 NetCDFParam.to_xarray() directly — no MODFLOW 6 simulation required.
@@ -62,8 +62,8 @@ def vertex_grid():
         [5, 200.0, 0.0],
     ]
     cell2d = [
-        [0, 50.0, 50.0, 0, 1, 4, 3],
-        [1, 150.0, 50.0, 1, 2, 5, 4],
+        [0, 50.0, 50.0, 4, 0, 1, 4, 3],
+        [1, 150.0, 50.0, 4, 1, 2, 5, 4],
     ]
     return VertexGrid(
         nlay=2,
@@ -99,7 +99,7 @@ def test_structured_dis_projection_variable(structured_grid, modeltime):
 
 
 def test_structured_dis_crs_wkt_attrs(structured_grid, modeltime):
-    """crs_wkt must be WKT2 and wkt must be WKT1 on the projection variable (CF-1.11)."""
+    """crs_wkt must be WKT2 and wkt must be WKT1 on the projection variable (CF-1.13)."""
     ds = structured_grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.STRUCTURED)
     proj = ds["projection"]
     assert "crs_wkt" in proj.attrs, "crs_wkt missing from projection variable"
@@ -135,7 +135,7 @@ def test_mesh_dis_topology_variable(structured_grid, modeltime):
 
 
 def test_mesh_dis_projection_crs_attrs(structured_grid, modeltime):
-    """crs_wkt must be WKT2 and wkt must be WKT1 on the projection variable (CF-1.11)."""
+    """crs_wkt must be WKT2 and wkt must be WKT1 on the projection variable (CF-1.13)."""
     ds = structured_grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.LAYERED_MESH)
     assert "projection" in ds
     proj = ds["projection"]
@@ -173,7 +173,7 @@ def test_mesh_dis_data_var_mesh_location_attrs_no_crs():
     var = "npf_k_l1"
     assert ds[var].attrs.get("mesh") == "mesh"
     assert ds[var].attrs.get("location") == "face"
-    assert ds[var].attrs.get("coordinates") == "mesh_face_x mesh_face_y"
+    assert ds[var].attrs.get("coordinates") == "mesh_face_x mesh_face_y z_l1"
 
 
 def test_mesh_disv_topology_variable(vertex_grid, modeltime):
@@ -184,7 +184,7 @@ def test_mesh_disv_topology_variable(vertex_grid, modeltime):
 
 
 def test_mesh_disv_projection_crs_attrs(vertex_grid, modeltime):
-    """crs_wkt must be WKT2 and wkt must be WKT1 on the projection variable (CF-1.11)."""
+    """crs_wkt must be WKT2 and wkt must be WKT1 on the projection variable (CF-1.13)."""
     ds = vertex_grid.to_xarray(modeltime=modeltime)
     assert "projection" in ds
     proj = ds["projection"]
@@ -242,6 +242,52 @@ def test_mesh_disv_layer_coord(vertex_grid, modeltime):
     assert ds["layer"].attrs.get("positive") == "down"
 
 
+def test_structured_dis_z_coord(structured_grid, modeltime):
+    """z must be a (layer, y, x) cell-center elevation auxiliary coordinate.
+
+    structured_grid: top=10.0, thickness=5.0 -> botm=[5.0, 0.0], so layer 1
+    z = (10+5)/2 = 7.5, layer 2 z = (5+0)/2 = 2.5.
+    """
+    ds = structured_grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.STRUCTURED)
+    assert "z" in ds, "z coordinate missing from structured dataset"
+    assert ds["z"].dims == ("layer", "y", "x")
+    assert ds["z"].attrs.get("standard_name") == "altitude"
+    assert ds["z"].attrs.get("positive") == "up"
+    assert ds["z"].attrs.get("long_name") == "cell center elevation"
+    assert np.allclose(ds["z"].isel(layer=0).values, 7.5)
+    assert np.allclose(ds["z"].isel(layer=1).values, 2.5)
+    assert ds["z"].encoding.get("_FillValue") is None
+
+
+def test_mesh_dis_z_coord(structured_grid, modeltime):
+    """z_l1, z_l2, ... must be per-layer cell-center elevation coordinates on nmesh_face."""
+    ds = structured_grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.LAYERED_MESH)
+    for k, expected in ((1, 7.5), (2, 2.5)):
+        var = f"z_l{k}"
+        assert var in ds, f"{var} missing from layered-mesh dataset"
+        assert ds[var].dims == ("nmesh_face",)
+        assert ds[var].attrs.get("standard_name") == "altitude"
+        assert ds[var].attrs.get("positive") == "up"
+        assert ds[var].attrs.get("long_name") == f"cell center elevation (layer {k})"
+        assert ds[var].attrs.get("layer") == k
+        assert np.allclose(ds[var].values, expected)
+        assert ds[var].encoding.get("_FillValue") is None
+
+
+def test_mesh_disv_z_coord(vertex_grid, modeltime):
+    """z_l1, z_l2, ... must be per-layer cell-center elevation coordinates for DISV too."""
+    ds = vertex_grid.to_xarray(modeltime=modeltime)
+    for k in (1, 2):
+        var = f"z_l{k}"
+        assert var in ds, f"{var} missing from layered-mesh DISV dataset"
+        assert ds[var].dims == ("nmesh_face",)
+        assert ds[var].attrs.get("standard_name") == "altitude"
+        assert ds[var].attrs.get("positive") == "up"
+        assert ds[var].attrs.get("long_name") == f"cell center elevation (layer {k})"
+        assert ds[var].attrs.get("layer") == k
+        assert ds[var].encoding.get("_FillValue") is None
+
+
 def test_mesh_dis_face_nodes_fill_in_encoding(structured_grid, modeltime):
     """mesh_face_nodes _FillValue must be in encoding, not attrs (DIS path).
 
@@ -381,7 +427,10 @@ def test_mesh_data_var_cf_attrs_survive_roundtrip(structured_grid, modeltime, tm
 
 
 def test_structured_dis_gdal_geotransform(structured_grid, modeltime, tmp_path):
-    """GDAL reads the correct GeoTransform from the projection variable."""
+    """GDAL derives the correct extent/GetGeoTransform() from crs_wkt plus
+    the projected x/y dimension coordinate arrays alone -- no explicit
+    GeoTransform attribute is written.
+    """
     gdal = pytest.importorskip("osgeo.gdal")
 
     import xarray as xr
@@ -472,8 +521,8 @@ def vertex_grid_no_crs():
         [5, 200.0, 0.0],
     ]
     cell2d = [
-        [0, 50.0, 50.0, 0, 1, 4, 3],
-        [1, 150.0, 50.0, 1, 2, 5, 4],
+        [0, 50.0, 50.0, 4, 0, 1, 4, 3],
+        [1, 150.0, 50.0, 4, 1, 2, 5, 4],
     ]
     return VertexGrid(
         nlay=2,
@@ -504,8 +553,8 @@ def vertex_grid_mixed_poly():
         [4, 2.0, 0.0],
     ]
     cell2d = [
-        [0, 0.667, 1.0, 0, 1, 2],  # triangle: 3 vertices (len-3=3)
-        [1, 1.5, 1.0, 1, 3, 4, 2],  # quad: 4 vertices (len-3=4)
+        [0, 0.667, 1.0, 3, 0, 1, 2],  # triangle: 3 vertices
+        [1, 1.5, 1.0, 4, 1, 3, 4, 2],  # quad: 4 vertices
     ]
     return VertexGrid(nlay=2, ncpl=2, vertices=vertices, cell2d=cell2d)
 
@@ -518,7 +567,7 @@ def _model_ds(grid, flopy4_time, gridtype, mesh=None):
     the grid.to_xarray() level.
     """
     dims = [2, 2, 3, 3] if gridtype == "structured" else [2, 2, 2]
-    attrs = {"mesh": "layered"} if mesh == "layered" else {}
+    attrs = {"modflow_mesh": "layered"} if mesh == "layered" else {}
     nc_model = NetCDFModel.from_dict(
         meta={
             "modeltype": "gwf6",
@@ -558,8 +607,10 @@ def test_structured_global_attrs(structured_grid, flopy4_time):
     ds = _model_ds(structured_grid, flopy4_time, gridtype="structured")
     assert ds.attrs.get("modflow_grid") == "structured"
     assert ds.attrs.get("modflow_model") == "gwf6: testmodel"
-    assert ds.attrs.get("Conventions") == "CF-1.11"
-    assert "mesh" not in ds.attrs, "mesh global attr must be absent for structured format"
+    assert ds.attrs.get("Conventions") == "CF-1.13"
+    assert "modflow_mesh" not in ds.attrs, (
+        "modflow_mesh global attr must be absent for structured format"
+    )
 
 
 def test_mesh_dis_global_attrs(structured_grid, flopy4_time):
@@ -567,16 +618,20 @@ def test_mesh_dis_global_attrs(structured_grid, flopy4_time):
     ds = _model_ds(structured_grid, flopy4_time, gridtype="structured", mesh="layered")
     assert ds.attrs.get("modflow_grid") == "structured"
     assert ds.attrs.get("modflow_model") == "gwf6: testmodel"
-    assert ds.attrs.get("Conventions") == "CF-1.11 UGRID-1.0"
-    assert ds.attrs.get("mesh") == "layered", "mesh global attr must be 'layered' (lowercase)"
+    assert ds.attrs.get("Conventions") == "CF-1.13 UGRID-1.0"
+    assert ds.attrs.get("modflow_mesh") == "layered", (
+        "modflow_mesh global attr must be 'layered' (lowercase)"
+    )
 
 
 def test_mesh_disv_global_attrs(vertex_grid, flopy4_time):
     """modflow_grid must be 'vertex' and Conventions must include UGRID for DISV layered mesh."""
     ds = _model_ds(vertex_grid, flopy4_time, gridtype="vertex", mesh="layered")
     assert ds.attrs.get("modflow_grid") == "vertex"
-    assert ds.attrs.get("Conventions") == "CF-1.11 UGRID-1.0"
-    assert ds.attrs.get("mesh") == "layered", "mesh global attr must be 'layered' (lowercase)"
+    assert ds.attrs.get("Conventions") == "CF-1.13 UGRID-1.0"
+    assert ds.attrs.get("modflow_mesh") == "layered", (
+        "modflow_mesh global attr must be 'layered' (lowercase)"
+    )
 
 
 # time coordinate attrs
@@ -697,21 +752,191 @@ def test_structured_bnds_fill_value_suppressed(structured_grid, modeltime):
 
 
 # projection variable: GeoTransform and spatial_ref
-def test_structured_projection_gdal_attrs(structured_grid, modeltime):
-    """GeoTransform and spatial_ref must be present on the structured projection variable.
-
-    GeoTransform is required for GDAL-based tools (QGIS, ArcGIS) to place the
-    raster correctly. spatial_ref is GDAL's companion WKT1 attribute.
+def test_structured_projection_no_geotransform(structured_grid, modeltime):
+    """GeoTransform and spatial_ref are never written; CF-standard numeric
+    grid_mapping parameters are written instead on the structured
+    projection variable for GDAL-based tools (QGIS, ArcGIS) to position
+    the projected raster correctly.
     """
     ds = structured_grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.STRUCTURED)
     proj = ds["projection"]
-    assert "GeoTransform" in proj.attrs, "GeoTransform missing from structured projection variable"
-    gt_parts = str(proj.attrs["GeoTransform"]).split()
-    assert len(gt_parts) == 6, (
-        f"GeoTransform must have 6 values, got: {proj.attrs['GeoTransform']!r}"
+    assert "GeoTransform" not in proj.attrs, "GeoTransform must never be written"
+    assert "spatial_ref" not in proj.attrs, "spatial_ref must never be written"
+
+    assert proj.attrs["grid_mapping_name"] == "transverse_mercator"
+    assert proj.attrs["longitude_of_central_meridian"] == pytest.approx(-117.0)
+    assert proj.attrs["latitude_of_projection_origin"] == pytest.approx(0.0)
+    assert proj.attrs["scale_factor_at_central_meridian"] == pytest.approx(0.9996)
+    assert proj.attrs["false_easting"] == pytest.approx(500_000.0)
+    assert proj.attrs["false_northing"] == pytest.approx(0.0)
+    assert proj.attrs["semi_major_axis"] == pytest.approx(6378137.0)
+    assert proj.attrs["inverse_flattening"] == pytest.approx(298.257223563)
+
+
+def test_structured_projection_crs_wkt_rotated(modeltime):
+    """A rotated grid's crs_wkt must be wrapped in a DerivedProjectedCRS
+    encoding the rotation via EPSG:9624 (affine parametric
+    transformation), with A0..B2 matching the world->local rotation
+    formula exactly; x/y must retain axis/standard_name since they are
+    correct native coordinates in the derived CRS.
+    """
+    import re
+
+    xorigin, yorigin, angrot = 573000.0, 4100000.0, 25.0
+    grid = StructuredGrid.uniform(
+        nlay=1,
+        nrow=3,
+        ncol=4,
+        delr=100.0,
+        delc=150.0,
+        top=10.0,
+        thickness=10.0,
+        xoff=xorigin,
+        yoff=yorigin,
+        angrot=angrot,
+        crs="EPSG:26911",
+    )
+    ds = grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.STRUCTURED)
+    proj = ds["projection"]
+    crs_wkt = proj.attrs["crs_wkt"]
+
+    assert crs_wkt.startswith("DERIVEDPROJCRS["), (
+        f"crs_wkt must be wrapped in a DerivedProjectedCRS: {crs_wkt[:80]}"
+    )
+    assert "Affine parametric transformation" in crs_wkt
+    assert 'ID["EPSG",9624]' in crs_wkt
+
+    params = {}
+    for pname in ("A0", "A1", "A2", "B0", "B1", "B2"):
+        m = re.search(rf'PARAMETER\["{pname}",([^,\]]+)', crs_wkt)
+        assert m, f"{pname} parameter not found in crs_wkt"
+        params[pname] = float(m.group(1))
+
+    ang = np.radians(angrot)
+    expected = {
+        "A0": -(xorigin * np.cos(ang) + yorigin * np.sin(ang)),
+        "A1": np.cos(ang),
+        "A2": np.sin(ang),
+        "B0": xorigin * np.sin(ang) - yorigin * np.cos(ang),
+        "B1": -np.sin(ang),
+        "B2": np.cos(ang),
+    }
+    for pname, value in expected.items():
+        assert params[pname] == pytest.approx(value, rel=1e-10), (
+            f"{pname}: {params[pname]} != expected {value}"
+        )
+
+    assert "axis" in ds["x"].attrs, "x axis attr must be present when the wrap succeeds"
+    assert "axis" in ds["y"].attrs, "y axis attr must be present when the wrap succeeds"
+    assert "GeoTransform" not in proj.attrs
+    assert "spatial_ref" not in proj.attrs
+
+    # The A0..B2 parameters alone don't prove x/y are written on the
+    # coordinate convention the wrap actually assumes (zero-based local,
+    # not xorigin/yorigin-offset) -- transform the *stored* x/y values
+    # themselves through the wrapped CRS and compare against the true
+    # world position computed independently.
+    pyproj = pytest.importorskip("pyproj")
+    derived_crs = pyproj.CRS.from_wkt(crs_wkt)
+    base_crs = pyproj.CRS.from_epsg(26911)
+    transformer = pyproj.Transformer.from_crs(derived_crs, base_crs, always_xy=True)
+    for lx, ly in zip(ds["x"].values, ds["y"].values):
+        world_x, world_y = transformer.transform(float(lx), float(ly))
+        expected_x = xorigin + lx * np.cos(ang) - ly * np.sin(ang)
+        expected_y = yorigin + lx * np.sin(ang) + ly * np.cos(ang)
+        assert np.allclose((world_x, world_y), (expected_x, expected_y), rtol=1e-9), (
+            f"stored (x,y)=({lx},{ly}) resolves to ({world_x},{world_y}), "
+            f"expected ({expected_x},{expected_y})"
+        )
+
+
+def test_structured_projection_lcc_2sp(modeltime):
+    """lambert_conformal_conic (2SP) CRS: grid_mapping_name plus the
+    standard_parallel array and the 6 shared CF numeric parameters are
+    all written."""
+    grid = StructuredGrid.uniform(
+        nlay=1,
+        nrow=3,
+        ncol=4,
+        delr=100.0,
+        delc=150.0,
+        top=10.0,
+        thickness=10.0,
+        xoff=1_005_000.0,
+        yoff=1_005_000.0,
+        crs="EPSG:3081",  # NAD83 / Texas State Mapping System
+    )
+    ds = grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.STRUCTURED)
+    proj = ds["projection"]
+    assert proj.attrs["grid_mapping_name"] == "lambert_conformal_conic"
+    assert np.allclose(proj.attrs["standard_parallel"], [27.4166666666667, 34.9166666666667])
+    assert proj.attrs["longitude_of_central_meridian"] == pytest.approx(-100.0)
+    assert proj.attrs["latitude_of_projection_origin"] == pytest.approx(31.1666666666667)
+    assert proj.attrs["false_easting"] == pytest.approx(1_000_000.0)
+    assert proj.attrs["false_northing"] == pytest.approx(1_000_000.0)
+    assert proj.attrs["semi_major_axis"] == pytest.approx(6378137.0)
+    assert proj.attrs["inverse_flattening"] == pytest.approx(298.257222101)
+
+
+def test_structured_projection_albers(modeltime):
+    """albers_conical_equal_area CRS: grid_mapping_name plus the
+    standard_parallel array and the 6 shared CF numeric parameters are
+    all written."""
+    grid = StructuredGrid.uniform(
+        nlay=1,
+        nrow=3,
+        ncol=4,
+        delr=100.0,
+        delc=150.0,
+        top=10.0,
+        thickness=10.0,
+        xoff=100_000.0,
+        yoff=100_000.0,
+        crs="EPSG:5070",  # NAD83 / Conus Albers
+    )
+    ds = grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.STRUCTURED)
+    proj = ds["projection"]
+    assert proj.attrs["grid_mapping_name"] == "albers_conical_equal_area"
+    assert np.allclose(proj.attrs["standard_parallel"], [29.5, 45.5])
+    assert proj.attrs["longitude_of_central_meridian"] == pytest.approx(-96.0)
+    assert proj.attrs["latitude_of_projection_origin"] == pytest.approx(23.0)
+    assert proj.attrs["false_easting"] == pytest.approx(0.0)
+    assert proj.attrs["false_northing"] == pytest.approx(0.0)
+    assert proj.attrs["semi_major_axis"] == pytest.approx(6378137.0)
+    assert proj.attrs["inverse_flattening"] == pytest.approx(298.257222101)
+
+
+def test_structured_projection_lcc_1sp_warns(modeltime):
+    """lambert_conformal_conic (1SP, scale_factor-based) CRS: CF has no
+    scale_factor attribute for this method, so numeric grid_mapping
+    parameter extraction must warn and write nothing, while still
+    writing grid_mapping_name from the CRS's projection method alone."""
+    grid = StructuredGrid.uniform(
+        nlay=1,
+        nrow=3,
+        ncol=4,
+        delr=100.0,
+        delc=150.0,
+        top=10.0,
+        thickness=10.0,
+        xoff=600_100.0,
+        yoff=600_100.0,
+        crs="EPSG:2062",  # Madrid 1870 (Madrid) / Spain LCC
     )
-    assert "spatial_ref" in proj.attrs, "spatial_ref missing from structured projection variable"
-    assert proj.attrs["spatial_ref"].startswith("PROJCS["), "spatial_ref must be WKT1"
+    with pytest.warns(UserWarning, match="1SP"):
+        ds = grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.STRUCTURED)
+    proj = ds["projection"]
+    assert proj.attrs["grid_mapping_name"] == "lambert_conformal_conic"
+    for attr in (
+        "standard_parallel",
+        "longitude_of_central_meridian",
+        "latitude_of_projection_origin",
+        "false_easting",
+        "false_northing",
+        "semi_major_axis",
+        "inverse_flattening",
+    ):
+        assert attr not in proj.attrs, f"{attr} must not be written for an unsupported 1SP variant"
 
 
 def test_mesh_dis_no_gdal_attrs(structured_grid, modeltime):
@@ -830,18 +1055,17 @@ def test_mesh_dis_face_coord_attrs(structured_grid, modeltime):
     assert ds["mesh_face_y"].attrs.get("grid_mapping") == "projection"
 
 
-def test_mesh_dis_face_bounds_fill_in_encoding(structured_grid, modeltime):
-    """mesh_face_xbnds/ybnds _FillValue must be FILL_FLOAT64 in encoding (DIS path).
+def test_mesh_dis_face_bounds_no_fill_attr(structured_grid, modeltime):
+    """mesh_face_xbnds/ybnds must not declare a _FillValue attribute (DIS path).
 
-    These are padded float arrays where padding slots must use NF90_FILL_DOUBLE
-    so CF validators and masking tools correctly identify them. xarray's default
-    for float64 is NaN which is not a valid CF fill value.
+    Matches MF6 (MeshNCModel.f90): padding slots hold NF90_FILL_DOUBLE in the
+    data itself, but the attribute is never declared on these variables.
     """
     ds = structured_grid.to_xarray(modeltime=modeltime, netcdf_format=NetCDFFormat.LAYERED_MESH)
     assert "_FillValue" not in ds["mesh_face_xbnds"].attrs
     assert "_FillValue" not in ds["mesh_face_ybnds"].attrs
-    assert ds["mesh_face_xbnds"].encoding.get("_FillValue") == FILL_FLOAT64
-    assert ds["mesh_face_ybnds"].encoding.get("_FillValue") == FILL_FLOAT64
+    assert ds["mesh_face_xbnds"].encoding.get("_FillValue") is None
+    assert ds["mesh_face_ybnds"].encoding.get("_FillValue") is None
 
 
 def test_mesh_disv_node_coord_attrs(vertex_grid, modeltime):
@@ -867,13 +1091,41 @@ def test_mesh_disv_face_coord_attrs(vertex_grid, modeltime):
     assert ds["mesh_face_y"].attrs.get("grid_mapping") == "projection"
 
 
-def test_mesh_disv_face_bounds_fill_in_encoding(vertex_grid, modeltime):
-    """mesh_face_xbnds/ybnds _FillValue must be FILL_FLOAT64 in encoding (DISV path)."""
+def test_mesh_disv_face_bounds_no_fill_attr(vertex_grid, modeltime):
+    """mesh_face_xbnds/ybnds must not declare a _FillValue attribute (DISV path)."""
     ds = vertex_grid.to_xarray(modeltime=modeltime)
     assert "_FillValue" not in ds["mesh_face_xbnds"].attrs
     assert "_FillValue" not in ds["mesh_face_ybnds"].attrs
-    assert ds["mesh_face_xbnds"].encoding.get("_FillValue") == FILL_FLOAT64
-    assert ds["mesh_face_ybnds"].encoding.get("_FillValue") == FILL_FLOAT64
+    assert ds["mesh_face_xbnds"].encoding.get("_FillValue") is None
+    assert ds["mesh_face_ybnds"].encoding.get("_FillValue") is None
+
+
+def test_mesh_disv_face_bounds_padding_value(modeltime):
+    """Padding slots (cells with fewer vertices than max) must hold FILL_FLOAT64
+    in the actual data, matching MF6's NF90_FILL_DOUBLE -- even though neither
+    system declares a _FillValue attribute for it (see the no_fill_attr tests).
+    """
+    # 1 quad (4 verts) + 1 triangle (3 verts): exercises real padding.
+    vertices = [
+        [0, 0.0, 100.0],
+        [1, 100.0, 100.0],
+        [2, 200.0, 100.0],
+        [3, 0.0, 0.0],
+        [4, 100.0, 0.0],
+        [5, 200.0, 0.0],
+    ]
+    cell2d = [
+        [0, 50.0, 50.0, 4, 0, 1, 4, 3],
+        [1, 150.0, 66.7, 3, 1, 2, 5],
+    ]
+    grid = VertexGrid(nlay=1, ncpl=2, vertices=vertices, cell2d=cell2d, crs=CRS)
+    ds = grid.to_xarray(modeltime=modeltime)
+    assert ds.sizes["max_nmesh_face_nodes"] == 4
+    # triangle (face index 1) has one padding slot
+    assert ds["mesh_face_xbnds"].values[1, -1] == FILL_FLOAT64
+    assert ds["mesh_face_ybnds"].values[1, -1] == FILL_FLOAT64
+    # quad (face index 0) is fully populated, no padding
+    assert not np.any(ds["mesh_face_xbnds"].values[0] == FILL_FLOAT64)
 
 
 def test_mesh_disv_face_bounds_padding_data_values(vertex_grid_mixed_poly, modeltime):
@@ -976,6 +1228,71 @@ def test_mesh_dis_param_long_name(structured_grid):
     assert ds["npf_k_l1"].attrs.get("long_name") == "hydraulic conductivity (L/T) layer 1"
 
 
+def test_from_model_int_param_dtype():
+    """NetCDFModel.from_model() (the real end-to-end path, not the lower-level
+    NetCDFParam helper above) must preserve integer dtype for integer GRIDDATA
+    params, not silently cast them to float64.
+    """
+    from flopy4.mf6.gwf import Dis, Gwf, Ic, Npf
+    from flopy4.mf6.utils.time import Time
+
+    dis = Dis(nlay=1, nrow=1, ncol=3, delr=1.0, delc=1.0, top=1.0, botm=0.0)
+    gwf = Gwf(dis=dis, name="m")
+    Ic(parent=gwf, strt=1.0)
+    Npf(parent=gwf, icelltype=1, k=1.0)
+    time = Time(perlen=[1.0], nstp=[1], tsmult=[1.0])
+
+    nc_model = NetCDFModel.from_model(gwf, netcdf_format=NetCDFFormat.STRUCTURED, time=time)
+    ds = nc_model.to_xarray()
+    assert np.issubdtype(ds["npf_icelltype"].dtype, np.integer), (
+        f"icelltype must be integer-typed, got {ds['npf_icelltype'].dtype}"
+    )
+    assert np.issubdtype(ds["npf_k"].dtype, np.floating), (
+        f"k must remain float-typed, got {ds['npf_k'].dtype}"
+    )
+
+
+def test_from_model_structured_layer_coordinates_z():
+    """Structured per-layer (layer, y, x) fields must link to the z auxiliary
+    coordinate via a `coordinates="z"` attribute, matching MF6's own
+    DisNCStructured.f90 output; a non-layered (y, x) field like dis_top must
+    not, since it has no per-layer elevation to link to.
+    """
+    from flopy4.mf6.gwf import Dis, Gwf, Ic, Npf
+    from flopy4.mf6.utils.time import Time
+
+    dis = Dis(nlay=2, nrow=1, ncol=3, delr=1.0, delc=1.0, top=10.0, botm=[5.0, 0.0])
+    gwf = Gwf(dis=dis, name="rch")
+    Ic(parent=gwf, strt=1.0)
+    Npf(parent=gwf, icelltype=1, k=1.0)
+    time = Time(perlen=[1.0], nstp=[1], tsmult=[1.0])
+
+    nc_model = NetCDFModel.from_model(gwf, netcdf_format=NetCDFFormat.STRUCTURED, time=time)
+    ds = nc_model.to_xarray()
+    assert ds["npf_k"].attrs.get("coordinates") == "z"
+    assert ds["npf_icelltype"].attrs.get("coordinates") == "z"
+    assert ds["dis_botm"].attrs.get("coordinates") == "z"
+    assert "coordinates" not in ds["dis_top"].attrs
+
+
+def test_from_model_title():
+    """title must match MF6's own model-type-specific convention (NCModel.f90),
+    preserving the model name's original casing -- not blanket-lowercased.
+    """
+    from flopy4.mf6.gwf import Dis, Gwf, Ic, Npf
+    from flopy4.mf6.utils.time import Time
+
+    dis = Dis(nlay=1, nrow=1, ncol=3, delr=1.0, delc=1.0, top=1.0, botm=0.0)
+    gwf = Gwf(dis=dis, name="rch")
+    Ic(parent=gwf, strt=1.0)
+    Npf(parent=gwf, icelltype=1, k=1.0)
+    time = Time(perlen=[1.0], nstp=[1], tsmult=[1.0])
+
+    nc_model = NetCDFModel.from_model(gwf, netcdf_format=NetCDFFormat.STRUCTURED, time=time)
+    ds = nc_model.to_xarray()
+    assert ds.attrs.get("title") == "RCH hydraulic head array input"
+
+
 # modflow_input attribute format
 def test_structured_param_modflow_input(structured_grid):
     """modflow_input must follow 'modelname/packagename/paramname' format (all lowercase)."""
diff --git a/test/mf6/test_mf6_netcdf_postprocess.py b/test/mf6/test_mf6_netcdf_postprocess.py
index 19315cf4..6f0210ef 100644
--- a/test/mf6/test_mf6_netcdf_postprocess.py
+++ b/test/mf6/test_mf6_netcdf_postprocess.py
@@ -1,5 +1,8 @@
 """Tests for flopy4.mf6.utils.netcdf_postprocess."""
 
+import subprocess
+import sys
+
 import numpy as np
 import pytest
 import xarray as xr
@@ -149,25 +152,15 @@ def test_postprocess_structured_nc_adds_attrs(tmp_path):
     assert p.attrs["wkt"].startswith("PROJCS["), "wkt must be WKT1"
     assert "grid_mapping_name" in p.attrs
     assert p.attrs["grid_mapping_name"] == "transverse_mercator"
-    assert "GeoTransform" in p.attrs
-    assert "spatial_ref" in p.attrs
-    assert p.attrs["spatial_ref"].startswith("PROJCS["), "spatial_ref must be WKT1 for GDAL"
-    ds.close()
-
-
-def test_postprocess_structured_nc_geotransform_values(tmp_path):
-    src = tmp_path / "out.nc"
-    _write_minimal_structured_nc(src)
-    postprocess_structured_nc(src)
-
-    ds = xr.open_dataset(src, mask_and_scale=False)
-    gt = [float(v) for v in ds["projection"].attrs["GeoTransform"].split()]
-    # x_left  = x_bnds[0,0] = 250,  x_right = x_bnds[-1,1] = 2250
-    # y_top   = y_bnds[0,1] = 3500, y_bot   = y_bnds[-1,0] = 500
-    assert gt[0] == pytest.approx(250.0)  # x origin (left edge)
-    assert gt[1] == pytest.approx(500.0)  # dx = (2250-250)/4
-    assert gt[3] == pytest.approx(3500.0)  # y origin (top edge)
-    assert gt[5] == pytest.approx(-1000.0)  # dy = (500-3500)/3
+    assert "GeoTransform" not in p.attrs
+    assert "spatial_ref" not in p.attrs
+    assert p.attrs["longitude_of_central_meridian"] == pytest.approx(-117.0)
+    assert p.attrs["latitude_of_projection_origin"] == pytest.approx(0.0)
+    assert p.attrs["scale_factor_at_central_meridian"] == pytest.approx(0.9996)
+    assert p.attrs["false_easting"] == pytest.approx(500_000.0)
+    assert p.attrs["false_northing"] == pytest.approx(0.0)
+    assert p.attrs["semi_major_axis"] == pytest.approx(6378137.0)
+    assert p.attrs["inverse_flattening"] == pytest.approx(298.257222101)
     ds.close()
 
 
@@ -179,7 +172,7 @@ def test_postprocess_structured_nc_idempotent(tmp_path):
 
     ds = xr.open_dataset(src, mask_and_scale=False)
     assert "wkt" in ds["projection"].attrs
-    assert "GeoTransform" in ds["projection"].attrs
+    assert ds["projection"].attrs["false_easting"] == pytest.approx(500_000.0)
     ds.close()
 
 
@@ -195,8 +188,9 @@ def test_postprocess_structured_nc_no_projection(tmp_path):
     ds2.close()
 
 
-def test_postprocess_structured_nc_fallback_no_bnds(tmp_path):
-    """GeoTransform is derived from cell-centre spacing when x_bnds/y_bnds are absent."""
+def test_postprocess_structured_nc_no_bnds(tmp_path):
+    """CF numeric grid_mapping parameters only need crs_wkt -- they are
+    added whether or not x_bnds/y_bnds are present."""
     x = np.array([100.0, 200.0, 300.0])
     y = np.array([600.0, 500.0, 400.0])
     src = tmp_path / "no_bnds.nc"
@@ -212,10 +206,89 @@ def test_postprocess_structured_nc_fallback_no_bnds(tmp_path):
 
     postprocess_structured_nc(src)
     ds2 = xr.open_dataset(src, mask_and_scale=False)
-    gt = [float(v) for v in ds2["projection"].attrs["GeoTransform"].split()]
-    # dx=100, dy=-100; x_left = x[0] - 0.5*dx = 50; y_top = y[0] - 0.5*dy = 650
-    assert gt[0] == pytest.approx(50.0)  # x origin (left edge)
-    assert gt[1] == pytest.approx(100.0)  # dx
-    assert gt[3] == pytest.approx(650.0)  # y origin (top edge)
-    assert gt[5] == pytest.approx(-100.0)  # dy
+    assert ds2["projection"].attrs["grid_mapping_name"] == "transverse_mercator"
+    assert ds2["projection"].attrs["false_easting"] == pytest.approx(500_000.0)
     ds2.close()
+
+
+def test_postprocess_structured_nc_already_wrapped_rotated(tmp_path):
+    """A rotated grid's crs_wkt may already be a wrapped DerivedProjectedCRS
+    (current-release MF6 output) -- postprocessing it must be a safe no-op,
+    not crash. WKT1 has no derived-CRS syntax, so wkt must come from the
+    unwrapped base CRS, not from converting the wrapped CRS directly."""
+    pyproj = pytest.importorskip("pyproj")
+    from flopy4.mf6.utils.crs import wrap_rotated_crs
+
+    base = pyproj.CRS.from_epsg(26918)
+    wrapped = wrap_rotated_crs(base, 500000.0, 4500000.0, 15.0)
+    wrapped_wkt = wrapped.to_wkt(pyproj.enums.WktVersion.WKT2_2019)
+    assert wrapped_wkt.startswith("DERIVEDPROJCRS[")
+
+    src = tmp_path / "rotated.nc"
+    ds = xr.Dataset(
+        {"projection": xr.Variable([], np.int32(1), attrs={"crs_wkt": wrapped_wkt})},
+        coords={
+            "x": xr.Variable(["x"], np.array([50.0, 150.0, 250.0])),
+            "y": xr.Variable(["y"], np.array([100.0, 0.0])),
+        },
+    )
+    ds.to_netcdf(src)
+    ds.close()
+
+    postprocess_structured_nc(src)
+    ds2 = xr.open_dataset(src, mask_and_scale=False)
+    p = ds2["projection"]
+    assert p.attrs["crs_wkt"].startswith("DERIVEDPROJCRS["), "wrap must be preserved"
+    assert p.attrs["wkt"].startswith("PROJCS["), "wkt must be the unwrapped base CRS"
+    assert p.attrs["grid_mapping_name"] == "transverse_mercator"
+    assert p.attrs["false_easting"] == pytest.approx(500_000.0)
+    ds2.close()
+
+
+def _run_ncfix(*args):
+    """Invoke ncfix via the module entry point so it works before pip install."""
+    return subprocess.run(
+        [sys.executable, "-m", "flopy4.mf6.utils.netcdf_postprocess", *args],
+        capture_output=True,
+        text=True,
+    )
+
+
+def test_ncfix_cli_structured(tmp_path):
+    src = tmp_path / "structured.nc"
+    _write_minimal_structured_nc(src)
+    dst = tmp_path / "fixed.nc"
+    result = _run_ncfix(str(src), "-o", str(dst), "--mode", "structured")
+    assert result.returncode == 0, result.stderr
+
+    ds = xr.open_dataset(dst, mask_and_scale=False)
+    assert "GeoTransform" not in ds["projection"].attrs
+    assert "crs_wkt" in ds["projection"].attrs
+    assert ds["projection"].attrs["false_easting"] == 500_000.0
+    ds.close()
+
+
+def test_ncfix_cli_mesh_explicit_mode(tmp_path):
+    src = tmp_path / "out.nc"
+    _write_minimal_mesh_nc(src)
+    result = _run_ncfix(str(src), "--mode", "mesh", "--verbose")
+    assert result.returncode == 0, result.stderr
+    assert "wrote:" in result.stdout
+
+    ds = xr.open_dataset(src, mask_and_scale=False)
+    assert "crs_wkt" in ds["projection"].attrs
+    ds.close()
+
+
+def test_ncfix_cli_auto_detect_structured(tmp_path):
+    """Auto-detect falls back to structured when no UGRID markers are present."""
+    src = tmp_path / "out.nc"
+    _write_minimal_structured_nc(src)
+    result = _run_ncfix(str(src), "--verbose")
+    assert result.returncode == 0, result.stderr
+    assert "wrote:" in result.stdout
+
+    ds = xr.open_dataset(src, mask_and_scale=False)
+    assert "GeoTransform" not in ds["projection"].attrs
+    assert ds["projection"].attrs["false_easting"] == 500_000.0
+    ds.close()
diff --git a/test/mf6/test_mf6_output_readers.py b/test/mf6/test_mf6_output_readers.py
index 5d1598d9..50ade517 100644
--- a/test/mf6/test_mf6_output_readers.py
+++ b/test/mf6/test_mf6_output_readers.py
@@ -322,8 +322,11 @@ def test_open_hds_disv_netcdf_ugrid_layered(function_tmpdir, disv_model_output):
     _open_hds_netcdf — CF-UGRID layered format for DISV.
 
     Writes a synthetic NetCDF file with ``head_l1``, ``head_l2`` variables
-    on ``(time, nmesh_face)`` dimensions and a ``mesh`` global attribute,
-    matching what MODFLOW 6 writes for DISV grids in layered mesh mode.
+    on ``(time, nmesh_face)`` dimensions and a ``modflow_mesh`` global
+    attribute, matching what MODFLOW 6 currently writes for DISV grids in
+    layered mesh mode (renamed from the bare ``mesh`` attribute, which
+    collided with UGRID's reserved per-variable ``mesh`` attribute in some
+    CF readers).
     """
     from pathlib import Path
 
@@ -340,7 +343,7 @@ def test_open_hds_disv_netcdf_ugrid_layered(function_tmpdir, disv_model_output):
     ds = xr.Dataset(
         data_vars,
         coords={"time": np.array([1.0, 2.0])},
-        attrs={"mesh": 1, "modflow_grid": "VERTEX"},
+        attrs={"modflow_mesh": 1, "modflow_grid": "VERTEX"},
     )
     nc_path = Path(function_tmpdir) / "head_ugrid.nc"
     ds.to_netcdf(nc_path)
@@ -358,6 +361,40 @@ def test_open_hds_disv_netcdf_ugrid_layered(function_tmpdir, disv_model_output):
         np.testing.assert_allclose(head.values[:, k, :], layer_data[k])
 
 
+def test_open_hds_disv_netcdf_ugrid_layered_legacy_mesh_attr(function_tmpdir, disv_model_output):
+    """
+    _open_hds_netcdf still recognizes the legacy bare ``mesh`` global
+    attribute (pre-rename) as a fallback, for files written before MF6
+    renamed it to ``modflow_mesh``.
+    """
+    from pathlib import Path
+
+    paths = disv_model_output
+    nlay, ncpl = paths["nlay"], paths["ncpl"]
+    ntime = 2
+
+    rng = np.random.default_rng(1)
+    layer_data = [rng.random((ntime, ncpl)) for _ in range(nlay)]
+    data_vars = {
+        f"head_l{k + 1}": xr.DataArray(layer_data[k], dims=("time", "nmesh_face"))
+        for k in range(nlay)
+    }
+    ds = xr.Dataset(
+        data_vars,
+        coords={"time": np.array([1.0, 2.0])},
+        attrs={"mesh": 1, "modflow_grid": "VERTEX"},
+    )
+    nc_path = Path(function_tmpdir) / "head_ugrid_legacy.nc"
+    ds.to_netcdf(nc_path)
+
+    head = open_hds(nc_path, paths["grb"])
+
+    assert isinstance(head, xu.UgridDataArray)
+    assert head.shape[0] == ntime
+    assert head.shape[1] == nlay
+    assert head.shape[2] == ncpl
+
+
 def test_open_hds_netcdf_grid_type_mismatch(function_tmpdir, dis_model_output):
     """
     _open_hds_netcdf raises ValueError when the NetCDF modflow_grid attribute