diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 96708027..90b71641 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -28,6 +28,7 @@ Upcoming Version
 * Add ``format_labels()`` on ``Constraints``/``Variables`` and ``format_infeasibilities()`` on ``Model`` that return strings instead of printing to stdout, allowing usage with logging, storage, or custom output handling. Deprecate ``print_labels()`` and ``print_infeasibilities()``.
 * Add ``fix()``, ``unfix()``, and ``fixed`` to ``Variable`` and ``Variables`` for fixing variables to values via equality constraints. Supports automatic rounding for integer/binary variables.
 * Add ``relax()``, ``unrelax()``, and ``relaxed`` to ``Variable`` and ``Variables`` for LP relaxation of integer/binary variables. Supports partial relaxation via filtered views (e.g. ``m.variables.integers.relax()``). Semi-continuous variables raise ``NotImplementedError``.
+* Add compatibility to latest xarray version.
 
 
 Version 0.6.6
diff --git a/linopy/expressions.py b/linopy/expressions.py
index 37d3e6bd..bda0b896 100644
--- a/linopy/expressions.py
+++ b/linopy/expressions.py
@@ -284,7 +284,7 @@ def sum(self, use_fallback: bool = False, **kwargs: Any) -> LinearExpression:
                 index.names = [str(col) for col in orig_group.columns]
                 index.name = GROUP_DIM
                 new_coords = Coordinates.from_pandas_multiindex(index, GROUP_DIM)
-                ds = xr.Dataset(ds.assign_coords(new_coords))
+                ds = ds.assign_coords(new_coords)
 
             ds = ds.rename({GROUP_DIM: final_group_name})
             return LinearExpression(ds, self.model)
@@ -391,8 +391,7 @@ def __init__(self, data: Dataset | Any | None, model: Model) -> None:
         coeffs_vars_dict = {str(k): v for k, v in coeffs_vars.items()}
         data = assign_multiindex_safe(data, **coeffs_vars_dict)
 
-        # transpose with new Dataset to really ensure correct order
-        data = Dataset(data.transpose(..., TERM_DIM))
+        data = data.transpose(..., TERM_DIM)
 
         # ensure helper dimensions are not set as coordinates
         if drop_dims := set(HELPER_DIMS).intersection(data.coords):
@@ -2098,7 +2097,7 @@ def __init__(self, data: Dataset | None, model: Model) -> None:
             raise ValueError(f"Size of dimension {FACTOR_DIM} must be 2.")
 
         # transpose data to have _term as last dimension and _factor as second last
-        data = xr.Dataset(data.transpose(..., FACTOR_DIM, TERM_DIM))
+        data = data.transpose(..., FACTOR_DIM, TERM_DIM)
         self._data = data
 
     @property
diff --git a/linopy/model.py b/linopy/model.py
index 2a635680..4748db6f 100644
--- a/linopy/model.py
+++ b/linopy/model.py
@@ -359,7 +359,9 @@ def parameters(self, value: Dataset | Mapping) -> None:
         """
         Set the parameters of the model.
         """
-        self._parameters = Dataset(value)
+        self._parameters = (
+            value.copy() if isinstance(value, Dataset) else Dataset(value)
+        )
 
     @property
     def solution(self) -> Dataset:
diff --git a/pyproject.toml b/pyproject.toml
index 87e79ab2..f012ebc6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -31,7 +31,7 @@ dependencies = [
     "bottleneck",
     "toolz",
     "numexpr",
-    "xarray>=2024.2.0,<2026.4",
+    "xarray>=2024.2.0",
     "dask>=0.18.0",
     "polars>=1.31.1",
     "tqdm",