Depressions, flow directions and flow accumulation

In this tutorial you take a real elevation model of a forested, mountainous catchment in south-western Oregon and ask the oldest question of terrain hydrology: where does the water go? You will remove the small depressions that trap water in every elevation model, by filling them and by breaching them, give every cell the direction its water flows to, and add up, for every cell, the area of the hillslopes that drain through it. At the end you will have a map of flow accumulation in which the valleys and streams appear on their own, you will know how the five flow methods of Terraflux differ and when to choose each, and you will have weighted the accumulation to ask how much cleared land lies upstream of every reach of the streams. Parts 2 and 3 build streams, watersheds, the height above the nearest drainage and the wetness indices on top of this model.

The idea#

A digital elevation model (DEM) is a raster: a grid of square cells, each holding the height of the ground. Water on the ground runs downhill, so the heights alone say where rain that falls on a cell will go: to a lower neighbour, then to its lower neighbour, and so on, until it leaves the map or stops in a hollow. Following those paths cell by cell is what terrain hydrology does. O'Callaghan & Mark (1984) showed how to extract a whole drainage network from a DEM this way, and the method is now the first step of almost every analysis of catchments, streams, soil moisture and erosion (Moore et al., 1991).

It takes three steps:

1. Condition the DEM. Real elevation models contain depressions: cells, or groups of cells, lower than all their neighbours, from which water cannot leave. A few are real (a lake, a sinkhole, a quarry), but most are errors: noise of the sensor, a bridge or a road embankment across a valley, trees that hide the valley floor. Each one breaks the flow paths, so they are removed first.
2. Give each cell a flow direction: which neighbour (or neighbours) its water goes to.
3. Accumulate: for each cell, count the cells whose water passes through it, itself included. That count, times the area of a cell, is the cell's contributing area (or catchment area): the size of the hillslope and valley above it.

The figure shows the three steps on a toy valley of 5 × 5 cells that drains south, out of the bottom edge. The cell in orange is a pit, 9 m high, lower than all its neighbours.

Enlarge

On the toy, the pit (code 0) collects the water of 20 cells and the outlet receives only 5. One false pit can cut a river in two. That is why every analysis starts by removing depressions, and there are two ways to do it.

Filling raises each depression to the height at which it would overflow, its spill point, as if it filled with water. Barnes, Lehman & Mulla (2014a) gave the fastest exact method, Priority-Flood: it starts from the cells of the map's edge (where water can leave), keeps the cells still to visit in a queue ordered by height, and always takes the lowest one; any neighbour lower than the cell it is reached from is raised to that cell's height. Every cell is visited once. Filling leaves flat areas where the depressions were, and a flat area has no downhill direction; Barnes, Lehman & Mulla (2014b) also gave the method Terraflux uses to drain flats, towards the lower terrain at their outlet and away from the higher terrain around them.

Breaching cuts a channel through the rim of a depression instead, down to a cell lower than its bottom, as a river cuts through a road embankment through its culvert. Lindsay & Dhun (2015) and Lindsay (2016) chose, for each depression, the channel of least cost, the one that changes the elevation least, and showed that breaching alters the DEM less than filling and keeps the flow paths closer to the real ones in landscapes full of roads and embankments.

The profile below follows the main stream of the study area for 3.8 km. The elevation model (black) is full of dips: where the stream crosses them, filling (blue) draws flat lakes over them, while breaching (orange) cuts the rims and keeps the stream close to the original ground.

Enlarge

When to use it#

Flow directions and accumulation are rarely the end of an ecological analysis, but almost everything about water in a landscape is built on them:

• Streams and catchments. Cells with a large accumulation are streams; the cells draining to a point are its catchment (part 2). Where no good map of small streams exists, these are often the best one available.
• What happens upstream. Stream ecologists relate the fish, invertebrates and water quality of a reach to the land use of its whole catchment, not only of its banks (Allan, 2004). A weighted accumulation sums any quantity over the catchment of every cell: the area of cleared land, of a soil type, of farms.
• Soil moisture, wetness and erosion. The wetness, stream power and sediment transport indices (part 3) combine the contributing area of each cell with its slope.
• The routing of models. InVEST's sediment and water models, and SWAT+ (part 4), route water and sediment along the same directions.

