diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index b76256b7..a11013d9 100644 --- a/CONTRIBUTORS.md +++ b/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) | diff --git a/README.md b/README.md index e9795ea1..24bf72de 100644 --- a/README.md +++ b/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 diff --git a/ROADMAP.md b/ROADMAP.md index b69a2dc7..443ab3ef 100644 --- a/ROADMAP.md +++ b/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. diff --git a/docs/writing-a-probe.md b/docs/writing-a-probe.md index 93be2903..58b80dba 100644 --- a/docs/writing-a-probe.md +++ b/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` | diff --git a/models/reference_recharge_blind/ht_adapter.py b/models/reference_recharge_blind/ht_adapter.py new file mode 100644 index 00000000..9dc6261f --- /dev/null +++ b/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() diff --git a/models/reference_recharge_blind/model.yaml b/models/reference_recharge_blind/model.yaml new file mode 100644 index 00000000..6a13423a --- /dev/null +++ b/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} diff --git a/models/reference_recharge_order_exact/ht_adapter.py b/models/reference_recharge_order_exact/ht_adapter.py new file mode 100644 index 00000000..f570f3be --- /dev/null +++ b/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() diff --git a/models/reference_recharge_order_exact/model.yaml b/models/reference_recharge_order_exact/model.yaml new file mode 100644 index 00000000..0f9f4679 --- /dev/null +++ b/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} diff --git a/models/reference_recharge_overshoot/ht_adapter.py b/models/reference_recharge_overshoot/ht_adapter.py new file mode 100644 index 00000000..9cb39d2d --- /dev/null +++ b/models/reference_recharge_overshoot/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 = "overshoot" +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() diff --git a/models/reference_recharge_overshoot/model.yaml b/models/reference_recharge_overshoot/model.yaml new file mode 100644 index 00000000..37e39966 --- /dev/null +++ b/models/reference_recharge_overshoot/model.yaml @@ -0,0 +1,16 @@ +name: reference_recharge_overshoot +version: "0.1.0" +description: "Recharge-ordering reference that sends 120% of a recharge day to the river, so its storage falls below the control's." +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} diff --git a/models/result.csv b/models/result.csv index e1827160..78238657 100644 --- a/models/result.csv +++ b/models/result.csv @@ -780,3 +780,13 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-10-01,modflow6,6.7.0-adapter.2,0.1.0,docker,momentum/routing-lag-consistency,N/A,INCOMPLETE,not run,,"does not report mrro; model does not declare that it consumes forcing pr; model does not declare that it consumes static main_channel_length_km, centroid_channel_length_km" 2026-10-01,modflow6,6.7.0-adapter.2,0.1.0,docker,momentum/stage-discharge-monotonic,N/A,INCOMPLETE,not run,,"does not report dis, channel, stage; model does not declare that it consumes static bed_elevation_m" 2026-10-01,modflow6,6.7.0-adapter.2,0.1.0,docker,momentum/uniform-flow-friction-consistency,N/A,INCOMPLETE,not run,,"does not report dis, stage; model does not declare that it consumes static width_m, cross_section_shape, bed_elevation_m, slope, manning_n" +2026-10-07,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/aquifer-recharge-ordering,PASS,OK,full record,500220947 1007173061 1514125175,aquifer_recharge_ordering: added recharge preserves storage and exchange ordering (I=10 mm) +2026-10-07,cwatm,1.11-5baaadd.4,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; 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-07,dhbv2,0.5.4-hbv2ep100.3,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; 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-07,flex_lumped,1.2.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,"does not report gw_sw_exchange, gw_to_sw, sw_to_gw; model does not declare that it consumes forcing gw_recharge, sw_stage_m; model does not declare that it consumes static aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day" +2026-10-07,flex_topo,1.1.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,"does not report gw_sw_exchange, gw_to_sw, sw_to_gw; model does not declare that it consumes forcing gw_recharge, sw_stage_m; model does not declare that it consumes static aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day" +2026-10-07,google_flood_forecast,0.1.0-828dfc5,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-07,lisflood,5.0.0-onecell.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; 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-07,sacsma_snow17,1.1.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,"does not report gw_sw_exchange, gw_to_sw, sw_to_gw; model does not declare that it consumes forcing gw_recharge, sw_stage_m; model does not declare that it consumes static aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day" +2026-10-07,summa,4.0.0-f787fa5.4,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; 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-07,wflow_sbm,1.0.4-ht.4,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" diff --git a/probes/mass/aquifer-recharge-ordering/README.md b/probes/mass/aquifer-recharge-ordering/README.md new file mode 100644 index 00000000..e83ee8bb --- /dev/null +++ b/probes/mass/aquifer-recharge-ordering/README.md @@ -0,0 +1,78 @@ +# mass/aquifer-recharge-ordering + +Two runs of the same passive aquifer under the same seeded river stage differ +only by a direct `gw_recharge` pulse of 5 mm/day on the first two scored days. +The question is the groundwater comparison principle (#135): water added to an +aquifer whose outflow does not fall as its storage rises can never leave it +with less water, can never make the exchange less favourable to the river, and +has to be accounted for. + +## The case + +One aquifer under 1 km² with specific yield 0.2, its head at the river's +stage, draining to the river through a conductance of 1000 m²/day, the value +`mass/groundwater-datum-invariance` uses. The record is 18 scored days after +2 days of spinup. At this conductance the pulse raises the exchange toward the +river by about 0.8 mm over the record, eight times the tolerance, so a reversal +of it is visible; at 100 m²/day it was 0.08 mm and no reversal could be seen. + +## `aquifer_recharge_ordering` + +Let `I_n` be the cumulative added recharge, `D_n` the cumulative incremental +exchange toward the river (control minus perturbed `gw_sw_exchange`, which is +positive into the aquifer), `B_n` the cumulative change in a declared +`gw_boundary` (positive into the aquifer) and `ΔS_n` the paired storage +difference. Three one-sided errors are scored: + +* `E_S = max_n [-ΔS_n]_+`, the largest fall of storage below the control; +* `E_Q = Σ_i [-d_i]_+ Δt`, the exchange increment sent the wrong way; and +* `E_B = max_n |ΔS_n + D_n - B_n - I_n|`, the worst running partition residual. + +`E_B` largely restates `groundwater_balance`, which this probe also checks on +both runs: if each run closes its own budget, the paired one closes to within +the two runs' tolerances, which can add to slightly more than `ε`. The probe's +independent content is `E_S` and `E_Q`. + +All three share one tolerance, `ε = max(0.01 · I_*, 1e-3 mm, min(1e-5 · S_max, +0.05 mm))`, where `I_*` is the total added recharge (10 mm, so ε = 0.1 mm) and +`S_max` the largest absolute storage either run reports. The scale comes from +the intervention rather than the background storage; the capped third term +absorbs rounding in a large absolute `gw` without letting a datum buy unlimited +slack, as in `groundwater_balance`. + +## A declared boundary + +A model with a boundary on the aquifer, a GHB or a well, reports it as +`gw_boundary`. Both `groundwater_balance` and the paired budget credit it into +the aquifer. A steady boundary cancels between the runs; a head-dependent one +drains part of the pulse, and the order still holds because that flux too does +not fall as the head rises. + +## What this probe does not catch + +The order is non-strict. A zero exchange increment passes, as it must for an +exchange branch that is genuinely insensitive, so a model whose exchange +ignores its own aquifer head stores the whole pulse, never releases it, and +passes. Asking that the water drain back out needs the aquifer's response +time. For a single store that is `S_y · A / C`, from the case's own +attributes, but for a distributed aquifer it depends on a transmissivity and a +geometry the case does not supply: MODFLOW 6's 10 × 10 grid at K = 1 m/day +drains on a time scale of thousands of days, and a floor built on +`S_y · A / C` failed it. `tests/test_recharge_ordering.py` pins the gap, so a +recovery check added later shows up there as a change of scope. + +## Baselines + +The four standard catchment references have no direct aquifer recharge, so +`must_pass` is the exact reference alone. MODFLOW 6, a physical model with this +process, is archived as passing in `models/result.csv`. + +| Reference | What it does | What catches it | +| --- | --- | --- | +| `reference_recharge_order_exact` | the exact daily solution of one aquifer cell draining to the river | nothing, it must pass | +| `reference_recharge_blind` | ignores `gw_recharge`, so its own budget no longer closes | `aquifer_recharge_ordering` (`E_B`), and `groundwater_balance` | +| `reference_recharge_overshoot` | sends 120% of a recharge day to the river | `aquifer_recharge_ordering`: `E_S` at about 20 times the tolerance, and `E_Q` at about 1.5 as the river refills the deficit | + +A reversed exchange increment that keeps each run's budget closed trips `E_Q` +alone. No reference fakes it; `tests/test_recharge_ordering.py` builds it from +the exact reference and checks that `E_Q` is the arm that fails. diff --git a/probes/mass/aquifer-recharge-ordering/generate.py b/probes/mass/aquifer-recharge-ordering/generate.py new file mode 100644 index 00000000..d4974f74 --- /dev/null +++ b/probes/mass/aquifer-recharge-ordering/generate.py @@ -0,0 +1,47 @@ +"""A passive aquifer under a fixed, seeded river stage, with and without a +direct recharge pulse on the first two scored days. + +The river conductance is the one mass/groundwater-datum-invariance uses. At +100 m2/day the pulse raised the exchange by about 0.08 mm over the record, +under the 0.1 mm tolerance, so no reversal of it could be seen; at 1000 m2/day +a reversed increment is 0.8 mm, eight times the tolerance (#152). +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +N_STEPS = 20 +PULSE_ROWS = slice(2, 4) +PULSE_MM_PER_DAY = 5.0 +VARIANTS = ("control", "recharge_added") +STATIC = { + "area_km2": 1.0, + "aquifer_specific_yield": 0.20, + "aquifer_storage_coefficient": 0.002, + "aquifer_initial_head_m": 100.0, + "aquifer_top_m": 110.0, + "aquifer_bottom_m": 90.0, + "river_conductance_m2_per_day": 1000.0, + "river_bottom_offset_m": 1.0, +} + + +def generate(seed: int, variant: str = "control"): + if variant not in VARIANTS: + raise ValueError(f"unknown variant {variant!r}; expected {VARIANTS}") + rng = np.random.default_rng(seed) + # Drawn before the variant matters, so both runs see the same river. + stage = 100.0 + 0.01 * rng.normal(size=N_STEPS) + recharge = np.zeros(N_STEPS) + if variant == "recharge_added": + recharge[PULSE_ROWS] = PULSE_MM_PER_DAY + forcing = pd.DataFrame( + { + "time": pd.date_range("2000-01-01", periods=N_STEPS, freq="D").strftime("%Y-%m-%d"), + "gw_recharge": recharge, + "sw_stage_m": np.round(stage, 6), + } + ) + return forcing, dict(STATIC) diff --git a/probes/mass/aquifer-recharge-ordering/probe.yaml b/probes/mass/aquifer-recharge-ordering/probe.yaml new file mode 100644 index 00000000..fdb336c1 --- /dev/null +++ b/probes/mass/aquifer-recharge-ordering/probe.yaml @@ -0,0 +1,66 @@ +id: mass/aquifer-recharge-ordering +title: Aquifer recharge preserves storage and exchange ordering +law: mass +track: synthetic +version: 1 +authors: + - name: Zhenjiang Wu + affiliation: Institute of Global Environmental Change, Xi’an Jiaotong University, Xi’an 710049, China + orcid: 0000-0002-5329-4917 + github: JustinLoiii +description: > + Identical passive aquifer runs differ only by a nonnegative direct recharge + pulse. The added recharge must never lower groundwater storage, must not make + net exchange less favorable to the river, and must be partitioned between + storage and cumulative incremental exchange without loss or creation. +requires: + fluxes: [gw_sw_exchange, gw_to_sw, sw_to_gw] + states: [gw] + forcing: [gw_recharge, sw_stage_m] + static: [area_km2, aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day] +case: + generator: generate.py + n_seeds: 3 + timestep: PT1D + period_days: 18 + spinup_days: 2 + min_window_days: 10 + max_output_mb: 2.0 + max_runtime_s: 60 + variants: [control, recharge_added] +criteria: + - aquifer_recharge_ordering: + control: control + perturbed: recharge_added + recharge: gw_recharge + exchange: gw_sw_exchange + storage: gw + rel_tol: 0.01 + abs_tol_mm: 1.0e-3 + sources: [gw_boundary] + - groundwater_balance: + recharge: gw_recharge + exchange: gw_sw_exchange + storage: gw + sources: [gw_boundary] + rel_tol: 0.01 + abs_tol_mm: 2.0e-4 + all_variants: true + - exchange_components: + rel_tol: 1.0e-6 + abs_tol: 1.0e-4 + all_variants: true + - state_bounds: + gw: [0.0, null] + all_variants: true +baselines: + # The standard catchment references do not expose a direct aquifer source. + # MODFLOW 6 is the independent physical baseline and is run with paired + # perturbations for this process. + must_pass: [reference_recharge_order_exact] + must_fail: + reference_recharge_blind: aquifer_recharge_ordering + reference_recharge_overshoot: aquifer_recharge_ordering +provenance: > + Synthetic deterministic forcing based on the accepted proposal in issue + https://github.com/Flood-Lab/HydroTuring/issues/135. diff --git a/site/index.html b/site/index.html index 48298bd3..3be98eef 100644 --- a/site/index.html +++ b/site/index.html @@ -350,58 +350,60 @@

How it works

HydroTuring does it conserve what physics says it must? - + MASS water budget closes closure - - ungauged - - event budget - - resolution - - warming - - causality - - dry-down - - steady state - - spin-up cycle - - drift - - extreme rain - - bounds - - area - - nonnegative - - antecedent - - phase - - time origin - - rain ±20% - - withdrawal - - snowpack - - gw exchange - - exchange - - datum - + + ungauged + + event budget + + resolution + + warming + + causality + + dry-down + + steady state + + spin-up cycle + + drift + + extreme rain + + bounds + + area + + nonnegative + + antecedent + + phase + + time origin + + rain ±20% + + withdrawal + + snowpack + + gw exchange + + exchange + + datum + + recharge + ENERGY energy budget closes @@ -422,7 +424,7 @@

How it works

snowmelt probes wanted - + MOMENTUM routing conserves momentum @@ -438,7 +440,7 @@

How it works

subcritical probes wanted - seed → case → run → one bit + seed → case → run → one bit seeded generator @@ -706,7 +708,7 @@

Credit