[PROBE: momentum] Add state-dependent flood-wave celerity check

[mimosapudical]

What this probe asserts

Adds momentum/wave-celerity-bounds, which tests whether a model's local transient response around a hydraulic state obeys the same physics as that state.

Existing probes can verify mass conservation, causality, increasing lag with reach length, and steady Manning consistency while still accepting a router whose transient propagation speed is fixed.

For example, a model can satisfy

longer reach -> larger lag

while propagating low, medium, and high hydraulic states at the same fixed speed.

At low, medium, and high prescribed steady river inflows, this probe runs paired short and long reaches with identical q_in and identical static inputs except reach_length_m. A small transient pulse is applied around each operating state.

Because both responses are measured at the outlet, the paired centroid-time difference identifies the observed propagation speed:

c_obs = (L_long - L_short) / (t_centroid,long - t_centroid,short)

The result is compared with the local rectangular-channel Manning/kinematic expectation:

c_kin = dQ/dA

The observed celerity must:

• remain positive;
• stay within the pre-declared 5% allowance of c_kin;
• increase from low to medium to high hydraulic state.

The temporal centroid is used instead of peak lag because dispersive routing can move the hydrograph peak even when the propagation timescale is correct. A Hayami regression test verifies that the first temporal moment still recovers the intended celerity under substantial diffusion.

The PR also adds a narrow Wflow routing path:

q_in -> native river routing -> dis

This isolates channel-wave propagation from rainfall and land-storage dynamics. Wflow's native kinematic-wave solver produced:

0.704 < 0.927 < 1.225 m/s

across the three hydraulic states, without changing Wflow's routing equations or relaxing the probe tolerance.

Discrimination

Reference model	Expected	Criterion that catches it
reference_saint_venant	PASS	n/a
reference_fixed_celerity	FAIL	wave_celerity_bounds

reference_saint_venant advances one-dimensional continuity and momentum with an independent finite-volume dynamic-wave implementation.

reference_fixed_celerity deliberately routes every hydraulic state at the same fixed speed. It remains causal and preserves length-dependent lag, but fails the state-dependent celerity criterion.

The gate additionally requires the negative control to fail only wave_celerity_bounds, preventing an unrelated failure from producing false discrimination.

Checklist

• There is an accepted proposal issue and this PR closes it
• authors in probe.yaml includes name, affiliation, GitHub handle, and ORCID; CONTRIBUTORS.md and CITATION.cff are in sync
• ht validate passes on the current PR head
• ht gate --probe momentum/wave-celerity-bounds passes on the current PR head
• The generator is deterministic given a seed and commits no data
• The focused probe gate runs in under one minute
• Tolerance and denominator are justified in probe.yaml
• A deliberately unphysical fixed-celerity model is caught specifically by wave_celerity_bounds

[chrimerss]

/review

[github-actions]

Automated review by Claude (changes-suggested), not a maintainer approval.

Adds momentum/wave-celerity-bounds: a paired short (4 km) / long (20 km) reach driven by prescribed q_in at 8/16/32 m3/s, each perturbed by a 6-hour +5% pulse. The new wave_celerity_bounds criterion differences the two outlet response centroids, so any timing common to both variants cancels, and requires the implied celerity to be positive, within 5% of exact rectangular Manning dQ/dA evaluated at q_in, and increasing with flow. Ships a fixed-1 m/s negative control, a sub-daily transient path for reference_saint_venant on a length-independent ~80 m mesh, a two-cell native river-routing path for wflow_sbm, and the sum_inflow m3/s→mm conversion in closure.

The physics is right: for a diffusion wave the first temporal moment is exactly L/c, which is why the centroid beats the peak, and the pre-pulse settling precondition is genuinely enforced rather than documented. Docs bookkeeping (README/ROADMAP/CONTRIBUTORS/CITATION/site/result.csv, wflow 19 of 22) all check out against the in-repo tests.

Two things for a maintainer. First, the estimator validates the record before the pulse rigorously but not after: the response is truncated at 72 h with no completeness check, so a slow-tailed router gets a silently biased celerity (finding 1). Second, the only independent physical evidence is Wflow; the README reports LISFLOOD implying 12.9–63.1 m/s against 0.66–1.00 m/s expected, and that 40–60× disagreement is parked as N/A rather than diagnosed. Also, ht validate and ht gate are unticked in the checklist — CI decides those.

_{Posted by the pr-review workflow. Run}

[cehw]

