diff --git a/imod/tests/fixtures/flow_transport_simulation_fixture.py b/imod/tests/fixtures/flow_transport_simulation_fixture.py
index 3f64cf97f..9b40812bd 100644
--- a/imod/tests/fixtures/flow_transport_simulation_fixture.py
+++ b/imod/tests/fixtures/flow_transport_simulation_fixture.py
@@ -207,10 +207,10 @@ def flow_transport_simulation():
         relaxation_factor=0.9,
     )
 
-    duration = pd.to_timedelta("2000d")
+    duration = pd.to_timedelta("20d")
     start = pd.to_datetime("2000-01-01")
     simulation.create_time_discretization(additional_times=[start, start + duration])
-    simulation["time_discretization"]["n_timesteps"] = 100
+    simulation["time_discretization"]["n_timesteps"] = 10
 
     return simulation
     # %%
diff --git a/imod/tests/fixtures/mf6_circle_fixture.py b/imod/tests/fixtures/mf6_circle_fixture.py
index 42a0f1e5d..eeaa6a856 100644
--- a/imod/tests/fixtures/mf6_circle_fixture.py
+++ b/imod/tests/fixtures/mf6_circle_fixture.py
@@ -75,7 +75,7 @@ def make_circle_model():
         reordering_method=None,
         relaxation_factor=0.97,
     )
-    simtimes = pd.date_range(start="2000-01-01", end="2001-01-01", freq="W")
+    simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
     simulation.create_time_discretization(additional_times=simtimes)
     return simulation
 
@@ -158,7 +158,7 @@ def make_circle_model_flow_with_transport_data(species: list[str]):
         reordering_method=None,
         relaxation_factor=0.97,
     )
-    simtimes = pd.date_range(start="2000-01-01", end="2001-01-01", freq="W")
+    simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
     simulation.create_time_discretization(simtimes)
     return simulation
 
@@ -357,7 +357,7 @@ def circle_model_transport():
         reordering_method=None,
         relaxation_factor=0.97,
     )
-    simtimes = pd.date_range(start="2000-01-01", end="2001-01-01", freq="W")
+    simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
     simulation.create_time_discretization(additional_times=simtimes)
     return simulation
 
@@ -437,7 +437,7 @@ def circle_model_transport_vsc():
         reordering_method=None,
         relaxation_factor=0.97,
     )
-    simtimes = pd.date_range(start="2000-01-01", end="2001-01-01", freq="W")
+    simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
     simulation.create_time_discretization(additional_times=simtimes)
     return simulation
 
@@ -514,6 +514,6 @@ def circle_model_transport_multispecies_variable_density():
         reordering_method=None,
         relaxation_factor=0.97,
     )
-    simtimes = pd.date_range(start="2000-01-01", end="2001-01-01", freq="W")
+    simtimes = pd.date_range(start="2000-01-01", end="2000-01-03")
     simulation.create_time_discretization(additional_times=simtimes)
     return simulation
diff --git a/imod/tests/fixtures/mf6_twri_fixture.py b/imod/tests/fixtures/mf6_twri_fixture.py
index 717282173..279bae2b3 100644
--- a/imod/tests/fixtures/mf6_twri_fixture.py
+++ b/imod/tests/fixtures/mf6_twri_fixture.py
@@ -232,7 +232,7 @@ def transient_twri_model():
         save_flows=True,
     )
     simulation.create_time_discretization(
-        additional_times=pd.date_range("2000-01-01", " 2000-01-31")
+        additional_times=pd.date_range("2000-01-01", " 2000-01-03")
     )
     return simulation
 
@@ -282,7 +282,7 @@ def transient_unconfined_twri_model():
     # Write specific discharges
     gwf_model["npf"]["save_specific_discharge"] = True
     simulation.create_time_discretization(
-        additional_times=pd.date_range("2000-01-01", " 2000-01-31")
+        additional_times=pd.date_range("2000-01-01", " 2000-01-03")
     )
     return simulation
 
diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py
index 4c761214c..03883f5d1 100644
--- a/imod/tests/test_mf6/test_circle.py
+++ b/imod/tests/test_mf6/test_circle.py
@@ -3,6 +3,7 @@
 from datetime import datetime
 
 import numpy as np
+import pandas as pd
 import pytest
 import xugrid as xu
 
@@ -24,17 +25,18 @@ def test_simulation_write_and_run(circle_model, tmp_path):
     with pytest.raises(
         RuntimeError, match="Simulation circle has not been written yet."
     ):
-        circle_model.run()
+        simulation.run()
 
     modeldir = tmp_path / "circle"
 
-    mask = deepcopy(circle_model["GWF_1"].domain)
+    mask = deepcopy(simulation["GWF_1"].domain)
     mask.values[:, 133] = -1
     simulation.mask_all_models(mask)
 
     simulation.write(modeldir, binary=False, use_absolute_paths=True)
     simulation.run()
-
+    # Test opening with different time unit, see if time unit is correctly
+    # converted to seconds
     head = imod.mf6.open_hds(
         modeldir / "GWF_1/GWF_1.hds",
         modeldir / "GWF_1/disv.disv.grb",
@@ -44,10 +46,10 @@ def test_simulation_write_and_run(circle_model, tmp_path):
     assert isinstance(head, xu.UgridDataArray)
     assert np.all(
         head["time"].values[0]
-        == np.array(datetime(1999, 1, 1, 0, 0, 7), dtype="datetime64[ns]")
+        == np.array(datetime(1999, 1, 1, 0, 0, 1), dtype="datetime64[ns]")
     )
     assert head.dims == ("time", "layer", "mesh2d_nFaces")
-    assert head.shape == (52, 2, 216)
+    assert head.shape == (2, 2, 216)
 
 
 def test_gwfmodel_render(circle_model, tmp_path):
@@ -87,7 +89,7 @@ def test_simulation_write_and_run_evt(circle_model_evt, tmp_path):
     with pytest.raises(
         RuntimeError, match="Simulation circle has not been written yet."
     ):
-        circle_model_evt.run()
+        simulation.run()
 
     modeldir = tmp_path / "circle_evt"
     simulation.write(modeldir, binary=True)
@@ -98,22 +100,22 @@ def test_simulation_write_and_run_evt(circle_model_evt, tmp_path):
     )
     assert isinstance(head, xu.UgridDataArray)
     assert head.dims == ("time", "layer", "mesh2d_nFaces")
-    assert head.shape == (52, 2, 216)
+    assert head.shape == (2, 2, 216)
 
 
 def test_simulation_write_and_run_evt__no_segments(circle_model_evt, tmp_path):
     simulation = circle_model_evt
 
-    evt = circle_model_evt["GWF_1"].pop("evt")
+    evt = simulation["GWF_1"].pop("evt")
     evt_ds = evt.dataset
     evt_ds = evt_ds.drop_vars(["proportion_rate", "proportion_depth"])
     evt_dict = {key: evt_ds[key] for key in evt_ds.keys()}
-    circle_model_evt["GWF_1"]["evt"] = imod.mf6.Evapotranspiration(**evt_dict)
+    simulation["GWF_1"]["evt"] = imod.mf6.Evapotranspiration(**evt_dict)
 
     with pytest.raises(
         RuntimeError, match="Simulation circle has not been written yet."
     ):
-        circle_model_evt.run()
+        simulation.run()
 
     modeldir = tmp_path / "circle_evt"
     simulation.write(modeldir, binary=True)
@@ -124,7 +126,7 @@ def test_simulation_write_and_run_evt__no_segments(circle_model_evt, tmp_path):
     )
     assert isinstance(head, xu.UgridDataArray)
     assert head.dims == ("time", "layer", "mesh2d_nFaces")
