diff --git a/flixopt/calculation.py b/flixopt/calculation.py
index c912b083b..a695b285b 100644
--- a/flixopt/calculation.py
+++ b/flixopt/calculation.py
@@ -164,7 +164,7 @@ def solve(self, solver: _Solver, log_file: pathlib.Path | None = None, log_main_
             from .io import document_linopy_model
 
             document_linopy_model(self.model, paths.model_documentation)
-            self.flow_system.to_netcdf(paths.flow_system)
+            self.flow_system.to_netcdf(paths.flow_system, engine='h5netcdf')
             raise RuntimeError(
                 f'Model was infeasible. Please check {paths.model_documentation=} and {paths.flow_system=} for more information.'
             )
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 63b6c6c52..4fedbddfd 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -211,7 +211,13 @@ def as_dataset(self, constants_in_dataset: bool = False) -> xr.Dataset:
         ds.attrs = self.as_dict(data_mode='name')
         return ds
 
-    def to_netcdf(self, path: str | pathlib.Path, compression: int = 0, constants_in_dataset: bool = True):
+    def to_netcdf(
+        self,
+        path: str | pathlib.Path,
+        compression: int = 0,
+        constants_in_dataset: bool = True,
+        engine: str = 'h5netcdf',
+    ):
         """
         Saves the FlowSystem to a netCDF file.
         Args:
@@ -220,7 +226,7 @@ def to_netcdf(self, path: str | pathlib.Path, compression: int = 0, constants_in
             constants_in_dataset: If True, constants are included as Dataset variables.
         """
         ds = self.as_dataset(constants_in_dataset=constants_in_dataset)
-        fx_io.save_dataset_to_netcdf(ds, path, compression=compression)
+        fx_io.save_dataset_to_netcdf(ds, path, compression=compression, engine='h5netcdf')
         logger.info(f'Saved FlowSystem to {path}')
 
     def plot_network(
diff --git a/flixopt/io.py b/flixopt/io.py
index 35304634d..07e990c45 100644
--- a/flixopt/io.py
+++ b/flixopt/io.py
@@ -208,6 +208,7 @@ def save_dataset_to_netcdf(
     ds: xr.Dataset,
     path: str | pathlib.Path,
     compression: int = 0,
+    engine: str = 'h5netcdf',
 ) -> None:
     """
     Save a dataset to a netcdf file. Store the attrs as a json string in the 'attrs' attribute.
@@ -240,6 +241,7 @@ def save_dataset_to_netcdf(
         encoding=None
         if not apply_encoding
         else {data_var: {'zlib': True, 'complevel': compression} for data_var in ds.data_vars},
+        engine='h5netcdf',
     )
 
 
diff --git a/flixopt/results.py b/flixopt/results.py
index 5a76e3e0a..a34875381 100644
--- a/flixopt/results.py
+++ b/flixopt/results.py
@@ -328,8 +328,8 @@ def to_file(
 
         paths = fx_io.CalculationResultsPaths(folder, name)
 
-        fx_io.save_dataset_to_netcdf(self.solution, paths.solution, compression=compression)
-        fx_io.save_dataset_to_netcdf(self.flow_system, paths.flow_system, compression=compression)
+        fx_io.save_dataset_to_netcdf(self.solution, paths.solution, compression=compression, engine='h5netcdf')
+        fx_io.save_dataset_to_netcdf(self.flow_system, paths.flow_system, compression=compression, engine='h5netcdf')
 
         with open(paths.summary, 'w', encoding='utf-8') as f:
             yaml.dump(self.summary, f, allow_unicode=True, sort_keys=False, indent=4, width=1000)
@@ -338,7 +338,7 @@ def to_file(
             if self.model is None:
                 logger.critical('No model in the CalculationResults. Saving the model is not possible.')
             else:
-                self.model.to_netcdf(paths.linopy_model)
+                self.model.to_netcdf(paths.linopy_model, engine='h5netcdf')
 
         if document_model:
             if self.model is None:
diff --git a/pyproject.toml b/pyproject.toml
index dcc7ffa12..1a3ece7ca 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -37,7 +37,7 @@ dependencies = [
     "xarray >= 2024.2.0, < 2026.0",  # CalVer: allow through next calendar year
     # Optimization and data handling
     "linopy >= 0.5.1, < 0.6",  # Widened from patch pin to minor range
-    "netcdf4 >= 1.6.1, < 2",
+    "h5netcdf>=1.0.0, < 2",
     # Utilities
     "pyyaml >= 6.0.0, < 7",
     "rich >= 13.0.0, < 15",