🧭 Concept Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 2 of 20 Updated

Flow Direction Explained: D8, D-Infinity and MFD Compared

Problem statement

Every DEM hydrology workflow turns elevation into flow direction: for each cell, where does its water go? D8 sends all of it to the single steepest of eight neighbours. D-infinity (D∞) sends it along the steepest downslope direction at any angle, split between the two neighbours that bracket it. Multiple flow direction (MFD) shares it among every lower neighbour. The choice barely moves a large river and completely changes a hillslope.

Measured with pysheds on the conditioned 10 m DEM of the Esopus Creek catchment in New York (2,406 × 3,470 cells), and on a synthetic tilted plane:

  • Along channels draining more than 1 km², D∞ accumulation was within 10% of D8 for 82.6% of cells; MFD for 50.2%.
  • On hillslopes, D8 left 10.0% of the basin's cells receiving no flow at all, against 2.1% for D∞ and 0.1% for MFD; median accumulation was 10, 17.8 and 25.3 cells.
  • On a smooth plane sloping 22.5° off the grid axes, D8 accumulation along one row ranged from 1 to 868 cells — parallel streaks — while D∞ and MFD varied by 5% and 1.6%.
  • MFD cost eight times the memory: a 534 MB flow-direction array against 67 MB, and 1.87 s of accumulation against 0.26 s.

Quick answer 

import numpy as np

if not hasattr(np, "in1d"):          # pysheds 0.5 on NumPy 2.4+
    np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
from pysheds.grid import Grid

grid = Grid.from_raster("esopus_breached.tif")
dem = grid.read_raster("esopus_breached.tif")
for routing in ("d8", "dinf", "mfd"):
    fdir = grid.flowdir(dem, routing=routing)
    acc = grid.accumulation(fdir, routing=routing)
    print(routing, np.asarray(fdir).shape, float(acc.max()) / 1e4, "km2")
d8 (2406, 3470) 513.4838 km2
dinf (2406, 3470) 513.4389212565115 km2
mfd (8, 2406, 3470) 497.0677272964447 km2

Use D8 for stream networks, watershed delineation and anything that needs one downstream cell. Use D∞ or MFD for hillslope quantities — wetness indices, erosion, contributing area on slopes — where dispersion is real.

Table comparing how D8 sends flow to one neighbour, D-infinity splits it between two neighbours by angle, and MFD shares it among all lower neighbours, with the memory each needs.
One neighbour, two neighbours, or all lower neighbours: the rest of the differences follow from that.

Step-by-step solution

1. Condition the DEM first

Every method needs a downslope neighbour for every cell. Run the comparison on a depression-free DEM — here a WhiteboxTools least-cost breach of the 3DEP 10 m model — or flow will stop at pits whatever the method; see breaching or filling a DEM.

2. D8: one neighbour, eight choices

D8 compares the drop to each of the eight neighbours, dividing diagonal drops by √2, and sends everything to the steepest. The result is an integer code per cell — in pysheds 64, 128, 1, 2, 4, 8, 16, 32 for north, north-east, east and so on clockwise — and a flow network that never splits. On the Esopus DEM it took 0.11 s and produced a 67 MB array; accumulation took 0.26 s.

3. D∞: any angle, split between two neighbours

D∞ fits planar triangles between the cell and pairs of neighbours, finds the steepest downslope angle, and divides flow between the two neighbours on either side of it in proportion to angular distance. Its output is an angle per cell (float64, 67 MB); accumulation took 1.18 s. Flow can split but never spreads to more than two cells from any one.

4. MFD: all lower neighbours

MFD gives each lower neighbour a share of the flow proportional to a power of its slope. pysheds stores the eight proportions for every cell, an array of shape (8, 2406, 3470) — 534 MB in float64. Accumulation took 1.87 s. Flow spreads widely, which suits divergent hillslopes and exaggerates spreading on valley floors.

5. Compare where it matters: channels

On cells with at least 1 km² of D8 drainage, D∞ accumulation had a median ratio to D8 of 0.999; MFD 0.908. At the Coldbrook gauge the largest accumulation within 150 m was 492.44 km² with D8, 491.98 km² with D∞ and only 408.26 km² with MFD, because MFD spreads the river's flow across neighbouring floodplain cells instead of concentrating it in one.