What it is not: a hydrological model. Flow accumulation is a property of the terrain only. It assumes that every cell produces the same water, that water follows the surface (not the bedrock, not pipes or ditches), and it says nothing about how much water flows or when. For volumes and timing, use a watershed model such as SWAT+ (part 4).

The alternatives: for the shape of the terrain itself (slope, curvature, topographic position), see Terrain from an elevation model. For a mapped stream network you trust, you may not need to extract one; you will still need the directions to delineate catchments.

What you need#

• An elevation model, as a GeoTIFF or any format GDAL reads. Here: orelevation.tif, the SRTM elevation (1 arc-second, about 30 m) of a forested area in the Oregon Coast Range foothills, the example data of landschaft (in Terraflux's repository, extern/landschaft/src/landschaft/surface/data/). Save it as elevation.tif next to your model. It is the same file as in Terrain from an elevation model.
• Projected coordinates, ideally. The hydro.* operators take the size of the cells of a map in longitude and latitude in metres at its centre latitude, so they run on it, but streams snapped to outlets, the index of connectivity of part 3 and every distance read from the results are simpler in metres. The model therefore starts by reprojecting the DEM to UTM zone 10 north (EPSG:32610) in square cells of 30 m, 369 rows by 270 columns, from 440 m to 1,388 m above sea level.
• Optionally a land cover map for the weighted accumulation. Here: the summer NDVI (a vegetation index, near 1 over dense green canopy) of the same area, orforest.tif in the same folder, saved as ndvi.tif.
• A margin. A catchment cut by the edge of the map is incomplete: its accumulation is too small. Read a DEM that covers the whole catchment of the streams you care about, with a margin.

Enlarge

Build the model#

Read, reproject and fill#

format = "terraflux/1"

[[nodes]]
id = "dem"
op = "raster.read"
params = { path = "elevation.tif" }
note = "SRTM elevation, 1 arc-second."

[[nodes]]
id = "elevation"
op = "raster.reproject"
inputs = { raster = "dem.raster" }
params = { crs = "EPSG:32610", resolution = 30, extent = [500010, 4760850, 508110, 4771920], resampling = "bilinear" }
note = "Square cells of 30 m."

[[nodes]]
id = "filled"
op = "hydro.fill"
inputs = { dem = "elevation.raster" }
note = "Every depression raised to its spill point."

[[nodes]]
id = "save_filled"
op = "raster.write"
inputs = { value = "filled.map" }
params = { path = "filled.tif", type = "f64" }

The reprojection uses bilinear resampling, a weighted mean of the four nearest cells, because elevation is continuous: the nearest-cell method would create small terraces, and each terrace edge would become a false step (see Terrain from an elevation model for the extent and the reasons).

hydro.fill is Priority-Flood. It has one setting, drain_flats. Left false, a filled depression is perfectly flat, and the directions on it are worked out later by the flat-drainage method of Barnes et al. (2014b), which every direction and accumulation operator applies by default (their flats setting). Set to true, filling raises each cell by the smallest step a number can take above the cell it drains to (the Priority-Flood+ε of Barnes et al., 2014a), so the filled DEM itself slopes, imperceptibly, towards each depression's outlet. Use it when another program, one without a flat-drainage method, will read the filled DEM. Both give every cell a way out: on this map, directions computed from the plain filled DEM with flats set to false leave 1,457 cells without a direction; with flats left on, or with drain_flats on, none.

Filling here raised 1,456 cells (1.5 % of the map), by up to 21.5 m: 224 separate depressions, 66 of them a single cell, the largest 241 cells. On the raw DEM, 251 cells away from the edge had no lower neighbour.

Flow directions and accumulation#

format = "terraflux/1"

[[nodes]]
id = "dem"
op = "raster.read"
params = { path = "elevation.tif" }

[[nodes]]
id = "elevation"
op = "raster.reproject"
inputs = { raster = "dem.raster" }
params = { crs = "EPSG:32610", resolution = 30, extent = [500010, 4760850, 508110, 4771920], resampling = "bilinear" }

[[nodes]]
id = "filled"
op = "hydro.fill"
inputs = { dem = "elevation.raster" }

[[nodes]]
id = "directions"
op = "hydro.flow_direction"
inputs = { dem = "filled.map" }
params = { method = "d8" }
note = "ESRI codes: 1 E, 2 SE, 4 S ... 128 NE."

[[nodes]]
id = "accumulation"
op = "hydro.flow_accumulation"
inputs = { dem = "filled.map" }
params = { method = "d8", units = "cells" }
note = "Cells draining through each cell."

Notice that both operators take the DEM, not each other's output: hydro.flow_accumulation computes the directions it needs itself, by its own method. The directions node is there to see and save the directions; you need it only if another program reads them.

hydro.flow_direction has three methods:

• d8 (the default): ESRI's codes, as in the key idea above.
• dinf: the D-infinity direction of Tarboton (1997), an angle in degrees clockwise from north (0 north, 90 east), any angle, not only multiples of 45°; −1 where the cell has no downhill direction. Tarboton fitted planes to the eight triangles formed by a cell and each pair of adjacent neighbours and took the steepest downhill direction among them.
• mfd: the multiple flow directions of InVEST's routing (its library pygeoprocessing), each cell's share towards each of its eight neighbours packed into one whole number of 32 bits, 4 bits per neighbour. It is a format for other programs, not a map to read.

Enlarge

hydro.flow_accumulation has five methods, the core of this tutorial:

Single-direction D8 concentrates water in lines one cell wide, which is right for a channel but wrong on an open hillslope, where water spreads as a sheet; on planar slopes it draws artificial parallel lines. The multiple flow methods (fd8, quinn, mfd) spread water on hillslopes, and D-infinity sits between them. The results of d8, dinf and fd8 are WhiteboxTools' (D8FlowAccumulation, DInfFlowAccumulation, FD8FlowAccumulation) cell for cell, and mfd is pygeoprocessing's, so that numbers can be compared with those programs.

Two more settings:

• units: cells (the default), area (the cells times the area of a cell, in square map units: m² here) or sca, the specific contributing area, the area divided by the width of a cell (m²/m), the area that drains through each metre of contour. The specific contributing area is what the wetness and stream power indices of part 3 use. At the outlet of the study area's largest catchment: 37,974 cells, 34.18 km², and 1,139,220 m²/m.
• weights: a map of what each cell contributes instead of 1 (an input port; missing cells count as 0). With weights the accumulation is the sum of the weights upslope, as the next step shows.

What lies upstream: a weighted accumulation#

Which reaches of the streams drain cleared land? Treat each cell of open land as contributing 1 and each forest cell 0: the weighted accumulation is then the number of open cells upslope, and dividing it by the plain accumulation gives the share of the catchment that is open.

format = "terraflux/1"

[[nodes]]
id = "dem"
op = "raster.read"
params = { path = "elevation.tif" }

[[nodes]]
id = "elevation"
op = "raster.reproject"
inputs = { raster = "dem.raster" }
params = { crs = "EPSG:32610", resolution = 30, extent = [500010, 4760850, 508110, 4771920], resampling = "bilinear" }

[[nodes]]
id = "filled"
op = "hydro.fill"
inputs = { dem = "elevation.raster" }

[[nodes]]
id = "accumulation"
op = "hydro.flow_accumulation"
inputs = { dem = "filled.map" }

[[nodes]]
id = "ndvi_read"
op = "raster.read"
params = { path = "ndvi.tif" }

[[nodes]]
id = "ndvi"
op = "raster.reproject"
inputs = { raster = "ndvi_read.raster" }
params = { crs = "EPSG:32610", resolution = 30, extent = [500010, 4760850, 508110, 4771920], resampling = "bilinear" }

[[nodes]]
id = "open"
op = "raster.calc"
inputs = { n = "ndvi.raster" }
params = { expression = "1 if n <= 0.55 else 0" }
note = "1 open land, 0 forest."

[[nodes]]
id = "open_acc"
op = "hydro.flow_accumulation"
inputs = { dem = "filled.map", weights = "open.map" }
note = "Open cells upslope."

[[nodes]]
id = "open_share"
op = "raster.calc"
inputs = { o = "open_acc.map", a = "accumulation.map" }
params = { expression = "100 * o / a if a >= 500 else None" }
note = "% of the catchment that is open, on the streams."

The weights can be anything that adds up: the area of a land use, a runoff coefficient, the number of houses or farms, an estimate of nitrogen applied per cell. The result is always "how much of it lies upstream".

terraflux check flow.toml
terraflux run flow.toml --out output

Read the results#

The accumulation map. Flow accumulation ranges over four orders of magnitude, from 1 on ridges to 37,974 cells at the outlet, so it is always drawn on a logarithmic scale. Drag the handle to compare it on the raw DEM and on the filled one, with the same colours:

Raw DEM

Filled DEM

What to look for:

• A connected network. After conditioning, every valley should carry a continuous bright line down to the edge of the map. A bright line that stops suddenly means a depression that was not removed, or a hole in the DEM (see the pitfalls).
• Straight parallel lines on hillslopes are an artefact of D8 on even slopes, not features of the land.
• Where water leaves the map. With D8, 47 cells have code 0, all on the edge: the water of the whole map leaves through them. The largest outlet, in the south-east corner, drains 34.2 km², 38 % of the map.

The weighted accumulation. Open land covers 3.8 % of this forested map, and the streams reflect it: half the stream cells have less than 3.1 % of open land upstream, and 90 % less than 7.8 %. A few short reaches stand out, with up to 39 % open land upstream: small streams that drain clear-cuts. In a study of stream invertebrates or fish, this is the variable you would relate to each sampling site, and the map tells you which reaches to sample to cover the gradient.

Enlarge

Settings that change the answer#

The flow method changes the accumulation of every hillslope cell, and little on the large streams. Move the slider: the five methods drawn with the same colours.

Enlarge

PlayFlow method

The catchment of the outlet is the same within 0.3 %, whatever the method: on the large streams all methods agree. The differences are on the hillslopes. The chart shows, for each method, the share of the map's cells with at least a given accumulation:

A close look at one hillslope makes the difference plain:

Which method to choose:

• D8 for anything that needs one path per cell: streams, stream links, watersheds and HAND (part 2) all use D8 directions in Terraflux. It is also the most common method in the literature, which helps comparison.
• FD8, Quinn or MFD for indices of hillslope wetness and erosion (part 3), where the spreading of water on slopes matters. In temperate forest plots, a wetness index computed with D8 explained less than half of the variation in measured soil moisture and plant species composition that the same index computed with FD8 explained (Kopecký et al., 2021); Sørensen et al. (2006) found that the best method varied with the variable and the site. Part 3 returns to this.
• D-infinity as a compromise: it spreads water, but only between two neighbours, so it disperses less than FD8.
• MFD when your results must match InVEST's sediment and water models.

The exponent of FD8 (1.1 by default, Freeman's value) controls the spreading: larger values send more water to the steepest neighbour and, as they grow, approach D8.

Common pitfalls#

• Skipping the conditioning. Directions and accumulation on a raw DEM stop at every depression (the comparison above). Always fill or breach first; every hydro.* operator says "best with depressions filled" in its description of the DEM.
• Holes in the DEM. Cells without a value are treated as places where water leaves the map, as the edge is: water flows into a hole and disappears. Fill small holes first (raster.fill), or mask large water bodies on purpose.
• Catchments cut by the edge. If the DEM ends before the watershed does, the accumulation of every cell below the cut is too small, and breaching may cut channels to the edge: here a depression near the east edge was breached by a 40 m deep cut through the valley side, straight to the edge. Read a DEM larger than the catchments you study, and cut the results afterwards.
• Real depressions. Lakes, sinkholes, karst and prairie potholes really do hold water. Filling them all connects them to streams that do not exist. Mask them or treat them apart.
• Roads and dams. An embankment across a valley is a dam in the DEM. Filling makes a lake behind it; breaching cuts through it, which is usually what the culvert does in reality (Lindsay & Dhun, 2015).
• Surface models. In forests, radar DEMs such as SRTM partly measure the top of the canopy. Clearings become pits, forest edges become walls, and streams under the canopy are hard to see. Prefer a bare-earth DEM (lidar) for small streams.
• Units. cells and area differ by the area of a cell (900 m² here); sca divides by its width as well. Say which you report, and do not compare counts of cells between DEMs of different resolutions.
• Resolution. A DEM of 30 m cannot see a ditch 2 m wide, and one of 1 m sees every wheel rut. Choose the resolution for the flow paths you care about.

Going further#

• Part 2: Streams, watersheds and HAND: turn the accumulation into a stream network, split it into links with their Strahler order, delineate the watershed of a gauge and measure the height of every cell above the nearest stream.
• Part 3: Wetness, stream power and sediment connectivity: the topographic wetness index and its relatives, and the index of connectivity of InVEST's sediment model.
• Part 4: Watershed simulation with SWAT+: from the terrain to daily water and sediment budgets.
• Terrain from an elevation model: slope, aspect, curvature and topographic position.
• The operators' reference: hydro.fill, hydro.breach, hydro.flow_direction, hydro.flow_accumulation.

References#

• Allan, J. D. (2004). Landscapes and riverscapes: The influence of land use on stream ecosystems. Annual Review of Ecology, Evolution, and Systematics, 35, 257–284. https://doi.org/10.1146/annurev.ecolsys.35.120202.110122
• Barnes, R., Lehman, C., & Mulla, D. (2014a). Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Computers & Geosciences, 62, 117–127. https://doi.org/10.1016/j.cageo.2013.04.024
• Barnes, R., Lehman, C., & Mulla, D. (2014b). An efficient assignment of drainage direction over flat surfaces in raster digital elevation models. Computers & Geosciences, 62, 128–135. https://doi.org/10.1016/j.cageo.2013.01.009
• Freeman, T. G. (1991). Calculating catchment area with divergent flow based on a regular grid. Computers & Geosciences, 17(3), 413–422. https://doi.org/10.1016/0098-3004(91)90048-I
• Kopecký, M., Macek, M., & Wild, J. (2021). Topographic Wetness Index calculation guidelines based on measured soil moisture and plant species composition. Science of the Total Environment, 757, 143785. https://doi.org/10.1016/j.scitotenv.2020.143785
• Lindsay, J. B. (2016). Efficient hybrid breaching-filling sink removal methods for flow path enforcement in digital elevation models. Hydrological Processes, 30(6), 846–857. https://doi.org/10.1002/hyp.10648
• Lindsay, J. B., & Dhun, K. (2015). Modelling surface drainage patterns in altered landscapes using LiDAR. International Journal of Geographical Information Science, 29(3), 397–411. https://doi.org/10.1080/13658816.2014.975715
• Moore, I. D., Grayson, R. B., & Ladson, A. R. (1991). Digital terrain modelling: A review of hydrological, geomorphological, and biological applications. Hydrological Processes, 5(1), 3–30. https://doi.org/10.1002/hyp.3360050103
• O'Callaghan, J. F., & Mark, D. M. (1984). The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing, 28(3), 323–344. https://doi.org/10.1016/S0734-189X(84)80011-0
• Quinn, P., Beven, K., Chevallier, P., & Planchon, O. (1991). The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models. Hydrological Processes, 5(1), 59–79. https://doi.org/10.1002/hyp.3360050106
• Sørensen, R., Zinko, U., & Seibert, J. (2006). On the calculation of the topographic wetness index: Evaluation of different methods based on field observations. Hydrology and Earth System Sciences, 10(1), 101–112. https://doi.org/10.5194/hess-10-101-2006
• Tarboton, D. G. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research, 33(2), 309–319. https://doi.org/10.1029/96WR03137