From f5ba6069514e05694d7a6c0857376a5490ef3513 Mon Sep 17 00:00:00 2001 From: Zhenjiang Wu <852356187@qq.com> Date: Mon, 28 Sep 2026 12:02:38 +0800 Subject: [PATCH 1/5] PROBE: aquifer recharge ordering --- CITATION.cff | 4 ++ CONTRIBUTORS.md | 1 + README.md | 1 + ROADMAP.md | 4 ++ docs/writing-a-probe.md | 1 + models/reference_recharge_blind/ht_adapter.py | 20 ++++++ models/reference_recharge_blind/model.yaml | 16 +++++ .../ht_adapter.py | 20 ++++++ .../reference_recharge_order_exact/model.yaml | 16 +++++ .../ht_adapter.py | 20 ++++++ .../reference_recharge_overshoot/model.yaml | 16 +++++ models/result.csv | 10 +++ .../mass/aquifer-recharge-ordering/README.md | 7 +++ .../aquifer-recharge-ordering/generate.py | 15 +++++ .../mass/aquifer-recharge-ordering/probe.yaml | 61 +++++++++++++++++++ site/index.html | 7 +++ src/hydroturing/criteria/__init__.py | 1 + src/hydroturing/criteria/recharge_ordering.py | 25 ++++++++ src/hydroturing/spec.py | 3 + tests/test_report_detail.py | 1 + 20 files changed, 249 insertions(+) create mode 100644 models/reference_recharge_blind/ht_adapter.py create mode 100644 models/reference_recharge_blind/model.yaml create mode 100644 models/reference_recharge_order_exact/ht_adapter.py create mode 100644 models/reference_recharge_order_exact/model.yaml create mode 100644 models/reference_recharge_overshoot/ht_adapter.py create mode 100644 models/reference_recharge_overshoot/model.yaml create mode 100644 probes/mass/aquifer-recharge-ordering/README.md create mode 100644 probes/mass/aquifer-recharge-ordering/generate.py create mode 100644 probes/mass/aquifer-recharge-ordering/probe.yaml create mode 100644 src/hydroturing/criteria/recharge_ordering.py diff --git a/CITATION.cff b/CITATION.cff index 57b001bb..d0b18bef 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -79,6 +79,10 @@ authors: family-names: Wang affiliation: University of Illinois Urbana-Champaign orcid: "https://orcid.org/0009-0001-7485-4696" + - given-names: Zhenjiang + family-names: Wu + affiliation: Institute of Global Environmental Change, Xi’an Jiaotong University, Xi’an 710049, China + orcid: "https://orcid.org/0000-0002-5329-4917" - given-names: Mofan family-names: Zhang affiliation: Department of Civil and Environmental Engineering, Stanford University, Stanford, CA, USA diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 25716747..53dd0974 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -26,6 +26,7 @@ long. | `mass/multi-decadal-drift` | Bing Li (Independent Researcher) | | `mass/extreme-rain` | Zhi Li (CU Boulder) | | `mass/runoff-bounds` | Zhi Li (CU Boulder) | +| `mass/aquifer-recharge-ordering` | Zhenjiang Wu (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 e07a5903..58b84501 100644 --- a/README.md +++ b/README.md @@ -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 | Added direct aquifer recharge preserves storage and incremental exchange ordering. | `reference_recharge_blind`, `reference_recharge_overshoot` | | [`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` | | [`energy/latent-heat-et-consistency`](probes/energy/latent-heat-et-consistency) | energy | The evaporation a model reports as water and the evaporation implied by the latent heat it reports: are they the same evaporation? | `reference_two_head`, `reference_constant_lambda`, `reference_sublimation_blind`, `reference_energy_leak` | diff --git a/ROADMAP.md b/ROADMAP.md index 7ff3f787..3c2f3698 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -97,6 +97,10 @@ 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. + ### `mass/area-invariance` · **merged** The same weather on a ten times larger catchment: every depth identical. diff --git a/docs/writing-a-probe.md b/docs/writing-a-probe.md index eeb7876d..cf690376 100644 --- a/docs/writing-a-probe.md +++ b/docs/writing-a-probe.md @@ -106,6 +106,7 @@ Every one is binary. | `regime_transfer` | closure holds out of range as well as in range | labelled stretches | | `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 | +| `aquifer_recharge_ordering` | added direct 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 | diff --git a/models/reference_recharge_blind/ht_adapter.py b/models/reference_recharge_blind/ht_adapter.py new file mode 100644 index 00000000..119fc79d --- /dev/null +++ b/models/reference_recharge_blind/ht_adapter.py @@ -0,0 +1,20 @@ +from __future__ import annotations +import argparse,csv,json +from pathlib import Path +COLUMNS=["time","gw_sw_exchange","gw_to_sw","sw_to_gw","gw"] +def simulate(forcing,static,mode="exact"): + sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] + for row in forcing: + stage=float(row["sw_stage_m"]); r=float(row["gw_recharge"]); q=k*(stage-head); new=head+(r+q)/b; q=b*(new-head)-r + if mode=="blind": r_effect=0.; new=head+(r_effect+q)/b; q=b*(new-head)-r_effect + elif mode=="overshoot" and r>0: q=-1.2*r; new=head+(r+q)/b + head=new; gw=b*head + rows.append({"time":row["time"],"gw_sw_exchange":q,"gw_to_sw":min(q,0.),"sw_to_gw":max(q,0.),"gw":gw}) + return rows +def main(): + p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent + with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) + static=json.loads((d/req["input"]["static"]).read_text()); name=Path(__file__).parent.name; mode="blind" if "blind" in name else ("overshoot" if "overshoot" in name else "exact"); rows=simulate(forcing,static,mode) + with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) + (d/req["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..fc00e11d --- /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. +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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] +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..119fc79d --- /dev/null +++ b/models/reference_recharge_order_exact/ht_adapter.py @@ -0,0 +1,20 @@ +from __future__ import annotations +import argparse,csv,json +from pathlib import Path +COLUMNS=["time","gw_sw_exchange","gw_to_sw","sw_to_gw","gw"] +def simulate(forcing,static,mode="exact"): + sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] + for row in forcing: + stage=float(row["sw_stage_m"]); r=float(row["gw_recharge"]); q=k*(stage-head); new=head+(r+q)/b; q=b*(new-head)-r + if mode=="blind": r_effect=0.; new=head+(r_effect+q)/b; q=b*(new-head)-r_effect + elif mode=="overshoot" and r>0: q=-1.2*r; new=head+(r+q)/b + head=new; gw=b*head + rows.append({"time":row["time"],"gw_sw_exchange":q,"gw_to_sw":min(q,0.),"sw_to_gw":max(q,0.),"gw":gw}) + return rows +def main(): + p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent + with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) + static=json.loads((d/req["input"]["static"]).read_text()); name=Path(__file__).parent.name; mode="blind" if "blind" in name else ("overshoot" if "overshoot" in name else "exact"); rows=simulate(forcing,static,mode) + with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) + (d/req["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..8253ca90 --- /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. +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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] +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..119fc79d --- /dev/null +++ b/models/reference_recharge_overshoot/ht_adapter.py @@ -0,0 +1,20 @@ +from __future__ import annotations +import argparse,csv,json +from pathlib import Path +COLUMNS=["time","gw_sw_exchange","gw_to_sw","sw_to_gw","gw"] +def simulate(forcing,static,mode="exact"): + sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] + for row in forcing: + stage=float(row["sw_stage_m"]); r=float(row["gw_recharge"]); q=k*(stage-head); new=head+(r+q)/b; q=b*(new-head)-r + if mode=="blind": r_effect=0.; new=head+(r_effect+q)/b; q=b*(new-head)-r_effect + elif mode=="overshoot" and r>0: q=-1.2*r; new=head+(r+q)/b + head=new; gw=b*head + rows.append({"time":row["time"],"gw_sw_exchange":q,"gw_to_sw":min(q,0.),"sw_to_gw":max(q,0.),"gw":gw}) + return rows +def main(): + p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent + with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) + static=json.loads((d/req["input"]["static"]).read_text()); name=Path(__file__).parent.name; mode="blind" if "blind" in name else ("overshoot" if "overshoot" in name else "exact"); rows=simulate(forcing,static,mode) + with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) + (d/req["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..d15c2ae5 --- /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. +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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] +supports: + perturbation: true +resources: {cpu: 1, memory_gb: 1, gpu: false} diff --git a/models/result.csv b/models/result.csv index b10d9eb0..1e04d306 100644 --- a/models/result.csv +++ b/models/result.csv @@ -726,3 +726,13 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-09-23,flex_lumped,1.2.0,0.1.0,subprocess,mass/spinup-cycle-invariance,PASS,OK,full record,107008012 613960126 1120912240,spinup_cycle_invariance: the repeated forcing reaches the same evaluation-year response for all 3 selected spin-up cycles (worst short<->long:gw departure 0.00%) 2026-09-23,flex_topo,1.1.0,0.1.0,subprocess,mass/spinup-cycle-invariance,PASS,OK,full record,107008012 613960126 1120912240,spinup_cycle_invariance: the repeated forcing reaches the same evaluation-year response for all 3 selected spin-up cycles (worst short<->plus3:channel departure 0.00%) 2026-09-23,sacsma_snow17,1.1.0,0.1.0,subprocess,mass/spinup-cycle-invariance,PASS,OK,full record,107008012 613960126 1120912240,spinup_cycle_invariance: the repeated forcing reaches the same evaluation-year response for all 3 selected spin-up cycles (worst short<->long:gw departure 0.00%) +2026-09-28,cwatm,1.11-5baaadd.4,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,dhbv2,0.5.4-hbv2ep100.3,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,flex_lumped,1.2.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,flex_topo,1.1.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,google_flood_forecast,0.1.0-828dfc5,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,lisflood,5.0.0-onecell.5,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,modflow6,6.7.0-adapter.1,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,sacsma_snow17,1.1.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,summa,4.0.0-f787fa5.4,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-28,wflow_sbm,1.0.4-ht.4,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases diff --git a/probes/mass/aquifer-recharge-ordering/README.md b/probes/mass/aquifer-recharge-ordering/README.md new file mode 100644 index 00000000..1f83444c --- /dev/null +++ b/probes/mass/aquifer-recharge-ordering/README.md @@ -0,0 +1,7 @@ +# mass/aquifer-recharge-ordering + +The paired case adds only direct `gw_recharge` to the same passive aquifer. The +`aquifer_recharge_ordering` criterion evaluates storage ordering, the sign of +the incremental exchange, and the cumulative partition identity using the +added recharge volume as its scale. A plateau with zero incremental exchange +is valid; a recharge-blind or nonmonotone update is not. diff --git a/probes/mass/aquifer-recharge-ordering/generate.py b/probes/mass/aquifer-recharge-ordering/generate.py new file mode 100644 index 00000000..9b121e4d --- /dev/null +++ b/probes/mass/aquifer-recharge-ordering/generate.py @@ -0,0 +1,15 @@ +from __future__ import annotations +import numpy as np +import pandas as pd +N_STEPS=20 +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,"river_conductance_m2_per_day":100.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) + stage=100.0+0.01*rng.normal(size=N_STEPS) + recharge=np.zeros(N_STEPS) + if variant=="recharge_added": recharge[2:4]=5.0 + 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..ac0aaa7d --- /dev/null +++ b/probes/mass/aquifer-recharge-ordering/probe.yaml @@ -0,0 +1,61 @@ +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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] +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 + - groundwater_balance: + recharge: gw_recharge + exchange: gw_sw_exchange + storage: gw + 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: + 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 1d852076..baafc2af 100644 --- a/site/index.html +++ b/site/index.html @@ -317,6 +317,7 @@

How it works

gw exchange exchange + recharge order ENERGY energy budget @@ -660,6 +661,7 @@

Credit

"infra.probe.snowpack":"snowpack", "infra.probe.ungauged":"ungauged", "infra.probe.exchange":"exchange", + "infra.probe.recharge":"recharge order", "infra.probe.gwexchange":"gw exchange", "infra.probe.pet":"PET demand", "infra.probe.latent":"latent heat", @@ -742,6 +744,7 @@

Credit

mass/extreme-rainmassThe largest storm scaled to ten times: runoff cannot fall, nor exceed the rain added.merged mass/runoff-boundsmassOver ten years, is the runoff possible at all: at least rain minus demand minus storage, at most rain plus storage? Runoff only.merged mass/area-invariancemassThe same weather on a ten times larger catchment: every depth identical.merged +mass/aquifer-recharge-orderingmassThe same weather on a ten times larger catchment: every depth identical.merged mass/response-nonnegativitymassOne storm added: from that day on, runoff may never be lower than without it.merged mass/antecedent-monotonicitymassThe same storm after a dry month and a wet one: more runoff from the wet one, no more than the extra water.merged mass/phase-counterfactualmassThe same water as rain instead of snow: timing moves, the integrated volumes do not.merged @@ -883,6 +886,7 @@

Credit

"infra.probe.snowpack":"manto nival", "infra.probe.ungauged":"sin aforo", "infra.probe.exchange":"intercambio", + "infra.probe.recharge":"orden de recarga", "infra.probe.gwexchange":"intercambio subt.", "infra.probe.pet":"demanda", "infra.probe.latent":"calor latente", @@ -965,6 +969,7 @@

Credit

mass/extreme-rainmasaLa mayor tormenta multiplicada por diez: la escorrentía no puede bajar ni superar la lluvia añadida.integrada mass/runoff-boundsmasaEn diez años, ¿es posible esa escorrentía? Al menos lluvia menos demanda menos almacenamiento; a lo sumo lluvia más almacenamiento. Solo con caudal.integrada mass/area-invariancemasaEl mismo tiempo sobre una cuenca diez veces mayor: cada lámina idéntica.integrada +mass/aquifer-recharge-orderingmasaEl mismo tiempo sobre una cuenca diez veces mayor: cada lámina idéntica.integrada mass/response-nonnegativitymasaUna tormenta añadida: desde ese día la escorrentía nunca puede ser menor que sin ella.integrada mass/antecedent-monotonicitymasaLa misma tormenta tras un mes seco y tras uno húmedo: más escorrentía en el húmedo, y no más que el agua extra.integrada mass/phase-counterfactualmasaLa misma agua como lluvia en vez de nieve: cambia el momento, no los volúmenes integrados.integrada @@ -1106,6 +1111,7 @@

Credit

"infra.probe.snowpack":"积雪收支", "infra.probe.ungauged":"无资料", "infra.probe.exchange":"对外交换", + "infra.probe.recharge":"补给顺序", "infra.probe.gwexchange":"地下水交换", "infra.probe.pet":"蒸散需求", "infra.probe.latent":"潜热一致", @@ -1188,6 +1194,7 @@

Credit

mass/extreme-rain质量最大暴雨放大十倍：径流不得减少，也不得超过所加降雨。已合并 mass/runoff-bounds质量十年之内，这样的径流是否可能：至少为降雨减去需求再减去蓄水，至多为降雨加蓄水？仅凭径流即可判定。已合并 mass/area-invariance质量同样的天气，流域面积放大十倍：每一项深度都必须完全相同。已合并 +mass/aquifer-recharge-ordering质量直接补给不能降低地下水储量或逆转交换增量。已合并 mass/response-nonnegativity质量加入一场暴雨：自那天起，径流在任何一天都不得低于未加雨的情形。已合并 mass/antecedent-monotonicity质量同一场暴雨，分别落在干月与湿月之后：湿的流域径流更多，且不超过多出的水量。已合并 mass/phase-counterfactual质量同样的水以雨而非雪的形式落下：时间会变，积分体积不变。已合并 diff --git a/src/hydroturing/criteria/__init__.py b/src/hydroturing/criteria/__init__.py index 7167ed31..cbc624c3 100644 --- a/src/hydroturing/criteria/__init__.py +++ b/src/hydroturing/criteria/__init__.py @@ -25,6 +25,7 @@ human, limits, radiation, + recharge_ordering, rating, regime, response, diff --git a/src/hydroturing/criteria/recharge_ordering.py b/src/hydroturing/criteria/recharge_ordering.py new file mode 100644 index 00000000..d37d3a66 --- /dev/null +++ b/src/hydroturing/criteria/recharge_ordering.py @@ -0,0 +1,25 @@ +from __future__ import annotations +import numpy as np +from hydroturing.criteria.base import FAIL,PASS,CriterionResult,criterion,make_window +from hydroturing.protocol import RunResult +from hydroturing.spec import ProbeSpec +@criterion("aquifer_recharge_ordering",paired=True) +def aquifer_recharge_ordering(runs:dict[str,RunResult],probe:ProbeSpec,params:dict)->CriterionResult: + c=str(params.get("control","control")); p=str(params.get("perturbed","recharge_added")) + if c not in runs or p not in runs: raise ValueError("aquifer_recharge_ordering needs control and perturbed variants") + wc=make_window(runs[c],probe); wp=make_window(runs[p],probe) + recharge=str(params.get("recharge","gw_recharge")); exchange=str(params.get("exchange","gw_sw_exchange")); storage=str(params.get("storage","gw")) + for w,n in ((wc,"control"),(wp,"perturbed")): + if recharge not in w.forcing.columns: raise ValueError(f"{n} forcing needs '{recharge}'") + for v in (exchange,storage): + if v not in w.table.columns: raise ValueError(f"{n} result needs '{v}'") + if len(wc.table)!=len(wp.table): raise ValueError("paired recharge runs have different scored lengths") + added=wp.volume(wp.forcing[recharge])-wc.volume(wc.forcing[recharge]); I=np.cumsum(added); Istar=float(I[-1]) + if Istar<=0: raise ValueError("aquifer_recharge_ordering needs positive added recharge") + ds=wp.table[storage].to_numpy(float)-wc.table[storage].to_numpy(float) + ds-=float(wp.state0[storage])-float(wc.state0[storage]) + d=wc.volume(wc.table[exchange])-wp.volume(wp.table[exchange]); D=np.cumsum(d) + es=float(np.max(np.maximum(-ds,0.0))); eq=float(np.sum(np.maximum(-d,0.0))); eb=float(np.max(np.abs(ds+D-I))) + eps=max(float(params.get("rel_tol",0.01))*Istar,float(params.get("abs_tol_mm",1e-3))) + ratios={"storage":es/eps,"exchange":eq/eps,"budget":eb/eps}; worst_name,worst=max(ratios.items(),key=lambda x:x[1]); ok=worst<=1.0 + return CriterionResult(name="aquifer_recharge_ordering",status=PASS if ok else FAIL,value=worst,threshold=1.0,message=(f"added recharge preserves storage and exchange ordering (I={Istar:.6g} mm)" if ok else f"recharge ordering fails {worst_name}: {max(es,eq,eb):.6g} mm exceeds {eps:.6g} mm"),diagnostics={"added_recharge_mm":Istar,"epsilon_mm":eps,"E_S_mm":es,"E_Q_mm":eq,"E_B_mm":eb,"ratios":ratios}) diff --git a/src/hydroturing/spec.py b/src/hydroturing/spec.py index e4ca6dd4..6ba49341 100644 --- a/src/hydroturing/spec.py +++ b/src/hydroturing/spec.py @@ -177,6 +177,9 @@ "reference_driven_exchange", "reference_evolving_exchange", "reference_recharge_exchange", + "reference_recharge_order_exact", + "reference_recharge_blind", + "reference_recharge_overshoot", "reference_noise_sink", "reference_token_exchange", "reference_steady_sink", diff --git a/tests/test_report_detail.py b/tests/test_report_detail.py index 54138a3e..88794525 100644 --- a/tests/test_report_detail.py +++ b/tests/test_report_detail.py @@ -47,6 +47,7 @@ "mass/extreme-rain": ("monotone_response",), "mass/exchange-response": ("exchange_response",), "mass/gw-sw-exchange-consistency": ("exchange_components",), + "mass/aquifer-recharge-ordering": ("aquifer_recharge_ordering",), "mass/human-abstraction": ("human_abstraction",), "mass/multi-decadal-drift": ("state_bounds", "total_storage_drift"), "mass/phase-counterfactual": ("phase_invariance",), From 6abeb976574b634363bccba171ae27870bbae26c Mon Sep 17 00:00:00 2001 From: Zhenjiang Wu <852356187@qq.com> Date: Tue, 29 Sep 2026 11:16:42 +0800 Subject: [PATCH 2/5] Address review feedback for recharge ordering probe --- README.md | 5 +- docs/writing-a-probe.md | 3 + models/modflow6/ht_adapter.py | 10 +-- models/modflow6/model.yaml | 6 +- models/reference_recharge_blind/ht_adapter.py | 3 +- models/reference_recharge_blind/model.yaml | 2 +- .../ht_adapter.py | 3 +- .../reference_recharge_order_exact/model.yaml | 2 +- .../ht_adapter.py | 3 +- .../reference_recharge_overshoot/model.yaml | 2 +- models/result.csv | 54 +++++++++++++--- .../mass/aquifer-recharge-ordering/README.md | 2 + .../mass/aquifer-recharge-ordering/probe.yaml | 6 +- site/index.html | 13 ++-- src/hydroturing/criteria/recharge_ordering.py | 62 ++++++++++++------- tests/test_recharge_ordering.py | 26 ++++++++ 16 files changed, 149 insertions(+), 53 deletions(-) create mode 100644 tests/test_recharge_ordering.py diff --git a/README.md b/README.md index 58b84501..f639c657 100644 --- a/README.md +++ b/README.md @@ -168,9 +168,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**, 1 of 1 probe scored | +| [`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. The updated datum/recharge paired cases are currently N/A because the local Docker runner is unavailable. | **N/A (INCOMPLETE)**, 0 of 0 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 groundwater-exchange probe 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 | passive aquifer reference that partitions direct recharge between storage and river exchange | must pass `mass/aquifer-recharge-ordering` | +| `reference_recharge_blind` | broken | ignores the added recharge while keeping an internally closed budget | caught by `aquifer_recharge_ordering` | +| `reference_recharge_overshoot` | broken | explicit update overshoots after the recharge pulse | caught by `aquifer_recharge_ordering` | | [`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**, 24 of 24 | | [`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**, 23 of 23 | | [`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**, 24 of 24 | diff --git a/docs/writing-a-probe.md b/docs/writing-a-probe.md index cf690376..e54c8c64 100644 --- a/docs/writing-a-probe.md +++ b/docs/writing-a-probe.md @@ -372,6 +372,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` | passive aquifer reference that partitions recharge exactly | nothing, it must pass | +| `reference_recharge_blind` | ignores the added recharge while closing its own budget | `aquifer_recharge_ordering` | +| `reference_recharge_overshoot` | explicit update overshoots after the pulse | `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/modflow6/ht_adapter.py b/models/modflow6/ht_adapter.py index 5c6965b2..6ed1aead 100644 --- a/models/modflow6/ht_adapter.py +++ b/models/modflow6/ht_adapter.py @@ -14,12 +14,10 @@ import flopy -MODEL = {"name": "modflow6", "version": "6.7.0-adapter.1"} +MODEL = {"name": "modflow6", "version": "6.7.0-adapter.2"} COLUMNS = ["time", "gw_sw_exchange", "gw_to_sw", "sw_to_gw", "gw"] NROW = 10 NCOL = 10 -TOP = 20.0 -BOTM = 0.0 K = 1.0 RIVER_ROW = 5 RIVER_COL = 5 @@ -86,6 +84,8 @@ def run_modflow(forcing: list[dict], static: dict) -> list[dict]: river_conductance = float(static["river_conductance_m2_per_day"]) river_bottom_offset = float(static["river_bottom_offset_m"]) initial_head = float(static.get("aquifer_initial_head_m", 10.0)) + top = float(static["aquifer_top_m"]) + botm = float(static["aquifer_bottom_m"]) with tempfile.TemporaryDirectory(prefix="hydroturing-mf6-") as workdir: sim = flopy.mf6.MFSimulation(sim_name="ht_gw", sim_ws=workdir, exe_name=mf6_exe) @@ -101,7 +101,7 @@ def run_modflow(forcing: list[dict], static: dict) -> list[dict]: gwf = flopy.mf6.ModflowGwf(sim, modelname="ht_gw", save_flows=True) flopy.mf6.ModflowGwfdis( gwf, nlay=1, nrow=NROW, ncol=NCOL, delr=delr, delc=delc, - top=TOP, botm=BOTM, + top=top, botm=botm, ) flopy.mf6.ModflowGwfic(gwf, strt=initial_head) flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=K) @@ -135,7 +135,7 @@ def run_modflow(forcing: list[dict], static: dict) -> list[dict]: budget_file = gwf.output.budget() rows = [] - gw_mm = specific_yield * max(initial_head - BOTM, 0.0) * 1000.0 + gw_mm = specific_yield * max(initial_head - botm, 0.0) * 1000.0 for period, forcing_row in enumerate(forcing): q_records = budget_file.get_data(text="RIV", kstpkper=(0, period))[0] # MODFLOW RIV q is positive into the aquifer, i.e. the river diff --git a/models/modflow6/model.yaml b/models/modflow6/model.yaml index 68e01b24..b96c613a 100644 --- a/models/modflow6/model.yaml +++ b/models/modflow6/model.yaml @@ -1,5 +1,5 @@ name: modflow6 -version: "6.7.0-adapter.1" +version: "6.7.0-adapter.2" description: > MODFLOW 6.7.0, the USGS's modular hydrologic model, run as a transient groundwater model with one river cell (MODFLOW-ORG/modflow6, release @@ -18,8 +18,8 @@ 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_storage_coefficient, aquifer_specific_yield, river_conductance_m2_per_day, river_bottom_offset_m] +needs_static: [area_km2, aquifer_storage_coefficient, aquifer_specific_yield, aquifer_initial_head_m, aquifer_top_m, aquifer_bottom_m, river_conductance_m2_per_day, river_bottom_offset_m] supports: - perturbation: false + perturbation: true resources: {cpu: 1, memory_gb: 2, gpu: false} window_days: full diff --git a/models/reference_recharge_blind/ht_adapter.py b/models/reference_recharge_blind/ht_adapter.py index 119fc79d..8f8d226b 100644 --- a/models/reference_recharge_blind/ht_adapter.py +++ b/models/reference_recharge_blind/ht_adapter.py @@ -1,6 +1,7 @@ from __future__ import annotations import argparse,csv,json from pathlib import Path +MODE = "blind" COLUMNS=["time","gw_sw_exchange","gw_to_sw","sw_to_gw","gw"] def simulate(forcing,static,mode="exact"): sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] @@ -14,7 +15,7 @@ def simulate(forcing,static,mode="exact"): def main(): p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) - static=json.loads((d/req["input"]["static"]).read_text()); name=Path(__file__).parent.name; mode="blind" if "blind" in name else ("overshoot" if "overshoot" in name else "exact"); rows=simulate(forcing,static,mode) + static=json.loads((d/req["input"]["static"]).read_text()); rows=simulate(forcing,static,MODE) with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) (d/req["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 index fc00e11d..74c4dce2 100644 --- a/models/reference_recharge_blind/model.yaml +++ b/models/reference_recharge_blind/model.yaml @@ -10,7 +10,7 @@ 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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] +needs_static: [area_km2, aquifer_specific_yield, aquifer_initial_head_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 index 119fc79d..042b46fe 100644 --- a/models/reference_recharge_order_exact/ht_adapter.py +++ b/models/reference_recharge_order_exact/ht_adapter.py @@ -1,6 +1,7 @@ from __future__ import annotations import argparse,csv,json from pathlib import Path +MODE = "exact" COLUMNS=["time","gw_sw_exchange","gw_to_sw","sw_to_gw","gw"] def simulate(forcing,static,mode="exact"): sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] @@ -14,7 +15,7 @@ def simulate(forcing,static,mode="exact"): def main(): p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) - static=json.loads((d/req["input"]["static"]).read_text()); name=Path(__file__).parent.name; mode="blind" if "blind" in name else ("overshoot" if "overshoot" in name else "exact"); rows=simulate(forcing,static,mode) + static=json.loads((d/req["input"]["static"]).read_text()); rows=simulate(forcing,static,MODE) with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) (d/req["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 index 8253ca90..0455544e 100644 --- a/models/reference_recharge_order_exact/model.yaml +++ b/models/reference_recharge_order_exact/model.yaml @@ -10,7 +10,7 @@ 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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] +needs_static: [area_km2, aquifer_specific_yield, aquifer_initial_head_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 index 119fc79d..5c0fe3c2 100644 --- a/models/reference_recharge_overshoot/ht_adapter.py +++ b/models/reference_recharge_overshoot/ht_adapter.py @@ -1,6 +1,7 @@ from __future__ import annotations import argparse,csv,json from pathlib import Path +MODE = "overshoot" COLUMNS=["time","gw_sw_exchange","gw_to_sw","sw_to_gw","gw"] def simulate(forcing,static,mode="exact"): sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] @@ -14,7 +15,7 @@ def simulate(forcing,static,mode="exact"): def main(): p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) - static=json.loads((d/req["input"]["static"]).read_text()); name=Path(__file__).parent.name; mode="blind" if "blind" in name else ("overshoot" if "overshoot" in name else "exact"); rows=simulate(forcing,static,mode) + static=json.loads((d/req["input"]["static"]).read_text()); rows=simulate(forcing,static,MODE) with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) (d/req["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 index d15c2ae5..81e2f27b 100644 --- a/models/reference_recharge_overshoot/model.yaml +++ b/models/reference_recharge_overshoot/model.yaml @@ -10,7 +10,7 @@ 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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] +needs_static: [area_km2, aquifer_specific_yield, aquifer_initial_head_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 1e04d306..62457ea2 100644 --- a/models/result.csv +++ b/models/result.csv @@ -726,13 +726,47 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-09-23,flex_lumped,1.2.0,0.1.0,subprocess,mass/spinup-cycle-invariance,PASS,OK,full record,107008012 613960126 1120912240,spinup_cycle_invariance: the repeated forcing reaches the same evaluation-year response for all 3 selected spin-up cycles (worst short<->long:gw departure 0.00%) 2026-09-23,flex_topo,1.1.0,0.1.0,subprocess,mass/spinup-cycle-invariance,PASS,OK,full record,107008012 613960126 1120912240,spinup_cycle_invariance: the repeated forcing reaches the same evaluation-year response for all 3 selected spin-up cycles (worst short<->plus3:channel departure 0.00%) 2026-09-23,sacsma_snow17,1.1.0,0.1.0,subprocess,mass/spinup-cycle-invariance,PASS,OK,full record,107008012 613960126 1120912240,spinup_cycle_invariance: the repeated forcing reaches the same evaluation-year response for all 3 selected spin-up cycles (worst short<->long:gw departure 0.00%) -2026-09-28,cwatm,1.11-5baaadd.4,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,dhbv2,0.5.4-hbv2ep100.3,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,flex_lumped,1.2.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,flex_topo,1.1.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,google_flood_forecast,0.1.0-828dfc5,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,lisflood,5.0.0-onecell.5,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,modflow6,6.7.0-adapter.1,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,sacsma_snow17,1.1.0,0.1.0,subprocess,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,summa,4.0.0-f787fa5.4,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases -2026-09-28,wflow_sbm,1.0.4-ht.4,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPLETE,not run,,does not declare support for paired groundwater recharge cases +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPATIBLE,not run,500220947 1007173061 1514125175,"static does not provide aquifer_top_m, aquifer_bottom_m" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,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, river_conductance_m2_per_day" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/evaporative-partition,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, hfls, hfss, hfg, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/latent-heat-et-consistency,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, hfls, hfss, hfg, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/pet-consistency,N/A,INCOMPLETE,not run,,"does not report evspsbl, mrro, mrso" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/radiation-consistency,N/A,INCOMPLETE,not run,,"does not report rlus, ts; model timestep PT1D does not cover the probe's PT1H; model does not declare that it consumes forcing rlds; model does not declare that it consumes static eps" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/snowmelt-energy-water,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, sbl, hfls, hfss, hfg, mrso, snw, canopy, lwsnl, csnow; model does not declare that it consumes forcing pr, rn" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/soil-heat-storage-consistency,N/A,INCOMPLETE,not run,,"does not report hfg, hfg_bottom, tsoil_layer; model timestep PT1D does not cover the probe's PT1H; model does not declare that it consumes forcing pr, tas, rsds, rlds; model does not declare that it consumes static soil_heat_capacity_areal, soil_layer_depth_m, soil_temperature_initial" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/surface-energy-closure,N/A,INCOMPLETE,not run,,"does not report hfls, hfss, hfg; model timestep PT1D does not cover the probe's PT1H" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/antecedent-monotonicity,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPATIBLE,not run,500220947 1007173061 1514125175,"static does not provide aquifer_top_m, aquifer_bottom_m" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/area-invariance,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/catchment-closure,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/causality,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/dry-down,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/exchange-response,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, gwex, mrso, snw, canopy; model does not declare that it consumes forcing gwh" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/extreme-event-closure,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/extreme-rain,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/gw-sw-exchange-consistency,N/A,INCOMPATIBLE,not run,479803588 986755702 1493707816 2000659930 360128397,"static does not provide aquifer_top_m, aquifer_bottom_m" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/human-abstraction,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/multi-decadal-drift,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/phase-counterfactual,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/precipitation-counterfactual,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/resolution-invariance,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/response-nonnegativity,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/runoff-bounds,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/snowpack-mass-closure,N/A,INCOMPLETE,not run,,"does not report snm, snw" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/spinup-cycle-invariance,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/steady-state,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/time-origin-invariance,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/ungauged-basin-closure,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/warming-response,N/A,INCOMPLETE,not run,,does not report mrro +2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,momentum/routing-conservation,N/A,INCOMPLETE,not run,,"does not report mrro, channel" +2026-09-29,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-09-29,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-09-29,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" diff --git a/probes/mass/aquifer-recharge-ordering/README.md b/probes/mass/aquifer-recharge-ordering/README.md index 1f83444c..0930593b 100644 --- a/probes/mass/aquifer-recharge-ordering/README.md +++ b/probes/mass/aquifer-recharge-ordering/README.md @@ -5,3 +5,5 @@ The paired case adds only direct `gw_recharge` to the same passive aquifer. The the incremental exchange, and the cumulative partition identity using the added recharge volume as its scale. A plateau with zero incremental exchange is valid; a recharge-blind or nonmonotone update is not. + +The exact reference supplies the deterministic gate. The independent MODFLOW 6 baseline is enabled for paired perturbations and is archived separately; the four standard catchment references do not expose a direct aquifer recharge interface. diff --git a/probes/mass/aquifer-recharge-ordering/probe.yaml b/probes/mass/aquifer-recharge-ordering/probe.yaml index ac0aaa7d..889e5bd4 100644 --- a/probes/mass/aquifer-recharge-ordering/probe.yaml +++ b/probes/mass/aquifer-recharge-ordering/probe.yaml @@ -17,7 +17,7 @@ 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_storage_coefficient, river_conductance_m2_per_day, river_bottom_offset_m] + static: [area_km2, aquifer_specific_yield, aquifer_initial_head_m, river_conductance_m2_per_day] case: generator: generate.py n_seeds: 3 @@ -37,6 +37,7 @@ criteria: storage: gw rel_tol: 0.01 abs_tol_mm: 1.0e-3 + sources: [gw_boundary] - groundwater_balance: recharge: gw_recharge exchange: gw_sw_exchange @@ -52,6 +53,9 @@ criteria: 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 diff --git a/site/index.html b/site/index.html index baafc2af..42b602da 100644 --- a/site/index.html +++ b/site/index.html @@ -268,7 +268,7 @@

How it works

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

How it works

gw exchange exchange + recharge order ENERGY @@ -722,7 +723,7 @@

Credit

summasubmittedSUMMA 4.0.0, a land model solving the water and energy balances of canopy, snow, soil and aquifer with one implicit solver; one lumped HRU, with the radiation and humidity it needs mocked from the forcing.9 / 24FAIL (VIOLATION) cwatmsubmittedCWatM 1.11, IIASA's Community Water Model: a process-based global hydrological model run on one grid cell, reporting its snow, soil, canopy, groundwater and runoff-concentration stores.17 / 21FAIL (VIOLATION) lisfloodsubmittedLISFLOOD 5.0.0, the EC JRC's distributed model behind EFAS and GloFAS, stepped through its own Python framework on one representative 5 km cell.18 / 21FAIL (ERROR) -modflow6submittedMODFLOW 6.7.0, the USGS's modular groundwater model: a transient one-river-cell aquifer run through flopy, reporting signed exchange and groundwater storage from its own RIV and STO budget terms.1 / 1PASS +modflow6submittedMODFLOW 6.7.0, the USGS's modular groundwater model with datum-aware paired inputs.0 / 0N/A (INCOMPLETE) flex_lumpedphysicalLumped FLEX/HBV from chrimerss/HydrologicModels: interception, beta-partitioned soil, fast and slow reservoirs, a geometry-aware lag.24 / 24PASS flex_topophysicalFLEX-Topo: plateau, hillslope and wetland units on the Wark catchment's real fractions, sharing one groundwater store.23 / 23PASS sacsma_snow17physicalThe NWS operational SAC-SMA with Snow-17 and a gamma unit hydrograph, ported from the legacy Fortran and checked against it.24 / 24PASS @@ -744,7 +745,7 @@

Credit

mass/extreme-rainmassThe largest storm scaled to ten times: runoff cannot fall, nor exceed the rain added.merged mass/runoff-boundsmassOver ten years, is the runoff possible at all: at least rain minus demand minus storage, at most rain plus storage? Runoff only.merged mass/area-invariancemassThe same weather on a ten times larger catchment: every depth identical.merged -mass/aquifer-recharge-orderingmassThe same weather on a ten times larger catchment: every depth identical.merged +mass/aquifer-recharge-orderingmassAdded direct aquifer recharge preserves storage and incremental exchange ordering.merged mass/response-nonnegativitymassOne storm added: from that day on, runoff may never be lower than without it.merged mass/antecedent-monotonicitymassThe same storm after a dry month and a wet one: more runoff from the wet one, no more than the extra water.merged mass/phase-counterfactualmassThe same water as rain instead of snow: timing moves, the integrated volumes do not.merged @@ -947,7 +948,7 @@

Credit

summaenviadoSUMMA 4.0.0, un modelo de superficie terrestre que resuelve los balances de agua y energía del dosel, la nieve, el suelo y el acuífero con un único esquema implícito; una HRU agregada, con la radiación y la humedad que necesita simuladas a partir del forzamiento.9 / 24FALLA (VIOLACIÓN) cwatmenviadoCWatM 1.11, el Community Water Model del IIASA: un modelo hidrológico global basado en procesos, ejecutado sobre una sola celda, que reporta sus almacenes de nieve, suelo, dosel, agua subterránea y concentración de escorrentía.17 / 21FALLA (VIOLACIÓN) lisfloodenviadoLISFLOOD 5.0.0, el modelo distribuido del JRC de la Comisión Europea detrás de EFAS y GloFAS, avanzado paso a paso con su propio marco de Python en una celda representativa de 5 km.18 / 21FALLA (ERROR) -modflow6enviadoMODFLOW 6.7.0, el modelo modular de agua subterránea del USGS: un acuífero transitorio de una celda de río ejecutado con flopy, que reporta el intercambio señalado y el almacenamiento subterráneo a partir de sus propios términos de balance RIV y STO.1 / 1PASA +modflow6enviadoMODFLOW 6.7.0, el modelo modular de agua subterránea del USGS: un acuífero transitorio de una celda de río ejecutado con flopy, que reporta el intercambio señalado y el almacenamiento subterráneo a partir de sus propios términos de balance RIV y STO.0 / 0N/A (INCOMPLETE) flex_lumpedfísicoFLEX/HBV agregado de chrimerss/HydrologicModels: intercepción, suelo con partición beta, embalses rápido y lento, un retardo basado en la geometría.24 / 24PASA flex_topofísicoFLEX-Topo: unidades de meseta, ladera y humedal con las fracciones reales de la cuenca Wark, un solo almacén subterráneo.23 / 23PASA sacsma_snow17físicoEl SAC-SMA operativo del NWS con Snow-17 y un hidrograma unitario gamma, portado del Fortran heredado y verificado contra él.24 / 24PASA @@ -969,7 +970,7 @@

Credit

mass/extreme-rainmasaLa mayor tormenta multiplicada por diez: la escorrentía no puede bajar ni superar la lluvia añadida.integrada mass/runoff-boundsmasaEn diez años, ¿es posible esa escorrentía? Al menos lluvia menos demanda menos almacenamiento; a lo sumo lluvia más almacenamiento. Solo con caudal.integrada mass/area-invariancemasaEl mismo tiempo sobre una cuenca diez veces mayor: cada lámina idéntica.integrada -mass/aquifer-recharge-orderingmasaEl mismo tiempo sobre una cuenca diez veces mayor: cada lámina idéntica.integrada +mass/aquifer-recharge-orderingmasaLa recarga directa conserva el almacenamiento y el orden del intercambio incremental.integrada mass/response-nonnegativitymasaUna tormenta añadida: desde ese día la escorrentía nunca puede ser menor que sin ella.integrada mass/antecedent-monotonicitymasaLa misma tormenta tras un mes seco y tras uno húmedo: más escorrentía en el húmedo, y no más que el agua extra.integrada mass/phase-counterfactualmasaLa misma agua como lluvia en vez de nieve: cambia el momento, no los volúmenes integrados.integrada @@ -1172,7 +1173,7 @@

Credit

summa提交模型SUMMA 4.0.0，陆面模型：以同一隐式求解器求解冠层、积雪、土壤与含水层的水量和能量平衡；按单个集总 HRU 运行，所需辐射与湿度由强迫数据模拟。9 / 24未通过（违反） cwatm提交模型CWatM 1.11，IIASA 的社区水模型（Community Water Model）：基于过程的全球水文模型，在单个网格上运行，报告积雪、土壤、冠层、地下水与产流汇集滞蓄。17 / 21未通过（违反） lisflood提交模型LISFLOOD 5.0.0，欧盟委员会联合研究中心支撑 EFAS 与 GloFAS 的分布式模型，通过其自身的 Python 框架在单个 5 公里代表性网格上逐步运行。18 / 21未通过（错误） -modflow6提交模型MODFLOW 6.7.0，美国地质调查局的模块化地下水模型：通过 flopy 运行的单河流网格瞬态含水层，依据其自身的 RIV 与 STO 收支项报告带符号交换量与地下水储量。1 / 1通过 +modflow6提交模型MODFLOW 6.7.0，美国地质调查局的模块化地下水模型：通过 flopy 运行的单河流网格瞬态含水层，依据其自身的 RIV 与 STO 收支项报告带符号交换量与地下水储量。0 / 0N/A (INCOMPLETE) flex_lumped物理模型来自 chrimerss/HydrologicModels 的集总 FLEX/HBV：截留、beta 分配的土壤、快慢水库与基于几何的滞时。24 / 24通过 flex_topo物理模型FLEX-Topo：台地、坡地与湿地单元按 Wark 流域的真实面积比共享一个地下水库。23 / 23通过 sacsma_snow17物理模型美国国家气象局业务用 SAC-SMA + Snow-17 + 伽马单位线，由遗留 Fortran 移植并与之比对。24 / 24通过 diff --git a/src/hydroturing/criteria/recharge_ordering.py b/src/hydroturing/criteria/recharge_ordering.py index d37d3a66..85fcfaa6 100644 --- a/src/hydroturing/criteria/recharge_ordering.py +++ b/src/hydroturing/criteria/recharge_ordering.py @@ -1,25 +1,45 @@ +"""Order preservation for a nonnegative direct aquifer-recharge pulse.""" from __future__ import annotations import numpy as np -from hydroturing.criteria.base import FAIL,PASS,CriterionResult,criterion,make_window +from hydroturing.criteria.base import FAIL, PASS, CriterionResult, criterion, make_window from hydroturing.protocol import RunResult from hydroturing.spec import ProbeSpec -@criterion("aquifer_recharge_ordering",paired=True) -def aquifer_recharge_ordering(runs:dict[str,RunResult],probe:ProbeSpec,params:dict)->CriterionResult: - c=str(params.get("control","control")); p=str(params.get("perturbed","recharge_added")) - if c not in runs or p not in runs: raise ValueError("aquifer_recharge_ordering needs control and perturbed variants") - wc=make_window(runs[c],probe); wp=make_window(runs[p],probe) - recharge=str(params.get("recharge","gw_recharge")); exchange=str(params.get("exchange","gw_sw_exchange")); storage=str(params.get("storage","gw")) - for w,n in ((wc,"control"),(wp,"perturbed")): - if recharge not in w.forcing.columns: raise ValueError(f"{n} forcing needs '{recharge}'") - for v in (exchange,storage): - if v not in w.table.columns: raise ValueError(f"{n} result needs '{v}'") - if len(wc.table)!=len(wp.table): raise ValueError("paired recharge runs have different scored lengths") - added=wp.volume(wp.forcing[recharge])-wc.volume(wc.forcing[recharge]); I=np.cumsum(added); Istar=float(I[-1]) - if Istar<=0: raise ValueError("aquifer_recharge_ordering needs positive added recharge") - ds=wp.table[storage].to_numpy(float)-wc.table[storage].to_numpy(float) - ds-=float(wp.state0[storage])-float(wc.state0[storage]) - d=wc.volume(wc.table[exchange])-wp.volume(wp.table[exchange]); D=np.cumsum(d) - es=float(np.max(np.maximum(-ds,0.0))); eq=float(np.sum(np.maximum(-d,0.0))); eb=float(np.max(np.abs(ds+D-I))) - eps=max(float(params.get("rel_tol",0.01))*Istar,float(params.get("abs_tol_mm",1e-3))) - ratios={"storage":es/eps,"exchange":eq/eps,"budget":eb/eps}; worst_name,worst=max(ratios.items(),key=lambda x:x[1]); ok=worst<=1.0 - return CriterionResult(name="aquifer_recharge_ordering",status=PASS if ok else FAIL,value=worst,threshold=1.0,message=(f"added recharge preserves storage and exchange ordering (I={Istar:.6g} mm)" if ok else f"recharge ordering fails {worst_name}: {max(es,eq,eb):.6g} mm exceeds {eps:.6g} mm"),diagnostics={"added_recharge_mm":Istar,"epsilon_mm":eps,"E_S_mm":es,"E_Q_mm":eq,"E_B_mm":eb,"ratios":ratios}) + +@criterion("aquifer_recharge_ordering", paired=True) +def aquifer_recharge_ordering(runs: dict[str, RunResult], probe: ProbeSpec, params: dict) -> CriterionResult: + """Check storage ordering, exchange ordering and cumulative partitioning. + + E_S is a maximum because one negative excursion violates comparison; E_Q + is a sum because any reversed increment is water sent in the wrong + direction; E_B is a maximum so local losses cannot cancel later. The + storage floor is capped to protect absolute-storage rounding without + letting a large datum buy unlimited tolerance. + """ + control_name = str(params.get("control", "control")); perturbed_name = str(params.get("perturbed", "recharge_added")) + if control_name not in runs or perturbed_name not in runs: raise ValueError("aquifer_recharge_ordering needs control and perturbed variants") + control = make_window(runs[control_name], probe); perturbed = make_window(runs[perturbed_name], probe) + recharge = str(params.get("recharge", "gw_recharge")); exchange = str(params.get("exchange", "gw_sw_exchange")); storage = str(params.get("storage", "gw")) + for window, label in ((control, "control"), (perturbed, "perturbed")): + if recharge not in window.forcing.columns: raise ValueError(f"{label} forcing needs '{recharge}'") + for name in (exchange, storage): + if name not in window.table.columns: raise ValueError(f"{label} result needs '{name}'") + if len(control.table) != len(perturbed.table): raise ValueError("paired recharge runs have different scored lengths") + added = perturbed.volume(perturbed.forcing[recharge]) - control.volume(control.forcing[recharge]) + tol = float(params.get("input_tol_mm", 1.0e-10)) + if np.any(added < -tol): raise ValueError("aquifer_recharge_ordering requires nonnegative added recharge at every step") + cumulative_added = np.cumsum(added); added_total = float(cumulative_added[-1]) + if added_total <= tol: raise ValueError("aquifer_recharge_ordering needs positive added recharge") + delta_storage = perturbed.table[storage].to_numpy(float) - control.table[storage].to_numpy(float) + delta_storage -= float(perturbed.state0[storage]) - float(control.state0[storage]) + d_exchange = control.volume(control.table[exchange]) - perturbed.volume(perturbed.table[exchange]); cumulative_exchange = np.cumsum(d_exchange) + source_delta = np.zeros(len(control.table)) + for source in params.get("sources", []): + if source in control.table.columns and source in perturbed.table.columns: + source_delta += perturbed.volume(perturbed.table[source]) - control.volume(control.table[source]) + storage_error = float(np.max(np.maximum(-delta_storage, 0.0))); exchange_error = float(np.sum(np.maximum(-d_exchange, 0.0))) + budget_error = float(np.max(np.abs(delta_storage + cumulative_exchange + np.cumsum(source_delta) - cumulative_added))) + relative = float(params.get("rel_tol", 0.01)) * added_total; absolute = float(params.get("abs_tol_mm", 1.0e-3)) + storage_max = max(float(np.max(np.abs(control.table[storage]))), float(np.max(np.abs(perturbed.table[storage])))); storage_floor = min(1.0e-5 * storage_max, 0.05) + epsilon = max(relative, absolute, storage_floor) + ratios = {"storage": storage_error / epsilon, "exchange": exchange_error / epsilon, "budget": budget_error / epsilon}; worst_name, worst = max(ratios.items(), key=lambda item: item[1]); ok = worst <= 1.0 + return CriterionResult(name="aquifer_recharge_ordering", status=PASS if ok else FAIL, value=worst, threshold=1.0, message=(f"added recharge preserves storage and exchange ordering (I={added_total:.6g} mm)" if ok else f"recharge ordering fails {worst_name}: {max(storage_error, exchange_error, budget_error):.6g} mm exceeds {epsilon:.6g} mm"), diagnostics={"added_recharge_mm": added_total, "epsilon_mm": epsilon, "E_S_mm": storage_error, "E_Q_mm": exchange_error, "E_B_mm": budget_error, "ratios": ratios}) diff --git a/tests/test_recharge_ordering.py b/tests/test_recharge_ordering.py new file mode 100644 index 00000000..2e5f28c5 --- /dev/null +++ b/tests/test_recharge_ordering.py @@ -0,0 +1,26 @@ +from __future__ import annotations +import numpy as np, pandas as pd, pytest +from dataclasses import replace +from hydroturing import registry +from hydroturing.criteria.recharge_ordering import aquifer_recharge_ordering +from hydroturing.harness import build_case +from hydroturing.protocol import RunResult + +def make(qc, qp, sc, sp, rc=None, rp=None): + probe=registry.find_probe('mass/aquifer-recharge-ordering'); case=build_case(probe,7,'control'); control_case=case; n=len(case.forcing) + def run(q,s,r,run_case): + forcing=run_case.forcing.copy(); forcing['gw_recharge']=np.asarray(r if r is not None else np.zeros(n)); table=pd.DataFrame({'gw_sw_exchange':q,'gw':s,'gw_to_sw':np.minimum(q,0),'sw_to_gw':np.maximum(q,0)}); return RunResult(case=replace(run_case, forcing=forcing),table=table,meta={},wall_seconds=0.) + control=run(np.asarray(qc),np.asarray(sc),rc if rc is not None else np.zeros(n),control_case) + # Rebuild the case so the two forcing frames do not alias. + case=build_case(probe,7,'recharge_added') + perturbed=run(np.asarray(qp),np.asarray(sp),rp if rp is not None else np.zeros(n),case) + return {'control':control,'recharge_added':perturbed},probe + +def params(probe): return dict(next(c.params for c in probe.criteria if c.name=='aquifer_recharge_ordering')) +def test_recharge_ordering_rejects_exchange_reversal(): + n=len(build_case(registry.find_probe('mass/aquifer-recharge-ordering'),7,'control').forcing); r=np.zeros(n); r[2:4]=5; q=np.zeros(n); qp=q.copy(); qp[2]=-1; s=np.full(n,2000.); sp=s.copy(); sp[2:]=sp[2:]+np.cumsum(r)[2:] + runs,p=make(q,qp,s,sp,np.zeros(n),r); assert not aquifer_recharge_ordering(runs,p,params(p)).passed +def test_recharge_ordering_rejects_negative_input_increment(): + n=len(build_case(registry.find_probe('mass/aquifer-recharge-ordering'),7,'control').forcing); r=np.zeros(n); r[2]=5; r[3]=-1; q=np.zeros(n); s=np.full(n,2000.); runs,p=make(q,q,s,s,r,np.zeros(n)); + with pytest.raises(ValueError,match='nonnegative'): aquifer_recharge_ordering(runs,p,params(p)) + From 197e214d978687a2ebba8bbb6a215999c1b49ccb Mon Sep 17 00:00:00 2001 From: Zhenjiang Wu <852356187@qq.com> Date: Tue, 29 Sep 2026 11:24:01 +0800 Subject: [PATCH 3/5] Document recharge ordering tolerance derivation --- .../mass/aquifer-recharge-ordering/README.md | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/probes/mass/aquifer-recharge-ordering/README.md b/probes/mass/aquifer-recharge-ordering/README.md index 0930593b..45180892 100644 --- a/probes/mass/aquifer-recharge-ordering/README.md +++ b/probes/mass/aquifer-recharge-ordering/README.md @@ -7,3 +7,22 @@ added recharge volume as its scale. A plateau with zero incremental exchange is valid; a recharge-blind or nonmonotone update is not. The exact reference supplies the deterministic gate. The independent MODFLOW 6 baseline is enabled for paired perturbations and is archived separately; the four standard catchment references do not expose a direct aquifer recharge interface. + +## Criterion and tolerance + +For each scored step, let `I_n` be cumulative added recharge and `D_n` be +cumulative incremental exchange (positive when the added recharge sends more +water to the river). The criterion records three one-sided errors: + +* `E_S = max([-ΔS_n]_+)`, the largest storage decrease; +* `E_Q = sum([-d_i]_+ · Δt_i)`, the total exchange increment in the wrong + direction; and +* `E_B = max(|ΔS_n + D_n + B_n - I_n|)`, the worst cumulative partition + residual, where `B_n` is the cumulative declared boundary-source change. + +The common tolerance is `ε = max(0.01 · I_* , 1e-3 mm, f · S_max)` where +`I_*` is total added recharge, `S_max` is the largest absolute storage in +either run, and `f · S_max` is a capped storage precision floor (`f = 1e-5`, +cap `0.05 mm`). Every error must be no larger than `ε`. This makes the scale +come from the perturbation itself while accounting for storage rounding without +allowing a large absolute datum to buy unlimited tolerance. From eda8355206dda171f6383882cdb0c0f27033409a Mon Sep 17 00:00:00 2001 From: Zhenjiang Wu <852356187@qq.com> Date: Tue, 29 Sep 2026 15:54:34 +0800 Subject: [PATCH 4/5] Run recharge ordering probe against MODFLOW 6 --- README.md | 2 +- models/reference_recharge_blind/model.yaml | 2 +- models/reference_recharge_order_exact/model.yaml | 2 +- models/reference_recharge_overshoot/model.yaml | 2 +- models/result.csv | 3 +-- probes/mass/aquifer-recharge-ordering/generate.py | 2 +- probes/mass/aquifer-recharge-ordering/probe.yaml | 2 +- site/index.html | 6 +++--- 8 files changed, 10 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index f639c657..f957235b 100644 --- a/README.md +++ b/README.md @@ -168,7 +168,7 @@ 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. The updated datum/recharge paired cases are currently N/A because the local Docker runner is unavailable. | **N/A (INCOMPLETE)**, 0 of 0 probes scored | +| [`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. The updated recharge paired case is covered by the local MODFLOW 6 container. | **PASS**, 1 of 1 probe 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 groundwater-exchange probe 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 | passive aquifer reference that partitions direct recharge between storage and river exchange | must pass `mass/aquifer-recharge-ordering` | diff --git a/models/reference_recharge_blind/model.yaml b/models/reference_recharge_blind/model.yaml index 74c4dce2..07818d42 100644 --- a/models/reference_recharge_blind/model.yaml +++ b/models/reference_recharge_blind/model.yaml @@ -10,7 +10,7 @@ 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, river_conductance_m2_per_day] +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/model.yaml b/models/reference_recharge_order_exact/model.yaml index 0455544e..b7d01ec0 100644 --- a/models/reference_recharge_order_exact/model.yaml +++ b/models/reference_recharge_order_exact/model.yaml @@ -10,7 +10,7 @@ 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, river_conductance_m2_per_day] +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/model.yaml b/models/reference_recharge_overshoot/model.yaml index 81e2f27b..e719313b 100644 --- a/models/reference_recharge_overshoot/model.yaml +++ b/models/reference_recharge_overshoot/model.yaml @@ -10,7 +10,7 @@ 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, river_conductance_m2_per_day] +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 62457ea2..5c08e250 100644 --- a/models/result.csv +++ b/models/result.csv @@ -732,7 +732,6 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-09-29,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, river_conductance_m2_per_day" 2026-09-29,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, river_conductance_m2_per_day" 2026-09-29,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, river_conductance_m2_per_day" -2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPATIBLE,not run,500220947 1007173061 1514125175,"static does not provide aquifer_top_m, aquifer_bottom_m" 2026-09-29,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, river_conductance_m2_per_day" 2026-09-29,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, river_conductance_m2_per_day" 2026-09-29,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, river_conductance_m2_per_day" @@ -744,7 +743,6 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/soil-heat-storage-consistency,N/A,INCOMPLETE,not run,,"does not report hfg, hfg_bottom, tsoil_layer; model timestep PT1D does not cover the probe's PT1H; model does not declare that it consumes forcing pr, tas, rsds, rlds; model does not declare that it consumes static soil_heat_capacity_areal, soil_layer_depth_m, soil_temperature_initial" 2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,energy/surface-energy-closure,N/A,INCOMPLETE,not run,,"does not report hfls, hfss, hfg; model timestep PT1D does not cover the probe's PT1H" 2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/antecedent-monotonicity,N/A,INCOMPLETE,not run,,does not report mrro -2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/aquifer-recharge-ordering,N/A,INCOMPATIBLE,not run,500220947 1007173061 1514125175,"static does not provide aquifer_top_m, aquifer_bottom_m" 2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/area-invariance,N/A,INCOMPLETE,not run,,does not report mrro 2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/catchment-closure,N/A,INCOMPLETE,not run,,"does not report pr, evspsbl, mrro, mrso, snw, canopy" 2026-09-29,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/causality,N/A,INCOMPLETE,not run,,does not report mrro @@ -770,3 +768,4 @@ run_date,model,version,suite_version,runner,probe,verdict,reason,window,seeds,de 2026-09-29,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-09-29,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-09-29,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-09-29,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) diff --git a/probes/mass/aquifer-recharge-ordering/generate.py b/probes/mass/aquifer-recharge-ordering/generate.py index 9b121e4d..62ad878c 100644 --- a/probes/mass/aquifer-recharge-ordering/generate.py +++ b/probes/mass/aquifer-recharge-ordering/generate.py @@ -3,7 +3,7 @@ import pandas as pd N_STEPS=20 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,"river_conductance_m2_per_day":100.0,"river_bottom_offset_m":1.0} +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":100.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) diff --git a/probes/mass/aquifer-recharge-ordering/probe.yaml b/probes/mass/aquifer-recharge-ordering/probe.yaml index 889e5bd4..d0c0dbc7 100644 --- a/probes/mass/aquifer-recharge-ordering/probe.yaml +++ b/probes/mass/aquifer-recharge-ordering/probe.yaml @@ -17,7 +17,7 @@ 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, river_conductance_m2_per_day] + 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 diff --git a/site/index.html b/site/index.html index 42b602da..5280e21d 100644 --- a/site/index.html +++ b/site/index.html @@ -723,7 +723,7 @@

Credit

summasubmittedSUMMA 4.0.0, a land model solving the water and energy balances of canopy, snow, soil and aquifer with one implicit solver; one lumped HRU, with the radiation and humidity it needs mocked from the forcing.9 / 24FAIL (VIOLATION) cwatmsubmittedCWatM 1.11, IIASA's Community Water Model: a process-based global hydrological model run on one grid cell, reporting its snow, soil, canopy, groundwater and runoff-concentration stores.17 / 21FAIL (VIOLATION) lisfloodsubmittedLISFLOOD 5.0.0, the EC JRC's distributed model behind EFAS and GloFAS, stepped through its own Python framework on one representative 5 km cell.18 / 21FAIL (ERROR) -modflow6submittedMODFLOW 6.7.0, the USGS's modular groundwater model with datum-aware paired inputs.0 / 0N/A (INCOMPLETE) +modflow6submittedMODFLOW 6.7.0, the USGS's modular groundwater model with datum-aware paired inputs.1 / 1PASS flex_lumpedphysicalLumped FLEX/HBV from chrimerss/HydrologicModels: interception, beta-partitioned soil, fast and slow reservoirs, a geometry-aware lag.24 / 24PASS flex_topophysicalFLEX-Topo: plateau, hillslope and wetland units on the Wark catchment's real fractions, sharing one groundwater store.23 / 23PASS sacsma_snow17physicalThe NWS operational SAC-SMA with Snow-17 and a gamma unit hydrograph, ported from the legacy Fortran and checked against it.24 / 24PASS @@ -948,7 +948,7 @@

Credit

summaenviadoSUMMA 4.0.0, un modelo de superficie terrestre que resuelve los balances de agua y energía del dosel, la nieve, el suelo y el acuífero con un único esquema implícito; una HRU agregada, con la radiación y la humedad que necesita simuladas a partir del forzamiento.9 / 24FALLA (VIOLACIÓN) cwatmenviadoCWatM 1.11, el Community Water Model del IIASA: un modelo hidrológico global basado en procesos, ejecutado sobre una sola celda, que reporta sus almacenes de nieve, suelo, dosel, agua subterránea y concentración de escorrentía.17 / 21FALLA (VIOLACIÓN) lisfloodenviadoLISFLOOD 5.0.0, el modelo distribuido del JRC de la Comisión Europea detrás de EFAS y GloFAS, avanzado paso a paso con su propio marco de Python en una celda representativa de 5 km.18 / 21FALLA (ERROR) -modflow6enviadoMODFLOW 6.7.0, el modelo modular de agua subterránea del USGS: un acuífero transitorio de una celda de río ejecutado con flopy, que reporta el intercambio señalado y el almacenamiento subterráneo a partir de sus propios términos de balance RIV y STO.0 / 0N/A (INCOMPLETE) +modflow6enviadoMODFLOW 6.7.0, el modelo modular de agua subterránea del USGS: un acuífero transitorio de una celda de río ejecutado con flopy, que reporta el intercambio señalado y el almacenamiento subterráneo a partir de sus propios términos de balance RIV y STO.1 / 1PASS flex_lumpedfísicoFLEX/HBV agregado de chrimerss/HydrologicModels: intercepción, suelo con partición beta, embalses rápido y lento, un retardo basado en la geometría.24 / 24PASA flex_topofísicoFLEX-Topo: unidades de meseta, ladera y humedal con las fracciones reales de la cuenca Wark, un solo almacén subterráneo.23 / 23PASA sacsma_snow17físicoEl SAC-SMA operativo del NWS con Snow-17 y un hidrograma unitario gamma, portado del Fortran heredado y verificado contra él.24 / 24PASA @@ -1173,7 +1173,7 @@

Credit

summa提交模型SUMMA 4.0.0，陆面模型：以同一隐式求解器求解冠层、积雪、土壤与含水层的水量和能量平衡；按单个集总 HRU 运行，所需辐射与湿度由强迫数据模拟。9 / 24未通过（违反） cwatm提交模型CWatM 1.11，IIASA 的社区水模型（Community Water Model）：基于过程的全球水文模型，在单个网格上运行，报告积雪、土壤、冠层、地下水与产流汇集滞蓄。17 / 21未通过（违反） lisflood提交模型LISFLOOD 5.0.0，欧盟委员会联合研究中心支撑 EFAS 与 GloFAS 的分布式模型，通过其自身的 Python 框架在单个 5 公里代表性网格上逐步运行。18 / 21未通过（错误） -modflow6提交模型MODFLOW 6.7.0，美国地质调查局的模块化地下水模型：通过 flopy 运行的单河流网格瞬态含水层，依据其自身的 RIV 与 STO 收支项报告带符号交换量与地下水储量。0 / 0N/A (INCOMPLETE) +modflow6提交模型MODFLOW 6.7.0，美国地质调查局的模块化地下水模型：通过 flopy 运行的单河流网格瞬态含水层，依据其自身的 RIV 与 STO 收支项报告带符号交换量与地下水储量。1 / 1PASS flex_lumped物理模型来自 chrimerss/HydrologicModels 的集总 FLEX/HBV：截留、beta 分配的土壤、快慢水库与基于几何的滞时。24 / 24通过 flex_topo物理模型FLEX-Topo：台地、坡地与湿地单元按 Wark 流域的真实面积比共享一个地下水库。23 / 23通过 sacsma_snow17物理模型美国国家气象局业务用 SAC-SMA + Snow-17 + 伽马单位线，由遗留 Fortran 移植并与之比对。24 / 24通过 From dc363ee9da0dde53262d385b992c65ceccabf393 Mon Sep 17 00:00:00 2001 From: chrimerss Date: Wed, 7 Oct 2026 11:12:53 -0600 Subject: [PATCH 5/5] Credit gw_boundary into the aquifer, pin each ordering arm, restore the site - aquifer_recharge_ordering credited a declared gw_boundary out of the aquifer (+B) where groundwater_balance has it in (dS = R + Q + B), so a correct model with a head-dependent boundary was off by twice the boundary's response: 1.7 times the tolerance already at 0.1 mm/day per m. It is now -B and closes to 1e-11 mm. The probe's groundwater_balance declares the same source. - The exchange arm was untested: the "reversal" test sent the increment the right way and failed on E_B, and deleting E_Q left the tests and the gate green. tests/test_recharge_ordering.py now builds every case on the probe's generator and checks that each arm is the one that catches its fault: a reversed increment that keeps each budget closed trips E_Q alone, and E_B is pinned as a running maximum. Removing any arm, or flipping the boundary sign, fails a test. - The criterion fails a non-finite output itself; it picked the worst arm with max() over the ratios, which passes a NaN exchange. The protocol already rejects NaN columns, so this was latent. - The criterion and the three reference adapters are rewritten in the repository's style; the exact reference now integrates the cell exactly over each day. Same case, same verdicts. - Probe README: the arms, the tolerance, the boundary, which content is independent (E_B largely restates groundwater_balance), which arms each reference trips, and what the probe does not catch. The order is non-strict, so a model that banks the pulse passes; a recovery check needs the aquifer's response time, which for MODFLOW 6's 10 x 10 grid at K = 1 m/day is thousands of days and is not derivable from the case's attributes. A test pins the gap. - site/index.html: the modflow6 models row had been replaced by "$13 / 3" in all three languages; restored as 3 / 3. The flowchart's new entry overprinted "datum"; it now follows it, the three pillars are 6 px taller, and the 24 MASS entries are evenly spaced, which also clears the spin-up/drift overlap main had. - README: thirty-six probes, twenty-four mass, three groundwater-exchange probes; reference descriptions match what the adapters do. - models/result.csv: the ten rows re-archived on the gate seeds; modflow6 still passes (E_B 7e-10 mm), the other nine are N/A as before. Co-Authored-By: Claude Opus 5.5 --- README.md | 18 +- docs/writing-a-probe.md | 6 +- models/reference_recharge_blind/ht_adapter.py | 79 +++++-- models/reference_recharge_blind/model.yaml | 2 +- .../ht_adapter.py | 79 +++++-- .../reference_recharge_order_exact/model.yaml | 2 +- .../ht_adapter.py | 79 +++++-- .../reference_recharge_overshoot/model.yaml | 2 +- models/result.csv | 20 +- .../mass/aquifer-recharge-ordering/README.md | 102 ++++++--- .../aquifer-recharge-ordering/generate.py | 54 ++++- .../mass/aquifer-recharge-ordering/probe.yaml | 1 + site/index.html | 110 ++++----- src/hydroturing/criteria/recharge_ordering.py | 160 ++++++++++--- tests/test_recharge_ordering.py | 214 ++++++++++++++++-- 15 files changed, 706 insertions(+), 222 deletions(-) diff --git a/README.md b/README.md index 447b8adb..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,7 +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 | Added direct aquifer recharge preserves storage and incremental exchange ordering; the paired budget arm restates `groundwater_balance`, while the independent content is the storage and exchange ordering. | `reference_recharge_blind`, `reference_recharge_overshoot` | +| [`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` | @@ -171,11 +171,11 @@ the probe cannot ask the declared model interface this question. | [`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**, 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 two groundwater-exchange probes that need outputs it does not report | +| `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 | passive aquifer reference that partitions direct recharge between storage and river exchange | must pass `mass/aquifer-recharge-ordering` | -| `reference_recharge_blind` | broken | ignores the added recharge while keeping an internally closed budget | caught by `aquifer_recharge_ordering` | -| `reference_recharge_overshoot` | broken | explicit update overshoots after the recharge pulse | caught by `aquifer_recharge_ordering` | +| `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` | @@ -413,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/docs/writing-a-probe.md b/docs/writing-a-probe.md index 9b28e143..58b80dba 100644 --- a/docs/writing-a-probe.md +++ b/docs/writing-a-probe.md @@ -377,9 +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` | passive aquifer reference that partitions recharge exactly | nothing, it must pass | -| `reference_recharge_blind` | ignores the added recharge while closing its own budget | `aquifer_recharge_ordering` | -| `reference_recharge_overshoot` | explicit update overshoots after the pulse | `aquifer_recharge_ordering` | +| `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 index 8f8d226b..9dc6261f 100644 --- a/models/reference_recharge_blind/ht_adapter.py +++ b/models/reference_recharge_blind/ht_adapter.py @@ -1,21 +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,csv,json + +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="exact"): - sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] - for row in forcing: - stage=float(row["sw_stage_m"]); r=float(row["gw_recharge"]); q=k*(stage-head); new=head+(r+q)/b; q=b*(new-head)-r - if mode=="blind": r_effect=0.; new=head+(r_effect+q)/b; q=b*(new-head)-r_effect - elif mode=="overshoot" and r>0: q=-1.2*r; new=head+(r+q)/b - head=new; gw=b*head - rows.append({"time":row["time"],"gw_sw_exchange":q,"gw_to_sw":min(q,0.),"sw_to_gw":max(q,0.),"gw":gw}) - return rows +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(): - p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent - with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) - static=json.loads((d/req["input"]["static"]).read_text()); rows=simulate(forcing,static,MODE) - with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) - (d/req["output"]["run"]).write_text(json.dumps({"status":"ok"})) -if __name__=="__main__": 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 index 07818d42..6a13423a 100644 --- a/models/reference_recharge_blind/model.yaml +++ b/models/reference_recharge_blind/model.yaml @@ -1,6 +1,6 @@ name: reference_recharge_blind version: "0.1.0" -description: Recharge ordering reference. +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 diff --git a/models/reference_recharge_order_exact/ht_adapter.py b/models/reference_recharge_order_exact/ht_adapter.py index 042b46fe..f570f3be 100644 --- a/models/reference_recharge_order_exact/ht_adapter.py +++ b/models/reference_recharge_order_exact/ht_adapter.py @@ -1,21 +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,csv,json + +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="exact"): - sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] - for row in forcing: - stage=float(row["sw_stage_m"]); r=float(row["gw_recharge"]); q=k*(stage-head); new=head+(r+q)/b; q=b*(new-head)-r - if mode=="blind": r_effect=0.; new=head+(r_effect+q)/b; q=b*(new-head)-r_effect - elif mode=="overshoot" and r>0: q=-1.2*r; new=head+(r+q)/b - head=new; gw=b*head - rows.append({"time":row["time"],"gw_sw_exchange":q,"gw_to_sw":min(q,0.),"sw_to_gw":max(q,0.),"gw":gw}) - return rows +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(): - p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent - with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) - static=json.loads((d/req["input"]["static"]).read_text()); rows=simulate(forcing,static,MODE) - with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) - (d/req["output"]["run"]).write_text(json.dumps({"status":"ok"})) -if __name__=="__main__": 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 index b7d01ec0..0f9f4679 100644 --- a/models/reference_recharge_order_exact/model.yaml +++ b/models/reference_recharge_order_exact/model.yaml @@ -1,6 +1,6 @@ name: reference_recharge_order_exact version: "0.1.0" -description: Recharge ordering reference. +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 diff --git a/models/reference_recharge_overshoot/ht_adapter.py b/models/reference_recharge_overshoot/ht_adapter.py index 5c0fe3c2..9cb39d2d 100644 --- a/models/reference_recharge_overshoot/ht_adapter.py +++ b/models/reference_recharge_overshoot/ht_adapter.py @@ -1,21 +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,csv,json + +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="exact"): - sy=float(static["aquifer_specific_yield"]); area=float(static["area_km2"])*1e6; b=sy*1000.; k=float(static["river_conductance_m2_per_day"])/area*1000.; head=float(static["aquifer_initial_head_m"]); gw=b*head; rows=[] - for row in forcing: - stage=float(row["sw_stage_m"]); r=float(row["gw_recharge"]); q=k*(stage-head); new=head+(r+q)/b; q=b*(new-head)-r - if mode=="blind": r_effect=0.; new=head+(r_effect+q)/b; q=b*(new-head)-r_effect - elif mode=="overshoot" and r>0: q=-1.2*r; new=head+(r+q)/b - head=new; gw=b*head - rows.append({"time":row["time"],"gw_sw_exchange":q,"gw_to_sw":min(q,0.),"sw_to_gw":max(q,0.),"gw":gw}) - return rows +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(): - p=argparse.ArgumentParser(); p.add_argument("--request",required=True); a=p.parse_args(); path=Path(a.request).resolve(); req=json.loads(path.read_text()); d=path.parent - with open(d/req["input"]["forcing"],newline="") as f: forcing=list(csv.DictReader(f)) - static=json.loads((d/req["input"]["static"]).read_text()); rows=simulate(forcing,static,MODE) - with open(d/req["output"]["table"],"w",newline="") as f: w=csv.DictWriter(f,fieldnames=COLUMNS); w.writeheader(); w.writerows(rows) - (d/req["output"]["run"]).write_text(json.dumps({"status":"ok"})) -if __name__=="__main__": 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 index e719313b..37e39966 100644 --- a/models/reference_recharge_overshoot/model.yaml +++ b/models/reference_recharge_overshoot/model.yaml @@ -1,6 +1,6 @@ name: reference_recharge_overshoot version: "0.1.0" -description: Recharge ordering reference. +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 diff --git a/models/result.csv b/models/result.csv index 84c035e0..78238657 100644 --- a/models/result.csv +++ b/models/result.csv @@ -780,13 +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-05,modflow6,6.7.0-adapter.2,0.1.0,docker,mass/aquifer-recharge-ordering,PASS,OK,full record,428911210 732642421 1064022182,aquifer_recharge_ordering: added recharge preserves storage and exchange ordering (I=10 mm) -2026-10-05,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-05,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-05,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-05,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-05,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-05,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-05,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-05,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-05,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" +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 index 45180892..e83ee8bb 100644 --- a/probes/mass/aquifer-recharge-ordering/README.md +++ b/probes/mass/aquifer-recharge-ordering/README.md @@ -1,28 +1,78 @@ # mass/aquifer-recharge-ordering -The paired case adds only direct `gw_recharge` to the same passive aquifer. The -`aquifer_recharge_ordering` criterion evaluates storage ordering, the sign of -the incremental exchange, and the cumulative partition identity using the -added recharge volume as its scale. A plateau with zero incremental exchange -is valid; a recharge-blind or nonmonotone update is not. - -The exact reference supplies the deterministic gate. The independent MODFLOW 6 baseline is enabled for paired perturbations and is archived separately; the four standard catchment references do not expose a direct aquifer recharge interface. - -## Criterion and tolerance - -For each scored step, let `I_n` be cumulative added recharge and `D_n` be -cumulative incremental exchange (positive when the added recharge sends more -water to the river). The criterion records three one-sided errors: - -* `E_S = max([-ΔS_n]_+)`, the largest storage decrease; -* `E_Q = sum([-d_i]_+ · Δt_i)`, the total exchange increment in the wrong - direction; and -* `E_B = max(|ΔS_n + D_n + B_n - I_n|)`, the worst cumulative partition - residual, where `B_n` is the cumulative declared boundary-source change. - -The common tolerance is `ε = max(0.01 · I_* , 1e-3 mm, f · S_max)` where -`I_*` is total added recharge, `S_max` is the largest absolute storage in -either run, and `f · S_max` is a capped storage precision floor (`f = 1e-5`, -cap `0.05 mm`). Every error must be no larger than `ε`. This makes the scale -come from the perturbation itself while accounting for storage rounding without -allowing a large absolute datum to buy unlimited tolerance. +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 index 63f31984..d4974f74 100644 --- a/probes/mass/aquifer-recharge-ordering/generate.py +++ b/probes/mass/aquifer-recharge-ordering/generate.py @@ -1,15 +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 -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) - stage=100.0+0.01*rng.normal(size=N_STEPS) - recharge=np.zeros(N_STEPS) - if variant=="recharge_added": recharge[2:4]=5.0 - 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) +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 index d0c0dbc7..fdb336c1 100644 --- a/probes/mass/aquifer-recharge-ordering/probe.yaml +++ b/probes/mass/aquifer-recharge-ordering/probe.yaml @@ -42,6 +42,7 @@ criteria: 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 diff --git a/site/index.html b/site/index.html index 8776d16e..3be98eef 100644 --- a/site/index.html +++ b/site/index.html @@ -350,60 +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 - - recharge order - - 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 @@ -424,7 +424,7 @@

