arXiv CSopen access

Water-network decisions share one hydraulic gradient, and it can now be computed exactly

· arxiv_cs

arXiv CS · Papers · License: Open Access

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

Water-network decisions share one hydraulic gradient, and it can now be computed exactly Tianwei Mu1,2,3 , Yue Wang3 , Mingzhe Yuan1,4,∗ , Wenhong Wang1 , Qing Luo2 , Min Xiao2 , Jun Li3 , Hui Yang3 , Manhong Huang5 1

Guangzhou Institute of Industrial Intelligence, Guangzhou 510000, China.

Key Laboratory of Ecological

arXiv:2609.06323v1 [cs.DC] 6 Sep 2026

Restoration of Regional Contaminated Environment, Ministry of Education, College of Environment, Shenyang University, Shenyang 110044, China. 3 School of Municipal Engineering and Environment, Shenyang Jianzhu University, Shenyang 110168, China. 4 Shenyang Institute of Automation, Chinese Academy of Sciences, Shenyang 110169, China. 5 College of Environmental Science and Engineering, Donghua University, Shanghai 201620, China. ∗

Corresponding author: [email protected].

Abstract Calibration, leak localisation and sensor placement on water distribution networks (WDNs) are decisions about continuous parameters, yet the hydraulic engine that defines the physics returns a solution and no derivatives, so practice falls back on derivative-free search or on surrogates whose error the answer inherits. We make the global gradient algorithm itself exactly differentiable: the forward pass reproduces the reference engine’s discrete devices, status switching and low-flow linearisation included, and the backward pass solves the implicit adjoint by reusing the forward pass’s terminal factorisation, so one extra sparse solve returns every parameter’s gradient at once, batched over scenarios on one graphics processor. Across 52 public, synthetic and operational networks and 8,140 simulation frames, every network meets the acceptance criterion, the largest head deviation from EPANET 2.2 is 1.137 × 10−13 ft and 25 agree exactly. One adjoint solve replaces the 906 simulations a finite-difference roughness Jacobian costs on the 905-pipe L-TOWN benchmark, and a leak-inversion training loop runs at 463–470 ms per optimiser step for 256 scenarios, 191 times the prior pipeline. Gradient calibration reaches its endpoint within a median 595 model calls, where the strongest of five tuned metaheuristics needs 8,060 to match it on the training loss and two never do within 20,000. On a 554-link operating network, one adjoint pass audits, pipe by pipe, which roughness parameters the installed sensors can constrain and which sensors to add, on the model the utility already operates. Keywords: Water distribution network; Hydraulic model; Automatic differentiation; Model calibration; Leak diagnosability; Sensor placement

Introduction

Most quantitative work on a water distribution network (WDN) is a search over continuous parameters: which roughness values reproduce the observed heads, which junction is leaking and by how much, where the next pressure sensor should go. Each runs against the same physics, the global gradient algorithm (GGA) of Todini and Pilati (1988), Todini and Rossman (2013), whose reference implementation EPANET 2.2 (Rossman et al., 2020) the field treats as ground truth; each is gradient-shaped, yet each is solved without gradients, because the engine returns a solution and nothing else. The workaround is derivative-free search (Maier et al., 2014), whose cost grows with parameter count and, with leakage management under regulatory and economic pressure worldwide (Sousa et al., 2026), is paid again at network scale whenever the model changes. Two lines of work point at the missing ingredient. One makes models differentiable by replacing them: learned surrogates support state estimation and calibration (Truong et al., 2024, Kerimov et al., 2023), and gradient-based control of urban drainage became tractable because a neural internal model supplied the derivatives the physics engine could not (Zhang et al., 2026). A surrogate is fast, but its gradient is the surrogate’s and its approximation error enters every downstream decision unbounded. The other keeps the physics: Ulusoy et al. (2022) showed on an operational network that gradient-based optimisation outperforms evolutionary search on large continuous design-for-control problems. What that argument needs, and what no published tool supplies, is the exact gradient of the deployed hydraulic engine. A separate acceleration literature makes the forward solve faster, by topological reduction, domain decomposition and parallel computing (Guo et al., 2024) or by replacing the linear solver. All of it optimises one network for one demand frame; this work optimises over a batch of operating conditions and requires the result to be differentiable, so the two are complementary (Section 4). We therefore take the reference numerical scheme exactly as it ships and make it differentiable. Analytic sensitivities of the idealised steady state are classical (Piller et al., 2017) and we claim no priority over them, but a solver in production switches link status, clamps options at parse time and linearises the head-loss law below a flow threshold, and the object differentiated here includes those devices. Three contributions follow. (i) The gradient object exists and can be trusted. The backward pass is the implicit adjoint of the very iteration the utility runs, discrete devices included, not of a smoothed stand-in, and the forward pass it hangs on is verified against the reference binary across 52 networks and 8,140 frames, largest head deviation 1.137 × 10−13 ft, none excluded (Section 3.1), so transfer to the utility’s own model carries no surrogate residual. (ii) The object is affordable at decision scale. Batching the whole solver, status machine included, on one graphics processor and reusing each scenario’s terminal factorisation for the adjoint brings one forward-plus-backward L-TOWN scenario from 333 ms to 1.90 ms at batch size 256; per single solve it is not a faster simulator, and Section 3.2 keeps the qualifications in view. (iii) Three decisions, one object. Calibration is least squares on the adjoint Jacobian, the success or failure of leak search is the coherence of its columns, and sensor placement is its Fisher information, so coherence can be read before a search and the Cramér–Rao lower bound (CRLB) prices the identifiability added sensors buy. The chain runs end to end on public benchmarks and on a utility’s operating network, with leak candidates drawn from its own work orders; placement keeps the sensors the utility has and prescribes, in simulation, additions that restore parameter identifiability but do not recover the leaks noise defeated, while a second 2

objective, against coherence itself, recovers one of the three as coherence predicts (Section 3.6).

Materials and methods

2.1

Networks and data

The verification suite has 52 networks: 23 public benchmark models, 21 fetched from their upstream sources and checked against recorded SHA-256 digests plus two from the Kentucky dataset release; 3 EPANET distribution examples; 23 randomly generated networks released with the code; and 3 operational models: City D (542 nodes, 554 links, 79 throttle-control valves), its emitter variant, and City H (921 nodes, 1,038 links, 6 pumps). Figure S1 draws nine of them. The mainline public network is L-TOWN, the BattLeDIM benchmark (Vrachimis et al., 2020, 2022, CC BY 4.0, Zenodo record 4017659), a cleaned copy differing from the published file by one removed default-pattern line and a line-ending change, both digests recorded. Public calibration studies use a synthetic 25-frame diurnal profile; City D the utility’s own 24-hour pattern. Both operational models are released with the permission of the operating utility, in anonymised form: coordinates carry a rigid transform that removes georeferencing, and identifiers are replaced. In the text they are reported by counts and positional labels, and the three injected leaks are labelled L1 –L3 .

2.2

The global gradient algorithm, and a bit-faithful replica

For a network with unknown junction heads H and link flows Q, one GGA iteration solves the symmetric positive definite Schur system AH(k+1) = F(k) with A = A21 D−1 A12 , where D = diag(dhL,k /dQk ) collects the head-loss derivatives, then corrects every flow by an independent scalar update satisfying nodal continuity exactly at every iterate. EPANET 2.2 adds the engineering devices a replica must reproduce: guard constants (108 for closed links and active pressure-reducing-valve (PRV) penalties, a 10−7 floor linearising the head-loss law at vanishing flow), a discrete status machine for check valves, pumps, PRVs and tanks with hysteresis and throttled checking, and convergence requiring status-consistency, not only a small flow change (Elhay and Simpson, 2011). The implementation is built on an automatic-differentiation tensor library (Paszke et al., 2019) over standard numerical Python, and carries two paths: a replica path for bit-level comparison against the reference engine, and a batched differentiable path (Mu et al., 2026) (Fig. 1). Transcribing formulas is not what makes the replica agree to 10−14 ft; four things invisible in the mathematics are. The symbolic factorisation’s elimination order and left-looking Cholesky loop were ported verbatim (George and Liu, 1981), since any other ordering changes the sequence of floating-point additions; the power and logarithm routines are taken from the reference binary’s own C runtime, which is why bit-level claims are platformspecific; literal constants are reproduced, down to that build’s truncated fallback π; and so are parse-time clamps, down to its silent clamping of the engine’s convergence-tolerance option. Backward error analysis says nothing finer than the arithmetic order is verifiable: with κ2 (A) = 1.64 × 1011 on L-TOWN with all three PRVs active, two implementations differing only in reduction order may disagree at the 10−4 ft level, and the same batched solve on CPU and GPU differs by 1.023 × 10−5 ft. A re-implementation agreeing at 10−6 ft is therefore at its noise floor (Supplementary Note 1). 3

a Network and parameters H0

b Batched forward (GGA)

pump

c Implicit adjoint

assemble A(Q; s), F

tank

∂L=∂µ = ¡¸ > ∂F=∂µ

[B, nnz] CSR values; pattern fixed

reservoir

factorise + solve A H = F

not converged

stateful GPU solver (cuDSS) active PRV: constraint row

PRV

leak (emitter): rank-one diag + RHS

adjoint solve A ¸ = ∂L=∂H

reuse

demand

reuses the terminal factorisation zero factorisations in backward

flow update Q mass balance by construction Woodbury row fix at active PRVs

∂L=∂H

status machine (CV, PRV)

naive reuse: O(1) error

active PRV: diagonal penalty row

all scenarios converged and status-consistent?

H¤

loss L(H ¤ )

∂L=∂H0 ∂L=∂r ∂L=∂d ∂L=∂Ke

∂L=∂µ

Figure 1: The differentiable solver. a, Parameter classes: reservoir head, pipe resistance, demand, emitter coefficient; a leak is a rank-one diagonal update. b, Batched forward pass with the full status machine. c, The implicit adjoint, reusing the terminal factorisation, returning all gradient classes at once.

2.3

Exact gradients: the implicit adjoint

Two reverse-mode routes are provided. A truncated K-step unrolling differentiates the solver itself, the route for learning inside the iteration; it is checked against a looser threshold and ships with a health check, because its truncated gradient does not converge in K on every network. The workhorse is the implicit-function adjoint: writing the converged state as the root of the fixed-point residual assembled from EPANET’s own update rules, clamped head-loss branches included, block elimination returns exactly the Schur matrix A of Section 2.2 augmented by the emitter diagonal, so the adjoint system has the same sparsity, conditioning and factorisation as the forward solve. One backward solve and one contraction return the gradient of a scalar loss for every parameter of four closed-form classes (demand, emitter coefficient, fixed head, pipe resistance) at memory cost independent of the iteration count, the standard implicit construction of differentiable optimisation layers (Amos and Kolter, 2017) applied to a new object. Finite-differencing one loss against all 905 L-TOWN pipe roughnesses costs 906 simulations; one forward-plus-adjoint pair returns the same gradient, and a full sensor-by-pipe Jacobian one adjoint solve per sensor (Section 3.2). Where the forward pass factorised a PRV penalty matrix, naive reuse corrupts the demand gradient at the valve’s downstream node by order one while nothing visibly fails; a Woodbury row replacement over the p replaced rows restores relative error to 1.4 × 10−10 (Supplementary Notes 2 and 3). On clamped branches the derivative is exactly zero because the engine’s own guard is the definition; status switching is combinatorial; the backward pass freezes the configuration the forward pass converged to.

2.4

The leak perturbation operator and the signature dictionary

Model-based leak localisation was posed on the steady-state equations by Pudar and Liggett (1992) and carried into practice by sensitivity-matrix methods (Pérez et al., 2011); what follows changes its price and its diagnosability, not its formulation. EPANET represents a leak three ways: an emitter qE,i = Ci pγi , pressure-driven demand (Wagner et al., 1988), or a known discharge added to nodal demand; none touches 4

an off-diagonal entry of A. Inserting a leak at junction i is A 7→ A + βi ei e⊤i plus a right-hand-side component: a rank-one diagonal update leaving the ordering, fill-in and symbolic factorisation invariant, so a candidate sweep over M nodes needs M solves, not M factorisations (Sherman and Morrison, 1950), and the adjoint returns ∂L/∂Ci for all of them simultaneously. Stacking the normalised sensor responses of each candidate gives the signature dictionary any sparse localisation method implicitly works with; its mutual coherence (Tropp and Gilbert, 2007) decides whether sparse recovery can succeed; Section 3.5 measures it, and grouping candidates into coherent, network-adjacent clusters converts it into a group-sparse district search (Supplementary Note 11). The exponent γ is held fixed and only the coefficient estimated (Supplementary Note 4).

2.5

Batched scenarios on one GPU

The batched port runs EPANET’s two cadences of discrete logic for the whole batch by masks: status branches are all evaluated and combined by one-hot selection, an active PRV becomes a penalty on the diagonal, and the loop ends when every scenario is simultaneously converged and status-consistent. On the networks carried in this paper’s mainline and operating-network results the assembled matrix is bit-wise symmetric and one factorisation serves both passes; on two large benchmarks outside the batched path, parallel links in mixed directions perturb A − A⊤ at the last bit, so the reuse there is symmetric to rounding rather than bit-exactly. The sparse route builds the compressed-row pattern once (a leak, a closed link or an active PRV changes only values) and feeds a stateful GPU sparse direct solver whose symbolic analysis is planned once per batch shape and cached; an undersized cache silently re-plans every step at 4.0–4.5× the cost. The backward pass adds no factorisation, verified by independent counters. Unsupported features raise at construction rather than degrade. Timings come from two RTX 5090 nodes measured independently in double precision, ranges being minimum–maximum; memory is peak device residency in one fresh process per cell (Supplementary Notes 5–7).

2.6

Calibration problem, baselines and fairness protocol

Calibration estimates Hazen–Williams roughness coefficients from noisy junction heads: on Hanoi, 34 free pipes from 25 synthetic-diurnal frames at 31 junctions (25 training sensors, 20 training frames, 5 validation frames), with σ = 0.1 ft Gaussian head noise per seed; on City D, 432 free pipes after freezing 43 structurally unidentifiable ones. The gradient calibrator runs Adam, then L-BFGS, then a Levenberg–Marquardt polish whose Jacobian comes from the adjoint. Five tuned derivative-free baselines, differential evolution (DE), particle swarm (PSO), CMA-ES, simulated annealing (SA) and a DE→Levenberg–Marquardt hybrid, were each tuned over four configurations and three dedicated seeds (≈240,000 model calls each on Hanoi; SA 60,576, reduced because it is serial, and declared; 3,168–3,600 on City D), then evaluated on 30 (Hanoi) or 5 (City D) fresh seeds at matched model-call budgets, with backward passes costed at their measured ≈1 % of a forward call. Baselines received batched forward evaluation, a more converged forward solve, ground-truth-prior initialisation and generous snapshots; effect sizes are Vargha–Delaney A12 with two-sided Wilcoxon tests, cross-checked against brute-force implementations (Supplementary Note 8). Enhanced gradient arms, two preconditioners and batched multi-start, share the same pipeline (Supplementary Note 13). 5

2.7

Sensor placement objective

Placement selects k pressure-sensor junctions to make roughness identifiable. From the adjoint-built sensitivity matrix of every candidate junction head to every pipe parameter, a Bayesian D-optimal objective (log-determinant of the prior-regularised Fisher information, σprior = 15, σnoise = 0.1 ft) is maximised by lazy greedy search (Krause et al., 2008), with an optimality certificate and a submodularity check. An augmentation mode holds an installed set fixed and adds k sensors greedily, for the same objective or for coverage (pipes crossing the census threshold), asserting at every step that no pipe loses identifiability; every added sensor is virtual (Supplementary Note 10).

2.8

Verification protocol

Forward verification compares every network frame by frame against the double-precision EPANET 2.2 dynamic library obtained through WNTR (Klise et al., 2017), requiring max |∆H| < 10−6 ft, max |∆Q| < 10−6 cfs and per-frame equality of time steps, statuses, settings and Newton iteration counts. Every input field class, reservoir heads and valve settings included, is read from the model file’s text rather than through a unit round-trip (Section 3.1). Gradient verification is four-fold: Richardson-extrapolated central differences, an unrolled-vs-implicit cross-check, a 22-item audit of degenerate cases (six on the operating networks; all pass, the worst reaching 0.57 of its threshold), and external checks driving finite differences through the compiled reference library. A 54-item regression suite (53 numerical checks, all passing, including seven on the operational models, plus one packaging guard) covers alignment, replay, extended-period simulation and the symmetry and status-schedule guards. The full protocol is Supplementary Note 9; per-network and per-coordinate results are Tables S1–S5.

Results

Every result below reads one object: the Jacobian of the deployed engine’s converged state with respect to its parameters, delivered whole by the implicit adjoint of Section 2.3.

3.1

Exact gradients of the solver as deployed

Fidelity is the warrant for transfer, and it is scoped: bit-level figures belong to the serial replica path, the batched path that runs the decision chain agrees within the conditioning envelope of Section 2.2, and “no surrogate residual” means nothing beyond the reference engine’s own rounding, not GPU bit-identity. Over the 52-network suite and 8,140 frames, all 52 meet the acceptance criterion; the largest head deviation is 1.137 × 10−13 ft (flow 1.421 × 10−14 cfs), 25 agree exactly, and per-frame Newton iteration counts are identical on all 52 (Table 1, Fig. 2). One input path deserves a caution. Every field class the solver reads can be rebuilt bit for bit from the input text, but a parser that takes reservoir heads and valve settings through a metre representation and back loses one unit in the last place on models written in US customary units. A reservoir head is a boundary condition, and the suite’s largest network (12,527 nodes) amplifies one: perturbing a single roughness value inside the reference library by one unit in the last place moves its own answer by 1.3 × 10−6 to 6.9 ft. Reading those two field classes from the input text is what brings

five networks to zero over every frame, one of them a 609-frame model, and it changes neither a field nor a digit on eight control networks (Table S13, Supplementary Note 1). The 1.137 × 10−13 ft that remains has another cause. Table 1: Forward fidelity against EPANET 2.2 over the 52-network suite by provenance. Deviation columns are maxima over every frame of every network, none excluded, with reservoir heads and valve settings read from the input text (Section 3.1); five rows reach these values only that way (Table S13). Per-network rows: Table S1. Networks

Frames

≤ 10−12 ft

Exact zero

Worst |∆H| (ft)

Public benchmark models

7,964

1.137 × 10−13

EPANET distribution examples

5.684 × 10−14

Randomly generated

2.842 × 10−14

Operational (City D, variant, City H)

1.421 × 10−14

All

8,140

1.137 × 10−13

Provenance

public (23)

synthetic (23)

example (3)

operating (3)

metre round-trip path (5) one unit in the last place on a reservoir head, amplified

max |ΔH| vs reference engine (ft)

100

10−4 acceptance threshold 10−6 ft

10−8 worst remaining: 1.1 × 10−13 ft 10−12

exact zero: 25 networks 0 (exact) 101

102

103

104

network size (nodes)

Figure 2: Forward fidelity: per-network maximum head deviation from the reference engine against network size, all 52 networks. Twenty-five agree exactly (floor row) and the worst is 1.1 × 10−13 ft; open circles mark the five networks that reach the floor only when reservoir heads and valve settings are read from the input text (Table S13).

The backward pass meets the same standard. Against Richardson-extrapolated central differences the implicit adjoint reaches a worst relative error of 5.9 × 10−8 over 61 coordinates and four parameter classes per network, the operating network included, against a 10−6 threshold. Externally, with finite differences driven through the compiled EPANET library, eight coordinates on the public Hanoi model with five emitters agree to a worst 5.56 × 10−5 against a 10−4 threshold, and 24 non-clamped City D pipes spanning 1.7 decades of sensitivity to 4.78 × 10−5 ; a second check on the emitter variant covers four parameter classes to 2.39 × 10−5 . The six clamped pipes are proven signal-free: the analytic gradient is exactly zero while the finite-difference residual, at most 6.45 × 10−7 , matches the independently predicted noise floor and sits a factor ≥ 1.6 × 104 below the smallest non-clamped gradient.

3.2

What batching buys

On L-TOWN (Fig. 3, Table S6), one forward-plus-backward scenario falls from 333 ms to 1.90 ms at B = 256 (175–176×), peak memory at B = 1,024 is 22.6× below the dense path’s, and the sparse-overdense crossover sits near 300 junctions. In a real inversion loop the steady state costs 463–470 ms per optimiser step for 256 scenarios, 190.7–191.6× the prior pipeline. A whole City D sensitivity matrix is

forward, ms per scenario

dense assembly sparse CSR + cuDSS

102

101

100 21

batch size B (scenarios)

forward + backward, ms per scenario

prior serial CPU adjoint GPU adjoint, dense GPU adjoint, sparse

103

102

out of memory in the timing job

101

100 21

peak device memory, forward + backward (MiB)

20.0–21.4× faster than sequential backward passes, 1.22–1.32× end to end. c dense sparse CSR + cuDSS

105

32 GiB card

104

103

102

batch size B (scenarios)

a sparse gain 1.27× at B=1 to 6.42× at B=1024 b at B=256 the factor splits in two: F1=50× adjoint onto the GPU, F2=3.5× sparse vs dense, product 175× vs a serial baseline c measured one fresh process per point: the dense forward first fails at B=1280; sparse was scanned to B=16384 (18.5 GiB) without reaching a boundary

Figure 3: Cost on the public 782-junction L-TOWN benchmark; bands are minimum–maximum over two GPU nodes, double precision. a, Forward wall time per scenario against batch size. b, Forward-plus-backward time, the end-to-end factor a product of two measured factors against a serial CPU-adjoint baseline. c, Peak device memory.

Four qualifications travel with these factors. First, the end-to-end factor is two things never to be quoted as one: the GPU adjoint contributes F1 = 2.8× at B = 1 rising to 57× at B = 512, sparse-versusdense linear algebra F2 = 3.1–3.8× for B ≥ 8. Second, the baseline is our own prior serial CPU-adjoint pipeline, not EPANET, which supplies no gradient at any price; a perfectly eight-way-parallel CPU adjoint would divide these factors by up to eight, arithmetically rather than by measurement. Third, the win is not universal: below the crossover the sparse route is slower than the dense one, down to 0.06× at B = 1,024 on a nine-junction network. Fourth, per single solve the dense default path is slower than the reference engine, by 26 to 73× across a 14-network public size sweep (Table S7).

3.3

Calibration: least squares on the gradient object, at matched budget

Calibration under uncertainty is a live problem here (Kerimov et al., 2023, Du et al., 2026). On Hanoi at σ = 0.1 ft over 30 seeds, the gradient calibrator reaches its endpoint in a median 595 model calls. At every budget from 200 to 5,000 calls all five tuned baselines are worse on every metric, the training-loss separation total (Vargha–Delaney A12 = 0.00, Wilcoxon p = 1.86×10−9 ; Fig. 4). Read as run length, the DE→LM hybrid reaches within 5 % of the gradient endpoint on 96.7 % of seeds at a median 8,060 calls; CMA-ES on 30 %; DE and PSO never within 20,000; and SA (serial, with a reduced, declared budget) not within the 5,000 calls it was run to. At 10,000 and 20,000 calls the hybrid closes the gap: the effect size collapses to negligible (A12 = 0.46 and 0.47) while the difference stays detectable (p = 3.05 × 10−5 and 8

0.031), the gradient endpoint still marginally ahead. Its refinement stage uses our own adjoint Jacobian, so the strongest baseline is itself a consumer of the instrument under test, and the equal-budget claim is bounded: it holds for 200–5,000 calls and should not be quoted beyond its range. DE → LM hybrid

CMA-ES

PSO

b 100 seeds within 5% of gradient endpoint loss (%)

training loss (median)

102 101 100 10−1 10−2 gradient endpoint, median 595 calls SA run to 5,000 only (serial; reduced, declared budget) 3

10 model-call budget

gradient

80 60 40 20 0

103 model-call budget

104

Figure 4: Calibration at matched model-call budget (Hanoi, σ = 0.1 ft, 30 seeds). a, Median training loss against budget for five tuned baselines; the star marks the gradient endpoint (median 595 calls). b, Fraction of seeds reaching within 5 % of the gradient endpoint; SA is right-censored at its declared 5,000-call ceiling.

On City D, on the utility’s own diurnal pattern, 432 free pipes recover with median final training error 8.4 × 10−4 , 9.4 × 10−3 and 8.5 × 10−2 under head noise σ ∈ {0.03, 0.1, 0.3} ft at five seeds each (Table S8), one full 432-pipe calibration costing 220–231 s of wall clock. With truth demands perturbed by ±15 % and the inversion run on nominal demands the median degrades only from 9.11 × 10−3 to 9.51 × 10−3 , though one of four repeats reaches 1.62 × 10−2 ; 5 % sensor bias yields a median 9.86 × 10−3 . Eight Latin-hypercube starts land within a 1.0087× spread of final training error, so the endpoint is not a local-minimum artefact. In the one configuration where a baseline wins, the DE→LM hybrid at a budget of 2,000 calls edges the gradient endpoint (0.0094 vs 0.00941, A12 = 0.60, p = 0.0625, the attainable floor at n = 5), reported as a baseline win. Two controls bound that comparison. Tuning budget: the City D baselines are tuned at 2,000 calls, equal to their evaluation budget, over four configurations and three seeds, and each of the four algorithms selects the same configuration as under the 3,168–3,600-call tuning of Section 2.6; differential evolution at ten times the budget still ends at 1.4 × 10−2 . Wall clock: 631 s for the gradient configuration’s 200 calls against 4,756 s for the hybrid’s 1,997 (one seed); divided by the eight-way parallelism the fairness protocol grants the baseline, the hybrid’s figure becomes 595 s, a tie, so the equal-budget claim rests on the call count. Enhancements to the gradient search do not extend it. A Schur-complement diagonal preconditioner improves the roughness error on the design’s own subspace, 17.71 to 12.40, on 10 of 10 seeds (p = 0.0020), but on the layout-independent ruler of Section 3.6 the same comparison is 7 of 10 (p = 0.34); an exact Gauss–Newton diagonal is worse than the cheap approximation, and batching B starts buys at most 2.8× per model call, never B×, at a cost in accuracy. No arm beats the 200-call gradient configuration, which leads on both rulers (Supplementary Note 13, Table S16).

3.4

Leak search with the same object’s columns

Because a leak is a rank-one update (Section 2.4), each candidate contributes one column to the same Jacobian and one adjoint prices them all at once. The test is candidates a crew would dig: 49 junctions distilled from a year of the utility’s work orders, 40 pressure sensors and three injected leaks of 3.0, 1.5 and 2.2 L s−1 on City D. The noiseless inversion recovers the exact three-leak support with a maximum relative discharge error of 5.9 × 10−7 and residual leakage exactly zero on all 46 other candidates. That recovery is an inverse crime by construction, the leaks being injected into the same model the inversion searches, and it does not survive noise: under σ = 0.1 ft the same inversion recovers none of the three, returning two candidates that are neither of them (Fig. 5a). a

leak discharge (L s−1)

noisy run: all three at zero; the two candidates returned are neither

two denominators, two medians: all 1,770 pairs: 0.826 (dashed) 1,495 within-zone: 0.871 (dotted)

candidate pairs per bin

injected recovered, noiseless recovered, σ=0.1 ft

nearest rival of each true leak: 0.987–0.99997

300 250

cross-zone: 275 pairs, all exactly 0

200 150

within-zone: 1,495 pairs

100 50

0.0

0.2

0.4

0.6

0.8

1.0

signature-dictionary coherence

Figure 5: Success and failure of leak search, and the mechanism. a, Work-order search on the operating network: injected versus recovered discharge for L1 –L3 , exact without noise and empty with it. b, Signature-dictionary coherence on the L-TOWN configuration, 1,770 candidate pairs, both denominators marked; arrows mark each true leak’s nearest rival.

The bracket is not specific to that network. On L-TOWN the test uses three self-injected leaks rather than the BattLeDIM competition scenarios, so published localisation results (Daniel et al., 2022) set context without being numerically comparable; a 60-candidate, 33-sensor inversion returns a true leak node as its largest coefficient but not the full support, the other two finishing 19th and 32nd, beneath their coherent neighbours. The failure is not an optimiser pathology but a measurable property of the columns.

3.5

Where the search stops: the coherence of the columns

The signature dictionary’s mutual coherence on this L-TOWN configuration has maximum 0.9999999981 and median 0.8258 over all 1,770 candidate pairs, including 275 made exactly orthogonal by a PRV, or 0.8713 over the 1,495 within-zone pairs: two denominators, both reported (Fig. 5b). Forty-two pairs, all within-zone, exceed 0.999, and each true leak has an inseparable rival at 0.999965, 0.999231 and 0.987023; the spurious recoveries are the high-coherence neighbours of the one node found. Coherence this close to one violates every sufficient condition for sparse recovery, so the inversion stops where the dictionary says it must, and the whole check costs under a minute of single-process CPU time. The same boundary in a crew’s units: the candidates coherence above 0.99 makes indistinguishable 10

from L-TOWN’s two lost leaks carry 984 and 512 m of connected pipe, while the leak the inversion returns has no rival above that threshold and 96 m at stake. The work-order dictionary (49 candidates, 40 sensors) reads the same way, with a median coherence of 0.934, 51 pairs above 0.999 and a rival at 0.99999 for the 1.5 L s−1 leak: both failures are diagnosed post hoc by one quantity. Ambiguity of that size is a unit of work. Grouping candidates mutually coherent above a threshold τ and adjacent in the network’s Voronoi partition, then replacing the sparse penalty by its group form, searches districts, not junctions (Supplementary Note 11 and Table S14). The trade-off, not the hit rate, is the result. On City D the single-point inversion ranks the three injected leaks 49th, 15th and 16th of 49; at τ = 0.999 two of the three fall inside the top three clusters, which hold seven candidates and 7.42 km of main, 8.8 % of the network, the largest 750 m in radius. On L-TOWN τ = 0.995 returns all three inside 22 clusters averaging 98 m. Neither figure means anything without the radius: at τ = 0 City D is one 9,043 m cluster with a hit rate of 3 of 3, the hit rate of naming the whole network. Thresholds are chosen after seeing the tabulated sweep; against 20 size-matched random groupings that keep the cluster sizes and destroy topology and coherence the design is never significant (p ≥ 0.095) and wins only on radius, 179 against 1,700 m. The leak the amplitude limit defeats stays defeated: signal-to-noise 0.50 against 166 and 151, a singleton at every τ ≥ 0.5, and 49th to 36th even when forced in with its twelve most coherent neighbours. Coherence is a property of the sensor set as much as of the network, which makes it a design variable rather than a fixed limit.

3.6

Sensor placement: a prescription, and how far it closes the loop

Pressure-sensor placement is a standing topic here (Zhou et al., 2024, Cheng et al., 2024) and informationbased sampling design is older still (Kapelan et al., 2005); what differs is the input, exact adjoint sensitivities of the deployed engine rather than a graph proxy, and the certificate. (No graph-spectral baseline is included, since only a proxy of the published method would be available for comparison.) On City D, a gradient-attribution census (one adjoint pass over 25 frames) partitions all 554 links into 279 whose roughness the designated 40-sensor data can identify, 149 it cannot, 41 dead branches, 6 clamped at the low-flow threshold and 79 non-pipe links (Fig. 6a); the zero classes sum exactly to the 196 all-zero sensitivity columns, on each of which central differences also return zero. The census precedes any calibration and reveals, rather than causes, an unidentifiability the sensor configuration already determines. On Hanoi the discriminating number is the identifiability bound: at k = 10 the D-optimal design attains an identifiable-subspace CRLB trace of 1.4746 against a random median of 13.642 (Table 2, Fig. 6b), with a worst-case optimality-certificate ratio of 0.763, zero submodularity violations in 200 samples and an 8-ms search. Hanoi has no structurally unobservable pipe, so “zero unobservable” claims would be vacuous there. On City D the question is which sensors to add to the 40 the utility has, not where 40 would go if none existed. Reselecting from scratch answers the latter: a k = 40 D-optimal layout over 541 candidate junctions recovers 52 of the 149 unobservable pipes against a random median of 31.5 but loses 38 the designated sensors could see, so the unobservable total falls only to 135, and no from-scratch layout on the k ≤ 80 grid reaches 80 % recovery (Fig. 6c). The prescription answers the former: the 40 stay and k are added greedily, for the D-optimal objective or for coverage under the census’s threshold, with 11

a clamped 6 identifiable 279

unobservable 149 dead branch 41

100

200

300

non-pipe 79

400

500

links (554 total), one adjoint pass

100

coverage D-optimal reselect 40 + k

9.25× tighter

recovered of the same 149 (virtual additions, none lost)

101

−1

d D-optimal random (median)

random, median of 30 random, interquartile D-optimal

recovered of the census's 149 unobservable pipes

CRLB trace, identifiable subspace

160 80% target: never reached on grid

140 120

k=40: 52; random needs k ≈ 78

100 80 60 40 20

140 120

all 149 80%

100

+32

+56

80 60 40 reselection loses 15 to 38 pipes

0 2

160

0 0

sensors k (Hanoi)

sensors k, reselected from scratch

sensors k added to the 40 designated

Figure 6: Placement as a prescription, by counts only. a, Gradient-attribution census of the operating network’s 554 links. b, Identifiable-subspace CRLB trace against sensor count on Hanoi, D-optimal versus 30 random layouts. c, Recovery of the 149 unobservable pipes by reselection from scratch. d, The same with the 40 designated sensors kept and k added.

Table 2: Sensor placement for roughness identifiability. Hanoi: identifiable-subspace CRLB trace, Bayesian D-optimal greedy against the median of 30 random layouts. Operating network: of the 149 pipes unobservable under the 40 designated sensors, how many a layout recovers and loses. Augmentation keeps the 40 and adds k; every design is virtual. Hanoi (34 pipes, 31 candidate sites, 1 frame): CRLB trace k

D-optimal

0.034

0.153

0.454

0.832

1.475

Random (median)

0.077

0.696

2.509

5.164

13.642

Operating network (475 pipes, 541 candidate sites, 25 frames): of the 149 unobservable k sensors added to the 40

Augment, coverage objective: recovered

133

149

Augment, D-optimal: recovered

Augment, either objective: lost

Reselect 40 + k from scratch: recovered / lost Still unobservable (cover. / D-opt. / resel.)

59/38

68/24

73/21

94/15

112/128/128

87/126/128

54/114/105

16/97/97

0/73/70

Reselect 40 from scratch (k = 0): 52 recovered, 38 lost, 135 still unobservable; random median 31.5 recovered

the identifiable set asserted never to shrink (zero violations in 80 steps). Coverage restores 37 of the 149 at k = 5 and all 149 at k = 80, D-optimal 21 and 76, neither losing a pipe at any budget (Table 2, Fig. 6d, Table S9). On L-TOWN, adding to the 33 sensors restores all 62 pipes they could not inform at 36 coverage additions, where reselecting 33 from scratch loses 99 (Supplementary Note 10). The downstream gain is graded on a ruler that does not move with the design. Roughness error over a layout’s own informative pipes is not comparable across layouts, that set growing with the sensors added from 279 pipes to 412, so designs are scored on the roughness error projected onto a reference subspace fixed before any design was chosen: the 98 directions a fully instrumented candidate pool would identify, against a prior of 28.92 (Supplementary Note 12). Head misfit ranks nothing here: at 0.1 ft it is 96 % noise variance, and over eight restarts on identical data it moves 0.4 % where this roughness error moves 25 %. On that ruler coverage +20 scores 17.45 and D-optimal +20 18.04 against a random median of 21.73, no draw among 100 random additions from the same pool at the same budget matching either (p = 0.0099, the floor at R = 100); all twelve budget-by-noise cells agree, five only to the coarser floor R = 20 allows (Tables S10 and S15). Two costs travel with it: in the 334 discarded directions the designs do not beat random (60 of 100 at least as good), and in 22 of 48 cell-design-metric combinations the design leads every random layout on roughness while trailing most on heads, so choosing a layout by head misfit chooses backwards. Repeating the work-order search with the augmented sets, same injected triple, noise and inverse-crime setting, does not recover the leaks 40 sensors lost: on L-TOWN all ten augmented sets return the same one of three, and on City D every set returns L2 with 1–4 % discharge error, a hit a control attributes to an added sensor on L2 ’s own junction. L1 moves no head by more than 0.09 ft, below the noise; L3 is recovered by three of five random +20 sets and none of the six designed (Table S11). Placement for roughness identifiability does not close the leak-search loop. A second objective, additions maximising the same dictionary’s coherence, does better: pooling every designed City D set against 40 random draws, L2 recovers 12/14 against 3/40 (p = 1.05 × 10−7 ), and what sorts recovery is each set’s own rival-coherence gap, which separates every recovering configuration from every losing one across all 55 sets whatever objective chose it. L3 ’s recovery is indistinguishable from random (p = 0.51) and L1 by neither, its column norm an order of magnitude below L2 ’s: an amplitude limit, not a coherence one. Both conclusions hold at the coarser stage the L-TOWN inversion shares, where the same mechanism bounds the loop: of 746 positions only one and 14 separate its two unrecovered leaks from their nearest rival (Table S12).

Discussion

Calibration, leak search and sensor siting were three separate black-box searches because the engine defining them exposed no derivatives. With derivatives exposed they are three readings of one object. The same coherence that forbids a junction-level answer sizes the district-level one that survives it: on the operating network the ambiguity a crew must walk is 7.42 km of main, not 84. The placement reading is a prescription, not a redesign: the utility’s sensors stay, the adjoint names the additions, and the census’s unidentifiable pipes become identifiable, though only additions against coherence itself recover any leak noise defeated. Because the forward pass reproduces the reference engine to 1.137 × 10−13 ft, every link

transfers to the model a utility already operates with no surrogate residual, and the census and coherence reading arrive as by-products of one pass, not assumptions. The forward-acceleration literature is orthogonal and the borrowing bounded. The batch path is limited by per-matrix host dispatch, not arithmetic; removing that floor is worth an estimated 4–9×, an estimate rather than a measurement. Three incompatibilities bound what the decomposition line (Guo et al., 2024) can lend: it multiplies the per-matrix bottleneck count, early stopping breaks the implicit adjoint’s premise, and any partition changes the elimination order the bit-faithful path reproduces. Adoption asks a utility for a hydraulic model inside the documented feature subset, the pressure records it logs and one CUDA GPU. The gradients are those of steady-state frames: extended-period simulation is forward-verified, but storage coupling between frames is not differentiated through, so multi-day control synthesis is out of scope. The binding limits bear repeating. Bit-level agreement is a claim about one reference binary on one platform; elsewhere it degrades to ≈ 10−6 ft. Gradients are gradients of a frozen status configuration, one-sided where a perturbation would flip a status. The batched differentiable path covers a documented feature subset, refusing what it does not implement. The speed factors carry the four qualifications of Section 3.2, including that the dense default path is slower per solve than the reference engine. The equal-budget calibration claim is bounded at 5,000 calls and rests on the call count: on the wall clock the strongest baseline, given the eight-way parallelism it is promised, draws level. The leak demonstrations are one network and one injected triple each, with the noiseless operating-network recovery an inverse crime whose noisy counterpart fails; identifiability-driven additions did not change that outcome, coherencedriven additions changed it for one of the three, for the reason coherence predicts; what generalises is the diagnosis, not the recovery. The cluster-level reading buys hits with inspection radius, at a threshold chosen after the sweep, and never separates from size-matched random grouping on hit rate alone. Placement gains are established on a projected roughness error under a single noise realisation, and not in the directions that projection discards. The sensor additions are virtual, designed and verified in simulation, none installed. This is research code, not an operational product.

Conclusions

An exactly differentiable global gradient algorithm, verified against EPANET 2.2 over 52 networks (to 1.137 × 10−13 ft at worst, none excluded) and batched on one GPU, hands three standing WDN decisions one gradient object. Calibration reaches its endpoint in a median 595 model calls where the best of five tuned metaheuristics needs 8,060; leak search runs over a utility’s own work-order candidates, and where noise defeats a junction-level answer the same object’s column coherence says so in advance and sizes the district-level answer that survives; sensor siting keeps an operating network’s 40 sensors and adds the few that restore, in simulation, the identifiability of the pipes their data could not inform, without recovering the leaks noise defeated. The same adjoint prices its own limits, so a utility can decide where sensors, crews and trust should go before committing them.

CRediT authorship contribution statement T. Mu: Conceptualization, Methodology, Software, Validation, Investigation, Visualization, Writing – original draft. Y. Wang: Data curation, Investigation, Writing – review & editing. M. Yuan: Conceptualization, Supervision, Funding acquisition, Writing – review & editing. W. Wang: Software, Formal analysis, Writing – review & editing. Q. Luo: Resources, Funding acquisition, Writing – review & editing. M. Xiao: Validation, Formal analysis. J. Li: Investigation, Resources. H. Yang: Supervision, Project administration. M. Huang: Supervision, Writing – review & editing.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability The two operational network models, and the year of leak repair records drawn from one of them, are released with the permission of the operating utility, in anonymised form under CC BY 4.0. The public networks are fetched from their upstream sources and verified against recorded SHA-256 digests rather than redistributed; sources and licences are listed with the code. The randomly generated networks, their generator, the emitter-augmented Hanoi configuration and every measurement report cited here are released with the code (Mu et al., 2026): development is public at https://github.com/ mutianwei521/wdsgpu and the submitted version is archived with a DOI at Zenodo. L-TOWN is available under CC BY 4.0 as Zenodo record 4017659 (Vrachimis et al., 2020). HydroGrad is the Python package dgga, MIT licence; the sparse route requires CUDA and reference cross-checks the EPANET 2.2 library obtained through WNTR.

Acknowledgements This work was funded by the open fund of the Key Laboratory of Ecological Restoration of Regional Contaminated Environment (Shenyang University), Ministry of Education (grant KF-26-11), and by the Guangdong Province Natural Science Foundation General Project (2026A1515011817). HydroGrad re-implements the EPANET 2.2 hydraulic engine from its MIT-licensed C sources; it is not endorsed by the US Environmental Protection Agency or the OpenWaterAnalytics community. WNTR supplied input parsing and the reference library.

Supplementary Information Water-network decisions share one hydraulic gradient, and it can now be computed exactly Tianwei Mu, Yue Wang, Mingzhe Yuan, Wenhong Wang, Qing Luo, Min Xiao, Jun Li, Hui Yang, Manhong Huang

Contents • Supplementary Notes 1–13: the implementation, derivation and protocol detail that Section 2 of the main text compresses: reproducing the reference arithmetic, including the input round-trip that costs one unit in the last place on a boundary condition (Note 1), the two differentiation routes (Note 2), the big-M pressure-reducing-valve penalty and its Woodbury correction (Note 3), the leak perturbation operator and the signature dictionary (Note 4), the batched status machine (Note 5), sparsity-pattern invariance and the stateful sparse solver (Note 6), the cost decomposition (Note 7), the calibration comparison protocol (Note 8), the verification protocol in full (Note 9), the sensor augmentation protocol with calibration and leak search under the augmented sets, control experiments and the coherence-driven placement test (Note 10), cluster-level leak diagnosability and its trade-off curve (Note 11), the layout-independent metric that scores sensor designs (Note 12) and the optimiser enhancements, the baseline tuning budget and the wall-clock reading (Note 13). • Supplementary Discussion: reference solutions that must not be used as training labels; the relation to learned surrogates; and the complete list of limitations. • Figure S1: the drawn topology of nine of the verified networks, ordered by node count, three orders of magnitude of size in one plate. • Tables S1–S16: the per-network fidelity suite (Table S1), truncated-unrolling health (Table S2), the external per-coordinate gradient checks (Table S3), the 22-item degenerate-case audit (Table S4), the 54-item regression suite (Table S5), L-TOWN cost and memory (Table S6), the 16-network size sweep (Table S7), the calibration noise ladders (Table S8), the sensor-augmentation account (Table S9), calibration with the augmented sets (Table S10), leak search with the augmented sets and its controls (Table S11), the coherence-driven placement test against random controls (Table S12), the exact-input control and its negative control (Table S13), the cluster-level trade-off curve (Table S14), the layout-independent placement metric across every budget-by-noise cell (Table S15) and the optimiser arms with the wall-clock reading (Table S16). Provenance and anonymisation. Every table and every measured value quoted in the prose is generated by a single script from the recorded measurement files it names, and the two operational distribution models are carried at the full 52-network calibre of the main text: 52 networks, 8,140 frames. Throughout, the operational models City D and City H are reported by counts and positional labels only; no model node, link or sensor identifier of an operating network appears in any table or sentence, and rows summarised by sensitivity band rather than listed per coordinate (Table S3b, c) are summarised for exactly that reason. Each model is released with the permission of the operating utility.

Supplementary Notes Supplementary Note 1: reproducing the reference arithmetic Two forward paths, one parser. The implementation offers an EPANET-replica path, which follows the reference control flow statement by statement and carries every bit-level claim in this work, and a dense batched path, which assembles the Schur complement as a dense tensor by an index-scatter accumulation and solves many scenarios at once. They share one parser and one set of coefficient formulas. Transcribing the head-loss laws, the pump curves and the valve coefficients is mechanical; what separates 10−6 ft agreement from 10−14 ft agreement is five things that are invisible in the mathematics. (i) Operation ordering in the linear solve. EPANET performs one symbolic factorisation at load time. It builds a de-duplicated adjacency list (parallel pipes share a single off-diagonal slot and their conductances are added), reorders the junction subgraph by multiple minimum degree (Liu, 1985), simulates the elimination graph to allocate fill-in slots, and stores the strictly lower triangle in compressed form. The numerical phase is √ a left-looking Cholesky with · adapted from George and Liu (1981). Both the elimination order and the loop structure of the column algorithm are ported verbatim, because any other ordering, and any other loop structure over the same sparsity pattern, produces a different sequence of floating-point additions, and at the conditioning quantified below that is not a difference in the last bits of the head. (ii) The transcendental library. The reference engine’s Windows build bundled with WNTR (Klise et al., 2017) is a MinGW-family build whose statically linked power routine agrees bit for bit with the one exported by the Windows C runtime and differs from other implementations, on a small fraction of arguments, by one unit in the last place. Because the Hazen–Williams resistance r = 4.727 L C −1.852 d−4.871 is evaluated once per pipe and then multiplies every downstream quantity, a one-unit difference there is amplified by the conditioning of the assembled system. Every call to the power and logarithm routines is therefore routed through the same C runtime as the reference build; a vectorised power function from an array library is not sufficient, and the residual it leaves grows with the iteration count. This is the reason the bit-level claims are stated for one platform and one binary, and for nothing else. (iii) Literal constants, including the inexact ones. The reference build was compiled without the C runtime’s π constant defined and uses the fallback literal π ≈ 3.141592654 from its own header rather than the nearest double. The initial flow assignment πd2 /4 then differs in its last bits, and the difference is amplified through the first iterations. The literal is reproduced. P (iv) Accumulation fold order. The stopping statistic is the relative flow change, k |∆Qk | divided by P k |Qk |. When its value sits near the convergence tolerance (the engine’s ACCURACY option), a different summation order flips the iteration count, and a different iteration count means the two implementations are no longer answering the same question. The left-fold order of the reference C loop is reproduced exactly. (v) Parsing semantics, including the clamps. EPANET does not necessarily solve the problem the input file describes. The ACCURACY option read from an input file is clamped to [10−5 , 10−1 ] during parsing, while the same option set through the programmatic interface is not, so a file requesting 10−8 is silently given 10−5 . A re-implementation that honours the file literally converges to a different stopping point with every formula correct. The clamp is reproduced.

Why care with the formulas cannot substitute. Let Ĥ be the head produced by a backward-stable double-precision solve. Standard backward error analysis gives ∥Ĥ − H∥ ≲ c(N ) ε κ2 (A), ∥H∥

ε = 2−53 ,

(1)

with c(N ) a modest polynomial in the dimension. The conditioning is set by the guard constants of the main text’s Section 2.2: a closed link contributes a conductance of 10−8 , a fully open valve 106 , and an active pressure-reducing valve (PRV) a penalty of 108 on the diagonal, the engine’s large-coefficient constant (CBIG), so a single assembled matrix spans sixteen orders of magnitude. On L-TOWN with all three PRVs active the converged Schur complement has κ2 (A) = 1.64 × 1011 , and (1) with c(N ) = 1 then bounds the relative disagreement between two implementations differing only in the order of their arithmetic by εκ2 ≈ 1.8 × 10−5 (arithmetic on the two constants just quoted, not a measurement). The directly measured figure sits inside that bound: the same batched solve on CPU and on GPU, where nothing differs but the reduction order of the underlying kernels, disagrees by 1.023 × 10−5 ft. Two consequences follow. Agreement at 10−14 ft is not obtainable by being careful with formulas; it is obtainable only by reproducing the arithmetic. And a re-implementation that agrees at 10−5 –10−6 ft is not defective (it is at the noise floor implied by (1)), but it cannot verify anything finer, and in particular cannot distinguish a modelling error from reduction-order noise. The bit-level target is the only one at which a failed comparison is unambiguously a modelling discrepancy rather than reduction-order noise. (vi) The unit round-trip on the boundary conditions. The five items above concern the arithmetic. A sixth concerns the inputs. Every field class the solver reads can be rebuilt bit for bit from the input file’s own text (diameters, lengths, resistances, minor-loss coefficients, demands, emitter coefficients); reservoir heads and valve settings are the two a parser is most likely to take through an internal metre representation and back, and on a model written in US customary units that round trip costs one unit in the last place. Five networks of the 52 are sensitive to it, four of them above the acceptance threshold and a fifth above 10−12 ft. The sensitivity is a property neither of the reference engine nor of those models; it is a property of the input path. The amplification is the conditioning of (1) applied to a boundary condition. On the suite’s largest network (12,527 nodes), perturbing a single roughness value inside the reference library by one unit in the last place moves the reference library’s own answer by 1.3 × 10−6 to 6.9 ft, and a one-unit error in a fixed head enters the same amplification. Reading those two field classes from the input text, through a non-default parser entry point so that the default numerical path is unchanged to the last bit, drives all five networks to bit-level agreement over every frame (Table S13a): 6.836, 2.30 × 10−5 , 6.086 × 10−6 , 2.146 × 10−6 and 6.08 × 10−12 ft all to zero, with max |∆Q| ≤ 7.11 × 10−15 cfs and with statuses, settings, frame times and Newton iteration counts equal frame by frame. The 609-frame extended-period model among them is covered across all its frames. The negative control establishes the cause. On eight further networks the same input path, reading the same two field classes from the input text and, in a second arm, the elevations as well, changes no field, and every measured maximum is identical (Table S13b). The effect reaches exactly the models a unit round trip can damage, the ones written in US customary units. The 1.14 × 10−13 ft that remains as the suite’s worst deviation is on two metric models the correction does not touch, and its cause is the reduction-order noise of (1), not the inputs. Nothing in this note licenses the stronger claim that all 52 networks agree bit for bit: 25 of them do, against 20 under the default input path, and the rest sit between 10−14 and 1.14 × 10−13 ft.

Supplementary Note 2: the two differentiation routes Route A: truncated unrolling. The dense forward path executes a fixed number K of global-gradientalgorithm (GGA) iterations with out-of-place tensor operations, and ordinary reverse-mode autodifferentiation is applied to the whole graph. What is returned is the exact derivative of the K-step map, at O(K) memory. Because K is finite it is not the derivative of the fixed point, and the gap is measured rather than assumed: this is why the route is checked against a 10−4 threshold where the adjoint is checked at 10−6 , and why the implementation carries a health check for the truncated gradient. Over 7 public networks the truncated demand gradient converges in K on four, is borderline on Pescara and does not converge on 2 (ky4, Net3; Table S2). This is a list of measured networks, not a conditioning criterion, and the health check should be run before the route is trusted on a new one. The route additionally refuses every PRV network, L-TOWN included, at solver construction. Route B: the implicit-function adjoint. Collect the unknowns as z = (Q, qE , H) and write the converged state as the root of F(z; θ) = 0, whose three blocks are the fixed-point form of EPANET’s own update rules, link energy: emitter: continuity:

rk = φk (Qk ) − (Hn1(k) − Hn2(k) ) = 0, 1/γ

rE,i = Ke,i qE,i − (Hi − zi ) = 0, P P mi = k: n2(k)=i Qk − k: n1(k)=i Qk − qE,i − di = 0,

(2) (3) (4)

where φk is exactly the head-loss branch the solver used (including the branch clamped at the engine’s low-flow resistance floor, RQtol) and zi is the node elevation. The Jacobian has the block form   D 0 −A12   J= 0 (5) E′ S , A21 Pe 0 with D = diag(φ′k ) equal to the head-loss derivative the solver’s own Newton update uses (on a clamped branch, the derivative of the clamped expression, so that the Jacobian is consistent with the residual actually driven to zero), E′ the diagonal of emitter gradients and S = Pe = −I. Eliminating ∆Q and ∆qE from (5) returns exactly the Schur complement A of the forward iteration, augmented by the emitter diagonal: the adjoint system has the same structure, the same sparsity pattern and the same conditioning as the forward solve. For a scalar loss L(z) the backward pass solves J⊤ λ = ∂L/∂z once and contracts ∂L/∂θ = −λ⊤ ∂F/∂θ at a memory cost independent of the iteration count. The four closed forms.

Each supported parameter class has ∂F/∂θ available analytically:

• nodal demand di : continuity rows, ∂mi /∂di = −1; • emitter coefficient Ke,i : emitter rows, ∂rE,i /∂Ke,i = |qE,i |1/γ−1 qE,i , and exactly zero on the clamped branch; • fixed head H0,f : link rows through A10 , contributing −1 when the fixed-head node is the link’s upstream end and +1 when it is the downstream end; • pipe resistance rk : link rows, ∂rk /∂r = sgn(Qk )|Qk |n , and exactly zero on the RQtol-clamped branch, on valves and on closed links. Pump speed ω and the pump curve coefficients h0 and R are supported the same way, branch by branch of the pump coefficient routine. For a constant-power pump on the clamped branch, EPANET’s iteration uses +g where the derivative of the residual is −g; that asymmetry is a half-step device serving the Newton iteration, not a Jacobian entry, and the implicit differentiation uses the derivative of the residual. 19

Newton polish, and why it is not optional. In route B the state is first converged with the GGA and then polished with a few full Newton steps on F(z; θ) = 0 using a sparse LU factorisation of the assembled J. The requirement is quantitative: on an ill-conditioned network the GGA’s Schur-complement half-iteration stalls at a relative error near 10−7 , whereas the assembled Newton step drives ∥F∥∞ to ∼ 10−14 ft (on the emitter-augmented Hanoi configuration, 4.263 × 10−14 ft). Without the polish, finite-difference verification measures the residual of the forward solve rather than the error of the gradient. Three sources of non-smoothness, which are not equivalent. Clamped branches (the RQtol floor on gk , the CBIG/RQtol bracketing of constant-power pumps, and the emitter’s one-sided behaviour at Ke = 0) give a piecewise-defined but continuous map: exact inside a branch, one-sided at a boundary, with EPANET’s own guard used as the definition. In particular ∂ ·/∂Ke,i at Ke,i = 0 is exactly zero, matching the engine’s short circuit. Pressure-driven demand would add the three-segment Wagner function (Wagner et al., 1988) with barrier penalties; it is not implemented, so the class does not arise, at the cost of not supporting the model. Status switching is a genuine combinatorial event: a switch changes the effective topology or replaces an equation. The backward pass freezes the configuration the forward pass converged to, so the reported derivative is the derivative of the smooth branch the solution sits on, and with respect to a parameter whose perturbation would flip a status it is one-sided at best. Nothing is smoothed; the 22-item degenerate-case audit (Table S4), 6 of whose cases run on the operating network’s closed-link and throttle-valve neighbourhoods, characterises exactly this regime.

Supplementary Note 3: the big-M PRV penalty matrix, and the Woodbury correction This part of the construction is given in full, because a direct reuse of the forward factorisation at an active PRV returns a wrong gradient with no visible failure. The reduction. Eliminating the link and emitter rows of J⊤ λ = (gQ , gE , gH ) leaves a node-level system whose coefficient matrix is exactly Ā, the converged GGA matrix the forward pass assembled last. Writing P = D̄−1 for the per-branch conductance, B for the junction incidence and h′E for the emitter derivative, the backward pass is Ā λm = B⊤ P gQ − em ⊙ gE /h′E − gH ,

λe = (gE + λm )/h′E , (6) so the whole backward pass is one triangular solve against the forward pass’s terminal factorisation plus sparse matrix–vector products. No new factorisation is required, and independently instrumented counters confirm that none is performed. λl = P gQ − (λm [j2 ] − λm [j1 ]) ,

Where the reuse breaks. At an active PRV the true Jacobian row is a constraint row, Hn2 = hset with Dkk = 0; the valve’s flow column then gives the constraint λm [j2 ]−λm [j1 ] = gQ [prv], and the corresponding λl is a free quantity that enters no parameter gradient (the resistance derivative is zero, there is no fixed-head coupling, and P = 0, so the scatter vanishes identically). But the matrix the forward pass actually factorised is the big-M penalty matrix Â = Ā + CBIG

p X

ej2a e⊤j2a ,

CBIG = 108 ,

(7)

a=1

for p active PRVs. Solving the adjoint system naively with Â−1 drives the adjoint variable at the valve’s downstream row to ≈ rj2 /CBIG, which is numerically zero. On L-TOWN this puts the demand gradient at

that node at −8.9 × 10−8 where the true value is 2.38 × 102 : wrong by order one, with a small residual for the wrong system and nothing visibly failing. The correction. Â and the true node-level matrix M differ by a replacement of the p downstream rows, M = Â +

p X

ej2a ma − âa

⊤

(8)

a=1

where m⊤a is the constraint row and â⊤a the penalty row it replaces: a rank-p update, so the Woodbury identity applies. Writing x0 = Â−1 r for the naive solve and S = Â−1 [ej21 , . . . , ej2p ] for the p correction directions, which cost p triangular solves against the factorisation already in hand, the corrected solve is M−1 r = x0 − S C−1 x0 [j2 ] − x0 [j1 ] − r[j2 ] , Cab = (sb )j2a − (sb )j1a , (9) with C ∈ Rp×p inverted directly. The bracket is exactly the residual of the p valve constraints at the naive solution, so (9) reads as: solve once with the penalty matrix, measure how far the constraints are violated, and remove that violation along the p directions the factorisation already supplies. One outer step of iterative refinement against the true M then follows, which at κ2 ≈ 1.64 × 1011 is a guard rather than a nicety. Accuracy of the corrected adjoint. With the correction, all four parameter classes agree with a reference computed by a sparse LU factorisation of J⊤ to a relative error of 1.4 × 10−10 ; the GPU adjoint agrees with the CPU adjoint to 7.3 × 10−11 over mixed batches whose scenarios hold different valve states (active, open and closed PRVs in one batch); and it agrees with central finite differences (each perturbed point solved through the full status machine, switch-crossing coordinates discarded) to 7.0 × 10−7 at worst. A terminal refactorisation is therefore required: at κ ∼ 1010 –1012 , with a factorisation from a stale status configuration the contraction factor of the iterative refinement exceeds one and the refinement diverges (BWSN Network 1). Finally, what the GPU adjoint cannot do it refuses loudly rather than degrading: pump-parameter gradients, Darcy–Weisbach, the clamped branch of constant-power pumps, cascaded PRVs sharing a downstream node, float32 and second derivatives all raise at construction or call time.

Supplementary Note 4: the leak perturbation operator and the signature dictionary Three leak paths, one sparsity pattern. EPANET 2.2 represents a leak in three superposable ways. (a) An emitter, a pressure-dependent orifice discharge qE,i = Ci pγi with pi = Hi − zi and γ a global exponent, 1/γ internally inverted into a virtual link to a virtual reservoir at the node elevation with head loss h−zi = Ke,i qE −1/γ and Ke,i = u Ci for a unit-conversion factor u; assembly adds 1/gE,i to Aii and (hL,E + zi )/gE,i to Fi . (b) Pressure-driven demand, the Wagner function inverted into a virtual link to a virtual reservoir at zi + pmin (Germanopoulos, 1985, Wagner et al., 1988), structurally identical to (a). (c) A known discharge added to nodal demand, di ← di + Qleak,i , which changes only the right-hand side. None of the three touches an off-diagonal entry of A. The exponent γ is held fixed throughout this work and only the coefficient Ci is 1/γ estimated: γ enters the emitter residual through qE,i , so treating it as a second unknown per candidate would double the parameter count while the data of one operating window constrain mainly the product’s magnitude. A leak is a diagonal rank-one update.

Inserting a leak at junction i therefore acts as

A 7−→ A + βi ei e⊤i ,

F 7−→ F + ϕi ei ,

(10)

with βi = 1/gE,i ≥ 0 and ϕi = (hL,E,i + zi )/gE,i for path (a), the analogous quantities for (b), and βi = 0, ϕi = −Qleak,i for (c): a rank-one, diagonal, sign-definite update of a symmetric positive definite matrix, plus a single-component update of its right-hand side. Three consequences follow. 21

For a known discharge the matrix is unchanged, so ∂H/∂Qleak,i = −A−1 ei : the i-th column of −A−1 , obtained by one solve against an existing factorisation, and by symmetry also the sensitivity of node i to a unit demand anywhere. Because the update in (10) is rank one, Sherman–Morrison applies to the emitter case (Sherman and Morrison, 1950), (A + βi ei e⊤i )−1 = A−1 −

βi A−1 ei e⊤i A−1 , 1 + βi (A−1 )ii

(11)

so a candidate sweep over M nodes replaces M factorisations by M solves, and a batch of candidates is a low-rank rather than a rank-one update with the same conclusion. And at fixed pressure ∂qE,i /∂Ci = pγi , so through the coupled system the adjoint of Note 2 returns ∂L/∂Ci for all candidates simultaneously with one solve rather than one forward simulation per candidate. That is the operational difference between derivative-free leak search and gradient-based leak inversion. The signature dictionary. Fix a sensor set S and a candidate set C, and define the signature of candidate i as the normalised column ∂hS . ∂hS ∈ R|S|·T , (12) di = ∂Ci ∂Ci 2 stacked over T operating frames. The matrix Φ = [d1 , . . . , d|C| ] is the dictionary that any sparse localisation method (ℓ1 minimisation, matching pursuit (Tropp and Gilbert, 2007), or Bayesian selection in a surrogate’s latent space (Mücke et al., 2023)) implicitly works with, and its mutual coherence µ = maxi̸=j |⟨di , dj ⟩| controls whether sparse recovery can succeed at all. Without differentiability Φ must be built by finite differences, one simulation per candidate per frame; here it is a by-product of the adjoint, and the 60-candidate dictionary of the main text costs 29 s of single-process CPU time. Two medians, two denominators. The coherence statistics of the main text’s Section 3.5 use two denominators, and neither may be quoted against the other’s population. The median over all 1,770 candidate pairs is 0.825840; that population includes 275 cross-zone pairs which are exactly orthogonal, because 5 candidates sit behind a PRV and a pressure signal does not cross the closed boundary. The median over the 1,495 within-zone pairs is 0.871. Both are true; they answer different questions, and the 42 pairs above 0.999 are 42 of 1,770 overall and 42 of 1,495 within-zone: all of them within-zone. The dictionary is not an artefact of the stopping rule. The L-TOWN input file ships with an ACCURACY setting of 10−2 , and both the base and the probe solves behind the reported dictionary are tightened well past it. As a control, the whole dictionary is built twice, at the shipped setting and tightened. The coherence statistics are insensitive to the choice: the maximum is identical to nine significant figures, the count of pairs above 0.999 is identical at 42, the median over all 1,770 pairs moves by 6.7 × 10−4 , and each true leak’s closest-rival coherence agrees to five significant figures. One thing does change, and it is the reason for tightening: at the shipped setting the solve leaves spurious non-zero sensor responses across the valve boundary, so the exact orthogonality of the cross-zone pairs is destroyed and the pair count is inflated. The coherence conclusion is therefore robust to the stopping rule; the zone decomposition quoted above is not, and requires the tightened solve.

Supplementary Note 5: the batched status machine EPANET’s hydraulic loop interleaves two cadences of discrete logic: valve statuses are re-examined on every iteration, whereas general link status checks run every CHECKFREQ iterations (the engine’s status-check frequency option) and once more at convergence. The batched port reproduces both cadences literally, and 22

does not bucket scenarios by status. For each controlled link the coefficient contributions of all status branches are evaluated for the whole batch and combined by masks; the transition rules are applied as one-hot selections over the same four cases as the reference routine; and the loop terminates only when every scenario in the batch is simultaneously converged and status-consistent. An active PRV replaces its downstream row by the penalty construction of Note 3, and the port writes the penalty on the diagonal only, so that A remains bitwise symmetric. That is not an assumption but a checked invariant: the regression suite (Table S5) tests max |A − A⊤ | on every Newton round of 86 network–path configurations (the operating network included) and finds it exactly zero, and four seeded one-sided-assembly variants of the port, among them “write the penalty row one-sided”, are each detected. The invariant the guard checks is a property of the assembly code, not of any network; the two large public benchmarks on which parallel links in mixed directions perturb A − A⊤ at the last bit (main text, Section 2.5) lie outside the batched path by capability, not by exemption. Agreement with the serial replica is behavioural, not merely terminal. On 134 valve-forcing scenarios of L-TOWN, chosen to drive each PRV through all three of its states, the batched and serial paths reach element-identical terminal states in 134 of 134, and their non-empty status-transition subsequences agree item by item in 131 of 134; the three divergences are stopping-time-only, one path declaring convergence one iteration earlier, and are not transition disagreements. Across five public PRV networks every divergence is attributable either to a status decision sitting within rounding distance of its threshold or to accumulated drift on penalty-degenerate branches, and the count of cases in which the two paths saw the same decision inputs and transitioned differently is zero. With the status machine the batched path covers 18 of the 21-network public suite (11 without it). The remaining three need valve types or head-loss laws the batched path does not implement: Net6 and BWSN Network 2 are outside it because of a capability gate and not because of a memory limit: Net6 combines check-valve pipes (the engine’s CVPIPE link type) with a PRV, and BWSN Network 2 combines them with pressure-sustaining and flow-control valves (PSV and FCV).

Supplementary Note 6: sparsity-pattern invariance and the stateful sparse solver The invariance, and why it holds. The Schur complement is A = A21 D−1 A12 with A12 = A⊤21 the junction incidence and D = diag(gk ). Entrywise, for i ̸= j, X X Aij = − gk−1 , Aii = gk−1 + (emitter term), (13) k: {n1(k),n2(k)}={i,j}

k∋i

so the off-diagonal pattern of A is exactly the adjacency of the junction subgraph and does not depend on the values gk . Every device the status machine can apply changes only those values: a closed link is given gk−1 = CBIG−1 = 10−8 rather than being deleted, a fully open valve 106 , an active PRV a CBIG term written on the diagonal. By (10) the same is true of a leak on any of its three paths. Hence the pattern, the elimination ordering, the fill-in pattern and the symbolic factorisation are invariant under the insertion, removal or modification of any leak and under any status change, and the invariance is structural rather than approximate. The one thing it requires of the implementation is that closed links keep their slots, which is why the reference engine’s choice not to delete them matters here. What the sparse route does with it. The sparse route (compressed-row assembly handed to the vendor sparse direct solver; deliberately not the default) builds the compressed-row structure once at construction, asserting at that point that every slot a PRV can write already exists, and thereafter updates values alone. The assembled values are not approximately but exactly those of the dense path: on L-TOWN all 17 iterations of the nominal frame give max |∆A| = 0 against the dense assembly, and the terminal heads, flows, statuses and iteration counts are bit-identical when both are handed to the same linear solver. 23

The plan cache and its capacity rule. The vendor sparse direct solver is used statefully: the symbolic analysis is planned once per batch shape and cached under a least-recently-used cache keyed by batch size, floating-point precision, device and gradient slot. The rule that matters in training is (distinct batch sizes passing through the solver within one optimiser step) × (gradient slots) ≤ cache cap. (14) Below that, the cache evicts a plan that is about to be needed and silently re-plans every step. Measured over six training shapes on two public networks, the rule is broken by more than the gradient-accumulation case usually described: an ordinary bucketed loader that runs the backward pass immediately after each forward pass still triggers it as soon as the number of buckets exceeds the cap, costing 3.99–4.45× the step time on ky4 and 2.33–2.43× on Modena, with the counter signature unambiguous (16 plans rebuilt per epoch at the default cap, zero at the correct one). Raising the gradient-slot count is not the remedy: it multiplies the number of cache keys, and in every measured cell its best outcome was inside node-to-node noise while its worst was 24.3–25.8× slower and 2438 MiB more resident. The sparse route refuses rather than degrades on float32, CPU execution, second derivatives and a missing vendor solver library.

Supplementary Note 7: the cost crossover, and why the end-to-end factor must be split The crossover is real and it is not one number. The sparse-versus-dense ratio was measured on 9 public networks at six batch sizes, forward-only and forward-plus-backward, on three independently measured RTX 5090 nodes (the node-to-node spread of any cell is at most 10%). Three readings follow, and the third is the one an implementer needs. First, for the forward solve alone the crossover at B ≥ 64 lies between Nj = 268 (Modena, 0.71–0.73×, sparse losing) and Nj = 959 (ky4, 6.6–6.9× at B = 1,024, sparse winning); at B = 8 it has already moved below 268 (Modena 1.21–1.24×). Second, for forward-plus-backward the crossover is lower, between Nj = 92 and Nj = 268: Net3 is level at B = 8 (0.97–1.04×) and loses above it, Modena wins at every batch size. No public network in the dense-capable set has 92 < Nj < 268, so an interval rather than a point is reported. (On the scale axis of pure single-frame forward cost, a different calibre, the main text’s Section 3.2 places the crossover near 300 junctions, with the 541-junction operating network the first measured point above it.) Third, on small networks the sparse route loses badly and the loss grows with B: down to 0.06× at B = 1,024 on the nine-junction Net1. The mechanism is visible in the raw timings (the dense per-scenario cost keeps amortising with B while the sparse direct solve is a floor that rises with B and is nearly independent of Nj ), so this is a structural property of the route, not a tuning artefact. The end-to-end factor is a product of two independent things. Table S6a reports the L-TOWN wall time per scenario, forward and forward-plus-backward, with the decomposition F1 F2 : F1 is the gain from moving the serial CPU adjoint onto the GPU at the same dense linear algebra, F2 the further gain from sparse-versus-dense linear algebra. To confirm that the split is a property of the linear algebra rather than of the pipeline, one round of linear algebra was additionally timed three ways on the same assembled system: (A) dense with generic reverse-mode autodifferentiation through a dense Cholesky factorisation, which is what the shipped dense path does; (B) dense with a hand-written adjoint reusing the same factor, which exists only in the measurement script and deliberately not in the package; and (C) the sparse route through the vendor direct solver. All three compute the same gradient, agreeing to between 4.3 × 10−16 and 7.1 × 10−9 , so the comparison is meaningful. The decomposition is unambiguous. On ky4 the sparse linear algebra genuinely wins, 4.11–4.36× at B = 64. On Modena it does not: 0.21–0.22× at B = 1,024, that is, the sparse route loses by a factor of nearly five. The end-to-end numbers those two networks show (50.4–54.0× on ky4 at B = 64 and about 2× on Modena) therefore come from the second factor almost entirely: the debt the dense path’s generic autodifferentiation carries, 10.6–10.9× on Modena and 12.1–12.5× on ky4. That debt grows 24

monotonically with Nj and with B, from a fixed overhead of about 1.97–2.11× on small networks, because the vector–Jacobian product of a Cholesky factorisation costs O(B Nj3 ) triangular solves and matrix products where a hand-written adjoint is one back-substitution. The correct statement for Modena is therefore “the dense path’s backward pass is expensive”, not “sparse is faster”. The product is accordingly never quoted as a single number, and the same applies to the F1 F2 split of Table S6. Memory and admission. Table S6b gives peak device memory, measured as the tensor library’s reserved peak plus device residency outside the tensor library, with the CUDA context subtracted, in one fresh process per cell; the two measuring nodes agree bitwise on every cell. The admission consequence is the one that changes what can be run: on L-TOWN at B = 1,024 the sparse forward-plus-backward path holds 1 286 MiB against the dense path’s 29 096 MiB, and where on a 32-GiB card the dense forward first fails at B = 1,280 in a fresh process, the sparse forward-plus-backward path was scanned to B = 16,384 without reaching a boundary. Two cautions accompany this. Peak memory is a fresh-process figure: within a long-running process, allocator fragmentation at large B can exhaust the device below the structural bound the memory table reports. And a ratio of peak memories depends on which layer is being compared: on ky4 at B = 64 the sparsity-pattern activation ratio Nj2 /nnz is 286, the tensor library’s peak-memory ratio is 57.9× and the ratio that actually decides admission is 33.6×. Absolute cost against the reference engine. Table S7 times one steady-state frame on 16 networks: 14 public models spanning two and a half orders of magnitude in size, plus the two operating networks. On the public rows the replica path costs 26× (Net3) to 73× (Fossolo) the reference library’s, with fitted exponents 1.01 for our solver and 0.99 for the library, so the gap is a constant factor rather than an asymptotic one; the smallest factor anywhere in the table is the operating network City D at 21×, with City H at 47×. What that factor buys is a derivative, which the reference library does not supply at any price. These are single-machine timings and are not a portable performance claim.

Supplementary Note 8: the calibration comparison protocol Problems. Calibration estimates Hazen–Williams roughness coefficients from noisy junction heads by minimising a mean-squared head misfit with box constraints C ∈ [40, 160] and prior regularisation. On Hanoi, 34 free pipes are estimated from 25 synthetic-diurnal frames (20 training, 5 validation) at 25 training sensors, with σ = 0.1 ft Gaussian head noise regenerated per seed; on City D, 432 free pipes (of 475; 41 dead branches and 2 all-frame-clamped pipes are frozen as structurally unidentifiable) from the utility’s own 24-hour diurnal pattern, 32 training sensors of the 40 designated, and the same noise model. Table S8 reports the full noise ladders. The gradient calibrator is Adam, then L-BFGS, then a Levenberg–Marquardt polish whose Jacobian comes from the adjoint of Note 2. The Modena rows. The public Modena rows of Table S8 are the two σ = 0 runs (per-pipe and grouped truth) and nothing else: the truth was generated by the same model the inversion searches, so both are inverse crimes by construction and are labelled so in the table. They document that the machinery runs on a third network; they are not evidence of calibration under uncertainty, which is why the main text’s calibration evidence is Hanoi and City D, whose ladders carry measured noise arms (σ ∈ {0.03, 0.1, 0.3} ft). The σ = 0 rows of Hanoi and City D are inverse crimes for the same reason and are labelled identically. Baselines and tuning. Five derivative-free baselines, differential evolution (DE), particle swarm (PSO), CMA-ES, simulated annealing (SA) and a DE→Levenberg–Marquardt hybrid, were each tuned over four

configurations and three dedicated seeds before evaluation: ≈240,000 model calls of tuning each on Hanoi, except SA, which is serial, cannot use the batched forward, and received 60,576 (reduced and declared); on City D the tuning budgets were 3,168–3,600 calls. Evaluation then used 30 fresh seeds (Hanoi) or 5 (City D) at matched model-call budgets {200, 600, 1,000, 2,000, 5,000}, extended to {10,000, 20,000} for DE, PSO, CMA-ES and the hybrid; backward passes are costed at their measured ≈1 % of a forward call. Baselines received batched forward evaluation, a more converged forward solve, ground-truth-prior initialisation and generous snapshots. On City D the baselines are additionally tuned at 2,000 calls, equal to their evaluation budget (Note 13); the one configuration where a baseline wins (City D, hybrid at 2,000 calls) is reported in the main text as a baseline win. Censoring, stated per method. SA’s run-length ceiling is its own declared budget: its hit-rate row is right-censored at 5,000 calls, and the statement “never reaches the endpoint within 20,000” is measured only for DE and PSO. Hit rates are metric-specific: on the training loss DE, PSO and SA are at 0 % where on the identifiable-subspace parameter-error RMSE (the roughness-error projection onto the well-conditioned subspace of the sensitivity matrix, a parameter-space quantity distinct from any pressure RMSE) the same table reads hybrid 73.3 %, CMA-ES 26.7 %, DE 10.0 %, PSO 6.7 %, SA 6.7 %, so every hit-rate quoted in the main text names its metric. Statistics. Effect sizes are Vargha–Delaney A12 with two-sided Wilcoxon rank-sum tests. Both implementations were cross-checked against brute-force re-implementations (direct pairwise counting for A12 ; exact rank-sum enumeration at these sample sizes) and agree exactly; the attainable two-sided floor at n = 5 per arm is p = 0.0625, which is why the City D baseline win is quoted with exactly that p-value. Robustness arms on City D (demand model error, sensor bias, and eight-start multi-start) are in Table S8b, including the one-in-four demand-uncertainty repeat that degrades to 1.62 × 10−2 .

Supplementary Note 9: the verification protocol in full Forward fidelity. The sweep comprises 52 networks: 23 public benchmark models, 3 EPANET distribution examples, 23 randomly generated networks and 3 operational rows (City D, its emitter variant and City H), all compared against the double-precision EPANET 2.2 dynamic library obtained through WNTR, frame by frame over the full extended-period simulation, 8,140 frames in all. Acceptance requires max |∆H| < 10−6 ft, max |∆Q| < 10−6 cfs and per-frame equality of the integer time sequence, the link statuses, the valve and pump settings and the Newton iteration counts. Table S1 lists every network and every measured maximum, and Figure S1 draws nine of them, three orders of magnitude of node count in one plate: the public panels are the verified files themselves, resolved through the same index the alignment harness reads, and the two operating panels are drawn from the released model files, whose coordinates are metres under a rigid-body transform carrying no georeference, so the plate can be redrawn from the published data alone. With inputs read exactly (Note 1), all 52 networks meet the criterion: every one agrees to ≤ 10−12 ft and 25 agree exactly, against 20 under the default input path; the largest deviation is 1.14 × 10−13 ft in head and 1.42 × 10−14 cfs in flow, and the three operational rows sit at 1.421 × 10−14 ft. Under the default input path the corresponding figures are 47 networks at ≤ 10−12 ft, 20 exact, and a worst head deviation of 6.836 ft on the largest network (Table S13). The suite spans both implemented head-loss laws (51 Hazen–Williams, 1 Darcy–Weisbach) and three flow-unit systems (37 LPS, 13 GPM, 2 CMH); the other seven unit systems EPANET supports are unit-tested only. Every value quoted in this document, Table S1 included, is the maximum over all frames.

Gradient correctness. Four independent checks are used, and they differ in what they would catch. Against Richardson-extrapolated central finite differences the implicit adjoint reaches a worst relative error of 5.85 × 10−8 over 61 coordinates and four parameter classes per network, the operating network included, at the same threshold as the public one. The Darcy–Weisbach path is differentiated analytically through all three flow regimes and checked separately over 24 coordinates: 1.15 × 10−8 in the fully turbulent regime and 2.85 × 10−6 on the laminar and transitional branch, where the gradient magnitudes are themselves of order 10−8 and the absolute agreement is ∼ 10−13 . The external checks (Table S3) are independent of the replica’s own arithmetic: a loss is computed by EPANET’s own compiled library, with demands and settings written through the programmatic interface and heads read back, and differenced against our analytic gradient. On the public emitter-augmented Hanoi configuration, eight coordinates spanning emitter, demand, reservoir-head and roughness classes agree to a worst 5.56 × 10−5 against a 1.00 × 10−4 threshold (Table S3a); on City D the roughness class is checked in breadth (24 non-clamped pipes spanning 1.7 decades of sensitivity to a worst 4.78 × 10−5 , plus six clamped pipes proven signal-free, their finite-difference residuals of at most 6.45 × 10−7 matching the independently predicted noise floor; Table S3b), and the City D emitter variant covers four parameter classes to a worst 2.39 × 10−5 (Table S3c). The finite-difference scheme takes symmetric steps at several fractions of a per-coordinate base and fits the odd model dL(δ) = 2gδ + bδ 3 by least squares, with denominators taken from measured applied values rather than nominal ones. Finally, 22 hand-constructed degenerate cases (Table S4), 6 of them on the operating network’s closed-link neighbourhoods, throttle-valve endpoints and mixed batches, are each checked against their own threshold; all pass, the worst reaching 0.57 of its own threshold, and 7 are exactly zero because the correct gradient is structurally absent: the low-flow linearisation, closed links and valves make the derivative of head loss with respect to resistance analytically zero, and naive autodifferentiation of a re-implementation reproduces those zeros without warning. The regression suite. Table S5 lists the suite item by item: 54 items, 53 numerical checks and one packaging guard. The 53 numerical checks all pass (391 s for the full suite on the reference machine). 7 of the numerical items exercise the two operational models (alignment, replay, extended-period simulation, physical crossvalidation, batch consistency and the three-way gradient cross-check) and are marked in the table. The packaging guard scans file contents, paths and version-control metadata for thirteen name variants of the two utilities and requires zero hits, the mechanical form of the anonymisation rule stated at the head of this document; over every file, path and compiled PDF of this submission package it reports zero.

Supplementary Note 10: sensor augmentation, calibration and leak search with the augmented sets, and the control experiments Every design here is virtual. New sensors are simulated on the calibrated model; none has been installed, and field installation is a later stage. The augmentation mode of the main text (Section 2.7) holds an installed set S0 fixed and adds k junctions one at a time from the candidate pool, either for the Bayesian D-optimal objective of the placement link (the same σprior = 15 and σnoise = 0.1 ft) or for coverage, the number of target pipes whose maximum sensitivity over sensors and frames crosses the census threshold; threshold, frames and structural mask are those of the census, so “restored” means exactly the census’s transition from an all-zero to a non-zero sensitivity column. At every step an assertion checks that the identifiable set has not shrunk, that no per-pipe posterior variance has risen and that no pipe visible to S0 has been lost; the curves carry zero violations. On City D, S0 is the 40 designated sensors, the target set is the census’s 149 unobservable pipes, and the candidate pool is all 541 junctions, which includes the three work-order leak junctions (a fact that matters in the leak search below). On L-TOWN, S0 is the 33 sensors of the main-text inversion, the target set is the 62 pipes they cannot inform, and the pool excludes the three injected leak

junctions, as the original sensor draw did. On Hanoi, S0 is ten random sensors over the 25 synthetic frames of the calibration ladder, with one unobservable pipe. Curves (Table S9a, b). On City D coverage restores 37/62/95/133/149 of the 149 at k = 5/10/20/40/80, 80 % at +32 and all at +56 (the coverage sequence saturates at step 57 and reverts to the D-optimal objective thereafter); D-optimal restores 21/23/35/52/76, and its well-conditioned subspace CRLB trace falls from 742 to 222 at +20 while that subspace widens from 20 to 36 dimensions. Reselecting 40 + k sensors from scratch recovers more of the 149 at small k but loses 38/38/24/21/15 pipes the designated sensors could see; augmentation loses none. The eps-rank CRLB of the main text’s Section 3.6 is dominated on City D by directions at the rank threshold and moves by up to 10−6 relative between platforms, so only counts and the well-conditioned trace are quoted for the augmented sets. On L-TOWN coverage restores all 62 at +36 (80 % at +24) where reselecting 33 from scratch loses 99; on Hanoi the single pipe is restored at +1 (coverage) or +4 (D-optimal), and the two objectives select identical sets at k = 5, 10 and 20. An independent re-computation of the City D curve, with a separate Fisher-matrix and census implementation on two hosts, reproduces every count, the well-conditioned traces agreeing to 10−9 relative. Calibration with the augmented sets (Table S10). The calibrator, truth, noise realisation and 20 % sensor hold-out of Note 8 are kept and only the sensor set changes. “Restored” pipes are the target pipes that cross the census threshold under the augmented set; their error is |Ĉ − Ctrue | and the prior’s error is |130 − Ctrue |. At σ = 0.1 ft the restored pipes move off the prior (median 11.5 against 21.9 under D-optimal +20; 15.5 against 19.1 under coverage +20); the RMSE over all informative pipes moves from 25.6 to 23.2–24.3; at σ = 0.3 ft the +20 sets are worse than S0 alone (29.5–30.3 against 26.3) and coverage’s restored pipes sit above their prior. One noise seed per cell on City D; the Hanoi ladder (three seeds) moves from 18.7 to 7.0 at +20 and σ = 0.1 ft. Leak search with the augmented sets (Table S11). The work-order inversion of the main text (49 candidates, three injected leaks, σ = 0.1 ft, the same regularisation and step budget, the same inverse-crime setting) is repeated with S0 ∪ Sk . Two facts about the baseline fix how the comparison is read. First, the 40-sensor set of the main-text demonstration and the census set S0 are different draws sharing 22 junctions; both return none of the three under this experiment’s noise realisation. Second, that noise field is generated over all junctions and read at the sensors, so every set sees one realisation, distinct from the realisation of the main-text demonstration; the “none of three to one of three” comparison is within this experiment. Each of the six augmented sets returns L2 with a discharge error of 1–4 % and nothing else. The control experiments attribute that hit: the coverage sequence places a sensor on L2 ’s own junction at its third step and the D-optimal sequence at its fifth; removing that one sensor from coverage +20 returns none of three; S0 plus that one sensor alone returns the same one of three (3 % error); twenty random additions from the same pool recover L2 in none of five draws, and forty in none of one. L1 (3.0 L s−1 ) changes no junction head by more than 0.09 ft under any set, below the 0.1 ft noise, so no sensor set can recover it. L3 is recovered by 3 of the five random +20 sets and by D-optimal +20 once the L2 -junction sensor is removed, but by none of the six designed sets, nor by S0 plus a sensor on its own junction; its nearest rival’s coherence is 0.9994–0.9997 in every set, and the designed sequences never place a sensor within three hops of it in 80 steps. Coherence medians are over all candidate pairs; on City D no pair is exactly orthogonal, so the two denominators of Note 4 coincide. Augmentation lowers the median from 0.928 to 0.79–0.90 and the pairs above 0.999 from 48 to 18–34, without lowering the rivals that matter. On L-TOWN the ten augmented sets, passed through the main-text inversion in a noiseless and a 0.1 ft group, all return the same one of three as S0 , with the two lost leaks at rank 19 and between ranks 27 and 32 of 60 and their rivals above 0.9985 in every set (Table S11b). 28

Coherence-driven placement against random controls (Table S12). A second placement objective, greedily adding sensors to minimise the same signature dictionary’s coherence (the diagnostic of Note 4, not the census’s Fisher information), was run on City D at k = 20 and k = 40 and compared against budget-matched random draws over the same candidate pool and rejection rule as Note 10’s augmentation curves: 20 fair seeds at each k (seeds that draw a leak junction are discarded and redrawn, so every random set is a valid candidate set), plus eight sensor sets chosen purely to minimise coherence (the greedy coherence sequence at 5, 10, 20, 40 and 80 added sensors; its worst-pair variant, which minimises the largest pairwise coherence, at 20; and its full-pool variant, drawn without the leak-junction rejection, at 20 and 40), for 40 random draws in total. The pooled significance tests below go further and ask a stricter question: is L2 ’s recovery a property of any purposefully chosen sensor set, or specifically of a coherence-minimising one? They pool those eight coherence-minimising sets with the six identifiability-driven augmented sets above (the coverage augmentations at 20, 40 and 80 added sensors and the D-optimal augmentations at 20, 40 and 80), 14 designed sets in total, against the 40 random draws – a design pool that already recovers L2 in every one of its six identifiability-driven members (Table S11a), so a coherence effect would have to survive being tested alongside that fact, not instead of it. The result reads differently leak by leak (Table S12a). L2 is recovered by 12 of the 14 pooled designed sets against 3 of the 40 random draws, a one-sided Fisher exact test of 1.05 × 10−7 ; of the eight coherence-minimising sets alone, six recover L2 (the greedy sequence’s 5- and 80-sensor sets do not). What actually sorts recovery is not which objective chose the set but each set’s own rival-coherence gap: across all 55 reconstructible configurations (every pooled designed set, every random draw, and every other City D sensor set this note or the main text records) the gap of every configuration that recovers L2 exceeds the gap of every configuration that does not, with no overlap – true of the identifiabilitydriven recoveries as much as the coherence-minimising ones, which is why the gap, not the label, is the mechanism reported in the main text. At the matched k = 20 budget alone (the worst-pair variant with 20 added sensors against 20 random k = 20 draws) the same conclusion holds by an exact permutation test (4.76 × 10−2 , none of 20 random draws recovering L2 ); at k = 40 (the greedy coherence sequence at 40 added sensors against 20 random k = 40 draws) the matched-budget test alone does not reach significance (1.90 × 10−1 , 3 of 20 random draws recovering L2 ), and the shortfall is not one of sample size: the random group’s 95 % Clopper–Pearson interval on its own hit rate is [0.032, 0.379] (pooled over both budgets, [0.016, 0.204], Table S12b), wide enough that additional random draws at k = 40 alone would not settle it either way. L3 ’s recovery is statistically indistinguishable from random placement (3 of 14 pooled designed sets against 7 of 40 random draws, p = 0.51) and its rival-coherence gap ranges overlap (Table S12a); it is not a coherence effect, and no claim in the main text credits it as one. L1 is recovered by neither objective at any budget: its dictionary column norm (0.15–0.64) sits an order of magnitude below L2 ’s (5.1–17.1), so the 0.1 ft noise floor swamps its signature however far the rival coherence is pushed down (to 0.121 at best on this pool); it is an amplitude limit, not a coherence one. The two inversion stages. The work-order inversion above adds a discrete support-refinement stage (nonlinear orthogonal matching pursuit and a support-swap polish) after its Adam-plus-ℓ1 coarse stage; the L-TOWN inversion routine of the main text’s Section 3.4 (Adam with an ℓ1 penalty) has only the coarse stage, ranking candidates by coefficient magnitude with no discrete search. Re-scoring every one of the 55 City D configurations on the coarse stage alone – whether the truth is among the top three by magnitude, the criterion the L-TOWN inversion routine itself reports – and repeating the Fisher and gap tests on that coarser judgement leaves both conclusions unchanged: L2 stays separated (2.10 × 10−6 , 13 of 14 coherence-selected sets against 8 of 40 random draws, rival-coherence gap still non-overlapping) and L3 stays statistically indistinguishable from random (8.92 × 10−1 , 5 of 14 against 20 of 40). The discrete support-refinement step changes zero L2 configurations from lost to recovered across the 55 and six from recovered to lost, so it does not produce the 29

coherence result; it only removes coarse-stage hits the flow error then fails to confirm. One further asymmetry appears only at the coarse stage: coherence-selected sets rank L1 in the top three by magnitude significantly more often than random draws do (2.39 × 10−3 , 6 of 14 against 2 of 40), yet not one of those coarse hits survives the discrete refinement step under the final flow-error criterion, where both groups still recover L1 zero times. Coherence can move L1 up the coarse ranking without ever producing an accurate recovery, which is consistent with, not contrary to, its being amplitude-limited rather than coherence-limited. Why no placement closes the L-TOWN loop (Table S12c). The same candidate-pool geometry that explains City D’s graded outcome bounds what any placement on L-TOWN could achieve for the two leaks Section 3.4’s inversion never recovers (truth ranks 19 and 32 of 60). Of the full 746-position fair candidate pool, only 1 and 14 positions, respectively, can push that leak’s coherence with its nearest rival below the 0.999 threshold Note 4 uses; an oracle greedy search dedicated to just that one pair, run for 80 steps and free to ignore every other leak and every budget constraint, drives the pair’s own coherence down to 0.996775 and 0.990741 but leaves the next-closest rival at 0.999400 and 0.997689, both still inside the inseparable regime. The one leak L-TOWN’s inversion always recovers has 746 of 746 fair positions that separate it from its rival, and an oracle search on that pair alone reaches 0.903043. No placement objective, on this pool, buys the other two leaks anywhere close to that margin. The conclusion carried to the main text is therefore: placement for roughness identifiability does not close the leak-search loop, but a second objective, placement against coherence itself, closes it for one of the three leaks for the reason coherence predicts, is indistinguishable from random placement for a second, and cannot close it for the third because that leak’s limit is amplitude, not coherence; the asymmetry between the two inversions used in this note does not drive the coherence result, and the candidate-pool geometry that bounds City D’s coherence effect also bounds what placement can do on L-TOWN.

Supplementary Note 11: from single-point localisation to cluster-level diagnosability The construction. Note 4’s dictionary makes the failure of single-point localisation measurable; it also says what is resolvable, because two candidates whose signatures are coherent above a threshold are exactly the pair the data cannot separate. Group them, and the inversion answers the question the dictionary can answer. Candidates are clustered by constrained agglomerative complete linkage: two clusters may merge only if they are adjacent in the network’s own Voronoi partition (multi-source Dijkstra over pipe lengths, no tuning parameter), and the merge is accepted only while the minimum cross-coherence stays at or above τ , so at termination every pair inside a cluster has µij ≥ τ . The inversion then replaces the non-negative Lasso by its P p group form, minx≥0 12 ∥Ax − y∥2 + λ g |g| ∥xg ∥2 , solved by FISTA with a non-negative group proximal step; with every group a singleton it reduces exactly to the single-point inversion, which is the control reported beside it. λ is chosen by a discrepancy rule that never sees the truth: sweeping λ downwards with warm starts, √ √ the first value whose residual sum of squares falls below max{σ 2 (n + 2 2n), min RSS + 2σ 2 2n}, the second term a floor because a linearised dictionary does not in general reach the pure-noise level. Clusters are ranked by ∥xg ∥2 and reported with the radius, diameter and connected pipe length of the district they name. On the operating network (542 nodes, 554 links, 49 candidates, 1,000 observations) the whole pipeline costs 0.06 s for the pairwise pipe distances, 0.002 s for the clustering and 1.94 s for a 14-point warm-started λ path on one CPU core. What it buys, and what it costs. Table S14 is the trade-off curve, and the curve is the result: a hit rate quoted without the radius that produced it is not one. At τ = 0 the operating network is a single cluster of radius 9,043 m covering all 83.9 km of main, and the top-3 hit rate is 3 of 3, the hit rate of saying “somewhere

in the network”. At the other end, at τ > 1 every cluster is a singleton and the reading is the single-point inversion, 0 of 3. Between them, at τ = 0.999, 25 clusters of mean radius 179 m put two of the three injected leaks inside the top three; those three clusters hold 7 of the 49 candidates and 7.42 km of main, 8.8 % of the network, the largest of them 750 m in radius. Panel (c) gives that inspection burden at five thresholds, because the mean radius over all clusters is dominated by singletons and understates what a crew would walk. On L-TOWN the corresponding point is τ = 0.995: 22 clusters, mean radius 98 m, 1.96 km of district, all three leaks inside the top three. Four things this result is not. First, the thresholds quoted are chosen after seeing the whole sweep; both sweeps are printed in full so that the choice can be discounted. Second, the design does not beat chance on hit rate. Against 20 size-matched random groupings that preserve the cluster-size distribution while destroying topology and coherence, the one-sided exchangeability p never falls below 0.09 on either network (panel d); what the design wins is the radius, 179 m against 1,700 m on the operating network and 40 m against 298 m on L-TOWN at the same hit rate. Third, the equal-budget single-point control is not a formality: given the pooled candidate count of the three ranked clusters, the single-point ranking recovers two of three at τ ≤ 0.9 and none at τ ≥ 0.95, so the cluster reading wins only in the tight-radius regime and loses in the loose one. Fourth, the leak that noise defeats is not rescued. Its own signal-to-noise ratio ∥Ca∥2 /σ 2 is 0.50 against 166 and 151 for the other two; it is a singleton at every τ ≥ 0.5 because its most coherent rival is not Voronoi-adjacent to it; and forcing it into a group with its twelve most coherent neighbours moves it only from 49th to 36th, the group’s non-centrality parameter and its degrees of freedom rising together so that the standardised score falls from 177.9 to 159.0. Aggregation redistributes the burden of saying which; it does not create signal. On L-TOWN, where the same leak label has signal-to-noise 5.6 × 104 , the same construction moves it from 23rd to 3rd.

Supplementary Note 12: scoring a sensor layout on a ruler that does not move with it Why a layout-independent ruler is needed. The natural score for an augmented layout, roughness error over the pipes that layout makes informative or over its own identifiable subspace, is not comparable across layouts, because the set being averaged over is itself a function of the design. Over the designs compared here the informative set runs from 279 to 412 pipes and the design-specific subspace rank from 20 to 56; a design that adds sensors is scored on a larger and differently conditioned set than the one it is being compared with. Head-space misfit is not an alternative. At σ = 0.1 ft it is 96 % noise variance, and across eight optimiser restarts on identical data and an identical layout the training head misfit moves 0.4 % and the held-out-frame misfit 1.4 % while the roughness error in the reference subspace moves 25 %: the head metrics cannot resolve what is being compared. The ruler. A reference subspace is fixed once, before any design is considered, from the adjoint sensitivity spectrum of the fully instrumented candidate pool (541 junctions, 25 frames): direction vj is admitted when its posterior standard deviation would fall to σprior /γ, an absolute criterion involving no layout. At γ = 2 this √ ⊤ (Ĉ − C admits kref = 98 of the 432 directions and the score is RMSEident = ∥Vref true )∥2 / kref , against a prior reference of 28.92 (every pipe left at its initial value). The same subspace scores every layout, designed or random, and γ = 3 and γ = 10 are reported beside it. What the designs are worth on it. Table S15a is the full grid: two budgets, three noise levels, two objectives, each against random additions drawn from the same pool at the same budget. In all twelve combinations no random layout matched the design, so each cell’s one-sided exact randomisation p sits at its attainable floor,

1/(R + 1): 0.0099 in the cell run to R = 100 and 0.0476 in the five run to R = 20. That the five coarser cells stop at 0.0476 is a sample-size limit and nothing else; reaching p < 0.01 there needs R ≥ 100 and no change of statistic would substitute. R was fixed per cell in advance of looking at any p, and the whole grid rests on a single noise realisation, so the p values are conditional on it. The costs. Panel (b) reads the main cell on five rulers. The conclusion survives γ = 3, γ = 10 and the unprojected 432-pipe RMSE, but in the 334 reference directions the projection discards, 60 of 100 random layouts are at least as good as the design (p = 0.60): concentrating a finite observation budget on identifiable directions is a choice to give up the rest, and the number says so. Panel (c) prices the alternative practice: selecting the best of 100 random additions by training head MSE returns a layout ranked 25th of 100 on the identifiability ruler, by held-out-frame RMSE one ranked 56th, a coin flip. Across the whole grid there are 22 of 48 (cell, design, head-metric) combinations in which the design is ahead of every random layout on parameters and behind most of them on heads. Panel (d) relates the design-time prediction to the achieved error: the linear-Gaussian prediction correlates with the achieved error at Spearman ρ = 0.38 within a budget (n = 100) and 0.55 when budgets are pooled, so the prediction separates budget levels reliably and orders layouts within a budget only weakly. Because the D-optimal objective and this evaluation subspace both descend from the same sensitivity matrix, the coverage row, whose objective is criterion-independent, is the one that carries the general claim.

Supplementary Note 13: optimiser enhancements, the baseline tuning budget and the wall clock Device memory. The experiment driver releases the terminal Cholesky factor after each backward pass; device memory holds at 67.5 MiB across the arms. The noise floor, and which metric can resolve anything. At the City D configuration the training loss of the true parameter vector is 1.010 × 10−2 , and the plain first-order arm already reaches 9.49 × 10−3 : below the floor, fitting noise. Near the floor the loss and the parameter error are anti-correlated across arms (ρ = −0.96), so arms are compared on parameter error, and on two of them: the design’s own identifiable subspace (the placement objective’s own ruler, of rank 41 and prior 30.41 for all 77 runs, so internally consistent) and the layout-independent ruler of Note 12. What the enhancements are worth. Table S16a lists every arm. On the internal ruler the Schur-complement diagonal preconditioner is a clear win (17.71 to 12.40, better on 10 of 10 seeds, sign test p = 0.0020) while making the training loss worse on 10 of 10, which is what a noise floor implies. On the layout-independent ruler the same paired comparison is 7 of 10 with p = 0.3438 and the median gap shrinks from 4.40 to 0.41: the gain is real in the subspace the design itself defines and not established outside it, and both readings are reported because neither alone is the answer. The exact Gauss–Newton diagonal is worse than the cheap Jacobi approximation (1 of 10, p = 0.0215). Batched multi-start buys at most 2.8× per model call at B = 32, not B× (panel b), because 20 frames already saturate the device; on the parameter error, added to the preconditioner, it is worse on 10 of 10 (p = 0.0020). No arm beats the 200-call gradient configuration, which leads on both rulers. Baseline tuning budget. The City D baselines are tuned at 2,000 model calls, equal to their evaluation budget and matching the Hanoi protocol, over four configurations and three dedicated seeds; each of the four algorithms selects the same configuration as under the 3,168–3,600-call tuning of Note 8. Differential

evolution at ten times the budget (20,000 calls) reaches a parameter error of 16.66, still worse than the gradient configuration’s 200-call result. The wall clock. Panel (c) gives the other axis. The gradient configuration takes 631 s for its 200 calls; the DE→LM hybrid takes 4,756 s for 1,997, on the same node, and that row rests on one seed, not five. Per model call the baselines are the cheaper side (0.081 s against 0.113 s), so the accounting gives them no timing advantage; but dividing the hybrid’s wall clock by the eight-way parallelism it is promised in the fairness protocol gives 595 s against 631 s, a tie. The equal-budget claim therefore stands on the model-call axis, 200 against 1,997, and not on the wall clock.

Supplementary Discussion Some EPANET reference solutions must not be used as labels A bit-level re-implementation is an unusually sensitive instrument for auditing the reference engine, because every disagreement has to be explained rather than absorbed. Two of the candidate explanations are properties of the reference solution rather than of the replica, and the finding is useful independently of anything else in this paper. On ky5, controls close two pumps and isolate a dead-end branch, which then retains a stale through-flow at the input file’s own ACCURACY setting; tightening the setting reduces the residual by more than two orders of magnitude. On BWSN Network 2, three junctions form dead-end pockets that reach the main network only through pumps and flow-control valves closed for 31 or 32 of the 32 frames. They carry a stale through-flow whose head is close to arbitrary: perturbing a single roughness value inside the reference library by one unit in the last place moves the library’s own answer at those junctions by 6.857 ft. Both BWSN Network 1 and BWSN Network 2 continue to oscillate without converging when the accuracy is tightened to 10−8 over 1,000 trials. The replica reproduces these states faithfully, as a replica should, but they are simulator artefacts. The implication is general and cheap to act on. Any work that generates supervision by running EPANET over benchmark networks (Hajgátó et al., 2021, Kerimov and Taormina, 2024, Ashraf et al., 2024, 2025), and any use of large simulated corpora (Truong et al., 2025), will train on those artefacts unless it compares each frame’s residual against the input file’s own stopping criterion and flags the frames that stopped rather than converged. The same applies to the mainline benchmark used here: L-TOWN ships with an ACCURACY setting of 10−2 , whose stopping solution sits 1.7 × 10−4 ft from the converged one, and every reference solution used in this work is tightened first.

Relation to the surrogate literature Learned surrogates and a differentiable replica answer different questions and are complementary rather than competing. A surrogate is asked to be fast, can be orders of magnitude faster than EPANET, and can generalise across topologies. A replica is asked to be correct; per solve it is slower than the reference engine (Note 7), and it generalises exactly as far as its feature coverage does: to any network built from the elements the package documents as supported, a strict subset of what EPANET reads, which the sweep exercises on 52 networks rather than on all of them. Two consequences are worth stating. First, a differentiable replica is a supervision source that supplies not only H(θ) but ∂H/∂θ, so a surrogate could be trained to match the sensitivity field and not only the pressure field. No such surrogate is trained here, and nothing is claimed about what it would buy. Second, it supplies the reference against which a surrogate’s gradients, and not merely its predictions, can be audited. At present that literature reports prediction error and is silent on gradient error, even where the gradient is what the surrogate is used for (Strotherm et al., 2025). Extending EPANET without modifying it is a recurring strategy (language bindings (Eck, 2016, Arandia and Eck, 2018), plugin architectures (Sela et al., 2019), domain decomposition (Diao et al., 2014) and digital twins built on it (Park et al., 2026)), and treating the engine as a callable black box is the right decision when the quantity of interest is the output, and the wrong one when it is a derivative of the output.

Limitations, in full The main text states the binding limitations where they arise; this is the complete set, in the order in which a reader is likely to be constrained by them.

(1) Unimplemented features, and coverage is per path. Pressure-breaker and general-purpose valves, pressure-driven demand analysis, tanks with volume curves and the entire water-quality module are not implemented, and are rejected at parse or construction time with an explicit error rather than approximated. A positive damping limit (the engine’s DAMPLIMIT option, DAMPLIMIT > 0) and non-zero headerror or flow-change limits (its HEADERROR and FLOWCHANGE options) are likewise rejected. Chezy–Manning is implemented but unvalidated, because no network in the suite uses it. The batched differentiable path is narrower than the replica path: Darcy–Weisbach networks and PSV, FCV, PBV and GPV valves are construction-time rejections there, which is why Net6 and BWSN Network 2 lie outside it: a capability gate, not a memory limit. With the status machine of Note 5 the batched path reaches 18 of the 21 public benchmark networks (11 without it). (2) Bit-level agreement is platform-dependent and is a claim about one reference binary. The 10−14 ft figures were obtained on Windows x64 calling the power and logarithm routines of the Windows C runtime, against the reference engine’s Windows build shipped inside WNTR. On another platform, against a differently compiled EPANET, or against a WNTR release carrying a different binary, a single one-unit difference in the power routine is enough to move the answer and agreement degrades to roughly 10−6 ft, far inside engineering tolerance but no longer bit-level. The 10−14 figures must not be quoted for another build without re-measuring. (3) The exact-input control is what makes the whole suite pass, and it is not the default path. Reading reservoir heads and valve settings from the input text is a non-default parser entry point (Note 1). The default path, on which every other result in this package is computed, is unchanged to the last bit, and on a model written in US customary units it still loses one unit in the last place on those two field classes. Five rows of Table S1 carry that loss. (4) “All 52 pass” is not “all 52 agree bit for bit.” 25 networks agree exactly; the other 27 sit between 10−14 and 1.14 × 10−13 ft, and the exact-input control does not move them, because their residual comes from reduction order rather than from the inputs. (5) Two reference solutions are unconverged transients. ky5 and BWSN Network 2 stop rather than converge at their input files’ own ACCURACY setting, and neither may be used as a training label; the replica reproduces them faithfully, which is not the same as their being right. Net6 is verified over all 609 of its frames, with frame times, statuses and Newton iteration counts equal frame by frame (Table S13a). (6) Gradients are gradients of a frozen status configuration. Derivatives with respect to parameters that would flip a valve or pump status are one-sided at best (Note 2). (7) The unrolled route is not usable everywhere, and the refusals are explicit. The unrolled-iteration path refuses every PRV network including L-TOWN, raising at construction or at the call rather than degrading; every PRV gradient in this paper comes from the implicit adjoint. Its truncated demand gradient does not converge in K on ky4, Net3 and is borderline on Pescara (Table S2). (8) The GPU adjoint and the sparse route refuse rather than degrade. The GPU adjoint raises on pump-parameter gradients, Darcy–Weisbach, the clamped branch of constant-power pumps, cascaded PRVs sharing a downstream node, float32 and second derivatives; the sparse route raises on float32, CPU execution and a missing vendor solver library; float32 with a PRV network is refused at construction on every path. None of these fall back silently. (9) The dense and unrolled paths do not have the accuracy of the replica path. The bit-level claims apply to the replica path only. GPU double precision differs from CPU double precision by an amount set 35

mainly by the network rather than by the card, from 1.2 × 10−10 ft on a pipes-only model upward, while repeating one measurement eight times on a single card already spreads it by 1.73×; single precision is unusable, drifting by 126.6 ft with the iteration count rising from 4–5 to 15–23. The batched paths additionally carry run-to-run jitter on identical inputs of up to 5.7 × 10−6 ft in head, on both the dense and the sparse route, from the non-deterministic reduction order of the index-scatter assembly; acceptance criteria on these paths are therefore set relative to each path’s own run-to-run jitter and never as a fixed threshold. (10) Coverage is 52 networks and a documented feature subset, not “all of EPANET”. Three of the ten supported flow-unit systems are exercised bit-exactly; the other seven are unit-tested only. (11) The dense default path is slower than the reference engine, by a factor of 26 to 73 per single steadystate frame across the public rows of the 16-network size sweep, 21 to 47× on the operating rows (Table S7), and at B = 1 the GPU loses to the CPU. The large factors reported for L-TOWN are against our own prior serial CPU-adjoint pipeline at 1.0-core occupancy, not against EPANET, which supplies no gradient at any price; a perfectly eight-way parallel CPU adjoint would divide them by up to eight, and that division is arithmetic rather than a measurement. Below the crossover the sparse route is slower than the dense one, down to 0.06× at B = 1,024 on a nine-junction network. (12) The leak demonstrations are one network and one injected triple each, and the two inversions used differ in exactly one stage. The operating-network work-order inversion injects its three leaks into the same model the inversion searches (an inverse crime by construction, stated in the main text), and its noisy counterpart recovers none of them; the L-TOWN inversion is one candidate set and one injected triple, with a single-stage ℓ1 objective deliberately not augmented with the support-refinement step City D’s inversion has. Repeating both with the identifiability-driven augmented sensor sets of Note 10 recovers no further leak that a control does not attribute to a sensor placed on the leak’s own junction. A second, coherence-driven placement objective does recover one of the three leaks the noisy work-order inversion loses, for a reason the signature coherence predicts, and is indistinguishable from random placement for the other, and re-scoring that result on the coarser stage the two inversions share leaves both conclusions unchanged (Note 10): the stage the two inversions do not share is not what drives them. What generalises is the diagnosis (identifiability census and signature coherence), not the recovery, with one qualified exception (Note 10). (13) Every added sensor is virtual. The augmentation prescriptions of Note 10 are designed and verified in simulation on the calibrated model; none has been installed, and installation, re-calibration on field data and the field-side check of the restored identifiability are later work. (14) This is research code, not validated for operational use.

Supplementary Figures pipe

Hanoi 32 nodes, 34 links

reservoir

1,000 a.u.

2 km

valve

L-TOWN 785 nodes, 909 links

Modena 272 nodes, 317 links

1 km

500 m

ky4 964 nodes, 1,158 links

2 km

Net3 97 nodes, 119 links

pump

10 a.u.

City D 542 nodes, 554 links

tank

City H 921 nodes, 1,038 links

10 km

Net6 3,356 nodes, 3,892 links

50 a.u.

BWSN Network 2 12,527 nodes, 14,831 links

10 km

Figure S1. Nine of the 52 verified networks. a–i, Pipes, reservoirs, tanks, pumps and valves; junctions are not drawn. Panels are ordered by node count and independently scaled, each with its own bar: metric where coordinates carry the model’s length unit, arbitrary units (a.u.) otherwise. City D and City H are drawn from the released models; their coordinates are local, rigidly transformed.

Supplementary Tables

Table S1. Complete forward verification suite: all 52 networks against the EPANET 2.2 library, maxima over every frame. Newton iteration counts are identical on all 52 (column omitted). p constant-power pump; *custom curve; ‡ the deviation shown is the default input path’s, and this row reaches the acceptance criterion only under the exact-input control of Note 1 (Table S13). Network

Source

Anytown Anytown (WNTR) Balerma BWSN Network 1 BWSN Network 2

Public Public Public Public Public

L Loss Units Features

22 41 H-W GPM 25 46 H-W GPM 447 454 D-W LPS 129 178 H-W GPM 12527 14831 H-W GPM

pump1* pump3* tank2 pipes only pump2 PRV8 tank2 ctrl1 rule4 pump4* PSV1 FCV4 tank2 ctrl1067

Frames max|∆H| max|∆Q| Outcome (ft) (cfs) 9 6.082e-12 4.750e-12 pass‡ 1441 0.000e+0 1.776e-15 pass 1 2.842e-14 8.882e-16 pass 207 6.086e-6 7.046e-6 pass‡ 32 6.836e+0 1.406e-5 pass‡

C-Town (BATADAL) Public D-Town Public Fossolo (poly1) Public Hanoi Public ky10 Public ky3 Public ky4 Public ky5 Public L-TOWN Public L-TOWN (Real) Public Modena Public Net1 Public Net2 Public Net3 Public Net6 Public Pescara Public Richmond (skeleton) Public Richmond (standard) Public

396 407 37 32 935 279 964 430 785 785 272 11 36 97 3356 71 48 872

444 H-W LPS pump11 PRV3 FCV1 tank7 ctrl20 850 5.684e-14 8.882e-16 pass 459 H-W LPS pump11 PRV4 TCV1 tank7 ctrl24 819 5.684e-14 1.776e-15 pass 58 H-W LPS pipes only 1 0.000e+0 6.939e-18 pass 34 H-W LPS pipes only 1 1.421e-14 1.421e-14 pass 1061 H-W GPM pump13p PRV5 tank13 ctrl6 1 0.000e+0 2.220e-16 pass 354 H-W GPM pump5p tank3 ctrl2 25 0.000e+0 8.882e-16 pass 1158 H-W GPM pump2p tank4 ctrl2 1 0.000e+0 4.441e-16 pass 507 H-W GPM pump11p tank3 ctrl8 31 0.000e+0 3.553e-15 pass 909 H-W CMH pump1 PRV3 tank1 ctrl2 2031 2.842e-14 1.110e-16 pass 909 H-W CMH pump1 PRV3 tank1 ctrl2 2030 2.842e-14 1.110e-16 pass 317 H-W LPS pipes only 1 2.842e-14 8.882e-16 pass 13 H-W GPM pump1 tank1 ctrl2 27 0.000e+0 4.441e-16 pass 40 H-W GPM tank1 56 0.000e+0 2.220e-16 pass 119 H-W GPM pump2 tank3 ctrl18 183 2.146e-6 5.684e-7 pass‡ 3892 H-W GPM pump61p PRV2 tank32 ctrl124 70/609 1.525e-5 4.739e-6 pass‡ 99 H-W LPS pipes only 1 1.421e-14 8.882e-16 pass 51 H-W LPS pump7* tank6 ctrl14 91 1.137e-13 2.220e-16 pass 957 H-W LPS pump7* PRV1 tank6 ctrl16 55 1.137e-13 2.220e-16 pass

EXA4 EXA5 EXA6

Example Example Example

396 193 111

444 H-W LPS 273 H-W LPS 123 H-W GPM

City D City D (emitter) City H

Operating Operating Operating

542 542 921

554 H-W LPS 554 H-W LPS 1038 H-W LPS

TCV79

Synthetic-M0 Synthetic-M1 Synthetic-M10 Synthetic-M11 Synthetic-M12 Synthetic-M13 Synthetic-M14 Synthetic-M15 Synthetic-M16 Synthetic-M17 Synthetic-M18 Synthetic-M19 Synthetic-M2 Synthetic-M3 Synthetic-M4 Synthetic-M5 Synthetic-M6 Synthetic-M7 Synthetic-M8 Synthetic-M9 Synthetic-S0 Synthetic-S1 Synthetic-S2

Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic Synthetic

104 46 39 72 89 72 106 47 98 89 117 85 101 36 81 87 37 65 74 39 43 117 73

176 H-W LPS 67 H-W LPS 63 H-W LPS 133 H-W LPS 162 H-W LPS 105 H-W LPS 220 H-W LPS 79 H-W LPS 142 H-W LPS 178 H-W LPS 235 H-W LPS 168 H-W LPS 135 H-W LPS 52 H-W LPS 143 H-W LPS 141 H-W LPS 61 H-W LPS 105 H-W LPS 148 H-W LPS 50 H-W LPS 68 H-W LPS 181 H-W LPS 122 H-W LPS

pipes only

pump11 PRV3 TCV1 tank7 ctrl22 pump2 PRV3 tank1 ctrl4 rule1 pump1 tank1

TCV79 emit5 pump6

pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only pipes only

1 5.684e-14 8.882e-16 pass 52 5.684e-14 8.882e-16 pass 25 0.000e+0 8.882e-16 pass 25 1.421e-14 1.776e-15 pass 25 1.421e-14 3.553e-15 pass 25 1.421e-14 7.105e-15 pass 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

0.000e+0 5.551e-17 pass 0.000e+0 1.110e-16 pass 0.000e+0 5.551e-17 pass 2.842e-14 2.776e-17 pass 2.842e-14 2.220e-16 pass 2.842e-14 2.220e-16 pass 2.842e-14 1.110e-16 pass 0.000e+0 1.388e-17 pass 2.842e-14 2.220e-16 pass 2.842e-14 5.551e-17 pass 0.000e+0 5.551e-17 pass 2.842e-14 5.551e-17 pass 0.000e+0 5.551e-17 pass 0.000e+0 1.388e-17 pass 0.000e+0 1.388e-17 pass 0.000e+0 4.441e-16 pass 2.842e-14 5.551e-17 pass 2.842e-14 5.551e-17 pass 2.842e-14 2.220e-16 pass 0.000e+0 6.939e-18 pass 1.421e-14 1.110e-16 pass 2.842e-14 1.110e-16 pass 0.000e+0 5.551e-17 pass

Table S2. Health of the truncated-unrolling demand gradient on the seven public networks it was measured on: relative error of the K-step gradient against the implicit adjoint. The route refuses PRV networks at construction.

Network

rel. error at K

outcome

Net1 Hanoi Net2 Modena Pescara Net3 ky4

1.25 × 10−13 1.52 × 10−11 4.03 × 10−10 5.47 × 10−11 1.10 × 10−7 3.98 × 10−3 1.80 × 10−4

usable usable usable usable conditional not recommended not recommended

Table S3. External gradient checks against the compiled EPANET 2.2 library. (a) The public emitter-augmented Hanoi configuration, per coordinate. (b) City D roughness class by sensitivity band and (c) City D emitter variant by parameter class, summarised without identifiers per the anonymisation rule. (a) The public emitter-augmented Hanoi configuration (Hanoi + 5 emitters), per coordinate; threshold 10−4 Class

Coordinate

emitter coefficient C

12 17 30 21 22 13 1 21

nodal demand d

reservoir head H0 pipe roughness CHW

EPANET-FD

Adjoint

Rel. error

1.2 -0.305558467 -0.305558468 1 0.221943309 0.22194331 1.5 0.115063779 0.115063777 9.1228 9.35450703 9.35448105 4.7576 12.5026841 12.5026823 9.221 -4.47928896 -4.47903978 328.08 -1.21584257 -1.2158337 130 -0.809110601 -0.80911228

8.60 × 10−10 6.35 × 10−9 1.79 × 10−8 2.78 × 10−6 1.43 × 10−7 5.56 × 10−5 7.30 × 10−6 2.07 × 10−6

(b) City D, roughness class in breadth; 30 pipes by sensitivity band, no identifiers Sensitivity band

Pipes

|dL/dC| ≥ 3 × 10−2 (12 pipes sampled of 22) [10−2 , 3 × 10−2 ) (12 sampled of 29) clamped at the low-flow floor (all 6; analytic gradient exactly 0)

Worst rel. error Outcome

12 1.26 × 10−5 PASS 12 4.78 × 10−5 PASS 6 |gFD | ≤ 6.45 × 10−7 PASS

Coords Worst rel. error Outcome

emitter coefficient C nodal demand d reservoir head H0 pipe roughness CHW

3 3 1 1

9.43 × 10−6 2.39 × 10−5 3.93 × 10−7 3.04 × 10−6

PASS PASS PASS PASS

Table S4. Gradient audit on degenerate cases: 22 hand-constructed cases, each against its own threshold; o marks the six on the operating network. “0 (structural)” means the correct derivative is analytically zero and the computed one is exactly zero. Degenerate case

Worst

near-clamp step, dead-end pipe, resistance gradient near-clamp step, dead-end junction, demand gradient near-clamp step, finite-difference check, 9 coordinates true clamp, clamped low-flow link, resistance gradient must be exactly zero true clamp, dead-end junction, demand gradient true clamp, finite-difference check, 9 coordinates operating network, closed-link neighbourhood: step inside the clamped branch, |∆G| against its noise ceilingo operating network, base load added: implicit vs unrolled cross-checko operating network, base load added: finite differences, micro-switch averagedo operating network, throttle-valve endpoints: demand gradient vs finite differenceso operating network, throttle valves and closed links: resistance gradient must be exactly zeroo implicit adjoint at a zero-emitter node: finite and exactly zero truncated unrolling at a zero-emitter node: finite and exactly zero control: finite differences at a live emitter coordinate zero-demand junction, demand gradient two-reservoir head gradient against finite differences head-shift identity direct reservoir-head term GPU against CPU on the same batch, synthetic 0009 operating network, GPU against CPU on the same batcho the tensor library’s own finite-difference gradient check, synthetic 0015 the tensor library’s own finite-difference gradient check, synthetic 0006

Threshold Fraction Outcome

0 (structural) 2.08 × 10−11 3.57 × 10−11 0 (structural)

exactly 0 – PASS 1.00 × 10−6 < 0.001 PASS 1.00 × 10−6 < 0.001 PASS exactly 0 – PASS

4.87 × 10−8 2.49 × 10−10 1.74 × 10−1

1.00 × 10−6 0.049 PASS 1.00 × 10−6 < 0.001 PASS 1.00 × 100 0.174 PASS

2.73 × 10−7

1.00 × 10−4

0.003 PASS

1.22 × 10−5

1.00 × 10−3

0.012 PASS

5.75 × 10−7

1.00 × 10−6

0.575 PASS

0 (structural)

exactly 0

– PASS

0 (structural) 0 (structural)

exactly 0 exactly 0

– PASS – PASS

8.02 × 10−11 1.00 × 10−6 1.56 × 10−9 1.00 × 10−6 4.03 × 10−12 1.00 × 10−6 1.32 × 10−15 1.00 × 10−9 1.11 × 10−16 1.00 × 10−12 3.48 × 10−16 1.00 × 10−8 2.76 × 10−16 1.00 × 10−8 0 (structural) exactly 0 0 (structural)

exactly 0

< 0.001 0.002 < 0.001 < 0.001 < 0.001 < 0.001 < 0.001 –

PASS PASS PASS PASS PASS PASS PASS PASS

– PASS

Table S5. The regression suite, item by item: 53 numerical checks, all passing, plus one packaging guard, clean over this submission package (Note 9). o operating-network item; ‡ two networks, City D and one public; § deliberately seeded mutants that must fail. Item

Measured (every threshold met; ft and cfs)

Steady-state alignment against the reference library 23 synthetic networks worst max |∆H| = 2.84 × 10−14 City Do ∆H 1.42 × 10−14 , ∆Q 1.78 × 10−15 , ∆Qcl 4.60 × 10−8 o City D (emitter) ∆H 1.42 × 10−14 , ∆Q 3.55 × 10−15 , ∆Qcl 4.65 × 10−8 , ∆e 2.78 × 10−16 Snapshot replay replay EXA4 replay EXA6 replay City Ho replay ky3 replay ky5

∆H 5.68 × 10−14 , ∆Q 8.88 × 10−16 ∆H 0, ∆Q 8.88 × 10−16 ∆H 1.42 × 10−14 , ∆Q 7.11 × 10−15 ∆H 0, ∆Q 8.88 × 10−16 ∆H 0, ∆Q 3.55 × 10−15

Autonomous extended-period simulation EPS EXA4 EPS EXA5 EPS EXA6 EPS City Ho EPS ky3 EPS ky5 EPS Anytown

∆H 5.68 × 10−14 , ∆Q 8.88 × 10−16 ∆H 5.68 × 10−14 , ∆Q 8.88 × 10−16 ∆H 0, ∆Q 8.88 × 10−16 ∆H 1.42 × 10−14 , ∆Q 7.11 × 10−15 ∆H 0, ∆Q 8.88 × 10−16 ∆H 0, ∆Q 3.55 × 10−15 ∆H 6.08 × 10−12 , ∆Q 4.75 × 10−12

Physical cross-validation cross-validation EXA4 cross-validation EXA5 cross-validation EXA6 cross-validation City Ho cross-validation ky3 cross-validation ky5 cross-validation City Do

mass 1.14 × 10−6 , dem 6.94 × 10−18 , H–W 1.34 × 10−8 mass 1.62 × 10−6 , dem 6.94 × 10−18 , H–W 6.43 × 10−15 mass 1.28 × 10−6 , dem 5.55 × 10−17 , H–W 3.73 × 10−10 mass 6.17 × 10−7 , dem 2.22 × 10−16 , H–W 3.55 × 10−15 mass 3.07 × 10−9 , dem 5.55 × 10−17 , H–W 5.86 × 10−15 mass 5.99 × 10−6 , dem 1.39 × 10−17 , H–W 1.79 × 10−13 mass 6.76 × 10−7 , dem 4.44 × 10−16 , H–W 1.43 × 10−14

Gradient cross-checks 2 networks × 4 parameter classes the tensor library’s own finite-difference gradient check, synthetic 0009 City D, batch-of-8 vs single-scenario consistencyo Degenerate-case gradient audit 22 degenerate cases, own thresholds Capability parity D–W, Chezy–Manning and FCV vs the library

worst implicit 5.85 × 10−8 , worst unrolled 1.21 × 10−6 , City D included‡ pass demand 0, resistance 0 (exact)

22/22 pass (6 on City D; Table S4) ∆H ≤ 2.842 × 10−14 ft on all three

Darcy–Weisbach gradient adjoint vs central differences, Balerma

worst 6.142 × 10−7

Symmetry guard max |A − A⊤ | per Newton round

86/86 bitwise symmetric; 4/4 mutants red§

Status-schedule parity status schedule and PRV subsequences

0 disagreements; 3/3 mutants red

Bit-level reference comparison L-TOWN bit for bit; float32 refused

785/785 heads, 909/909 flows

Packaging guard tracked contents, paths, commit message

42variants; 0 hits over every file, path and compiled 13 name PDF of this submission package (Note 9)

Table S6. L-TOWN cost against batch size on two RTX 5090 nodes (minimum–maximum), double precision. (a) Wall time per scenario with the two-factor decomposition F1 (adjoint onto GPU) and F2 (sparse versus dense). (b) Peak device memory; fwd, forward only; f+b, forward plus backward. (a) Wall time per scenario (ms); forward and forward-plus-backward forward B

dense

forward + backward

sparse

d÷s

CPU adj.

GPU dense

GPU sparse

F1 F2

1 85.98–100.51 56.97–79.01 1.27–1.51 420.06–422.51 152.41–152.90 107.80–109.03 2.75–2.77 1.40–1.42 3.88–3.90 8 20.57–22.98 7.71–10.76 2.13–2.67 349.51–350.87 47.07–47.13 15.18–15.22 7.42–7.44 3.09–3.11 22.96–23.12 64 5.93–6.19 1.38–2.01 3.08–4.31 334.26–335.33 12.01–12.08 3.195–3.205 27.76–27.83 3.75–3.78 104.29–104.95 256 3.78–3.85 0.69–1.03 3.75–5.48 332.90–334.39 6.642–6.670 1.895–1.909 50.12–50.13 3.49–3.50 175.19–175.66 512 3.48–3.50 0.57–0.85 4.12–6.10 333.40–333.49 5.816–5.824 1.671–1.679 57.26–57.33 3.46–3.49 198.52–199.60 1024 3.33–3.35 0.52–0.77 4.36–6.42 OOM OOM 1.612–1.617 – – 198.06–201.23 (b) Peak device memory (MiB), one fresh process per cell B fwd dense fwd sparse f+b dense f+b sparse dense ÷ sparse 1 8 64 256 512 1024

112 292 1 626 7 338 14 564 29 068

46 48 124 296 590 1 212

194 300 1 634 7 346 14 572 29 096

214 86 160 332 626 1 286

0.91 3.49 10.21 22.13 23.28 22.63

Table S7. Cost per steady-state frame across the size sweep: 14 public networks plus the two operating networks, replica path against the reference library, single machine. The factor buys a derivative the reference engine does not supply. Network

L replica (ms) reference (ms) replica ÷ ref.

Hanoi 32 34 synthetic 0003 36 52 Fossolo (poly1) 37 58 Pescara 71 99 Net3 97 119 Modena 272 317 EXA4 396 444 Balerma 447 454 City D (operating) 542 554 L-TOWN 785 909 Richmond (standard) 872 957 City H (operating) 921 1038 ky10 935 1061 ky4 964 1158 Net6 3356 3892 BWSN Network 2 12527 14831

0.498 1.052 1.590 2.824 1.831 6.489 8.181 5.473 6.371 36.199 13.982 15.007 24.195 23.036 79.354 478.214

0.0184 0.0294 0.0219 0.0594 0.0715 0.1871 0.2140 0.1181 0.3038 0.7390 0.3449 0.3200 0.5110 0.5339 2.6096 9.5961

27.0 35.8 72.6 47.5 25.6 34.7 38.2 46.3 21.0 49.0 40.5 46.9 47.3 43.1 30.4 49.8

Table S8. Calibration noise ladders (a) and City D robustness arms (b). ic inverse crime: truth generated by the model being inverted; Modena has only such rows, which is why it appears here and not in Results (Supplementary Note 8). (a) Noise ladders; median final training MSE over the seeds of each arm Network

Free pipes σ (ft)

Seeds Median MSE Note 1 3.00 × 10−11 inverse crimeic 1 9.73 × 10−28 inverse crimeic 3 8.85 × 10−4 3 9.83 × 10−3 3 8.85 × 10−2

Hanoi (public) Hanoi (public) Hanoi (public) Hanoi (public) Hanoi (public)

34 34 34 34 34

0 (truth per-pipe) 0 (truth grouped) 0.03 0.1 0.3

Modena (public) Modena (public)

317 0 (truth per-pipe) 317 0 (truth grouped)

1 1

3.25 × 10−6 inverse crimeic 1.96 × 10−5 inverse crimeic

City D (operating) City D (operating) City D (operating) City D (operating) City D (operating)

432 432 432 432 432

1 1 5 5 5

2.98 × 10−6 inverse crimeic 3.39 × 10−6 inverse crimeic 8.43 × 10−4 9.40 × 10−3 8.48 × 10−2

0 (truth per-pipe) 0 (truth grouped) 0.03 0.1 0.3

(b) City D robustness arms at σ = 0.1 ft (median final training MSE) Arm demand model error: truth × U [0.85, 1.15], inverted on nominal demand paired control, demand known sensor bias 5 % (σb = 0.05) multi-start: 8 Latin-hypercube starts, C ∈ [80, 145]

Runs

MSE

9.51 × 10−3 (worst 1.62 × 10−2 )

4 9.11 × 10−3 5 9.86 × 10−3 8 9.34 × 10−3 –9.42 × 10−3 (1.0087× spread)

Table S9. Sensor augmentation, all virtual (designed and verified in simulation; none installed). (a) City D: recovery of the 149 unobservable pipes when k sensors are added to the 40 designated, both objectives, against reselection from scratch. (b) The same prescription on L-TOWN and Hanoi (Supplementary Note 10). (a) City D: 149 pipes unobservable under the 40 designated sensors; 541 candidates, 25 frames, census threshold, σnoise = 0.1 ft k Sensors Recovered Still unobs. Lost Rank ksub

Design

CRLBsub

S0 (designated) 0 cover +k 5 cover +k 10 cover +k 20 cover +k 40 cover +k 80

40 45 50 60 80 120

0 37 62 95 133 149

149 112 87 54 16 0

0 0 0 0 0 0

110 122 132 144 163 202

20 22 28 33 36 59

742 489 686 656 347 486

D-opt +k D-opt +k D-opt +k D-opt +k D-opt +k

5 10 20 40 80

45 50 60 80 120

21 23 35 52 76

128 126 114 97 73

0 0 0 0 0

118 123 133 157 196

24 28 36 56 79

457 359 222 512 1430

reselect 40+k reselect 40+k reselect 40+k reselect 40+k reselect 40+k reselect 40

5 10 20 40 80 0

45 50 60 80 120 40

59 59 68 73 94 52

128 128 105 97 70 135

38 38 24 21 15 38

128 134 149 166 201 122

– – – – – –

(b) Public networks: recovered of the unobservable pipes at k added (nothing lost at any k under either objective) +40

+80

L-TOWN, 33 fixed, 62 unobservable add, coverage 21 33 46 62 add, D-optimal 2 2 3 7 reselect 33 + k: recovered / lost 6/99 6/99 8/99 13/98

62 10 15/2

Network

Hanoi, 10 fixed, 1 unobservable

+5 +10 +20

Design

add, coverage add, D-optimal

1 1

– –

Table S10. Calibration with the augmented sensor sets, the calibrator, truth, noise and hold-out of Supplementary Note 8 unchanged. (a) City D: error of the restored pipes against their prior, and the RMSE over all informative pipes. (b) Hanoi (Supplementary Note 10). (a) City D, 432 free pipes, one noise seed per cell, 20 % of sensors held out Set

Sensors σ (ft) Inform. Restored Median |∆C| Prior

S0 S0 S0

40 0.03 40 0.1 40 0.3

279 279 279

0 0 0

D-opt +20 D-opt +20 D-opt +20

60 0.03 60 0.1 60 0.3

314 314 314

D-opt +40 D-opt +40 D-opt +40

80 0.03 80 0.1 80 0.3

cover +20 cover +20 cover +20 cover +40 cover +40 cover +40

– – –

Better RMSE, inform.

– – –

24.4 25.6 26.3

35 35 35

11.5 21.9 11.5 21.9 12.2 21.9

21/35 24/35 20/35

22.9 23.2 30.3

331 331 331

52 52 52

12.1 18.3 13.2 18.3 11.9 18.3

26/52 26/52 23/52

22.4 23.6 23.9

60 0.03 60 0.1 60 0.3

374 374 374

95 95 95

14.8 19.1 15.5 19.1 19.7 19.1

45/95 44/95 40/95

23.5 24.3 29.5

80 0.03 80 0.1 80 0.3

412 412 412

133 133 133

16.9 19.5 58/133 18.4 19.5 56/133 22.0 19.5 49/133

23.0 24.2 30.0

(b) Hanoi, 34 free pipes, median RMSE over the informative pipes across three noise seeds Set S0 (10 random) D-opt +5 D-opt +10 D-opt +20

Sensors Unobservable σ = 0.03 σ = 0.1 10 15 20 30

1 0 0 0

18.3 9.9 7.8 6.2

18.7 11.0 10.2 7.0

σ = 0.3 19.4 14.4 17.0 9.8

Table S11. Leak search with the augmented sensor sets. (a) The work-order search on City D with the designated set, the six augmented sets and the control experiments; † an added sensor sits on L2 ’s junction, ‡ the main-text 40-sensor set, sharing 22 junctions with S0 . (b) The L-TOWN inversion with its ten augmented sets (Supplementary Note 10). (a) Work-order search with the augmented sets, City D: 49 candidates, σ = 0.1 ft, one noise realisation for every set Sensors Hits Which Error Coh. median > 0.999 Rival L2 Rival L3

Sensor set ‡

40 of main-text Section 3.4 S0 (designated) cover +20† D-opt +20† cover +40† D-opt +40† cover +80† D-opt +80†

40 40 60 60 80 80 120 120

0/3 none 0/3 none 1/3 L2 1/3 L2 1/3 L2 1/3 L2 1/3 L2 1/3 L2

– – 2.3 % 3.7 % 2.3 % 1.5 % 2.1 % 1.1 %

0.934 0.928 0.796 0.897 0.789 0.846 0.806 0.791

51 48 22 34 18 27 20 23

1.00000 1.00000 0.99898 0.99893 0.99893 0.99900 0.99909 0.99934

0.9993 0.9994 0.9995 0.9995 0.9996 0.9996 0.9997 0.9997

cover +20 minus the L2 -junction sensor D-opt +20 minus that sensor S0 plus that sensor alone S0 plus an L3 -junction sensor random +20, draw 1 random +20, draw 2 random +20, draw 3 random +20, draw 4 random +20, draw 5 random +40

59 59 41 41 60 60 60 60 60 80

0/3 none 1/3 L3 1/3 L2 0/3 none 1/3 L3 1/3 L3 0/3 none 1/3 L3 0/3 none 0/3 none

– 7.0 % 3.0 % – 3.9 % 4.2 % – 4.2 % – –

0.796 0.902 0.928 0.920 0.933 0.919 0.930 0.913 0.940 0.913

27 39 44 47 38 37 37 35 46 33

1.00000 1.00000 0.99905 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000

0.9995 0.9995 0.9994 0.9995 0.9995 0.9995 0.9995 0.9995 0.9995 0.9997

(b) L-TOWN inversion with the augmented sets: 60 candidates, 1,770 pairs (275 orthogonal); hits of 3 in the noiseless and the σ = 0.1 ft group Sensor set S0 (33 of main-text Section 3.4) D-opt +5 D-opt +10 D-opt +20 D-opt +40 D-opt +80 cover +5 cover +10 cover +20 cover +40 cover +80

Sensors Clean Noisy Med. all Med. zone > 0.999 33 38 43 53 73 113 38 43 53 73 113

1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3

0.8258 0.8248 0.8260 0.8262 0.8281 0.8243 0.8252 0.8304 0.8200 0.8132 0.8161

0.8713 0.8678 0.8716 0.8704 0.8682 0.8662 0.8737 0.8807 0.8627 0.8567 0.8618

Rival A

Rival B

42 0.999231 40 0.998821 40 0.998797 41 0.998912 43 0.998878 44 0.998560 41 0.999204 37 0.999289 38 0.999362 31 0.999175 30 0.999057

0.999965 0.999965 0.999966 0.999967 0.999966 0.999917 0.999966 0.999966 0.999860 0.999849 0.999859

Table S12. Coherence-driven placement on City D against random controls, the same test re-scored on the coarse stage L-TOWN’s inversion shares, and the L-TOWN candidate-pool bound. (a) Fisher test and rival-coherence gap by leak and inversion stage. (b) L2 ’s matched-budget exact test and the random group’s 95 % Clopper–Pearson interval. (c) Separable positions and the oracle-pair-greedy floor for L-TOWN’s two never-recovered leaks (Supplementary Note 10). (a) City D: 14 designed sets (8 coherence-min., 6 identifiability-driven) vs 40 random draws, by leak / stage Leak Stage

Design Random

L1 L1 L2 L2 L3 L3

0/14 6/14 12/14 13/14 3/14 5/14

final stage1 final stage1 final stage1

0/40 2/40 3/40 8/40 7/40 20/40

Fisher p Gap, recovered 1.00 × 100 2.39 × 10−3 1.05 × 10−7 2.10 × 10−6 5.12 × 10−1 8.92 × 10−1

Gap, lost Separated

– 0.667 no 0.55 0.624 no 3.12e-05 2.8e-05 yes 3.01e-07 3.01e-07 yes 0.000327 0.000606 no 0.000327 0.00059 no

(b) L2 , same-budget exact permutation test and the random group’s 95% Clopper–Pearson CI Budget

Random draws Random hits

p (exact) Random rate, 95% CI

k = 20 (worst-pair variant, +20) k = 40 (greedy coherence sequence, +40)

20 20

0/20 4.76 × 10−2 0.000–0.168 3/20 1.90 × 10−1 0.032–0.379

pooled random rate, L1 pooled random rate, L2 pooled random rate, L3

40 40 40

0/40 3/40 7/40

– 0.000–0.088 – 0.016–0.204 – 0.073–0.328

Rival coh. at S0 Separable of Pool Oracle floor

rank 19 of 60 rank 32 of 60

0.999965 0.999231

1 14

746 746

0.996775 0.990741

Rival after 0.999400 0.997689

Table S13. The exact-input control: reservoir heads and valve settings read from the input text instead of through a unit round-trip (Supplementary Note 1). (a) Every network the control changes: the metre round-trip path against the input-text path, over all frames. (b) The negative control: the same code path on eight networks whose fields it leaves untouched. default input path Network

Frames Outcome

(a) every network the control changes Net3 183 fail BWSN Network 1 207 fail BWSN Network 2 32 fail Net6 609 fail Anytown 9 pass

max|∆H| (ft)

fields

max|∆Q| res./valve (cfs) changed

exact-input control max|∆H| (ft)

max|∆Q| (cfs)

2.15 × 10−6 5.68 × 10−7 6.09 × 10−6 7.05 × 10−6 6.84 × 100 1.41 × 10−5 2.30 × 10−5 4.74 × 10−6 6.08 × 10−12 4.75 × 10−12

1/0 1/3 2/1 1/0 2/0

0 0 0 0 0

3.55 × 10−15 8.88 × 10−16 7.11 × 10−15 7.11 × 10−15 1.78 × 10−15

(b) controls: the same code path on eight networks it leaves alone EXA5 52 pass 5.68 × 10−14 8.88 × 10−16 City D 25 pass 1.42 × 10−14 1.78 × 10−15 City H 25 pass 1.42 × 10−14 7.11 × 10−15 C-Town (BATADAL) 850 pass 5.68 × 10−14 8.88 × 10−16 D-Town 819 pass 5.68 × 10−14 1.78 × 10−15 L-TOWN 2,031 pass 2.84 × 10−14 1.11 × 10−16 Richmond (skeleton) 91 pass 1.14 × 10−13 2.22 × 10−16 Richmond (standard) 55 pass 1.14 × 10−13 2.22 × 10−16

0/0 0/0 0/0 0/0 0/0 0/0 0/0 0/0

5.68 × 10−14 1.42 × 10−14 1.42 × 10−14 5.68 × 10−14 5.68 × 10−14 2.84 × 10−14 1.14 × 10−13 1.14 × 10−13

8.88 × 10−16 1.78 × 10−15 7.11 × 10−15 8.88 × 10−16 1.78 × 10−15 1.11 × 10−16 2.22 × 10−16 2.22 × 10−16

Table S14. Cluster-level leak diagnosability: the hit-rate against inspection-radius trade-off (Supplementary Note 11). (a), (b) Full threshold sweeps with the geometry-only, coherence-only and equal-budget single-point controls. (c) The inspection burden of the three highest-ranked clusters. (d) Size-matched random grouping. τ

clusters

mean mean cluster level geometry coherence equal-budget radius district top-1 top-3 only only single point (m) (km) top-3 top-3 top-3

(a) operating network, 49 work-order candidates, 40 sensors, σ = 0.1 ft 0 1 9043 83.92 3/3 3/3 3/3 3/3 0.5 4 1902 20.98 1/3 2/3 2/3 3/3 0.8 5 1506 16.78 1/3 2/3 1/3 3/3 0.9 6 1427 13.99 1/3 1/3 2/3 1/3 0.95 8 1086 10.49 1/3 2/3 2/3 2/3 0.98 12 788 6.99 1/3 2/3 1/3 2/3 0.99 13 729 6.46 1/3 2/3 1/3 2/3 0.995 17 379 4.94 1/3 2/3 1/3 2/3 0.999 25 179 3.36 1/3 2/3 1/3 2/3 0.9995 28 140 3.00 1/3 1/3 1/3 1/3 0.9999 30 114 2.80 1/3 1/3 1/3 1/3 >1 49 0 1.71 0/3 0/3 0/3 0/3

3/3 2/3 2/3 2/3 0/3 0/3 0/3 0/3 0/3 0/3 0/3 0/3

(b) L-TOWN, 60 candidates, 33 sensors, 256 scenarios, σ = 0.1 ft 0 1 1575 43.16 3/3 3/3 3/3 3/3 0.5 2 934 21.58 3/3 3/3 3/3 3/3 0.8 4 542 10.79 2/3 3/3 3/3 3/3 0.9 6 399 7.19 2/3 3/3 3/3 3/3 0.95 9 295 4.80 1/3 2/3 2/3 3/3 0.98 15 185 2.88 1/3 2/3 3/3 2/3 0.99 17 160 2.54 1/3 3/3 3/3 3/3 0.995 22 98 1.96 1/3 3/3 2/3 3/3 0.999 36 40 1.20 1/3 2/3 1/3 2/3 0.9995 41 30 1.05 0/3 1/3 1/3 1/3 0.9999 53 7 0.81 0/3 1/3 1/3 1/3 >1 60 0 0.72 0/3 1/3 1/3 1/3

– – – – – – – – – – – –

(c) what a crew would actually walk: the three highest-ranked clusters, operating network τ clusters candidates pipe in share of largest of the mean radius mean district in top 3 top 3 (km) network three radii (m) all clusters (m) all clusters (km) 0.9 0.95 0.99 0.995 0.999

6 8 13 17 25

23 14 10 8 7

34.70 18.53 10.43 8.66 7.42

41.3 % 22.1 % 12.4 % 10.3 % 8.8 %

2344 2344 1041 750 750

1427 1086 729 379 179

13.99 10.49 6.46 4.94 3.36

(d) size-matched random grouping: same cluster size distribution, topology and coherence destroyed Network

operating network 0.9 0.99 0.999 L-TOWN 0.9 0.99 0.999

design radius random radius ≥ design top-3 (m) top-3 (m) 1/3 2/3 2/3 3/3 3/3 2/3

1427 729 179 399 160 40

2.10/3 0.65/3 0.75/3 2.30/3 1.70/3 1.50/3

4415 3544 1700 1223 698 298

20/20 1/20 2/20 8/20 3/20 8/20

1.000 0.095 0.143 0.429 0.190 0.429

Table S15. Sensor layouts scored on a reference subspace fixed before any design is chosen (Supplementary Note 12). (a) All twelve budget-by-noise cells against same-budget random additions. (b) The main cell on five rulers, including the discarded directions. (c) The cost of selecting a layout by head misfit. (d) Design-time prediction against achieved error. (a) every cell: the designed additions against random additions from the same pool at the same budget Budget σ (ft)

+20

0.03

+20

0.1

+20

0.3

+40

0.03

+40

0.1

+40

0.3

objective

R identifiable-subspace roughness RMSE

coverage 20 D-optimal 20 coverage 100 D-optimal 100 coverage 20 D-optimal 20 coverage 20 D-optimal 20 coverage 20 D-optimal 20 coverage 20 D-optimal 20

design

best median

worst ≥ design

15.87 17.50 17.45 18.04 18.13 21.85 14.47 15.64 16.24 18.02 19.50 19.99

18.58 18.58 19.70 19.70 21.87 21.87 16.60 16.60 18.04 18.04 20.06 20.06

21.18 21.18 23.81 23.81 26.56 26.56 21.13 21.13 22.01 22.01 24.45 24.45

20.48 20.48 21.73 21.73 24.15 24.15 19.22 19.22 20.11 20.11 22.48 22.48

random

0/20 0/20 0/100 0/100 0/20 0/20 0/20 0/20 0/20 0/20 0/20 0/20

0.0476 0.0476 0.0099 0.0099 0.0476 0.0476 0.0476 0.0476 0.0476 0.0476 0.0476 0.0476

(b) the same cell (+20, σ = 0.1 ft, R = 100, coverage objective) read on five different rulers Ruler

design random median random ≥ design

identifiable subspace, γ = 2 (k = 98) identifiable subspace, γ = 3 (k = 70) identifiable subspace, γ = 10 (k = 29) all 432 pipes, no projection reference null space (432 − 98 directions)

17.45 15.83 9.86 25.10 26.93

21.73 21.83 20.05 25.79 26.86

0/100 0/100 0/100 0/100 60/100

p 0.0099 0.0099 0.0099 0.0099 0.6040

(c) the price of choosing a layout by head misfit: pick the best of 100 random additions by each head metric Selection criterion

training head MSE held-out-frame head RMSE held-out-sensor head RMSE

its ident. RMSE

rank percentile best in pool

21.05 25/100 21.80 56/100 20.65 10/100

24 % 55 % 9%

regret

19.70 19.70 19.70

+1.36 +2.10 +0.95

(d) how well the design-time linear-Gaussian prediction orders the achieved error Pool +20 only (n = 100) +40 only (n = 20) both budgets pooled (n = 120)

Spearman ρ

p predicted spread achieved spread

0.38 1.3e−04 0.69 6.0e−04 0.55 6.3e−11

5.6 % 4.7 % 8.2 %

19.0 % 19.6 % 26.9 %

Table S16. Optimiser enhancements and the wall-clock reading on City D (Supplementary Note 13). (a) Every arm at a 2,000-call evaluation budget. (b) What batching multiple starts buys per model call. (c) The derivative-free baselines on the same node, with the eight-way-parallel boundary the fairness protocol promises them. (a) City D, 432 free pipes, evaluation budget 2,000 model calls, medians over noise seeds Arm

seeds

reference gradient configuration first order, no preconditioner first order, Schur-complement diagonal first order, exact Gauss–Newton diagonal first order, 8 starts batched first order, 8 starts + Schur first order, 16 starts + Schur Levenberg–Marquardt tail LM tail + Schur LM tail + 8 starts + Schur

4 10 10 10 10 10 10 5 4 4

calls training loss roughness RMSE wall (s) 200 1976 1976 1976 1976 1976 1976 1997 1998 1997

9.109e-03 9.492e-03 9.541e-03 9.459e-03 9.575e-03 9.542e-03 9.580e-03 9.459e-03 9.565e-03 9.052e-03

11.32 17.71 12.40 16.04 15.00 13.63 14.09 15.17 11.37 13.69

631 223 221 232 139 141 127 2331 2326 2383

(b) what batching many starts buys: it is not a factor of B starts B s per step s per model call 1 2 4 8 16 32

0.1123 0.1507 0.2326 0.3773 0.6861 1.2770

speed-up per call

0.1123 0.0754 0.0581 0.0472 0.0429 0.0399

1.00× 1.49× 1.93× 2.38× 2.62× 2.81×

(c) the derivative-free baselines on the same node and the same budget: the wall-clock reading Algorithm differential evolution particle swarm CMA-ES DE → LM hybrid differential evolution, 10× budget

seeds

calls training loss wall (s)

5 1984 5 1984 5 1980 1 1997 3 20000

1.449e-01 8.000e-02 2.488e-01 9.400e-03 1.433e-02

156 156 156 4756 1691

wall / 8 (s) 20 20 20 595 211

References Brandon Amos and J. Zico Kolter. OptNet: Differentiable optimization as a layer in neural networks. In Proceedings of the 34th International Conference on Machine Learning (ICML), pages 136–145, 2017. arXiv:1703.00443. Ernesto Arandia and Bradley J. Eck. An r package for epanet simulations. Environmental Modelling & Software, 107:59–63, 2018. doi: 10.1016/j.envsoft.2018.05.016. Inaam Ashraf, Janine Strotherm, Luca Hermes, and Barbara Hammer. Physics-informed graph neural networks for water distribution systems. In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI-24), volume 38, pages 21905–21913, 2024. doi: 10.1609/aaai.v38i20.30192. arXiv:2403.18570. Inaam Ashraf, André Artelt, and Barbara Hammer. Scalable and robust physics-informed graph neural networks for water distribution systems. In 2025 International Joint Conference on Neural Networks (IJCNN), pages 1–8. IEEE, 2025. doi: 10.1109/ijcnn64981.2025.11229349. Menglong Cheng, Juan Li, Chunyue Wang, Chaoxiong Ye, and Zheng Chang. Graph Laplace regularizationbased pressure sensor placement strategy for leak localization in the water distribution networks under joint hydraulic and topological feature spaces. Water Research, 257:121666, 2024. doi: 10.1016/j.watres.2024. 121666. Ivo Daniel, Jorge Pesantez, Simon Letzgus, Mohammad Ali Khaksar Fasaee, Faisal Alghamdi, Emily Berglund, G. Mahinthakumar, and Andrea Cominola. A sequential pressure-based algorithm for data-driven leakage identification and model-based localization in water distribution networks. Journal of Water Resources Planning and Management, 148(6):04022025, 2022. doi: 10.1061/(ASCE)WR.1943-5452. 0001535. Kegong Diao, Zhengji Wang, Gregor Burger, Chien-Hsun Chen, Wolfgang Rauch, and Yuwen Zhou. Speedup of water distribution simulation by domain decomposition. Environmental Modelling & Software, 52: 253–263, 2014. doi: 10.1016/j.envsoft.2013.09.025. K. Du, J. Yu, F. Zheng, Z. Kapelan, and D. Savic. Decoupling elevation errors from pipe roughness calibration in hydraulic network models. Water Research, 290:125058, 2026. doi: 10.1016/j.watres.2025.125058. Bradley J. Eck. An r package for reading epanet files. Environmental Modelling & Software, 84:149–154, 2016. doi: 10.1016/j.envsoft.2016.06.027. Sylvan Elhay and Angus R. Simpson. Dealing with zero flows in solving the nonlinear equations for water distribution systems. Journal of Hydraulic Engineering, 137(10):1216–1224, 2011. doi: 10.1061/(ASCE) HY.1943-7900.0000411. Alan George and Joseph W. H. Liu. Computer Solution of Large Sparse Positive Definite Systems. PrenticeHall, 1981. George Germanopoulos. A technical note on the inclusion of pressure dependent demand and leakage terms in water supply network models. Civil Engineering Systems, 2(3):171–179, 1985. doi: 10.1080/ 02630258508970401. S. Guo, K. Xin, T. Tao, and H. Yan. A deep-level decomposed model to accelerate hydraulic simulations in large water distribution networks. Water Research, 266:122318, 2024. doi: 10.1016/j.watres.2024.122318.

Gergely Hajgátó, Bálint Gyires-Tóth, and György Paál. Reconstructing nodal pressures in water distribution systems with graph neural networks. arXiv:2104.13619, 2021. Zoran S. Kapelan, Dragan A. Savic, and Godfrey A. Walters. Optimal sampling design methodologies for water distribution model calibration. Journal of Hydraulic Engineering, 131(3):190–200, 2005. doi: 10.1061/(ASCE)0733-9429(2005)131:3(190). Bulat Kerimov and Riccardo Taormina. Towards transferable metamodels for water distribution systems with edge-based graph neural networks. Water Research, 261:121933, 2024. doi: 10.1016/j.watres.2024.121933. Bulat Kerimov, Roberto Bentivoglio, Alexander Garzón, Elvin Isufi, Franz Tscheikner-Gratl, David B. Steffelbauer, and Riccardo Taormina. Assessing the performances and transferability of graph neural network metamodels for water distribution systems. Journal of Hydroinformatics, 25(6):2223–2234, 2023. doi: 10.2166/hydro.2023.031. Katherine A. Klise, David Hart, Dylan Moriarty, Michael Bynum, Regan Murray, Jonathan Burkhardt, and Terranna Haxton. Water network tool for resilience (WNTR) user manual. Technical report, U.S. Environmental Protection Agency, 2017. Andreas Krause, Jure Leskovec, Carlos Guestrin, Jeanne VanBriesen, and Christos Faloutsos. Efficient sensor placement optimization for securing large water distribution networks. Journal of Water Resources Planning and Management, 134(6):516–526, 2008. doi: 10.1061/(ASCE)0733-9496(2008)134:6(516). Joseph W. H. Liu. Modification of the minimum-degree algorithm by multiple elimination. ACM Transactions on Mathematical Software, 11(2):141–153, 1985. doi: 10.1145/214392.214398. H.R. Maier, Z. Kapelan, J. Kasprzyk, J. Kollat, L.S. Matott, M.C. Cunha, G.C. Dandy, M.S. Gibbs, E. Keedwell, A. Marchi, A. Ostfeld, D. Savic, D.P. Solomatine, J.A. Vrugt, A.C. Zecchin, B.S. Minsker, E.J. Barbour, G. Kuczera, F. Pasha, A. Castelletti, M. Giuliani, and P.M. Reed. Evolutionary algorithms and other metaheuristics in water resources: Current status, research challenges and future directions. Environmental Modelling & Software, 62:271–299, 2014. doi: 10.1016/j.envsoft.2014.09.013. Tianwei Mu, Yue Wang, Mingzhe Yuan, Wenhong Wang, Qing Luo, Min Xiao, Jun Li, and Hui Yang. HydroGrad: an exactly differentiable global gradient algorithm (Python package dgga). Software, Zenodo, 2026. Version and DOI to be assigned at the archived release accompanying submission. Nikolaj T. Mücke, Prerna Pandey, Shashi Jain, Sander M. Bohté, and Cornelis W. Oosterlee. A probabilistic digital twin for leak localization in water distribution networks using generative deep learning. Sensors, 23 (13):6179, 2023. doi: 10.3390/s23136179. Ji-Ye Park, Kwang-Ju Kim, Minhyuk Jeung, In-Su Jang, Jung-Won Yu, Mi-Seon Kang, Hyun-Su Bae, Changyoon Jeong, and Sang-Soo Baek. Digital-twin tool for a drinking water distribution system using augmented reality and epanet. Environmental Modelling & Software, 197:106829, 2026. doi: 10.1016/j. envsoft.2025.106829. Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems (NeurIPS), pages 8024–8035, 2019.

Ramon Pérez, Vicenç Puig, Josep Pascual, Joseba Quevedo, Edson Landeros, and Antonio Peralta. Methodology for leakage isolation using pressure sensitivity analysis in water distribution networks. Control Engineering Practice, 19(10):1157–1167, 2011. doi: 10.1016/j.conengprac.2011.06.004. Olivier Piller, Sylvan Elhay, Jochen Deuerlein, and Angus R. Simpson. Local sensitivity of pressure-driven modeling and demand-driven modeling steady-state solutions to variations in parameters. Journal of Water Resources Planning and Management, 143(2):04016074, 2017. doi: 10.1061/(ASCE)WR.1943-5452. 0000729. See also the erratum, ibid. 143(8), 2017, doi 10.1061/(ASCE)WR.1943-5452.0000813. Ranko S. Pudar and James A. Liggett. Leaks in pipe networks. Journal of Hydraulic Engineering, 118(7): 1031–1046, 1992. doi: 10.1061/(ASCE)0733-9429(1992)118:7(1031). Lewis A. Rossman, Hyoungmin Woo, Michael Tryby, Feng Shang, Robert Janke, and Terranna Haxton. EPANET 2.2 user manual. Technical report, U.S. Environmental Protection Agency, Cincinnati, OH, 2020. Water Infrastructure Division, Center for Environmental Solutions and Emergency Response. Lina Sela, Elad Salomons, and Mashor Housh. Plugin prototyping for the epanet software. Environmental Modelling & Software, 119:49–56, 2019. doi: 10.1016/j.envsoft.2019.05.010. Jack Sherman and Winifred J. Morrison. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. The Annals of Mathematical Statistics, 21(1):124–127, 1950. doi: 10.1214/aoms/1177729893. Ana Luís Sousa, Angela F. Brochado, Eugénio Rocha, and António Andrade-Campos. A state-of-the-art review of scientific, industrial, and commercial developments in water leakage management for water supply systems. Water Research, 299:125853, 2026. doi: 10.1016/j.watres.2026.125853. Janine Strotherm, Luca Hermes, Inaam Ashraf, and Barbara Hammer. Go with the flow: Leveraging physicsinformed gradients to solve real-world problems in water distribution systems. In Lecture Notes in Computer Science, pages 41–59. Springer, 2025. doi: 10.1007/978-3-032-06129-4_3. Ezio Todini and S. Pilati. A gradient algorithm for the analysis of pipe networks. In Bryan Coulbeck and Chun-Hou Orr, editors, Computer Applications in Water Supply, Volume 1: Systems Analysis and Simulation. Research Studies Press, 1988. Ezio Todini and Lewis A. Rossman. Unified framework for deriving simultaneous equation algorithms for water distribution networks. Journal of Hydraulic Engineering, 139(5):511–526, 2013. doi: 10.1061/ (ASCE)HY.1943-7900.0000703. Joel A. Tropp and Anna C. Gilbert. Signal recovery from random measurements via orthogonal matching pursuit. IEEE Transactions on Information Theory, 53(12):4655–4666, 2007. doi: 10.1109/tit.2007.909108. Huy Truong, Andrés Tello, Alexander Lazovik, and Victoria Degeler. Graph neural networks for pressure estimation in water distribution systems. Water Resources Research, 60, 2024. doi: 10.1029/2023WR036741. Huy Truong, Andrés Tello, Alexander Lazovik, and Victoria Degeler. DiTEC-WDN: A large-scale dataset of hydraulic scenarios across multiple water distribution networks. Scientific Data, 12(1):1733, 2025. doi: 10.1038/s41597-025-06026-0. arXiv:2503.17167. Aly-Joy Ulusoy, Herman A. Mahmoud, Filippo Pecci, Edward C. Keedwell, and Ivan Stoianov. Bi-objective design-for-control for improving the pressure management and resilience of water distribution networks. Water Research, 222:118914, 2022. doi: 10.1016/j.watres.2022.118914. 54

Stelios G. Vrachimis, Demetrios G. Eliades, Riccardo Taormina, Avi Ostfeld, Zoran Kapelan, Shuming Liu, Marios Kyriakou, Pavlos Pavlou, Mengning Qiu, and Marios M. Polycarpou. BattLeDIM: Battle of the leakage detection and isolation methods. In Proc. 2nd Int. CCWI/WDSA Joint Conference, 2020. Dataset: Zenodo record 4017659, https://doi.org/10.5281/zenodo.4017659, CC BY 4.0. Stelios G. Vrachimis, Demetrios G. Eliades, Riccardo Taormina, Zoran Kapelan, Avi Ostfeld, Shuming Liu, Marios Kyriakou, Pavlos Pavlou, Mengning Qiu, and Marios M. Polycarpou. Battle of the leakage detection and isolation methods. Journal of Water Resources Planning and Management, 148(12):04022068, 2022. doi: 10.1061/(ASCE)WR.1943-5452.0001601. Janet M. Wagner, Uri Shamir, and David H. Marks. Water distribution reliability: Simulation methods. Journal of Water Resources Planning and Management, 114(3):276–294, 1988. doi: 10.1061/(ASCE) 0733-9496(1988)114:3(276). Not to be confused with the companion paper in the same issue, “Water Distribution Reliability: Analytical Methods“, 114(3):253–275. Zhiyu Zhang, Wenchong Tian, Zhenliang Liao, and Zhiguo Yuan. Differentiable neural network-based models enable gradient-based optimization for model predictive control of urban drainage networks. Water Research, 291:125188, 2026. doi: 10.1016/j.watres.2025.125188. X. Zhou, X. Wan, S. Liu, K. Su, W. Wang, and R. Farmani. An all-purpose method for optimal pressure sensor placement in water distribution networks based on graph signal analysis. Water Research, 266:122354, 2024. doi: 10.1016/j.watres.2024.122354.

Record · ID 668011 · SHA-256 4af699863d226f12
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.