# Petrel Flow Accumulation Layer v0.3 — Methodology ## Data pipeline ``` GEDTM30 v1.2 bare-earth DEM (© OpenGeoHub / Ho et al. 2025, CC BY 4.0) │ ▼ read each 10° land tile with a 1° (~110 km) halo of neighbour DEM land mask: GSHHG-f · closed basins: depression-depth scan │ ▼ pysheds: fill_pits → fill_depressions → resolve_flats → flowdir │ ▼ per-tile D8 accumulation (truncated at every tile edge — an intermediate, never a product) │ ▼ CROSS-TILE PROPAGATION — iterate edge exports between tiles until the transfers stop changing; drainage in AREA units across every seam │ reservoir nodes (Volga/Kama): inflow terminates at the polygon, outflow re-injected at the dam, exit carved along the mapped river │ ▼ propagated contributing area ──────────────────► accumulation (km²) │ ┌───────────────┴───────────────┐ ▼ ▼ channels = contributing area > 1 km² TWI = ln( area·1e6 / tan slope ) │ │ ▼ ▼ Euclidean distance transform → metres ── twi │ ▼ ── distance_to_drainage │ ▼ mosaic to global COG → coast_clip (gshhg_f, all_touched) → uint16 reencode with PETREL_SCALE (cap distance 50 km, TWI 35) → overviews (distance: NEAREST — see below) ``` ## Flow derivation — pysheds D8 Per DEM window, the drainage grid is computed with the standard pysheds depression-handling chain: **`fill_pits`** (single-cell sinks) → **`fill_depressions`** (flood multi-cell depressions to their spill level) → **`resolve_flats`** (impose a gradient across flat areas) → **`flowdir`** (D8: each cell drains to its steepest of 8 neighbours) → **`accumulation`**. Two masks are applied before routing: - **Land.** The GSHHG-f shoreline masks the ocean out of the flow grid. Without it the sea is routed as flat terrain and coastal cells inherit an invented upstream area. - **Closed basins.** Terminal sinks (Caspian, Aral, Chad, Eyre, Tarim, Saharan sebkhas) are found by depression-fill depth — a filled depth over 10 m on a basin of at least 1,000 km² — not from a hard-coded list. Water reaching one ends there, as it does in nature. Those cells carry no drainage and are **NoData** in both shipped rasters. ## Cross-tile propagation Flow accumulation is inherently global: water crosses tile boundaries, and a tile computed in isolation under-counts everything entering from outside. Measured at reference gauges, the per-tile field captures **2.5–25%** of true drainage. It is an intermediate and never ships. v0.3 reconciles the tiles by iteration. Each round, every tile is accumulated with its pending inflow applied to its edges, and what leaves each edge is pushed to the neighbour; rounds repeat until the transfers stop changing. Tiles are visited in an approximate upstream-first order so water crosses a continent in few rounds rather than one tile per round. **Drainage is carried in area, not in cells.** A grid cell is a different area at every latitude — 23% larger at 45°N than at 55°N — so exporting a cell count across a north–south seam and re-counting it at the neighbour's cell size inflates or deflates the river. Measured before the fix: the Mississippi read 108.6% of its published basin at Vicksburg (it flows south, into larger cells) and the Nile 90.0% at Aswan (north, into smaller cells). Each row's weight is its own cell area relative to the tile centre, and every seam transfer is rescaled by the ratio of the two tiles' cell areas. **Known residual.** A drainage divide running along a tile boundary can make two tiles feed each other, and a handful of island tiles never fully settle; they creep by less than 0.1% of their inflow per round. The build accepts a residual of a few such tiles once no reservoir node is still changing. The structural fix — solving the boundary-cell graph as a DAG in one pass rather than iterating tile transfers (Barnes et al.) — is the flagged upgrade. ## Reservoirs as flow nodes Flat water has no gradient, so inside a large reservoir the DEM cannot say which way the river goes, and D8 picks arbitrarily. On the Volga this severed the river: the Volga/Kama confluence lies inside the Kuybyshev Reservoir on a tile seam, and the mainstem was routed east across it, leaving Volgograd with **15.5%** of its published basin. For the eight Volga/Kama reservoirs (Rybinsk, Cheboksary, Kama, Votkinsk, Nizhnekamsk, Kuybyshev, Saratov, Volgograd) v0.3 replaces DEM routing with an explicit node: 1. Everything draining into the polygon is summed — the reservoir's inflow. 2. That inflow plus the reservoir's own surface is re-injected at one cell: the HydroLAKES pour point, moved onto the first real river cell below the dam. 3. The exit is traced along the **Natural Earth** mapped centerline, across tile boundaries, until it reaches the next reservoir (a handoff, credited to that reservoir's inflow), the sea, or open river. Where no centerline exists, the corridor follows the DEM's own downhill path; failing that, the shortest land path to the nearest strictly lower lake. 4. The reservoir's own cells are then filled with the contributing area passing through it, so the river reads continuously across the water body rather than vanishing into it. Guards, each added after a measured failure: no node may feed itself, directly or after a path that leaves and re-enters; the node graph is checked for cycles and each is broken at the edge into the smallest watershed; a lake whose rim already carries far more than its own published watershed is a river passing by, not a lake to terminate, and keeps DEM routing. **All other lakes keep DEM routing.** A global version of this treatment was built and measured: it fixed the Volga and broke four other basins at floodplain lakes. See [`validation_report.md`](validation_report.md). ## Contributing area, distance to drainage, TWI **Contributing area** is the propagated upstream area of each cell in **km²**, including the cell itself. **Distance to drainage** is the Euclidean distance transform (`scipy.ndimage.distance_transform_edt`) of the complement of the channel mask `contributing area > 1 km²`, converted from pixels to **metres** with a latitude-scaled metres-per-pixel factor (`res° × 111 320 × √cos(lat)`). Cells on a channel read 0. The far tail is capped at **50 km**, the top of the uint16 encoding. Measured over the whole raster, the cap is reached on 1,683 cells (118 km²): islands and cays too small to contain a cell with 1 km² of contributing area, so no channel exists on them at all. It is not reached in deserts, which carry organised drainage within 50 km. The threshold is an **area**, not a cell count: 1000 cells is 77 km² at the equator and 39 km² at 60°N, so a cell-count threshold makes the definition of "a river" drift with latitude. **TWI** is `ln( contributing_area_m² / tan(slope) )`, with slope the local gradient magnitude of the DEM floored at 1e-4 rad to avoid dividing by zero on flats. It is unitless, typically spans ~0–30, and is capped at **35** for encoding. High TWI marks flat, water-gathering, saturation-prone terrain; low TWI marks steep, well-drained slopes. Because TWI is a logarithm of area, it compresses hard: a 23% error in drainage moves it by 0.2 units, and on large rivers it saturates at the cap. ## Overview resampling Overviews are built bilinear for continuous bands, with one exception. `distance_to_drainage` builds **nearest** overviews, because its 50 km cap encodes to 65534 and NoData is 65535: blending the two across the coastline produced values that landed on the cap, so zoomed-out views (previews, and the API's tile server, both of which read overviews) showed spurious "50 km" pixels over open ocean. Nearest removes the blend entirely. TWI keeps bilinear; its capped cells are inland and it shows no such artefact at any level. ## Coastline masking The canonical `coast_clip` step (GSHHG-f, `all_touched=True`) removes ocean cells after the global raster is assembled. Any cell touching real land keeps its value; open ocean and — by Petrel standard — Antarctica (GSHHG L1 only) return NoData. The product ships at its native derivation grid with no coarse-to-fine downscale, so there is no Carolina-Beach-class coastal-hole regression here; the audit still verifies coastal land cells resolve to data and ocean / Antarctica resolve to NoData. ## Sanity checks the COG must pass 1. **PETREL_COAST_CLIP + PETREL_SCALE tags present** — the canonical mask and encoding steps ran. 2. **Basin capture at gauges** — measured against published upstream areas at nine reference gauges before the bundle is built; see [`validation_report.md`](validation_report.md). An over-count is treated as proof of a specific defect, not as noise. 3. **Structure** — distance-to-drainage reads 0 along major river lines, including through reservoirs, and high on ridgelines; TWI reads high on valley floors and low on steep slopes. 4. **Mask correctness** — Antarctica, open ocean and closed-basin sinks are NoData; coastal land cells resolve to data. 5. **Range** — distance within [0, 50 000] m; TWI within [~0, 35]. ## Versioning - `v0.1` — pysheds D8, halo-tiled at 9 arcsec, no cross-tile reconciliation. Captured 2.5–19.3% of true drainage at reference gauges. - `v0.3` — cross-tile propagation in area units; GEDTM30 bare earth; GSHHG-f land mask; closed-basin termination; channel threshold as an area; the Volga/Kama reservoirs as flow nodes; contributing area shipped alongside the derived layers. - Next — boundary-cell DAG in place of iteration; a treatment for floodplain lakes that does not depend on a mapped centerline; 3 arcsec (90 m) regional derivations where the DEM and downstream hazards justify it.