-    assert head.shape == (52, 2, 216)
+    assert head.shape == (2, 2, 216)
 
 
 def test_gwfmodel_render_evt(circle_model_evt, tmp_path):
@@ -161,11 +163,10 @@ def test_gwfmodel_render_evt(circle_model_evt, tmp_path):
 
 def test_simulation_write_and_run_transport(circle_model_transport, tmp_path):
     simulation = circle_model_transport
-
     with pytest.raises(
         RuntimeError, match="Simulation circle has not been written yet."
     ):
-        circle_model_transport.run()
+        simulation.run()
 
     modeldir = tmp_path / "circle_transport"
     simulation.write(modeldir)
@@ -175,11 +176,11 @@ def test_simulation_write_and_run_transport(circle_model_transport, tmp_path):
 
     assert isinstance(head, xu.UgridDataArray)
     assert head.dims == ("time", "layer", "mesh2d_nFaces")
-    assert head.shape == (52, 2, 216)
+    assert head.shape == (2, 2, 216)
 
     assert isinstance(concentration, xu.UgridDataArray)
     assert concentration.dims == ("time", "layer", "mesh2d_nFaces")
-    assert concentration.shape == (52, 2, 216)
+    assert concentration.shape == (2, 2, 216)
 
 
 def test_simulation_write_and_run_transport_vsc(circle_model_transport_vsc, tmp_path):
@@ -198,11 +199,11 @@ def test_simulation_write_and_run_transport_vsc(circle_model_transport_vsc, tmp_
 
     assert isinstance(head, xu.UgridDataArray)
     assert head.dims == ("time", "layer", "mesh2d_nFaces")
-    assert head.shape == (52, 2, 216)
+    assert head.shape == (2, 2, 216)
 
     assert isinstance(concentration, xu.UgridDataArray)
     assert concentration.dims == ("time", "layer", "mesh2d_nFaces")
-    assert concentration.shape == (52, 2, 216)
+    assert concentration.shape == (2, 2, 216)
 
 
 def run_simulation(simulation, modeldir):
@@ -251,8 +252,8 @@ def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path)
 
     # Test if model runs
     head_half, concentration_half = run_simulation(half_simulation, tmp_path / "half")
