Hydrological Flow as Optimization
modelling Level 4

Hydrological Flow as Optimization

Water always flows downhill, following the path of steepest descent. Given a digital elevation model, we can trace every raindrop from where it lands to where it exits the landscape. This model shows how local gradients generate global drainage networks — a beautiful example of optimization in nature.

Prerequisites: optimization, graph traversal, local minima, flow accumulation

Updated 13 min read

1. The Question

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.

2. The Conceptual Model

Flow as a Graph Problem

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, Global Pattern

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).

3. Building the Mathematical Model

Step 1: Flow Direction (D8 Algorithm)

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.

Step 2: Flow Accumulation

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.

Step 3: Watershed Delineation

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).

Step 4: Stream Network Extraction

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.

4. Worked Example by Hand

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) Determine the flow direction for cell (2, 2) using D8.
(b) If every cell contributes 1 unit of flow, what is the flow accumulation at cell (4, 4)?

Solution

(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.)

5. Computational Implementation

Below is an interactive flow routing visualization.

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.

6. Interpretation

Stream Network Threshold

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).

Watershed Area

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\]

Sinks and Depressions

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).

7. What Could Go Wrong?

Flat Areas

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

Parallel Flow on Uniform Slopes

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.

Spurious Sinks from Noise

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).

Edge Effects

Flow at the edge of a DEM exits the grid. These cells need special handling (treat as outlets or mark as no-data).

8. Extension: Real-World Applications

Flood Prediction

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)

Erosion and Sediment Transport

Flow accumulation is a proxy for stream power — the ability to erode and transport sediment.

High accumulation + steep slope = high erosion potential.

Pollution Routing

Track contaminants (agricultural runoff, industrial spills) as they move downslope through the drainage network.

Ecosystem Connectivity

Flow networks define hydrologic connectivity — which wetlands, streams, and riparian zones are linked. Important for fish migration, amphibian breeding, nutrient cycling.

9. Math Refresher: Graphs and Traversal

Graph Terminology

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.

Traversal Algorithms

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.

Connectivity

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.

Summary

  • Water flows downhill along the path of steepest descent
  • D8 algorithm: Each cell flows to the steepest of 8 neighbors
  • Flow accumulation: Count how many cells drain through each location
  • High accumulation → streams; low accumulation → ridges
  • Watersheds: Areas draining to the same outlet
  • Sinks and flat areas require special handling
  • Flow routing underlies flood modelling, erosion prediction, and watershed delineation

References