Short answer

Watershed delineation extracts the upslope contributing area that drains to a chosen outlet by chaining five hydrology operations: fill sinks → flow direction → flow accumulation → stream definition → pour-point snap → watershed. Skip the fill, and your flow network fragments at every spurious pit. Skip the pour-point snap, and your basin will be a handful of cells. Get the order and the conditioning right and the result is reproducible from any hydrologically conditioned DEM.

The pipeline, step by step

1. Hydrologically condition the DEM (fill or breach)

Raw DEMs contain depressions — some real (sinkholes, quarries), most artefacts of interpolation or sensor noise. A single one-cell pit will trap flow and stop the drainage network from being continuous. Two strategies:

  • Filling raises pit cells to the level of their lowest pour point (ArcGIS Fill, QGIS/GRASS r.fill.dir, WhiteboxTools FillDepressions).
  • Breaching carves a channel out of the depression, which alters fewer cells and avoids creating large flat artificial lakes; BreachDepressionsLeastCost in WhiteboxTools is often preferable on flat or LiDAR-derived terrain.

A pure fill on a flat coastal DEM can manufacture enormous flat pools; breaching usually gives a more physical result there.

2. Flow direction

This assigns each cell the direction water leaves it. The default is D8: flow goes entirely to the single steepest of the eight neighbours, encoded as a power-of-two bitfield (1=E, 2=SE, 4=S, 8=SW, 16=W, 32=NW, 64=N, 128=NE in ESRI's scheme). D8 is simple and gives crisp channels but suffers parallel-flow striping on smooth hillslopes. D-infinity (Tarboton) splits flow proportionally between the two neighbours that bracket the steepest descent direction, which represents divergence on hillslopes far better. For catchment outlines D8 is standard; for erosion or wetness indices, D-infinity is usually the better physics.

3. Flow accumulation

Each cell receives the count (or weighted sum) of all cells draining through it. Channel cells light up as high-accumulation lines; ridgelines stay near zero. This raster is both your stream-extraction input and the surface you snap pour points to.

4. Stream definition

Threshold the accumulation raster: cells above a chosen accumulation are "stream." The threshold sets drainage density and is a real decision, not a default — on a 10 m DEM a threshold of 1,000 cells means roughly 0.1 km² of contributing area initiates a channel. Calibrate it against a known hydrography layer (e.g. national stream data) rather than guessing.

5. Snap the pour point and delineate

The outlet you digitise rarely lands exactly on a channel cell. If it sits one cell off the channel, the watershed tool returns the tiny local catchment of that off-stream cell. Snap the pour point to the nearest high-accumulation cell within a tolerance (a few cells) first, then run the catchment tool.

Worked example: QGIS + GRASS (r.watershed) and WhiteboxTools

QGIS ships the GRASS hydrology tools in its Processing toolbox. The cleanest single-pass approach is r.watershed, which conditions internally:

Processing ▸ GRASS ▸ r.watershed
  Elevation: dem_10m
  Outputs: drainage (direction), accumulation, stream
  Threshold: 1000     # min accumulated cells to start a stream

Then for a specific outlet:

Processing ▸ GRASS ▸ r.water.outlet
  Drainage direction: drainage (from r.watershed)
  Coordinates of outlet: 482300,5630100
  → basin raster, then Raster ▸ polygonize

A discrete, explicit ESRI-style pipeline (e.g. in WhiteboxTools, callable from QGIS) looks like:

FillDepressions  --dem=dem_10m.tif --output=filled.tif
D8Pointer        --dem=filled.tif  --output=fdir.tif
D8FlowAccumulation --input=fdir.tif --pntr --output=facc.tif
ExtractStreams   --flow_accum=facc.tif --threshold=1000 --output=streams.tif
JensonSnapPourPoints --pour_pts=outlet.shp --streams=streams.tif \
                     --snap_dist=50 --output=snapped.shp
Watershed        --d8_pntr=fdir.tif --pour_pts=snapped.shp --output=basin.tif

JensonSnapPourPoints with a 50 m snap distance is the step that prevents the empty-basin problem. The ArcGIS equivalent is the Spatial Analyst sequence Fill → Flow Direction → Flow Accumulation → Snap Pour Point → Watershed.

Common pitfalls and why they happen

  • No fill/breach. The most common error: flow direction on an unconditioned DEM produces a network full of internal sinks, and downstream accumulation is broken into disconnected pieces.
  • Pour point off the channel. Because the outlet was digitised by eye over a basemap, it lands beside the modelled stream. Always snap.
  • Geographic CRS. Running hydrology in EPSG:4326 means cell sizes are in degrees and accumulation areas are distorted by latitude. Reproject to a projected CRS first: gdalwarp -t_srs EPSG:32631 dem_ll.tif dem_utm.tif.
  • Over-filling flat terrain. A blanket fill on low-relief DEMs creates artificial flat lakes that route flow incorrectly; prefer breaching there.
  • Edge basins. Catchments touching the DEM border are truncated — their true contributing area extends off-tile. Flag or clip them; never report a border basin's area as complete.
  • Threshold mistaken for truth. The stream threshold controls how many channels appear. Two analysts with different thresholds get different basins; document the value.

QA and validation

  • Overlay the modelled streams on reference hydrography. Tune the accumulation threshold until drainage density matches the known network.
  • Check the divides against contours or a hillshade. Basin boundaries should follow ridgelines; a divide cutting across a valley signals an unfilled sink or a CRS problem.
  • Sum the area and sanity-check it against gauged drainage area if a stream gauge exists at your outlet — agreement within a few percent is a strong sign the conditioning is sound.
  • Confirm units and CRS with gdalinfo before starting; hydrology only behaves in projected, metric space.

Bathyl perspective

A delineated watershed is the product of a specific conditioning recipe — fill versus breach, D8 versus D-infinity, the accumulation threshold, the snap distance. We record every one of those values alongside the basin polygon, because the boundary is meaningless to a reviewer who cannot see how it was derived. Reproducibility here is not bureaucracy; it is the difference between a defensible contributing area and a pretty blue line.

Related reading

Sources