Issue #1381 average well rates at steady state (#1382)

Fixes #1381

# Description

Upon receiving simulation ``times`` with string "steady-state"
``Well.from_imod5_data`` now assumes the simulation is steady-state and
well rates should be averaged.
This brings well rates of packages a lot closer to what iMOD5 produces
in the water balance for steady-state simulations.
``GroundwaterFlowModel.from_imod5_data`` sets ``times`` arg for
``Well.from_imod5_data`` to ``"steady-state"`` if no Storage package is
present in the projectfile.

I opted for this for now as we currently lack of a clear API to
construct steady-state simulations. See also #1308.

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` now also accept the argument
+  ``times = "steady-state"``, for the simulation is assumed to be "steady-state"
+  and well timeseries are averaged.
 
 Fixed
 ~~~~~

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
 
 
@@ -271,7 +271,8 @@ def from_imod5_data(
         result["npf"] = npf_pkg
 
         # import sto
-        if "sto" in imod5_data.keys():
+        is_transient = "sto" in imod5_data.keys()
+        if is_transient:
             result["sto"] = StorageCoefficient.from_imod5_data(
                 imod5_data,
                 grid,
@@ -308,6 +309,7 @@ def from_imod5_data(
         imod5_keys = list(imod5_data.keys())
 
         # import wells
+        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]
@@ -322,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, times
-                )
-            else:
-                result[wel_key_truncated] = LayeredWell.from_imod5_data(
-                    wel_key, imod5_data, times
-                )
+
+            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] = (
+                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

imod/mf6/wel.py

@@ -4,8 +4,9 @@
 import itertools
 import textwrap
 import warnings
+from collections.abc import Iterable
 from datetime import datetime
-from typing import Any, Optional, Tuple, Union
+from typing import Any, Optional, Tuple, Union, cast
 
 import cftime
 import numpy as np
@@ -45,9 +46,9 @@
     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 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"
@@ -108,25 +109,33 @@ 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: StressPeriodTimesType,
 ) -> 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 = [
-            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: pd.api.typing.DataFrameGroupBy, start_times: StressPeriodTimesType
+):
+    if _is_steady_state(start_times):
+        return average_timeseries(df_group)
+    else:
+        return resample_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,
@@ -155,7 +164,7 @@ def _prepare_df_ipf_associated(
 
 
 def _prepare_df_ipf_unassociated(
-    pkg_data: dict, start_times: list[datetime]
+    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"])
@@ -185,6 +194,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.
@@ -197,13 +210,12 @@ 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: StressPeriodTimesType, 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, 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)
 
@@ -287,6 +299,33 @@ def derive_cellid_from_points(
     return cellid.astype(int)
 
 
+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 _is_steady_state(times):
+        return times
+    elif _is_iterable_of_datetimes(times):
+        return cast(list[datetime], 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
@@ -539,7 +578,7 @@ def from_imod5_data(
         cls,
         key: str,
         imod5_data: dict[str, dict[str, GridDataArray]],
-        times: list[datetime],
+        times: StressPeriodTimesType,
         minimum_k: float = 0.1,
         minimum_thickness: float = 0.05,
     ) -> "GridAgnosticWell":
@@ -571,10 +610,15 @@ 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 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
           "steady-state" timestamp, the resulting dataset still uses simulation
           times.
@@ -618,8 +662,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
@@ -629,10 +674,12 @@ 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)
 
-        df = _unpack_package_data(pkg_data, times, all_well_times)
+        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)
 
         # Groupby unique wells, to get dataframes per time.
@@ -650,7 +697,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

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):

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):

imod/tests/test_mf6/test_utilities/test_resampling.py

@@ -1,8 +1,9 @@
 from datetime import datetime
 
+import numpy as np
 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 +161,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 = 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]
+
+    pd.testing.assert_frame_equal(
+        mean_timeseries, expected_timeseries, check_dtype=False
+    )

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):