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,526 cells (0.43%) sat in 2,829 depressions; 1,408 were single cells and 2,645 were shallower than 0.5 m.
• Filling changed 21,587 cells, raising them by up to 10.01 m. Least-cost breaching with a 20-cell search changed 10,645 cells, raising 7,311 and lowering 3,334.
• Unconstrained breaching cut one cell 190.51 m — from 183.65 m to −6.86 m — and lowered 5,786 cells by more than 5 m.
• Every sensible option gave a watershed between 491.9 and 492.3 km²; breaching limited to 2 m deep and 50 cells long raised a third of the basin's cells by up to 212 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_max_procs(1)                      # one thread: parallel runs give slightly different results
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.

fill                   0.44 s
breach_lc20            3.25 s
breach_lc100_cap50    15.71 s
breach_unlimited       1.62 s
breach_2m_50cells      2.11 s

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, 1.00% of cells were in depressions, with a median depth of 1.27 m and a maximum of 20.47 m. Most were small: 2,586 of the 2,829 depressions inside the NLDI basin covered 10 cells or fewer. The few large ones are the ones to look at — here a depression bottoming out at 256.03 m and covering 4,150 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.44 s and changed 21,587 cells inside the basin by up to 10.01 m, adding 1.077 million m³ of "terrain". The watershed at the gauge was 492.21 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,645 cells and the deepest cut was 16.32 m; with dist=100, 8,087 cells and 23.50 m, because a longer search finds routes through higher barriers. On the full 12.2-million-cell DEM a 100-cell search took 377 s on one thread, against 0.74 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 3.32 m, no cell was lowered by more than 5 m, and the watershed stayed at 492.33 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,933 cells in the basin, 5,786 of them by more than 5 m, and removed 35.2 million m³. Its deepest cut, 190.51 m, drove a trench below sea level from a flat valley floor near the Ashokan Reservoir. The watershed area was still 492.22 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,649,388 cells inside the basin — a third of them — were raised by up to 212 m, flooding whole valleys to the level of their outlets, and the watershed shrank to 448.15 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 20 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 10,017 cells, 1.00 km², though their areas were only 0.28 km² apart; their stream cells above 1 km² overlapped for 27,652 of about 28,500 cells.

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_max_procs(1)                      # one thread: parallel runs give slightly different results
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.44 s
breach_lc20            3.25 s
breach_lc100_cap50    15.71 s
breach_unlimited       1.62 s
breach_2m_50cells      2.11 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,587  raised    21,587 (max  +10.01 m)  lowered       0 (max    +0.00 m, >5 m: 0)  net    +1.077 Mm3
breach_lc20          changed    10,645  raised     7,311 (max   +1.79 m)  lowered   3,334 (max   -16.32 m, >5 m: 25)  net    -0.058 Mm3
breach_lc100_cap50   changed    13,035  raised     9,733 (max  +10.01 m)  lowered   3,302 (max    -3.32 m, >5 m: 0)  net    +0.661 Mm3
breach_unlimited     changed    46,085  raised       152 (max   +1.64 m)  lowered  45,933 (max  -190.51 m, >5 m: 5,786)  net   -35.199 Mm3
breach_2m_50cells    changed 1,652,937  raised 1,649,388 (max +211.87 m)  lowered   3,549 (max    -2.00 m, >5 m: 0)  net +11978.253 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.21 km2 | deepest cut     0.00 m at 42.20717, -74.53064 (690.40 m -> 690.40 m)
breach_lc20          watershed  491.93 km2 | deepest cut   -16.32 m at 42.10444, -74.22214 (578.92 m -> 562.60 m)
breach_lc100_cap50   watershed  492.33 km2 | deepest cut    -3.32 m at 42.10198, -74.45160 (376.61 m -> 373.28 m)
breach_unlimited     watershed  492.22 km2 | deepest cut  -190.51 m at 42.00116, -74.26675 (183.65 m -> -6.86 m)
breach_2m_50cells    watershed  448.15 km2 | deepest cut    -2.00 m at 42.14518, -74.21305 (579.21 m -> 577.21 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.40 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.

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.
• Run WhiteboxTools on one thread (wbt.set_max_procs(1)): multithreaded runs here differed in 10,471 cells.
• Breach runtime grows with dist; on the full grid a 100-cell least-cost search took 377 s on one thread.
• 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, run once, left a watershed of 459.0 km² on this DEM, 7.7% short; a second fill and resolve brought it to 491.0 km². 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.

Internal links

• DEM hydrology explained: from elevation to where water goes — where conditioning fits
• Flow direction explained: D8, D-infinity and MFD compared — the step after conditioning
• How to fill sinks and derive flow direction from a DEM — the basic fill workflow
• Fixing streams that stop at roads, bridges and embankments — barriers in the DEM
• How to burn a known river network into a DEM — forcing known channels
• Fixing streams that run in straight parallel lines across flat ground — what large fills cause
• How to delineate a watershed from a pour point in Python — using the conditioned DEM
• Digital elevation models explained: DEM, DSM and DTM — why bare-earth models need less conditioning

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 3.32 m.

Does breaching or filling change the watershed area?

Very little when done sensibly: fill, least-cost breaching and capped breaching all gave 491.9–492.3 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.51 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,587 cells (0.44%) and least-cost breaching with a 20-cell search 10,645.

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.2-million-cell DEM a 100-cell search took 377 s on one thread against 0.74 s for filling; on the basin crop, with a cost cap, 15.7 s.