Fractional open-water coverage for the recharge mask and polder DRN

src/nhflotools/geoconverter/utils.py

@@ -119,7 +119,7 @@ def read_tabular(
     y_column: str | None,
     wkt_column: str | None,
     target_crs: str = "EPSG:28992",
-    **kwargs
+    **kwargs,
 ) -> gpd.GeoDataFrame:
     """Read tabular data with geometry information.
 

src/nhflotools/lakes.py

@@ -314,10 +314,12 @@ def recharge_pond_mask(ds, panden_riv=None, *, fractional=False):
     fractional : bool, optional
         When False (default) return the boolean mask: True at carved lake cells and panden
         RIV cells. When True return a float fraction per cell instead: 1.0 at carved lake
-        cells (the whole cell is modeled as lake bed held by the RIV) and at panden RIV
-        cells, and the open-water coverage fraction at cells below the carve threshold —
-        precipitation on such a partial lake sliver feeds the lake (whose balance the
-        prescribed stage absorbs), not the aquifer, so the caller can scale recharge by
+        cells (the whole cell is modeled as lake bed held by the RIV), the panden-covered
+        fraction from ``ds['panden_coverage']`` (set by
+        :func:`nhflotools.panden.riv_from_oppervlakte_pwn`; 1.0 at panden RIV cells when
+        that variable is absent), and the open-water coverage fraction at cells below the
+        carve threshold — precipitation on such a partial lake or pand sliver feeds the
+        stage boundary, not the aquifer, so the caller can scale recharge by
         ``1 - fraction`` there instead of keeping the full land P-E.
 
     Returns
@@ -332,4 +334,10 @@ def recharge_pond_mask(ds, panden_riv=None, *, fractional=False):
         mask.loc[{"icell2d": cells}] = True
     if not fractional:
         return mask
-    return xr.where(mask, 1.0, ds["lake_coverage"], keep_attrs=True)
+    if "panden_coverage" in ds:
+        coverage = ds["lake_coverage"] + ds["panden_coverage"]
+    elif panden_riv is not None:
+        coverage = xr.where(mask & ~ds["lake_cell"], 1.0, ds["lake_coverage"])
+    else:
+        coverage = ds["lake_coverage"]
+    return xr.where(ds["lake_cell"], 1.0, coverage.clip(max=1.0), keep_attrs=True)

src/nhflotools/panden.py

@@ -5,9 +5,11 @@
 import geopandas as gpd
 import nlmod
 import numpy as np
+import xarray as xr
 
 logger = logging.getLogger(__name__)
 
+
 def get_oppervlakte_pwn_shapes(data_path_panden):
     """Get oppervlakte shapes for PWN area.
 
@@ -59,14 +61,23 @@ def riv_from_oppervlakte_pwn(ds, gwf, data_path_panden):
     Returns
     -------
     flopy.mf6.ModflowGwfriv or None