Thanks, Jingzhi — the paired-reach design helps separate propagation from a shared upstream delay, and the fixed-speed control makes the added state-response check clear.

I ran 80fff8b on Linux and found one reproducibility blocker: probes/momentum/wave-celerity-bounds/probe.yaml supplies the ORCID as a URL, but the schema requires the bare identifier. This also prevents other probes from loading through the registry. Please change that field to 0009-0001-4031-8545; the URL in CITATION.cff can stay. With only that correction in my review copy, ht validate, 91 focused tests and the three-seed probe gate pass.

I'm requesting changes for this loading failure. I haven't rerun native Wflow, and the response-tail concern already raised remains open.

[Barbhuiya12]

Reviewed at 480c818. ht validate (34 probes, 77 models), ht gate --probe momentum/wave-celerity-bounds and the full pytest -q (1256 passed) all pass here. The ORCID fix clears the registry loading failure @cehw found. The two stale points from the automated pass are also resolved: the 256-cell comment now matches the length-scaled mesh, and the unresolved-timing branch now fails through _ResponseFailure instead of writing inf and -1.0 into the diagnostics.

Margin

I ran both references over seeds 0–19:

model	verdict	residual against c_kin
reference_saint_venant	PASS 20/20	+0.85% to +0.92% on every state
reference_fixed_celerity	FAIL 20/20	−35.1% to +68.9%

The Saint-Venant offset is steady at about +0.9%. That is about what you would expect from evaluating dQ/dA at q_in while the +5% pulse rides slightly faster. So the reference sits at about a fifth of the 5% allowance, and the offset is systematic rather than noise.

The bound does work the ordering check cannot

The fixed-celerity control fails on both magnitude and ordering, so the gate alone doesn't show that the 5% bound is needed. To check that, I replaced the control's dis in the runner with a pure translation at a state-dependent speed. This router is causal, length-dependent and increases with flow, so it passes the ordering check. Over the three gate seeds:

translation speed	verdict	residual
dQ/dA at the Manning normal depth	PASS	−0.2% to +1.1%
mean velocity Q/A at the same depth	FAIL	−39.7% to −40.0%

Routing at the water velocity instead of the wave celerity is a common real mistake, off by roughly 5/3. The probe catches it through the magnitude bound alone. The first row also shows that the criterion recovers a known celerity to about 1% at an hourly step. You could add this router as a second must_fail, but I don't think it's required.

Smaller points, none blocking

1. The tail check reads one sample. tail_excess is the excess at response_stop - 1 only. That fixes the case the automated review described: a K ≈ 40 h reservoir leaves about 16% of peak at 72 h and is now refused. But a response that crosses base at the last step and is still material a few steps earlier gets through. Bounding the truncated share of the response volume would be the direct form of the argument the README makes. Neither reference comes near this, so a follow-up is fine.
2. closure's sum_inflow. On main this denominator already existed and treated q_in as mm/day. Converting it through area_km2 is right, but no probe in the suite, this one included, uses sum_inflow, so nothing here exercises it except test_routing_inflow_contract.py. It's worth one line in the PR description, so the shared-criterion change isn't missed at merge.
3. Wflow version. 480c818 changes WflowSbmAdapter.jl after the 1.0.4-ht.5 archive rows were written. The change only adds an error on the q_in path when pr or pet is non-zero. This probe's generator sets both to zero, and the other probes never enter that path, so every archived row should reproduce. I can't run the native image here, so please confirm that, or bump to ht.6 and regenerate.
4. Proposal link. #148 is accepted and assigned to you, but the PR doesn't close it and the first checklist box is unticked. Adding Closes #148 would tidy that.

The physics is right, and the paired-centroid design is the correct way to cancel the shared upstream response. I also found the LISFLOOD note more convincing for being left as a disagreement rather than tuned away. Approving. Items 1–4 can be handled here or after merge.

[cehw]

Thanks, Jingzhi — the new tail check makes an unfinished response distinguishable from a measured wave speed, which makes the result much clearer.

I reran 480c818 on climet3: ht validate, 93 focused tests and the three-seed wave-celerity gate pass. The ORCID fix clears my earlier loading blocker, and the new tail and strict-JSON regressions pass as well.

I also ran Wflow 1.0.4 through the PR's Julia adapter on all three original gate seeds: all pass, with wave-speed deviations of 0.34–0.93% against the 5% allowance. The new pr/pet rejection checks work, and the ordinary no-q_in path still runs. This was native Julia execution, not a Docker-image check.