How it works

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

How it works

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

Credit

"infra.probe.snowpack":"snowpack", "infra.probe.ungauged":"ungauged", "infra.probe.exchange":"exchange", - "infra.probe.recharge":"recharge order", + "infra.probe.recharge":"recharge", "infra.probe.datum":"datum shift", "infra.probe.gwexchange":"gw exchange", "infra.probe.pet":"PET demand", @@ -869,7 +869,7 @@

Credit

summasubmittedSUMMA 4.0.0, a land model solving the water and energy balances of canopy, snow, soil and aquifer with one implicit solver; one lumped HRU, with the radiation and humidity it needs mocked from the forcing.9 / 24FAIL (VIOLATION) cwatmsubmittedCWatM 1.11, IIASA's Community Water Model: a process-based global hydrological model run on one grid cell, reporting its snow, soil, canopy, groundwater and runoff-concentration stores.17 / 21FAIL (VIOLATION) lisfloodsubmittedLISFLOOD 5.0.0, the EC JRC's distributed model behind EFAS and GloFAS, stepped through its own Python framework on one representative 5 km cell.18 / 21FAIL (ERROR) -$13 / 3PASS +modflow6submittedMODFLOW 6.7.0, the USGS's modular groundwater model: a transient one-river-cell aquifer run through flopy, reporting signed exchange and groundwater storage from its own RIV and STO budget terms.3 / 3PASS flex_lumpedphysicalLumped FLEX/HBV from chrimerss/HydrologicModels: interception, beta-partitioned soil, fast and slow reservoirs, a geometry-aware lag.25 / 25PASS flex_topophysicalFLEX-Topo: plateau, hillslope and wetland units on the Wark catchment's real fractions, sharing one groundwater store.24 / 24PASS sacsma_snow17physicalThe NWS operational SAC-SMA with Snow-17 and a gamma unit hydrograph, ported from the legacy Fortran and checked against it.25 / 25PASS @@ -1034,7 +1034,7 @@