-        RIV package, or None when no data intersects the model extent.
+        RIV package, or None when no data intersects the model extent. As a side effect
+        ``ds['panden_coverage']`` is set: the panden-covered fraction of each cell, used by
+        :func:`nhflotools.lakes.recharge_pond_mask` for the fractional recharge mask.
     """
     panden_shp = get_oppervlakte_pwn_shapes(data_path_panden=data_path_panden)
     rivdata = nlmod.grid.gdf_to_grid(panden_shp, ds)
     if rivdata.empty:
         logger.warning("No RIV data found within the provided extent.")
         return None
 
+    piece_area = rivdata.groupby("cellid")["area"].sum()
+    coverage = xr.zeros_like(ds["area"], dtype=float)
+    coverage.loc[{"icell2d": piece_area.index.to_numpy()}] = (
+        piece_area.to_numpy() / ds["area"].sel(icell2d=piece_area.index.to_numpy()).to_numpy()
+    ).clip(max=1.0)
+    ds["panden_coverage"] = coverage
+
     rivdata["cond"] = rivdata["area"] / rivdata["c"]
     agg = nlmod.grid.aggregate_vector_per_cell(
         rivdata,

src/nhflotools/polder.py

@@ -34,12 +34,13 @@ def drn_from_waterboard_data(ds, gwf, wb="Hollands Noorderkwartier", cbot=1.0, e
         Bottom resistance of the drains [days], by default 1.0. The per-cell conductance
         is ``cell_area / cbot``.
     exclude : xarray.DataArray or None, optional
-        Boolean mask over ``icell2d`` marking cells at which no DRN reach should be
-        emitted. Both the conductance and the elevation are nulled at these cells, so a
-        reach is dropped whether its stage comes from HHNK peilgebied data or from the
-        maaiveld fallback (``nlmod.gwf.drn`` masks on ``cond > 0``). The default None
-        applies no exclusion and is fully backward-compatible. Used to hand carved lake
-        cells over to a dedicated stage boundary.
+        Boolean mask or open-water fraction (0-1) over ``icell2d``. The conductance is
+        scaled by the drainable land share ``1 - exclude`` (a boolean ``True`` scales by
+        0), so a reach is dropped at fully open-water cells — whether its stage comes from
+        HHNK peilgebied data or from the maaiveld fallback (``nlmod.gwf.drn`` masks on
+        ``cond > 0``) — and only drains the land share elsewhere. The default None applies
+        no exclusion and is fully backward-compatible. Used to hand carved lake cells over
+        to a dedicated stage boundary and to keep open water out of the drained area.
 
     Returns
     -------
@@ -91,10 +92,11 @@ def drn_from_waterboard_data(ds, gwf, wb="Hollands Noorderkwartier", cbot=1.0, e
     ds["drn_cond"] = drn_cond
 
     if exclude is not None:
-        # Null both fields at excluded cells. nlmod.gwf.drn builds a reach wherever
-        # cond > 0, so nulling drn_cond drops the reach whether its stage came from the
-        # celldata assignment or the maaiveld fallback above.
-        ds["drn_cond"] = ds["drn_cond"].where(~exclude)
-        ds["drn_elev"] = ds["drn_elev"].where(~exclude)
+        # Scale the conductance by the drainable land share; nlmod.gwf.drn builds a reach
+        # wherever cond > 0, so a zero land share drops the reach whether its stage came
+        # from the celldata assignment or the maaiveld fallback above.
+        land = (1.0 - exclude.astype(float)).clip(min=0.0)
+        ds["drn_cond"] *= land
+        ds["drn_elev"] = ds["drn_elev"].where(land > 0)
 
     return nlmod.gwf.drn(ds, gwf, elev="drn_elev", cond="drn_cond")

tests/test_lakes.py

@@ -384,3 +384,23 @@ def test_lak_gdf_from_lakes_pwn_returns_none_when_no_lake_cells(disv_grid):
     sliver = _lake_gdf([_piece(CELL_40, box(0, 400, 20, 420), 2.0, 1.0)])
     ds_carved, _ = lakes.carve_lake_cells(ds, sliver, min_area_fraction=0.5)
     assert lakes.lak_gdf_from_lakes_pwn(ds_carved, sliver, min_area_fraction=0.5) is None
+
+
+def test_pond_mask_fractional_uses_panden_coverage(disv_grid):
+    """With ``ds['panden_coverage']`` present the panden share is fractional and sums with the lake share."""
+    ds, _gwf, geoms = disv_grid()
+    rows = [
+        _piece(CELL_40, box(0, 400, 80, 600), 2.0, 1.0),  # 40% lake -> not carved
+        _piece(CELL_FULL, geoms[CELL_FULL], 2.0, 1.0),  # carved
+    ]
+    ds_carved, _ = lakes.carve_lake_cells(ds, _lake_gdf(rows), min_area_fraction=0.5)
+    panden_coverage = xr.zeros_like(ds_carved["lake_coverage"])
+    panden_coverage.loc[{"icell2d": PANDEN_A}] = 0.3
+    panden_coverage.loc[{"icell2d": CELL_40}] = 0.9
+    ds_carved["panden_coverage"] = panden_coverage
+
+    frac = lakes.recharge_pond_mask(ds_carved, None, fractional=True)
+
+    assert frac.sel(icell2d=PANDEN_A).item() == pytest.approx(0.3)  # fraction, not 1.0
+    assert frac.sel(icell2d=CELL_40).item() == pytest.approx(1.0)  # 0.4 lake + 0.9 panden, clipped
+    assert frac.sel(icell2d=CELL_FULL).item() == pytest.approx(1.0)  # carved stays fully covered

tests/test_panden.py

@@ -182,3 +182,14 @@ def test_riv_ssm_registration_is_transport_gated_and_idempotent(tmp_path, transp
         assert ds.attrs["ssm_sources"] == [riv.package_name]
     else:
         assert "ssm_sources" not in ds.attrs
+
+
+def test_riv_stores_panden_coverage(tmp_path, panden_ds_gwf):
+    """The panden-covered fraction of each cell is stored for the fractional recharge mask."""
+    ds, gwf = panden_ds_gwf
+    # 50 m x 50 m pand wholly inside cell 0 of the 2x2 grid of 100 m x 100 m cells.
+    path = _write_panden_shp(tmp_path, [box(10, 110, 60, 160)], ["ICAS-noord"])
+
+    riv_from_oppervlakte_pwn(ds, gwf, data_path_panden=path)
+
+    np.testing.assert_allclose(ds["panden_coverage"].values, [2500.0 / 10000.0, 0.0, 0.0, 0.0])

tests/test_polder.py

@@ -17,6 +17,7 @@
 import geopandas as gpd
 import nlmod
 import numpy as np
+import pytest
 import xarray as xr
 from shapely.geometry import box
 
@@ -187,3 +188,20 @@ def test_drn_fallback_excludes_lake_cells(disv_grid, monkeypatch):
     assert CELL_REAL_STAGE not in excl_cells
     assert CELL_FALLBACK in excl_cells  # a fallback (maaiveld) cell is unaffected
     assert excl_cells == base_cells - {CELL_REAL_STAGE}
+
+
+def test_drn_fractional_exclude_scales_land_share(disv_grid, monkeypatch):
+    """A fractional exclude scales the conductance by ``1 - fraction`` and drops fully covered cells."""
+    ds, gwf, _geoms = disv_grid()
+    frac = xr.zeros_like(ds["top"]).astype(float)
+    frac.loc[{"icell2d": CELL_REAL_STAGE}] = 0.4  # 40% open water at the real-stage cell
+    frac.loc[{"icell2d": CELL_FALLBACK}] = 1.0  # fully open water -> reach dropped
+
+    drn = _run_drn(monkeypatch, ds, gwf, exclude=frac)
+    conds = {rec["cellid"][-1]: float(rec["cond"]) for rec in drn.stress_period_data.data[0]}
+
+    # cbot=1.0, so the unscaled conductance equals the cell area.
+    assert conds[CELL_REAL_STAGE] == pytest.approx(0.6 * float(ds["area"].sel(icell2d=CELL_REAL_STAGE)))
+    assert CELL_FALLBACK not in conds
+    untouched = next(c for c in conds if c not in {CELL_REAL_STAGE, CELL_FALLBACK})
+    assert conds[untouched] == pytest.approx(float(ds["area"].sel(icell2d=untouched)))