What is flow direction in GIS and how is the D8 method computed?

Flow direction is a raster that records, for each cell, the direction water leaves it. The common D8 method points each cell to the single steepest of its eight neighbors and encodes the result as a power-of-two code (1=E, 2=SE, 4=S, up to 128=NE in the Esri scheme). It is the input flow accumulation is built from.

What does flow accumulation actually show?

Flow accumulation counts, for each cell, how many upstream cells drain into it through the flow-direction grid. Cells with high counts trace valley bottoms and channels; thresholding the accumulation grid extracts a stream network. The count times cell area gives contributing (catchment) area.

Why do you have to fill pits before computing flow accumulation?

A pit is a cell or group of cells with no lower neighbor, often a DEM artifact. Flow routed into a pit has nowhere to go, so accumulation stalls and the drainage network fragments. Filling (or breaching) raises pits to their pour point so water can continue downslope.

What is the difference between D8 and D-infinity flow direction?

D8 sends all flow from a cell to one neighbor, which is simple and robust but produces parallel artificial lines on hillslopes. D-infinity splits flow between the two neighbors that bracket the steepest downslope angle, giving smoother, more realistic dispersion on slopes at higher computational cost.

Flow Direction and Flow Accumulation Explained

Short answer

Flow direction and flow accumulation are the two foundational rasters of hydrologic terrain analysis. Flow direction records, for every DEM cell, where water would leave it; flow accumulation counts how many upslope cells ultimately drain through each cell. Threshold the accumulation grid and you get a synthetic stream network; trace upslope from a point and you get its watershed. Both are extremely sensitive to errors in the input elevation surface, which is why pit filling and CRS/unit checks come before anything else.

The pipeline, in order

Hydrologic derivatives are computed as a chain, and skipping a step corrupts everything downstream. The canonical sequence is:

Condition the DEM (fill or breach pits and flat areas).
Compute flow direction from the conditioned surface.
Compute flow accumulation from flow direction.
Extract streams by thresholding accumulation.
Delineate watersheds / catchments from flow direction and pour points.

You cannot reorder these. Flow accumulation reads the flow-direction grid; flow direction reads the conditioned DEM. Feed a raw, pitted DEM into step 2 and the drainage network will be fragmented no matter how good the later tools are.

How flow direction is computed (D8)

The most widely used algorithm is D8 (deterministic eight-neighbor), introduced by O'Callaghan and Mark in 1984. For each cell, D8 examines its eight neighbors, computes the drop to each (elevation difference divided by distance, where diagonal neighbors are 1.414 cell-widths away), and assigns flow to the single neighbor with the steepest descent.

Esri encodes the result as powers of two so directions can be combined bitwise: 1 = east, 2 = southeast, 4 = south, 8 = southwest, 16 = west, 32 = northwest, 64 = north, 128 = northeast. So a flow-direction value of 4 means "this cell drains south." Other toolchains (GRASS, TauDEM, WhiteboxTools) use their own encodings, so never assume the codes are portable between tools.

D8's strength is simplicity and robustness; its weakness is that all flow from a cell goes to exactly one neighbor. On planar hillslopes this produces visible parallel "flow lines" that look artificial. D-infinity (Tarboton, 1997) addresses this by representing flow direction as a continuous angle and splitting flow proportionally between the two neighbors that bracket that angle. The result disperses flow more realistically across slopes and is preferred for soil-erosion and wetness-index work, at higher cost. MFD (multiple flow direction) methods spread flow to all downslope neighbors.

How flow accumulation is computed

Once every cell has a flow direction, accumulation is a routing problem: follow the directions and tally how many cells drain through each. By convention, most tools report the number of upstream cells excluding the cell itself, so a ridge-top cell has accumulation 0 and a major valley cell can have hundreds of thousands.

To turn the count into physical contributing area, multiply by cell area. For a 10 m DEM, each cell is 100 m², so an accumulation value of 50,000 corresponds to a contributing area of 5,000,000 m² = 5 km². This is the number that matters for hydrologic modelling, not the raw count.

Extracting a stream network is then a threshold on accumulation. A common starting point is "channels begin where contributing area exceeds some value," for example a few hundred cells on a fine DEM. The right threshold is landscape-dependent; in arid terrain channels require more contributing area than in humid uplands, so calibrate against mapped streams rather than copying a default.