I have no further blocking findings. The one-sample tail limitation is already covered in @Barbhuiya12's follow-up.

Approving.

[mimosapudical]

Thanks both for taking another careful look — really appreciate it :)

I’ve wrapped up the last tail-completeness cleanup as well, including the regression for that edge case. That should close out the remaining technical follow-up on my side.

[Barbhuiya12]

Re-reviewed at 29baebe, which is current with main. ht validate (34 probes, 77 models), ht gate --probe momentum/wave-celerity-bounds and the focused tests (29) pass. On seeds 0–19 nothing moved: reference_saint_venant passes 20/20 at +0.85 % to +0.92 %, reference_fixed_celerity fails 20/20, and the Q/A velocity router still fails at −40 %.

The tail check is now the direct form. It takes the share of the captured response volume in the final pulse-length block, rather than one edge sample. I tested it with a router that translates at the exact dQ/dA and then passes the outlet through a linear reservoir of residence time K, on two gate seeds:

K	verdict
5 h	PASS, residuals unchanged (−0.2 % to +1.1 %)
15 h	refused: final 6 h hold 2.5 % of the captured volume
40 h	refused: final 6 h hold 3.6 %

So a slow tail can no longer pass by crossing base on the last row, which was my point 1.

One observation, not a request. A refused response is scored FAIL, not N/A. For K = 15 h the truncation bias would largely cancel between the two reaches, so that model is failed on the window length rather than on its celerity. The message says "not contained", which is honest, and failing is the conservative side of the asymmetric rule, so I'm fine with it as it is. If a submitted model with real floodplain storage lands there, a longer response_hours is the lever.

Two small things from my last review are still open, neither blocking: the PR body still doesn't carry Closes #148, and wflow_sbm stays at 1.0.4-ht.5 across the 480c818 adapter change (outputs unchanged on every path the suite runs, as far as I can see without the native image).

Approving again.

[chrimerss]

Thanks, Jingzhi. The paired-reach design is right, and the first moment recovers dQ/dA for every honest router I tried. I'm asking for one change before merge, in how the centroid is taken. The other points below are small.

What I checked, at 29baebe

• pytest -q: 1257 passed. ht validate: 34 probes, 77 models. Full ht gate: passes. On this probe, reference_saint_venant passes and reference_fixed_celerity trips wave_celerity_bounds only.
• wflow_sbm in Docker at this head, gate seeds: PASS, 0.704 < 0.927 < 1.225 m/s. The row is identical to the archived one apart from the date, so the Wflow evidence survives the 480c818 adapter change and the 29baebe tail check.
• Synthetic routers on the real generated cases, calling the criterion directly (gate seeds plus seeds 1–5):
  • pure translation at c_kin: within 1.2%;
  • translation at the wide-channel 5/3 u: within 2.8%;
  • a single linear reservoir with K = L/c_kin: within 0.4%;
  • a 10-cell cascade: within 0.8%;
  • translation at 0.9 × 5/3 u: fails at −10%;
  • 1.06 × c_kin: fails;
  • a length-independent outlet store with K = 15 h or 40 h: refused by the tail check (0.025 and 0.036 of the captured volume).

The change: take the first moment of the signed excess

_centroid clips the excess at zero before taking the moment (wave_celerity.py:198). The first-moment identity the README rests on holds for the signed response. A router whose outflow first dips below base therefore loses that dip from the moment, and the dip is read as a celerity error.

The router that shows this is textbook Muskingum, run as one subreach at the hourly step with the exact K = L/c_kin. On the 20 km reach, K ≈ 8 h, so C0 = (Δt − 2KX)/(2K(1−X) + Δt) is negative and the outflow dips first. Its signed first moment is exactly K, so its celerity is right:

Muskingum, one subreach, K = L/c_kin	as merged (clipped)	signed first moment
X = 0.2	FAIL on 1 of 8 seeds (+5.0%)	PASS 8/8, ≤ 0.3%
X = 0.3	FAIL 8/8 (+11.5%)	PASS 8/8, ≤ 0.5%
X = 0.45	FAIL 8/8 (+25%)	PASS 8/8, ≤ 0.8%
X = 0.2, 2 km subreaches (no dip)	PASS 8/8	PASS 8/8

