River Network Formation: How Water Carves Fractal Drainage Basins

Eight directions, one rule

A digital elevation model is a grid of heights. To turn it into a river network you first need a rule for which way water flows off every cell, and the rule almost everyone uses is D8: each cell drains into whichever of its eight neighbours gives the steepest downhill slope. No hydraulics, no momentum, just the direction of steepest descent, computed once per cell.

slope(neighbour) = (elev[cell] - elev[neighbour]) / distance
distance = 1        for the 4 orthogonal neighbours
distance = sqrt(2)  for the 4 diagonal neighbours

flowDir[cell] = argmax over the 8 neighbours of slope(neighbour)

D8 is crude at the scale of a single cell — real water spreads across a hillslope, not down one thread — but averaged over a whole basin it reproduces the branching, tree-like pattern of real drainage remarkably well. Pits (local minima with no downhill neighbour) have to be filled or breached first, or the flow simply stops and pools.

Flow accumulation: counting upstream area

Once every cell has a flow direction, flow accumulation counts, for each cell, how many upstream cells eventually drain through it. Process cells from highest to lowest elevation (a topological sort along the flow-direction graph) and add each cell's accumulated flow to whatever it drains into:

accum[cell] = 1                      // itself
for cell in cells sorted by elevation, highest first:
    downstream = flowDir[cell]
    accum[downstream] += accum[cell]

// accum[cell] is now the contributing (drainage) area, in cell counts

High accumulation values trace exactly the shape of a river network — thresholding accum at some value (say 100 cells) and keeping only cells above it is the standard way to extract a channel network from a raw elevation grid. This is the same computation QGIS, ArcHydro and TauDEM run on real satellite terrain.

The stream power law: where and how fast it erodes

Flow accumulation gives geometry; the stream power law gives dynamics. It states that the rate a channel cuts down into rock is proportional to a power of the drainage area (a proxy for discharge) times a power of the local slope:

dz/dt = U - K * A^m * S^n

U = rock uplift rate         K = erodibility (rock type, climate)
A = upstream drainage area   S = local channel slope
m ≈ 0.5,  n ≈ 1              (typical calibrated exponents)

This one equation explains why headwater streams are steep and slow-eroding while a big river with a huge catchment can grind through bedrock on a gentle slope: area substitutes for the water and sediment doing the cutting. Rearranged for a channel in equilibrium with uplift, it also predicts the concave-up long profile every real river has.

Horton's laws and Strahler stream order

Robert Horton and later Arthur Strahler noticed the branching itself is statistically regular. Give every headwater stream order 1; when two order-n streams meet, the result is order n+1 (two streams of different order leave the higher one unchanged). Horton's laws then say the number, average length and average drainage area of streams change by a roughly constant ratio from one order to the next — a bifurcation ratio near 3–5 holds across rivers on every continent, which is a large part of why river networks look fractal at any zoom level.

Frequently asked questions

Try it live

Everything above runs in your browser — open River Network Formation and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.