[PROBE: mass] interception-response: a larger canopy must return the water as evaporation

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

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

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 |

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

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}

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

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}