-    assert head_half.shape == (52, 2, 108)
-    assert concentration_half.shape == (52, 2, 108)
+    assert head_half.shape == (2, 2, 108)
+    assert concentration_half.shape == (2, 2, 108)
 
 
 def test_simulation_clip_and_constant_state_at_boundary__transient_chd(
@@ -289,12 +290,12 @@ def test_simulation_clip_and_constant_state_at_boundary__transient_chd(
     assert conc_clipped.notnull().sum() == 20
     head_clipped = half_simulation["GWF_1"]["chd_clipped"]["head"]
     assert head_clipped.sizes["layer"] == 2
-    assert head_clipped.sizes["time"] == 52
-    assert head_clipped.notnull().sum() == 20 * 52
+    assert head_clipped.sizes["time"] == 2
+    assert head_clipped.notnull().sum() == 20 * 2
     # Test if model runs
     head_half, concentration_half = run_simulation(half_simulation, tmp_path / "half")
-    assert head_half.shape == (52, 2, 108)
-    assert concentration_half.shape == (52, 2, 108)
+    assert head_half.shape == (2, 2, 108)
+    assert concentration_half.shape == (2, 2, 108)
 
 
 def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
@@ -304,9 +305,17 @@ def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
     Test with a transient chd boundary condition and transient states for
     boundary conditions. The time steps of the chd package are outside of the
     time domain of the states_for_boundary.
+
+    Extend the time discretization as it was before
+    https://github.com/Deltares/imod-python/pull/1029, to allow clipping longer
+    timeseries.
     """
     # Arrange
     simulation = circle_model_transport
+    simulation.create_time_discretization(
+        additional_times=pd.date_range(start="2000-01-01", end="2001-01-01", freq="W")
+    )
+
     idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
     time = simulation["time_discretization"].dataset.coords["time"].values
 
diff --git a/imod/tests/test_mf6/test_ex01_twri.py b/imod/tests/test_mf6/test_ex01_twri.py
index 5709b8df5..71b24c480 100644
--- a/imod/tests/test_mf6/test_ex01_twri.py
+++ b/imod/tests/test_mf6/test_ex01_twri.py
@@ -4,6 +4,7 @@
 from pathlib import Path
 
 import numpy as np
+import pandas as pd
 import pytest
 import xarray as xr
 
@@ -15,8 +16,32 @@
 from imod.typing.grid import ones_like
 
 
-def test_dis_render(twri_model, tmp_path):
+@pytest.fixture(scope="function")
+def twri_model_extended(twri_model):
+    """Extend the time discretization as it was before
+    https://github.com/Deltares/imod-python/pull/1029, to make sure that
+    integration tests succeed here."""
     simulation = twri_model
+    simulation.create_time_discretization(
+        additional_times=pd.date_range("2000-01-01", " 2000-01-31")
+    )
+    return simulation
+
+
+@pytest.fixture(scope="function")
+def transient_twri_model_extended(transient_twri_model):
+    """Extend the time discretization as it was before
+    https://github.com/Deltares/imod-python/pull/1029, to make sure that
+    integration tests succeed here."""
+    simulation = transient_twri_model
+    simulation.create_time_discretization(
+        additional_times=pd.date_range("2000-01-01", " 2000-01-31")
+    )
+    return simulation
+
+
+def test_dis_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     dis = simulation["GWF_1"]["dis"]
     actual = dis._render(
         directory=tmp_path,
@@ -63,8 +88,8 @@ def test_dis_render(twri_model, tmp_path):
     assert (tmp_path / "dis" / "idomain.bin").is_file()
 
 
-def test_chd_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_chd_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     globaltimes = simulation["time_discretization"]["time"].values
     chd = simulation["GWF_1"]["chd"]
     actual = chd._render(
@@ -99,8 +124,8 @@ def test_chd_render(twri_model, tmp_path):
     assert (tmp_path / "chd" / "chd.bin").is_file()
 
 
-def test_drn_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_drn_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     globaltimes = simulation["time_discretization"]["time"].values
     drn = simulation["GWF_1"]["drn"]
     actual = drn._render(
@@ -135,8 +160,8 @@ def test_drn_render(twri_model, tmp_path):
     assert (tmp_path / "drn" / "drn.bin").is_file()
 
 
-def test_ic_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_ic_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     ic = simulation["GWF_1"]["ic"]
     actual = ic._render(
         directory=tmp_path,
@@ -161,8 +186,8 @@ def test_ic_render(twri_model, tmp_path):
     assert (tmp_path / "ic.ic").is_file()
 
 
-def test_npf_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_npf_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     npf = simulation["GWF_1"]["npf"]
     actual = npf._render(
         directory=tmp_path, pkgname="npf", globaltimes=None, binary=True
@@ -197,8 +222,8 @@ def test_npf_render(twri_model, tmp_path):
     assert (tmp_path / "npf.npf").is_file()
 
 
-def test_oc_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_oc_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     globaltimes = simulation["time_discretization"]["time"].values
     oc = simulation["GWF_1"]["oc"]
 
@@ -242,8 +267,8 @@ def test_oc_render(twri_model, tmp_path):
     assert actual == expected
 
 
-def test_rch_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_rch_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     globaltimes = simulation["time_discretization"]["time"].values
     rch = simulation["GWF_1"]["rch"]
     actual = rch._render(
@@ -275,8 +300,8 @@ def test_rch_render(twri_model, tmp_path):
     assert (tmp_path / "rch" / "rch.bin").is_file()
 
 
-def test_wel_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_wel_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     simulation.write(
         tmp_path,
     )
@@ -307,8 +332,8 @@ def test_wel_render(twri_model, tmp_path):
     assert (tmp_path / "GWF_1/wel" / "wel.bin").is_file()
 
 
-def test_solver_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_solver_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     solver = simulation["solver"]
     actual = solver._render(
         directory=tmp_path,
@@ -342,8 +367,8 @@ def test_solver_render(twri_model, tmp_path):
     assert (tmp_path / "solver.ims").is_file()
 
 
-def test_gwfmodel_render(twri_model, tmp_path):
-    simulation = twri_model
+def test_gwfmodel_render(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     globaltimes = simulation["time_discretization"]["time"].values
     gwfmodel = simulation["GWF_1"]
     path = Path(tmp_path.stem).as_posix()
@@ -375,8 +400,8 @@ def test_gwfmodel_render(twri_model, tmp_path):
     assert (tmp_path / "GWF_1").is_dir()
 
 
-def test_simulation_render(twri_model):
-    simulation = twri_model
+def test_simulation_render(twri_model_extended):
+    simulation = twri_model_extended
     write_context = WriteContext(".")
     actual = simulation._render(write_context)
 
@@ -406,13 +431,13 @@ def test_simulation_render(twri_model):
 
 
 @pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7")
-def test_simulation_write_and_run(twri_model, tmp_path):
-    simulation = twri_model
+def test_simulation_write_and_run(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
 
     with pytest.raises(
         RuntimeError, match="Simulation ex01-twri has not been written yet."
     ):
-        twri_model.run()
+        simulation.run()
 
     modeldir = tmp_path / "ex01-twri"
     simulation.write(
@@ -429,10 +454,10 @@ def test_simulation_write_and_run(twri_model, tmp_path):
     )
     assert isinstance(head, xr.DataArray)
     assert head.dims == ("time", "layer", "y", "x")
-    assert head.shape == (1, 3, 15, 15)
+    assert head.shape == (30, 3, 15, 15)
     assert np.all(
-        head["time"].values
-        == np.array("1999-01-02T00:00:00.000000", dtype="datetime64[ns]")
+        head["time"].values[-1]
+        == np.array("1999-01-31T00:00:00.000000000", dtype="datetime64[ns]")
     )
     meanhead_layer = head.groupby("layer").mean(dim=xr.ALL_DIMS)
     mean_answer = np.array([59.79181509, 30.44132373, 24.88576811])
@@ -440,8 +465,8 @@ def test_simulation_write_and_run(twri_model, tmp_path):
 
 
 @pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7")
-def test_simulation_write_storage(transient_twri_model, tmp_path):
-    simulation = transient_twri_model
+def test_simulation_write_storage(transient_twri_model_extended, tmp_path):
+    simulation = transient_twri_model_extended
     modeldir = tmp_path / "ex01-twri-transient"
     simulation.write(modeldir, binary=True)
     simulation.run()
@@ -492,8 +517,8 @@ def test_simulation_validate_false(transient_twri_model, tmp_path):
 
 
 @pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7")
-def test_simulation_write_errors(twri_model, tmp_path):
-    simulation = twri_model
+def test_simulation_write_errors(twri_model_extended, tmp_path):
+    simulation = twri_model_extended
     model = simulation["GWF_1"]
     model.pop("sto")
     modeldir = tmp_path / "ex01-twri"
@@ -523,16 +548,17 @@ def test_slice_and_run(transient_twri_model, tmp_path):
 
 
 @pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7")
-def test_slice_and_run_with_state(transient_twri_model, tmp_path):
+def test_slice_and_run_with_state(transient_twri_model_extended, tmp_path):
+    simulation = transient_twri_model_extended
     # TODO: bring back well once slicing is implemented...
-    transient_twri_model["GWF_1"].pop("wel")
-    transient_twri_model.write(tmp_path, False, True, False)
-    transient_twri_model.run()
-    head = transient_twri_model.open_head(simulation_start_time="2000-01-01")
+    simulation["GWF_1"].pop("wel")
+    simulation.write(tmp_path, False, True, False)
+    simulation.run()
+    head = simulation.open_head(simulation_start_time="2000-01-01")
 
     # set the boundary to  recognizable value, in this case 1.23
     state_for_boundary = {"GWF_1": 1.23 * ones_like(head)}
-    clipped_simulation = transient_twri_model.clip_box(
+    clipped_simulation = simulation.clip_box(
         time_min="2000-01-10",
         time_max="2000-01-20",
         layer_min=1,
diff --git a/imod/tests/test_mf6/test_mf6_out.py b/imod/tests/test_mf6/test_mf6_out.py
index e3b412b10..c61233e6e 100644
--- a/imod/tests/test_mf6/test_mf6_out.py
+++ b/imod/tests/test_mf6/test_mf6_out.py
@@ -415,7 +415,7 @@ def test_open_cbc__dis_transient(transient_twri_result):
             "wel_wel",
         ]
         for array in cbc.values():
-            assert array.shape == (30, 3, 15, 15)
+            assert array.shape == (2, 3, 15, 15)
             assert isinstance(array, xr.DataArray)
             assert isinstance(array.data, dask.array.Array)
 
@@ -457,7 +457,7 @@ def test_open_cbc__dis_transient_unconfined(transient_unconfined_twri_result):
             "wel_wel",
         ]
         for array in cbc.values():
-            assert array.shape == (30, 3, 15, 15)
+            assert array.shape == (2, 3, 15, 15)
             assert isinstance(array, xr.DataArray)
             assert isinstance(array.data, dask.array.Array)
 
@@ -479,10 +479,10 @@ def test_open_cbc__disv(circle_result):
         ]
         for key, array in cbc.items():
             if key in ("chd_chd", "flow-lower-face"):
-                assert array.shape == (52, 2, 216)
+                assert array.shape == (2, 2, 216)
                 assert array.dims[-1] == array.ugrid.grid.face_dimension
             else:
-                assert array.shape == (52, 2, 342)
+                assert array.shape == (2, 2, 342)
                 assert array.dims[-1] == array.ugrid.grid.edge_dimension
             assert isinstance(array, xu.UgridDataArray)
             assert isinstance(array.data, dask.array.Array)
@@ -505,10 +505,10 @@ def test_open_cbc__disv_offset(circle_result__offset_origins):
         ]
         for key, array in cbc.items():
             if key in ("chd_chd", "flow-lower-face"):
-                assert array.shape == (52, 2, 216)
+                assert array.shape == (2, 2, 216)
                 assert array.dims[-1] == array.ugrid.grid.face_dimension
             else:
-                assert array.shape == (52, 2, 342)
+                assert array.shape == (2, 2, 342)
                 assert array.dims[-1] == array.ugrid.grid.edge_dimension
             assert isinstance(array, xu.UgridDataArray)
             assert isinstance(array.data, dask.array.Array)
@@ -554,10 +554,10 @@ def test_open_cbc__disv_sto(circle_result_sto):
         ]
         for key, array in cbc.items():
             if key in ("chd_chd", "flow-lower-face", "sto-ss"):
-                assert array.shape == (52, 2, 216)
+                assert array.shape == (2, 2, 216)
                 assert array.dims[-1] == array.ugrid.grid.face_dimension
             else:
-                assert array.shape == (52, 2, 342)
+                assert array.shape == (2, 2, 342)
                 assert array.dims[-1] == array.ugrid.grid.edge_dimension
             assert isinstance(array, xu.UgridDataArray)
             assert isinstance(array.data, dask.array.Array)
diff --git a/imod/tests/test_mf6/test_mf6_regrid_transport.py b/imod/tests/test_mf6/test_mf6_regrid_transport.py
index 1f816fc64..aa5405eb3 100644
--- a/imod/tests/test_mf6/test_mf6_regrid_transport.py
+++ b/imod/tests/test_mf6/test_mf6_regrid_transport.py
@@ -87,11 +87,8 @@ def test_regrid_transport(
     # Compute and compare the total mass of the original model and the regridded mode.
     # Constants like porosity are not included.
     for species in ["species_a", "species_b", "species_c", "species_d"]:
-        original_mass = np.sum(
-            conc.sel(time=2000, species=species).values * cell_volume
-        )
+        original_mass = np.sum(conc.sel(time=20, species=species).values * cell_volume)
         regridded_mass = np.sum(
-            regridded_conc.sel(time=2000, species=species).values
-            * regridded_cell_volume
+            regridded_conc.sel(time=20, species=species).values * regridded_cell_volume
         )
         assert abs((original_mass - regridded_mass)) / original_mass < 3e-2
diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py
index 2fc5f42c9..888daaa0e 100644
--- a/imod/tests/test_mf6/test_mf6_simulation.py
+++ b/imod/tests/test_mf6/test_mf6_simulation.py
@@ -199,19 +199,21 @@ def test_simulation_open_head(circle_model, tmp_path):
     modeldir = tmp_path / "circle"
     simulation.write(modeldir)
     simulation.run()
-    head = simulation.open_head()
 
-    assert isinstance(head, xu.UgridDataArray)
-    assert head.dims == ("time", "layer", "mesh2d_nFaces")
-    assert head.shape == (52, 2, 216)
+    # open heads without time conversion
+    head_notime = simulation.open_head()
+
+    assert isinstance(head_notime, xu.UgridDataArray)
+    assert head_notime.dims == ("time", "layer", "mesh2d_nFaces")
+    assert head_notime.shape == (2, 2, 216)
 
     # open heads with time conversion.
     head = simulation.open_head(
         simulation_start_time=datetime(2013, 3, 11, 22, 0, 0), time_unit="w"
     )
     assert head.dims == ("time", "layer", "mesh2d_nFaces")
-    assert head.shape == (52, 2, 216)
-    assert str(head.coords["time"].values[()][0]) == "2013-04-29T22:00:00.000000"
+    assert head.shape == (2, 2, 216)
+    assert str(head.coords["time"].values[()][0]) == "2013-03-18T22:00:00.000000"
 
 
 class PathCases:
@@ -235,9 +237,10 @@ def test_simulation_write_run_open__different_paths(circle_model, tmp_path, path
         simulation.write(path)
         simulation.run()
         head = simulation.open_head()
-        # Assert not an empty array is returned
+
         assert isinstance(head, xu.UgridDataArray)
-        assert head.shape == (52, 2, 216)
+        assert head.dims == ("time", "layer", "mesh2d_nFaces")
+        assert head.shape == (2, 2, 216)
 
 
 def test_simulation_open_flow_budget(circle_model, tmp_path):
diff --git a/imod/tests/test_mf6/test_multimodel/test_mf6_modelsplitter_transport.py b/imod/tests/test_mf6/test_multimodel/test_mf6_modelsplitter_transport.py
index c834c5c81..83f7cea8b 100644
--- a/imod/tests/test_mf6/test_multimodel/test_mf6_modelsplitter_transport.py
+++ b/imod/tests/test_mf6/test_multimodel/test_mf6_modelsplitter_transport.py
@@ -133,28 +133,28 @@ def test_split_flow_and_transport_model_evaluate_output(
     submodel_labels = submodel_labels.sel(layer=0, drop=True)
 
     # for reference run the original model and load the results
-    simulation.write(tmp_path / "original", binary=False)
+    simulation.write(tmp_path / "original")
     simulation.run()
     original_conc = simulation.open_concentration()
     original_head = simulation.open_head()
 
     # split the model , run the split model and load the results
     new_simulation = simulation.split(submodel_labels)
-    new_simulation.write(tmp_path, binary=False)
+    new_simulation.write(tmp_path)
     new_simulation.run()
     conc = new_simulation.open_concentration()
     head = new_simulation.open_head()
 
     # Compare
     np.testing.assert_allclose(
-        head.sel(time=2000)["head"].values,
-        original_head.sel(time=2000).values,
+        head.sel(time=20)["head"].values,
+        original_head.sel(time=20).values,
         rtol=1e-4,
         atol=1e-6,
     )
     np.testing.assert_allclose(
-        conc.sel(time=2000)["concentration"].values,
-        original_conc.sel(time=2000).values,
+        conc.sel(time=20)["concentration"].values,
+        original_conc.sel(time=20).values,
         rtol=1e-4,
         atol=1e-6,
     )
@@ -191,8 +191,8 @@ def test_split_flow_and_transport_model_evaluate_output_with_species(
 
     # Compare
     np.testing.assert_allclose(
-        conc.sel(time=2000)["concentration"].values,
-        original_conc.sel(time=2000).values,
+        conc.sel(time=20)["concentration"].values,
+        original_conc.sel(time=20).values,
         rtol=1e-4,
         atol=1e-6,
     )
diff --git a/imod/tests/test_mf6/test_multimodel/test_mf6_partitioned_simulation_postprocessing.py b/imod/tests/test_mf6/test_multimodel/test_mf6_partitioned_simulation_postprocessing.py
index 82cac94fe..bb94caeda 100644
--- a/imod/tests/test_mf6/test_multimodel/test_mf6_partitioned_simulation_postprocessing.py
+++ b/imod/tests/test_mf6/test_multimodel/test_mf6_partitioned_simulation_postprocessing.py
@@ -60,38 +60,7 @@ def test_import_heads_structured(
     assert np.allclose(merged_heads.coords["layer"].values, [1, 2, 3])
     assert np.allclose(
         merged_heads.coords["time"].values,
-        [
-            1.0,
-            2.0,
-            3.0,
-            4.0,
-            5.0,
-            6.0,
-            7.0,
-            8.0,
-            9.0,
-            10.0,
-            11.0,
-            12.0,
-            13.0,
-            14.0,
-            15.0,
-            16.0,
-            17.0,
-            18.0,
-            19.0,
-            20.0,
-            21.0,
-            22.0,
-            23.0,
-            24.0,
-            25.0,
-            26.0,
-            27.0,
-            28.0,
-            29.0,
-            30.0,
-        ],
+        [1.0, 2.0],
     )
 
 
@@ -106,9 +75,7 @@ def test_import_heads_unstructured(tmp_path, circle_partitioned):
 
     # Assert
     assert np.allclose(merged_heads.coords["layer"].values, [1, 2])
-    assert np.allclose(
-        merged_heads.coords["time"].values, list(np.arange(7.0, 365.0, 7.0))
-    )
+    assert np.allclose(merged_heads.coords["time"].values, [1.0, 2.0])
     assert np.allclose(merged_heads.coords["mesh2d_nFaces"].values, list(range(216)))
 
 
diff --git a/imod/tests/test_mf6/test_multimodel/test_mf6_partitioning_unstructured.py b/imod/tests/test_mf6/test_multimodel/test_mf6_partitioning_unstructured.py
index 4edefa4d3..88eea96c2 100644
--- a/imod/tests/test_mf6/test_multimodel/test_mf6_partitioning_unstructured.py
+++ b/imod/tests/test_mf6/test_multimodel/test_mf6_partitioning_unstructured.py
@@ -444,14 +444,14 @@ def run_simulation(tmp_path, simulation, species=None):
     simulation.run()
     head = simulation.open_head()
     flow_budget = simulation.open_flow_budget()
-    flow_budget = flow_budget.sel(time=364)
+    flow_budget = flow_budget.isel(time=-1)
     concentration = None
     transport_budget = None
     transport_budget = None
     if has_transport:
         concentration = simulation.open_concentration()
         transport_budget = simulation.open_transport_budget(species)
-        transport_budget = transport_budget.sel(time=364)
+        transport_budget = transport_budget.isel(time=-1)
     return head, concentration, flow_budget, transport_budget