[PROBE: mass] add aquifer recharge-response ordering

CONTRIBUTORS.md

@@ -27,6 +27,7 @@ long.
 | `mass/extreme-rain` | Zhi Li (CU Boulder) |
 | `mass/runoff-bounds` | Zhi Li (CU Boulder) |
 | `mass/groundwater-datum-invariance` | Zhenjiang Wu (Xi’an Jiaotong University) |
+| `mass/aquifer-recharge-ordering` | Zhenjiang Wu (Institute of Global Environmental Change, Xi’an Jiaotong University) |
 | `mass/area-invariance` | Zhi Li (CU Boulder) |
 | `mass/response-nonnegativity` | Zhi Li (CU Boulder) |
 | `mass/antecedent-monotonicity` | Zhi Li (CU Boulder) |

README.md

@@ -68,17 +68,17 @@ not.
 
 ## The probes
 
-Thirty-five: twenty-three under mass, seven under energy and five under momentum.
+Thirty-six: twenty-four under mass, seven under energy and five 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
 hand-written FLEX models and the NWS's SAC-SMA with Snow-17, must pass
 every probe that can ask them anything; six of the seven energy probes and
-the two groundwater-exchange probes need outputs they do not report and are not
+the three groundwater-exchange probes need outputs they do not report and are not
 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. Twelve of the thirty-five require no model output
