Expand snm coverage for packaged snow models

[jinlong-hu]

Closes #137

Summary

This PR expands snm coverage for packaged models with an explicit snow module, so that mass/snowpack-mass-closure can be evaluated where a model-native snow-module outflow is available.

The mappings added are:

Model	snm mapping
CWatM	Rain + SnowMelt + IceMelt
δHBV2	RAIN + tosoil
LISFLOOD	Rain + SnowMelt
Wflow SBM	native snow.runoff
SUMMA	native scalarRainPlusMelt

The other relevant packaged models were also reviewed. flex_lumped and flex_topo do not expose a separate snow-module control volume with a distinct liquid outflow matching the snm contract. google_flood_forecast has no snow module, while modflow6 represents groundwater and surface-water exchange without snow processes, so they remain unchanged.

Implementation

For each supported model, this PR:

• maps the model-native snow-module outflow to snm;
• updates model.yaml to declare the new flux;
• documents the physical interpretation of the mapping;
• preserves the native model quantity rather than reconstructing snm from a residual;
• updates the archived probe results and model standings.

For δHBV2, snm is RAIN + tosoil. tosoil is read directly from the model, while RAIN is reconstructed using the same temperature-threshold rule as the pinned hydrodl2 version. Each adapter run cross-checks that reconstruction against the model's exposed snow stores using Δ(SNOWPACK + MELTWATER) = pr - snm. The check retains a 1e-4 mm absolute tolerance floor and expands it, when necessary, to four float32 ULPs at the native snow-state and flux scale, accounting for δHBV2's native numerical precision at deep snowpacks while retaining the reconstruction check. On the collaborator-provided 730-day deep-snow reproduction, the maximum step residual is 0.000195183 mm against a precision-aware tolerance of 0.00048828125 mm, and the adapter writes all 730 output rows.

For Wflow SBM, snm is read directly from the native snow.runoff variable. Its snowpack budget closes to machine precision. In the snowpack-closure probe canopy_capacity_mm = 0, mapped to Cmax = 0, so precipitation reaches the snow module without upstream interception.

For SUMMA, snm is mapped to scalarRainPlusMelt. With explicit snow layers this is the basal liquid flux leaving the snowpack; without an explicit snow layer it also carries rain and melt entering the soil. The adapter also reports signed sbl = -scalarSnowSublimation, preserving SUMMA's native sublimation/frost exchange: positive values remove snow by sublimation and negative values represent frost deposition.

Results

After rebasing onto the current main, I reran the adapter checks and the full gate-seed suite.

Model	Standing	Snowpack closure	Spin-up cycle invariance
CWatM	18 / 22	PASS	PASS
δHBV2	12 / 22	PASS	FAIL
LISFLOOD	20 / 22	PASS	PASS
Wflow SBM	19 / 22	PASS	PASS
SUMMA	10 / 25	PASS	PASS

The δHBV2 spin-up failure is the existing mrso state-bound issue: the reported soil moisture exceeds the declared capacity under all three tested spin-up lengths.

SUMMA now passes mass/snowpack-mass-closure, with a cumulative residual of 0.0000% of precipitation. The closure uses independent native quantities: scalarRainPlusMelt as snm, signed -scalarSnowSublimation as sbl, and scalarSWE as snw; it is not reconstructed from a residual.

Reporting signed sbl exposes a semantic difference between its use in the snowpack balance and in the current latent-heat split criterion. The snowpack balance accepts negative sbl as net frost deposition, while the latent-heat criterion treats the reported sublimating share as non-negative and therefore flags SUMMA's negative frost-deposition steps. The affected energy probes remain FAIL, while the snowpack mass balance closes to machine precision.

Validation

All five modified adapters pass ht verify-adapter.

The repository-level checks also pass:

pytest -q
1271 passed, 8 skipped

ht validate
✅ validation passed: 35 probe(s), 80 model(s), suite 0.1.0

ht gate
✅ GATE PASSED: every probe separates the reference models as declared.

The archived rows for these model versions were regenerated against the current 35-probe suite, including mass/groundwater-datum-invariance, mass/spinup-cycle-invariance, and momentum/froude-regime.

[chrimerss]

/review

[github-actions]

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

The PR maps snm (snow-module liquid outflow) for five packaged models so mass/snowpack-mass-closure can score them, bumps each model version, and archives a full 33-probe gate-seed row set per model, updating the top-level README, each model README and the site's three language blocks.

The mappings themselves look sound. CWatM (Rain + SnowMelt + IceMelt), LISFLOOD (Rain + SnowMelt) and Wflow SBM (native snow.runoff) are read from native terms and close the pack budget algebraically; the Julia column re-index at WflowSbmAdapter.jl:714-722 and the LISFLOOD record unpacking are both consistent. The archived standings (12/22, 19/22, 9/25, 18/22, 20/22) match models/result.csv exactly, and every new version has a row for all 33 probes, so tests/test_docs_in_sync.py should pass.

