Where does the silt come from? Soil-erosion risk in a reservoir's catchment¶

Take-home notebook · Workshop 1 · IC-TWEP 2026 · QuestFeed¶

Reservoirs in North Africa and the Middle East lose storage every year to sediment washed off their catchments. Knowing where soil erosion is most intense tells you where soil and water conservation will protect the reservoir best. Here we estimate erosion risk with the Revised Universal Soil Loss Equation (RUSLE):

$$A = R \times K \times LS \times C \times P$$

Factor Meaning Data used here
R rainfall erosivity the catchment's own CHIRPS rainfall record (1981 → today)
K soil erodibility SoilGrids 250 m: sand, silt, clay and organic carbon
LS slope length and steepness Copernicus DEM (90 m) with flow accumulation
C land cover and management ESA WorldCover 2021 (10 m)
P support practices (terraces, contours) set to 1 — no conservation assumed

The result, A, is an estimate of long-term average soil loss in tonnes per hectare per year. RUSLE is a screening tool: use it to rank areas, not as a precise measurement.

Ready.
In [2]:
# @title 1 · Choose a reservoir and set up a 90 m grid over its catchment
RESERVOIR = "Beni Haroun — Algeria"  # @param ["Beni Haroun — Algeria", "Hammam Boughrara (Tlemcen) — Algeria", "Beni Bahdel (Tlemcen) — Algeria", "Al Massira — Morocco", "King Fahd Dam — Saudi Arabia", "Wadi Dayqah — Oman"]
sites = {s["label"]: s for s in fetch_json("sites.json")}
site = sites[RESERVOIR]["id"]
info = fetch_json(f"{site}/summary.json")
catchment = fetch_vector(f"{site}/catchment.geojson")
lon0, lat0 = catchment.unary_union.centroid.coords[0]
crs = f"EPSG:{(32600 if lat0 >= 0 else 32700) + int((lon0 + 180) // 6) + 1}"
c_utm = catchment.to_crs(crs)
x0, y0, x1, y1 = c_utm.total_bounds
RES = 90
gbox = GeoBox.from_bbox((x0, y0, x1, y1), crs=crs, resolution=RES)
inside = features.rasterize(c_utm.geometry, out_shape=gbox.shape, transform=gbox.transform, fill=0, default_value=1).astype(bool)
print(f"{info['name']}: catchment {info['catchment_km2']:,.0f} km², grid {gbox.shape[1]} × {gbox.shape[0]} at {RES} m")

Beni Haroun Reservoir: catchment 7,574 km², grid 1584 × 954 at 90 m

R — rainfall erosivity¶

Erosivity grows faster than rainfall itself, because intense storms detach far more soil. We use the relationship of Renard & Freimund (1994) for mean annual precipitation P (mm): R = 0.0483 · P¹·⁶¹ (MJ mm ha⁻¹ h⁻¹ yr⁻¹), valid for P up to about 850 mm, with P from the catchment's CHIRPS record.

In [3]:
# @title 2 · Rainfall erosivity
rain = fetch_csv(f"{site}/rainfall_monthly.csv")
rain["year"] = rain["month"].str[:4].astype(int)
P = rain[(rain.year >= 1991) & (rain.year <= 2020)].groupby("year")["precip_mm"].sum().mean()
R = 0.0483 * P ** 1.61
print(f"Mean annual rainfall 1991–2020: {P:,.0f} mm → R = {R:,.0f} MJ mm ha⁻¹ h⁻¹ yr⁻¹")

Mean annual rainfall 1991–2020: 569 mm → R = 1,319 MJ mm ha⁻¹ h⁻¹ yr⁻¹

K — soil erodibility¶

Fine sands and silts erode most easily; clays and soils rich in organic matter resist. We compute K from SoilGrids topsoil properties with the EPIC equation (Williams, 1995), converted to SI units (t ha h ha⁻¹ MJ⁻¹ mm⁻¹).

In [4]:
# @title 3 · Soil erodibility from SoilGrids (about a minute)
def soilgrids(prop):
    url = f"https://files.isric.org/soilgrids/latest/data/{prop}/{prop}_0-5cm_mean.vrt"
    out = np.full(gbox.shape, np.nan, np.float32)
    with rasterio.open(url) as src:
        reproject(rasterio.band(src, 1), out, dst_transform=gbox.transform, dst_crs=crs, resampling=Resampling.average,
                  src_nodata=src.nodata, dst_nodata=np.nan)
    return out
sand, silt, clay, soc = (soilgrids(p) for p in ("sand", "silt", "clay", "soc"))
SAN, SIL, CLA = sand / 10, silt / 10, clay / 10          # g/kg → %
C_org = soc / 100                                          # dg/kg → %
f_csand = 0.2 + 0.3 * np.exp(-0.256 * SAN * (1 - SIL / 100))
f_clsi = (SIL / np.maximum(CLA + SIL, 1e-3)) ** 0.3
f_orgc = 1 - 0.25 * C_org / (C_org + np.exp(3.72 - 2.95 * C_org))
SN1 = 1 - SAN / 100
f_hisand = 1 - 0.7 * SN1 / (SN1 + np.exp(-5.51 + 22.9 * SN1))
K = 0.1317 * f_csand * f_clsi * f_orgc * f_hisand
print(f"K over the catchment: median {np.nanmedian(K[inside]):.3f}, range {np.nanpercentile(K[inside], 5):.3f}–{np.nanpercentile(K[inside], 95):.3f}")

K over the catchment: median 0.017, range 0.016–0.017

LS — slope length and steepness¶

Water gathers speed and volume as it flows downhill. LS combines the upslope contributing area (from flow accumulation on the DEM) with slope steepness (Moore & Burch, 1986): LS = (Aₛ / 22.13)^0.4 · (sin β / 0.0896)^1.3, with the slope length capped at 300 m, a common practice that keeps channels from dominating. Flow direction is the steepest of the eight neighbours (D8).

In [5]:
# @title 4 · Slope length and steepness (about a minute)
dem_items = list(Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["cop-dem-glo-90"], bbox=list(catchment.total_bounds)).items())
dem = load(dem_items, bands=["data"], geobox=gbox, resampling="bilinear")["data"].isel(time=0).values.astype("float64")
dem[(dem < -100) | ~np.isfinite(dem)] = np.nan
z = np.nan_to_num(dem, nan=np.nanmax(dem))
H, W = z.shape
# D8 steepest-descent receivers
offs = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]
dist = np.array([np.sqrt(2), 1, np.sqrt(2), 1, 1, np.sqrt(2), 1, np.sqrt(2)]) * RES
zp = np.pad(z, 1, mode="edge")
drops = np.stack([(z - zp[1 + dy:1 + dy + H, 1 + dx:1 + dx + W]) / d for (dy, dx), d in zip(offs, dist)])
k = drops.argmax(0)
rr = np.indices((H, W))[0] + np.array([o[0] for o in offs])[k]
cc = np.indices((H, W))[1] + np.array([o[1] for o in offs])[k]
ok = (drops.max(0) > 0) & (rr >= 0) & (rr < H) & (cc >= 0) & (cc < W)
recv = np.where(ok, rr * W + cc, -1).ravel()
# Flow accumulation by repeated donor summation; six passes reach the 300 m slope-length cap at 90 m cells.
acc = np.ones(H * W)
has = recv >= 0
for _ in range(6):
    acc = 1 + np.bincount(recv[has], weights=acc[has], minlength=H * W)
acc = acc.reshape(H, W)
gy, gx = np.gradient(z, RES)
beta = np.arctan(np.hypot(gx, gy))
As = np.minimum(acc * RES, 300.0)                     # specific contributing area, slope length capped at 300 m
LS = (As / 22.13) ** 0.4 * (np.sin(beta) / 0.0896) ** 1.3
LS[~inside] = np.nan
print(f"LS over the catchment: median {np.nanmedian(LS):.2f}, 95th percentile {np.nanpercentile(LS, 95):.1f}")

LS over the catchment: median 2.56, 95th percentile 12.1

C — cover and management¶

Dense vegetation protects soil; bare ground and cultivated land do not. C values per land-cover class follow ranges reported in the literature (e.g. Panagos et al., 2015). Edit the table to reflect local knowledge.

In [6]:
# @title 5 · Cover factor from ESA WorldCover 2021
C_TABLE = {10: 0.003,   # tree cover
           20: 0.03,    # shrubland
           30: 0.05,    # grassland
           40: 0.20,    # cropland
           50: 0.0,     # built-up (sealed surfaces do not erode)
           60: 0.35,    # bare / sparse vegetation
           70: 0.0, 80: 0.0, 90: 0.0, 95: 0.003, 100: 0.10}   # snow, water, wetland, mangroves, moss/lichen
catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=pc.sign_inplace)
wc_items = list(catalog.search(collections=["esa-worldcover"], bbox=list(catchment.total_bounds),
                               query={"esa_worldcover:product_version": {"eq": "2.0.0"}}).items())
lc = load(wc_items, bands=["map"], geobox=gbox, resampling="mode")["map"].isel(time=0).values
Cf = np.vectorize(lambda v: C_TABLE.get(int(v), np.nan))(lc).astype("float32")
Cf[~inside] = np.nan
shares = pd.Series(lc[inside]).value_counts(normalize=True) * 100
names = {10: "Tree cover", 20: "Shrubland", 30: "Grassland", 40: "Cropland", 50: "Built-up", 60: "Bare / sparse",
         80: "Water", 90: "Wetland"}
print("Land cover in the catchment: " + ", ".join(f"{names.get(k, k)} {v:.0f} %" for k, v in shares.items() if v >= 1))

Land cover in the catchment: Cropland 55 %, Grassland 34 %, Tree cover 5 %, Built-up 3 %, Shrubland 1 %

A — estimated soil loss, and where it concentrates¶

