Add nhflotools.lakes (Bergen pond/lake stage boundary) + polder DRN exclude

[bdestombe]

What

Adds nhflotools.lakes and a polder-DRN exclude option — the nhflotools half of the Bergen pond/lake stage-boundary fix (model-script side: NHFLO/models#112 and #119, wired up in NHFLO/models#127).

• nhflotools/lakes.py — a single source-of-truth aggregator _aggregate_lake_cells (one record per lake per cell, plus a per-cell coverage Series) shared by:

  • carve_lake_cells(ds, gdf_lake_grid, min_area_fraction=0.5) — lowers the model top to the deepest lake bottom in cells whose combined lake coverage clears the threshold, and stores two markers: boolean ds['lake_cell'] and float ds['lake_coverage'] (coverage fraction, also below the carve threshold);
  • riv_from_lakes_pwn(...) — holds the carved cells at their prescribed stage with one RIV reach per lake per cell (stage=strt, rbot= that lake's bed, cond=Σ piece_area/clake, boundname=identificatie). Lakes sharing a cell are never merged, so per-lake budgets stay attributable and a future MVR mover/weir setup can address individual lakes;
  • lak_gdf_from_lakes_pwn(...) — LAK variant for detailed studies, as data prep only: it aggregates the same aggregator output to the per-(cell, lake) frame that nlmod.gwf.lake_from_gdf consumes (the model script builds the package and its meteorology directly with nlmod), with an effective bed resistance whose bedleak * cell_area reproduces the RIV summed-piece conductance and one exact strt per lake. The total exchange is nevertheless somewhat weaker than the RIV's: MF6 places the connected cell's half-cell vertical resistance (0.5*Δz/k33) in series with the lakebed for VERTICAL connections (physics-review verified against MF6 to machine precision; 12–16% lower on the PWN layers, and exact equivalence via bedleak is impossible where clake < 0.5*Δz/k33, e.g. the clake=1 Vlotter lakes);
  • recharge_pond_mask(ds, panden_riv=None, *, fractional=False) — flags cells whose meteoric input is carried by a stage boundary so RCH can exclude them. fractional=True returns a per-cell fraction instead (1.0 at carved/panden cells, the open-water coverage at cells below the carve threshold) for recharge * (1 - fraction) scaling.

  Because the carved-cell set and the stage-boundary cell set come from the same aggregator, they are equal by construction — no carved cell is left without a stage, and no stage reach lands on an un-carved cell.

• nhflotools/polder.py — drn_from_waterboard_data(..., exclude=None): nulls drn_cond/drn_elev at excluded cells (backward-compatible default None), so carved lake cells are handed to the dedicated stage boundary instead of a polder drain sitting at the pre-carve maaiveld.

Why

In the default Bergen extent nothing maintains the pond/lake stages (strt 1.92–3.30 m NAP): the panden RIV lies outside the extent, the LAK block is commented out, and no RWS/HHNK boundary touches a lake cell — so carved cells behaved as dry confined depressions, recharge was double-counted over them, and the polder fallback drain sat at the pre-carve surface. See NHFLO/models#112 / #119.

Design decisions (reviewed)

• Per-lake RIV (default): rbot = the lake's own bed caps bed infiltration once the head drops below the lakebed (perched-pond behaviour) while draining freely above the stage; cond = piece_area/clake matches how LAK forms lakebed conductance. In a cell carved by a single lake (the usual case) rbot equals the carved top; where lakes share a cell, the cell is carved to the deepest bed and the shallower lake's reach keeps its own higher rbot — physically each lake perches at its own bed. Confirmed by independent physics review (exact flux equivalence with the previous blended reach for uniform clake; strictly more faithful when beds differ).
• LAK variant available for detailed studies; outlets ('-1' string, [lakes_pwn] External outlet encoded as string '-1' crashes nlmod.gwf.lake.lake_from_gdf data#72) remain out of scope. Known limitation (docstring-documented): on a cell shared by two lakes, MF6 applies each lake's RAINFALL over the full cell area while RCH excluded the cell once, over-applying the meteoric term there — no current lakes_pwn cell is shared.
• Recharge/EVT double-counting (review question): verified impossible in the current build — nlmod.read.knmi.get_recharge (default method='linear') nets Makkink evaporation into the single recharge variable, no EVT package is built anywhere in the script, and the boolean RCH mask drops a cell's entire meteoric term exactly once; the prescribed-stage boundary then carries the lake's own water balance.

Verification

• ruff format / ruff check clean.
• 31/31 unit tests on a synthetic disv grid (no solver, no PWN data, no network): coverage threshold (strict > on combined coverage), deepest-bed carve, unified dropna, carve-set == stage-set identity, per-lake aggregation (pieces merge within a lake, never across; regression test fails on the previous per-cell blend, including the per-lake rbot in a shared cell), top-active-layer placement, boolean + fractional recharge masks, LAK connection/conductance/from_ds-meteorology equivalence, DRN cond/elev nulling, and the empty/all-below-threshold edge.
• Test suite mutation-tested: 13 mutations of the numerical kernel (lake-merge, rbot-blend, threshold boundary, coverage sum→max, stage weighting, conductance, mask branches, LAK conductance) — all caught.
• Not run end-to-end against real PWN data / MF6 — the physical claim (heads at carved cells settle near strt) is unverified here. Pending.

Follow-ups (from review, non-blocking)

• Wire fractional=True recharge scaling into the model script (ds['recharge'] *= 1 - fraction; the per-cell transient KNMI recharge path honours per-cell scaling) — currently the boolean mask is used.
• Assert the lake RIV is registered in ds.attrs['ssm_sources'] (SSM path) and cover the transport=0 branch.
• LAK secondary branches: per-lake outlets and gwt→LKT are unexercised by tests.
• Slim lak_gdf_from_lakes_pwn once Aggregate lake pieces per cell in lake_from_gdf, like RIV/DRN celldata aggregation gwmod/nlmod#582 lands (piece→cell aggregation + strt float tolerance inside lake_from_gdf, mirroring the RIV/DRN surface_water.aggregate idiom): the effective-clake and strt-collapse lines move upstream and the helper reduces to threshold-filtering pieces from the shared aggregator.
• Add an end-to-end MF6 smoke test that a carved cell's head holds near strt.
• Consider factoring the shared RIV-build tail with panden.riv_from_oppervlakte_pwn.
• Generalize into a single consolidated adjust-top step: when more carve sources join (future bathymetry / panden opt-ins), merge all bed elevations into one top with a single set_model_top call (it is a lossy one-way ratchet; monotonic lowering is safe, raising is not) — see the architecture note on NHFLO/models#126.

Relationship

This must land before NHFLO/models#127 (the model script imports carve_lake_cells, riv_from_lakes_pwn, recharge_pond_mask, lak_gdf_from_lakes_pwn). Draft until the end-to-end verification above is done.

[bdestombe]

Cross-ref: the top-lifecycle architecture note on NHFLO/models#126 places nhflotools.lakes in the target lifecycle for 09pwnmodel2: it is the seed of the one consolidated adjust-top step between the layer model and the first flopy package. Verified properties relevant to this module:

• nlmod.layers.set_model_top is value-idempotent for equal input and safe for repeated monotonic lowering (zero-thickness botms track the top), but it is a one-way ratchet: remove_inactive_layers drops layers permanently and a later raise refills only the highest surviving layer (nlmod/dims/layers.py:863-871). When more carve sources join later (bathymetry, panden opt-in), merge them into one top and keep a single set_model_top call inside this module's successor rather than carving incrementally — added as a follow-up checkbox in the description.
• riv_from_lakes_pwn's rbot = carved top keeps the geometry compatible with the deferred LAK upgrade: lake_from_gdf connects claktype="VERTICAL" at the first active layer (nlmod/gwf/lake.py:155-196), i.e. the model top at lake cells is the lake bed either way.
• The decision rule from the note for future scope questions: carve only where the model top itself must be the water bed (LAK/stage-held large bodies, sea bathymetry); linear/small features (ditches, drains, panden) stay rbot-placed via build_spd(layer_method="lay_of_rbot"), which resolves the layer from the current geometry at build time (nlmod/gwf/surface_water.py:481-491) and needs no carve.

[bdestombe]

Does this mean that lake cells are merged? I think it might become more difficult to configure the mover package, as they are assigned to individual lakes. Same goes for weirs

[bdestombe]

Confirmed and fixed in 44a9b1a: aggregation is now per (cellid, identificatie), so lakes sharing a cell each keep their own reach, stage, rbot, cond and boundname — MVR/weir and LAK config can address individual lakes. Regression test fails on the old blend."

[bdestombe]

Not sure from where this function is called, but I presume from the function that created the top in nhflotools

[bdestombe]

Called from the model script directly after the layer model is finalized, before any package build — not from a top-creation function. Lifecycle contract now in the docstring; folding it into one consolidated adjust-top step stays the models#126 follow-up.

[bdestombe]

Thoroughly review MF6 and flopy to confirm that all recharge and evaporation is not already deducted from the EVT or RECH packages.

[bdestombe]

Note that fractional overlap should also be possible