The thing to decide before merging is SUMMA. Its new snm row is a FAIL at 18.85% of precipitation, and the README explains it as a contract limitation: sbl cannot be reported because HydroTuring's sbl must be non-negative. Nothing in AGENTS.md, spec.py or the criteria says that — coherence.py explicitly nets out deposition reported in sbl. So this publishes a snowpack conservation violation against SUMMA that is most likely an adapter omission, not the model.

Two secondary items: the LISFLOOD README now contradicts itself about the two ten-year probes (ERROR vs completed), and δHBV2's reconstructed RAIN term is not exercised by any archived case.

---

models/lisflood/README.md:389 medium: LISFLOOD README still documents ERROR rows the archive no longer has

The rewritten Result section (line 398) says "FAIL (VIOLATION), 20 of 22 probes passed" and models/result.csv:762-794 shows mass/precipitation-counterfactual PASS and mass/human-abstraction VIOLATION, both on 3650-day windows. But lines 389-394 still read "A ten-year daily record takes 97 s there. That is over the 60 s budget of mass/precipitation-counterfactual and mass/human-abstraction … so the archived ERROR rows record no wall time of their own", and lines 434 and 450 still say those probes are ERROR. The "Native re-run" section that explained how the provisional ERROR rows would be replaced was deleted in this PR, and model.yaml's closing comment (which also carried the "Default flood-event window" note) went with it.

Nothing in the adapter changed its speed and both probes still carry max_runtime_s: 60 (probes/mass/precipitation-counterfactual/probe.yaml:35, probes/mass/human-abstraction/probe.yaml:35), so the verdict moved from FAIL (ERROR) to FAIL (VIOLATION) purely because a different host ran it — and the README no longer records which host produced the archived rows.

Failure scenario: a reader checking why LISFLOOD's verdict changed finds the README asserting the two probes cannot finish inside their budget on the host, and no statement of the host these rows came from; the next person to re-run gets different rows and cannot tell whether that is a regression.

Suggested fix: rewrite lines 378-394 for the host that produced the .6 rows (or say plainly that the archive is now native x86-64 and the emulated timings are historical), and restore the window sentence to model.yaml.

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

[cehw]

Thanks, Jinlong — exposing the native snow-module outflow makes the existing snow-balance probe usable for more models. δHBV2’s mapping matched the native model in my cold-snow and rain-on-snow runs.

One issue to fix before merging: the new 1e-4 mm cross-check at 21a8308 stops δHBV2 on a 730-day synthetic seasonal case with peak snow water equivalent of 1,456 mm. The residual is 0.000199 mm, comparable to native float32 roundoff. With the same forcing and published ep100 weights, the previous adapter writes all 730 rows; the new one writes none.

Could the check account for native state precision while retaining the rain-reconstruction check? I ran this on Linux with the pinned model dependencies, using Python 3.10 outside Docker, and have included the reproduction input below. That would keep this added coverage usable for deeper snowpacks too.

Reproduction input and environment

Native CPU inference on Linux: Python 3.10.19, torch 2.5.1+cpu, numpy 1.26.4, pandas 2.3.3, dmg 1.4.3, hydrodl2 1.3.5, dhbv2 0.5.4; DHBV_THREADS=2 and OMP_NUM_THREADS=2. I used the published hbv_2_ep100.pt, config.yaml and normalization_statistics.json from the same official release bundle named in the Dockerfile. This was not a Docker-image test.

The input generator below is extracted from the diagnostic I ran. Its regenerated forcing and static parameters matched the recorded experiment exactly. It produces two synthetic seasons, with 180 cold days, 120 dry melt days and 65 warm days per season. Total precipitation is 2,973.091151 mm, with a maximum daily value of 81.878663 mm. No explicit elevation is supplied, so the adapter uses its training-mean elevation.

Save as reproduce_input.py in a scratch directory and run it with the stated numpy/pandas versions:

import json
from pathlib import Path
import numpy as np
import pandas as pd

STATIC = {
    "area_km2": 250.0,
    "soil_capacity_mm": 320.0,
    "canopy_capacity_mm": 0.0,
    "degree_day_factor_mm_per_C_day": 3.2,
    "baseflow_coefficient": 0.006,
    "snow_threshold_degC": 0.0,
    "latitude_deg": 44.0,
}

def frame(time, pr, tas):
    doy = pd.to_datetime(pd.Series(time)).dt.dayofyear.to_numpy()
    daylength = 1.0 + 0.35 * np.cos(2 * np.pi * (doy - 172) / 365)
    pet = np.maximum(0.0, 0.13 * (tas + 5.0)) * daylength  # the probe generator's PET form
    return pd.DataFrame({"time": list(time), "pr": np.round(pr, 6), "tas": np.round(tas, 6), "pet": np.round(pet, 6)})

def deep_snow_case(scale, years=2, seed=20261002):
    """Single deep winter per year: 180 cold days (all snow), 120 dry melt days, 65 warm days.
    Same draws at every scale, so only the amounts (and the SWE magnitude) change."""
    rng = np.random.default_rng(seed)
    n = 365 * years
    time = pd.date_range("1999-10-01", periods=n, freq="D").strftime("%Y-%m-%d")
    day = np.arange(n) % 365
    cold, melt = day < 180, (day >= 180) & (day < 300)
    noise, innov = np.zeros(n), rng.normal(0.0, 0.8, n)
    for t in range(1, n):
        noise[t] = 0.72 * noise[t - 1] + innov[t]
    tas = np.where(cold, np.minimum(-6 + noise, -3.0), np.where(melt, np.maximum(8 + noise, 4.0), 10 + noise))
    snow_draw, wet_cold = rng.gamma(0.8, 1.0, n), rng.random(n) < 0.5
    rain_draw, wet_late = rng.gamma(0.7, 8.0, n), rng.random(n) < 0.3
    pr = np.where(cold & wet_cold, scale * snow_draw, 0.0)
    pr = np.where(~cold & ~melt & wet_late, rain_draw, pr)
    return frame(time, pr, tas), dict(STATIC)

case = Path(__file__).resolve().parent / "publication_input"
(case / "input").mkdir(parents=True, exist_ok=True)
(case / "output").mkdir(exist_ok=True)
f, static = deep_snow_case(scale=18.0, years=2, seed=20261002)
f.to_csv(case / "input/forcing.csv", index=False)
(case / "input/static.json").write_text(json.dumps(static, indent=2))
request = {
    "case_id": "deep-snow-native-control", "seed": 20261002,
    "timestep": "PT1D", "n_steps": len(f),
    "input": {"forcing": "input/forcing.csv", "static": "input/static.json"},
    "output": {"table": "output/result.csv", "run": "output/run.json"},
}
(case / "request.json").write_text(json.dumps(request, indent=2))
print(case / "request.json")

The paired model runs used the base adapter, version .3 and the PR adapter, version .4, with the same model directory and input. To repeat the comparison, save these as ht_adapter_base.py and ht_adapter_pr.py, set DHBV_MODEL_DIR to the unpacked bundle directory, and create separate case copies before either run:

export DHBV_THREADS=2 OMP_NUM_THREADS=2
cp -R publication_input baseline_case
cp -R publication_input pr_case
python ht_adapter_base.py --request baseline_case/request.json
python ht_adapter_pr.py --request pr_case/request.json

Observed: base exit 0 and 730 output rows; PR exit 1 before writing the output table or run metadata. The error ends with:

RuntimeError: reconstructed HBV rain partition is inconsistent with SNOWPACK + MELTWATER: maximum step residual 0.000198998 mm

At the worst step, float32 spacing is 0.00012207 mm; the full-record signed snow-budget residual is only -1.65×10⁻⁸ of precipitation. The base snow-storage series matches the native states observed in the PR diagnostic exactly.

[jinlong-hu]

@cehw Thanks for the detailed reproduction. I reproduced the deep-snow case and updated the δHBV2 cross-check to account for the native float32 precision while retaining the rain-reconstruction validation.

The check now keeps the existing 1e-4 mm absolute floor and, where necessary, expands the per-step tolerance to four float32 ULPs at the native snow-state/flux scale.

On your 730-day reproduction case, the updated adapter completes successfully with:

• maximum step residual: 0.00019518285989761353 mm
• maximum precision-aware tolerance: 0.00048828125 mm

The full output contains all 730 rows.

I also re-ran:

• ht verify-adapter --model dhbv2 — PASS
• ht run --model dhbv2 --probe mass/snowpack-mass-closure --gate-seeds — PASS
• the full δHBV2 suite — unchanged at 12/22, with snowpack closure PASS
• repository tests — 1259 passed, 8 skipped
• ht validate — 34 probes / 77 models
• ht gate — PASS

The branch has also been rebased onto the current main, and the five modified model archives have been regenerated on the current 34-probe suite.

[cehw]

Thanks for the fix, Jinlong — the precision-aware tolerance keeps this check useful for deeper snowpacks. I re-ran 29e53da on Linux with the same official weights: both cases complete, and deliberately omitting rain from snm still triggers the check before output.

This resolves my δHBV2 finding. My checks cover that adapter’s native inference outside Docker.