Credit

"infra.probe.snowpack":"manto nival", "infra.probe.ungauged":"sin aforo", "infra.probe.exchange":"intercambio", - "infra.probe.recharge":"orden de recarga", + "infra.probe.recharge":"recarga", "infra.probe.datum":"datum vert.", "infra.probe.gwexchange":"intercambio subt.", "infra.probe.pet":"demanda", @@ -1115,7 +1115,7 @@

Credit

summaenviadoSUMMA 4.0.0, un modelo de superficie terrestre que resuelve los balances de agua y energía del dosel, la nieve, el suelo y el acuífero con un único esquema implícito; una HRU agregada, con la radiación y la humedad que necesita simuladas a partir del forzamiento.9 / 24FALLA (VIOLACIÓN) cwatmenviadoCWatM 1.11, el Community Water Model del IIASA: un modelo hidrológico global basado en procesos, ejecutado sobre una sola celda, que reporta sus almacenes de nieve, suelo, dosel, agua subterránea y concentración de escorrentía.17 / 21FALLA (VIOLACIÓN) lisfloodenviadoLISFLOOD 5.0.0, el modelo distribuido del JRC de la Comisión Europea detrás de EFAS y GloFAS, avanzado paso a paso con su propio marco de Python en una celda representativa de 5 km.18 / 21FALLA (ERROR) -$13 / 3PASA +modflow6enviadoMODFLOW 6.7.0, el modelo modular de agua subterránea del USGS: un acuífero transitorio de una celda de río ejecutado con flopy, que reporta el intercambio señalado y el almacenamiento subterráneo a partir de sus propios términos de balance RIV y STO.3 / 3PASA flex_lumpedfísicoFLEX/HBV agregado de chrimerss/HydrologicModels: intercepción, suelo con partición beta, embalses rápido y lento, un retardo basado en la geometría.25 / 25PASA flex_topofísicoFLEX-Topo: unidades de meseta, ladera y humedal con las fracciones reales de la cuenca Wark, un solo almacén subterráneo.24 / 24PASA sacsma_snow17físicoEl SAC-SMA operativo del NWS con Snow-17 y un hidrograma unitario gamma, portado del Fortran heredado y verificado contra él.25 / 25PASA @@ -1361,7 +1361,7 @@

