diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 049ea850..e1ee40c6 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -105,6 +105,7 @@ Upcoming Version
 * ``linopy.merge`` with ``join="outer"``, ``"left"`` or ``"right"`` no longer raises an xarray ``MergeError`` when only one operand carries an auxiliary coordinate on a joined dimension. The constant was reindexed with a fill value of ``0`` that was also written into the auxiliary coordinate, so it no longer matched the coefficients' copy. The coordinate now follows its rows onto the joined grid and is NaN where no operand provides it, under both semantics. The same holds for ``.add`` / ``.sub`` / ``.mul`` / ``.div`` with a constant and a reindexing ``join=``.
 * A frozen constraint caches the label-to-position mapping of its matrix columns by weak reference and only rebuilds it when the constraint or the set of variables changes. While a persistent snapshot holds the arrays, repeated matrix assembly on an unchanged model returns the same objects, so the snapshot diff can again skip the comparison of untouched frozen constraints by object identity; one-off exports such as ``to_file`` retain no extra memory. (`#933 <https://github.com/PyPSA/linopy/issues/933>`__)
 * ``Model.copy`` (and the ``copy.copy``/``copy.deepcopy`` protocols) keeps frozen constraints as ``CSRConstraint`` instead of rebuilding them as dense ``Constraint``; ``deep=True`` copies their arrays, and duals are copied only with ``include_solution=True``. A CSR-backed objective is now deep-copied with ``deep=True`` as well. The auxiliary coordinates of a CSR grid are stored read-only, so a copy can share them safely. (`#981 <https://github.com/PyPSA/linopy/issues/981>`__)
+* ``Model.to_netcdf``/``linopy.read_netcdf`` now keep the auxiliary (non-index) coordinates of dense variables, constraints, expressions and ``Model.parameters``, such as the levels on a flat dimension that v1 recommends instead of a ``MultiIndex``, and scalar coordinates left by a scalar selection. ``read_netcdf`` dropped every non-index coordinate except ``MultiIndex`` levels, to keep out the scalar coordinates of other containers that xarray attaches on selection; it now drops only those. Frozen constraints already kept theirs.
 * ``Solver.close()`` no longer leaves dangling native handles behind. The solver model is now dropped before the environment that owns it, instead of after. And the COPT and MindOpt file interfaces no longer hand back a model they already disposed: after a file-based COPT or MindOpt solve, ``model.solver_model`` is ``None`` rather than a handle into freed memory. (`#899 <https://github.com/PyPSA/linopy/pull/899>`__)
 
 **Breaking Changes**
diff --git a/linopy/io.py b/linopy/io.py
index 80568f18..a5055d3f 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -1230,16 +1230,24 @@ def parse_multiindex_attr(value: str | Iterable[str]) -> list[str]:
             return [str(n) for n in json.loads(value)]
         return [str(n) for n in value]
 
+    # Selecting a container's data variables attaches the scalar coords of all
+    # containers, so each is mapped to its owner: the longest matching prefix,
+    # as names may contain dashes. Coords with dimensions only attach to the
+    # container whose dimensions they share.
+    prefixes = {str(k).rsplit("-", 1)[0] for k in ds}
+    scalar_owner = {
+        c: max(
+            (p for p in prefixes if str(c).startswith(f"{p}-")), key=len, default=None
+        )
+        for c in ds.coords
+        if ds[c].ndim == 0
+    }
+
     def get_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset:
         ds = ds[[k for k in ds if has_prefix(str(k), prefix)]]
-        multiindexes = []
-        for dim in ds.dims:
-            attr = ds.attrs.get(f"{dim}_multiindex")
-            if attr is None:
-                continue
-            for name in parse_multiindex_attr(attr):
-                multiindexes.append(prefix + "-" + name)
-        ds = ds.drop_vars(set(ds.coords) - set(ds.dims) - set(multiindexes))
+        ds = ds.drop_vars(
+            [c for c in ds.coords if scalar_owner.get(c, prefix) != prefix]
+        )
         to_rename = set([*ds.dims, *ds.coords, *ds])
         ds = ds.rename({d: remove_prefix(d, prefix) for d in to_rename})
         ds.attrs = {
diff --git a/test/test_io.py b/test/test_io.py
index dabee6f5..9b07fe05 100644
--- a/test/test_io.py
+++ b/test/test_io.py
@@ -416,6 +416,32 @@ def test_model_to_netcdf_with_multiindex_scipy_engine(
     assert_model_equal(m, read_netcdf(fn))
 
 
+def test_model_to_netcdf_aux_coords(tmp_path: Path) -> None:
+    # Flat dimension with levels as auxiliary coordinates (the v1 alternative
+    # to a MultiIndex), on a variable, a constraint, a named expression and a
+    # parameter.
+    m = Model()
+    coords = xr.Coordinates(
+        {
+            "snapshot": np.arange(4),
+            "time": ("snapshot", pd.date_range("2020-01-01", periods=4, freq="h")),
+            "period": ("snapshot", [2020, 2020, 2030, 2030]),
+        }
+    )
+    x = m.add_variables(lower=0, coords=coords, name="x")
+    m.add_constraints(x >= 1, name="c")
+    m.add_expressions(2 * x, name="e")
+    # Scalar coords are attached to every container on read; the dash name
+    # makes this container's prefix extend the prefix of "c".
+    m.add_constraints(x.sel(snapshot=0) >= 2, name="c-first")
+    m.add_objective(x.sum())
+    m.parameters["demand"] = xr.DataArray(np.ones(4), coords=coords)
+
+    fn = tmp_path / "test_aux_coords.nc"
+    m.to_netcdf(fn)
+    assert_model_equal(m, read_netcdf(fn))
+
+
 def test_model_to_netcdf_with_expressions(
     model_with_expressions: Model, tmp_path: Path
 ) -> None: