← Back to the map

How it works

The end-to-end methodology: live discharge telemetry routed through the Indus basin's barrage network in one dimension, coupled to precomputed 2D LISFLOOD-FP physics — extent, depth, channel stress, pondage and exposure, validated against radar observation.

01What it does

The platform models flooding across Pakistan's Indus river system as one connected, regulated machine. It opens on the country's current state: Flood Forecasting Division discharge readings, scraped every 15 minutes from 27 stations, are routed through the full water network — six rivers, five confluences, eight link canals and twelve regulated barrages — so the first map you see is where today's water actually is and where it is going. From there an operator intervenes: raise a headwork discharge to a hypothetical flood, or cap the release at any barrage and watch held water pond behind the gates while the entire basin re-routes in under a second.

Two layers make that possible. Layer one is a live one-dimensional network router that traverses the river-and-canal graph with Muskingum lag and attenuation, sums flows at confluences, applies each barrage's release cap, and flags spill wherever routed discharge exceeds the published bankfull capacity. Layer two is a library of two-dimensional LISFLOOD-FP shallow-water simulations — real physics, precomputed offline as 20-day transient runs across seven discharge bands for each river valley — that turns the router's per-node outflows into inundation extent and depth on the ground. The coupling is what keeps the platform interactive: the physics is genuine but already solved, so moving a slider costs a re-route and a tile fetch, not a model run.

This is the difference from a static global hazard map. Products like the JRC return-period layers give a single fixed picture and cannot answer "what happens at this discharge, today, here — and what changes if Trimmu holds?" The platform is interactive, scenario-driven and locally conditioned, and it is validated against observation: the September 2014 Chenab flood (Marala peak 562,848 cusecs, ~15,930 m³/s) is reproduced inside the observed 10–20 km corridor, and the 2025 super-flood is scored cell-by-cell against a Sentinel-1 radar mask (POD ≈ 0.78). Section 14 reports every number, including the unflattering ones.

outflow = min(inflow, cap); held = inflow − outflow; spill = max(Q − Q_bankfull, 0)
Regulation and overtopping in the network router: each of the 12 regulated barrages passes at most its user-set release cap (cap = 0 holds everything), and routed discharge above a node's FFD bankfull capacity is flagged as spill.
Two-layer flood platform architecture Flowchart: live FFD telemetry, headwork discharges and barrage release caps feed a one-dimensional network router (Muskingum reaches, confluence sums, twelve barrage caps, spill versus bankfull, sub-second re-route). Offline, FABDEM terrain with a JRC-burned channel and WorldCover roughness, Gaussian flood waves and paired barrage-wall runs feed precomputed LISFLOOD-FP valley libraries stored as maximum-depth and pondage-pool COGs. Ratio-banded serving joins the two layers and drives the map layers — extent, depth, channel stress, pondage pools, exposure — validated against Sentinel-1 2025, the observed 2014 extent and FFD hydrographs. LIVE INPUTS FFD live telemetry 27 stations · every 15 min Headwork discharge 5 rim stations · sliders Barrage release caps 12 barrages · cap 0 = hold LAYER 1 · LIVE 1-D network router Muskingum lag + attenuation 6 rivers · 5 confluences · 8 link canals 12 barrage caps · spill vs bankfull re-routes the whole basin < 1 s OFFLINE INPUTS FABDEM terrain 120 m JRC channel · Manning roughness Gaussian flood waves 7 bands · 50k–1M cusecs Barrage wall runs paired free vs walled LAYER 2 · PRECOMPUTED LISFLOOD-FP 8.0.3 2-D shallow water · run offline 5 valley libraries · 20-day transients max-depth COGs · pondage pools 36 pondage libs · ≤2% mass closure SERVING Ratio-banded serving latitude slabs where outflow / baseline shifts > 2% blend bracketing bands held volume → pool state OUTPUT MAP Flood map Extent Depth Channel stress Pondage pools Exposure Validation — observed vs modelled Sentinel-1 2025: CSI ≈0.40 · POD ≈0.78 | 2014 extent 13.8 km in 10–20 km | routing NSE 0.72–0.86 routed outflow per node band + pool COGs
The two-layer architecture: live FFD telemetry and barrage controls drive the 1-D network router; its per-node outflows select and blend precomputed 2-D LISFLOOD-FP depth and pondage libraries, served to the map as extent, depth, channel stress, pool and exposure layers — validated against Sentinel-1 (2025), the observed 2014 extent, and FFD hydrographs.

02The live network router — regulation in one dimension

Layer one is a numpy implementation of classical flood routing (routing_core.py) over the full Flood Forecasting Division water network: six rivers, five confluences — Trimmu (Chenab + Jhelum), Panjnad (Chenab + Ravi + Sutlej), Attock (Indus + Kabul) and Guddu (Indus + Panjnad) among them — and the eight link canals that transfer water between rivers. Every reach carries a Muskingum lag-and-attenuation model, so a flood wave takes realistic time to travel and flattens as it goes, and at each confluence the routed outflows simply sum. The router is not a cartoon of the basin; it is the basin's plumbing, traversed whole on every interaction.

Twelve barrages are regulated: Khanki, Qadirabad, Trimmu and Panjnad on the Chenab axis; Rasul on the Jhelum; Balloki and Sidhnai on the Ravi; Islam on the Sutlej; Kalabagh, Chashma, Taunsa and Guddu on the Indus. Each has a single control — a release cap in cusecs. Outflow is min(inflow, cap); the difference is held upstream as a growing pool (Section 11); cap 0 holds everything. This mirrors the operational reality of Pakistan's rivers, where a barrage under threat holds or spills to protect its headworks, and it is why the network view matters: cap Rasul on the Jhelum and the Chenab below Trimmu drops with it, because the confluence sum carries the deficit across rivers.

The router recomputes the entire basin — every reach, confluence, canal and barrage — in well under a second, so the controls are live (input is throttled at 350 ms and the map never waits on the physics). At each node the routed discharge is compared with the FFD's published bankfull capacity, and the excess, max(Q − Q_bankfull, 0), is flagged as spill: the one-dimensional signal that the two-dimensional layers translate into flooded ground. Against the observed 2025 monsoon the router reproduces the FFD gauge hydrographs with Nash–Sutcliffe efficiencies of 0.72–0.86 on the upper reaches (Section 14) — honest skill for a lag-and-attenuation model, and the reason we trust it as the backbone.

O_t = C1·I_t + C2·I_(t−1) + C3·O_(t−1)
Muskingum routing per reach: outflow O at time t from the current and lagged inflow I and the previous outflow, with coefficients fixed by each reach's travel time K and attenuation weight X.

03Discharge data — live telemetry and the historical record

The router's boundary conditions are live. A scraper polls the Flood Forecasting Division's published discharge readings every 15 minutes across 27 stations, including all five rim headworks where water enters Pakistan's regulated system — Tarbela on the Indus, Mangla on the Jhelum, Marala on the Chenab, Jassar on the Ravi and Sulemanki on the Sutlej. These readings are what the app opens on, and they are the anchor every experiment perturbs and returns to (Section 15).

Beneath the live feed sits a long record. We assemble a continuous daily series at Marala Barrage spanning 1997–2025 — 28 years — by reconciling the inflow components (Jammu Tawi, Manawar Tawi and the Chenab mainstem) against the regulated outflows (the Upper Chenab and MR-Link canals plus the downstream release). A dedicated ETL stage unifies the heterogeneous source spreadsheets into a single, gap-checked daily series and derives a flood-event catalogue, tagging each event with its peak, rise rate, duration and flood class. The September 2014 flood — the largest in the record — anchors the historical hindcast: its observed peak of 562,848 cusecs (~15,930 m³/s) agrees to within ~4% of the independently published Marala peak of ~15,300 m³/s, a cross-check that confirms both the record and the unit conversion.

Events are classified with the Federal Flood Commission's Marala thresholds, giving the sliders and the scenario bands a shared vocabulary: Low (100,000 cusecs), Medium (150,000), High (200,000), Very High (400,000) and Exceptionally High (600,000 cusecs). A selected discharge is no longer a single point-source hydrograph: it becomes a headwork boundary condition of the network router (Section 2), which turns it into per-node outflows everywhere downstream, and the two-dimensional layer reads those routed outflows — not the raw gauge — when it decides what to draw (Section 12).

Q_peak(Sep 2014) = 562,848 cusecs ≈ 15,930 m³/s (≈ 15,300 m³/s reported at Marala; within ~4%); 1 m³/s ≈ 35.31 cusecs
2014 benchmark peak and the cusecs–m³/s conversion

04Terrain

The terrain base is the FABDEM V1-2 bare-earth digital elevation model (a Copernicus GLO-30 derivative with buildings and forests removed): a globally consistent 30 m product that is the runnable foundation for the whole pipeline. The production 2-D libraries run on a 120 m resampling of this surface — the resolution at which five whole river valleys (Chenab, Jhelum, Ravi, Sutlej, Indus) can be simulated as 20-day transients within a practical compute budget — while targeted validation hindcasts, such as the 2025 Chenab super-flood (Section 14), run at the native 30 m. Both grids are the same bare-earth surface; only the cell size differs.

We are explicit about what a DEM of this class can and cannot deliver in this setting. The Punjab floodplain is flat alluvium, where the model's vertical error is roughly 1 to 4 m — the same order of magnitude as the flood depths we are simulating. That reality sets our product contract: we report flood extent and discrete hazard bands, not sub-0.5 m depth precision. Quoting depth to the centimetre over terrain this flat would imply an accuracy the elevation data does not support. Consistent with that contract, the pipeline is validated on flood footprint (Section 14) rather than on calibrated absolute depths.

Before any hydraulics run, the DEM is hydro-conditioned so that water routes physically rather than pooling in artefacts. We pit-fill with a priority-flood algorithm, removing the spurious sinks that would otherwise trap flow and break the 2D solver's mass routing. The conditioned surface is then prepared for the engine: the bank-to-bank channel from the JRC surface-water extent (Section 5) is burned into it, and floodplain friction is mapped from ESA WorldCover land cover to per-cell Manning roughness rather than a single basin constant. That triplet — bare-earth FABDEM, JRC channel, WorldCover roughness — is the terrain every valley library is built on.

05River extent from satellite

The hydraulic model needs the full bank-to-bank footprint of each river, not just the water visible on any single day. The Chenab and its neighbours are braided, sand-bed rivers: the active channel migrates across a wide belt of bars and secondary threads that are dry at low flow but carry water at high stage. To capture this true channel, we derive extent from a multi-decadal record of where water has ever been, rather than from a single snapshot.

The primary product is the JRC Global Surface Water layer, which encodes the per-pixel water occurrence frequency from the full Landsat archive (Pekel et al. 2016). We take all pixels with occurrence of at least 10 percent as the channel mask. Because occurrence integrates over more than three decades and every season, the mask includes the bars and secondary threads that are wet only at high flow, capturing the braided channel that defines each river's conveyance corridor.

A single dry-season optical image, by contrast, shows only the narrow active thread and would systematically understate channel width. To track the current active water surface we add a Sentinel-2 L2A product: a multi-scene MNDWI composite that takes the union of water across acquisition dates, flagging water where the Modified Normalized Difference Water Index exceeds zero (Xu 2006). The JRC occurrence mask is the geometry used downstream: it is the channel burned into the DEM for every valley library (Section 9) and the polygon behind the width, depth and channel-stress chain (Sections 6, 7 and 13). The Sentinel-2 composite supplies the present-day active extent for display and cross-checking.

MNDWI = (Green - SWIR) / (Green + SWIR), water where MNDWI > 0
Modified Normalized Difference Water Index (Sentinel-2 active-water mask)

06Channel width

Channel width is measured directly from the bank-to-bank channel polygon rather than assumed from a regional hydraulic-geometry relationship. The channel extent is rasterized and a Euclidean distance transform is computed: every wetted cell is assigned its distance to the nearest channel boundary. This distance field rises to a ridge along the medial axis of the channel, where the value equals the perpendicular distance to the nearest bank — i.e. the local half-width.

To turn that distance field into a width defined consistently across each cross-section, a maximum filter propagates the centre-line (ridge) value outward across the channel. This carries the bank-to-bank measurement from the medial axis to the surrounding cells, so each reach is characterised by its full cross-sectional width rather than the smaller distance recorded at an off-centre pixel. Full width is twice the propagated half-width.

Across the Chenab corridor this yields a median channel width of about 536 m, with maxima near 1,680 m in the widest braided reaches. Because the measurement is taken from the multi-year channel polygon (JRC Global Surface Water 35-year Landsat occurrence), it captures the full braided cross-section — including bars that are wet only at high flow — rather than just the active thread. As an independent check, the resulting widths are consistent in magnitude with the SWORD river database for these reaches. The per-reach width field feeds the depth inversion (Section 7) and the channel-stress screening layer (Section 13).

W = 2 · max(EDT) per cross-section
Channel width from the Euclidean distance transform (EDT) of the channel polygon

07Channel depth — hydraulic inversion

Satellites see only the water surface, and a DEM over a wet channel measures that surface, not the bed beneath it. For the Chenab no measured bathymetry exists, and the global "depth" products that might fill the gap are themselves empirical formulas rather than observations. Our answer is to recover channel depth physically rather than assume it — and the key asset that makes this tractable is that, unlike most remote-sensing studies, we have a measured discharge record at Marala.

Knowing Q collapses the hard inverse problem into a direct calculation. Treating the wide braided Chenab as a channel where hydraulic radius approximates mean flow depth (R ≈ d), Manning's law combined with continuity (Q = A·V, A = W·d) inverts to a closed form for depth. We evaluate it with the measured channel width W from the extent step, the measured bed slope S = 0.33 m/km (0.00033, surveyed Marala→Khanki), and a braided sand-bed roughness n ≈ 0.025–0.035 — n being the dominant remaining uncertainty, which we bound rather than fit.

Evaluated at the bankfull "within-channel" threshold (Q ≈ 100,000 cusecs ≈ 2,830 m³/s), the inversion yields a mean bankfull depth of about 3.6 m. This is independently corroborated by the Andreadis et al. (2013) global bankfull relation d = 0.27·Q^0.3, which gives roughly 3 m for the same discharge. Applied per reach, continuity does the physically correct thing automatically: narrow necks run deeper and wide braided reaches run shallower, producing a spatially varying depth field rather than a single constant. By construction this is a reach-averaged mean depth; the thalweg runs deeper (on the order of 1.5–2×), and because no in-situ bathymetry exists the field is a physics-based estimate whose roughness parameter is bounded from the literature rather than fitted.

The approach sits in a well-established lineage: at-station hydraulic geometry (Leopold & Maddock 1953), mass-conserved flow-law inversion and AMHG (Gleason & Smith 2014), the BAM family (Hagemann & Durand 2017) and geoBAMr (Brinkerhoff et al. 2020) — run here in the well-posed direction, discharge known, depth solved. Two things in the current platform stand on this field. It is the depth term inside the channel-stress screening layer (Section 13), and it is the burn used by the original 30 m Chenab corridor library whose 2014 hindcast still ships as the platform's historical tab: that conditioned run lands the modelled 2014 corridor inside the observed 10–20 km envelope (Section 14). The production 120 m valley libraries burn the channel from the JRC extent directly (Section 9); the inversion remains the Chenab-specific, physics-first version of that step.

d = [ n · Q / ( W · √S ) ] ^ (3/5)
Manning–continuity depth inversion (R ≈ d for a wide channel): Q = discharge (m³/s), W = measured width (m), S = bed slope, n = Manning's roughness
Manning continuity depth inversion: river channel cross-section A river channel cross-section with floodplain on both sides, a wetted water surface of width W, mean depth d below it, bed slope S, and the equation d equals n Q over W root S, raised to the three fifths power. Depth is inferred from measured discharge, width and slope because a satellite cannot see the river bed. River Channel Cross-Section — Manning Continuity Depth Inversion Satellite measures from above: width W, discharge Q, slope S — never the bed floodplain floodplain W wetted surface width d mean depth (inferred) bed slope S d = [ n·Q W · S ] 3/5 n Manning roughness Q measured discharge W surface width S bed slope Solve for: d mean depth bed is not observed Depth d is recovered from measured discharge Q, surface width W, and slope S — a satellite cannot see the bed directly.
Channel cross-section behind the Manning–continuity inversion: satellite-measurable width W, discharge Q and slope S recover the unobserved mean depth d — about 3.6 m at bankfull, deeper at narrow necks, shallower across wide braids.

08Channel conditioning — the 2014 Chenab hindcast

A satellite-era DEM measures the water surface, not the channel bed, so the raw terrain has no incised channel to carry low and moderate flows: an inflow spreads laterally as soon as it enters the domain, producing an unphysically wide, shallow sheet. Channel conditioning fixes this by carving a channel into the terrain before it reaches the hydraulic engine. The original Chenab corridor build — the library that still powers the platform's September 2014 hindcast tab — is the fullest expression of the idea.

There, we lower every DEM cell that falls inside the river polygon (Section 5) by the local Manning-continuity depth field (Section 7), producing a channel-conditioned DEM in which the bed sits a physically estimated distance below the surrounding banks. Because the depth field is reach-resolved, continuity makes narrow necks deeper and wide braided reaches shallower along the corridor. LISFLOOD-FP then keeps flow confined to this incised channel and spills onto the floodplain only where discharge genuinely overtops the banks.

Conditioning is applied as a discrete, versioned step so a no-burn baseline and the conditioned terrain can be compared directly. On the September 2014 event it is the main correction to floodplain over-spread: the modelled perpendicular extent falls from roughly 16.8 km in the no-burn baseline to ~14 km with the physics-based burn, inside the 10–20 km reported for the flood (Khan et al.), and every contingency score against the observed mask improves (Section 14). The production whole-valley libraries inherit the lesson in simplified form — the channel is incised from the JRC occurrence extent across all five valleys (Section 9) — while the reach-resolved inversion burn remains the Chenab 2014 story and the depth term of the screening layer.

z_conditioned(x, y) = z_DEM(x, y) − d(x, y), for cells inside the river polygon
Channel burn: subtract the reach-resolved Manning-continuity depth field d (Section 7) from the DEM within the river polygon.

092D hydraulics — LISFLOOD-FP

The hydraulic engine is Bristol LISFLOOD-FP 8.0.3, built from the published Zenodo source and run on its local-inertial (acceleration) solver. Rather than solving the full Saint-Venant equations, this solver advances a simplified form of the shallow-water equations across the raster floodplain, retaining the local-inertial term while dropping only the convective acceleration. This gives stable, mass-conserving 2D flow at a fraction of the cost of a full dynamic-wave or CFD solver — which is what makes it tractable to simulate five whole river valleys (Chenab, Jhelum, Ravi, Sutlej and the braided Indus) end to end at 120 m, rather than a single corridor.

The channel is carried in the terrain itself: the JRC bank-to-bank extent (Section 5) is burned into the FABDEM surface so flow stays in-channel until it genuinely overtops, and floodplain friction comes from ESA WorldCover land cover mapped to per-cell Manning roughness. The model is driven by text inputs — the ASCII DEM, a .par configuration file and paired .bci/.bdy boundary files — and the forcing is a calibrated Gaussian flood wave, peaking at day 5.0 with a standard deviation of 2.2 days, scaled to each discharge band and run as a 20-day transient: long enough for the wave to traverse the valley and the floodplain to fill and begin to drain. One domain decision matters downstream: the Jhelum and Ravi domains continue past their confluences into the Chenab channel, so water arriving at a junction flows on downstream instead of pooling artificially at a domain edge.

The solver integrates these shallow-water dynamics forward in time and writes water depth on the same grid; what the platform stores is the maximum-depth envelope of each 20-day run — the deepest water every cell sees over the whole event — exported as Cloud-Optimized GeoTIFFs. It is worth being explicit that the physics is two-dimensional, depth-averaged flow over terrain — not a 3D fluid simulation. The 3D appearance in the viewer (deck.gl buildings, extruded river and channel layers) is visualisation of those 2D depth grids, not an additional dimension of the hydraulic computation.

∂q/∂t + g·h·∂(h+z)/∂x + g·n²·q·|q| / h^(7/3) = 0
Local-inertial (acceleration) form of the shallow-water momentum equation solved per cell; q = unit-width discharge, h = water depth, z = bed elevation, n = Manning's roughness.

10The scenario libraries — transient physics, precomputed

The hydraulic engine is not run live. For each of the five valleys, LISFLOOD-FP is executed ahead of time at seven discharge bands spanning 50,000 to 1,000,000 cusecs — from below the FFC "Low" threshold to well past the observed record. Every band is a full 20-day transient driven by the calibrated Gaussian wave (Section 9), and what is stored is its maximum-depth envelope as a Cloud-Optimized GeoTIFF, so the library is an indexed stack of unsteady-physics inundation surfaces keyed by river and discharge band — not steady-state stills.

Seven bands per valley is deliberately few, because the serving layer interpolates: at request time the platform blends the two bracketing bands around each target discharge, so the effective resolution in Q is far finer than seven steps, and the band a given reach receives is set by the live router rather than by a single global slider (Section 12 describes the coupling). The upper bands (800,000 and 1,000,000 cusecs) sit above the record peak of 562,848 cusecs (September 2014); they are physically modelled but presented as hypothetical extremes for stress planning rather than validated events.

The September 2014 event is carried separately from the band matrix: the original 120 m Chenab-corridor library, driven by the full 2014 Marala hydrograph over the Manning-inversion-conditioned terrain (Section 8), powers the platform's dedicated 2014 hindcast tab. It is the platform's history — the corridor-scale build the whole-basin architecture grew out of — and it remains the cleanest like-for-like comparison against the best-documented flood in the record.

11Pondage — what holding water costs

Holding water at a barrage is not free: it ponds. When an operator caps a release, the held volume accumulates upstream of the gates, and a model that let water vanish there would make "hold everything" look like a magic wand. So the platform models the pool with the same physics as the flood. For each regulated barrage and inflow band we run a paired LISFLOOD-FP experiment — a free reference run, and a walled run identical in every respect except a two-row wall across the channel at the barrage — and snapshot the growing pool every 0.25 days.

The pool is defined as the walled-minus-free depth difference, which cleanly removes water that is merely in transit down the channel and leaves only what the closure impounds. Integrating that difference gives the pool volume at each snapshot, and the snapshots are indexed into eight log-spaced volume bands from 10⁷ to 10⁹ m³, stored as multi-band COGs. At serving time the router integrates each barrage's held hydrograph — cusecs held, over time, into cubic metres — and the matching real solver pool state renders behind the gates: a shallow-water solution shaped by the actual terrain, not a cartoon ellipse.

The numbers behind the claim: 36 pondage libraries (12 barrages × 3 inflow bands), every one passing a mass-closure gate that requires the integrated pool volume to agree with the held inflow to within 2% — the median closure error is ~0.01%. Deep, near-static pools exposed a known weakness of local-inertial schemes — oscillation over still water — which we damp with the solver's θ parameter set to 0.7 for the pondage runs; that trades a little numerical diffusion for stable pools, and we say so rather than hide it. Pools beyond ~2.9×10⁸ m³ reach the edge of some simulation domains and clip visually; the volume accounting remains correct (Section 16).

V_pool(t) = Σ_cells [ h_walled(t) − h_free(t) ] · A_cell
Pool volume as the walled-minus-free depth difference summed over the domain — transit water cancels, impounded water remains. Snapshots every 0.25 d, indexed into 8 log-spaced volume bands (10⁷–10⁹ m³).

12Coupling the layers — ratio-banded serving

The join between the 1-D router and the 2-D libraries is what makes an intervention visible on the ground. When the routed state changes, the inundation API compares each river's per-node outflow against its baseline (the uncontrolled routed state) and splits the river into latitude slab groups wherever the outflow-to-baseline ratio changes by more than 2%. Each slab is then served at its headwork-equivalent discharge — the headwork Q scaled by the slab's ratio — blended between the two bracketing library bands. The result is a single seamless raster in which different stretches of the same river legitimately show different floods.

The slabs capture both directions of causality. Cap Khanki and the Chenab below Khanki drops band by band while the reach above holds its level — an own-river effect. Cap Rasul on the Jhelum and the Chenab below Trimmu re-bands too, because the router's confluence sum carries the deficit across rivers — a cross-river effect no per-river lookup could produce. The same request resolves each barrage's held volume into its pondage pool state (Section 11), so the map shows, simultaneously, less water downstream and more water behind the gates.

Two approximations are worth naming. Between band levels the depth field is a blend, not a re-solve: the physics is exact at the seven library bands and interpolated between them. And where two rivers' rasters meet near a confluence, the served extent is their union rather than a coupled joint solution; the Jhelum and Ravi domains continuing into the Chenab channel (Section 9) keeps water moving through the junctions, but true joint-backwater Y-domains are the next build (Section 16).

Q_slab = Q_headwork × ( Q_routed(slab) / Q_baseline(slab) ), slab split where the ratio changes > 2%
Headwork-equivalent discharge per latitude slab: the library band is chosen by scaling the headwork discharge by the slab's routed-to-baseline ratio, then blending the bracketing bands.

13Channel conveyance & overtopping

Alongside the 2D inundation grids, the platform reports a complementary screening diagnostic: how close the channel is to its conveyance limit at any given discharge. For each reach we evaluate the river's carrying capacity Q_cap from Manning's equation, using the bank-to-bank channel polygon (Section 5), the measured cross-section width (Section 6), the physics-based channel depth (Section 7), and the measured bed slope S = 0.00033 (0.33 m/km, Marala-Khanki reach, Awan 2003).

The computation runs as a raster over the river polygon, masked to the channel. Flow area is A = W*d and, for these wide alluvial cross-sections, the hydraulic radius reduces to the mean depth, R ~= d. With Manning's n bounded to the braided sand-bed range (0.025-0.035, the dominant uncertainty), this yields a per-reach Q_cap field: continuity makes narrow necks lower-capacity and wide reaches higher-capacity, so the diagnostic resolves where conveyance is relatively constrained rather than assuming a single section governs. Because the depth term is a reach-mean estimate, Q_cap is a reach-mean conveyance (the thalweg is locally deeper), so the result is best read as a relative ranking of overtopping susceptibility along the corridor.

The map colours the ratio of the current discharge to local capacity, discharge / Q_cap, and recolours it live as the routed state changes. Reaches well below capacity render green (flow within banks), reaches approaching unity render amber (near overtopping), and reaches at or above unity render red (channel overtops). This screening raster is distinct from the router's spill flag (Section 2), which compares routed discharge with the FFD's published bankfull capacity at network nodes: Q_cap resolves relative susceptibility continuously along the channel between those nodes. Like the depth field it builds on, it is a physics-based diagnostic whose roughness is bounded rather than fitted, so the colours should be read as relative susceptibility, not absolute thresholds.

Q_cap = (1/n) * A * R^(2/3) * sqrt(S), with A = W*d and R ~= d
Manning per-reach channel carrying capacity
River reach coloured by channel stress (discharge / carrying capacity) A plan-view river reach. Wide segments are green (within banks), narrowing to amber near capacity, and red at narrow necks where flow overtops. The Manning capacity equation is labelled. Channel stress along a river reach Stress = discharge Q / carrying-capacity Q_cap. Flow direction: left to right. flow Q narrow neck — overtopping narrow neck — overtopping wide — within banks wide — within banks Channel stress S = Q / Q_cap Green S < 0.8 — within banks Amber 0.8 ≤ S < 1.0 — near capacity Red S ≥ 1.0 — overtopping banks Stress rises where the channel narrows (Q_cap falls). Carrying capacity — Manning Qcap = (1/n) · A · R2/3 · √S n Manning roughness A flow area R hydraulic radius A/P S bed slope
Per-reach Q_cap raster on the river polygon, coloured live by discharge / Q_cap: green within banks, amber near capacity, red overtopping.

14Validation

The test we hold the platform to is cell-by-cell agreement with radar. For the 2025 Chenab super-flood — the largest event since 2014 — we ran a 30 m hindcast and scored it against an observed water mask built by Sentinel-1 change detection: a VV-backscatter decrease criterion for open water plus a VH criterion for flooded vegetation, computed on radiometrically terrain-corrected (RTC) scenes from the Microsoft Planetary Computer. Radar sees through monsoon cloud, so unlike optical benchmarks this reference actually sees the flood. The scores: CSI ≈ 0.40 over the full footprint, rising to 0.44 where modelled depth exceeds 1 m, with POD ≈ 0.78 — the model finds roughly four-fifths of the observed flooding. Single-scene placement checks score POD 0.86–0.92. The production 120 m whole-basin library, scored against the same reference, gives POD 0.82 / CSI 0.34 — coarser, as expected, but the same flood in the same places.

Agreement is quantified with the standard contingency-table skill scores over cells classified as hits (A, wet in both), false alarms (B, wet only in the model) and misses (C, wet only in the observation). CSI is the headline metric because it ignores correctly-dry cells, which otherwise dominate and inflate agreement; POD, FAR and bias resolve whether residual disagreement is under- or over-prediction. The September 2014 flood remains the historical benchmark: driven by the full 2014 Marala hydrograph, the channel-conditioned hindcast produces a mean cross-valley extent of 13.8 km, inside the 10–20 km corridor reported from field and remote-sensing analysis (Khan et al.). The earlier optical comparison for 2014 — against the Global Flood Database MODIS mask (DFO_4179) — survives as a footnote with its limitation stated: that monsoon flood was heavily cloud-covered, the MODIS mask captured on the order of 400 km² of a footprint well over 1,000 km², and the resulting CSI ≈ 0.11 says more about the benchmark than the model. That ceiling is exactly why the SAR reference now leads.

The 1-D layer is validated separately: routed hydrographs against the FFD's observed 2025 gauge records give Nash–Sutcliffe efficiencies of 0.72–0.86 on the upper reaches. As an order-of-magnitude overlay, the JRC CEMS-GloFAS RP50 return-period layer is served across the whole Punjab — a sanity benchmark produced with the same LISFLOOD-FP family, never a calibration target. And what is not yet validated is stated plainly in Section 16: the low discharge bands (50,000–150,000 cusecs) have no SAR reference at all, because radar change detection only yields a clean mask when there is a large flood to see.

CSI = A / (A + B + C); POD = A / (A + C); FAR = B / (A + B); bias = (A + B) / (A + C)
Extent-agreement contingency metrics (A = hits, B = false alarms, C = misses)
Flood-extent validation: 2x2 contingency / overlap diagram Modelled versus observed flooded area showing Hits A (overlap), False alarms B (model-only), Misses C (observed-only), with the Critical Success Index CSI equal to A divided by the sum of A, B and C. Flood-extent validation: modelled vs observed 2×2 overlap of flooded area B False alarms model only A Hits overlap C Misses observed only Modelled flooded Observed flooded CRITICAL SUCCESS INDEX CSI = A / (A + B + C) A Hits B False alarms C Misses
Contingency scoring of modelled vs observed extent: hits, false alarms, misses. Applied to the 2025 Chenab super-flood against the Sentinel-1 change-detection mask (CSI ≈ 0.40, POD ≈ 0.78 at 30 m; POD 0.82 / CSI 0.34 for the 120 m whole-basin library) and to the 2014 hindcast (13.8 km cross-valley extent vs the observed 10–20 km).

15Live now-cast

The platform opens on now. The 15-minute FFD scrape (Section 3) keeps 27 stations current — the five rim headworks plus the barrage and gauge chain downstream — and the router immediately routes the latest readings through the network, so the first paint is Pakistan's present hydraulic state, not a demo scenario. During a developing flood this turns the map into a what-if console on the actual event: the operator watches the observed wave move down the routed network and can test holds at the barrages ahead of it.

Every experiment starts from and returns to reality. "Set to LIVE conditions" snaps all headwork discharges and barrage caps back to the latest feed, so a hypothetical never silently masquerades as the current state. The distinction is enforced in the interface rather than left to the user's memory.

We are precise about what "live" buys. The now-cast is a routing forecast, not a weather forecast: lead time comes from hydraulic travel time — a wave observed at a rim station takes one to seven days to reach the lower basin, and the router propagates it — but nothing upstream of the telemetry is predicted. Rainfall-driven forecasting above the gauges is not yet built (Section 16).

16Honest limits

Validation coverage is uneven, and we say where. The SAR-scored skill (Section 14) belongs to the big-flood bands: the low bands (50,000–150,000 cusecs) have no radar reference, because change detection only produces a clean observed mask when there is a large flood to see, so at low flows the libraries are unscored physics on validated terrain. On the Indus, the low band over-spreads west of the braid belt and should be read as screening-grade there. And the depth contract of Section 4 stands everywhere: extent and hazard bands, not centimetre depths, with Manning's n bounded from the literature rather than fitted.

The architecture has known seams. Where two rivers meet, the served rasters are a per-river union rather than a coupled joint solution — the Jhelum and Ravi domains flowing on into the Chenab channel keeps water moving through the junctions, but genuine joint-backwater Y-domains at Trimmu and Panjnad are the next build. Between library bands, depth is blended rather than re-solved. Pondage pools beyond ~2.9×10⁸ m³ clip at some domain edges — the rendered pool stops at the boundary while the volume accounting stays correct.

Finally, the forecast horizon. Lead time today is hydraulic travel time downstream of the telemetry — one to seven days depending on where the wave is and where you care about — because rainfall-driven forecasting upstream of the gauges is not yet built. The platform tells you what the rivers will do with the water they already carry; it does not yet tell you what the sky will add. We publish these limits because the platform's value depends on being believed exactly as far as it deserves.

17Visualisation

The front end is a Next.js application built on MapLibre GL with a dark Carto basemap, chosen so that flood-depth and channel-stress colour ramps read clearly against an unlit background. deck.gl renders the interactive 3D layers on top of the map in a single shared WebGL context: extruded buildings for exposure, the river network, and the channel geometry. Map and deck.gl share one camera, so vector overlays, raster tiles, and 3D meshes stay registered as the user pans, zooms, and tilts.

Continuous raster fields — flood depth, pondage pools and channel stress — are served as Cloud-Optimised GeoTIFFs (COGs) through TiTiler, a dynamic tile server. Rather than pre-baking a colour image per scenario, TiTiler reads only the requested tiles from the COG and recolours them on the fly, applying a rescale window and a colour map (turbo for depth, a green-to-red ramp for stress) so a single source raster can be restyled instantly. Vector geometry — river network, both-bank channel extent, and the JRC hazard overlay — is delivered as PMTiles or GeoJSON, with the FastAPI plus PostGIS/SFCGAL backend supplying the underlying features.

Every data product is an independent, toggleable layer, and the controls are the router's. Per-barrage release-cap sliders and headwork discharges drive the 1-D re-route (throttled at 350 ms); the served slabs re-band and blend (Section 12); pondage pools render behind capped gates at their integrated volume (Section 11); and the channel-stress raster recolours — all without re-running the hydraulic engine. A depth-highlight slider drives the TiTiler rescale window, a dedicated tab replays the September 2014 hindcast from the original Chenab corridor library (Section 8), and "Set to LIVE conditions" (Section 15) returns everything to the current feed.

Data sources

Live telemetryFFD discharge feed, 27 stations incl. 5 rim headworks (15-min scrape)boundary conditions
Discharge recordMarala daily 1997–2025 (WAPDA / Irrigation)history & flood classes
DEMFABDEM V1-2 bare-earth (GLO-30 derivative)terrain — 120 m libraries, 30 m hindcasts
Water historyJRC Global Surface Water (35-yr Landsat)bank-to-bank channel burn
Land coverESA WorldCover 10 mManning roughness map
ImagerySentinel-2 L2A (Earth Search / AWS)active-water MNDWI
SAR referenceSentinel-1 RTC (Microsoft Planetary Computer)2025 observed flood mask
Slope / bankfullWMO/GWP Chenab case study (Awan 2003)S = 0.33 m/km
Hazard mapJRC CEMS-GloFAS return-period layersbenchmark overlay
Hydraulic engineLISFLOOD-FP 8.0.3 (Bristol)2D shallow-water

References

  1. Shaw et al. (2021) LISFLOOD-FP 8.0, Geosci. Model Dev. 14. https://doi.org/10.5194/gmd-14-3577-2021
  2. Pekel, Cottam, Gorelick & Belward (2016) High-resolution mapping of global surface water and its long-term changes, Nature 540.
  3. Awan (2003) Flood Management: River Chenab Marala-Khanki, WMO/GWP APFM case study.
  4. Khan et al., Extent of 2014 Flood Damages in Chenab Basin, IntechOpen ch. 62591.
  5. Dottori et al. / JRC CEMS-GloFAS Global Flood Hazard maps.
  6. McCarthy (1938) The unit hydrograph and flood routing, US Army Corps of Engineers (the Muskingum method); Cunge (1969) On the subject of a flood propagation computation method, J. Hydraulic Research 7.
  7. Flood Forecasting Division, Pakistan Meteorological Department — daily river & barrage discharge bulletins. https://ffd.pmd.gov.pk
  8. Khan et al., Extent of 2014 Flood Damages in Chenab Basin, IntechOpen ch.62591 (https://www.intechopen.com/chapters/62591)
  9. Awan (2003), Flood Management: River Chenab Marala–Khanki, WMO/GWP APFM case study (https://www.floodmanagement.info/publications/casestudies/cs_pakistan_chenab_full.pdf)
  10. Discharge record: WAPDA / Provincial Irrigation, Marala daily 1997–2025; flood classes per FFC Marala scheme
  11. Copernicus GLO-30 DEM, dataspace.copernicus.eu (open licence)
  12. FABDEM V1-2 (Hawker et al. 2022, Univ. Bristol; bare-earth 30 m, CC-BY-NC)
  13. Zanaga et al. (2021) ESA WorldCover 10 m land cover. https://worldcover2021.esa.int
  14. Pekel, Cottam, Gorelick & Belward (2016), High-resolution mapping of global surface water and its long-term changes, Nature 540.
  15. Xu (2006), Modification of normalised difference water index (MNDWI) to enhance open water features, Int. J. Remote Sensing 27(14).
  16. JRC Global Surface Water v1.4 (CC BY), global-surface-water.appspot.com; Sentinel-2 L2A via earth-search.aws.element84.com (open).
  17. Pekel, Cottam, Gorelick, Belward (2016) High-resolution mapping of global surface water and its long-term changes, Nature 540. https://global-surface-water.appspot.com
  18. Altenau et al. (2021) The Surface Water and Ocean Topography (SWOT) Mission River Database (SWORD), Water Resources Research. http://gaia.geosci.unc.edu/SWORD/
  19. Leopold & Maddock (1953), The hydraulic geometry of stream channels, USGS PP 252
  20. Andreadis et al. (2013), A simple global river bankfull width and depth database, Water Resources Research, doi:10.1002/wrcr.20440 (data: zenodo 61758)
  21. Gleason & Smith (2014), Toward global mapping of river discharge using inverted hydraulic geometry (AMHG), PNAS 111(13), doi:10.1073/pnas.1317606111
  22. Hagemann & Durand (2017), BAM (Bayesian AMHG–Manning), WRR; Brinkerhoff et al. (2020), geoBAMr, WRR (code: github.com/markwh/bamr, github.com/craigbrinkerhoff/geoBAMr)
  23. Awan (2003), Flood Management: River Chenab Marala–Khanki, WMO/GWP APFM case study (slope S = 0.33 m/km, bankfull threshold)
  24. Neal, Schumann & Bates (2012) A sub-grid channel model for simulating river hydraulics and floodplain inundation, Water Resources Research.
  25. Neal et al. (2021) Estimating river channel bathymetry in large-scale flood models, Water Resources Research. https://doi.org/10.1029/2020WR028301
  26. Shaw et al. (2021) LISFLOOD-FP 8.0, Geoscientific Model Development 14. https://doi.org/10.5194/gmd-14-3577-2021
  27. Khan et al. Extent of 2014 Flood Damages in Chenab Basin, IntechOpen ch. 62591. https://www.intechopen.com/chapters/62591
  28. Bates & De Roo (2000) A simple raster-based model for flood inundation, J. Hydrology 236.
  29. Bates, Horritt & Fewtrell (2010) A simple inertial formulation of the shallow water equations (acceleration solver), J. Hydrology 387.
  30. Shaw et al. (2021) LISFLOOD-FP 8.0, Geosci. Model Dev. 14, 3577. https://doi.org/10.5194/gmd-14-3577-2021 (code: Zenodo 4073011).
  31. Bates, Horritt & Fewtrell (2010) A simple inertial formulation of the shallow water equations (local-inertial/acceleration solver), J. Hydrology 387.
  32. Neal, Schumann & Bates (2012) A sub-grid channel model for simulating river hydraulics and floodplain inundation, Water Resources Research 48.
  33. de Almeida, Bates, Freer & Souvignet (2012) Improving the stability of a simple formulation of the shallow water equations for 2-D flood modeling, Water Resources Research 48 (the θ diffusion weighting).
  34. Manning's equation (open-channel uniform flow)
  35. Awan (2003), Flood Management: River Chenab Marala-Khanki, WMO/GWP APFM case study - S = 0.33 m/km
  36. Shaw et al. (2021), LISFLOOD-FP 8.0, GMD 14, https://doi.org/10.5194/gmd-14-3577-2021
  37. Khan et al., Extent of 2014 Flood Damages in Chenab Basin, IntechOpen ch. 62591, https://www.intechopen.com/chapters/62591
  38. Tellman et al. (2021), Satellite imaging reveals increased proportion of population exposed to floods, Nature 596 (Global Flood Database), https://doi.org/10.1038/s41586-021-03695-w
  39. Bates & De Roo (2000), A simple raster-based model for flood inundation simulation, J. Hydrology 236 (source of the F/CSI fit statistic)
  40. Dottori et al. / JRC CEMS-GloFAS Global Flood Hazard maps (return-period RP50/RP100 tiles)
  41. Sentinel-1 RTC (radiometrically terrain-corrected) scenes via the Microsoft Planetary Computer, planetarycomputer.microsoft.com
  42. MapLibre GL JS — https://maplibre.org
  43. deck.gl — https://deck.gl
  44. TiTiler dynamic tile server (Development Seed) — https://developmentseed.org/titiler
  45. Cloud-Optimised GeoTIFF (COG) — https://cogeo.org
  46. PMTiles (Protomaps) — https://protomaps.com/docs/pmtiles

Limitations: the DEM caps depth precision in flat terrain (hazard bands, not metres); low discharge bands (50k–150k cusecs) are not SAR-validated; the Indus low band is screening-grade west of the braid; confluence rasters are a per-river union; pondage pools beyond ~2.9×10⁸ m³ clip at domain edges (volume stays correct); lead time is hydraulic travel time downstream of telemetry (1–7 days) — rainfall-driven forecasting upstream of the gauges is not yet built. Section 16 gives the full list.