6. Compare where it matters: hillslopes

On cells with fewer than 100 cells of D8 drainage, D∞ accumulation had a median ratio of 1.38 to D8 and MFD 2.00. D8 left one in ten cells with nothing flowing into them — ridges of the grid rather than of the terrain. The 99th percentile of accumulation was 2,590 cells for D8, 2,917 for D∞ and 4,578 for MFD.

7. Look for straight parallel lines

On a plane whose true downslope bearing falls between two D8 directions, D8 alternates between them in long straight runs. On a 10% slope bearing 157.5°, D8 used only the south and south-east codes, and accumulation along a row across the flow ranged from 1 to 868 cells with a coefficient of variation of 1.694. D∞ ranged from 120 to 150 (CV 0.050) and MFD from 163 to 177 (CV 0.016). Real flats and gentle planes show the same streaks; see fixing streams that run in straight parallel lines.

8. Delineate watersheds with D8, or treat the others as fractions

A D8 catchment is a yes-or-no set of cells. With D∞ and MFD, a cell can send part of its flow into the catchment: pysheds counted 495.35 km² of cells contributing anything with D∞ and 515.49 km² with MFD, against 492.44 km² with D8. For boundaries, use D8; for contributing fractions, use the proportional methods deliberately.

Bar chart of the coefficient of variation of flow accumulation across a row of a smooth tilted plane for D8, D-infinity and MFD.
On a perfectly smooth plane, D8 invents channels; the proportional methods do not.

Code examples

Example 1 — the three methods, timed 

import time

import numpy as np

if not hasattr(np, "in1d"):
    np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
from pysheds.grid import Grid

grid = Grid.from_raster("esopus_breached.tif")
dem = grid.read_raster("esopus_breached.tif")
accumulations = {}
for routing in ("d8", "dinf", "mfd"):
    start = time.perf_counter()
    fdir = grid.flowdir(dem, routing=routing)
    middle = time.perf_counter()
    acc = grid.accumulation(fdir, routing=routing)
    end = time.perf_counter()
    accumulations[routing] = np.asarray(acc, dtype="float64")
    fd = np.asarray(fdir)
    print(f"{routing:4} flowdir {middle - start:5.2f} s, accumulation {end - middle:5.2f} s, "
          f"direction array {fd.shape} {fd.dtype} {fd.nbytes / 1e6:.0f} MB")
d8   flowdir  0.11 s, accumulation  0.26 s, direction array (2406, 3470) int64 67 MB
dinf flowdir  0.30 s, accumulation  1.18 s, direction array (2406, 3470) float64 67 MB
mfd  flowdir  0.31 s, accumulation  1.87 s, direction array (8, 2406, 3470) float64 534 MB

Timings vary between runs; in an earlier run on the same machine MFD accumulation took 4.73 s and D8 0.90 s.

Example 2 — how the methods differ inside the basin 

import geopandas as gpd
from rasterio import features

basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(grid.crs.srs)
inside = features.rasterize(basin.geometry, out_shape=dem.shape, transform=grid.affine).astype(bool)

d8 = accumulations["d8"][inside]
channels, hillslopes = d8 >= 1e4, d8 < 100
for routing, acc in accumulations.items():
    a = acc[inside]
    ratio = a / d8
    print(f"{routing:4} no inflow {np.isclose(a, 1).mean():5.1%}  median {np.median(a):6.2f} cells  "
          f"channels within 10% of D8 {np.mean(np.abs(ratio[channels] - 1) < 0.1):5.1%}  "
          f"hillslope median ratio {np.median(ratio[hillslopes]):.2f}")
d8   no inflow 10.0%  median  10.00 cells  channels within 10% of D8 100.0%  hillslope median ratio 1.00
dinf no inflow  2.1%  median  17.76 cells  channels within 10% of D8 82.6%  hillslope median ratio 1.38
mfd  no inflow  0.1%  median  25.30 cells  channels within 10% of D8 50.2%  hillslope median ratio 2.00

Example 3 — a smooth plane that D8 cannot represent 

from pysheds.view import Raster, ViewFinder

