⚙️ How-To Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 3 of 20 Updated

How to Condition a DEM for Hydrology: Breaching or Filling

Problem statement

Before a DEM can route water, its depressions — cells or groups of cells with no lower neighbour — have to go. Filling raises each depression to the level at which it would overflow. Breaching carves a path downward through the barrier that closes it. Both make every cell drain, and both change the terrain: filling builds flat lakes, breaching digs trenches. Which one, and with what limits, decides how much of the DEM you rewrite and whether a road embankment or a reservoir dam gets cut through.

Measured with WhiteboxTools on the USGS 3DEP 10 m DEM of the Esopus Creek catchment above Coldbrook, New York (published drainage area 497.3 km²):

  • Inside the basin, 21,745 cells (0.44%) sat in 2,855 depressions; 1,417 were single cells and 2,666 were shallower than 0.5 m.
  • Filling changed 21,775 cells, raising them by up to 10.14 m. Least-cost breaching with a 20-cell search changed 10,670 cells, raising 7,360 and lowering 3,310.
  • Unconstrained breaching cut one cell 190.42 m — from 183.74 m to −6.68 m — and lowered 5,792 cells by more than 5 m.
  • Every sensible option gave a watershed within 0.4 km² of 492.4 km²; breaching limited to 2 m deep and 50 cells long raised a third of the basin's cells by up to 211 m and cut the watershed to 448.2 km².

Quick answer

For most catchment work, breach short and shallow barriers and fill what remains, with a cap on cut depth:

import whitebox

wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.breach_depressions_least_cost("dem.tif", "conditioned.tif", dist=100, max_cost=50, fill=True)

dist is the longest breach channel in cells, max_cost caps the total elevation cut along it, and fill=True fills whatever cannot be breached within those limits. Then compare the conditioned DEM with the original: count changed cells, the deepest cut and the highest fill, and look at where they are.

Two panels comparing filling, which raises a depression to its spill level, with breaching, which cuts a channel through the barrier, with the cells each changed on the Esopus DEM.
Filling edits the depression; breaching edits the barrier. Both make the valley drain.

Step-by-step solution

1. Measure the depressions before removing them

depth_in_sink reports how deep each cell sits below its depression's spill level. On the basin DEM, 0.96% of cells were in depressions, with a median depth of 0.148 m and a maximum of 10.13 m. Most were noise: 2,613 of the 2,855 depressions inside the NLDI basin covered 10 cells or fewer. The few large ones are the ones to look at — here a flat area at 256.03 m covering 3,506 cells, and depressions behind road crossings.

2. Fill: simplest and flattest

fill_depressions raises every depression cell to its spill level; fix_flats=True then adds a tiny gradient so the flat still drains. It took 0.40 s and changed 21,775 cells inside the basin by up to 10.14 m, adding 1.149 million m³ of "terrain". The watershed at the gauge was 492.42 km², with an IoU of 0.985 against the USGS polygon. Filling never lowers terrain, so it cannot cut through a real dam — but it turns every depression, including those behind roads, into a flat lake.

3. Least-cost breaching: shortest, cheapest channels

breach_depressions_least_cost searches outward from each depression for the path to lower ground with the least total cut. With dist=20 it changed 10,670 cells and the deepest cut was 16.34 m; in an earlier run with dist=100, 8,217 cells and 70.00 m, because a longer search finds routes through higher barriers. Both left watersheds within 0.4 km² of the filled result. On the full 12.5-million-cell DEM a 100-cell search took 330 s, against 0.56 s for filling.

4. Cap the cost of a breach

max_cost stops breaches whose total cut exceeds a limit, and fill=True fills those depressions instead. With max_cost=50, the deepest cut inside the basin fell to 2.83 m, no cell was lowered by more than 5 m, and the watershed stayed at 492.41 km². This is the setting that keeps road embankments and dams from being carved into trenches while still removing the small artefacts by breaching.

5. Beware unconstrained breaching

breach_depressions with no limits lowered 45,919 cells in the basin, 5,792 of them by more than 5 m, and removed 34.9 million m³. Its deepest cut, 190.42 m, drove a trench below sea level from a flat valley floor near the Ashokan Reservoir. The watershed area was still 492.44 km²: a good area does not mean a sensible DEM.

6. Beware tight limits that force filling

Limiting breaches to 2 m deep and 50 cells long left most depressions unbreachable, and the fallback filled them. 1,648,677 cells inside the basin — a third of them — were raised by up to 211 m, flooding whole valleys to the level of their outlets, and the watershed shrank to 448.21 km², 9.9% short, because flow crossed the resulting plateaus in the wrong places.

7. Inspect where the DEM changed

Map the difference between conditioned and original DEMs, and list the largest changes with their coordinates (Example 3). Deep cuts should line up with real channels under bridges or culverts; large fills with real ponds or reservoirs. A 72 m cut through a hillside or a 200 m fill across a valley is a conditioning artefact that will reappear in every downstream product.

