Imported from lzwei196/KISS-Knowledge-Infrastructure-for-Scientific-Simulation (models/Landlab/SKILL.md). Install upstream with npx skills add lzwei196/KISS-Knowledge-Infrastructure-for-Scientific-Simulation --skill Landlab. Copyright stays with the author.

MANDATORY EXECUTION POLICY — READ BEFORE PROCEEDING

You MUST run the actual model binary or package described in this document. If the model fails to import, compile, or execute, you MUST:

1. Check diagnostics/triplets.yaml for a matching error pattern
2. Try the fix suggested in the triplet's remedy section
3. If still failing, report the error to the user with full details

You MUST NOT substitute a simplified Python formula, regression equation, or hand-coded approximation in place of the real model.

Before starting, run: python preflight_check.py (in this KI directory) to verify that the model binary/package and required data are available.

DEBUGGING PROTOCOL — When something goes wrong, follow this order:

1. Check triplets — diagnostics/triplets.yaml may already cover this error
2. Read official docs — The model's own documentation for expected formats/units
3. Find working examples — Check outputs/ or the model's shipped test data
4. Fix the tool — With knowledge of what "correct" looks like

Do NOT write custom debug scripts. The answers are in the docs and examples.

KI map — what to read, and when

Projected 2026-08-17 from the KI's actual contents — 9 components present. Refresh: python3 ki_tools_common/generate_skill_map.py --ki_dir <this KI>.

Executable tool index (projected — complete by construction)

Every public tool in this KI, by exact path. What each is FOR lives in the human-written Tool Inventory above; --help on any of these prints its arguments.

8 public tools; _-prefixed helpers and packaging files excluded.

Landlab v2.10 — Earth Surface Dynamics Modeling Toolkit

Package: hydrocraft-landlab v1.0.0 Model: Landlab 2.10.1 (Python library with Cython extensions) Domain: Geomorphology, hydrology, stratigraphy, glaciology Created by: Landlab Development Team (CU Boulder, U Washington, Tulane) Last updated: 2026-04-29 Stats: 7 tools | 5 skill documents | 22 diagnostic triplets | ~2,200 lines of validated code Validation status: production_validated (5 test cases: Whipple & Tucker 1999 analytical + SPACE binary Qs-Q + SPACE steady-state concavity + Loess Plateau SRTM real-DEM slope-area + Loess Plateau stream power SY + P export)

Data Preparation

Forcing data

Data Sources: Use from ki_tools_common.load_forcing import load_daily_forcing for CMFD/MSWX/NASA POWER.

Data Validation Reference: See data_ki/CMFD/SKILL.md for meteorological forcing documentation. See data_ki/USGS_Sediment/SKILL.md for suspended sediment observations.

Applicability

This KI is for landscape evolution and sediment-yield benchmarks only. It does NOT produce a daily discharge hydrograph and is NOT applicable to gauge-discharge tests (e.g. Bengbu, Wangjiaba). The surface_water__discharge field from FlowAccumulator is steady-state Q = drainage_area * uniform_runoff_rate, not a time series. Validation cases supported: slope-area, concavity, sediment yield, Whipple-Tucker analytical, SPACE Qs-Q. For streamflow benchmarks use a hydrologic KI (PIHM, mHM, VIC, etc.).

Overview

Landlab is an open-source Python package for numerical modeling of Earth surface dynamics. Unlike compiled, monolithic landscape evolution models, Landlab is a component-based framework: users compose simulations by coupling reusable process components (erosion, diffusion, flow routing, weathering, etc.) on a shared model grid. The grid stores spatially distributed fields (elevation, drainage area, soil depth) at nodes, links, patches, cells, corners, and faces.

Landlab supports six grid types: RasterModelGrid (uniform rectangular), HexModelGrid (hexagonal cells), VoronoiDelaunayGrid (irregular Voronoi), FramedVoronoiGrid, RadialModelGrid, NetworkModelGrid, and IcosphereGlobalGrid. All grids implement the same field-storage API, so components are grid-agnostic.

There are 60+ process components covering: flow routing (D8, D-infinity, MFD, priority-flood), stream-power erosion, hillslope diffusion, SPACE bedrock- alluvium erosion, overland flow (kinematic wave, de Almeida), flexural isostasy, tectonic faulting, soil production (exponential weathering), fire, vegetation dynamics, species evolution, tidal flow, marine sediment transport, and more.

Installation

From PyPI (recommended)