The dip is a numerical artifact of a coarse Muskingum step. Whether a response goes negative is mass/response-nonnegativity's question, though, and this probe would report it as a wave-speed error of +11.5%. X = 0.3 is a common default, and a submitted Muskingum model running one subreach per reach at an hourly step would be failed here with the wrong reason in the archive.

The fix is four lines:

-    response[first:response_stop] = np.clip(
-        discharge[first:response_stop] - base, 0.0, None
-    )
+    response[first:response_stop] = discharge[first:response_stop] - base
     total = float(response.sum())
@@
-    tail_sum = float(response[tail_start:response_stop].sum())
-    tail_fraction = tail_sum / max(total, 1.0e-12)
+    magnitude = np.abs(response)
+    tail_sum = float(magnitude[tail_start:response_stop].sum())
+    tail_fraction = tail_sum / max(float(magnitude.sum()), 1.0e-12)

With this applied:

• the 73 focused tests (test_wave_celerity_bounds.py, test_routing_inflow_contract.py, test_report_detail.py) pass;
• the probe gate passes;
• the Saint-Venant residual is unchanged at +0.9%;
• the fixed-celerity control still fails;
• the K = 15 h and 40 h tails are still refused at the same fractions;
• all three Muskingum rows above pass.

The existing total <= 0 guard still catches a response with no net positive volume. Please add a regression test with a Muskingum (or any dipping) router so the clip can't come back. With X = 0.3 and one subreach, it fails at 29baebe and passes after the change.

Smaller points, none blocking

1. The nine N/A rows for this probe are stale. Rerunning cwatm, dhbv2, flex_lumped, flex_topo, google_flood_forecast, lisflood, modflow6, sacsma_snow17 and summa at this head gives the same verdict and reason for each, but a different detail. The rows were written on 09-25, before the probe required q_in. For example, modflow6's row says the probe needs forcing pr, and none of them mention q_in. Since the criterion changes anyway, please regenerate them with ht run --model <name> --probe momentum/wave-celerity-bounds --gate-seeds --csv models/result.csv, and the wflow_sbm row with them.
2. README count. "Twelve of the thirty-four require no model output beyond runoff" should be thirteen: this probe requires only dis, and routing-lag-consistency is already counted with its geometry inputs.
3. wflow_sbm version. It stays at 1.0.4-ht.5 across the 480c818 adapter change. Outputs are unchanged on every probe the suite runs, and the Docker run above confirms it for this one, so I'm fine keeping the label. Bumping to ht.6 when you regenerate the rows would be cleaner.
4. Closes #148 in the PR body, so the proposal closes on merge.
5. The LISFLOOD paragraph in the probe README reports 12.9–63.1 m/s from an adapter wiring that isn't in the repository. A 40–60× error looks more like the declared reach length not reaching LISFLOOD's channel storage than like LISFLOOD's routing, so the numbers may describe that wiring rather than the model. Please either say the audit wiring is unverified or drop the numbers, until someone can reproduce them.

Process

• Momentum-pool approvals: @Barbhuiya12's at 29baebe is the one current approval. @cehw's approval was on 480c818, was dismissed by the push, and cehw is in the energy pool. Once the change lands, this needs Barbhuiya12's re-approval plus mine or @kawh1111's.
• The probe workflow shows action_required on all three heads, so CI hasn't run yet. A maintainer needs to approve the run on the next push.
• #130 defines q_in as a reach-indexed model output, while this PR defines it as a prescribed forcing. Whichever merges second has to reconcile the name. Nothing to do here, but it's worth knowing.

[chrimerss]

This clip drops the part of the excess that goes below base, but the first-moment identity holds for the signed response. A one-subreach Muskingum with the exact K = L/c_kin and X = 0.3 dips first on the 20 km reach (C0 < 0). Its signed first moment is right to 0.5%, but clipped it reads +11.5% and fails. Suggest response[first:response_stop] = discharge[first:response_stop] - base, with the tail share taken on np.abs(response). Details and the verified patch are in the review body.

[chrimerss]

@mimosapudical please request review once you're ready

[mimosapudical]

@mimosapudical please request review once you're ready

Thanks — the requested changes are now in on the current head, including the signed first-moment fix and its Muskingum regression test. I’ve requested re-review. Appreciate you taking another look.

[mimosapudical]

@Barbhuiya12 @chrimerss thanks again for the detailed feedback. I’ve synced the branch with the latest main, cleaned up the remaining conflicts, and reran the checks. Everything is passing on my side now.

Would really appreciate another look when you get a chance — thanks!