From DEM to TWI
Why are some hillslope cells wetter than others? Walk through every step — from raw elevation data to a wetness map you can interact with on a real catchment.
The six-step pipeline — click to jump
1 · What makes a hillslope cell wet?
Water flows from high cells to low cells under gravity. A cell tends to stay wet when:
- Many upslope cells drain into it — the contributing area
ais large. - The local slope is gentle — water lingers rather than rushing away, so
tan(β)is small.
Tiny drainage area + steep slope = dry.
TWI turns that sentence into a number [1].
2 · A real catchment: USGS gauge 01031500
The Piscataquis River drainage in central Maine — a CAMELS benchmark catchment. Toggle TWI (qualitative) in the layers control to see a simulated wetness pattern overlaid on the watershed. Click the gauge marker for metadata.
delineate_watershed("01031500") →
fetch_dem(3DEP, 30m) → compute_twi(method="d8") →
returns a georeferenced raster. The boundary polygon is from that same pipeline.
3 · Flow direction — click any cell
A 5×5 synthetic DEM. Each cell routes water to its lowest 8-neighbour (D8 algorithm). Click any cell to trace its flow path to the outlet, highlight its upstream contributing cells, and compute an approximate TWI.
■ upstream contributors ■ flow path to outlet ■ outlet cell (66 m, row 4 col 2)
Under the hood: three things worth understanding
Why ln() — not just a / tan(β)?
The raw ratio a / tan(β) spans roughly 4 orders of magnitude across a real catchment —
a valley floor might have a/tan(β) ≈ 40,000 while a steep ridge has ≈ 5.
Plotting those values linearly produces a map where virtually every cell looks the same shade and only
the valley floor is visible.
Taking the natural log compresses the range to roughly 2–14 (a 12-unit spread), where gradients are perceptible everywhere on the hillslope. It also has a physical basis: Beven & Kirkby (1979) derived TWI directly from the steady-state soil-moisture deficit equation — the ln falls out of the algebra, it is not a cosmetic choice.
A practical bonus: at very low slopes (tan(β) → 0, common in flat valley floors),
a / tan(β) diverges to infinity; ln keeps the index finite and
interpretable with a small epsilon floor.
D8 vs D∞ routing algorithms
D8 (O'Callaghan & Mark, 1984) routes 100% of each cell's flow to the single steepest of its 8 cardinal/diagonal neighbours. It is fast and deterministic, but forces flow into only 8 discrete directions — producing "stair-step" channel artefacts and overestimating channelisation on smooth hillslopes.
D∞ (Tarboton, 1997) computes the steepest downslope direction as a continuous angle, then splits flow proportionally between the two neighbours bounding that angle. This allows any of 2π possible directions, producing smoother contributing-area fields and more realistic distributed wetness patterns — but at higher computational cost.
For channel network extraction D8 is usually preferred (discrete routing matches the physics);
for spatially continuous TWI fields, D∞ reduces artefacts. The AI-Hydro
compute_twi tool defaults to method="d8" but accepts
method="dinf".
Resolution dependence — why you must report DEM resolution
Specific contributing area is a = A_upslope / w where w is the contour
width (≈ cell size in raster DEMs). At coarser resolution each cell accumulates the same physical
drainage area but w is larger — so a is proportionally larger, and
TWI inflates.
Sørensen et al. (2006) [3] quantified this: mean TWI varied by 1–2 units between 25 m and 200 m grids of the same catchments — a shift large enough to change wetness-class assignments. This means TWI values from different studies are not directly comparable unless computed at the same resolution.
Practical guidance: 30 m (SRTM, 3DEP-30) is the standard for continental-scale CAMELS comparisons; 10 m (3DEP-10) for sub-catchment studies. Always cite the source DEM and its resolution when reporting TWI statistics.
4 · Live TWI calculator
Adjust the inputs — TWI and the distribution histogram update instantly [1,2].
The equation
a = specific contributing area (m²/m) ·
β = local slope (radians)
Synthetic catchment distribution
16×16 synthetic hillslope (same as the Python heatmap below). Marker shows your current value.
| TWI | Class |
|---|---|
| < 5 | Dry — ridge / steep hillslope |
| 5 – 8 | Moderate — mid-slope |
| 8 – 11 | Wet — lower hillslope / hollow |
| > 11 | Very wet — riparian / saturated |
5 · Python cell — point TWI calculation
Source is regenerated from JavaScript each time you move a slider. Click Run to execute in the AI-Hydro kernel.
6 · Synthetic catchment TWI heatmap
A 16×16 synthetic hillslope — contributing area peaks in the central valley, slope tapers downhill. Run to see where wetness concentrates.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams.update({
'figure.facecolor':'#0a0a15', 'axes.facecolor':'#0f0f1e',
'text.color':'#7dd3fc', 'axes.labelcolor':'#7dd3fc',
'xtick.color':'#94a3b8', 'ytick.color':'#94a3b8',
'axes.edgecolor':'#334155',
})
n = 16
y_idx, x_idx = np.mgrid[0:n, 0:n]
area_grid = 50 + 200 * np.exp(-((x_idx - n/2)**2) / 18) * (1 + y_idx / n)
slope_grid = 18 - 12 * (y_idx / n) + 4 * np.abs(x_idx - n/2) / n
twi_grid = np.log(area_grid / np.maximum(np.tan(np.deg2rad(slope_grid)), 1e-6))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3.6))
im = ax1.imshow(twi_grid, cmap='GnBu', origin='lower')
ax1.set_title('TWI grid', fontsize=11)
ax1.set_xlabel('x'); ax1.set_ylabel('y (↑ downslope)')
cb = fig.colorbar(im, ax=ax1, fraction=0.046); cb.set_label('TWI')
ax2.hist(twi_grid.ravel(), bins=20, color='#00A3FF', edgecolor='#0a0a15', alpha=0.85)
ax2.axvline(twi_grid.mean(), color='#00FFFF', linewidth=1.5, linestyle='--')
ax2.text(twi_grid.mean()+0.1, ax2.get_ylim()[1]*0.9,
f'mean={twi_grid.mean():.1f}', color='#00FFFF', fontsize=9)
ax2.set_title('TWI distribution', fontsize=11)
ax2.set_xlabel('TWI'); ax2.set_ylabel('count')
plt.tight_layout()
plt.show()
print(f"Grid mean TWI = {twi_grid.mean():.2f}, std = {twi_grid.std():.2f}")
print(f"Dry (<5): {(twi_grid<5).sum()} | Wet (>8): {(twi_grid>8).sum()}")
7 · Checkpoint
8 · What the agent would do
The agent runs the full pipeline — DEM fetch, fill-sinks, D8 routing, accumulation, slope, TWI — and returns a georeferenced raster. The simulated overlay you saw above is a qualitative preview; the agent produces the real thing.
References
- Beven, K.J. & Kirkby, M.J. (1979). A physically based, variable contributing area model of basin hydrology. Hydrological Sciences Bulletin, 24(1), 43–69. DOI: 10.1080/02626667909491834
- Quinn, P.F., 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. DOI: 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, 101–112. DOI: 10.5194/hess-10-101-2006
- Addor, N., Newman, A.J., Mizukami, N. & Clark, M.P. (2017). The CAMELS data set: catchment attributes and meteorology for large-sample studies. Hydrology and Earth System Sciences, 21, 5293–5313. DOI: 10.5194/hess-21-5293-2017
- Bartos, M. (2020). pysheds: simple and fast watershed delineation in Python. github.com/mdbartos/pysheds
- 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. DOI: 10.1029/96WR03137
compute_twi in the agent for the real raster.