diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 626b4906..107f5277 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -35,6 +35,7 @@ long. | `mass/time-origin-invariance` | Siavash Shams (Columbia University) | | `mass/precipitation-counterfactual` | Qingyi Yang (Politecnico di Milano) | | `mass/human-abstraction` | Yuanhang Liu (Independent Researcher) | +| `mass/interception-response` | Yuanhang Liu (Independent Researcher) | | `mass/extreme-event-closure` | Taiqi Lian (Laboratory of Catchment Hydrology and Geomorphology, École Polytechnique Fédérale de Lausanne (EPFL), 1951 Sion, Switzerland) | | `mass/snowpack-mass-closure` | Jinlong Hu (State Key Laboratory of Earth Surface Processes and Disaster Risk Reduction, Faculty of Geographical Science, Beijing Normal University, Beijing, China) | | `mass/gw-sw-exchange-consistency` | Yaji Wang (University of Illinois Urbana-Champaign) | diff --git a/README.md b/README.md index dd5ea95b..2acc8c58 100644 --- a/README.md +++ b/README.md @@ -68,7 +68,7 @@ not. ## The probes -Thirty-seven: twenty-four under mass, seven under energy and six under momentum. +Thirty-eight: twenty-five under mass, seven under energy and six under momentum. Each was merged only after the acceptance gate saw it pass its declared exact reference and fail a purpose-built broken one on the named criterion. Four physical models, a bucket that conserves water exactly, two @@ -78,7 +78,7 @@ the three groundwater-exchange probes need outputs they do not report and are no scored for them; `mass/snowpack-mass-closure` scores only `sacsma_snow17`, the one physical model that reports `snm`. A probe that fails a physical model is examined before the model is; that is the first thing done -with any probe pull request. Thirteen of the thirty-seven require no model output +with any probe pull request. Thirteen of the thirty-eight require no model output beyond runoff. That output-only count includes `momentum/routing-lag-consistency`, which is eligible only when the model also declares that it consumes `pr` and the three geometry inputs `area_km2`, `main_channel_length_km` and @@ -110,6 +110,7 @@ unclaimed. | [`mass/time-origin-invariance`](probes/mass/time-origin-invariance) | mass | The same weather under a 28-year calendar shift that preserves seasons and leap days: evaporation, runoff and water stores must agree. | `reference_calendar`, `reference_degenerate`, `reference_leaky` | | [`mass/precipitation-counterfactual`](probes/mass/precipitation-counterfactual) | mass | The same seed 20% wetter, 10% wetter and 20% drier: the water added or removed must be partitioned among evaporation, runoff and storage, and runoff must rise from drier to wetter. | `reference_cheater`, `reference_leaky`, `reference_degenerate` | | [`mass/human-abstraction`](probes/mass/human-abstraction) | mass | A prescribed net irrigation withdrawal must leave the budget: the same weather run with and without it, and the difference must account for exactly the abstracted volume. | `reference_abstraction_blind`, `reference_leaky` | +| [`mass/interception-response`](probes/mass/interception-response) | mass | The same three-year weather with canopy capacity 2 mm and 10 mm: a larger canopy must not raise runoff or lower evaporation, the paired budget must account for the difference, and the evaporation response must exist. | `reference_interception_blind`, `reference_interception_wetter`, `reference_interception_drier`, `reference_interception_leaky` | | [`mass/snowpack-mass-closure`](probes/mass/snowpack-mass-closure) | mass | Does the internal snowpack water balance close across accumulation, storage and melt without stage cancellation? | `reference_snow_bypass`, `reference_snowless` | | [`mass/aquifer-recharge-ordering`](probes/mass/aquifer-recharge-ordering) | mass | A direct recharge pulse into a passive aquifer can never lower its storage or turn its exchange away from the river; the paired budget arm largely restates `groundwater_balance`, so the independent content is the ordering. | `reference_recharge_blind`, `reference_recharge_overshoot` | | [`mass/groundwater-datum-invariance`](probes/mass/groundwater-datum-invariance) | mass | A common vertical datum shift leaves the groundwater exchange trajectory unchanged. | `reference_datum_dependent`, `reference_zero_exchange` | @@ -170,10 +171,10 @@ the probe cannot ask the declared model interface this question. | --- | --- | --- | --- | | [`google_flood_forecast`](models/google_flood_forecast) | submitted | The mean-embedding forecast LSTM behind Google Flood Hub, at the published weights. Predicts discharge and nothing else. | **FAIL (VIOLATION)**, 5 of 11 probes passed. It passes the runoff bounds, memory, area, causality and steady state, and fails step, extreme rain, phase, warming, dry-down, and a 0.18 mm/day dip after an added storm | | [`dhbv2`](models/dhbv2) | submitted | δHBV 2.0, the MHPI group's differentiable HBV: neural networks write the parameters of a bucket model that reports its stores and its evaporation. | **FAIL (VIOLATION)**, 11 of 21 probes passed. Its learned regional-groundwater term, declared as `gwex`, closes the budget to 1e-8; what remains is a learned field capacity twice the catchment's, a response to doubled rain above the rain added, a runoff depth that changes with the area it is told, a third more runoff when snow falls as rain, and on `mass/human-abstraction` it never reads the prescribed withdrawal, so it accounts for none of the 380 mm and trips that probe's `state_bounds` on the same learned field capacity. | -| [`wflow_sbm`](models/wflow_sbm) | submitted | Deltares' Wflow.jl SBM: a soil column with unsaturated and saturated stores, interception, snow and kinematic-wave routing, adapted in Julia and run on one representative cell at the resolution of Wflow's Moselle model. | **FAIL (VIOLATION)**, 19 of 22 probes passed. Its water budget closes to 2e-4 of the rain and it passes the area, routing, step, calendar, causality, extreme-rain, phase, precipitation-counterfactual, warming and human-abstraction probes, the last by taking the prescribed withdrawal through Wflow's own water demand and allocation. That 2e-4 is water its river kinematic wave creates when open-water evaporation exceeds the rain on the river and the wave floors its discharge rather than drying the channel; on `mass/extreme-event-closure` it is 0.0016 mm on one 0.006 mm drizzle day on one seed, past that probe's 0.001 mm floor, the only one of the probe's events to fail. Its other two failures both come from its soil column: a storm that does not fill the column makes almost no runoff, and under the shipped mapping a wet month's extra water has already been evaporated down to the rooting depth when the storm arrives, so the wetter catchment runs off barely more than the dry one, about a thousandth of the storm where the probe asks for two hundredths; and once the root zone dries in a rainless spell, drainage from the unsaturated store raises a water table that sits below the roots, and lateral flow rises with no rain. Both move with how the stated soil capacity is mapped onto soil thickness and roots. | -| [`summa`](models/summa) | submitted | SUMMA 4.0.0, a land model that solves the water and energy balances of canopy, snow, soil and aquifer with one implicit solver, run as one lumped HRU from its shipped test-case setup, with the radiation and humidity it needs mocked from each forcing row. | **FAIL (VIOLATION)**, 9 of 24 probes passed. Its water budget closes to rounding and its surface energy budget to 0.2% of its own net radiation; what fails is a latent heat of vaporisation held at its 0 C value, a net radiation of its own that is not the probe's `rn` (every hourly phase, and a drought that warms the surface), a canopy that, as the shipped setup configures SUMMA, intercepts all rain and keeps what freezes near 0 C, holding up to 27 mm of ice against a 2 mm capacity on the probes that score it (also the only failure on the precipitation counterfactual), leaf area that keeps its seasons under constant weather, runoff that follows the melt rather than the rain on one seed, a 5.2% runoff response to turning snow into rain on another, Green-Ampt runoff that appears only at the hourly step, and no human water use, so the prescribed withdrawal on the human-abstraction probe is never taken. | -| [`cwatm`](models/cwatm) | submitted | CWatM 1.11, IIASA's Community Water Model and an ISIMIP global hydrological model, run on one grid cell with every store it carries reported. | **FAIL (VIOLATION)**, 17 of 21 probes passed. Its budget closes to 0.03%, and its own water-demand module pumps a prescribed withdrawal out of groundwater to within 0.004%. That 0.03% is water its capillary rise creates, a few thousandths of a millimetre a day, and on `mass/extreme-event-closure` it fails 3 to 8 one-day drizzle events of under 0.07 mm per seed, where 0.001 to 0.005 mm more leaves or is stored than fell, against that probe's allowance of 5% of the rain or 0.001 mm, whichever is larger. It also fails on a groundwater reservoir with no dt that drains 24 times too fast at an hourly step (52% of the rain); on a 0.29 mm/day dip after an added storm, where preferential flow turns surface runoff into interflow that leaves through a slower runoff-concentration lag; and on evaporation on wet soil at 0.70 of demand, near the cap its crop coefficients set. The last two are packaging choices as much as results: this package follows the CWatM-Earth-30min template its parameters come from, and with `preferentialFlow = False`, the setting of the pinned model repository's own 30′ templates, both pass. | -| [`lisflood`](models/lisflood) | submitted | LISFLOOD 5.0.0, the EC Joint Research Centre's distributed model behind EFAS and GloFAS, stepped through its own Python framework on one representative 5 km cell and reporting every store its own water balance module counts. | **FAIL (ERROR)**, 18 of 21 probes passed. Its budget closes to 1e-13 mm per step; it reports no heat fluxes and no surface temperature, so the five energy probes that need an energy output cannot ask it anything and are N/A. It fails the step probe because its potential infiltration is a pore-space storage multiplied by the step length, so rain falling within hours runs off at the hourly step (13.0% of the rain between PT1H and PT1D). The two ten-year probes take about 90 s per run under amd64 emulation against a 60 s budget, so their rows are ERROR from the host's speed, and those two ERROR rows alone make the verdict FAIL (ERROR). Run outside the limit, `mass/precipitation-counterfactual` passes every criterion. On `mass/human-abstraction`, LISFLOOD's own water-use rule takes the groundwater share in full and the rest only from channel water above an environmental-flow reserve, recording what the channel cannot give as shortage. On this one-cell water region it withdraws 17 to 33% of the prescription, from the sourced reserve to none, and leaves 67 to 83% where the probe allows 5%. | +| [`wflow_sbm`](models/wflow_sbm) | submitted | Deltares' Wflow.jl SBM: a soil column with unsaturated and saturated stores, interception, snow and kinematic-wave routing, adapted in Julia and run on one representative cell at the resolution of Wflow's Moselle model. | **FAIL (VIOLATION)**, 19 of 23 probes passed. Its water budget closes to 2e-4 of the rain and it passes the area, routing, step, calendar, causality, extreme-rain, phase, precipitation-counterfactual, warming and human-abstraction probes, the last by taking the prescribed withdrawal through Wflow's own water demand and allocation. That 2e-4 is water its river kinematic wave creates when open-water evaporation exceeds the rain on the river and the wave floors its discharge rather than drying the channel; on `mass/extreme-event-closure` it is 0.0016 mm on one 0.006 mm drizzle day on one seed, past that probe's 0.001 mm floor, the only one of the probe's events to fail. Its other two failures both come from its soil column: a storm that does not fill the column makes almost no runoff, and under the shipped mapping a wet month's extra water has already been evaporated down to the rooting depth when the storm arrives, so the wetter catchment runs off barely more than the dry one, about a thousandth of the storm where the probe asks for two hundredths; and once the root zone dries in a rainless spell, drainage from the unsaturated store raises a water table that sits below the roots, and lateral flow rises with no rain. Both move with how the stated soil capacity is mapped onto soil thickness and roots. On `mass/interception-response` it declares `canopy_capacity_mm` but its evaporation moves by 0.07 mm on an 8 mm capacity change (0.8%, below the 50% response floor), so the declared capacity is not reaching its interception. | +| [`summa`](models/summa) | submitted | SUMMA 4.0.0, a land model that solves the water and energy balances of canopy, snow, soil and aquifer with one implicit solver, run as one lumped HRU from its shipped test-case setup, with the radiation and humidity it needs mocked from each forcing row. | **FAIL (VIOLATION)**, 10 of 25 probes passed. Its water budget closes to rounding and its surface energy budget to 0.2% of its own net radiation; what fails is a latent heat of vaporisation held at its 0 C value, a net radiation of its own that is not the probe's `rn` (every hourly phase, and a drought that warms the surface), a canopy that, as the shipped setup configures SUMMA, intercepts all rain and keeps what freezes near 0 C, holding up to 27 mm of ice against a 2 mm capacity on the probes that score it (also the only failure on the precipitation counterfactual), leaf area that keeps its seasons under constant weather, runoff that follows the melt rather than the rain on one seed, a 5.2% runoff response to turning snow into rain on another, Green-Ampt runoff that appears only at the hourly step, and no human water use, so the prescribed withdrawal on the human-abstraction probe is never taken. It passes `mass/interception-response`: the larger canopy moves water into evaporation and out of runoff, and the paired budget closes. | +| [`cwatm`](models/cwatm) | submitted | CWatM 1.11, IIASA's Community Water Model and an ISIMIP global hydrological model, run on one grid cell with every store it carries reported. | **FAIL (VIOLATION)**, 17 of 22 probes passed. Its budget closes to 0.03%, and its own water-demand module pumps a prescribed withdrawal out of groundwater to within 0.004%. That 0.03% is water its capillary rise creates, a few thousandths of a millimetre a day, and on `mass/extreme-event-closure` it fails 3 to 8 one-day drizzle events of under 0.07 mm per seed, where 0.001 to 0.005 mm more leaves or is stored than fell, against that probe's allowance of 5% of the rain or 0.001 mm, whichever is larger. It also fails on a groundwater reservoir with no dt that drains 24 times too fast at an hourly step (52% of the rain); on a 0.29 mm/day dip after an added storm, where preferential flow turns surface runoff into interflow that leaves through a slower runoff-concentration lag; and on evaporation on wet soil at 0.70 of demand, near the cap its crop coefficients set. The last two are packaging choices as much as results: this package follows the CWatM-Earth-30min template its parameters come from, and with `preferentialFlow = False`, the setting of the pinned model repository's own 30′ templates, both pass. On `mass/interception-response` the paired residual is 1.54 mm on an 8 mm capacity change (19.2%, over the 12.5% effective limit), so the larger canopy does not fully account for the extra interception. | +| [`lisflood`](models/lisflood) | submitted | LISFLOOD 5.0.0, the EC Joint Research Centre's distributed model behind EFAS and GloFAS, stepped through its own Python framework on one representative 5 km cell and reporting every store its own water balance module counts. | **FAIL (ERROR)**, 19 of 22 probes passed. Its budget closes to 1e-13 mm per step; it reports no heat fluxes and no surface temperature, so the five energy probes that need an energy output cannot ask it anything and are N/A. It fails the step probe because its potential infiltration is a pore-space storage multiplied by the step length, so rain falling within hours runs off at the hourly step (13.0% of the rain between PT1H and PT1D). The two ten-year probes take about 90 s per run under amd64 emulation against a 60 s budget, so their rows are ERROR from the host's speed, and those two ERROR rows alone make the verdict FAIL (ERROR). Run outside the limit, `mass/precipitation-counterfactual` passes every criterion. On `mass/human-abstraction`, LISFLOOD's own water-use rule takes the groundwater share in full and the rest only from channel water above an environmental-flow reserve, recording what the channel cannot give as shortage. On this one-cell water region it withdraws 17 to 33% of the prescription, from the sourced reserve to none, and leaves 67 to 83% where the probe allows 5%. It passes `mass/interception-response`, where the paired partition closes and the evaporation response is present. | | [`modflow6`](models/modflow6) | submitted | MODFLOW 6.7.0, the USGS's modular groundwater model: a transient one-river-cell aquifer run through flopy, reporting the RIV and STO budget terms as signed exchange and groundwater storage. Archived as evidential, not gated on. | **PASS**, 3 of 3 probes scored | | `reference_bucket` | exact | conserves water exactly by construction | must pass every probe that can ask it anything; N/A on the six energy probes and the three groundwater-exchange probes that need outputs it does not report | | `reference_exchange_exact` | exact | a lagged-head bookkeeping reference that closes its own groundwater balance exactly and reports genuinely bidirectional exchange | must pass `mass/gw-sw-exchange-consistency`, the only probe that can ask it anything | @@ -183,8 +184,8 @@ the probe cannot ask the declared model interface this question. | `reference_datum_exact` | exact | linear single-aquifer reference with all absolute heads and elevations shifted consistently | must pass `mass/groundwater-datum-invariance` | | `reference_datum_dependent` | broken | positive exchange whose conductance depends on the numerical initial head | caught by `datum_flux_invariance` | | `reference_zero_exchange` | broken | frozen aquifer reporting no exchange | caught by `datum_flux_invariance` | -| [`flex_lumped`](models/flex_lumped) | physical | lumped FLEX/HBV: interception, beta-partitioned unsaturated store, fast and slow reservoirs, geometry-aware triangular lag | must pass every probe that can ask it anything; **PASS**, 25 of 25 | -| [`flex_topo`](models/flex_topo) | physical | FLEX-Topo: plateau, hillslope and wetland units on real Wark fractions sharing one groundwater store | must pass every probe that can ask it anything; **PASS**, 24 of 24 | +| [`flex_lumped`](models/flex_lumped) | physical | lumped FLEX/HBV: interception, beta-partitioned unsaturated store, fast and slow reservoirs, geometry-aware triangular lag | must pass every probe that can ask it anything; **PASS**, 26 of 26 | +| [`flex_topo`](models/flex_topo) | physical | FLEX-Topo: plateau, hillslope and wetland units on real Wark fractions sharing one groundwater store | must pass every probe that can ask it anything; **PASS**, 25 of 25 | | [`sacsma_snow17`](models/sacsma_snow17) | physical | the NWS's SAC-SMA with Snow-17 and a gamma unit hydrograph, ported from the legacy Fortran and checked against it | must pass every probe that can ask it anything; **PASS**, 25 of 25 | | `reference_coupled` | exact | the bucket with snow sublimation and a surface energy budget: every kilogram converted at the latent heat of the phase it actually underwent | must pass every criterion of the three energy-flux probes; supports daily and hourly steps | | `reference_snow_energy` | exact | an energy-balance snowpack: melt bought with `max(0, Rn - H - LE - G) / lambda_f`, liquid held in the pore space and reported as `lwsnl`, cold content carried as an energy deficit and reported as `csnow` | must pass every criterion of `energy/snowmelt-energy-water`; supports daily and hourly steps | @@ -240,6 +241,10 @@ the probe cannot ask the declared model interface this question. | `reference_rating_drift` | broken | derives its stage from a slowly decaying running maximum of discharge (`peak = max(q, 0.997 * peak)` per day), so the gauge ratchets up with each flood far faster than it relaxes, stepping down only a fraction of a percent a day | caught by `rating_monotonic` | | `reference_rating_inverted` | broken | reads the loop backwards, high while the flood is arriving and low once it is leaving: monotone in discharge, so only the loop sees it | caught by `rating_loop` | | `reference_flat_stage` | broken | reports a constant stage, so there is no rating and no loop, only a number that does not vary | caught by `non_degenerate` | +| `reference_interception_blind` | broken | declares `canopy_capacity_mm` and ignores it, so the two capacity variants are bit-identical | caught by `interception_response` (evaporation response floor) | +| `reference_interception_wetter` | broken | a larger capacity raises runoff and pays for it out of a separate store | caught by `interception_response` (runoff may not rise) | +| `reference_interception_drier` | broken | a larger capacity lowers evaporation and parks the water in storage | caught by `interception_response` (evaporation may not fall) | +| `reference_interception_leaky` | broken | returns part of the paired difference as evaporation and drops the rest | caught by `interception_response` (paired closure) | | `reference_uniform_flow` | exact | converts constant effective rainfall to discharge and solves exact rectangular Manning normal depth above the declared bed | must pass `uniform_flow_friction` | | `reference_saint_venant` | exact | advances one-dimensional continuity and momentum with a finite-volume solver from a non-equilibrium state; no normal-depth lookup | must pass `uniform_flow_friction` | | `reference_wrong_roughness` | broken | computes stage with Manning roughness 3% above the declared value, a 5.74% near-boundary residual | caught by `uniform_flow_friction` | @@ -417,7 +422,7 @@ where that conversation happens, before and alongside the issues. ## Status -Suite `0.1.0`, pre-release. Thirty-seven probes, twenty-four mass, seven energy, six momentum, synthetic track only. More +Suite `0.1.0`, pre-release. Thirty-eight probes, twenty-five mass, seven energy, six momentum, synthetic track only. More energy and momentum probes, and the real-data track, are next. The harness runs paired cases and scores labelled regimes; spatial and temporal closure, counterfactual response and invariance are represented in the suite. The roadmap lists the remaining diff --git a/ROADMAP.md b/ROADMAP.md index 9c2b29f0..172e37ad 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -146,6 +146,15 @@ for exactly the abstracted volume (net of return flow). Contributed by Yuanhang Liu. +### `mass/interception-response` · **merged** +The same three-year weather with canopy capacity 2 mm and 10 mm: a larger +canopy must not raise runoff or lower evaporation, the paired budget must +account for the difference, and the evaporation response must exist. +*Discriminates:* models that declare `canopy_capacity_mm` and ignore it, or +allocate the intercepted water the wrong way. + +Contributed by Yuanhang Liu. + ### `mass/exchange-response` · **merged** A declared head-driven exchange must respond to its external head. Once a model closes its budget with `gwex`, that flux is identically minus the diff --git a/docs/writing-a-probe.md b/docs/writing-a-probe.md index 9277a8bf..673cb5a4 100644 --- a/docs/writing-a-probe.md +++ b/docs/writing-a-probe.md @@ -146,6 +146,7 @@ Every one is binary. | `routing_conservation` | the channel store is non-negative and never exceeds `max_lag_days` of the largest recent runoff, plus `min_allowance_mm` | one run | | `lag_time_bounds` | rainfall-to-runoff peak lag lies inside a broad duration-corrected Snyder envelope derived from public catchment geometry | paired runs, a geometry ladder | | `scaling_monotonicity` | peak lag does not materially reverse and grows by a resolvable amount across catchment scales | paired runs, a geometry ladder | +| `interception_response` | a declared static capacity change must not raise runoff or lower evaporation, the paired budget must account for the difference, and the evaporation response must exist; `Delta cap` is read from the two runs' `case.static`, not from a number in the probe | paired runs, two capacity variants | | `rating_monotonic` | stage does not fall against its running maximum as the abscissa rises: equal-count bin medians are taken over the abscissa and the summed running-maximum deficit is compared with an explicit `tolerance` when one is given, and otherwise with 8% of the rating's span. The share is that large because a stage read off a store is hysteretic by construction, and its binned rating dips below its own running maximum by a visible fraction of the span for that reason alone | one run | | `rating_loop` | where the gauge loops against the reach's store, the loop must be small enough to be noise or run the right way: at the same storage the rising limb sits lower than the falling one. A single-valued rating, or a loop below `min_loop_m` with an inconsistent sign across bins, is read as "no loop" and passes | one run | | `uniform_flow_friction` | on each labelled low, medium and high steady plateau, the reported discharge and stage must make Manning friction slope agree with the declared bed slope for the explicit rectangular section; CV and first-to-last-quarter trend gates reject blocks that have not converged | one run, three labelled plateaus | diff --git a/models/reference_interception_blind/ht_adapter.py b/models/reference_interception_blind/ht_adapter.py new file mode 100644 index 00000000..ed68f9bf --- /dev/null +++ b/models/reference_interception_blind/ht_adapter.py @@ -0,0 +1,118 @@ +"""Synthetic negative control for mass/interception-response. + +MODE = "blind". The control declares `canopy_capacity_mm` and consumes it, +but allocates the paired difference the wrong way; it is not physical, and it +exists only to show that `interception_response` catches the named failure. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +MODE = "blind" +MODEL = {"name": "reference_interception_blind", "version": "1.0.0"} +COLUMNS = ["time", "pr", "evspsbl", "mrro", "mrso", "canopy", "gw"] + + +def simulate(forcing, static): + declared = float(static["canopy_capacity_mm"]) + # The base bucket runs at a fixed 2 mm capacity, so the declared capacity + # produces no physical response on its own. The paired difference is then + # exactly the synthetic allocation below, which is what the control is + # built to get wrong. `blind` never reads the declared value at all. + cap = 2.0 + + # This control is specific to mass/interception-response: its scored + # window starts after 365 spinup days, and the synthetic allocation is + # placed in the first 100 scored days so it cannot cancel in state0. + spinup = 365 + adjust_steps = 100 + + soil = 250.0 + canopy = 0.0 + gw = 100.0 + rows = [] + for i, step in enumerate(forcing): + pr = max(float(step["pr"]), 0.0) + pet = max(float(step["pet"]), 0.0) + + intercepted = min(max(cap - canopy, 0.0), pr) + canopy += intercepted + through = pr - intercepted + + canopy_evap = min(canopy, pet) + canopy -= canopy_evap + pet_left = max(pet - canopy_evap, 0.0) + + soil += through + soil_evap = min(soil, pet_left * 0.5 * min(1.0, soil / 100.0)) + soil -= soil_evap + runoff = max(0.0, soil - 200.0) + soil -= runoff + et = canopy_evap + soil_evap + + if spinup <= i < spinup + adjust_steps: + if MODE == "wetter": + # A larger canopy increases runoff, and pays for it out of + # storage. The paired budget still closes; criterion 1 fails. + runoff += 0.008 * declared + et += 0.006 * declared + gw -= 0.014 * declared + elif MODE == "drier": + # A larger canopy suppresses evaporation and parks the water + # in the separate groundwater store. Criterion 2 fails. + cut = 0.006 * declared + take = min(et, cut) + et -= take + gw += take + elif MODE == "leaky": + # Return only part of the water as evaporation and make the + # rest disappear from the separate store. Criterion 3 fails; + # the evaporation response is present, so criterion 4 passes. + et += 0.006 * declared + gw -= 0.008 * declared + + rows.append({ + "time": step["time"], + "pr": float(step["pr"]), + "evspsbl": max(et, 0.0), + "mrro": max(runoff, 0.0), + "mrso": max(soil, 0.0), + "canopy": max(canopy, 0.0), + "gw": max(gw, 0.0), + }) + return rows + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--request", required=True) + args = parser.parse_args() + + request_path = Path(args.request).resolve() + request = json.loads(request_path.read_text()) + io_dir = request_path.parent + + with open(io_dir / request["input"]["forcing"], newline="") as fh: + forcing = list(csv.DictReader(fh)) + static = json.loads((io_dir / request["input"]["static"]).read_text()) + + rows = simulate(forcing, static) + out_path = io_dir / request["output"]["table"] + out_path.parent.mkdir(parents=True, exist_ok=True) + with open(out_path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=COLUMNS) + writer.writeheader() + writer.writerows(rows) + (io_dir / request["output"]["run"]).write_text( + json.dumps({"status": "ok", "model": MODEL, "n_steps": len(rows)}, indent=2) + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/models/reference_interception_blind/model.yaml b/models/reference_interception_blind/model.yaml new file mode 100644 index 00000000..3ab5d18f --- /dev/null +++ b/models/reference_interception_blind/model.yaml @@ -0,0 +1,19 @@ +name: reference_interception_blind +version: "1.0.0" +description: > + Synthetic negative control for mass/interception-response (blind). + It is not physical and exists only to pin the criterion it is declared + against. +authors: ["Yuanhang Liu"] +license: PolyForm-Noncommercial-1.0.0 +runner: subprocess +entrypoint: ["python3", "ht_adapter.py"] +timestep: PT1D +emits: + fluxes: [pr, evspsbl, mrro] + states: [mrso, canopy, gw] +needs_forcing: [pr, pet] +needs_static: [canopy_capacity_mm] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/reference_interception_drier/ht_adapter.py b/models/reference_interception_drier/ht_adapter.py new file mode 100644 index 00000000..14122bd5 --- /dev/null +++ b/models/reference_interception_drier/ht_adapter.py @@ -0,0 +1,118 @@ +"""Synthetic negative control for mass/interception-response. + +MODE = "drier". The control declares `canopy_capacity_mm` and consumes it, +but allocates the paired difference the wrong way; it is not physical, and it +exists only to show that `interception_response` catches the named failure. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +MODE = "drier" +MODEL = {"name": "reference_interception_drier", "version": "1.0.0"} +COLUMNS = ["time", "pr", "evspsbl", "mrro", "mrso", "canopy", "gw"] + + +def simulate(forcing, static): + declared = float(static["canopy_capacity_mm"]) + # The base bucket runs at a fixed 2 mm capacity, so the declared capacity + # produces no physical response on its own. The paired difference is then + # exactly the synthetic allocation below, which is what the control is + # built to get wrong. `blind` never reads the declared value at all. + cap = 2.0 + + # This control is specific to mass/interception-response: its scored + # window starts after 365 spinup days, and the synthetic allocation is + # placed in the first 100 scored days so it cannot cancel in state0. + spinup = 365 + adjust_steps = 100 + + soil = 250.0 + canopy = 0.0 + gw = 100.0 + rows = [] + for i, step in enumerate(forcing): + pr = max(float(step["pr"]), 0.0) + pet = max(float(step["pet"]), 0.0) + + intercepted = min(max(cap - canopy, 0.0), pr) + canopy += intercepted + through = pr - intercepted + + canopy_evap = min(canopy, pet) + canopy -= canopy_evap + pet_left = max(pet - canopy_evap, 0.0) + + soil += through + soil_evap = min(soil, pet_left * 0.5 * min(1.0, soil / 100.0)) + soil -= soil_evap + runoff = max(0.0, soil - 200.0) + soil -= runoff + et = canopy_evap + soil_evap + + if spinup <= i < spinup + adjust_steps: + if MODE == "wetter": + # A larger canopy increases runoff, and pays for it out of + # storage. The paired budget still closes; criterion 1 fails. + runoff += 0.008 * declared + et += 0.006 * declared + gw -= 0.014 * declared + elif MODE == "drier": + # A larger canopy suppresses evaporation and parks the water + # in the separate groundwater store. Criterion 2 fails. + cut = 0.006 * declared + take = min(et, cut) + et -= take + gw += take + elif MODE == "leaky": + # Return only part of the water as evaporation and make the + # rest disappear from the separate store. Criterion 3 fails; + # the evaporation response is present, so criterion 4 passes. + et += 0.006 * declared + gw -= 0.008 * declared + + rows.append({ + "time": step["time"], + "pr": float(step["pr"]), + "evspsbl": max(et, 0.0), + "mrro": max(runoff, 0.0), + "mrso": max(soil, 0.0), + "canopy": max(canopy, 0.0), + "gw": max(gw, 0.0), + }) + return rows + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--request", required=True) + args = parser.parse_args() + + request_path = Path(args.request).resolve() + request = json.loads(request_path.read_text()) + io_dir = request_path.parent + + with open(io_dir / request["input"]["forcing"], newline="") as fh: + forcing = list(csv.DictReader(fh)) + static = json.loads((io_dir / request["input"]["static"]).read_text()) + + rows = simulate(forcing, static) + out_path = io_dir / request["output"]["table"] + out_path.parent.mkdir(parents=True, exist_ok=True) + with open(out_path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=COLUMNS) + writer.writeheader() + writer.writerows(rows) + (io_dir / request["output"]["run"]).write_text( + json.dumps({"status": "ok", "model": MODEL, "n_steps": len(rows)}, indent=2) + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/models/reference_interception_drier/model.yaml b/models/reference_interception_drier/model.yaml new file mode 100644 index 00000000..68d1f6e4 --- /dev/null +++ b/models/reference_interception_drier/model.yaml @@ -0,0 +1,19 @@ +name: reference_interception_drier +version: "1.0.0" +description: > + Synthetic negative control for mass/interception-response (drier). + It is not physical and exists only to pin the criterion it is declared + against. +authors: ["Yuanhang Liu"] +license: PolyForm-Noncommercial-1.0.0 +runner: subprocess +entrypoint: ["python3", "ht_adapter.py"] +timestep: PT1D +emits: + fluxes: [pr, evspsbl, mrro] + states: [mrso, canopy, gw] +needs_forcing: [pr, pet] +needs_static: [canopy_capacity_mm] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/reference_interception_leaky/ht_adapter.py b/models/reference_interception_leaky/ht_adapter.py new file mode 100644 index 00000000..ca24d4cd --- /dev/null +++ b/models/reference_interception_leaky/ht_adapter.py @@ -0,0 +1,118 @@ +"""Synthetic negative control for mass/interception-response. + +MODE = "leaky". The control declares `canopy_capacity_mm` and consumes it, +but allocates the paired difference the wrong way; it is not physical, and it +exists only to show that `interception_response` catches the named failure. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +MODE = "leaky" +MODEL = {"name": "reference_interception_leaky", "version": "1.0.0"} +COLUMNS = ["time", "pr", "evspsbl", "mrro", "mrso", "canopy", "gw"] + + +def simulate(forcing, static): + declared = float(static["canopy_capacity_mm"]) + # The base bucket runs at a fixed 2 mm capacity, so the declared capacity + # produces no physical response on its own. The paired difference is then + # exactly the synthetic allocation below, which is what the control is + # built to get wrong. `blind` never reads the declared value at all. + cap = 2.0 + + # This control is specific to mass/interception-response: its scored + # window starts after 365 spinup days, and the synthetic allocation is + # placed in the first 100 scored days so it cannot cancel in state0. + spinup = 365 + adjust_steps = 100 + + soil = 250.0 + canopy = 0.0 + gw = 100.0 + rows = [] + for i, step in enumerate(forcing): + pr = max(float(step["pr"]), 0.0) + pet = max(float(step["pet"]), 0.0) + + intercepted = min(max(cap - canopy, 0.0), pr) + canopy += intercepted + through = pr - intercepted + + canopy_evap = min(canopy, pet) + canopy -= canopy_evap + pet_left = max(pet - canopy_evap, 0.0) + + soil += through + soil_evap = min(soil, pet_left * 0.5 * min(1.0, soil / 100.0)) + soil -= soil_evap + runoff = max(0.0, soil - 200.0) + soil -= runoff + et = canopy_evap + soil_evap + + if spinup <= i < spinup + adjust_steps: + if MODE == "wetter": + # A larger canopy increases runoff, and pays for it out of + # storage. The paired budget still closes; criterion 1 fails. + runoff += 0.008 * declared + et += 0.006 * declared + gw -= 0.014 * declared + elif MODE == "drier": + # A larger canopy suppresses evaporation and parks the water + # in the separate groundwater store. Criterion 2 fails. + cut = 0.006 * declared + take = min(et, cut) + et -= take + gw += take + elif MODE == "leaky": + # Return only part of the water as evaporation and make the + # rest disappear from the separate store. Criterion 3 fails; + # the evaporation response is present, so criterion 4 passes. + et += 0.006 * declared + gw -= 0.008 * declared + + rows.append({ + "time": step["time"], + "pr": float(step["pr"]), + "evspsbl": max(et, 0.0), + "mrro": max(runoff, 0.0), + "mrso": max(soil, 0.0), + "canopy": max(canopy, 0.0), + "gw": max(gw, 0.0), + }) + return rows + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--request", required=True) + args = parser.parse_args() + + request_path = Path(args.request).resolve() + request = json.loads(request_path.read_text()) + io_dir = request_path.parent + + with open(io_dir / request["input"]["forcing"], newline="") as fh: + forcing = list(csv.DictReader(fh)) + static = json.loads((io_dir / request["input"]["static"]).read_text()) + + rows = simulate(forcing, static) + out_path = io_dir / request["output"]["table"] + out_path.parent.mkdir(parents=True, exist_ok=True) + with open(out_path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=COLUMNS) + writer.writeheader() + writer.writerows(rows) + (io_dir / request["output"]["run"]).write_text( + json.dumps({"status": "ok", "model": MODEL, "n_steps": len(rows)}, indent=2) + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/models/reference_interception_leaky/model.yaml b/models/reference_interception_leaky/model.yaml new file mode 100644 index 00000000..b0bad937 --- /dev/null +++ b/models/reference_interception_leaky/model.yaml @@ -0,0 +1,19 @@ +name: reference_interception_leaky +version: "1.0.0" +description: > + Synthetic negative control for mass/interception-response (leaky). + It is not physical and exists only to pin the criterion it is declared + against. +authors: ["Yuanhang Liu"] +license: PolyForm-Noncommercial-1.0.0 +runner: subprocess +entrypoint: ["python3", "ht_adapter.py"] +timestep: PT1D +emits: + fluxes: [pr, evspsbl, mrro] + states: [mrso, canopy, gw] +needs_forcing: [pr, pet] +needs_static: [canopy_capacity_mm] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/reference_interception_wetter/ht_adapter.py b/models/reference_interception_wetter/ht_adapter.py new file mode 100644 index 00000000..bcd9ba0e --- /dev/null +++ b/models/reference_interception_wetter/ht_adapter.py @@ -0,0 +1,118 @@ +"""Synthetic negative control for mass/interception-response. + +MODE = "wetter". The control declares `canopy_capacity_mm` and consumes it, +but allocates the paired difference the wrong way; it is not physical, and it +exists only to show that `interception_response` catches the named failure. +""" + +from __future__ import annotations + +import argparse +import csv +import json +import sys +from pathlib import Path + +MODE = "wetter" +MODEL = {"name": "reference_interception_wetter", "version": "1.0.0"} +COLUMNS = ["time", "pr", "evspsbl", "mrro", "mrso", "canopy", "gw"] + + +def simulate(forcing, static): + declared = float(static["canopy_capacity_mm"]) + # The base bucket runs at a fixed 2 mm capacity, so the declared capacity + # produces no physical response on its own. The paired difference is then + # exactly the synthetic allocation below, which is what the control is + # built to get wrong. `blind` never reads the declared value at all. + cap = 2.0 + + # This control is specific to mass/interception-response: its scored + # window starts after 365 spinup days, and the synthetic allocation is + # placed in the first 100 scored days so it cannot cancel in state0. + spinup = 365 + adjust_steps = 100 + + soil = 250.0 + canopy = 0.0 + gw = 100.0 + rows = [] + for i, step in enumerate(forcing): + pr = max(float(step["pr"]), 0.0) + pet = max(float(step["pet"]), 0.0) + + intercepted = min(max(cap - canopy, 0.0), pr) + canopy += intercepted + through = pr - intercepted + + canopy_evap = min(canopy, pet) + canopy -= canopy_evap + pet_left = max(pet - canopy_evap, 0.0) + + soil += through + soil_evap = min(soil, pet_left * 0.5 * min(1.0, soil / 100.0)) + soil -= soil_evap + runoff = max(0.0, soil - 200.0) + soil -= runoff + et = canopy_evap + soil_evap + + if spinup <= i < spinup + adjust_steps: + if MODE == "wetter": + # A larger canopy increases runoff, and pays for it out of + # storage. The paired budget still closes; criterion 1 fails. + runoff += 0.008 * declared + et += 0.006 * declared + gw -= 0.014 * declared + elif MODE == "drier": + # A larger canopy suppresses evaporation and parks the water + # in the separate groundwater store. Criterion 2 fails. + cut = 0.006 * declared + take = min(et, cut) + et -= take + gw += take + elif MODE == "leaky": + # Return only part of the water as evaporation and make the + # rest disappear from the separate store. Criterion 3 fails; + # the evaporation response is present, so criterion 4 passes. + et += 0.006 * declared + gw -= 0.008 * declared + + rows.append({ + "time": step["time"], + "pr": float(step["pr"]), + "evspsbl": max(et, 0.0), + "mrro": max(runoff, 0.0), + "mrso": max(soil, 0.0), + "canopy": max(canopy, 0.0), + "gw": max(gw, 0.0), + }) + return rows + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("--request", required=True) + args = parser.parse_args() + + request_path = Path(args.request).resolve() + request = json.loads(request_path.read_text()) + io_dir = request_path.parent + + with open(io_dir / request["input"]["forcing"], newline="") as fh: + forcing = list(csv.DictReader(fh)) + static = json.loads((io_dir / request["input"]["static"]).read_text()) + + rows = simulate(forcing, static) + out_path = io_dir / request["output"]["table"] + out_path.parent.mkdir(parents=True, exist_ok=True) + with open(out_path, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=COLUMNS) + writer.writeheader() + writer.writerows(rows) + (io_dir / request["output"]["run"]).write_text( + json.dumps({"status": "ok", "model": MODEL, "n_steps": len(rows)}, indent=2) + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/models/reference_interception_wetter/model.yaml b/models/reference_interception_wetter/model.yaml new file mode 100644 index 00000000..856630d5 --- /dev/null +++ b/models/reference_interception_wetter/model.yaml @@ -0,0 +1,19 @@ +name: reference_interception_wetter +version: "1.0.0" +description: > + Synthetic negative control for mass/interception-response (wetter). + It is not physical and exists only to pin the criterion it is declared + against. +authors: ["Yuanhang Liu"] +license: PolyForm-Noncommercial-1.0.0 +runner: subprocess +entrypoint: ["python3", "ht_adapter.py"] +timestep: PT1D +emits: + fluxes: [pr, evspsbl, mrro] + states: [mrso, canopy, gw] +needs_forcing: [pr, pet] +needs_static: [canopy_capacity_mm] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/result.csv b/models/result.csv index 30717a1e..6dc053b3 100644 --- a/models/result.csv +++ b/models/result.csv @@ -836,3 +836,13 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-10-05,flex_lumped,1.2.0,0.1.0,subprocess,momentum/wave-celerity-bounds,N/A,INCOMPATIBLE,not run,,model does not declare that it consumes forcing q_in; model does not declare that it consumes static reach_length_m 2026-10-05,sacsma_snow17,1.1.0,0.1.0,subprocess,momentum/wave-celerity-bounds,N/A,INCOMPATIBLE,not run,,model does not declare that it consumes forcing q_in; model does not declare that it consumes static reach_length_m 2026-10-16,wflow_sbm,1.0.4-ht.5,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,"does not report gw_sw_exchange, gw_to_sw, sw_to_gw, gw; model does not declare that it consumes forcing gw_recharge, sw_stage_m; model does not declare that it consumes static area_km2, aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day" +2026-10-23,flex_lumped,1.2.0,0.1.0,subprocess,mass/interception-response,PASS,OK,full record,1135680359 1642632473 2100940,"interception_response: capacity 2 -> 10 mm: evaporation +154.76 mm, runoff -143.44 mm, storage -11.32 mm, residual 0.00 mm of 8 mm (0.0%)" +2026-10-23,flex_topo,1.1.0,0.1.0,subprocess,mass/interception-response,PASS,OK,full record,1135680359 1642632473 2100940,"interception_response: capacity 2 -> 10 mm: evaporation +188.87 mm, runoff -179.03 mm, storage -9.84 mm, residual -0.00 mm of 8 mm (0.0%)" +2026-10-23,dhbv2,0.5.4-hbv2ep100.3,0.1.0,docker,mass/interception-response,N/A,INCOMPATIBLE,not run,,model does not declare that it consumes static canopy_capacity_mm +2026-10-23,google_flood_forecast,0.1.0-828dfc5,0.1.0,docker,mass/interception-response,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrso, canopy; model does not declare that it consumes static canopy_capacity_mm" +2026-10-23,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/interception-response,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, canopy; model does not declare that it consumes static canopy_capacity_mm" +2026-10-23,sacsma_snow17,1.1.0,0.1.0,subprocess,mass/interception-response,N/A,INCOMPATIBLE,not run,,model does not declare that it consumes static canopy_capacity_mm +2026-10-23,cwatm,1.11-5baaadd.4,0.1.0,docker,mass/interception-response,FAIL,VIOLATION,full record,1135680359 1642632473 2100940,"interception_response: capacity 2 -> 10 mm: the paired budget does not account for the difference (-1.54 mm, 19.2%, limit 12.5%)" +2026-10-23,lisflood,5.0.0-onecell.5,0.1.0,docker,mass/interception-response,PASS,OK,1095-day flood event,1135680359 1642632473 2100940,"interception_response: capacity 2 -> 10 mm: evaporation +148.76 mm, runoff -157.24 mm, storage +8.49 mm, residual -0.00 mm of 8 mm (0.0%)" +2026-10-23,summa,4.0.0-f787fa5.4,0.1.0,docker,mass/interception-response,PASS,OK,full record,1135680359 1642632473 2100940,"interception_response: capacity 2 -> 10 mm: evaporation +157.80 mm, runoff -177.11 mm, storage +19.32 mm, residual 0.00 mm of 8 mm (0.0%)" +2026-10-23,wflow_sbm,1.0.4-ht.5,0.1.0,docker,mass/interception-response,FAIL,VIOLATION,full record,1135680359 1642632473 2100940,"interception_response: capacity 2 -> 10 mm: the evaporation response is too small or absent (+0.07 mm, 0.8% of the capacity change, floor 50%)" diff --git a/probes/mass/interception-response/README.md b/probes/mass/interception-response/README.md new file mode 100644 index 00000000..07620d31 --- /dev/null +++ b/probes/mass/interception-response/README.md @@ -0,0 +1,106 @@ +# mass/interception-response + +A larger canopy intercepts more water. That water has to go somewhere: less +runoff, more evaporation, or a reported store. The probe runs the same +three-year weather twice, changing only `static.canopy_capacity_mm` from 2 mm +to 10 mm, and scores where the difference went. + +## The residual + +Write `Delta X = X(capacity_high) - X(capacity_low)` over the scored window, +with fluxes integrated to depths in mm: + +``` +R = Delta evspsbl + Delta mrro + Delta gwex + + Delta (mrso + snw + canopy + gw + channel) [mm] +``` + +A model that conserves water has `R = 0`. The capacity change is +`Delta cap = 8 mm`, and the four checks are: + +1. `max(0, Delta mrro) / Delta cap <= 0.05` — a larger canopy may not raise + runoff; +2. `max(0, -Delta evspsbl) / Delta cap <= 0.05` — it may not lower + evaporation; +3. `abs(R) / Delta cap <= 0.05` — the paired budget must account for the + difference, using every store the model reports; +4. `Delta evspsbl / Delta cap >= 0.5` — the response must exist and must be an + evaporation response. + +The absolute residual floor is 1 mm, so the effective limit is +`max(0.05, 1 mm / Delta cap) = 12.5%` at `Delta cap = 8 mm`. The floor is +there so a numerical residual on a multi-year record does not bind; the +measured exact bucket closes to better than 0.01 mm. + +Criterion 4 is what makes this more than a closure check. A model that +declares `canopy_capacity_mm` and ignores it produces two identical runs: +criteria 1-3 pass trivially and 4 fails. The floor is on evaporation rather +than on the canopy-store peak because a model that routes interception loss +directly to evaporation is physically sound and should not be failed for not +exposing a separate store. + +## Why 2 -> 10 mm + +10 mm is a closed-canopy capacity and is representable by every model that +currently declares the static; LISFLOOD's `lai_for_canopy_capacity` +conversion tops out near 11.7 mm, and 20 mm raises before the native run. The +change is still a factor of five. Native CWatM 1.11, run outside Docker at +2 -> 20 mm on an earlier dry-down seed, moved +63.32 mm of ET and -62.84 mm of +runoff with a 0.784 mm paired residual (reported by @cehw in #161), which is +independent evidence that the design asks a real question. Its row on this +probe is generated through the Docker adapter, not by hand. + +## Measured on the gate seed + +Seed `1135680359`, 1,095 scored days (365 spinup days in front of them): + +| model | Delta sum evspsbl | Delta sum mrro | paired residual | +| --- | ---: | ---: | ---: | +| `reference_bucket` | +82.18 mm | -104.91 mm | 0.00 mm | +| `flex_lumped` | +144.39 mm | -151.57 mm | 0.00 mm | +| `flex_topo` | +180.58 mm | -180.64 mm | 0.00 mm | + +The response magnitudes differ by more than a factor of two, which is why +criterion 4 is a floor and not a magnitude test. All three pass on every gate +seed. + +## Baselines + +`must_pass`: `reference_bucket`, `flex_lumped`, `flex_topo`. + +`must_fail`, one control per criterion: + +| control | fault | caught by | +| --- | --- | --- | +| `reference_interception_blind` | declares the capacity and ignores it | criterion 4 (zero response) | +| `reference_interception_wetter` | larger capacity raises runoff | criterion 1 | +| `reference_interception_drier` | larger capacity lowers evaporation | criterion 2 | +| `reference_interception_leaky` | returns part of the difference as evaporation and drops the rest | criterion 3 | + +Each is listed under `baselines.must_fail_only`, so the gate requires it to +trip `interception_response` alone rather than a contract accident. + +`cwatm` and `summa` also declare `canopy_capacity_mm` and report `canopy`; +they are candidates for `must_pass` once their Docker rows pass the gate. +`sacsma_snow17` and `dhbv2` report `canopy` as an identically zero state but +do not declare the static, so they are `N/A (INCOMPATIBLE)` and their +standing does not move. + +## What this probe does not catch + +- A model that does not declare `canopy_capacity_mm` is not asked; the static + is the opt-in. +- A model whose conduction of the parameter is physically wrong but still + raises evaporation and closes the budget can pass. The probe scores + allocation, not the internal interception formulation. +- The case fixes 2 -> 10 mm. A model that cannot represent 10 mm is out of + scope until a second pair is agreed; the first gate does not silently score + it at a different range. + +## Reproducing + +```bash +ht run --model reference_bucket --probe mass/interception-response +ht run --model reference_interception_leaky --probe mass/interception-response +ht gate --probe mass/interception-response +``` diff --git a/probes/mass/interception-response/generate.py b/probes/mass/interception-response/generate.py new file mode 100644 index 00000000..18e47a4b --- /dev/null +++ b/probes/mass/interception-response/generate.py @@ -0,0 +1,91 @@ +"""Forcing generator for mass/interception-response. + +Two variants of the same three-year temperate record, differing only in the +declared canopy capacity: 2 mm (a sparse canopy) and 10 mm (a closed one). +The weather is drawn once, before the variant branch, so both runs see +byte-identical rain, temperature and potential evaporation. The record ends +with a dry spell long enough that both canopies drain, which means the extra +water a larger canopy intercepts has to show up in the event peaks rather +than in a final store. + +10 mm is chosen as the upper end because it is representable by every model +that currently declares `canopy_capacity_mm`; LISFLOOD's conversion tops out +near 11.7 mm. The change is still a factor of five. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +PERIOD_YEARS = 3 +SPINUP_DAYS = 365 +DRY_DAYS = 90 +N_STEPS = PERIOD_YEARS * 365 + SPINUP_DAYS + +STATIC = { + "area_km2": 250.0, + "soil_capacity_mm": 320.0, + "canopy_capacity_mm": 2.0, + "degree_day_factor_mm_per_C_day": 3.2, + "baseflow_coefficient": 0.006, + "snow_threshold_degC": 0.0, + "latitude_deg": 38.0, +} + +CAPACITY_LOW_MM = 2.0 +CAPACITY_HIGH_MM = 10.0 + + +def _weather(rng: np.random.Generator, n: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Temperate, rain-dominated daily weather with a terminal dry spell.""" + day = np.arange(n) + doy = day % 365 + p_wet = 0.27 + 0.12 * np.cos(2 * np.pi * (doy - 30) / 365) + wet = rng.random(n) < p_wet + pr = np.where(wet, rng.gamma(shape=0.7, scale=14.0, size=n), 0.0) + + seasonal = 15.0 - 7.0 * np.cos(2 * np.pi * (doy - 15) / 365) + noise = np.zeros(n) + innovation = rng.normal(0.0, 2.0, n) + for t in range(1, n): + noise[t] = 0.72 * noise[t - 1] + innovation[t] + tas = seasonal + noise + daylength = 1.0 + 0.35 * np.cos(2 * np.pi * (doy - 172) / 365) + pet = np.maximum(0.0, 0.13 * (tas + 5.0)) * daylength + + # The last dry spell is scored: both canopies have to empty before the + # window ends, so the paired response cannot hide in a final store. + pr = pr.copy() + pr[-DRY_DAYS:] = 0.0 + return pr, tas, pet + + +def _frame(pr, tas, pet) -> pd.DataFrame: + time = pd.date_range("2000-01-01", periods=len(pr), freq="D") + return pd.DataFrame({ + "time": time.strftime("%Y-%m-%d"), + "pr": np.round(pr, 6), + "tas": np.round(tas, 6), + "pet": np.round(pet, 6), + }) + + +def generate(seed: int, variant: str | None = None) -> tuple[pd.DataFrame, dict]: + rng = np.random.default_rng(seed) + pr, tas, pet = _weather(rng, N_STEPS) + static = dict(STATIC) + if variant == "capacity_low": + static["canopy_capacity_mm"] = CAPACITY_LOW_MM + elif variant == "capacity_high": + static["canopy_capacity_mm"] = CAPACITY_HIGH_MM + elif variant is not None: + raise ValueError(f"unknown variant {variant!r}") + return _frame(pr, tas, pet), static + + +if __name__ == "__main__": + frame, static = generate(20261021) + print(f"{N_STEPS} steps; {frame['pr'][:SPINUP_DAYS].sum():.0f} mm rain in spinup, " + f"{frame['pr'][SPINUP_DAYS:].sum():.0f} mm in the scored years; " + f"canopy capacity {static['canopy_capacity_mm']} mm") diff --git a/probes/mass/interception-response/probe.yaml b/probes/mass/interception-response/probe.yaml new file mode 100644 index 00000000..2f847952 --- /dev/null +++ b/probes/mass/interception-response/probe.yaml @@ -0,0 +1,74 @@ +id: mass/interception-response + +title: A larger canopy must intercept more water and give it back as evaporation + +law: mass +track: synthetic +version: 1 + +authors: + - name: Yuanhang Liu + affiliation: Independent Researcher + github: kawh1111 + orcid: 0009-0005-9633-6004 + +description: > + The same three-year weather twice, differing only in the declared canopy + capacity (2 mm against 10 mm). A larger canopy intercepts more water, so + runoff may not rise, evaporation may not fall, and the paired budget must + account for the difference: the extra water has to leave as evaporation or + remain in a reported store. A model that declares the capacity and ignores + it produces two identical runs and fails the response floor, which is what + makes this more than a closure check. + +requires: + fluxes: [pr, evspsbl, mrro] + states: [mrso, canopy] + static: [canopy_capacity_mm] + +case: + generator: generate.py + n_seeds: 3 + timestep: PT1D + period_years: 3 + spinup_days: 365 + max_output_mb: 2.0 + max_runtime_s: 120 + # The budget is taken over the whole scored record, including the terminal + # dry spell that drains both canopies. + min_window_days: 1095 + # The first name is the control. Both variants see identical weather and + # differ only in `static.canopy_capacity_mm`. + variants: [capacity_low, capacity_high] + +criteria: + - interception_response: + parameter: canopy_capacity_mm + low: capacity_low + high: capacity_high + tolerance: 0.05 + response_floor: 0.5 + # A 1 mm absolute floor on the paired residual. The measured exact + # bucket closes to 1e-6 mm; this only keeps a numerical residual from + # binding at the 8 mm capacity change. + floor_mm: 1.0 + +baselines: + must_pass: [reference_bucket, flex_lumped, flex_topo, lisflood, summa] + must_fail: + reference_interception_blind: interception_response + reference_interception_wetter: interception_response + reference_interception_drier: interception_response + reference_interception_leaky: interception_response + # Each control is meant to fail through `interception_response` alone, not + # through a contract accident. + must_fail_only: + - reference_interception_blind + - reference_interception_wetter + - reference_interception_drier + - reference_interception_leaky + +provenance: > + Synthetic. A three-year temperate record with a terminal dry spell, + generated at run time. The two variants share the weather and differ only + in the declared canopy capacity. diff --git a/site/index.html b/site/index.html index 0638bc13..3474578b 100644 --- a/site/index.html +++ b/site/index.html @@ -350,7 +350,7 @@

How it works

HydroTuring does it conserve what physics says it must? - + MASS water budget closes @@ -403,6 +403,8 @@

How it works

datum recharge + + interception ENERGY energy budget @@ -710,18 +712,18 @@

Credit