diff --git a/benchmarks/benchmark_snippets.py b/benchmarks/benchmark_snippets.py
index 3bfd4d2..9762106 100644
--- a/benchmarks/benchmark_snippets.py
+++ b/benchmarks/benchmark_snippets.py
@@ -22,6 +22,11 @@
 num_cells_params = [10_000_000]
 num_cells_ids = ["10M"]
 
+maximum_shapes = {
+    "1M": (1_000, 1_000),
+    "10M": (2_000, 5_000),
+}
+
 
 def make_array_shape(num_cells: int) -> tuple[int, int]:
     side_length = int(math.sqrt(num_cells))
@@ -44,6 +49,59 @@ def annotate_math_benchmark(benchmark, *, num_cells: int, formula: str, math_pat
     benchmark.extra_info["math_path"] = math_path
 
 
+def maximum_numpy(arr_maximum: np.ndarray, arr_values: np.ndarray) -> None:
+    np.maximum(arr_maximum, arr_values, out=arr_maximum)
+
+
+def setup_maximum_args(
+    shape: tuple[int, int], padded_domain: bool
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    rng = np.random.default_rng(72)
+    if padded_domain:
+        storage_shape = (shape[0] + 2, shape[1] + 2)
+        initial = rng.random(size=storage_shape, dtype=np.float32)
+        values = rng.random(size=storage_shape, dtype=np.float32)
+        arr_maximum = np.empty(storage_shape, dtype=np.float32)
+    else:
+        initial = rng.random(size=shape, dtype=np.float32)
+        values = rng.random(size=shape, dtype=np.float32)
+        arr_maximum = np.empty(shape, dtype=np.float32)
+    return initial, values, arr_maximum
+
+
+def run_maximum_benchmark(
+    benchmark,
+    implementation,
+    shape: tuple[int, int],
+    padded_domain: bool,
+    use_interior_view: bool,
+) -> None:
+    initial, values, arr_maximum = setup_maximum_args(shape, padded_domain)
+    if padded_domain and use_interior_view:
+        initial = initial[1:-1, 1:-1]
+        values = values[1:-1, 1:-1]
+        arr_maximum = arr_maximum[1:-1, 1:-1]
+    expected = np.maximum(initial, values)
+
+    np.copyto(arr_maximum, initial)
+    implementation(arr_maximum, values)
+    np.testing.assert_array_equal(arr_maximum, expected)
+
+    def setup():
+        np.copyto(arr_maximum, initial)
+        return (arr_maximum, values), {}
+
+    benchmark.pedantic(implementation, setup=setup, rounds=7, iterations=1)
+    benchmark.extra_info["lattice_updates"] = shape[0] * shape[1]
+    benchmark.extra_info["domain_layout"] = "padded" if padded_domain else "packed"
+    if not padded_domain:
+        benchmark.extra_info["call_layout"] = "packed"
+    elif use_interior_view:
+        benchmark.extra_info["call_layout"] = "interior_view"
+    else:
+        benchmark.extra_info["call_layout"] = "padded_storage"
+
+
 def setup_almeida_args(num_cells: int) -> tuple[object, ...]:
     arr_shape = make_array_shape(num_cells)
     rng = np.random.default_rng(42)
@@ -85,6 +143,39 @@ def setup_velocity_diagnostics_args(num_cells: int) -> tuple[object, ...]:
     return arr_qx, arr_qy, arr_h, arr_v, arr_vdir, arr_fr, g
 
 
+## element-wise maximum ##
+
+
+@pytest.mark.parametrize("shape", maximum_shapes.values(), ids=maximum_shapes.keys())
+@pytest.mark.parametrize("padded_domain", [False, True], ids=["packed", "padded_domain"])
+def test_benchmark_maximum_numpy(benchmark, shape, padded_domain):
+    run_maximum_benchmark(benchmark, maximum_numpy, shape, padded_domain, use_interior_view=True)
+
+
+@pytest.mark.parametrize("shape", maximum_shapes.values(), ids=maximum_shapes.keys())
+@pytest.mark.parametrize("padded_domain", [False, True], ids=["packed", "padded_domain"])
+def test_benchmark_maximum_cython_serial(benchmark, shape, padded_domain):
+    run_maximum_benchmark(
+        benchmark,
+        snippets.arr_maximum_serial,
+        shape,
+        padded_domain,
+        use_interior_view=False,
+    )
+
+
+@pytest.mark.parametrize("shape", maximum_shapes.values(), ids=maximum_shapes.keys())
+@pytest.mark.parametrize("padded_domain", [False, True], ids=["packed", "padded_domain"])
+def test_benchmark_maximum_cython_parallel(benchmark, shape, padded_domain):
+    run_maximum_benchmark(
+        benchmark,
+        snippets.arr_maximum_parallel,
+        shape,
+        padded_domain,
+        use_interior_view=False,
+    )
+
+
 ## velocity ##
 
 
diff --git a/src/itzi_core/array_definitions.py b/src/itzi_core/array_definitions.py
index cdc08d9..3471909 100644
--- a/src/itzi_core/array_definitions.py
+++ b/src/itzi_core/array_definitions.py
@@ -293,7 +293,7 @@ class ArrayDefinition:
         key="hmax",
         csdms_name="land_surface_water__max_of_depth",
         cf_name="",
-        category=[ArrayCategory.INTERNAL],  # output together with depth arrays
+        category=[ArrayCategory.INTERNAL, ArrayCategory.OUTPUT],
         description="Maximum water depth reached since the beginning of the simulation.",
         unit="m",
         cf_unit="",
@@ -323,7 +323,7 @@ class ArrayDefinition:
         key="vmax",
         csdms_name="land_surface_water_flow__max_of_speed",
         cf_name="",
-        category=[ArrayCategory.INTERNAL],  # output together with speed arrays
+        category=[ArrayCategory.INTERNAL, ArrayCategory.OUTPUT],
         description="Maximum water speed reached since the beginning of the simulation.",
         unit="m s-1",
         cf_unit="",
diff --git a/src/itzi_core/compute/snippets.pyx b/src/itzi_core/compute/snippets.pyx
index 92268a6..781b7de 100644
--- a/src/itzi_core/compute/snippets.pyx
+++ b/src/itzi_core/compute/snippets.pyx
@@ -44,6 +44,50 @@ cdef float DEG_360_F = 360.0
 cdef float EPS_F = 1e-12
 
 
+@cython.wraparound(False)
+@cython.boundscheck(False)
+@cython.initializedcheck(False)
+@cython.nonecheck(False)
+def arr_maximum_serial(
+    DTYPE_t[:, ::1] arr_maximum,
+    DTYPE_t[:, ::1] arr_values,
+):
+    """Update an element-wise maximum with a serial loop for benchmarking."""
+    cdef Py_ssize_t rows = arr_maximum.shape[0]
+    cdef Py_ssize_t cols = arr_maximum.shape[1]
+    cdef Py_ssize_t row_idx, col_idx
+
+    if rows != arr_values.shape[0] or cols != arr_values.shape[1]:
+        raise ValueError("Maximum arrays must have the same shape")
+
+    for row_idx in range(rows):
+        for col_idx in range(cols):
+            if arr_values[row_idx, col_idx] > arr_maximum[row_idx, col_idx]:
+                arr_maximum[row_idx, col_idx] = arr_values[row_idx, col_idx]
+
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+@cython.initializedcheck(False)
+@cython.nonecheck(False)
+def arr_maximum_parallel(
+    DTYPE_t[:, ::1] arr_maximum,
+    DTYPE_t[:, ::1] arr_values,
+):
+    """Update an element-wise maximum with OpenMP for benchmarking."""
+    cdef Py_ssize_t rows = arr_maximum.shape[0]
+    cdef Py_ssize_t cols = arr_maximum.shape[1]
+    cdef Py_ssize_t row_idx, col_idx
+
+    if rows != arr_values.shape[0] or cols != arr_values.shape[1]:
+        raise ValueError("Maximum arrays must have the same shape")
+
+    for row_idx in prange(rows, nogil=True, schedule="static"):
+        for col_idx in range(cols):
+            if arr_values[row_idx, col_idx] > arr_maximum[row_idx, col_idx]:
+                arr_maximum[row_idx, col_idx] = arr_values[row_idx, col_idx]
+
+
 @cython.wraparound(False)  # Disable negative index check
 @cython.cdivision(True)  # Don't check division by zero
 @cython.boundscheck(False)  # turn off bounds-checking for entire function
diff --git a/src/itzi_core/providers/base.py b/src/itzi_core/providers/base.py
index 64b2f8a..e05a2b6 100644
--- a/src/itzi_core/providers/base.py
+++ b/src/itzi_core/providers/base.py
@@ -23,7 +23,7 @@
 
     import numpy as np
 
-    from itzi_core.data_containers import DrainageNetworkData, MassBalanceData, SimulationData
+    from itzi_core.data_containers import DrainageNetworkData, MassBalanceData
     from itzi_core.providers.domain_data import DomainData
 
 
@@ -61,9 +61,8 @@ def write_arrays(
     ) -> None:
         """Write all arrays for the current time step."""
 
-    @abstractmethod
-    def finalize(self, final_data: SimulationData) -> None:
-        """Finalize outputs and cleanup."""
+    def finalize(self) -> None:
+        """Flush and close provider resources."""
 
 
 class VectorOutputProvider(ABC):
@@ -75,9 +74,8 @@ def write_vector(
     ) -> None:
         """Write simulation data for current time step."""
 
-    @abstractmethod
-    def finalize(self, drainage_data: DrainageNetworkData) -> None:
-        """Finalize outputs and cleanup."""
+    def finalize(self) -> None:
+        """Flush and close provider resources."""
 
 
 class MassBalanceOutputProvider(ABC):
diff --git a/src/itzi_core/providers/csv_output.py b/src/itzi_core/providers/csv_output.py
index 3731170..f11ff05 100644
--- a/src/itzi_core/providers/csv_output.py
+++ b/src/itzi_core/providers/csv_output.py
@@ -14,17 +14,21 @@
 
 from __future__ import annotations
 
+import csv
 from datetime import datetime, timedelta
-from typing import TypedDict, TYPE_CHECKING, Tuple, List
 from io import StringIO
 from pathlib import PurePosixPath, PureWindowsPath
-import csv
+from typing import TYPE_CHECKING, TypedDict
 
 import pandas as pd
 
+from itzi_core.data_containers import (
+    DrainageLinkAttributes,
+    DrainageLinkData,
+    DrainageNodeAttributes,
+    DrainageNodeData,
+)
 from itzi_core.providers.base import VectorOutputProvider
-from itzi_core.data_containers import DrainageLinkData, DrainageLinkAttributes
-from itzi_core.data_containers import DrainageNodeData, DrainageNodeAttributes
 
 if TYPE_CHECKING:
     from itzi_core.data_containers import DrainageNetworkData
@@ -117,9 +121,6 @@ def write_vector(
         self._update_csv(sim_time_str, "link", drainage_data.links)
         self.number_of_writes["link"] += 1
 
-    def finalize(self, drainage_data: DrainageNetworkData) -> None:
-        """Finalize outputs and cleanup."""
-
     def _check_existing_csv(self, geom_type: str):
         """In order to be compatible, an existing CSV should have:
             - Same headers
@@ -206,18 +207,16 @@ def _update_csv(
         self,
         sim_time_str: str,
         geom_type: str,
-        drainage_elements: Tuple[DrainageNodeData | DrainageLinkData, ...],
+        drainage_elements: tuple[DrainageNodeData | DrainageLinkData, ...],
     ):
         """Update adequate CSV in object store"""
         # Check compatibility on first write
         if 0 == self.number_of_writes[geom_type] and self.existing_ids[geom_type]:
             # IDs must match
-            new_ids = set(
-                [
-                    drainage_elem.attributes.model_dump()[f"{geom_type}_id"]
-                    for drainage_elem in drainage_elements
-                ]
-            )
+            new_ids = {
+                drainage_elem.attributes.model_dump()[f"{geom_type}_id"]
+                for drainage_elem in drainage_elements
+            }
             if not new_ids == self.existing_ids[geom_type]:
                 raise ValueError(
                     f"Object ids mismatch for {geom_type}: "
@@ -236,7 +235,7 @@ def _update_csv(
         updated_csv = existing_csv + new_rows
         obstore.put(self.store, self.file_paths[geom_type], file=updated_csv.encode("utf-8"))
 
-    def _attrs_line(self, drainage_element: DrainageNodeData | DrainageLinkData) -> List[str, ...]:
+    def _attrs_line(self, drainage_element: DrainageNodeData | DrainageLinkData) -> list[str, ...]:
         """Return a list of attributes"""
         # Convert attributes to list
         attributes = [str(a) for a in drainage_element.attributes.model_dump().values()]
diff --git a/src/itzi_core/providers/icechunk_output.py b/src/itzi_core/providers/icechunk_output.py
index e06866f..50a56e5 100644
--- a/src/itzi_core/providers/icechunk_output.py
+++ b/src/itzi_core/providers/icechunk_output.py
@@ -15,7 +15,7 @@
 from collections.abc import Mapping
 from datetime import datetime, timedelta
 from importlib.metadata import version
-from typing import TYPE_CHECKING, TypedDict
+from typing import TypedDict
 
 import numpy as np
 
@@ -25,6 +25,8 @@
     import pyproj
     import xarray as xr
     import zarr
+    from pyproj.exceptions import CRSError
+    from zarr.errors import BaseZarrError, GroupNotFoundError
 except ImportError:
     raise ImportError(
         "To use the Icechunk backend, install itzi with: "
@@ -35,9 +37,6 @@
 from itzi_core.array_definitions import ARRAY_DEFINITIONS
 from itzi_core.providers.base import RasterOutputProvider
 
-if TYPE_CHECKING:
-    from itzi_core.data_containers import SimulationData
-
 
 class IcechunkRasterOutputConfig(TypedDict):
     # A list of var names to be written
@@ -68,6 +67,7 @@ def __init__(self, config: IcechunkRasterOutputConfig) -> None:
                 raise
         self.spatial_coordinates = self._get_spatial_coordinates()
         self.append_mode = False
+        self.time_is_datetime: bool | None = None
         # Make sure new data matches existing one
         if self.has_existing_data():
             self.check_repo_match()
@@ -100,106 +100,102 @@ def _get_spatial_coordinates(self) -> list[tuple[str, np.ndarray, dict[str, str]
 
     def has_existing_data(self) -> bool:
         """Check if the repository already contains data."""
+        session = self.repo.readonly_session("main")
         try:
-            session = self.repo.readonly_session("main")
             existing_ds = xr.open_zarr(session.store)
-            # Check if there are any data variables and if time dimension has records
-            has_vars = len(existing_ds.data_vars) > 0
-            has_time_records = "time" in existing_ds.coords and len(existing_ds.coords["time"]) > 0
-            return has_vars and has_time_records
-        except Exception:
+        except GroupNotFoundError:
             return False
+        return len(existing_ds.data_vars) > 0
 
     def get_latest_timestamp(self) -> datetime | timedelta | None:
         """Get the latest timestamp from existing data, or None if no data exists."""
         session = self.repo.readonly_session("main")
         try:
             existing_ds = xr.open_zarr(session.store)
-            if "time" in existing_ds.coords and len(existing_ds.coords["time"]) > 0:
-                latest_time = existing_ds.coords["time"][-1]
-                # Convert numpy datetime64/timedelta64 back to Python types
-                if np.issubdtype(latest_time.dtype, np.datetime64):
-                    return latest_time.values.astype("datetime64[ms]").astype(datetime)
-                elif np.issubdtype(latest_time.dtype, np.timedelta64):
-                    return timedelta(
-                        seconds=int(latest_time.values.astype("timedelta64[s]").astype(int))
-                    )
-                else:
-                    return None
-            else:
-                return None
-        except Exception:
+        except GroupNotFoundError:
             return None
+        if "time" in existing_ds.coords and len(existing_ds.coords["time"]) > 0:
+            latest_time = existing_ds.coords["time"][-1]
+            # Convert numpy datetime64/timedelta64 back to Python types
+            if np.issubdtype(latest_time.dtype, np.datetime64):
+                return latest_time.values.astype("datetime64[ms]").astype(datetime)
+            elif np.issubdtype(latest_time.dtype, np.timedelta64):
+                return timedelta(
+                    milliseconds=int(latest_time.values.astype("timedelta64[ms]").astype(int))
+                )
+        return None
 
-    def check_repo_match(self):
+    def check_repo_match(self) -> None:
         """Raises ValueError if entry data does not match the existing repo."""
-        crs_match = False
-        vars_match = False
-        repo_match = False
-
-        # Get crs, variables, dims and coordinates from existing repo
         session = self.repo.readonly_session("main")
         try:
             existing_ds = xr.open_zarr(session.store)
-        except Exception as e:
+        except BaseZarrError as e:
             raise ValueError(f"Existing {session.store} is not a valid zarr store: {e}")
 
         try:
-            existing_crs = pyproj.CRS.from_wkt(existing_ds.attrs["crs_wkt"])
-        except Exception:
-            crs_match = False
-
-        existing_vars = existing_ds.data_vars
-        existing_num_vars = len(existing_vars)
-        # Check if existing variables coordinates match the new ones
-        # This implies coherence in variables names and dims
-        if existing_crs == self.crs:
-            crs_match = True
-        new_num_vars = len(self.out_map_names)
-        vars_match_dict = {}
-        for existing_var_name in existing_vars:
-            existing_var = existing_ds[existing_var_name]
-            try:
-                # incompatible dimension names will fail here
-                existing_var_x_coords = existing_var.coords["x"].values
-                existing_var_y_coords = existing_var.coords["y"].values
-            except Exception:
-                vars_match_dict[existing_var_name] = False
-                continue
-            try:  # allclose raises ValueError if not same len
-                var_x_coords_match = np.allclose(existing_var_x_coords, self.x_coords)
-                var_y_coords_match = np.allclose(existing_var_y_coords, self.y_coords)
-                if var_x_coords_match and var_y_coords_match:
-                    vars_match_dict[existing_var_name] = True
-                else:
-                    vars_match_dict[existing_var_name] = False
-            except ValueError:
-                vars_match_dict[existing_var_name] = False
-            if existing_var_name not in self.out_map_names.values():
-                vars_match_dict[existing_var_name] = False
-
-        if all(list(vars_match_dict.values())) and existing_num_vars == new_num_vars:
-            vars_match = True
-        # Raise if not full match
-        repo_match = crs_match and vars_match
-        if not repo_match:
+            existing_crs_wkt = existing_ds.attrs["crs_wkt"]
+        except KeyError as e:
+            raise KeyError("Existing repository has no 'crs_wkt' attribute") from e
+        try:
+            existing_crs = pyproj.CRS.from_wkt(existing_crs_wkt)
+        except CRSError as e:
+            raise ValueError("Existing repository 'crs_wkt' attribute is not valid WKT") from e
+        if existing_crs != self.crs:
             raise ValueError(
-                f"Provided data does not match existing icechunk repository {self.repo}: "
-                f"CRS match: {crs_match}, "
-                f"Variables match: {vars_match}, "
-                f"Existing CRS: {existing_crs.to_epsg()}, "
-                f"New CRS: {self.crs.to_epsg()}, "
-                f"Existing numbers of variables in the store: {existing_num_vars}, "
-                f"Numbers of variables to be written: {new_num_vars}, "
-                f"Matching variables coordinates: {vars_match_dict}, "
-                f"Existing vars: {[v for v in existing_vars]}, "
-                f"New vars: {list(self.out_map_names.values())}."
+                "Provided CRS does not match existing icechunk repository: "
+                f"existing={existing_crs.to_epsg()}, configured={self.crs.to_epsg()}"
             )
 
+        existing_names = set(existing_ds.data_vars)
+        expected_names = set(self.out_map_names.values())
+        if existing_names != expected_names:
+            raise ValueError(
+                "Configured output names do not match existing icechunk repository: "
+                f"existing={sorted(existing_names)}, configured={sorted(expected_names)}"
+            )
+
+        expected_dims = ("time", "y", "x")
+        for name, variable in existing_ds.data_vars.items():
+            if variable.dims != expected_dims:
+                raise ValueError(
+                    f"Existing variable {name!r} has dimensions {variable.dims}; expected "
+                    f"{expected_dims}. Create a new repository or explicitly migrate the "
+                    "legacy repository before appending."
+                )
+            try:
+                x_matches = np.allclose(variable.coords["x"].values, self.x_coords)
+                y_matches = np.allclose(variable.coords["y"].values, self.y_coords)
+            except (KeyError, ValueError):
+                x_matches = y_matches = False
+            if not x_matches or not y_matches:
+                raise ValueError(
+                    f"Coordinates for existing variable {name!r} do not match the configured grid"
+                )
+
+        if "time" not in existing_ds.coords:
+            raise ValueError("Existing repository has no compatible temporal time coordinate")
+        time_dtype = existing_ds["time"].dtype
+        if np.issubdtype(time_dtype, np.datetime64):
+            self.time_is_datetime = True
+        elif np.issubdtype(time_dtype, np.timedelta64):
+            self.time_is_datetime = False
+        else:
+            raise ValueError("Existing repository has no compatible temporal time coordinate")
+
     def write_arrays(
         self, array_dict: Mapping[str, np.ndarray], sim_time: datetime | timedelta
     ) -> None:
         """Write results to an icechunk repository"""
+        incoming_time_is_datetime = isinstance(sim_time, datetime)
+        if (
+            self.time_is_datetime is not None
+            and self.time_is_datetime != incoming_time_is_datetime
+        ):
+            raise ValueError(
+                "Incoming time coordinate type does not match the existing icechunk repository"
+            )
+        self.time_is_datetime = incoming_time_is_datetime
         vars_to_write = list(array_dict.keys())
         expected_map_keys = set(self.out_map_names.keys())
         if not expected_map_keys == set(vars_to_write):
@@ -209,7 +205,7 @@ def write_arrays(
                 f"Received: {vars_to_write}"
             )
         # prepare the data
-        dataset = self.get_dataset_from_dict(array_dict, sim_time)
+        dataset = self._build_dataset(array_dict, sim_time)
         first_var_name = next(iter(self.out_map_names.values()))
         time_encoding = dataset[first_var_name].encoding["time"]
         # Commit to the repo
@@ -225,53 +221,31 @@ def write_arrays(
             self._zarr_append(icechunk_session.store, dataset)
         icechunk_session.commit(commit_message)
 
-    def write_maxs(self, array_dict):
-        # Validate that all provided variables are valid max variables
-        valid_max_vars = ["hmax", "vmax"]
-        vars_to_write = list(array_dict.keys())
-        invalid_vars = [var for var in vars_to_write if var not in valid_max_vars]
-        if invalid_vars:
-            raise ValueError(
-                f"Invalid max variables: {invalid_vars}. Valid options: {valid_max_vars}"
-            )
-        dataset = self.get_dataset_from_dict(array_dict)
-        # Commit to the repo
-        icechunk_session = self.repo.writable_session("main")
-        icechunk.xarray.to_icechunk(dataset, icechunk_session, mode="a")
-        icechunk_session.commit("Maximum values of itzi simulation")
-
-    def get_dataset_from_dict(self, array_dict, sim_time=None):
+    def _build_dataset(
+        self, array_dict: Mapping[str, np.ndarray], sim_time: datetime | timedelta
+    ) -> xr.Dataset:
         """From a dict of arrays, return an xarray dataset."""
         data_vars = {}
-        if sim_time is not None:
-            if isinstance(sim_time, datetime):
-                time_dtype = "datetime64[ms]"
-                sim_time_np = np.datetime64(sim_time, "ms")
-                time_unit = "milliseconds since 1970-01-01T00:00:00"
-            elif isinstance(sim_time, timedelta):
-                time_dtype = "timedelta64[s]"
-                sim_time_np = np.timedelta64(sim_time, "s")
-                time_unit = "seconds"
-            else:
-                raise ValueError(f"Unknown temporal type: {type(sim_time)}")
-
-            time_coordinate = np.array([sim_time_np], dtype=time_dtype)
-            # Don't put any encoding-related attrs since they conflict with encoding
-            time_attrs = {}
-            time_encoding = {
-                "units": time_unit,
-                "dtype": time_dtype,
-            }
-            coordinates = [("time", time_coordinate, time_attrs)] + self.spatial_coordinates
+        if isinstance(sim_time, datetime):
+            time_dtype = "datetime64[ms]"
+            sim_time_np = np.datetime64(sim_time, "ms")
+            time_unit = "milliseconds since 1970-01-01T00:00:00"
+        elif isinstance(sim_time, timedelta):
+            time_dtype = "timedelta64[ms]"
+            sim_time_np = np.timedelta64(sim_time, "ms")
+            time_unit = "milliseconds"
         else:
-            coordinates = self.spatial_coordinates
+            raise ValueError(f"Unknown temporal type: {type(sim_time)}")
+
+        time_coordinate = np.array([sim_time_np], dtype=time_dtype)
+        time_encoding = {
+            "units": time_unit,
+            "dtype": time_dtype,
+        }
+        coordinates = [("time", time_coordinate, {})] + self.spatial_coordinates
 
         for key, arr in array_dict.items():
-            if key in ["hmax", "vmax"]:
-                max_mapping = {"hmax": "water_depth", "vmax": "v"}
-                var_name = f"{self.out_map_names[max_mapping[key]]}_max"
-            else:
-                var_name = self.out_map_names[key]
+            var_name = self.out_map_names[key]
             coords_shape = (len(self.y_coords), len(self.x_coords))
             if arr.shape != coords_shape:
                 raise ValueError(
@@ -283,8 +257,7 @@ def get_dataset_from_dict(self, array_dict, sim_time=None):
                 "standard_name": self.cf_names[key],
                 "long_name": self.descriptions[key],
             }
-            if sim_time is not None:
-                arr = np.expand_dims(arr, axis=0)
+            arr = np.expand_dims(arr, axis=0)
 
             data_array = xr.DataArray(
                 data=arr,
@@ -292,9 +265,8 @@ def get_dataset_from_dict(self, array_dict, sim_time=None):
                 name=var_name,  # Write the requested name
                 attrs=var_attributes,
             )
-            if sim_time is not None:
-                assert data_array["time"].dtype == time_dtype
-                data_array.encoding["time"] = time_encoding
+            assert data_array["time"].dtype == time_dtype
+            data_array.encoding["time"] = time_encoding
             data_vars[var_name] = data_array
         dataset_attributes = {
             "crs_wkt": self.crs.to_wkt(),
@@ -302,9 +274,7 @@ def get_dataset_from_dict(self, array_dict, sim_time=None):
         }
         dataset = xr.Dataset(data_vars, attrs=dataset_attributes)
 
-        # Set encoding on the time coordinate of the dataset itself
-        if sim_time is not None:
-            dataset["time"].encoding.update(time_encoding)
+        dataset["time"].encoding.update(time_encoding)
 
         return dataset
 
@@ -328,15 +298,3 @@ def _zarr_append(self, store, dataset: xr.Dataset) -> None:
             z_group[var_name].resize(new_shape)
             # Use direct assignment
             z_group[var_name][current_shape[0]] = data_array.values[0]
-
-    def finalize(self, final_data: "SimulationData") -> None:
-        """Write max values."""
-        arr_dict = {}
-        # Only process arrays that are actually configured for output
-        if self.out_map_names.get("water_depth"):
-            arr_dict["hmax"] = final_data.raw_arrays["hmax"]
-        if self.out_map_names.get("v"):
-            arr_dict["vmax"] = final_data.raw_arrays["vmax"]
-        # Only write max values if there are arrays to write
-        if arr_dict:
-            self.write_maxs(arr_dict)
diff --git a/src/itzi_core/providers/memory_output.py b/src/itzi_core/providers/memory_output.py
index 34f1fa8..b681dce 100644
--- a/src/itzi_core/providers/memory_output.py
+++ b/src/itzi_core/providers/memory_output.py
@@ -18,7 +18,7 @@
 
 import numpy as np
 
-from itzi_core.data_containers import DrainageNetworkData, SimulationData
+from itzi_core.data_containers import DrainageNetworkData
 from itzi_core.providers.base import RasterOutputProvider, VectorOutputProvider
 
 
@@ -38,9 +38,6 @@ def write_arrays(
             if isinstance(arr, np.ndarray):
                 self.output_maps_dict[arr_key].append((deepcopy(sim_time), arr.copy()))
 
-    def finalize(self, final_data: SimulationData) -> None:
-        """Finalize outputs and cleanup."""
-
 
 class MemoryVectorOutputProvider(VectorOutputProvider):
     """Save drainage simulation outputs in memory."""
@@ -54,6 +51,3 @@ def write_vector(
     ) -> None:
         """Save simulation data for current time step."""
         self.drainage_data.append((deepcopy(sim_time), deepcopy(drainage_data)))
-
-    def finalize(self, drainage_data: DrainageNetworkData) -> None:
-        """Finalize outputs and cleanup."""
diff --git a/src/itzi_core/report.py b/src/itzi_core/report.py
index 9b18300..a857429 100644
--- a/src/itzi_core/report.py
+++ b/src/itzi_core/report.py
@@ -18,6 +18,8 @@
 from datetime import datetime, timedelta
 from typing import TYPE_CHECKING
 
+import numpy as np
+
 from itzi_core.array_definitions import ARRAY_DEFINITIONS, ArrayCategory
 from itzi_core.compute import rastermetrics
 from itzi_core.const import TemporalType
@@ -56,8 +58,6 @@ def __init__(
         self.mass_balance_output_provider = mass_balance_output_provider
         # a dict containing lists of maps written to gis to be registered
         self.output_maplist = {k: [] for k in self.out_map_names}
-        # a dict of array written at a given step. Keys are the same as out_map_names
-        self.output_arrays = {}
         self.dt = dt
         self.last_step = copy.copy(start_time)
 
@@ -68,10 +68,8 @@ def step(self, simulation_data: SimulationData):
             converted_sim_time = sim_time - self.start_time
         else:
             converted_sim_time = sim_time
-        self.get_output_arrays(simulation_data)
-        self.raster_provider.write_arrays(
-            array_dict=self.output_arrays, sim_time=converted_sim_time
-        )
+        output_arrays = self.get_output_arrays(simulation_data)
+        self.raster_provider.write_arrays(array_dict=output_arrays, sim_time=converted_sim_time)
         if self.mass_balance_output_provider is not None:
             self.write_mass_balance(simulation_data, converted_sim_time)
         drainage_data = simulation_data.drainage_network_data
@@ -81,17 +79,17 @@ def step(self, simulation_data: SimulationData):
         self.last_step = copy.copy(sim_time)
         return self
 
-    def end(self, final_data: SimulationData):
+    def end(self):
         """Finalize output providers after the last report has been written."""
-        self.raster_provider.finalize(final_data)
-        if final_data.drainage_network_data is not None:
-            self.vector_provider.finalize(final_data.drainage_network_data)
+        self.raster_provider.finalize()
+        self.vector_provider.finalize()
         if self.mass_balance_output_provider is not None:
             self.mass_balance_output_provider.finalize()
         return self
 
-    def get_output_arrays(self, data: SimulationData):
+    def get_output_arrays(self, data: SimulationData) -> dict[str, np.ndarray]:
         """Returns a dict of arrays to be written to the disk"""
+        output_arrays = {}
         raw = data.raw_arrays
         accum_arrays = data.accumulation_arrays
         interval_s = (data.sim_time - self.last_step).total_seconds()
@@ -107,20 +105,20 @@ def get_output_arrays(self, data: SimulationData):
             # --- Direct raw arrays ---
             if arr_key in ["water_depth", "v", "vdir", "froude", "hmax", "vmax"]:
                 if arr_key in raw:
-                    self.output_arrays[arr_key] = raw[arr_key]
+                    output_arrays[arr_key] = raw[arr_key]
                 continue  # go to next key
 
             # --- Calculated arrays ---
             if arr_key == "water_surface_elevation":
-                self.output_arrays[arr_key] = rastermetrics.calculate_wse(
+                output_arrays[arr_key] = rastermetrics.calculate_wse(
                     raw["water_depth"], raw["dem"]
                 )
             elif arr_key == "qx":
-                self.output_arrays[arr_key] = rastermetrics.calculate_flux(raw["qe_new"], cell_dy)
+                output_arrays[arr_key] = rastermetrics.calculate_flux(raw["qe_new"], cell_dy)
             elif arr_key == "qy":
-                self.output_arrays[arr_key] = rastermetrics.calculate_flux(raw["qs_new"], cell_dx)
+                output_arrays[arr_key] = rastermetrics.calculate_flux(raw["qs_new"], cell_dx)
             elif arr_key == "volume_error":  # Volume error
-                self.output_arrays[arr_key] = accum_arrays["error_depth_accum"] * cell_area
+                output_arrays[arr_key] = accum_arrays["error_depth_accum"] * cell_area
 
         # --- Averaged accumulation arrays ---
         if interval_s <= 0:
@@ -137,10 +135,10 @@ def get_output_arrays(self, data: SimulationData):
                     conversion_factor = 1000 * 3600  # m/s to mm/h
                 else:
                     conversion_factor = 1.0
-                self.output_arrays[output_name] = rastermetrics.calculate_average_rate_from_total(
+                output_arrays[output_name] = rastermetrics.calculate_average_rate_from_total(
                     accum_arrays[accum_key], interval_s, conversion_factor
                 )
-        return self
+        return output_arrays
 
     def write_mass_balance(self, data: SimulationData, converted_sim_time: datetime | timedelta):
         """Calculate mass balance and log it."""
diff --git a/src/itzi_core/simulation.py b/src/itzi_core/simulation.py
index a80da28..a688ef8 100644
--- a/src/itzi_core/simulation.py
+++ b/src/itzi_core/simulation.py
@@ -142,17 +142,8 @@ def initialize(self) -> Self:
             padded=True,
         )
 
-        # Calculate hmax and vmax of the initial state
-        np.maximum(
-            self.raster_domain.get_array("hmax"),
-            self.raster_domain.get_array("water_depth"),
-            out=self.raster_domain.get_array("hmax"),
-        )
-        np.maximum(
-            self.raster_domain.get_array("vmax"),
-            self.raster_domain.get_array("v"),
-            out=self.raster_domain.get_array("vmax"),
-        )
+        self._update_maximum("water_depth", "hmax")
+        self._update_maximum("v", "vmax")
 
         for arr_key in self.accum_mapping.keys():
             self._update_accum_array(arr_key, self.sim_time)
@@ -287,9 +278,7 @@ def update_until(self, then: timedelta) -> Self:
     def finalize(self) -> None:
         """Flush already-written results and close runtime resources."""
         # The last interval is reported by update() when its end lands on end_time.
-        if self.continuity_data is None:
-            self.continuity_data = self.get_continuity_data()
-        self.report.end(self._build_simulation_data(self.sim_time, 0))
+        self.report.end()
         if self.drainage_model:
             self.drainage_model.close()
 
@@ -348,6 +337,10 @@ def set_array(
         if arr_id in ["inflow", "rain"]:
             self._update_accum_array(arr_id, current_time)
         self.raster_domain.update_array(arr_id, arr)
+        if arr_id in {"water_depth", "water_surface_elevation"}:
+            self._update_maximum("water_depth", "hmax")
+        elif arr_id == "v":
+            self._update_maximum("v", "vmax")
         if arr_id == "dem":
             self.surface_flow.update_flow_dir()
         return self
@@ -404,6 +397,12 @@ def _update_accum_array(self, k: str, sim_time: datetime) -> None:
             rastermetrics.accumulate_rate_to_total(accum_array, rate_array, time_diff, padded=True)
             self.accum_update_time[ak] = sim_time
 
+    def _update_maximum(self, value_key: str, maximum_key: str) -> None:
+        """Synchronize a cumulative maximum with its current value array."""
+        values = self.raster_domain.get_array(value_key)
+        maximum = self.raster_domain.get_array(maximum_key)
+        np.maximum(maximum, values, out=maximum)
+
     def create_hotstart(self) -> io.BytesIO:
         """Create a hotstart file with the current state of the simulation.
 
diff --git a/tests/test_array_definitions.py b/tests/test_array_definitions.py
index 3ca8b84..ec31622 100644
--- a/tests/test_array_definitions.py
+++ b/tests/test_array_definitions.py
@@ -53,3 +53,10 @@ def test_array_definitions():
             if attr == "cf_name" and len(duplicates) == 1 and "" in duplicates:
                 continue
             assert False, f"Found duplicates in <{attr}>: {duplicates}"
+
+
+def test_maximum_arrays_are_internal_outputs():
+    definitions = {arr_def.key: arr_def for arr_def in ARRAY_DEFINITIONS}
+    for key in ("hmax", "vmax"):
+        assert ArrayCategory.INTERNAL in definitions[key].category
+        assert ArrayCategory.OUTPUT in definitions[key].category
diff --git a/tests/test_hotstart_integration.py b/tests/test_hotstart_integration.py
index 3ef0536..9b254fd 100644
--- a/tests/test_hotstart_integration.py
+++ b/tests/test_hotstart_integration.py
@@ -170,7 +170,7 @@ def run_with_hotstart_checkpoints(
 
 
 def assert_final_state_matches(simulation: Simulation, reference: Simulation) -> None:
-    for key in ["water_depth", "qe", "qs"]:
+    for key in ["water_depth", "hmax", "vmax", "qe", "qs"]:
         arr_resumed = simulation.raster_domain.get_array(key)
         arr_reference = reference.raster_domain.get_array(key)
         np.testing.assert_allclose(arr_resumed, arr_reference, err_msg=f"Final {key} mismatch")
@@ -396,7 +396,7 @@ def test_resume_allows_output_map_name_change(
 
     resumed_output_map_names = helpers.make_output_map_names(
         "out_resume",
-        ["water_depth", "qx", "qy", "volume_error"],
+        ["water_depth", "hmax", "qx", "qy", "volume_error"],
     )
     sim_b_config = sim_a_config.model_copy(update={"output_map_names": resumed_output_map_names})
     resumed_output = MemoryRasterOutputProvider(resumed_output_map_names)
@@ -412,6 +412,10 @@ def test_resume_allows_output_map_name_change(
     assert sim_b.report.out_map_names == resumed_output_map_names
     assert resumed_output.out_map_names == resumed_output_map_names
     assert resumed_output.output_maps_dict["water_depth"]
+    assert resumed_output.output_maps_dict["hmax"]
+    np.testing.assert_allclose(
+        resumed_output.output_maps_dict["hmax"][-1][1], sim_b.get_array("hmax")
+    )
     assert sim_a_config.output_map_names["water_depth"] != resumed_output_map_names["water_depth"]
     assert_final_state_matches(sim_b, uninterrupted_simulation)
 
diff --git a/tests/test_icechunk_output.py b/tests/test_icechunk_output.py
index aa0dd55..6ed1724 100644
--- a/tests/test_icechunk_output.py
+++ b/tests/test_icechunk_output.py
@@ -13,7 +13,7 @@
 """
 
 import tempfile
-from typing import Dict, Mapping
+from collections.abc import Mapping
 from datetime import datetime, timedelta
 
 import numpy as np
@@ -26,13 +26,11 @@
 pytest.importorskip("pyproj")
 
 import icechunk
-import xarray as xr
 import pyproj
+import xarray as xr
 
-from itzi_core.providers.icechunk_output import IcechunkRasterOutputProvider
 from itzi_core.array_definitions import ARRAY_DEFINITIONS, ArrayCategory
-from itzi_core.data_containers import SimulationData
-
+from itzi_core.providers.icechunk_output import IcechunkRasterOutputProvider
 
 # Mark all tests in this module as cloud tests
 pytestmark = pytest.mark.cloud
@@ -55,7 +53,7 @@ def maps_dict():
 
 
 @pytest.fixture(scope="module")
-def coordinates(maps_dict: Dict):
+def coordinates(maps_dict: dict):
     """Generate x and y coordinates for the test arrays"""
     arr_shape = next(iter(maps_dict.values())).shape
     y_coords = np.linspace(start=1234, stop=1234 + arr_shape[0], num=arr_shape[0])
@@ -70,14 +68,14 @@ def crs():
 
 
 @pytest.fixture(scope="module")
-def out_map_names(maps_dict: Dict):
+def out_map_names(maps_dict: dict):
     """Output map names mapping for the test arrays"""
     return {key: f"test_{key}" for key in maps_dict.keys()}
 
 
 @pytest.fixture
 def icechunk_provider(
-    temp_dir: tempfile.TemporaryDirectory, coordinates: Dict, crs: pyproj.CRS, out_map_names: list
+    temp_dir: tempfile.TemporaryDirectory, coordinates: dict, crs: pyproj.CRS, out_map_names: list
 ):
     storage = icechunk.local_filesystem_storage(temp_dir.name)
     provider_config = {
@@ -91,6 +89,68 @@ def icechunk_provider(
     return icechunk_p
 
 
+def test_missing_zarr_group_is_treated_as_empty(
+    icechunk_provider: IcechunkRasterOutputProvider,
+):
+    assert not icechunk_provider.has_existing_data()
+    assert icechunk_provider.get_latest_timestamp() is None
+    with pytest.raises(ValueError, match="not a valid zarr store"):
+        icechunk_provider.check_repo_match()
+
+
+def test_missing_crs_metadata_is_rejected(
+    icechunk_provider: IcechunkRasterOutputProvider,
+    monkeypatch: pytest.MonkeyPatch,
+):
+    monkeypatch.setattr(xr, "open_zarr", lambda store: xr.Dataset(attrs={}))
+
+    with pytest.raises(KeyError, match="Existing repository has no 'crs_wkt' attribute"):
+        icechunk_provider.check_repo_match()
+
+
+def test_invalid_crs_metadata_is_rejected(
+    icechunk_provider: IcechunkRasterOutputProvider,
+    monkeypatch: pytest.MonkeyPatch,
+):
+    monkeypatch.setattr(xr, "open_zarr", lambda store: xr.Dataset(attrs={"crs_wkt": "invalid"}))
+
+    with pytest.raises(
+        ValueError, match="Existing repository 'crs_wkt' attribute is not valid WKT"
+    ):
+        icechunk_provider.check_repo_match()
+
+
+@pytest.mark.parametrize(
+    "method_name", ["has_existing_data", "get_latest_timestamp", "check_repo_match"]
+)
+def test_unexpected_zarr_errors_propagate(
+    icechunk_provider: IcechunkRasterOutputProvider,
+    monkeypatch: pytest.MonkeyPatch,
+    method_name: str,
+):
+    def raise_unexpected_error(*args, **kwargs):
+        raise RuntimeError("unexpected error")
+
+    monkeypatch.setattr(xr, "open_zarr", raise_unexpected_error)
+
+    with pytest.raises(RuntimeError, match="unexpected error"):
+        getattr(icechunk_provider, method_name)()
+
+
+def test_unexpected_crs_errors_propagate(
+    icechunk_provider: IcechunkRasterOutputProvider,
+    monkeypatch: pytest.MonkeyPatch,
+):
+    def raise_unexpected_error(*args, **kwargs):
+        raise RuntimeError("unexpected error")
+
+    monkeypatch.setattr(xr, "open_zarr", lambda store: xr.Dataset(attrs={"crs_wkt": "valid"}))
+    monkeypatch.setattr(pyproj.CRS, "from_wkt", raise_unexpected_error)
+
+    with pytest.raises(RuntimeError, match="unexpected error"):
+        icechunk_provider.check_repo_match()
+
+
 @pytest.mark.parametrize("start_year", [1, 1978, 3456])
 @pytest.mark.parametrize("time_step_s", [1, 60, 300])
 def test_write_arrays_absolute(
@@ -98,7 +158,7 @@ def test_write_arrays_absolute(
     temp_dir: tempfile.TemporaryDirectory,
     start_year: int,
     time_step_s: int,
-    maps_dict: Dict,
+    maps_dict: dict,
 ):
     # Write timesteps
     time_steps_num = 3
@@ -145,7 +205,7 @@ def test_write_arrays_relative(
     temp_dir: tempfile.TemporaryDirectory,
     start_seconds: int,
     time_step_s: int,
-    maps_dict: Dict,
+    maps_dict: dict,
 ):
     """Test writing arrays with relative time (timedelta)"""
     # Write timesteps, with 1 minute in between
@@ -184,7 +244,7 @@ def test_data_consistency(
     icechunk_provider: IcechunkRasterOutputProvider,
     temp_dir: tempfile.TemporaryDirectory,
     maps_dict: Mapping[str, np.ndarray],
-    coordinates: Dict,
+    coordinates: dict,
     crs: pyproj.CRS,
     out_map_names: Mapping[str, str],
 ):
@@ -219,7 +279,7 @@ def test_data_consistency(
     # Assert that all expected data variables are present
     expected_var_names = set(out_map_names.values())
     actual_var_names = set(ds.data_vars.keys())
-    assert expected_var_names.issubset(actual_var_names), (
+    assert expected_var_names == actual_var_names, (
         f"Expected {expected_var_names}, actual {actual_var_names}"
     )
 
@@ -240,6 +300,7 @@ def test_data_consistency(
     # Assert that data values are preserved for each variable at both timesteps
     for internal_key, zarr_var_name in out_map_names.items():
         if zarr_var_name in ds.data_vars:
+            assert ds[zarr_var_name].dims == ("time", "y", "x")
             # Check first timestep data
             original_data_1 = maps_dict_1[internal_key]
             actual_data_1 = ds[zarr_var_name].isel(time=0).values
@@ -273,8 +334,8 @@ def test_data_consistency(
 
 def test_non_matching_shape(
     temp_dir: tempfile.TemporaryDirectory,
-    maps_dict: Dict,
-    coordinates: Dict,
+    maps_dict: dict,
+    coordinates: dict,
     crs: pyproj.CRS,
     out_map_names: Mapping[str, str],
 ):
@@ -319,8 +380,8 @@ def test_non_matching_shape(
 
 def test_non_matching_variable_names(
     temp_dir: tempfile.TemporaryDirectory,
-    maps_dict: Dict,
-    coordinates: Dict,
+    maps_dict: dict,
+    coordinates: dict,
     crs: pyproj.CRS,
     out_map_names: Mapping[str, str],
 ):
@@ -359,8 +420,8 @@ def test_non_matching_variable_names(
 
 def test_non_matching_number_of_variables(
     temp_dir: tempfile.TemporaryDirectory,
-    maps_dict: Dict,
-    coordinates: Dict,
+    maps_dict: dict,
+    coordinates: dict,
     crs: pyproj.CRS,
     out_map_names: Mapping[str, str],
 ):
@@ -400,8 +461,8 @@ def test_non_matching_number_of_variables(
 
 def test_non_matching_coordinates_same_dimensions(
     temp_dir: tempfile.TemporaryDirectory,
-    maps_dict: Dict,
-    coordinates: Dict,
+    maps_dict: dict,
+    coordinates: dict,
     crs: pyproj.CRS,
     out_map_names: Mapping[str, str],
 ):
@@ -443,8 +504,8 @@ def test_non_matching_coordinates_same_dimensions(
 
 def test_non_matching_crs(
     temp_dir: tempfile.TemporaryDirectory,
-    maps_dict: Dict,
-    coordinates: Dict,
+    maps_dict: dict,
+    coordinates: dict,
     crs: pyproj.CRS,
     out_map_names: Mapping[str, str],
 ):
@@ -483,8 +544,8 @@ def test_non_matching_crs(
 
 def test_multi_session_data_persistence(
     temp_dir: tempfile.TemporaryDirectory,
-    maps_dict: Dict,
-    coordinates: Dict,
+    maps_dict: dict,
+    coordinates: dict,
     crs: pyproj.CRS,
     out_map_names: Mapping[str, str],
 ):
@@ -605,129 +666,59 @@ def test_multi_session_data_persistence(
     assert crs == crs_actual
 
 
-def test_finalize_max_values(
-    icechunk_provider: IcechunkRasterOutputProvider,
+def test_maxima_use_configured_names_without_base_arrays(
     temp_dir: tempfile.TemporaryDirectory,
-    maps_dict: Dict,
-    coordinates: Dict,
+    coordinates: dict,
     crs: pyproj.CRS,
-    out_map_names: Mapping[str, str],
 ):
-    """Test that finalize() method properly writes max values arrays
-    and does not affect existing data."""
-
-    # Create initial simulation data with regular arrays
-    sim_time_1 = datetime(year=2023, month=1, day=1, hour=10)
-    sim_time_2 = datetime(year=2023, month=1, day=1, hour=11)
-
-    # Write initial timesteps with regular data
-    icechunk_provider.write_arrays(maps_dict, sim_time_1)
-    icechunk_provider.write_arrays(maps_dict, sim_time_2)
-
-    # Create max value arrays (should be different from regular data)
-    rng = np.random.default_rng(seed=999)
-    arr_shape = next(iter(maps_dict.values())).shape
-    max_arrays = {
-        "hmax": rng.random(size=arr_shape, dtype=np.float32) + 10.0,
-        "vmax": rng.random(size=arr_shape, dtype=np.float32) + 5.0,
+    storage = icechunk.local_filesystem_storage(temp_dir.name)
+    out_map_names = {"hmax": "maximum_water_depth", "vmax": "maximum_water_speed"}
+    provider = IcechunkRasterOutputProvider(
+        {
+            "out_map_names": out_map_names,
+            "crs": crs,
+            "x_coords": coordinates["x_coords"],
+            "y_coords": coordinates["y_coords"],
+            "icechunk_storage": storage,
+        }
+    )
+    arrays = {
+        "hmax": np.full((6, 9), 2.0, dtype=np.float32),
+        "vmax": np.full((6, 9), 3.0, dtype=np.float32),
     }
+    provider.write_arrays(arrays, timedelta(seconds=30))
+    provider.finalize()
 
-    # Create a SimulationData object for finalize()
-    final_sim_time = datetime(year=2023, month=1, day=1, hour=12)
-
-    final_data = SimulationData(
-        sim_time=final_sim_time,
-        time_step=3600.0,  # 1 hour in seconds
-        time_steps_counter=1,
-        continuity_data=None,  # Not used in finalize()
-        raw_arrays=max_arrays,
-        accumulation_arrays={},  # Not used in finalize()
-        cell_dx=1.0,
-        cell_dy=1.0,
-        drainage_network_data=None,  # Not used in finalize()
-    )
+    ds = xr.open_zarr(icechunk.Repository.open(storage).readonly_session("main").store)
+    assert set(ds.data_vars) == set(out_map_names.values())
+    for key, name in out_map_names.items():
+        assert ds[name].dims == ("time", "y", "x")
+        np.testing.assert_array_equal(ds[name].isel(time=0).values, arrays[key])
 
-    # Call finalize to write max values
-    icechunk_provider.finalize(final_data)
 
-    # Read the data back and verify
+def test_legacy_static_maximum_schema_is_rejected(
+    temp_dir: tempfile.TemporaryDirectory,
+    coordinates: dict,
+    crs: pyproj.CRS,
+):
     storage = icechunk.local_filesystem_storage(temp_dir.name)
-    repo = icechunk.Repository.open(storage)
-    session = repo.readonly_session("main")
-    ds = xr.open_zarr(session.store)
-    print(ds)
-
-    # Assert that we still have 2 timesteps
-    assert ds.sizes["time"] == 2
-
-    # Verify that original data is still intact
-    for internal_key, var_name in out_map_names.items():
-        if var_name in ds.data_vars:
-            # Check first timestep data (should be unchanged)
-            original_data = maps_dict[internal_key]
-            actual_data_t0 = ds[var_name].isel(time=0).values
-            actual_data_t1 = ds[var_name].isel(time=1).values
-            assert np.allclose(actual_data_t0, original_data), (
-                f"First timestep data was modified for {var_name}"
-            )
-            assert np.allclose(actual_data_t1, original_data), (
-                f"Second timestep data was modified for {var_name}"
-            )
-
-    # Verify that max values are correctly written
-    # Check for hmax
-    if "water_depth" in icechunk_provider.out_map_names:
-        assert "test_water_depth_max" in ds.data_vars, (
-            "test_water_depth_max should be present in dataset"
-        )
-        expected_max_data = max_arrays["hmax"]
-        actual_max_data = ds["test_water_depth_max"].values
-        assert np.allclose(actual_max_data, expected_max_data), (
-            "hmax data does not match expected values"
-        )
-
-        # Ensure max values are different from regular data
-        regular_water_depth = (
-            ds["water_depth"].isel(time=0).values if "water_depth" in ds.data_vars else None
-        )
-        if regular_water_depth is not None:
-            assert not np.allclose(actual_max_data, regular_water_depth), (
-                "Max values should be different from regular data"
-            )
-
-    # Check for vmax
-    if "v" in icechunk_provider.out_map_names:
-        assert "test_v_max" in ds.data_vars, "test_v_max should be present in dataset"
-        expected_max_data = max_arrays["vmax"]
-        actual_max_data = ds["test_v_max"].values
-        assert np.allclose(actual_max_data, expected_max_data), (
-            "vmax data does not match expected values"
+    repo = icechunk.Repository.create(storage)
+    legacy = xr.Dataset(
+        {"maximum_water_depth": (("y", "x"), np.zeros((6, 9), dtype=np.float32))},
+        coords={"x": coordinates["x_coords"], "y": coordinates["y_coords"]},
+        attrs={"crs_wkt": crs.to_wkt()},
+    )
+    session = repo.writable_session("main")
+    icechunk.xarray.to_icechunk(legacy, session, mode="w-")
+    session.commit("legacy maximum")
+
+    with pytest.raises(ValueError, match="dimensions.*new repository.*migrate"):
+        IcechunkRasterOutputProvider(
+            {
+                "out_map_names": {"hmax": "maximum_water_depth"},
+                "crs": crs,
+                "x_coords": coordinates["x_coords"],
+                "y_coords": coordinates["y_coords"],
+                "icechunk_storage": storage,
+            }
         )
-
-        # Ensure max values are different from regular data
-        regular_v = ds["v"].isel(time=0).values if "v" in ds.data_vars else None
-        if regular_v is not None:
-            assert not np.allclose(actual_max_data, regular_v), (
-                "Max values should be different from regular data"
-            )
-
-    # Verify that timestamps are correct
-    expected_times = [sim_time_1, sim_time_2]
-    actual_times = [pd.to_datetime(t).to_pydatetime() for t in ds["time"].values]
-    assert len(actual_times) == len(expected_times)
-    for actual, expected in zip(actual_times, expected_times):
-        assert actual == expected, f"Expected {expected}, got {actual}"
-
-    # Verify that spatial coordinates are preserved
-    assert "x" in ds.coords
-    assert "y" in ds.coords
-    expected_x = coordinates["x_coords"]
-    actual_x = ds.coords["x"].values
-    assert np.allclose(actual_x, expected_x)
-    expected_y = coordinates["y_coords"]
-    actual_y = ds.coords["y"].values
-    assert np.allclose(actual_y, expected_y)
-
-    # Assert that CRS information is preserved
-    crs_actual = pyproj.CRS.from_wkt(ds.attrs["crs_wkt"])
-    assert crs == crs_actual
diff --git a/tests/test_max_values.py b/tests/test_max_values.py
index 77b5606..df0a4fb 100644
--- a/tests/test_max_values.py
+++ b/tests/test_max_values.py
@@ -15,6 +15,7 @@
 from __future__ import annotations
 
 from datetime import datetime, timedelta
+from itertools import pairwise
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -36,7 +37,7 @@
 def sim_5by5_max_values(domain_5by5, helpers) -> Simulation:
     """Run a 5x5 simulation for 2s with 1s record step.
 
-    Outputs: water_depth, v
+    Outputs: water_depth, hmax, v, vmax
     Used for testing that max values are correctly computed.
     """
     # Build SimulationConfig
@@ -52,7 +53,7 @@ def sim_5by5_max_values(domain_5by5, helpers) -> Simulation:
         ),
         output_map_names=helpers.make_output_map_names(
             "out_5by5_max_values",
-            ["water_depth", "v"],
+            ["water_depth", "hmax", "v", "vmax"],
         ),
         # Same values as 5by5_max_values.ini
         surface_flow_parameters=SurfaceFlowParameters(hmin=0.000001, dtmax=1, cfl=0.8),
@@ -88,44 +89,53 @@ def sim_5by5_max_values(domain_5by5, helpers) -> Simulation:
 class TestMaxValues:
     """Test that the maximum values of h and v are properly calculated.
 
-    The simulation tracks hmax and vmax internally. These should equal
-    the element-wise maximum across all time-step outputs.
+    The simulation tracks hmax and vmax internally and reports them as
+    cumulative maximum arrays.
     """
 
     def test_water_depth_max(self, sim_5by5_max_values):
-        """The internal hmax array should match max of all water_depth outputs."""
+        """Reported hmax is nondecreasing and ends at the internal maximum."""
         output_dict = sim_5by5_max_values.report.raster_provider.output_maps_dict
 
-        # Get all water_depth arrays from output
+        h_max_arrays = [arr for _, arr in output_dict["hmax"]]
         h_arrays = [arr for _, arr in output_dict["water_depth"]]
-
-        # Compute element-wise maximum across all time steps
-        h_max_computed = np.maximum.reduce(h_arrays)
-
-        # Get the internal hmax array
         h_max_internal = sim_5by5_max_values.get_array("hmax")
 
-        # The overall max values should match
-        assert np.isclose(np.nanmax(h_max_computed), np.nanmax(h_max_internal)), (
-            f"hmax mismatch: computed max={np.nanmax(h_max_computed):.6f}, "
-            f"internal max={np.nanmax(h_max_internal):.6f}"
-        )
+        assert len(h_max_arrays) == 3
+        np.testing.assert_allclose(h_max_arrays[0], h_arrays[0])
+        for previous, current in pairwise(h_max_arrays):
+            assert np.all(current >= previous)
+        np.testing.assert_allclose(h_max_arrays[-1], h_max_internal)
 
     def test_velocity_max(self, sim_5by5_max_values):
-        """The internal vmax array should match max of all v outputs."""
+        """Reported vmax is nondecreasing and ends at the internal maximum."""
         output_dict = sim_5by5_max_values.report.raster_provider.output_maps_dict
 
-        # Get all v arrays from output
+        v_max_arrays = [arr for _, arr in output_dict["vmax"]]
         v_arrays = [arr for _, arr in output_dict["v"]]
+        v_max_internal = sim_5by5_max_values.get_array("vmax")
 
-        # Compute element-wise maximum across all time steps
-        v_max_computed = np.maximum.reduce(v_arrays)
+        assert len(v_max_arrays) == 3
+        np.testing.assert_allclose(v_max_arrays[0], v_arrays[0])
+        for previous, current in pairwise(v_max_arrays):
+            assert np.all(current >= previous)
+        np.testing.assert_allclose(v_max_arrays[-1], v_max_internal)
 
-        # Get the internal vmax array
-        v_max_internal = sim_5by5_max_values.get_array("vmax")
 
-        # The overall max values should match
-        assert np.isclose(np.nanmax(v_max_computed), np.nanmax(v_max_internal)), (
-            f"vmax mismatch: computed max={np.nanmax(v_max_computed):.6f}, "
-            f"internal max={np.nanmax(v_max_internal):.6f}"
-        )
+def test_set_array_synchronizes_maxima(sim_5by5_max_values: Simulation) -> None:
+    simulation = sim_5by5_max_values
+
+    larger_depth = simulation.get_array("hmax").copy() + 1.0
+    simulation.set_array("water_depth", larger_depth)
+    np.testing.assert_allclose(simulation.get_array("hmax"), larger_depth)
+
+    simulation.set_array("water_depth", np.zeros_like(larger_depth))
+    np.testing.assert_allclose(simulation.get_array("hmax"), larger_depth)
+
+    larger_wse = simulation.get_array("dem") + larger_depth + 1.0
+    simulation.set_array("water_surface_elevation", larger_wse)
+    np.testing.assert_allclose(simulation.get_array("hmax"), larger_depth + 1.0)
+
+    larger_speed = simulation.get_array("vmax").copy() + 1.0
+    simulation.set_array("v", larger_speed)
+    np.testing.assert_allclose(simulation.get_array("vmax"), larger_speed)
diff --git a/tests/test_report.py b/tests/test_report.py
new file mode 100644
index 0000000..37fc995
--- /dev/null
+++ b/tests/test_report.py
@@ -0,0 +1,102 @@
+"""
+Copyright (C) 2026 Laurent G. Courty
+
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public License
+as published by the Free Software Foundation; either version 2.1
+of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Lesser General Public License for more details.
+"""
+
+from datetime import UTC, datetime, timedelta
+
+import numpy as np
+
+from itzi_core.const import TemporalType
+from itzi_core.data_containers import SimulationData
+from itzi_core.providers.memory_output import (
+    MemoryRasterOutputProvider,
+    MemoryVectorOutputProvider,
+)
+from itzi_core.report import Report
+
+
+def test_get_output_arrays_returns_a_fresh_selection() -> None:
+    start_time = datetime(2000, 1, 1, tzinfo=UTC)
+    out_map_names = {"water_depth": "depth", "hmax": "depth_max"}
+    raster_provider = MemoryRasterOutputProvider(out_map_names)
+    report = Report(
+        start_time=start_time,
+        temporal_type=TemporalType.ABSOLUTE,
+        raster_output_provider=raster_provider,
+        vector_output_provider=MemoryVectorOutputProvider(),
+        mass_balance_output_provider=None,
+        out_map_names=out_map_names,
+        dt=timedelta(seconds=1),
+    )
+    data = SimulationData(
+        sim_time=start_time,
+        time_step=1.0,
+        time_steps_counter=0,
+        continuity_data=None,
+        raw_arrays={
+            "water_depth": np.array([[1.0]], dtype=np.float32),
+            "hmax": np.array([[2.0]], dtype=np.float32),
+        },
+        accumulation_arrays={},
+        cell_dx=1.0,
+        cell_dy=1.0,
+        drainage_network_data=None,
+    )
+
+    report.step(data)
+    del out_map_names["water_depth"]
+    report.step(data)
+
+    assert len(raster_provider.output_maps_dict["water_depth"]) == 1
+    assert len(raster_provider.output_maps_dict["hmax"]) == 2
+
+
+def test_maxima_are_selected_independently_of_base_arrays() -> None:
+    start_time = datetime(2000, 1, 1, tzinfo=UTC)
+    out_map_names = {"hmax": "depth_max", "vmax": "speed_max"}
+    raster_provider = MemoryRasterOutputProvider(out_map_names)
+    report = Report(
+        start_time=start_time,
+        temporal_type=TemporalType.ABSOLUTE,
+        raster_output_provider=raster_provider,
+        vector_output_provider=MemoryVectorOutputProvider(),
+        mass_balance_output_provider=None,
+        out_map_names=out_map_names,
+        dt=timedelta(seconds=1),
+    )
+    data = SimulationData(
+        sim_time=start_time,
+        time_step=1.0,
+        time_steps_counter=0,
+        continuity_data=None,
+        raw_arrays={
+            "water_depth": np.array([[1.0]], dtype=np.float32),
+            "hmax": np.array([[2.0]], dtype=np.float32),
+            "v": np.array([[3.0]], dtype=np.float32),
+            "vmax": np.array([[4.0]], dtype=np.float32),
+        },
+        accumulation_arrays={},
+        cell_dx=1.0,
+        cell_dy=1.0,
+        drainage_network_data=None,
+    )
+
+    report.step(data)
+
+    assert set(raster_provider.output_maps_dict) == {"hmax", "vmax"}
+    np.testing.assert_array_equal(
+        raster_provider.output_maps_dict["hmax"][0][1], data.raw_arrays["hmax"]
+    )
+    np.testing.assert_array_equal(
+        raster_provider.output_maps_dict["vmax"][0][1], data.raw_arrays["vmax"]
+    )
diff --git a/tests/test_wse.py b/tests/test_wse.py
index 09667ff..3b93cc6 100644
--- a/tests/test_wse.py
+++ b/tests/test_wse.py
@@ -153,7 +153,10 @@ def test_timed_memory_input_updates_water_depth_from_wse(domain_5by5) -> None:
             "friction": "friction",
             "water_surface_elevation": "water_surface_elevation",
         },
-        output_map_names={"water_depth": "out_5by5_wse_timed_memory_water_depth"},
+        output_map_names={
+            "water_depth": "out_5by5_wse_timed_memory_water_depth",
+            "hmax": "out_5by5_wse_timed_memory_hmax",
+        },
         surface_flow_parameters=SurfaceFlowParameters(hmin=0.0001, dtmax=0.3, cfl=0.2),
         infiltration_model=InfiltrationModelType.NULL,
     )
@@ -182,10 +185,11 @@ def test_timed_memory_input_updates_water_depth_from_wse(domain_5by5) -> None:
             },
         }
     )
+    raster_output = MemoryRasterOutputProvider(sim_config.output_map_names)
     simulation = (
         SimulationBuilder(sim_config, domain_5by5.arr_mask, np.float32)
         .with_input_provider(input_provider)
-        .with_raster_output_provider(MemoryRasterOutputProvider(sim_config.output_map_names))
+        .with_raster_output_provider(raster_output)
         .with_vector_output_provider(MemoryVectorOutputProvider())
         .build()
     )
@@ -197,6 +201,10 @@ def test_timed_memory_input_updates_water_depth_from_wse(domain_5by5) -> None:
         np.full(domain_5by5.domain_data.shape, 0.2, dtype=np.float32),
         atol=1e-5,
     )
+    np.testing.assert_allclose(
+        simulation.raster_domain.get_array("hmax"),
+        simulation.raster_domain.get_array("water_depth"),
+    )
     assert simulation.timed_arrays is not None
     wse_timed_array = simulation.timed_arrays["water_surface_elevation"]
     assert wse_timed_array.arr_start == start_time
@@ -209,5 +217,11 @@ def test_timed_memory_input_updates_water_depth_from_wse(domain_5by5) -> None:
         np.full(domain_5by5.domain_data.shape, 0.35, dtype=np.float32),
         atol=1e-5,
     )
+    np.testing.assert_allclose(
+        simulation.raster_domain.get_array("hmax"),
+        simulation.raster_domain.get_array("water_depth"),
+    )
+    reported_hmax = raster_output.output_maps_dict["hmax"][-1][1]
+    np.testing.assert_allclose(reported_hmax, simulation.raster_domain.get_array("hmax"))
     assert wse_timed_array.arr_start == boundary_time
     assert wse_timed_array.arr_end == end_time