From local slopes to global watersheds — routing water across a DEM
2026-02-26
The Bow River draining east from the Canadian Rockies through Calgary did not choose its course. Every metre of its channel is the accumulated record of billions of individual decisions made by individual raindrops: at each point on the landscape, water moves to the lowest available neighbour. The river network is what happens when that single rule is applied, by every drop, over geological time. The tributary structure, the meander pattern, the location of the watershed divide at the Continental Divide — all of it emerges from local optimisation repeated at enormous scale.
This insight turns watershed delineation from a problem requiring complex physical theory into a problem of following a simple rule on a grid. Given a digital elevation model, applying the rule — water flows to the steepest available downslope neighbour — from every cell generates a complete directed flow graph. From that graph you can read off flow accumulation (how many cells drain through each point), stream channels (where accumulation exceeds a threshold), watershed boundaries (the ridgelines that separate contributing areas), and outlet locations. The mathematics is graph theory; the geography is hydrology. This model traces the connection between local gradient and global drainage structure.
Where does the water go?
A raindrop lands on a hillslope. Gravity pulls it downhill. At each point, it moves in the direction of steepest descent — the direction where elevation drops fastest. It continues until it reaches a stream, a lake, or the edge of the DEM.
If we know the elevation at every point, can we predict: 1. Flow direction — which way does water move from each cell? 2. Flow accumulation — how much water passes through each cell? 3. Watersheds — which areas drain to the same outlet? 4. Stream networks — where do channels form?
All four questions are answered by the same fundamental algorithm: routing flow according to the gradient.
A digital elevation model is a grid of cells. Each cell can pass water to one or more downslope neighbors.
Graph representation: - Nodes: Grid cells - Edges: Flow connections from one cell to a downslope neighbor - Directed: Water flows in one direction (downhill)
Flow routing algorithms: 1. D8 (Deterministic 8-direction): Flow goes to the steepest of 8 neighbors 2. D-infinity: Flow direction can be any angle (not constrained to 8 directions) 3. Multiple flow direction (MFD): Flow splits among all downslope neighbors
This model focuses on D8 — the simplest and most widely used.
Local rule: At each cell, send water to the steepest downslope neighbor.
Global pattern: This creates a dendritic network — a branching tree structure where small tributaries merge into larger streams.
Key insight: Complex drainage networks emerge from a simple local optimization (follow the steepest descent).
For each cell (i, j), examine the 8 neighbors:
[i-1,j-1] [i-1,j] [i-1,j+1]
[i,j-1] [i,j] [i,j+1]
[i+1,j-1] [i+1,j] [i+1,j+1]For each neighbor n, compute the slope from (i, j) to n:
s_n = \frac{z_{i,j} - z_n}{d_n}
Where: - z_{i,j} is the elevation of the center cell - z_n is the elevation of neighbor n - d_n is the distance to neighbor n: - Cardinal neighbors (N, S, E, W): d = \Delta x (cell size) - Diagonal neighbors (NE, SE, SW, NW): d = \Delta x \sqrt{2}
Flow direction: The neighbor with the maximum positive slope (steepest descent).
If all neighbors are upslope or equal elevation, the cell is a sink (local minimum) — water pools here.
Once flow direction is determined for every cell, trace water upstream to count how many cells drain through each location.
Algorithm: 1. Initialize accumulation to 1 for every cell (each cell contributes itself) 2. Process cells in order from highest to lowest elevation 3. For each cell, add its accumulation to the accumulation of its downstream neighbor
Result: Cells with high accumulation values are streams — many upstream cells drain through them.
A watershed (or catchment) is the area that drains to a specific outlet.
Algorithm: 1. Choose an outlet cell (e.g., where a stream exits the DEM) 2. Trace upstream from the outlet, marking all cells that flow (directly or indirectly) to it 3. The marked cells form the watershed
Multiple watersheds: Divide the landscape by tracing from multiple outlets. Cells on ridges are watershed boundaries (divides).
Cells with accumulation above a threshold are classified as streams.
Threshold: e.g., accumulation > 100 cells (0.1 km² if cells are 30 m × 30 m).
Network topology: Streams connect where tributaries meet. The flow direction graph defines the network structure.
Problem: A 5×5 DEM with cell size 10 m:
Col 0 Col 1 Col 2 Col 3 Col 4
Row 0: 20 19 18 17 16
Row 1: 19 18 17 16 15
Row 2: 18 17 16 15 14
Row 3: 17 16 15 14 13
Row 4: 16 15 14 13 12(a) Flow direction from cell (2, 2)
Cell (2, 2) has elevation z = 16 m.
Eight neighbors:
| Neighbor | Row | Col | Elevation | Distance | Slope |
|---|---|---|---|---|---|
| NW | 1 | 1 | 18 | 14.14 m | (16-18)/14.14 = -0.141 (uphill) |
| N | 1 | 2 | 17 | 10 m | (16-17)/10 = -0.1 (uphill) |
| NE | 1 | 3 | 16 | 14.14 m | (16-16)/14.14 = 0 (flat) |
| E | 2 | 3 | 15 | 10 m | (16-15)/10 = 0.1 |
| SE | 3 | 3 | 14 | 14.14 m | (16-14)/14.14 = 0.141 ✓ |
| S | 3 | 2 | 15 | 10 m | (16-15)/10 = 0.1 |
| SW | 3 | 1 | 16 | 14.14 m | (16-16)/14.14 = 0 |
| W | 2 | 1 | 17 | 10 m | (16-17)/10 = -0.1 (uphill) |
Steepest descent: SE (slope = 0.141)
Flow direction from (2, 2): Southeast → cell (3, 3)
(b) Flow accumulation at (4, 4)
Cell (4, 4) is the lowest point (elevation 12 m) — all other cells ultimately drain here.
Trace upstream to count contributing cells:
Every cell in the grid flows toward (4, 4) along the gradient. The entire 5×5 grid = 25 cells contribute.
Flow accumulation at (4, 4): 25
(This assumes uniform diagonal flow pattern. In reality, D8 creates a specific tree structure, but for this uniformly sloping DEM, all paths lead to the corner.)
Below is an interactive flow routing visualization.
<label>
Terrain type:
<select id="terrain-flow-select">
<option value="plane">Tilted plane</option>
<option value="valley">Valley</option>
<option value="ridge">Ridge</option>
<option value="cone">Cone (radial)</option>
<option value="crater">Crater (sink)</option>
</select>
</label>
<label>
Display mode:
<select id="display-flow-select">
<option value="elevation">Elevation</option>
<option value="flow-direction">Flow direction (arrows)</option>
<option value="accumulation">Flow accumulation</option>
</select>
</label>
<label>
Grid size:
<input type="range" id="grid-flow-slider" min="20" max="50" step="10" value="30">
<span id="grid-flow-value">30</span> × 30
</label><p id="flow-stats"></p>Try this: - Plane: Uniform flow direction; accumulation increases downslope - Valley: Flow converges to valley axis; high accumulation along the channel - Ridge: Flow diverges from ridgeline; low accumulation on the crest - Cone: Radial flow from peak; accumulation low everywhere - Crater: Flow converges to the center (sink); high accumulation in the crater
Key insight: High accumulation → streams. Low accumulation → ridges or hillslopes.
Question: At what accumulation value does a “stream” begin?
Answer: Depends on climate, geology, vegetation. Typical thresholds: - Humid regions: 0.1–1 km² contributing area - Arid regions: 10–100 km² contributing area
If cell size = 30 m, then 1 km² = 1,111 cells.
Stream order: Classify streams by their position in the network (Strahler ordering, Shreve ordering).
Sum the number of cells in a watershed, multiply by cell area:
A_{\text{watershed}} = N_{\text{cells}} \times (\Delta x)^2
Example: 5000 cells, 30 m cell size:
A = 5000 \times (30)^2 = 4,500,000 \text{ m}^2 = 4.5 \text{ km}^2
Real DEMs often contain sinks — cells lower than all neighbors (artifacts or real depressions like lakes).
Problem: Water accumulates in sinks and can’t flow out, breaking the drainage network.
Solution: Fill sinks — raise sink elevations to the lowest pour point (the lowest cell on the perimeter of the depression).
If multiple cells have identical elevation, slope is zero in all directions — flow direction is ambiguous.
Solutions: - Add small random perturbations to elevations - Use flow direction from surrounding higher cells to infer flow across flats - Specialized flat-routing algorithms
D8 forces flow into discrete 45° sectors. On a uniformly tilted plane, flow should be perpendicular to contours, but D8 approximates this with stepped paths.
D-infinity (infinite-direction) algorithms handle this better but are more complex.
Low-quality DEMs have measurement errors that create artificial pits.
Solution: Smooth the DEM before flow routing (but not so much that real features disappear).
Flow at the edge of a DEM exits the grid. These cells need special handling (treat as outlets or mark as no-data).
Given rainfall intensity and duration, route water across the DEM to predict: - Peak discharge at stream gauges - Inundation extent (which areas flood) - Time to peak (how long until floodwaters arrive)
Flow accumulation is a proxy for stream power — the ability to erode and transport sediment.
High accumulation + steep slope = high erosion potential.
Track contaminants (agricultural runoff, industrial spills) as they move downslope through the drainage network.
Flow networks define hydrologic connectivity — which wetlands, streams, and riparian zones are linked. Important for fish migration, amphibian breeding, nutrient cycling.
Graph: A set of nodes (vertices) connected by edges.
Directed graph: Edges have direction (A → B, not B → A).
Tree: A connected graph with no cycles.
Depth-first search (DFS): Follow one path as far as possible before backtracking.
Breadth-first search (BFS): Explore all neighbors before moving to the next level.
For flow accumulation, we process cells in topological order (upstream to downstream), which requires sorting by elevation.
Two cells are hydrologically connected if there’s a flow path from one to the other.
A watershed is a connected component of the flow graph rooted at an outlet.