n, bearing = 201, np.deg2rad(22.5)
yy, xx = np.mgrid[0:n, 0:n] * 10.0
plane = 1000 - (yy * np.cos(bearing) + xx * np.sin(bearing)) * 0.1   # 10% slope, 22.5° east of south
view = ViewFinder(affine=grid.affine, shape=plane.shape, crs=grid.crs, nodata=-9999.0)
synthetic = Grid(viewfinder=view)
for routing in ("d8", "dinf", "mfd"):
    fd = synthetic.flowdir(Raster(plane, viewfinder=view), routing=routing)
    row = np.asarray(synthetic.accumulation(fd, routing=routing), dtype="float64")[150, 60:141]
    print(f"{routing:4} accumulation across the flow: min {row.min():6.1f}, max {row.max():6.1f}, CV {row.std() / row.mean():.3f}")
d8   accumulation across the flow: min    1.0, max  868.0, CV 1.694
dinf accumulation across the flow: min  120.0, max  150.0, CV 0.050
mfd  accumulation across the flow: min  163.4, max  176.7, CV 0.016

Explanation

Why D8 makes clean networks and poor hillslopes

A single receiver means every cell belongs to exactly one flow path, so accumulation concentrates into one-cell channels and catchments are unambiguous. On a hillslope, though, water spreads; forcing it into one of eight directions collects it into artificial rills, while the cells between them receive nothing. That is why one cell in ten had no inflow.

Why D∞ is a compromise

Splitting between only two neighbours keeps dispersion limited, so channels stay nearly as concentrated as D8 — within 10% for four in five channel cells — while flow direction is no longer restricted to multiples of 45°. It removes most of D8's streaking on planes and still produces reasonably narrow channels.

Why MFD under-concentrates rivers

MFD shares flow with every lower neighbour, including floodplain cells beside the channel that are only slightly lower. Along a valley bottom that spreads the river over several cells, so no single cell carries the full drainage area: 408 km² at the gauge instead of 492 km². Variants that switch to single-direction flow above a threshold avoid this.

Why memory differs so much

D8 and D∞ need one value per cell; MFD needs eight. On a 8.3-million-cell grid that is the difference between 67 MB and 534 MB before accumulation even starts, and it grows with the square of resolution.

Table of flow direction array memory and accumulation time for D8, D-infinity and MFD on the 10 metre Esopus DEM.
MFD's realism on hillslopes costs eight times the memory and seven times the accumulation time.

Edge cases or notes

  • Flow direction codes differ between tools. pysheds, WhiteboxTools and ESRI use different D8 numbering; pass a dirmap or convert before mixing them.
  • Flats have no steepest direction. Resolve them before computing flow direction, or D8 codes them −1.
  • Edge cells drain off the grid; accumulation near the edge of a clipped DEM is truncated.
  • MFD exponents control dispersion; larger exponents concentrate flow towards D8.
  • Stream burning forces directions along a mapped network; see burning a known river network into a DEM.
  • Topographic wetness indices are sensitive to the method on hillslopes; report which one you used.
  • Lower resolution reduces streaking but also removes real small channels.

FAQ

What is the D8 flow direction algorithm?

It sends all of a cell's flow to the steepest of its eight neighbours, dividing diagonal drops by √2. It gives a single downstream cell everywhere, which makes clean stream networks and unambiguous watersheds.

What is the difference between D8 and D-infinity?

D-infinity allows any flow angle and splits flow between the two neighbours bracketing it. On Esopus Creek channels the two agreed within 10% for 82.6% of cells; on hillslopes D-infinity accumulation was a median 1.38 times D8.

When should I use MFD flow direction?

For hillslope processes where flow genuinely spreads, such as wetness indices. It under-concentrates rivers — 408 km² at the gauge instead of 492 km² — and needs eight times the memory of D8.

Why does D8 create parallel lines on smooth slopes?

When the true downslope direction falls between two of the eight directions, D8 alternates between them in straight runs. On a smooth plane accumulation varied from 1 to 868 cells across a row.

Which flow direction method should I use for watershed delineation?

D8. Its catchments are sets of cells with a clear boundary; D-infinity and MFD give partial contributions that need a threshold to become a boundary.

Do pysheds and WhiteboxTools use the same D8 codes?

No. pysheds uses 64, 128, 1, 2, 4, 8, 16, 32 for north through north-west; WhiteboxTools uses powers of two starting from north-east. Convert or pass a dirmap when combining them.