8. Validate with the watershed, then with streams

Compare the watershed area and boundary with a reference, and the derived stream network with a mapped one. Filled and least-cost-breached watersheds differed by 3,955 cells, 0.40 km²; their stream cells above 1 km² overlapped for 44,983 of about 46,400 cells.

Table of conditioning options on the Esopus DEM with cells changed, deepest cut, highest fill and watershed area for fill, least-cost breaching at two distances, capped breaching, unconstrained breaching and tightly limited breaching.
Four settings gave almost the same watershed area; the DEM edits behind them differed by two orders of magnitude.

Code examples

Example 1 — condition the DEM several ways 

import os
import time

import whitebox

os.makedirs("work", exist_ok=True)
wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("work"))
DEM = os.path.abspath("esopus_3dep13_utm_basin.tif")

OPTIONS = {
    "fill": lambda out: wbt.fill_depressions(DEM, out, fix_flats=True),
    "breach_lc20": lambda out: wbt.breach_depressions_least_cost(DEM, out, dist=20, fill=True),
    "breach_lc100_cap50": lambda out: wbt.breach_depressions_least_cost(DEM, out, dist=100, max_cost=50, fill=True),
    "breach_unlimited": lambda out: wbt.breach_depressions(DEM, out, fill_pits=True),
    "breach_2m_50cells": lambda out: wbt.breach_depressions(DEM, out, max_depth=2.0, max_length=50, fill_pits=True),
}
for name, run in OPTIONS.items():
    start = time.perf_counter()
    run(f"{name}.tif")
    print(f"{name:20} {time.perf_counter() - start:6.2f} s")
fill                   0.40 s
breach_lc20            2.94 s
breach_lc100_cap50    14.78 s
breach_unlimited       1.91 s
breach_2m_50cells      2.05 s

Example 2 — what each option changed, inside the basin 

import geopandas as gpd
import numpy as np
import rasterio
from rasterio import features

with rasterio.open(DEM) as src:
    original = src.read(1).astype("float64")
    valid = original != src.nodata
    transform, crs, cell = src.transform, src.crs, src.res[0]
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(crs)
inside = features.rasterize(basin.geometry, out_shape=original.shape, transform=transform).astype(bool) & valid


def changes(name):
    with rasterio.open(f"work/{name}.tif") as src:
        dz = np.where(inside, src.read(1).astype("float64") - original, 0.0)
    raised, lowered = dz > 1e-3, dz < -1e-3
    print(f"{name:20} changed {int((raised | lowered).sum()):>9,}  raised {int(raised.sum()):>9,} (max {dz.max():+7.2f} m)  "
          f"lowered {int(lowered.sum()):>7,} (max {dz.min():+8.2f} m, >5 m: {int((dz < -5).sum()):,})  "
          f"net {dz.sum() * cell * cell / 1e6:+9.3f} Mm3")


for name in OPTIONS:
    changes(name)
fill                 changed    21,775  raised    21,775 (max  +10.14 m)  lowered       0 (max    +0.00 m, >5 m: 0)  net    +1.149 Mm3
breach_lc20          changed    10,670  raised     7,360 (max   +0.92 m)  lowered   3,310 (max   -16.34 m, >5 m: 27)  net    -0.050 Mm3
breach_lc100_cap50   changed    12,964  raised     9,702 (max  +10.14 m)  lowered   3,262 (max    -2.83 m, >5 m: 0)  net    +0.671 Mm3
breach_unlimited     changed    46,098  raised       179 (max   +0.59 m)  lowered  45,919 (max  -190.42 m, >5 m: 5,792)  net   -34.887 Mm3
breach_2m_50cells    changed 1,652,204  raised 1,648,677 (max +211.24 m)  lowered   3,527 (max    -2.00 m, >5 m: 0)  net +11951.359 Mm3

Example 3 — the watershed each option gives, and where the deepest cut is 

from pyproj import Transformer
from shapely.geometry import Point

x, y = Transformer.from_crs("EPSG:4269", crs, always_xy=True).transform(-74.2701944, 42.0144722)
gpd.GeoDataFrame(geometry=[Point(x, y)], crs=crs).to_file("work/gauge.shp")

for name in OPTIONS:
    wbt.d8_pointer(f"{name}.tif", f"{name}_d8.tif")
    wbt.d8_flow_accumulation(f"{name}_d8.tif", f"{name}_acc.tif", out_type="cells", pntr=True)
    wbt.snap_pour_points("gauge.shp", f"{name}_acc.tif", f"{name}_outlet.shp", snap_dist=150)
    wbt.watershed(f"{name}_d8.tif", f"{name}_outlet.shp", f"{name}_ws.tif")
    with rasterio.open(f"work/{name}_ws.tif") as src:
        area = (src.read(1) == 1).sum() * cell * cell / 1e6
    with rasterio.open(f"work/{name}.tif") as src:
        dz = np.where(inside, src.read(1).astype("float64") - original, 0.0)
    row, col = np.unravel_index(dz.argmin(), dz.shape)
    lon, lat = Transformer.from_crs(crs, "EPSG:4326", always_xy=True).transform(*(transform * (col + 0.5, row + 0.5)))
    print(f"{name:20} watershed {area:7.2f} km2 | deepest cut {dz.min():8.2f} m at {lat:.5f}, {lon:.5f} "
          f"({original[row, col]:.2f} m -> {original[row, col] + dz.min():.2f} m)")