In [7]:
# @title 6 · Soil-loss map and erosion classes
A = R * K * LS * Cf * 1.0                       # P = 1
A[~inside] = np.nan
classes = [(0, 2, "very low"), (2, 5, "low"), (5, 10, "moderate"), (10, 20, "high"), (20, np.inf, "very high")]
px_ha = (RES * RES) / 1e4
tab = pd.DataFrame([{"class": n, "t/ha/yr": f"{lo:g}–{hi:g}" if np.isfinite(hi) else f"> {lo:g}",
                     "area (km²)": round(np.nansum((A >= lo) & (A < hi)) * px_ha / 100),
                     "share of area (%)": round(100 * np.nansum((A >= lo) & (A < hi)) / inside.sum(), 1),
                     "share of soil loss (%)": round(100 * np.nansum(np.where((A >= lo) & (A < hi), A, 0)) / np.nansum(A), 1)}
                    for lo, hi, n in classes]).set_index("class")
display(tab)
fig, ax = plt.subplots(figsize=(10, 7))
bb = gbox.boundingbox
im = ax.imshow(np.log10(np.clip(A, 0.1, 200)), cmap="YlOrBr", vmin=-1, vmax=np.log10(200),
               extent=(bb.left, bb.right, bb.bottom, bb.top))
c_utm.boundary.plot(ax=ax, color="#0F172A", lw=1)
ax.set_title(f"{info['name']} catchment — estimated soil loss (RUSLE)"); ax.axis("off")
cb = fig.colorbar(im, ax=ax, shrink=0.75); cb.set_ticks([-1, 0, 1, 2]); cb.set_ticklabels(["0.1", "1", "10", "100"])
cb.set_label("t ha⁻¹ yr⁻¹ (log scale)")
plt.show()
a = np.sort(A[np.isfinite(A)])[::-1]                 # soil loss, largest first
share10 = 100 * a[: max(1, len(a) // 10)].sum() / a.sum()
print(f"The most erosion-prone 10 % of the catchment produces about {share10:.0f} % of the estimated soil loss: "
      f"that is where conservation pays off most.")
t/ha/yr area (km²) share of area (%) share of soil loss (%)
class
very low 0–2 1962 25.9 2.8
low 2–5 1673 22.1 10.4
moderate 5–10 1820 24.0 23.8
high 10–20 1427 18.8 36.1
very high > 20 540 7.1 26.9
No description has been provided for this image 
The most erosion-prone 10 % of the catchment produces about 34 % of the estimated soil loss: that is where conservation pays off most.

From soil loss to reservoir sediment (order of magnitude)¶

Only part of the eroded soil reaches the reservoir; the rest is deposited on slopes and floodplains. A simple sediment delivery ratio that decreases with catchment size (Vanoni, 1975) gives an order-of-magnitude estimate of the sediment reaching the dam each year.

In [8]:
# @title 7 · Indicative sediment yield
area_km2 = inside.sum() * RES * RES / 1e6
gross_t = np.nansum(A) * px_ha
sdr = 0.42 * (area_km2 / 2.59) ** -0.125         # Vanoni (1975), area in square miles
yield_t = gross_t * sdr
volume_m3 = yield_t / 1.3                         # assumed dry bulk density of deposited sediment, 1.3 t/m³
print(f"Gross erosion ≈ {gross_t / 1e6:,.1f} million t/yr · delivery ratio ≈ {sdr:.2f} · "
      f"sediment reaching the reservoir ≈ {yield_t / 1e6:,.1f} million t/yr ≈ {volume_m3 / 1e6:,.1f} million m³/yr")
print("Treat this as an order of magnitude: compare it with bathymetric surveys of the reservoir where they exist.")

Gross erosion ≈ 5.5 million t/yr · delivery ratio ≈ 0.15 · sediment reaching the reservoir ≈ 0.9 million t/yr ≈ 0.7 million m³/yr
Treat this as an order of magnitude: compare it with bathymetric surveys of the reservoir where they exist.

Limits and next steps¶

  • RUSLE estimates sheet and rill erosion; gullies, landslides and riverbank erosion — often large sources in the Maghreb — are not included.
  • A single catchment-average R ignores the spatial pattern of storms; K, C and P values carry their own uncertainty. Use the map to rank areas, then verify in the field.
  • Combine the hot spots with the AI change radar from the live notebook: land-cover change in high-risk zones is where siltation can accelerate.

Data: Copernicus DEM GLO-90 (© DLR e.V. and © Airbus Defence and Space GmbH, provided under COPERNICUS by the European Union and ESA) · SoilGrids 2.0 (ISRIC — World Soil Information, CC-BY 4.0) · ESA WorldCover 2021 (© ESA WorldCover project, CC-BY 4.0) · CHIRPS v2.0 (Funk et al. 2015) · HydroBASINS (Lehner & Grill 2013). Methods: Renard & Freimund (1994); Williams (1995); Moore & Burch (1986); Panagos et al. (2015); Vanoni (1975).

QuestFeed · IC-TWEP 2026 take-home materials · questfeed.com