+with any probe pull request. Twelve of the thirty-six 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
@@ -108,6 +108,7 @@ unclaimed.
 | [`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/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` |
 | [`mass/gw-sw-exchange-consistency`](probes/mass/gw-sw-exchange-consistency) | mass | A model reporting groundwater-river exchange: do recharge, the net exchange (`gw_sw_exchange`, distinct from `gwex`) and the two signed directional components, and aquifer storage all agree? | `reference_exchange_sign_error` |
 | [`energy/pet-consistency`](probes/energy/pet-consistency) | energy | Evaporation reaches demand when the model's own soil is wettest, stays below it, and falls when the soil is driest. | `reference_thirsty` |
@@ -169,9 +170,12 @@ the probe cannot ask the declared model interface this question.
 | [`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%. |
-| [`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**, 2 of 2 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 two groundwater-exchange probes that need outputs it does not report |
+| [`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 |
+| `reference_recharge_order_exact` | exact | the exact daily solution of one aquifer cell draining to a river, so added recharge is stored, then released | must pass `mass/aquifer-recharge-ordering` |
+| `reference_recharge_blind` | broken | ignores `gw_recharge`, so its own budget no longer closes | caught by `aquifer_recharge_ordering` |
+| `reference_recharge_overshoot` | broken | sends 120% of a recharge day to the river, so its storage falls below the control's | caught by `aquifer_recharge_ordering` |
 | `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` |
@@ -409,7 +413,7 @@ where that conversation happens, before and alongside the issues.
 
 ## Status
 
-Suite `0.1.0`, pre-release. Thirty-five probes, twenty-three mass, seven energy, five momentum, synthetic track only. More
+Suite `0.1.0`, pre-release. Thirty-six probes, twenty-four mass, seven energy, five 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

ROADMAP.md

@@ -97,6 +97,11 @@ The largest storm scaled to ten times: runoff cannot fall, nor exceed the rain a
 ### `mass/runoff-bounds` · **merged**
 Over ten years, is the runoff possible at all? The mass question a runoff-only model has to answer.
 
+### `mass/aquifer-recharge-ordering` · **merged**
+Direct aquifer recharge must preserve groundwater storage and incremental exchange ordering.
+
+Contributed by Zhenjiang Wu (Xi'an Jiaotong University).
+
 ### `mass/groundwater-datum-invariance` · **merged**
 
 A common translation of all absolute groundwater heads and river elevations must leave the exchange trajectory unchanged.

docs/writing-a-probe.md

@@ -107,6 +107,7 @@ Every one is binary.
 | `counterfactual_response` | added or removed water is partitioned, not absorbed; `perturbed` may name one variant or a list, each scored against the control | paired runs |
 | `invariance` | a transform the physics ignores changes nothing | paired runs |
 | `datum_flux_invariance` | a common vertical datum shift leaves groundwater exchange unchanged, with a nonzero activity guard | paired runs |
+| `aquifer_recharge_ordering` | added direct aquifer recharge cannot lower storage or reverse incremental exchange, and is conserved cumulatively | paired runs |
 | `resolution_invariance` | integrated volumes agree between the same weather at two steps | paired runs at different steps |
 | `response_sign` | perturb one driver both ways, hold the rest: each response must point the way physics says, by a real share of the change in demand | paired runs |
 | `causality` | nothing may change before an added storm, and runoff must answer it after | paired runs |
@@ -376,6 +377,9 @@ The reference models available today:
 | `reference_driven_exchange` | a head-driven boundary exchange against a constant catchment head; must pass `mass/exchange-response` | nothing, it must pass |
 | `reference_evolving_exchange` | the same boundary against a catchment head that moves, five-day time constant; must pass | nothing, it must pass |
 | `reference_recharge_exchange` | the evolving boundary on a losing catchment that drains half its runoff through it, storativity 1 mm/m; must pass | nothing, it must pass |
+| `reference_recharge_order_exact` | the exact daily solution of one aquifer cell draining to a river; must pass `mass/aquifer-recharge-ordering` | nothing, it must pass |
+| `reference_recharge_blind` | ignores `gw_recharge`, so its own budget no longer closes | `aquifer_recharge_ordering` |
+| `reference_recharge_overshoot` | sends 120% of a recharge day to the river, so its storage falls below the control's | `aquifer_recharge_ordering` |
 | `reference_noise_sink` | declares the prescribed head and never reads it; its exchange is the day's accounting error | `exchange_response` |
 | `reference_token_exchange` | the noise sink plus a head term a million times too small to be anything but a token | `exchange_response` |
 | `reference_steady_sink` | the same token over a steady sink rather than a noisy one, which a share of the exchange's variation cannot see | `exchange_response` |

models/reference_recharge_blind/ht_adapter.py

@@ -0,0 +1,66 @@
+"""Recharge-ordering reference: one aquifer cell draining to a river.
+
+The three reference_recharge_* models share this file and differ only in MODE:
+
+exact         the exact solution of S_y dh/dt = r + C/A (stage - h) over each
+              day, so the budget closes and the added water drains on
+              tau = S_y A / C, 200 days on this probe's case
+blind         ignores gw_recharge in its storage and exchange, so the run's
+              own budget no longer closes
+overshoot     on a recharge day sends 120% of it to the river, so storage
+              falls below the control's
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import math
+from pathlib import Path
+
+MODE = "blind"
+COLUMNS = ["time", "gw_sw_exchange", "gw_to_sw", "sw_to_gw", "gw"]
+
+
+def simulate(forcing, static, mode=MODE):
+    b = float(static["aquifer_specific_yield"]) * 1000.0  # mm of storage per m of head
+    k = float(static["river_conductance_m2_per_day"]) / (float(static["area_km2"]) * 1.0e6) * 1000.0  # mm/day per m
+    head = float(static["aquifer_initial_head_m"])
+    bottom = float(static["aquifer_bottom_m"])
+    rows = []
+    for row in forcing:
+        stage = float(row["sw_stage_m"])
+        r = float(row["gw_recharge"])
+        if mode == "blind":
+            r = 0.0
+        if mode == "overshoot" and r > 0.0:
+            q = -1.2 * r
+        else:
+            equilibrium = stage + r / k
+            new_head = equilibrium + (head - equilibrium) * math.exp(-k / b)
+            q = b * (new_head - head) - r
+        head += (r + q) / b
+        rows.append({"time": row["time"], "gw_sw_exchange": q, "gw_to_sw": min(q, 0.0), "sw_to_gw": max(q, 0.0), "gw": b * (head - bottom)})
+    return rows
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--request", required=True)
+    path = Path(parser.parse_args().request).resolve()
+    request = json.loads(path.read_text())
+    root = path.parent
+    with open(root / request["input"]["forcing"], newline="") as f:
+        forcing = list(csv.DictReader(f))
+    static = json.loads((root / request["input"]["static"]).read_text())
+    rows = simulate(forcing, static)
+    with open(root / request["output"]["table"], "w", newline="") as f:
+        writer = csv.DictWriter(f, fieldnames=COLUMNS)
+        writer.writeheader()
+        writer.writerows(rows)
+    (root / request["output"]["run"]).write_text(json.dumps({"status": "ok"}))
+
+
+if __name__ == "__main__":
+    main()

models/reference_recharge_blind/model.yaml

@@ -0,0 +1,16 @@
+name: reference_recharge_blind
+version: "0.1.0"
+description: "Recharge-ordering reference that ignores gw_recharge, so its own budget no longer closes."
+authors: [HydroTuring maintainers]
+license: PolyForm-Noncommercial-1.0.0
+runner: subprocess
+entrypoint: [python3, ht_adapter.py]
+timestep: PT1D
+emits:
+  fluxes: [gw_sw_exchange, gw_to_sw, sw_to_gw]
+  states: [gw]
+needs_forcing: [gw_recharge, sw_stage_m]
+needs_static: [area_km2, aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day]
+supports:
+  perturbation: true
+resources: {cpu: 1, memory_gb: 1, gpu: false}

models/reference_recharge_order_exact/ht_adapter.py

@@ -0,0 +1,66 @@
+"""Recharge-ordering reference: one aquifer cell draining to a river.
+
+The three reference_recharge_* models share this file and differ only in MODE:
+
+exact         the exact solution of S_y dh/dt = r + C/A (stage - h) over each
+              day, so the budget closes and the added water drains on
+              tau = S_y A / C, 200 days on this probe's case
+blind         ignores gw_recharge in its storage and exchange, so the run's
+              own budget no longer closes
+overshoot     on a recharge day sends 120% of it to the river, so storage
+              falls below the control's
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import math
+from pathlib import Path
+
+MODE = "exact"
+COLUMNS = ["time", "gw_sw_exchange", "gw_to_sw", "sw_to_gw", "gw"]
+
+
+def simulate(forcing, static, mode=MODE):
+    b = float(static["aquifer_specific_yield"]) * 1000.0  # mm of storage per m of head
+    k = float(static["river_conductance_m2_per_day"]) / (float(static["area_km2"]) * 1.0e6) * 1000.0  # mm/day per m
+    head = float(static["aquifer_initial_head_m"])
+    bottom = float(static["aquifer_bottom_m"])
+    rows = []
+    for row in forcing:
+        stage = float(row["sw_stage_m"])
+        r = float(row["gw_recharge"])
+        if mode == "blind":
+            r = 0.0
+        if mode == "overshoot" and r > 0.0:
+            q = -1.2 * r
+        else:
+            equilibrium = stage + r / k
+            new_head = equilibrium + (head - equilibrium) * math.exp(-k / b)
+            q = b * (new_head - head) - r
+        head += (r + q) / b
+        rows.append({"time": row["time"], "gw_sw_exchange": q, "gw_to_sw": min(q, 0.0), "sw_to_gw": max(q, 0.0), "gw": b * (head - bottom)})
+    return rows
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--request", required=True)
+    path = Path(parser.parse_args().request).resolve()
+    request = json.loads(path.read_text())
+    root = path.parent
+    with open(root / request["input"]["forcing"], newline="") as f:
+        forcing = list(csv.DictReader(f))
+    static = json.loads((root / request["input"]["static"]).read_text())
+    rows = simulate(forcing, static)
+    with open(root / request["output"]["table"], "w", newline="") as f:
+        writer = csv.DictWriter(f, fieldnames=COLUMNS)
+        writer.writeheader()
+        writer.writerows(rows)
+    (root / request["output"]["run"]).write_text(json.dumps({"status": "ok"}))
+
+
+if __name__ == "__main__":
+    main()

models/reference_recharge_order_exact/model.yaml

@@ -0,0 +1,16 @@
+name: reference_recharge_order_exact
+version: "0.1.0"
+description: "Recharge-ordering reference: the exact solution of one aquifer cell draining to a river, so added recharge is stored, then released on the cell's own time constant."
+authors: [HydroTuring maintainers]
+license: PolyForm-Noncommercial-1.0.0
+runner: subprocess
+entrypoint: [python3, ht_adapter.py]
+timestep: PT1D
+emits:
+  fluxes: [gw_sw_exchange, gw_to_sw, sw_to_gw]
+  states: [gw]
+needs_forcing: [gw_recharge, sw_stage_m]
+needs_static: [area_km2, aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day]
+supports:
+  perturbation: true
+resources: {cpu: 1, memory_gb: 1, gpu: false}