fill                 watershed  492.42 km2 | deepest cut     0.00 m at 42.20711, -74.53061 (690.50 m -> 690.50 m)
breach_lc20          watershed  492.03 km2 | deepest cut   -16.34 m at 42.10438, -74.22223 (578.99 m -> 562.65 m)
breach_lc100_cap50   watershed  492.41 km2 | deepest cut    -2.83 m at 42.14960, -74.26344 (546.12 m -> 543.29 m)
breach_unlimited     watershed  492.44 km2 | deepest cut  -190.42 m at 42.00119, -74.26660 (183.74 m -> -6.68 m)
breach_2m_50cells    watershed  448.21 km2 | deepest cut    -2.00 m at 42.09874, -74.21335 (532.19 m -> 530.20 m)

Filling lowers nothing, so its "deepest cut" line just reports the first cell of an all-zero difference.

Explanation

Why DEMs are full of depressions

Elevation models are measured and interpolated surfaces. Bridges and dense vegetation register as barriers across channels; interpolation leaves small pits; rounding creates flat patches; real features such as ponds, quarries and reservoirs are genuine depressions. In a 10 m lidar-derived DEM most depressions are a few cells and a few centimetres, and a handful are large.

Why filling and breaching give the same watershed

A watershed boundary follows ridges, which neither method touches much. Both methods only change where flow goes inside the depression and across its barrier, and the flow still ends up leaving by the same outlet. That is why four of the settings gave areas within 0.41 km² of each other, and why area alone cannot tell you whether conditioning was sensible.

Why unconstrained breaching digs so deep

Without a cost limit, the algorithm breaches every depression by lowering a path until it reaches a lower cell, however deep that requires. A large flat valley floor adjoining a reservoir whose surface is lower, or a DEM edge with nodata, can pull that path down by hundreds of metres. The trench is invisible in a hillshade at basin scale and obvious in a difference map.

Why tight limits make things worse, not safer

max_depth and max_length stop breaching, not conditioning: the depressions they refuse to breach are then filled. When the limits are small compared with real barriers, the filling floods valleys to their spill levels, producing large flats where flow direction is arbitrary. Capping the total cost and filling the remainder, as max_cost does, avoids that by breaching the many cheap barriers and filling only the few expensive ones.

Bar chart of the deepest cut inside the basin for least-cost breaching at 20 and 100 cells, capped breaching and unconstrained breaching.
Search distance and cost caps decide whether breaching nudges a culvert or digs a canyon.

Edge cases or notes

  • Real depressions — karst, kettle lakes, closed basins — should be masked, not removed.
  • Reservoirs appear as large flats; flow direction across them is arbitrary after flat resolution.
  • Breach runtime grows with dist; on the full grid a 100-cell least-cost search took 330 s.
  • NoData edges act as sinks; make sure nodata is tagged before conditioning.
  • Stream burning is an alternative for known channels; see burning a river network into a DEM.
  • pysheds' own sequence of fill_pits, fill_depressions and resolve_flats left pits on this DEM and a watershed 74% too small; see fixing a watershed that comes out as a few pixels.
  • Keep the conditioned DEM for routing only. Use the original elevations for slope, HAND and anything physical.

FAQ

Should I breach or fill a DEM?

Breach small, cheap barriers and fill the rest. Least-cost breaching with a cost cap changed far fewer cells than filling on the Esopus DEM and kept the deepest cut to 2.83 m.

Does breaching or filling change the watershed area?

Very little when done sensibly: fill, least-cost breaching and capped breaching all gave 492.0–492.4 km². Very tight breach limits that forced large fills cut it to 448.2 km².

Why did breaching lower my DEM by hundreds of metres?

Unconstrained breaching lowers a path until it reaches lower ground. On the Esopus DEM it cut one cell by 190.42 m; set max_cost, or limit dist, and fill the rest.

How many cells does conditioning change?

Usually under 1%. Inside the Esopus basin filling changed 21,775 cells (0.44%) and least-cost breaching with a 20-cell search 10,670.

What does fix_flats do in WhiteboxTools fill_depressions?

It adds a very small gradient across filled flat areas so that every cell has a downslope direction afterwards. Without it, filled depressions are perfectly flat and flow direction is undefined there.

How long does least-cost breaching take?

It depends on the search distance. On the full 12.5-million-cell DEM a 100-cell search took 330 s against 0.56 s for filling; on the basin crop, with a cost cap, 14.8 s.