Worked example with GDAL, WhiteboxTools, and ArcGIS

Starting from a 10 m DEM (dem.tif) in a projected CRS with meters as both horizontal and vertical units:

Condition the surface (WhiteboxTools breaches depressions, which is often gentler than filling):

whitebox_tools -r=BreachDepressions --dem=dem.tif --output=dem_fill.tif

Compute D8 pointer and accumulation:

whitebox_tools -r=D8Pointer --dem=dem_fill.tif --output=fdir.tif
whitebox_tools -r=D8FlowAccumulation --input=fdir.tif --pntr --output=facc.tif

In ArcGIS Pro the equivalent Spatial Analyst chain is Fill → Flow Direction → Flow Accumulation, then Con("facc" > 500, 1) and Stream to Feature to vectorize. In QGIS the Processing toolbox exposes the same steps through SAGA (Fill Sinks, Flow Accumulation) and GRASS (r.watershed, which conditions, routes, and accumulates in one pass and is a good cross-check).

GDAL's gdaldem does not compute hydrology (it handles slope, aspect, hillshade, roughness, TRI, TPI); use it earlier to confirm your DEM is sane, but reach for Whitebox, SAGA, GRASS, TauDEM, or ArcGIS for routing.

Why pits and NoData break everything

A pit (or sink) is a cell, or cluster of cells, surrounded entirely by higher cells, so it has no downslope neighbor. Some pits are real (karst, quarries), but most in a DEM are artifacts of interpolation, sensor noise, or vertical quantization. Flow routed into a pit terminates there: accumulation stops growing, and the channel downstream of the pit appears disconnected. A handful of single-cell pits can shatter a drainage network into dozens of fragments.

Filling raises each pit to the elevation of its lowest rim cell (its pour point) so flow can escape. Breaching instead carves a shallow channel through the obstruction, which alters fewer cells and avoids creating large artificial flat ponds; many practitioners prefer breaching or a hybrid on detailed LiDAR DEMs.

Flat areas are the second problem: after filling, large flats have no gradient, so flow direction is undefined. Good tools (GRASS r.watershed, Whitebox) impose a small consistent gradient across flats toward the outlet rather than leaving them ambiguous.

NoData holes, bridges, and culverts are the third. A NoData void interrupts routing entirely. A road embankment over a culvert reads as a continuous dam in the DEM and dams the modelled flow even though water passes through the culvert in reality. These require manual hydro-enforcement (burning streams in, or stamping culverts) before routing.

Common pitfalls and their causes

Fragmented stream network. Almost always unfilled pits or NoData. Run the fill/breach step and re-check.
Streams that stop at roads. Embankments acting as dams; hydro-enforce culverts.
Direction codes look wrong in another tool. You assumed Esri's power-of-two scheme but the grid came from GRASS or TauDEM. Check the encoding.
Accumulation values that mean nothing physical. You read raw cell counts instead of multiplying by cell area; and remember resampling a DEM changes cell area and therefore every contributing-area number.
Vertical/horizontal unit mismatch. Elevation in feet on a DEM with meters spacing distorts the slope used to pick steepest descent. Confirm both units before routing.
Over-filling. Aggressive fill on a noisy DEM drowns real basins and creates large false flats; prefer breaching where possible.

QA and validation

Overlay the extracted streams on a hillshade and on an independent mapped hydrography layer (national hydrography dataset, topographic map). Channels should sit in valley bottoms, not on slopes or ridges.
Confirm watershed outlines close and that the total contributing area at the basin outlet matches an independent estimate within a reasonable margin.
Inspect the filled-minus-original difference raster: large or widespread fill depths flag DEM quality problems, not just isolated pits.
Check the edges. Cells near the DEM boundary have truncated upslope areas, so basins clipped by the extent are underestimated; buffer the DEM beyond the area of interest.
Re-run with a second algorithm (D8 vs D-infinity, or GRASS vs Whitebox) on a sample and confirm the major channels agree.

Bathyl perspective

We treat the conditioned DEM, not the raw one, as the real input to hydrology, and we record exactly how it was conditioned (fill vs breach, flat-resolution method, any burned streams) alongside the outputs. A drainage map is only as trustworthy as the depression handling behind it, so we validate every extracted network against independent hydrography before it informs a decision.

Flow Direction and Flow Accumulation Basics