Start with Ribamod docs

docs/_quarto-manual.yml

@@ -24,6 +24,8 @@ book:
         - coupler_metamod_technical.qmd
         - coupler_metamod_example.ipynb
         - coupler_ribamod_config.qmd
+        - coupler_ribamod_technical.qmd
+        - coupler_ribamod_preprocessing.qmd
         - primod_api/index.qmd
         - coupler_architecture.qmd
         - coupler_release.qmd

docs/_quarto-website.yml

@@ -64,6 +64,8 @@ website:
         - section: "Ribasim - MODFLOW 6"
           contents:
             - href: coupler_ribamod_config.qmd
+            - href: coupler_ribamod_technical.qmd
+            - href: coupler_ribamod_preprocessing.qmd
         - href: primod_api/index.qmd
         - href: coupler_architecture.qmd
         - href: coupler_release.qmd

docs/coupler_ribamod_preprocessing.qmd

@@ -0,0 +1,162 @@
+---
+title: Pre-processing
+---
+
+This document describes how to setup coupled Ribasim-MODFLOW 6 simulations.
+
+## The primod Python package
+
+To aid in setting up coupled models, we have created the `primod` Python
+package. It takes both models, derives the exchange relationships between
+both, and writes: 
+
+* the Ribasim model
+* the MODFLOW 6 simulation
+* the exchange files
+* the iMOD coupler configuration file 
+
+The derivation of exchange connections between the models is automatic, and
+based on the spatial information provided for both models.
+
+As `primod` is a Python package, both models must be represented in Python:
+
+* The MODFLOW 6 simulation is represented by the
+  [Modflow6Simulation](https://deltares.github.io/imod-python/api/generated/mf6/imod.mf6.Modflow6Simulation.html)
+  class of the `imod` Python package.
+* The Ribasim model is represented by the
+  [Model](https://deltares.github.io/Ribasim/python/reference/#model) class of
+  the `ribasim` Python package.
+* The combination of both models is represented by the `primod`
+  [RibaMod](primod_api/RibaMod.html#primod.RibaMod)
+  class.
+
+Both Python packages support reading a model from a TOML file and associated
+files.
+
+## Abbreviated example
+
+The following abbreviated example:
+
+* Reads a Ribasim model
+* Reads a MODFLOW 6 simulation
+* Reads a basin definition associated with the Ribasim model
+* Defines a driver coupling: which river and drainage packages are coupled
+* Sets up a coupled model
+* Writes the coupled model
+
+
+``` python
+import ribasim
+import geopandas as gpd
+import imod
+
+from primod.ribamod import DriverCoupling, RibaMod
+
+
+ribasim_model = ribasim.Model.read("ribasim/ribasim.toml")
+mf6_simulation = imod.mf6.Modflow6Simulation.from_file("modflow/GWF_1.toml")
+basin_definition = gpd.read_file("ribasim_network.gpkg", layer="basin_areas")
+
+
+driver_coupling = DriverCoupling(
+    mf6_model="GWF_1",
+    mf6_active_river_packages=["riv-1"],
+    mf6_passive_drainage_packages=["drn-1", "drn-2"],
+)
+
+ribamod_model = RibaMod(
+    ribasim_model=ribasim_model,
+    mf6_simulation=mf6sim,
+    coupling_list=[driver_coupling],
+    basin_definition=basin_definition,
+)
+
+ribamod_model.write(
+   directory="coupled-ribamod-simulation",
+   modflow6_dll=r"c:\bin\imod_coupler\modflow6\libmf6.dll",
+   ribasim_dll=r"c:\bin\imod_coupler\ribasim\bin\libribasim.dll",
+   ribasim_dll_dependency=r"c:\bin\imod_coupler\ribasim\bin",
+)
+```
+
+## Requirements
+
+* The start and end times of the Ribasim and MODFLOW 6 simulations must align.
+  `primod` will raise an error otherwise.
+* Spatial extents of both models need not coincide. Part of the Ribasim basins
+  may be located outside of the MODFLOW 6 simulation window. Uncoupled basins
+  will proceed with the regular drainage and irrigation terms define in the
+  Basin / Static or Basin / Time tables.
+* Similarly, not every River and Drainage boundary needs to be linked with
+  Ribasim. Boundaries outside of any basin polygon will simply use the regular
+  file input.
+
+## Exchange derivation
+
+The exchanges are derived based on spatial overlap of Ribasim basins and the
+grid-based River and Drainage representation. For an active coupling, the x
+and y locations of the subgrid elements is used to link each actively coupled
+MODFLOW 6 boundary with a subgrid element.
+
+
+### Passive coupling
+
+The derivation of exchanges proceeds in the following steps:
+
+* Rasterize the basin definition polygons (provided as a
+  `geopandas.GeoDataFrame`) to the MODFLOW 6 model grid.
+* Overlay the conductance on the rasterized basin definition,
+  and derive for each boundary the basin index.
+* Identify the indices of the coupled MODFLOW 6 boundaries.
+* Store the basin indices and boundary indices in a table.
+
+### Active coupling
+
+The derivation of active coupling exchanges proceeds largely the same, but also
+locates the nearest subgrid elements:
+
+* Rasterize the basin definition polygons (provided as a
+  `geopandas.GeoDataFrame`) to the MODFLOW 6 model grid.
+* Overlay the conductance on the rasterized basin definition, and derive for
+  each boundary the basin index.
+* Identify the indices of the coupled MODFLOW 6 boundaries.
+* Using the `meta_x` and `meta_y` columns in the Ribasim Basin / subgrid table,
+  find the nearest subgrid element for each MODFLOW 6 boundary.
+* Store the basin indices, boundary indices, subgrid indices in a table.
+
+## Modifications to the Ribasim model
+
+The `primod.RibaMod` class makes the following alteration to the Ribasim input
+before writing the Ribasim model: for basins that are coupled to MODFLOW 6,
+the infiltration and drainage columns are set to `NaN` / `Null` (nodata) in
+the Basin / Static or Basin / Time tables.
+
+This ensures that Ribasim does not overwrite the exchange flows while running
+coupled with MODFLOW 6.
+
+Conceptually, it also means that when a basin is coupled, it should generally
+located inside of the MODFLOW 6 model; after all, when half of the basin is
+located outside of the MODFLOW 6 model, it will not receive drainage or lose
+water to infiltration in that half.
+
+## Modifications to the MODFLOW 6 simulation
+
+Currently, no modifications are made in the MODFLOW 6 input.
+
+### Consistency between MODFLOW 6 and Ribasim subgrid
+
+During the coupling, water levels should not be set below the bed elevation
+of the boundary. For drainage packages, this is the drainage elevation
+provided in the MODFLOW 6 input; for river packages, this is the bottom
+elevation provided in the MODFLOW 6 input.
+
+There is potential for inconsistency here, as Ribasim also describes a bed
+elevation: the lowest level of the subgrid piecewise interpolation table:
+
+* In case the MODFLOW 6 bed elevation is higher than the subgrid elevation,
+  infiltration will stop before the Ribasim basin is empty.
+* In case the MODFLOW 6 bed elevation is lower than the subgrid elevation,
+  infiltration will proceed even when the Ribasim basin is empty.
+
+The second is a more pressing problem, as it will results in a discrepancy
+between the water balances of both models.

docs/coupler_ribamod_technical.qmd

@@ -0,0 +1,157 @@
+---
+title: Technical Reference
+---
+
+This document describes how Ribasim and MODFLOW6 are coupled. It is intended
+for groundwater modellers, who need to know which variables are exchanged
+between computational kernels and at which moment. For details of the inner
+workings of the code, we refer to the docstrings in the code.
+
+The following sequence diagram show the current sequential coupling scheme for
+a single MODFLOW 6 stress period.
+
+```{mermaid}
+sequenceDiagram
+    autonumber
+    Ribasim ->> MODFLOW6: exchange stage [T-1]
+    Note over MODFLOW6: solve T
+    MODFLOW6 ->> Ribasim: RIV flux [T-1]
+    Note over Ribasim: solve T
+```
+
+# Requirements
+
+* Ribasim only couples to the River (RIV) and Drainage (DRN) packages of
+  MODFLOW6.
+* One or more River and Drainage boundaries can be connected with one Ribasim
+  basin.
+* The MODFLOW6 River and Drainage entity is the smallest entity: coupling
+  multiple Ribasim basins to a single MODFLOW6 boundary is not supported.
+* The active coupling, in which MODFLOW6 boundary levels are changed, always
+  requires a "subgrid" table in Ribasim.
+
+# Data exchanges
+
+## MODFLOW6 to Ribasim
+
+### Flows
+
+The computed drainage and infiltration flows to and from River and Drainage
+packages are aggregated on basin level and set as:
+
+* infiltration: a flow from the surface water to the groundwater (positive).
+* drainage: a flow from the groundwater to the surface (positive).
+
+These are set as the infiltration and drainage forcing on the Ribasim basins.
+
+## Ribasim to MODFLOW 6
+
+### Levels
+
+Ribasim sets the stage of River packages, and the drainage elevation of Drainage
+packages in MODFLOW6.
+
+The levels of the basins are not set directly in MODFLOW6. The levels are
+interpolated using the "subgrid" functionality of Ribasim, which can be used to
+e.g. create sloping water levels within a single basin.
+
+For simplicity of the coupling, the subgrid functionality is also used in
+case of a 1:1 basin-boundary coupling.
+
+# Files
+
+The following files are required to couple the two model codes.
+
+## MODFLOW6
+
+No specific files are required. The MODFLOW6 model must contain River or
+Drainage packages which represent the surface water.
+
+## Ribasim
+
+No specific files are required. The Basin / Subgrid table must be defined. For
+preprocessing with `primod`, the metadata columns `meta_x`, `meta_y` must be
+included in this table to indicate the spatial location of each subgrid
+element.
+
+## Coupler
+
+These files provide the mappings from the MODFLOW boundary indices to the
+Ribasim basin indices. For actively coupled system (see below), it also
+includes subgrid indices.
+
+The files are identified by the MODFLOW6 package name.
+
+Note that the exchange files are **tab-separated** and use 0-based indices.
+They also include a header.
+
+### Passive coupling
+
+#### [package-name].tsv
+
+
+```
+basin_index     bound_index
+{basin_index0}  {bound_index0}
+{basin_index1}  {bound_index1}
+...
+```
+
+### Active coupling
+
+#### [package-name].tsv
+
+```
+basin_index     bound_index      subgrid_index
+{basin_index0}  {bound_index0}   {subgrid_index0}
+{basin_index1}  {bound_index1}   {subgrid_index1}
+...
+```
+
+# MODFLOW 6 surface water representation
+
+Ribasim is linked to the surface water representation of the MODFLOW 6 models.
+Specifically, Ribasim is linked to either River (RIV) or Drainage (DRN)
+packages. In a coupled simulation, Ribasim computes a surface water mass
+balance and sets the water levels for the MODFLOW 6 boundaries, and MODFLOW 6
+provides the drainage and infiltration flows. These flows are collected in the
+Ribasim basins.
+
+The RibaMod coupling only exchanges between Ribasim basins and MODFLOW 6 River
+and Drainage packages. This means that for a successful coupling with Ribasim,
+the surface water should **only** be represented by either River or Drainage
+boundaries.
+
+More advanced surface water boundaries such as the Streamflow-Routing (SFR)
+package or the Lake (LAK) package compute their own surface water mass balance,
+and cannot be coupled to Ribasim. Such packages may be used freely in the parts
+of the model domain that are **not** coupled with Ribasim, but should not be
+used in the part of the model domain that is coupled to Ribasim.
+
+## Active versus passive coupling
+
+In coupling Ribasim to MODFLOW 6, we make a distinction between active and
+passive coupling.
+
+In passive coupling, MODFLOW 6 computes drainage and infiltration flows. These
+are added to the water balance of the Ribasim basins. The water levels computed
+by Ribasim are **not** set back into MODFLOW 6. This type of coupling is
+convenient for boundaries that show little variation in terms of drainage
+elevation over time (e.g. surface runoff, or ditches with negligible water
+depth). A passive coupling requires little in terms of parametrization of the
+Ribasim model: all flows are simply added as a sort of "lateral" flow into the
+basin.
+
+In the active coupling, Ribasim does set the water levels of the MODFLOW 6
+model.
+
+The passive coupling should generally only be applied to coupling to Drainage
+packages, not River packages: in case of "over-infiltration", the water levels
+in Ribasim will drop, potentially resulting in a dry basin. However, since the
+water levels are not set in MODFLOW 6, no feedback can occur, and the MODFLOW 6
+model will keep infiltrating water that is not available. This creates a
+discrepancy between the water balances of both models.
+
+In the active coupling, large infiltration flows would result in lower water
+levels, and infiltration stops when the basin dries up and the River stage
+reaches the bed elevation.