python -m venv venv && source venv/bin/activate
pip install landlab

From source (development)

git clone https://github.com/landlab/landlab.git
cd landlab
pip install -e ".[dev,testing]"

Dependencies

• Core: numpy>=1.20, scipy, matplotlib, xarray>=0.16, pandas, pyyaml, netcdf4, bmipy, rich-click
• Build: cython, setuptools, wheel (Cython extensions for performance)
• Optional: pyshp (shapefiles), statsmodels, gflex (flexure)

Verify installation

import landlab; print(landlab.__version__)
from landlab import RasterModelGrid
mg = RasterModelGrid((10, 10), xy_spacing=100.0)
print(f"Grid OK: {mg.number_of_nodes} nodes")

Pipeline (7 stages)

Parallelism

• Stages 2 and 3 can run in parallel (independent input preparation)
• Stage 4 depends on both 2 and 3
• Stages 5–7 are sequential

Execution Model

Landlab is a Python API, not a CLI binary. Simulations are Python scripts:

from landlab import RasterModelGrid
from landlab.components import FlowAccumulator, StreamPowerEroder, LinearDiffuser

# 1. Create grid
mg = RasterModelGrid((50, 50), xy_spacing=100.0)
z = mg.add_zeros("topographic__elevation", at="node")
z += mg.node_x * 0.01 + np.random.rand(mg.number_of_nodes)

# 2. Set boundary conditions
mg.set_closed_boundaries_at_grid_edges(True, False, True, False)

# 3. Instantiate components
fa = FlowAccumulator(mg, flow_director="FlowDirectorD8")
sp = StreamPowerEroder(mg, K_sp=1e-5, m_sp=0.5, n_sp=1.0)
ld = LinearDiffuser(mg, linear_diffusivity=0.01)

# 4. Time loop
dt = 1000.0  # years
for t in range(2000):
    fa.run_one_step()
    sp.run_one_step(dt)
    ld.run_one_step(dt)

# 5. Output
from landlab.io import write_esri_ascii
write_esri_ascii("output.asc", mg, "topographic__elevation")

BMI Interface

Components can also be driven via the Basic Modeling Interface (BMI):

from landlab.bmi import wrap_as_bmi
BmiSPE = wrap_as_bmi(StreamPowerEroder)
model = BmiSPE()
model.initialize("config.yaml")
model.update()
model.finalize()

6. Output Description

Source of truth: dag.yaml. The DAG wins over this body if they ever disagree.

Headline output (validation_rank: 1):

sediment_yield — Annual catchment sediment yield via detachment-limited stream power (E = K_sp A^m S^n converted by bulk density). (t/km^2/yr)

The other DAG outputs are: topographic__elevation, concavity (theta), steepness (ks), relief, sediment__outflux (SPACE Qs), soil__depth, and surface_water__discharge (steady-state).

8. Unit Table

Source of truth: dag.yaml and docs/format_spec.yaml. Landlab is unit-agnostic, so the tools and user scripts must keep the unit system consistent.

11. Validated Results

Source of truth for validation bars: docs/validation_convention.yaml. Null convention bands are written as no cited threshold; no uncited threshold is substituted.

Performance Metrics -- judged against the field's bar

Documented validation cases

Data replacement tracking

Grid Field System

Fields are stored at 7 locations on the grid:

Standard Variable Names and Units

Unit Trap Table

These are the most dangerous unit mismatches in Landlab workflows:

Key Components Reference

Flow Routing

Erosion & Sediment Transport

Weathering & Soil Production

Hydrology

Tectonics & Lithology

I/O Formats

Critical Domain Knowledge

1. Landlab is unit-agnostic — it does NOT enforce units internally. The user MUST ensure all inputs share a consistent unit system. If elevation is in meters, spacing must be meters, K must be in m^(1-2m)/yr, diffusivity in m^2/yr, and dt in years. Mixing m/s with m/yr is the #1 source of silent failures. → dt_001, dt_005

2. K_sp units depend on m_sp and n_sp exponents — for the stream power law E = K A^m S^n, when m=0.5 and n=1.0, K has units m^(1-2m)/time = m^0/yr = 1/yr. Changing exponents changes K units. Copying K values from papers with different exponents produces wrong erosion rates. → dt_003

3. Boundary conditions must be set before component instantiation — FlowAccumulator needs at least one open boundary node for drainage. If all boundaries are closed, flow has nowhere to go and drainage_area will be zero everywhere. → dt_011