Credit

summa提交模型SUMMA 4.0.0，陆面模型：以同一隐式求解器求解冠层、积雪、土壤与含水层的水量和能量平衡；按单个集总 HRU 运行，所需辐射与湿度由强迫数据模拟。9 / 24未通过（违反） cwatm提交模型CWatM 1.11，IIASA 的社区水模型（Community Water Model）：基于过程的全球水文模型，在单个网格上运行，报告积雪、土壤、冠层、地下水与产流汇集滞蓄。17 / 21未通过（违反） lisflood提交模型LISFLOOD 5.0.0，欧盟委员会联合研究中心支撑 EFAS 与 GloFAS 的分布式模型，通过其自身的 Python 框架在单个 5 公里代表性网格上逐步运行。18 / 21未通过（错误） -$13 / 3通过 +modflow6提交模型MODFLOW 6.7.0，美国地质调查局的模块化地下水模型：通过 flopy 运行的单河流网格瞬态含水层，依据其自身的 RIV 与 STO 收支项报告带符号交换量与地下水储量。3 / 3通过 flex_lumped物理模型来自 chrimerss/HydrologicModels 的集总 FLEX/HBV：截留、beta 分配的土壤、快慢水库与基于几何的滞时。25 / 25通过 flex_topo物理模型FLEX-Topo：台地、坡地与湿地单元按 Wark 流域的真实面积比共享一个地下水库。24 / 24通过 sacsma_snow17物理模型美国国家气象局业务用 SAC-SMA + Snow-17 + 伽马单位线，由遗留 Fortran 移植并与之比对。25 / 25通过 diff --git a/src/hydroturing/criteria/recharge_ordering.py b/src/hydroturing/criteria/recharge_ordering.py index 85fcfaa6..427f05ed 100644 --- a/src/hydroturing/criteria/recharge_ordering.py +++ b/src/hydroturing/criteria/recharge_ordering.py @@ -1,45 +1,135 @@ -"""Order preservation for a nonnegative direct aquifer-recharge pulse.""" +"""Ordering for a direct aquifer-recharge pulse. + +`aquifer_recharge_ordering` compares a control run with one that differs only +by a nonnegative pulse of `gw_recharge` (mass/aquifer-recharge-ordering, #135). +It is the groundwater comparison principle: for a passive aquifer whose +outflow does not fall as its storage rises, adding recharge can never leave +less water in it, and can never make the net exchange less favourable to the +river. 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 any declared aquifer boundary flux (positive into the aquifer, the +sign `groundwater_balance` gives it) and dS_n the paired storage difference. +Three one-sided errors are scored: + + E_S = max_n [-dS_n]_+ storage never falls below control + E_Q = sum_i [-d_i]_+ dt the exchange increment never reverses + E_B = max_n |dS_n + D_n - B_n - I_n| the added water is all accounted for + +E_S is a maximum because one excursion below the control breaks the order; +E_Q is a sum because every reversed increment is water sent the wrong way; +E_B is a running maximum so that a loss cannot be cancelled by a later gain. +If each run closes its own budget, which `groundwater_balance` checks on both +variants, E_B is zero to within the two tolerances: it restates that check as +a pair. The independent content of the criterion is E_S and E_Q. + +All three share one tolerance scaled to the intervention rather than to the +background storage, eps = max(rel_tol * I_*, abs_tol_mm, min(1e-5 * S_max, +0.05 mm)). The last term absorbs rounding in a large absolute `gw` without +letting a datum buy unlimited slack, the same capped floor `groundwater_balance` +uses. + +The order is non-strict: a zero exchange increment passes, so a model that +banks the added recharge and never releases it is not caught here. Asking +that the water drain back out needs the aquifer's response time, and for a +distributed aquifer that depends on a transmissivity and a geometry the case +does not supply. +""" + from __future__ import annotations + import numpy as np -from hydroturing.criteria.base import FAIL, PASS, CriterionResult, criterion, make_window + +from hydroturing.criteria.base import ( + FAIL, + PASS, + CriterionResult, + criterion, + make_window, +) from hydroturing.protocol import RunResult from hydroturing.spec import ProbeSpec + @criterion("aquifer_recharge_ordering", paired=True) def aquifer_recharge_ordering(runs: dict[str, RunResult], probe: ProbeSpec, params: dict) -> CriterionResult: - """Check storage ordering, exchange ordering and cumulative partitioning. - - E_S is a maximum because one negative excursion violates comparison; E_Q - is a sum because any reversed increment is water sent in the wrong - direction; E_B is a maximum so local losses cannot cancel later. The - storage floor is capped to protect absolute-storage rounding without - letting a large datum buy unlimited tolerance. - """ - control_name = str(params.get("control", "control")); perturbed_name = str(params.get("perturbed", "recharge_added")) - if control_name not in runs or perturbed_name not in runs: raise ValueError("aquifer_recharge_ordering needs control and perturbed variants") - control = make_window(runs[control_name], probe); perturbed = make_window(runs[perturbed_name], probe) - recharge = str(params.get("recharge", "gw_recharge")); exchange = str(params.get("exchange", "gw_sw_exchange")); storage = str(params.get("storage", "gw")) - for window, label in ((control, "control"), (perturbed, "perturbed")): - if recharge not in window.forcing.columns: raise ValueError(f"{label} forcing needs '{recharge}'") - for name in (exchange, storage): - if name not in window.table.columns: raise ValueError(f"{label} result needs '{name}'") - if len(control.table) != len(perturbed.table): raise ValueError("paired recharge runs have different scored lengths") + """Added recharge never lowers storage, never reverses the exchange + increment, and is partitioned between them without loss or creation.""" + name = "aquifer_recharge_ordering" + control_name = str(params.get("control", "control")) + perturbed_name = str(params.get("perturbed", "recharge_added")) + for variant in (control_name, perturbed_name): + if variant not in runs: + raise ValueError(f"{name} needs variant {variant!r}") + control = make_window(runs[control_name], probe) + perturbed = make_window(runs[perturbed_name], probe) + recharge = str(params.get("recharge", "gw_recharge")) + exchange = str(params.get("exchange", "gw_sw_exchange")) + storage = str(params.get("storage", "gw")) + for window, label in ((control, control_name), (perturbed, perturbed_name)): + if recharge not in window.forcing.columns: + raise ValueError(f"{name}: {label} forcing needs '{recharge}'") + for column in (exchange, storage): + if column not in window.table.columns: + raise ValueError(f"{name}: {label} result needs '{column}'") + if len(control.table) != len(perturbed.table): + raise ValueError(f"{name}: paired runs have different scored lengths") + added = perturbed.volume(perturbed.forcing[recharge]) - control.volume(control.forcing[recharge]) - tol = float(params.get("input_tol_mm", 1.0e-10)) - if np.any(added < -tol): raise ValueError("aquifer_recharge_ordering requires nonnegative added recharge at every step") - cumulative_added = np.cumsum(added); added_total = float(cumulative_added[-1]) - if added_total <= tol: raise ValueError("aquifer_recharge_ordering needs positive added recharge") - delta_storage = perturbed.table[storage].to_numpy(float) - control.table[storage].to_numpy(float) - delta_storage -= float(perturbed.state0[storage]) - float(control.state0[storage]) - d_exchange = control.volume(control.table[exchange]) - perturbed.volume(perturbed.table[exchange]); cumulative_exchange = np.cumsum(d_exchange) - source_delta = np.zeros(len(control.table)) + input_tol = float(params.get("input_tol_mm", 1.0e-10)) + # The order is only implied by physics when water is added at every step; + # a variant that removes some is a different experiment. + if np.any(added < -input_tol): + raise ValueError(f"{name} requires nonnegative added recharge at every step") + added_total = float(np.sum(added)) + if added_total <= input_tol: + raise ValueError(f"{name} needs positive added recharge") + + s_control = control.table[storage].to_numpy(dtype=float) + s_perturbed = perturbed.table[storage].to_numpy(dtype=float) + s0 = float(perturbed.state0[storage]) - float(control.state0[storage]) + delta_storage = (s_perturbed - s_control) - s0 + d_exchange = control.volume(control.table[exchange]) - perturbed.volume(perturbed.table[exchange]) + d_boundary = np.zeros(len(control.table)) for source in params.get("sources", []): if source in control.table.columns and source in perturbed.table.columns: - source_delta += perturbed.volume(perturbed.table[source]) - control.volume(control.table[source]) - storage_error = float(np.max(np.maximum(-delta_storage, 0.0))); exchange_error = float(np.sum(np.maximum(-d_exchange, 0.0))) - budget_error = float(np.max(np.abs(delta_storage + cumulative_exchange + np.cumsum(source_delta) - cumulative_added))) - relative = float(params.get("rel_tol", 0.01)) * added_total; absolute = float(params.get("abs_tol_mm", 1.0e-3)) - storage_max = max(float(np.max(np.abs(control.table[storage]))), float(np.max(np.abs(perturbed.table[storage])))); storage_floor = min(1.0e-5 * storage_max, 0.05) - epsilon = max(relative, absolute, storage_floor) - ratios = {"storage": storage_error / epsilon, "exchange": exchange_error / epsilon, "budget": budget_error / epsilon}; worst_name, worst = max(ratios.items(), key=lambda item: item[1]); ok = worst <= 1.0 - return CriterionResult(name="aquifer_recharge_ordering", status=PASS if ok else FAIL, value=worst, threshold=1.0, message=(f"added recharge preserves storage and exchange ordering (I={added_total:.6g} mm)" if ok else f"recharge ordering fails {worst_name}: {max(storage_error, exchange_error, budget_error):.6g} mm exceeds {epsilon:.6g} mm"), diagnostics={"added_recharge_mm": added_total, "epsilon_mm": epsilon, "E_S_mm": storage_error, "E_Q_mm": exchange_error, "E_B_mm": budget_error, "ratios": ratios}) + d_boundary += perturbed.volume(perturbed.table[source]) - control.volume(control.table[source]) + # Checked here rather than left to the comparisons below: the worst arm is + # picked with max(), which passes over a NaN ratio. + if not (np.isfinite(delta_storage).all() and np.isfinite(d_exchange).all() and np.isfinite(d_boundary).all()): + return CriterionResult(name=name, status=FAIL, threshold=1.0, message="non-finite storage, exchange or boundary flux") + + storage_max = float(np.max(np.abs(np.concatenate([s_control, s_perturbed])))) + epsilon = max( + float(params.get("rel_tol", 0.01)) * added_total, + float(params.get("abs_tol_mm", 1.0e-3)), + min(1.0e-5 * storage_max, 0.05), + ) + budget = delta_storage + np.cumsum(d_exchange) - np.cumsum(d_boundary) - np.cumsum(added) + errors = { + "storage": float(np.max(np.maximum(-delta_storage, 0.0))), + "exchange": float(np.sum(np.maximum(-d_exchange, 0.0))), + "budget": float(np.max(np.abs(budget))), + } + ratios = {arm: error / epsilon for arm, error in errors.items()} + worst_arm = max(ratios, key=ratios.get) + ok = ratios[worst_arm] <= 1.0 + return CriterionResult( + name=name, + status=PASS if ok else FAIL, + value=ratios[worst_arm], + threshold=1.0, + message=( + f"added recharge preserves storage and exchange ordering (I={added_total:.6g} mm)" + if ok + else f"recharge ordering fails {worst_arm}: {errors[worst_arm]:.6g} mm exceeds {epsilon:.6g} mm" + ), + diagnostics={ + "added_recharge_mm": added_total, + "epsilon_mm": epsilon, + "E_S_mm": errors["storage"], + "E_Q_mm": errors["exchange"], + "E_B_mm": errors["budget"], + "ratios": ratios, + }, + ) diff --git a/tests/test_recharge_ordering.py b/tests/test_recharge_ordering.py index 2e5f28c5..102b6a88 100644 --- a/tests/test_recharge_ordering.py +++ b/tests/test_recharge_ordering.py @@ -1,26 +1,202 @@ +"""mass/aquifer-recharge-ordering: each arm caught on its own. + +The gate shows that the criterion separates its references, but a must_fail +baseline names a criterion, not an arm of it, and no reference reverses the +exchange. These tests build runs on the probe's own generator and check that +each arm of `aquifer_recharge_ordering` is the one that catches the fault it is +there for; a reversed exchange increment is caught by E_Q alone. Deleting any +one arm, or flipping the boundary's sign, fails a test here. +""" + from __future__ import annotations -import numpy as np, pandas as pd, pytest + +import runpy from dataclasses import replace +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + from hydroturing import registry from hydroturing.criteria.recharge_ordering import aquifer_recharge_ordering -from hydroturing.harness import build_case +from hydroturing.harness import build_case, evaluate_criteria from hydroturing.protocol import RunResult +from hydroturing.seeds import gate_seeds + +REPO_ROOT = Path(__file__).resolve().parents[1] +PROBE = registry.find_probe("mass/aquifer-recharge-ordering") +SEED = gate_seeds(PROBE.id, 1)[0] + + +def _reference(mode: str): + # The three reference_recharge_* adapters are one file that differs only in MODE. + simulate = runpy.run_path(str(REPO_ROOT / "models/reference_recharge_order_exact/ht_adapter.py"))["simulate"] + return lambda forcing, static: simulate(forcing, static, mode) + + +def _tables(mode: str = "exact") -> dict[str, tuple]: + simulate = _reference(mode) + out = {} + for variant in PROBE.variants: + case = build_case(PROBE, SEED, variant) + table = pd.DataFrame(simulate(case.forcing.to_dict("records"), case.static)).drop(columns="time") + out[variant] = (case, table) + return out + + +def _runs(tables) -> dict[str, RunResult]: + runs = {} + for variant, (case, table) in tables.items(): + q = table["gw_sw_exchange"] + table = table.assign(gw_to_sw=np.minimum(q, 0.0), sw_to_gw=np.maximum(q, 0.0)) + runs[variant] = RunResult(case=case, table=table, meta={}, wall_seconds=0.0) + return runs + + +def _params(name: str) -> dict: + return dict(next(c.params for c in PROBE.criteria if c.name == name)) + + +def _ordering(runs): + return aquifer_recharge_ordering(runs, PROBE, _params("aquifer_recharge_ordering")) + + +def _others(runs) -> dict[str, bool]: + results = evaluate_criteria(runs, PROBE, control="control") + return {r.name: r.passed for r in results if r.name != "aquifer_recharge_ordering"} + + +def test_exact_reference_passes_every_criterion(): + runs = _runs(_tables()) + ordering = _ordering(runs) + assert ordering.passed + assert ordering.diagnostics["E_B_mm"] < 1e-9 + assert all(_others(runs).values()) + + +def test_storage_arm_catches_an_overshooting_update(): + runs = _runs(_tables("overshoot")) + result = _ordering(runs) + assert not result.passed + ratios = result.diagnostics["ratios"] + # E_S, about 20 times the tolerance. E_Q follows at about 1.5 as the river + # refills the deficit, which is a reversed increment too. + assert ratios["storage"] > 10.0 + assert ratios["storage"] == max(ratios.values()) + assert ratios["budget"] <= 1.0 + + +def test_exchange_arm_catches_a_reversed_increment_the_budget_cannot_see(): + tables = _tables() + control = tables["control"][1] + case, perturbed = tables["recharge_added"] + # Send the added water's exchange response the wrong way, and move gw by + # the same amount so this run still closes its own budget. + q = control["gw_sw_exchange"] - (perturbed["gw_sw_exchange"] - control["gw_sw_exchange"]) + gw = perturbed["gw"] + np.cumsum(q - perturbed["gw_sw_exchange"]) + tables["recharge_added"] = (case, perturbed.assign(gw_sw_exchange=q, gw=gw)) + runs = _runs(tables) + result = _ordering(runs) + assert not result.passed + ratios = result.diagnostics["ratios"] + assert ratios["exchange"] > 1.0 + assert ratios["storage"] <= 1.0 and ratios["budget"] <= 1.0 + assert all(_others(runs).values()) + + +def test_budget_arm_catches_a_recharge_blind_model(): + runs = _runs(_tables("blind")) + result = _ordering(runs) + assert not result.passed + ratios = result.diagnostics["ratios"] + assert ratios["budget"] > 1.0 + assert ratios["storage"] <= 1.0 and ratios["exchange"] <= 1.0 + + +def test_the_budget_arm_is_a_running_maximum(): + """A loss that a later gain cancels is still a loss: 1 mm missing from + the perturbed storage for five days, then back.""" + tables = _tables() + case, table = tables["recharge_added"] + gw = table["gw"].copy() + gw.iloc[case.spinup_steps + 5 : case.spinup_steps + 10] -= 1.0 + tables["recharge_added"] = (case, table.assign(gw=gw)) + result = _ordering(_runs(tables)) + assert not result.passed + ratios = result.diagnostics["ratios"] + assert ratios["budget"] > 1.0 + assert ratios["storage"] <= 1.0 and ratios["exchange"] <= 1.0 + + +def test_a_store_that_keeps_the_pulse_is_outside_this_probe(): + """The order is non-strict, so an exchange that ignores the aquifer's head + passes: it banks the pulse and never releases it. Pinned so that a + recovery check, when one is added, shows up here as a change of scope.""" + tables = {} + for variant in PROBE.variants: + case = build_case(PROBE, SEED, variant) + s = case.static + b = s["aquifer_specific_yield"] * 1000.0 + k = s["river_conductance_m2_per_day"] / (s["area_km2"] * 1.0e6) * 1000.0 + q = k * (case.forcing["sw_stage_m"].to_numpy() - s["aquifer_initial_head_m"]) + gw0 = b * (s["aquifer_initial_head_m"] - s["aquifer_bottom_m"]) + gw = gw0 + np.cumsum(case.forcing["gw_recharge"].to_numpy() + q) + tables[variant] = (case, pd.DataFrame({"gw_sw_exchange": q, "gw": gw})) + runs = _runs(tables) + result = _ordering(runs) + assert result.passed + assert result.diagnostics["E_Q_mm"] == 0.0 + assert all(_others(runs).values()) + + +def _with_boundary(sign: float): + """The exact cell plus a head-dependent boundary (a GHB) draining it.""" + tables = {} + for variant in PROBE.variants: + case = build_case(PROBE, SEED, variant) + s = case.static + b = s["aquifer_specific_yield"] * 1000.0 + k = s["river_conductance_m2_per_day"] / (s["area_km2"] * 1.0e6) * 1000.0 + g = 0.2 * k + head0 = head = s["aquifer_initial_head_m"] + rows = [] + for stage, r in zip(case.forcing["sw_stage_m"], case.forcing["gw_recharge"]): + new = (b * head + r + k * stage + g * head0) / (b + k + g) + q, boundary = k * (stage - new), g * (head0 - new) + head = new + rows.append({"gw_sw_exchange": q, "gw_boundary": sign * boundary, "gw": b * (head - s["aquifer_bottom_m"])}) + tables[variant] = (case, pd.DataFrame(rows)) + return _runs(tables) + + +def test_a_declared_boundary_is_credited_into_the_aquifer(): + runs = _with_boundary(+1.0) + assert _ordering(runs).passed + assert all(_others(runs).values()) + + +def test_a_boundary_reported_with_the_wrong_sign_fails_the_budget(): + result = _ordering(_with_boundary(-1.0)) + assert not result.passed + assert result.diagnostics["ratios"]["budget"] > 1.0 + + +def test_removing_recharge_is_refused_rather_than_judged(): + tables = _tables() + case, table = tables["recharge_added"] + forcing = case.forcing.copy() + forcing.loc[5, "gw_recharge"] = -1.0 + tables["recharge_added"] = (replace(case, forcing=forcing), table) + with pytest.raises(ValueError, match="nonnegative"): + _ordering(_runs(tables)) -def make(qc, qp, sc, sp, rc=None, rp=None): - probe=registry.find_probe('mass/aquifer-recharge-ordering'); case=build_case(probe,7,'control'); control_case=case; n=len(case.forcing) - def run(q,s,r,run_case): - forcing=run_case.forcing.copy(); forcing['gw_recharge']=np.asarray(r if r is not None else np.zeros(n)); table=pd.DataFrame({'gw_sw_exchange':q,'gw':s,'gw_to_sw':np.minimum(q,0),'sw_to_gw':np.maximum(q,0)}); return RunResult(case=replace(run_case, forcing=forcing),table=table,meta={},wall_seconds=0.) - control=run(np.asarray(qc),np.asarray(sc),rc if rc is not None else np.zeros(n),control_case) - # Rebuild the case so the two forcing frames do not alias. - case=build_case(probe,7,'recharge_added') - perturbed=run(np.asarray(qp),np.asarray(sp),rp if rp is not None else np.zeros(n),case) - return {'control':control,'recharge_added':perturbed},probe - -def params(probe): return dict(next(c.params for c in probe.criteria if c.name=='aquifer_recharge_ordering')) -def test_recharge_ordering_rejects_exchange_reversal(): - n=len(build_case(registry.find_probe('mass/aquifer-recharge-ordering'),7,'control').forcing); r=np.zeros(n); r[2:4]=5; q=np.zeros(n); qp=q.copy(); qp[2]=-1; s=np.full(n,2000.); sp=s.copy(); sp[2:]=sp[2:]+np.cumsum(r)[2:] - runs,p=make(q,qp,s,sp,np.zeros(n),r); assert not aquifer_recharge_ordering(runs,p,params(p)).passed -def test_recharge_ordering_rejects_negative_input_increment(): - n=len(build_case(registry.find_probe('mass/aquifer-recharge-ordering'),7,'control').forcing); r=np.zeros(n); r[2]=5; r[3]=-1; q=np.zeros(n); s=np.full(n,2000.); runs,p=make(q,q,s,s,r,np.zeros(n)); - with pytest.raises(ValueError,match='nonnegative'): aquifer_recharge_ordering(runs,p,params(p)) +@pytest.mark.parametrize("column", ["gw", "gw_sw_exchange"]) +def test_a_non_finite_output_fails(column): + tables = _tables() + case, table = tables["recharge_added"] + table = table.copy() + table.loc[10, column] = np.nan + tables["recharge_added"] = (case, table) + assert not _ordering(_runs(tables)).passed