Short answer

Drainage and terrain risk screening turns a digital elevation model into a small set of derivative rasters — flow accumulation, slope, and a wetness index — and thresholds them to flag where surface water concentrates, where it ponds, and where steep banks meet channels. The workflow is: condition the DEM so water has a continuous downslope path, route flow, derive the indices, and overlay the flagged cells against assets and proposed works.

It is deliberately a screening product. It tells you where to look and where to commission a hydrologist or geotechnical engineer. It does not replace a HEC-RAS hydraulic run, a rainfall-runoff model, or a site geotechnical assessment, and the deliverable should say so explicitly.

What the DEM actually encodes

Every output below is computed from the elevation grid and its immediate neighbours, so the source DEM governs everything. Three properties matter most.

  • Cell size. A 1 m lidar-derived DEM resolves a roadside ditch; a 30 m SRTM or Copernicus GLO-30 grid does not. Drainage features narrower than two or three cells are invisible. Match the grid to the feature you are screening for.
  • Vertical accuracy and noise. Lidar DEMs are typically accurate to ~0.1-0.3 m vertically over open ground; global radar DEMs are several metres. Noise propagates into slope and is amplified in curvature.
  • Surface type. A digital surface model (DSM) includes buildings and canopy; a digital terrain model (DTM) is the bare earth. Route flow on a DTM. Routing on a DSM sends water over rooftops and around tree crowns, producing nonsense channels.

Confirm the horizontal CRS is projected (metres), not geographic degrees. Flow and slope algorithms assume planimetric cell dimensions; running them on EPSG:4326 produces meaningless gradients because a degree of longitude is not a degree of latitude in ground distance.

Step 1 — Hydrologically condition the DEM

Raw DEMs contain depressions: real ones (kettle holes, quarries) and spurious ones from interpolation or noise. Left in place, they trap routed flow and break the accumulation network into disconnected fragments. You have two repair strategies:

  • Fill raises pit cells to the level of their lowest pour point. Simple, but it flattens real basins and can create long artificial flats.
  • Breach carves a shallow channel from the pit to the downstream cell. It preserves drainage volume better and is preferred where genuine depressions exist.

In QGIS this is the SAGA or WhiteboxTools provider — for example WhiteboxTools BreachDepressionsLeastCost, or the SAGA "Fill Sinks (Wang & Liu)" algorithm, which fills and assigns flow directions in one pass. In raw GDAL there is no single fill tool; most teams use RichDEM or Whitebox for the conditioning step, then continue in GDAL/QGIS.

# WhiteboxTools, command line
whitebox_tools -r=BreachDepressionsLeastCost \
  --dem=dtm.tif --output=dtm_breached.tif --dist=100

Keep both the raw and conditioned DEMs. The difference raster (gdal_calc.py -A dtm_breached.tif -B dtm.tif --calc="A-B") is itself a useful diagnostic: large positive values show where conditioning moved a lot of elevation, which is exactly where to be sceptical of downstream results.

Step 2 — Flow direction and accumulation

Flow direction assigns each cell the direction water leaves it. The classic D8 method sends all flow to the single steepest of eight neighbours; D-infinity (Tarboton) splits flow between the two downslope neighbours bracketing the steepest gradient, which represents divergent hillslopes more realistically. Use D8 for crisp channel networks, D-infinity for hillslope wetness.

Flow accumulation then counts, for each cell, how many cells drain through it (or sums their area). High-accumulation cells are channels; the value scales with contributing area.

# QGIS Processing: SAGA "Flow Accumulation (Top-Down)"
# or WhiteboxTools:
whitebox_tools -r=D8FlowAccumulation \
  --input=dtm_breached.tif --output=facc.tif --out_type=cells

To convert accumulation in cells to contributing area in square metres, multiply by cell area: a 1 m grid with accumulation 50,000 represents 5 ha upslope.

Step 3 — Slope, curvature, and the wetness index

Compute slope from the conditioned DEM in degrees:

gdaldem slope dtm_breached.tif slope_deg.tif -compute_edges

The Topographic Wetness Index combines contributing area and slope:

TWI = ln( a / tan(β) )

where a is the upslope contributing area per unit contour width (accumulation × cell size) and β is local slope in radians. High TWI marks cells that receive large drainage on shallow gradients — saturation-prone hollows, footslopes, and floodplain margins. Build it with gdal_calc.py, guarding against a zero-slope divide:

gdal_calc.py -A facc.tif -B slope_deg.tif \
  --calc="log((A*1.0+1)/tan(radians(maximum(B,0.5))))" \
  --outfile=twi.tif

The maximum(B,0.5) floor on slope prevents infinite TWI on perfectly flat cells, a common artefact of over-filled DEMs.

Step 4 — Threshold and overlay

Screening is the act of turning continuous rasters into flagged zones. Reasonable starting thresholds, to be calibrated locally:

  • Concentrated flow / potential channel: flow accumulation above a contributing-area threshold (often a few hectares; smaller in steep humid terrain, larger in arid flats).
  • Ponding / saturation: TWI in roughly the top decile of the catchment, and slope below ~2-3 degrees.
  • Erosion / bank instability: slope above ~25-30 degrees within one or two cells of a high-accumulation channel.

Reclassify each raster (QGIS "Reclassify by table" or gdal_calc.py) to 0/1 flags, then combine. Overlay the flagged raster against buildings, roads, intakes, or a proposed alignment and clip to the area of interest. The output is a heat map of "investigate here", not a design surface.

Worked example — culvert siting screen

A site team needs to flag where an access track crosses concentrated flow. On a 1 m DTM: breach depressions, run D8 accumulation, threshold accumulation at 10,000 cells (1 ha), and rasterise the track centreline. Intersect the buffered track (ST_Buffer in PostGIS, or QGIS Buffer at 5 m) with the channel flags. Each intersection is a candidate culvert location. Sort by upslope contributing area to rank them — the cell with 80 ha upslope needs a larger structure than one with 3 ha. None of this sizes the culvert; it tells the engineer where to size one.

Common pitfalls and why they happen

  • Routing on a DSM. Canopy and buildings divert flow. The channels look plausible but follow infrastructure, not ground. Always use a bare-earth DTM.
  • Skipping conditioning. Without filling/breaching, accumulation breaks into short disconnected segments at every spurious pit, and TWI spikes on artificial flats. The network simply looks wrong.
  • Over-filling. Aggressive fill creates large artificial flats with near-zero slope, producing TWI infinities and false ponding flags. Prefer breaching where real depressions exist.
  • Slope from a geographic CRS. Degrees of latitude and longitude differ in ground length, so gradients are distorted and direction-dependent. Reproject to a projected CRS first.
  • Treating thresholds as universal. A 1 ha channel-initiation threshold tuned for humid uplands flags far too much in arid plains. Calibrate against mapped streams or aerial imagery.

Quality checks

  • Overlay the derived channel network on aerial imagery or a national hydrography layer; they should broadly agree along main valleys.
  • Verify total catchment accumulation at the outlet equals the cell count of the basin — a sanity check that flow routing is connected.
  • Inspect the DEM edges: -compute_edges aside, accumulation truncates where the basin extends beyond the tile, so buffer the extent or note the edge as unreliable.
  • Confirm CRS, cell size, DTM-vs-DSM, and source date are recorded in the layer metadata.

Bathyl perspective

We treat drainage screening as a triage layer that points specialists at the right places fast. The deliverable always names its DEM source and resolution, the conditioning method, and the thresholds used, so a reviewer can reproduce or challenge any flag — and so no one mistakes a screening overlay for a flood model.

Related reading

Sources