4. Depression routing is not automatic — D8 flow routing stops at pits. Without DepressionFinderAndRouter or PriorityFloodFlowRouter, interior pits accumulate flow and block downstream erosion. → dt_012

5. CFL stability for explicit diffusion — LinearDiffuser uses an explicit scheme with internal sub-stepping (alpha=0.15 × dx²/D). If the user wraps it in their own sub-stepping, they may defeat the internal CFL guard. → dt_013

6. Grid spacing affects slope calculation — Slope = Δz/Δx. If xy_spacing is in degrees instead of meters, slopes will be orders of magnitude wrong, producing unrealistic erosion. → dt_004

7. FlowAccumulator must run before erosion components — StreamPowerEroder, ErosionDeposition, and Space all require drainage_area, flow__receiver_node, and flow__upstream_node_order to be populated. These are outputs of FlowAccumulator. Running erosion first yields zero drainage area. → dt_014

Validation: SPACE Binary — Atchafalaya Qs-Q + Concavity (2026-04-29)

Configuration — Test 1 (Qs-Q binary)

• Grid: 3×51 quasi-1D channel, dx=500 m, slope=5×10⁻⁴, soil_depth=2 m
• Component: SpaceLargeScaleEroder (K_sed=2.5e-5, m_sp=0.5, n_sp=1.0, v_s=5 m/yr, H*=1 m)
• Method: warmup 800 steps × 1 yr; then 7 runoff-rate probes (single SPACE step, no uplift)
• Fit: log(Qs) vs log(Q) — sediment mass flux vs water flux

Configuration — Test 2 (Steady-state concavity)

• Grid: 50×50 RasterModelGrid, dx=200 m, open south boundary
• Duration: 4000 × 500 yr = 2 Myr; uplift U=1e-3 m/yr

Results

Key Findings

1. b_sim=0.897 reflects transport-limited SPACE behavior (86.5% sediment-dominated, long travel distance) — physically correct for these parameters
2. Detachment-limited regime (b→m=0.5) requires bare bedrock (H→0) or very high v_s; transport-limited (b→1.0) matches Atchafalaya obs b=1.077
3. Fitting Qs vs Q (not SSC vs Q) is essential — SSC = Qs/Q ∝ Q^(m-1) = Q^(-0.5) always gives negative exponent regardless of model correctness
4. Steady-state θ=0.464 confirms SPACE correctly implements stream-power scaling (Whipple & Tucker 1999 theoretical θ=0.5)

Output

• Figure: outputs/landlab_atchafalaya_ssc_q/validation_figure.png
• Metrics: outputs/landlab_atchafalaya_ssc_q/metrics.json
• Run: python tools/dissect_atchafalaya_ssc_q_surrogate.py

Validation: Loess Plateau Real-DEM Slope-Area (2026-04-29)

Configuration

• DEM: SRTM 30m, tile N36E109 (~36.64°N 109.33°E, Shaanxi, China)
• Clip: 500×500 px (15 km × 15 km), dx=30 m, z=[953, 1369] m
• Tool: convert_dem_to_grid.py → FlowAccumulator(D8 + DepressionFinderAndRouter)
• Slope method: 2D spatial gradient |∇z| (binned; D8 receiver slope is noisy on SRTM)
• Channel threshold: A ≥ 0.1 km²; 25 log-spaced area bins

Results vs Published Loess Plateau Concavity

Key Findings

1. Binned slope-area R²=0.866 confirms Landlab correctly extracts channel scaling from real terrain
2. θ=0.163 is consistent with actively eroding Loess Plateau gullies (lower than global mountain-river average θ≈0.45)
3. D8 receiver slope (topographic__steepest_slope) gives θ≈0.13 and R²≈0.08 on SRTM — use spatial gradient instead (dt_018)
4. Single-outlet boundary (closed_all_but_outlet) depresses θ toward 0 — use open_all for DEM clips (dt_019)
5. SRTM integer elevations create discrete slope bands; log-spaced binning recovers R²=0.87 from a raw scatter of R²=0.07

Triplets exercised

• dt_018: D8 slope artifact on filled DEMs → use spatial gradient
• dt_019: single-outlet boundary depresses θ → use open_all

Validation: Loess Plateau Sediment Yield + Particulate-P Export (2026-04-29)

Configuration

