From 8eafdc56a4102d46c53127ad7bf8da8afd834e02 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 14 Jan 2025 18:02:23 +0100
Subject: [PATCH 01/12] Average timeseries upon receving only one start time,
 meaning only one stress period.

---
 imod/mf6/wel.py                               | 11 ++++++++--
 .../test_utilities/test_resampling.py         | 21 ++++++++++++++++++-
 imod/util/expand_repetitions.py               | 21 ++++++++++++++++++-
 3 files changed, 49 insertions(+), 4 deletions(-)

diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index abefe7c23..f9cadf096 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -47,7 +47,7 @@
 from imod.select.points import points_indices, points_values
 from imod.typing import GridDataArray, Imod5DataDict
 from imod.typing.grid import is_spatial_grid, ones_like
-from imod.util.expand_repetitions import resample_timeseries
+from imod.util.expand_repetitions import average_timeseries, resample_timeseries
 from imod.util.structured import values_within_range
 
 ABSTRACT_METH_ERROR_MSG = "Method in abstract base class called"
@@ -119,12 +119,19 @@ def _prepare_well_rates_from_groups(
     if has_associated:
         # Resample times per group
         unique_well_groups = [
-            resample_timeseries(df_group, start_times)
+            _process_timeseries(df_group, start_times)
             for df_group in unique_well_groups
         ]
     return _df_groups_to_da_rates(unique_well_groups)
 
 
+def _process_timeseries(df_group, start_times):
+    if len(start_times) > 1:
+        return resample_timeseries(df_group, start_times)
+    else:
+        return average_timeseries(df_group)
+
+
 def _prepare_df_ipf_associated(
     pkg_data: dict, start_times: list[datetime], all_well_times: list[datetime]
 ) -> pd.DataFrame:
diff --git a/imod/tests/test_mf6/test_utilities/test_resampling.py b/imod/tests/test_mf6/test_utilities/test_resampling.py
index e24a9b722..c5fe7cd5b 100644
--- a/imod/tests/test_mf6/test_utilities/test_resampling.py
+++ b/imod/tests/test_mf6/test_utilities/test_resampling.py
@@ -2,7 +2,7 @@
 
 import pandas as pd
 
-from imod.util.expand_repetitions import resample_timeseries
+from imod.util.expand_repetitions import average_timeseries, resample_timeseries
 
 
 def initialize_timeseries(times: list[datetime], rates: list[float]) -> pd.DataFrame:
@@ -160,3 +160,22 @@ def test_timeseries_resampling_refine_and_coarsen():
     pd.testing.assert_frame_equal(
         original_timeseries, re_coarsened_timeseries, check_dtype=False
     )
+
+
+def test_mean_timeseries():
+    # In this test, we compute the mean of a timeseries.
+    times = [datetime(1989, 1, i) for i in [1, 3, 4, 5, 6]]
+    rates = [i * 100 for i in range(1, 6)]
+    timeseries = initialize_timeseries(times, rates)
+
+    mean_timeseries = average_timeseries(timeseries)
+
+    dummy_times = [datetime(1989, 1, 1)]
+    expected_rates = [300.0]
+    expected_timeseries = initialize_timeseries(dummy_times, expected_rates)
+    col_order = ["x", "y", "id", "filt_top", "filt_bot", "rate"]
+    expected_timeseries = expected_timeseries[col_order]
+
+    pd.testing.assert_frame_equal(
+        mean_timeseries, expected_timeseries, check_dtype=False
+    )
diff --git a/imod/util/expand_repetitions.py b/imod/util/expand_repetitions.py
index 60c358b14..361b1ad31 100644
--- a/imod/util/expand_repetitions.py
+++ b/imod/util/expand_repetitions.py
@@ -44,6 +44,24 @@ def expand_repetitions(
     return expanded
 
 
+def average_timeseries(well_rate: pd.DataFrame) -> pd.DataFrame:
+    """
+    Compute the mean value of the timeseries for a single well. Time column is
+    dropped.
+
+    Parameters
+    ----------
+    well_rate:  pd.DataFrame
+        input timeseries for a single well
+    """
+    # Take first item
+    output_frame = well_rate.iloc[[0]].drop(columns=["time", "rate"])
+    # Compute mean
+    col_index = output_frame.shape[1]
+    output_frame.insert(col_index, "rate", well_rate["rate"].mean())
+    return output_frame
+
+
 def resample_timeseries(
     well_rate: pd.DataFrame, times: list[datetime] | pd.DatetimeIndex
 ) -> pd.DataFrame:
@@ -98,7 +116,8 @@ def resample_timeseries(
 
     # compute time difference from perious to current row
     time_diff_col = intermediate_df["time"].diff()
-    intermediate_df.insert(7, "time_to_next", time_diff_col.values)
+    col_index = intermediate_df.shape[1] - 1
+    intermediate_df.insert(col_index, "time_to_next", time_diff_col.values)
 
     # shift the new column 1 place down so that they become the time to the next row
     intermediate_df["time_to_next"] = intermediate_df["time_to_next"].shift(-1)

From 665c240267da31ddb865c5983e3bd188ac9bb8e4 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 14 Jan 2025 18:11:00 +0100
Subject: [PATCH 02/12] Update docstring

---
 imod/mf6/wel.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index f9cadf096..f0dd798a3 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -578,10 +578,13 @@ def from_imod5_data(
             * Multiplication and addition factors need to remain constant through time
             * Same associated well cannot be assigned to multiple layers
         - The dataframe of the first projectfile timestamp is selected
-        - Rate timeseries are resampled with a time weighted mean to the
-          simulation times.
-        - When simulation times fall outside well timeseries range, the last
-          rate is forward filled.
+        - Timeseries are processed as follows:
+            * If ``len(times) > 2``, rate timeseries are resampled with a time
+              weighted mean to the simulation times. When simulation times fall
+              outside well timeseries range, the last rate is forward filled.
+            * If ``len(times) == 2``, the simulation is assumed to be
+              "steady-state" and an average rate is computed from the
+              timeseries.
         - Projectfile timestamps are not used. Even if assigned to a
           "steady-state" timestamp, the resulting dataset still uses simulation
           times.

From 0a75707e15bd4daf0fe99a36a5a44af504585a23 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 14 Jan 2025 18:14:53 +0100
Subject: [PATCH 03/12] Update changelog

---
 docs/api/changelog.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index be7a0670f..f9971fc8e 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -27,6 +27,10 @@ Changed
   > 0 to indicate an active cell and <0 to indicate a vertical passthrough cell,
   consistent with MODFLOW6. Previously this could only be indicated with 1 and
   -1.
+- :meth:`imod.mf6.Well.from_imod5_data` and
+  :meth:`imod.mf6.LayeredWell.from_imod5_data` upon receiving ``len(times) <=
+  2``, the simulation is assumed to be "steady-state" and well timeseries are
+  averaged.
 
 Fixed
 ~~~~~

From a9dc44e8fa56a561b6e78450e366887f3f2f6e77 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 15 Jan 2025 10:34:17 +0100
Subject: [PATCH 04/12] Update unused argument in the signature

---
 imod/mf6/wel.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index f0dd798a3..6ec0d0fae 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -133,7 +133,7 @@ def _process_timeseries(df_group, start_times):
 
 
 def _prepare_df_ipf_associated(
-    pkg_data: dict, start_times: list[datetime], all_well_times: list[datetime]
+    pkg_data: dict, all_well_times: list[datetime]
 ) -> pd.DataFrame:
     """Prepare dataframe for an ipf with associated timeseries in a textfile."""
     # Validate if associated wells are assigned multiple layers, factors,
@@ -210,7 +210,7 @@ def _unpack_package_data(
     start_times = times[:-1]  # Starts stress periods.
     has_associated = pkg_data["has_associated"]
     if has_associated:
-        return _prepare_df_ipf_associated(pkg_data, start_times, all_well_times)
+        return _prepare_df_ipf_associated(pkg_data, all_well_times)
     else:
         return _prepare_df_ipf_unassociated(pkg_data, start_times)
 

From 574e35ad7934696fe2c587240fcf68daa019a752 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 15 Jan 2025 10:45:16 +0100
Subject: [PATCH 05/12] Minor reduction of code duplication

---
 imod/mf6/wel.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index 6ec0d0fae..5d053fa0b 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -108,14 +108,13 @@ def _df_groups_to_da_rates(
 def _prepare_well_rates_from_groups(
     pkg_data: dict,
     unique_well_groups: pd.api.typing.DataFrameGroupBy,
-    times: list[datetime],
+    start_times: list[datetime],
 ) -> xr.DataArray:
     """
     Prepare well rates from dataframe groups, grouped by unique well locations.
     Resample timeseries if ipf with associated text files.
     """
     has_associated = pkg_data["has_associated"]
-    start_times = times[:-1]  # Starts stress periods.
     if has_associated:
         # Resample times per group
         unique_well_groups = [
@@ -204,10 +203,9 @@ def _prepare_df_ipf_unassociated(
 
 
 def _unpack_package_data(
-    pkg_data: dict, times: list[datetime], all_well_times: list[datetime]
+    pkg_data: dict, start_times: list[datetime], all_well_times: list[datetime]
 ) -> pd.DataFrame:
     """Unpack package data to dataframe"""
-    start_times = times[:-1]  # Starts stress periods.
     has_associated = pkg_data["has_associated"]
     if has_associated:
         return _prepare_df_ipf_associated(pkg_data, all_well_times)
@@ -642,7 +640,8 @@ def from_imod5_data(
         pkg_data = imod5_data[key]
         all_well_times = get_all_imod5_prj_well_times(imod5_data)
 
-        df = _unpack_package_data(pkg_data, times, all_well_times)
+        start_times = times[:-1]  # Starts stress periods.
+        df = _unpack_package_data(pkg_data, start_times, all_well_times)
         cls._validate_imod5_depth_information(key, pkg_data, df)
 
         # Groupby unique wells, to get dataframes per time.
@@ -660,7 +659,7 @@ def from_imod5_data(
             for (var, dtype), value in zip(varnames, index_values)
         }
         cls_input["rate"] = _prepare_well_rates_from_groups(
-            pkg_data, unique_well_groups, times
+            pkg_data, unique_well_groups, start_times
         )
         cls_input["minimum_k"] = minimum_k
         cls_input["minimum_thickness"] = minimum_thickness

From febd4cd6e94ce74b23bfca8193b663afba18ff00 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 15 Jan 2025 14:26:41 +0100
Subject: [PATCH 06/12] For GridAgnosticWell.from_imod5_data the ``times`` arg
 now accepts "steady-state" for a steady-state simulation.

---
 docs/api/changelog.rst              |  6 ++--
 imod/mf6/wel.py                     | 51 +++++++++++++++++++++--------
 imod/tests/test_mf6/test_mf6_wel.py | 13 ++++++++
 3 files changed, 53 insertions(+), 17 deletions(-)

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index f9971fc8e..a9a78e20e 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -28,9 +28,9 @@ Changed
   consistent with MODFLOW6. Previously this could only be indicated with 1 and
   -1.
 - :meth:`imod.mf6.Well.from_imod5_data` and
-  :meth:`imod.mf6.LayeredWell.from_imod5_data` upon receiving ``len(times) <=
-  2``, the simulation is assumed to be "steady-state" and well timeseries are
-  averaged.
+  :meth:`imod.mf6.LayeredWell.from_imod5_data` now also accept the argument
+  ``times = "steady-state"``, for the simulation is assumed to be "steady-state"
+  and well timeseries are averaged.
 
 Fixed
 ~~~~~
diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index 5d053fa0b..bb729ced1 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -5,7 +5,7 @@
 import textwrap
 import warnings
 from datetime import datetime
-from typing import Any, Optional, Tuple, Union
+from typing import Any, Literal, Optional, Tuple, Union
 
 import cftime
 import numpy as np
@@ -125,10 +125,10 @@ def _prepare_well_rates_from_groups(
 
 
 def _process_timeseries(df_group, start_times):
-    if len(start_times) > 1:
-        return resample_timeseries(df_group, start_times)
-    else:
+    if _times_is_steady_state(start_times):
         return average_timeseries(df_group)
+    else:
+        return resample_timeseries(df_group, start_times)
 
 
 def _prepare_df_ipf_associated(
@@ -161,7 +161,7 @@ def _prepare_df_ipf_associated(
 
 
 def _prepare_df_ipf_unassociated(
-    pkg_data: dict, start_times: list[datetime]
+    pkg_data: dict, start_times: list[datetime] | Literal["steady-state"]
 ) -> pd.DataFrame:
     """Prepare dataframe for an ipf with no associated timeseries."""
     is_steady_state = any(t is None for t in pkg_data["time"])
@@ -191,6 +191,10 @@ def _prepare_df_ipf_unassociated(
     da_multi = df_multi.to_xarray()
     indexers = {"ipf_row": ipf_row_index}
     if not is_steady_state:
+        if start_times == "steady-state":
+            raise ValueError(
+                "start_times cannot be 'steady-state' for transient wells without associated timeseries."
+            )
         indexers["time"] = start_times
     # Multi-dimensional reindex, forward fill well locations, fill well rates
     # with 0.0.
@@ -292,6 +296,21 @@ def derive_cellid_from_points(
     return cellid.astype(int)
 
 
+def _times_is_steady_state(times: list[datetime] | Literal["steady-state"]) -> bool:
+    return isinstance(times, str) and times == "steady-state"
+
+
+def _get_starttimes(times: list[datetime] | Literal["steady-state"]) -> list[datetime] | Literal["steady-state"]:
+    if _times_is_steady_state(times):
+        return times
+    elif hasattr(times, '__iter__') and isinstance(times[0], (datetime, np.datetime64, pd.Timestamp)):
+        return times[:-1]
+    else:
+        raise ValueError(
+            "Only 'steady-state' or a list of datetimes are supported for ``times``."
+        )
+
+
 class GridAgnosticWell(BoundaryCondition, IPointDataPackage, abc.ABC):
     """
     Abstract base class for grid agnostic wells
@@ -544,7 +563,7 @@ def from_imod5_data(
         cls,
         key: str,
         imod5_data: dict[str, dict[str, GridDataArray]],
-        times: list[datetime],
+        times: list[datetime] | Literal["steady-state"],
         minimum_k: float = 0.1,
         minimum_thickness: float = 0.05,
     ) -> "GridAgnosticWell":
@@ -576,11 +595,13 @@ def from_imod5_data(
             * Multiplication and addition factors need to remain constant through time
             * Same associated well cannot be assigned to multiple layers
         - The dataframe of the first projectfile timestamp is selected
-        - Timeseries are processed as follows:
-            * If ``len(times) > 2``, rate timeseries are resampled with a time
-              weighted mean to the simulation times. When simulation times fall
-              outside well timeseries range, the last rate is forward filled.
-            * If ``len(times) == 2``, the simulation is assumed to be
+        - Timeseries are processed based on the ``times`` argument of this
+          method:
+            * If ``times`` is a list of datetimes, rate timeseries are resampled
+              with a time weighted mean to the simulation times. When simulation
+              times fall outside well timeseries range, the last rate is forward
+              filled.
+            * If ``times = "steady-state"``, the simulation is assumed to be
               "steady-state" and an average rate is computed from the
               timeseries.
         - Projectfile timestamps are not used. Even if assigned to a
@@ -626,8 +647,9 @@ def from_imod5_data(
         imod5_data: dict
             iMOD5 data loaded from a projectfile with
             :func:`imod.formats.prj.open_projectfile_data`
-        times: list
-            Simulation times
+        times: list[datetime] | Literal["steady-state"]
+            Simulation times, a list of datetimes for transient simulations. Or
+            the string ``"steady-state"`` for steady-state simulations.
         minimum_k: float, optional
             On creating point wells, no point wells will be placed in cells with
             a lower horizontal conductivity than this. Wells are placed when
@@ -637,10 +659,11 @@ def from_imod5_data(
             a lower thickness than this. Wells are placed when ``to_mf6_pkg`` is
             called.
         """
+
         pkg_data = imod5_data[key]
         all_well_times = get_all_imod5_prj_well_times(imod5_data)
 
-        start_times = times[:-1]  # Starts stress periods.
+        start_times = _get_starttimes(times) # Starts stress periods.
         df = _unpack_package_data(pkg_data, start_times, all_well_times)
         cls._validate_imod5_depth_information(key, pkg_data, df)
 
diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py
index 91475cef7..44436cabe 100644
--- a/imod/tests/test_mf6/test_mf6_wel.py
+++ b/imod/tests/test_mf6/test_mf6_wel.py
@@ -1019,6 +1019,19 @@ def test_import_and_convert_to_mf6(imod5_dataset, tmp_path, wel_class):
     mf6_well._write("wel", [], write_context)
 
 
+@parametrize("wel_class", [Well, LayeredWell])
+@pytest.mark.usefixtures("imod5_dataset")
+def test_import__as_steady_state(imod5_dataset, wel_class):
+    data = imod5_dataset[0]
+    times = "steady-state"
+    # Import grid-agnostic well from imod5 data (it contains 1 well)
+    wel = wel_class.from_imod5_data("wel-WELLS_L3", data, times)
+
+    assert "time" not in wel.dataset.coords
+    assert wel.dataset["rate"].shape == (1,)
+    np.testing.assert_almost_equal(wel.dataset["rate"].values, -323.89361702)
+
+
 @parametrize("wel_class", [Well])
 @pytest.mark.usefixtures("imod5_dataset")
 def test_import_and_cleanup(imod5_dataset, wel_class: Well):

From f094c4fef4491df376d6e79c36dc200dcd3a1af2 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 15 Jan 2025 14:32:37 +0100
Subject: [PATCH 07/12] If no storage package, assume steady-state and provide
 wel times as "steady-state"

---
 imod/mf6/model_gwf.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index ebc878876..b3cac1104 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -271,7 +271,8 @@ def from_imod5_data(
         result["npf"] = npf_pkg
 
         # import sto
-        if "sto" in imod5_data.keys():
+        steady_state = "sto" in imod5_data.keys()
+        if steady_state:
             result["sto"] = StorageCoefficient.from_imod5_data(
                 imod5_data,
                 grid,
@@ -308,6 +309,10 @@ def from_imod5_data(
         imod5_keys = list(imod5_data.keys())
 
         # import wells
+        if steady_state:
+            wel_times = "steady-state"
+        else:
+            wel_times = times
         wel_keys = [key for key in imod5_keys if key[0:3] == "wel"]
         for wel_key in wel_keys:
             wel_key_truncated = wel_key[:16]
@@ -325,11 +330,11 @@ def from_imod5_data(
             layer = np.array(imod5_data[wel_key]["layer"])
             if np.any(layer == 0):
                 result[wel_key_truncated] = Well.from_imod5_data(
-                    wel_key, imod5_data, times
+                    wel_key, imod5_data, wel_times
                 )
             else:
                 result[wel_key_truncated] = LayeredWell.from_imod5_data(
-                    wel_key, imod5_data, times
+                    wel_key, imod5_data, wel_times
                 )
 
         if "cap" in imod5_keys:

From 89c4badb89dfe6e52e2d23d73b5a86d5cb9c5e58 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 15 Jan 2025 14:32:44 +0100
Subject: [PATCH 08/12] Format

---
 imod/mf6/wel.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index bb729ced1..debbece8f 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -300,10 +300,14 @@ def _times_is_steady_state(times: list[datetime] | Literal["steady-state"]) -> b
     return isinstance(times, str) and times == "steady-state"
 
 
-def _get_starttimes(times: list[datetime] | Literal["steady-state"]) -> list[datetime] | Literal["steady-state"]:
+def _get_starttimes(
+    times: list[datetime] | Literal["steady-state"],
+) -> list[datetime] | Literal["steady-state"]:
     if _times_is_steady_state(times):
         return times
-    elif hasattr(times, '__iter__') and isinstance(times[0], (datetime, np.datetime64, pd.Timestamp)):
+    elif hasattr(times, "__iter__") and isinstance(
+        times[0], (datetime, np.datetime64, pd.Timestamp)
+    ):
         return times[:-1]
     else:
         raise ValueError(
@@ -663,7 +667,7 @@ def from_imod5_data(
         pkg_data = imod5_data[key]
         all_well_times = get_all_imod5_prj_well_times(imod5_data)
 
-        start_times = _get_starttimes(times) # Starts stress periods.
+        start_times = _get_starttimes(times)  # Starts stress periods.
         df = _unpack_package_data(pkg_data, start_times, all_well_times)
         cls._validate_imod5_depth_information(key, pkg_data, df)
 

From 8cbfe9fa9e67493f35ec262e4c0664e406423698 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 15 Jan 2025 14:53:09 +0100
Subject: [PATCH 09/12] Fix mypy errors

---
 imod/mf6/model_gwf.py   |  4 ++--
 imod/mf6/wel.py         | 24 +++++++++++++-----------
 imod/typing/__init__.py |  2 ++
 3 files changed, 17 insertions(+), 13 deletions(-)

diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index b3cac1104..b1d492ce6 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -40,7 +40,7 @@
     SimulationAllocationOptions,
     SimulationDistributingOptions,
 )
-from imod.typing import GridDataArray
+from imod.typing import GridDataArray, StressPeriodTimesType
 from imod.typing.grid import zeros_like
 
 
@@ -310,7 +310,7 @@ def from_imod5_data(
 
         # import wells
         if steady_state:
-            wel_times = "steady-state"
+            wel_times: StressPeriodTimesType = "steady-state"
         else:
             wel_times = times
         wel_keys = [key for key in imod5_keys if key[0:3] == "wel"]
diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index debbece8f..8a451191a 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -5,7 +5,7 @@
 import textwrap
 import warnings
 from datetime import datetime
-from typing import Any, Literal, Optional, Tuple, Union
+from typing import Any, Optional, Tuple, Union, cast
 
 import cftime
 import numpy as np
@@ -45,7 +45,7 @@
     ValidationError,
 )
 from imod.select.points import points_indices, points_values
-from imod.typing import GridDataArray, Imod5DataDict
+from imod.typing import GridDataArray, Imod5DataDict, StressPeriodTimesType
 from imod.typing.grid import is_spatial_grid, ones_like
 from imod.util.expand_repetitions import average_timeseries, resample_timeseries
 from imod.util.structured import values_within_range
@@ -108,7 +108,7 @@ def _df_groups_to_da_rates(
 def _prepare_well_rates_from_groups(
     pkg_data: dict,
     unique_well_groups: pd.api.typing.DataFrameGroupBy,
-    start_times: list[datetime],
+    start_times: StressPeriodTimesType,
 ) -> xr.DataArray:
     """
     Prepare well rates from dataframe groups, grouped by unique well locations.
@@ -124,7 +124,9 @@ def _prepare_well_rates_from_groups(
     return _df_groups_to_da_rates(unique_well_groups)
 
 
-def _process_timeseries(df_group, start_times):
+def _process_timeseries(
+    df_group: pd.api.typing.DataFrameGroupBy, start_times: StressPeriodTimesType
+):
     if _times_is_steady_state(start_times):
         return average_timeseries(df_group)
     else:
@@ -161,7 +163,7 @@ def _prepare_df_ipf_associated(
 
 
 def _prepare_df_ipf_unassociated(
-    pkg_data: dict, start_times: list[datetime] | Literal["steady-state"]
+    pkg_data: dict, start_times: StressPeriodTimesType
 ) -> pd.DataFrame:
     """Prepare dataframe for an ipf with no associated timeseries."""
     is_steady_state = any(t is None for t in pkg_data["time"])
@@ -207,7 +209,7 @@ def _prepare_df_ipf_unassociated(
 
 
 def _unpack_package_data(
-    pkg_data: dict, start_times: list[datetime], all_well_times: list[datetime]
+    pkg_data: dict, start_times: StressPeriodTimesType, all_well_times: list[datetime]
 ) -> pd.DataFrame:
     """Unpack package data to dataframe"""
     has_associated = pkg_data["has_associated"]
@@ -296,19 +298,19 @@ def derive_cellid_from_points(
     return cellid.astype(int)
 
 
-def _times_is_steady_state(times: list[datetime] | Literal["steady-state"]) -> bool:
+def _times_is_steady_state(times: StressPeriodTimesType) -> bool:
     return isinstance(times, str) and times == "steady-state"
 
 
 def _get_starttimes(
-    times: list[datetime] | Literal["steady-state"],
-) -> list[datetime] | Literal["steady-state"]:
+    times: StressPeriodTimesType,
+) -> StressPeriodTimesType:
     if _times_is_steady_state(times):
         return times
     elif hasattr(times, "__iter__") and isinstance(
         times[0], (datetime, np.datetime64, pd.Timestamp)
     ):
-        return times[:-1]
+        return cast(list[datetime], times[:-1])
     else:
         raise ValueError(
             "Only 'steady-state' or a list of datetimes are supported for ``times``."
@@ -567,7 +569,7 @@ def from_imod5_data(
         cls,
         key: str,
         imod5_data: dict[str, dict[str, GridDataArray]],
-        times: list[datetime] | Literal["steady-state"],
+        times: StressPeriodTimesType,
         minimum_k: float = 0.1,
         minimum_thickness: float = 0.05,
     ) -> "GridAgnosticWell":
diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py
index da8bb6c19..c44eb9838 100644
--- a/imod/typing/__init__.py
+++ b/imod/typing/__init__.py
@@ -2,6 +2,7 @@
 Module to define type aliases.
 """
 
+from datetime import datetime
 from typing import TYPE_CHECKING, Literal, TypeAlias, TypedDict, TypeVar, Union
 
 import numpy as np
@@ -17,6 +18,7 @@
 UnstructuredData: TypeAlias = Union[xu.UgridDataset, xu.UgridDataArray]
 FloatArray: TypeAlias = NDArray[np.floating]
 IntArray: TypeAlias = NDArray[np.int_]
+StressPeriodTimesType: TypeAlias = list[datetime] | Literal["steady-state"]
 
 
 class SelSettingsType(TypedDict, total=False):

From 7956c659972e200dbaee26bcbbd89f9de6bb9309 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 15 Jan 2025 15:15:09 +0100
Subject: [PATCH 10/12] Add test for steady state wells

---
 imod/mf6/model_gwf.py                      | 10 +++----
 imod/tests/test_mf6/test_mf6_simulation.py | 32 ++++++++++++++++++++++
 2 files changed, 37 insertions(+), 5 deletions(-)

diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index b1d492ce6..6acf88bff 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -271,8 +271,8 @@ def from_imod5_data(
         result["npf"] = npf_pkg
 
         # import sto
-        steady_state = "sto" in imod5_data.keys()
-        if steady_state:
+        transient = "sto" in imod5_data.keys()
+        if transient:
             result["sto"] = StorageCoefficient.from_imod5_data(
                 imod5_data,
                 grid,
@@ -309,10 +309,10 @@ def from_imod5_data(
         imod5_keys = list(imod5_data.keys())
 
         # import wells
-        if steady_state:
-            wel_times: StressPeriodTimesType = "steady-state"
+        if transient:
+            wel_times: StressPeriodTimesType = times
         else:
-            wel_times = times
+            wel_times = "steady-state"
         wel_keys = [key for key in imod5_keys if key[0:3] == "wel"]
         for wel_key in wel_keys:
             wel_key_truncated = wel_key[:16]
diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py
index 73ec1b31f..8af9de08e 100644
--- a/imod/tests/test_mf6/test_mf6_simulation.py
+++ b/imod/tests/test_mf6/test_mf6_simulation.py
@@ -653,6 +653,38 @@ def test_import_from_imod5__correct_well_type(imod5_dataset):
     assert isinstance(simulation["imported_model"]["wel-WELLS_L5"], LayeredWell)
 
 
+@pytest.mark.unittest_jit
+@pytest.mark.usefixtures("imod5_dataset")
+def test_import_from_imod5__well_steady_state(imod5_dataset):
+    # Unpack
+    imod5_data = imod5_dataset[0]
+    period_data = imod5_dataset[1]
+
+    sto = imod5_data.pop("sto")
+
+    # Other arrangement
+    default_simulation_allocation_options = SimulationAllocationOptions
+    default_simulation_distributing_options = SimulationDistributingOptions
+    datelist = pd.date_range(start="1/1/1989", end="1/1/2013", freq="W")
+
+    # Act
+    simulation = Modflow6Simulation.from_imod5_data(
+        imod5_data,
+        period_data,
+        datelist,
+        default_simulation_allocation_options,
+        default_simulation_distributing_options,
+    )
+    # Assert
+    gwf = simulation["imported_model"]
+    assert "time" not in gwf["wel-WELLS_L3"].dataset.coords
+    assert "time" not in gwf["wel-WELLS_L4"].dataset.coords
+    assert "time" not in gwf["wel-WELLS_L5"].dataset.coords
+    # Teardown
+    # Reassign storage package again
+    imod5_data["sto"] = sto
+
+
 @pytest.mark.unittest_jit
 @pytest.mark.usefixtures("imod5_dataset")
 def test_import_from_imod5__nonstandard_regridding(imod5_dataset, tmp_path):

From db42d9898110999b1712613e720f5ba5daa27150 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 16 Jan 2025 10:16:36 +0100
Subject: [PATCH 11/12] Improve after review

---
 imod/mf6/model_gwf.py                         | 27 +++++++++----------
 imod/mf6/wel.py                               | 23 +++++++++++-----
 .../test_utilities/test_resampling.py         |  3 ++-
 3 files changed, 30 insertions(+), 23 deletions(-)

diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index 6acf88bff..853d52308 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -271,8 +271,8 @@ def from_imod5_data(
         result["npf"] = npf_pkg
 
         # import sto
-        transient = "sto" in imod5_data.keys()
-        if transient:
+        is_transient = "sto" in imod5_data.keys()
+        if is_transient:
             result["sto"] = StorageCoefficient.from_imod5_data(
                 imod5_data,
                 grid,
@@ -309,10 +309,7 @@ def from_imod5_data(
         imod5_keys = list(imod5_data.keys())
 
         # import wells
-        if transient:
-            wel_times: StressPeriodTimesType = times
-        else:
-            wel_times = "steady-state"
+        wel_times: StressPeriodTimesType = times if is_transient else "steady-state"
         wel_keys = [key for key in imod5_keys if key[0:3] == "wel"]
         for wel_key in wel_keys:
             wel_key_truncated = wel_key[:16]
@@ -327,15 +324,15 @@ def from_imod5_data(
                     """
                 )
                 raise KeyError(msg)
-            layer = np.array(imod5_data[wel_key]["layer"])
-            if np.any(layer == 0):
-                result[wel_key_truncated] = Well.from_imod5_data(
-                    wel_key, imod5_data, wel_times
-                )
-            else:
-                result[wel_key_truncated] = LayeredWell.from_imod5_data(
-                    wel_key, imod5_data, wel_times
-                )
+
+            wel_layer = imod5_data[wel_key]["layer"].values
+            is_allocated = np.any(wel_layer == 0)
+            wel_args = (wel_key, imod5_data, wel_times)
+            result[wel_key_truncated] = (
+                Well.from_imod5_data(*wel_args)
+                if is_allocated
+                else LayeredWell.from_imod5_data(*wel_args)
+            )
 
         if "cap" in imod5_keys:
             result["msw-sprinkling"] = LayeredWell.from_imod5_cap_data(imod5_data)  # type: ignore
diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index 8a451191a..50b3fcea0 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -4,6 +4,7 @@
 import itertools
 import textwrap
 import warnings
+from collections.abc import Iterable
 from datetime import datetime
 from typing import Any, Optional, Tuple, Union, cast
 
@@ -127,7 +128,7 @@ def _prepare_well_rates_from_groups(
 def _process_timeseries(
     df_group: pd.api.typing.DataFrameGroupBy, start_times: StressPeriodTimesType
 ):
-    if _times_is_steady_state(start_times):
+    if _is_steady_state(start_times):
         return average_timeseries(df_group)
     else:
         return resample_timeseries(df_group, start_times)
@@ -195,7 +196,7 @@ def _prepare_df_ipf_unassociated(
     if not is_steady_state:
         if start_times == "steady-state":
             raise ValueError(
-                "start_times cannot be 'steady-state' for transient wells without associated timeseries."
+                "``start_times`` cannot be 'steady-state' for transient wells without associated timeseries."
             )
         indexers["time"] = start_times
     # Multi-dimensional reindex, forward fill well locations, fill well rates
@@ -298,18 +299,26 @@ def derive_cellid_from_points(
     return cellid.astype(int)
 
 
-def _times_is_steady_state(times: StressPeriodTimesType) -> bool:
+def _is_steady_state(times: StressPeriodTimesType) -> bool:
+    # Shortcut when not string, to avoid ambigious bitwise "and" operation when
+    # its not.
     return isinstance(times, str) and times == "steady-state"
 
 
+def _is_iterable_of_datetimes(times: StressPeriodTimesType) -> bool:
+    return (
+        isinstance(times, Iterable)
+        and (len(times) > 0)
+        and isinstance(times[0], (datetime, np.datetime64, pd.Timestamp))
+    )
+
+
 def _get_starttimes(
     times: StressPeriodTimesType,
 ) -> StressPeriodTimesType:
-    if _times_is_steady_state(times):
+    if _is_steady_state(times):
         return times
-    elif hasattr(times, "__iter__") and isinstance(
-        times[0], (datetime, np.datetime64, pd.Timestamp)
-    ):
+    elif _is_iterable_of_datetimes(times):
         return cast(list[datetime], times[:-1])
     else:
         raise ValueError(
diff --git a/imod/tests/test_mf6/test_utilities/test_resampling.py b/imod/tests/test_mf6/test_utilities/test_resampling.py
index c5fe7cd5b..ad8fc2be9 100644
--- a/imod/tests/test_mf6/test_utilities/test_resampling.py
+++ b/imod/tests/test_mf6/test_utilities/test_resampling.py
@@ -1,5 +1,6 @@
 from datetime import datetime
 
+import numpy as np
 import pandas as pd
 
 from imod.util.expand_repetitions import average_timeseries, resample_timeseries
@@ -171,7 +172,7 @@ def test_mean_timeseries():
     mean_timeseries = average_timeseries(timeseries)
 
     dummy_times = [datetime(1989, 1, 1)]
-    expected_rates = [300.0]
+    expected_rates = np.mean(rates)
     expected_timeseries = initialize_timeseries(dummy_times, expected_rates)
     col_order = ["x", "y", "id", "filt_top", "filt_bot", "rate"]
     expected_timeseries = expected_timeseries[col_order]

From 095b4c10e2f67356599c06d73b7a6f0b262efed4 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 16 Jan 2025 13:02:28 +0100
Subject: [PATCH 12/12] Fix unittests

---
 imod/mf6/model_gwf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index 853d52308..c7e9a2af0 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -325,7 +325,7 @@ def from_imod5_data(
                 )
                 raise KeyError(msg)
 
-            wel_layer = imod5_data[wel_key]["layer"].values
+            wel_layer = np.array(imod5_data[wel_key]["layer"])
             is_allocated = np.any(wel_layer == 0)
             wel_args = (wel_key, imod5_data, wel_times)
             result[wel_key_truncated] = (