• DEM: SRTM 30m, tile N36E109, clip 500×500 px (15 km × 15 km), dx=30 m, z=[953, 1369] m
• Model: Detachment-limited stream power law E = K_sp × A^0.5 × S (no SPACE; see dt_022)
• Boundary: closed_all_but_outlet (lowest boundary node = outlet, z=959 m, node 189499)
• Flow routing: FlowAccumulator(D8 + DepressionFinderAndRouter), runoff_rate=3.17e-9 m/s (100 mm/yr)
• K_sp calibration: K_sp = SY_target / (mean(A^0.5 × S) × ρ_bulk × 1e6) → 1.64e-4 m^0/yr
• Bulk density: 1400 kg/m³; Soil P: 800 mg/kg; P enrichment ratio: 2.0
• P export: TP = SY × 10 × P_soil_ppm × ER × 1e-6

Results vs Published Yellow River Tributary Yields

Key Findings

1. SPACE is unsuitable for annual SY from real DEMs — transport-limited regime gives ~2.7% export efficiency; use DL stream power instead (dt_022)
2. K_sp calibrated analytically so mean gross erosion = target SY; no iterative tuning needed
3. Spatial gradient |∇z| (np.gradient on 2D elevation array) gives physically meaningful slopes for SY; D8 receiver slope is too noisy at cell scale (dt_018)
4. closed_all_but_outlet is mandatory for watershed-scale SY (single catchment definition); open_all creates many sub-basins and underestimates outlet flux
5. P export (TP=134.25 kg/ha/yr) is gross potential; actual at watershed outlet is 0.3–0.7× lower due to within-basin P retention

Triplets exercised

• dt_022: SPACE near-zero SY on real DEMs → use DL stream power for annual SY

Output

• Figure: outputs/landlab_loess_sediment/validation_figure.png
• Metrics: outputs/landlab_loess_sediment/metrics.json
• Run: python tools/dissect_loess_plateau_sediment_yield.py

Validation: Whipple & Tucker (1999) Steady-State Test (2026-03-25)

Configuration

• Grid: 50×50 RasterModelGrid, dx=100 m
• Components: FlowAccumulator(D8) + StreamPowerEroder(K=1e-5, m=0.5, n=1.0) + LinearDiffuser(D=0.01 m²/yr)
• Uplift: 1e-3 m/yr uniform
• Duration: 2 Myr, dt=500 yr (4000 steps)
• Boundary: South edge open, other three closed

Results vs Analytical Solution

Key Findings

1. Slope-area correlation r=0.9978 confirms steady-state channel profiles
2. Concavity index θ=0.488 vs theoretical 0.500 (2.4% error, within discretization)
3. Relief stabilized at 810.2 m after 2 Myr of evolution
4. All three components (FlowAccumulator, StreamPowerEroder, LinearDiffuser) validated
5. Wall time 2.2s for 4000 steps on 50×50 grid demonstrates computational efficiency

Tools Reference

Calibration Parameters (Priority Order)

Data Requirements

Quick Start

# Install
pip install landlab

# Run input converter
python tools/convert_dem_to_grid.py --dem input.asc --output grid.nc --spacing 100

# Run soil parameter converter
python tools/convert_soil_params.py --hwsd soil.csv --grid grid.nc --output params.json

# Run simulation
python tools/run_landlab.py --config config.yaml --output results/

# Parse output
python tools/parse_landlab_output.py --input results/ --output results/summary.csv --metrics slope_area

Diagnostic Triplets Summary

Reference: diagnostics/triplets.yaml

File Structure

ki/
├── SKILL.md                          # This file
├── tools/
│   ├── convert_dem_to_grid.py                    # DEM → Landlab grid
│   ├── convert_soil_params.py                    # Soil data → grid fields
│   ├── run_landlab.py                            # Execution wrapper
│   ├── parse_landlab_output.py                   # Output → CSV + metrics
│   ├── dissect_atchafalaya_ssc_q_surrogate.py    # SPACE Qs-Q + concavity validation
│   ├── dissect_loess_plateau_slope_area.py       # Real-DEM slope-area validation
│   └── dissect_loess_plateau_sediment_yield.py   # DL stream power SY + P export validation
├── docs/
│   ├── s1_grid_setup.md              # Grid creation skill
│   ├── s2_input_preparation.md       # DEM/forcing loading
│   ├── s3_soil_parameters.md         # Soil/weathering params
│   ├── s4_component_assembly.md      # Coupling components
│   └── s5_execution_output.md        # Running and analyzing
└── diagnostics/
    └── triplets.yaml                 # 22 symptom→diagnosis→remedy