Landlab git:20261024.0b1595dC

Immutable. This exact content never changes and is served at /api/v1/blob/b4fcc8da984c8252.

> **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:BEGIN (projected by generate_skill_map.py — edit the KI, not this table) -->
## KI map — what to read, and when

| when you need | read | why |
|---|---|---|
| FIRST, always | `preflight_check.py` | run it (`python preflight_check.py`): proves env/binary/data are usable and emits a machine-readable `PREFLIGHT_REPORT=` line. Do not debug a run that never had a healthy environment. |
| to run the pipeline stages | `tools/` (8 tools) | the executable pipeline. Read each tool's argparse (`--help`) before composing a command; SKILL.md's stage table says which tool serves which stage. |
| before running a stage | `docs/s*_*.md` (5 stage docs) | per-stage procedure, verification and traps — the how-to that SKILL.md's overview compresses. |
| on ANY error, before debugging | `diagnostics/triplets.yaml` (22 entries) | symptom → diagnosis → remedy for this model's known failure modes. Check here FIRST; the answer usually exists. Never renumber or rewrite entries. |
| to know what an output IS | `dag.yaml` | the model's identity: every output's medium, units, `validation_rank` (1 = the headline variable) and observability. Scoring and obs-binding read THIS — when asked 'what does this model predict', the dag is the answer, not a guess. |
| when building inputs / parsing outputs | `docs/format_spec.yaml` | exact I/O shapes + `known_issues`, projected from dag + triplets. Regenerate with `ki_tools_common/generate_format_spec.py` after changing either — never hand-edit. |
| to judge a run's skill | `docs/validation_convention.yaml` | how this model's field judges it validated: per-`dag_variable` metrics, directions and CITED pass-bands. A run is graded against these, not against intuition. |
| for claims and thresholds | `docs/gathered_papers.json` (22 papers) + `docs/papers_index.md` | the literature this KI is judged by; each entry's `text_path` is fetched full text in the central paper cache. `role: benchmark` marks the model's own skill paper. |
| for a machine-readable summary | `knowledge_infrastructure.yaml` | the manifest (package, pipeline, validation tier, counts) — projected by `ki_tools_common/generate_ki_manifest.py`; regenerate after structural changes, never hand-edit. |

*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>`.*
<!-- KI-MAP:END -->

<!-- KI-TOOL-INDEX:BEGIN (projected by generate_skill_map.py — the discoverability contract: every public tool, exact path; PURPOSE stays human-authored elsewhere) -->
### 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.

| tool (exact path) | invocation |
|---|---|
| `tools/convert_dem_to_grid.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/convert_dem_to_grid.py --help` |
| `tools/convert_soil_params.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/convert_soil_params.py --help` |
| `tools/dissect_atchafalaya_ssc_q_surrogate.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/dissect_atchafalaya_ssc_q_surrogate.py --help` |
| `tools/dissect_loess_plateau_sediment_yield.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/dissect_loess_plateau_sediment_yield.py --help` |
| `tools/dissect_loess_plateau_slope_area.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/dissect_loess_plateau_slope_area.py --help` |
| `tools/dissect_space_ssc_q_rating.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/dissect_space_ssc_q_rating.py --help` |
| `tools/parse_landlab_output.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/parse_landlab_output.py --help` |
| `tools/run_landlab.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/run_landlab.py --help` |

*8 public tools; `_`-prefixed helpers and packaging files excluded.*
<!-- KI-TOOL-INDEX:END -->

# 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**: current tool, stage-doc and triplet counts: the KI map at the top of this file
**Validation status**: production_validated

---

## 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)
```bash
python -m venv venv && source venv/bin/activate
pip install landlab
```

### From source (development)
```bash
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
```python
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)

| Stage | ID | Tool | Description |
|-------|----|------|-------------|
| 1. Grid Setup | s1_grid | — | Create model grid with DEM or synthetic topography |
| 2. Input Preparation | s2_input | `convert_dem_to_grid.py` | Load DEM, set boundary conditions, add fields |
| 3. Soil/Parameter Setup | s3_params | `convert_soil_params.py` | Map HWSD or other soil data to grid fields |
| 4. Component Assembly | s4_assembly | — | Instantiate and couple Landlab components |
| 5. Execution | s5_run | `run_landlab.py` | Time-stepping loop with run_one_step() |
| 6. Output Extraction | s6_output | `parse_landlab_output.py` | Extract fields to CSV/NetCDF, compute metrics |
| 7. Diagnostics | s7_diag | — | Validate results against analytical solutions |

### 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:

```python
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):
```python
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`)

| Output variable (dag `var`) | Rank | Unit | Emitted in | Description |
|-----------------------------|------|------|------------|-------------|
| `sediment_yield` | 1 | `t/km^2/yr` | derived metric (detachment-limited stream-power export) | Annual catchment sediment yield via detachment-limited stream power (E = K_sp A^m S^n converted by bulk density). |
| `sediment__outflux (SPACE Qs)` | 2 | `m^3/s (or kg/s)` | Landlab grid field at_node (SPACE / ErosionDeposition) | Sediment flux out of a node; the Qs in the SPACE Qs-Q rating relationship. |
| `surface_water__discharge (steady-state)` | 3 | `m^3/s` | Landlab grid field at_node (FlowAccumulator) | Steady-state discharge Q = drainage_area * runoff_rate from FlowAccumulator; NOT a hydrograph time series. |
| `topographic__elevation` | 4 | `m` | Landlab grid field at_node; written to ESRI ASCII (.asc) / NetCDF (.nc) | Final land-surface elevation of the evolved landscape (node field). |
| `concavity (theta)` | 5 | `dimensionless` | derived metric from slope-area regression on the grid | Land-surface/channel-topography slope-area concavity index theta = -d(logS)/d(logA); at steady state expected to approach m_sp/n_sp. |
| `steepness (ks)` | 6 | `dimensionless (channel steepness index)` | derived metric from slope-area regression | Channel steepness index from the slope-area regression. |
| `relief` | 7 | `m` | derived metric over core nodes | Max-minus-min elevation over core nodes; landscape relief at steady state. |
| `soil__depth` | 8 | `m` | Landlab grid field at_node | Regolith/soil thickness above bedrock (node field). |

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.

| Variable / quantity | Source unit handled by this KI | Model or DAG unit | Conversion / invariant | Where checked or used |
|---------------------|--------------------------------|-------------------|------------------------|-----------------------|
| `sediment_yield` | derived from detachment-limited stream-power export | `t/km^2/yr` | E = K_sp A^m S^n converted by bulk density | DAG rank-1 output; `tools/dissect_loess_plateau_sediment_yield.py` |
| `topographic__elevation` | DEM elevation in meters; feet are a known trap | `m` | feet to meters: multiply by 0.3048 | `convert_dem_to_grid.py`; dt_007 |
| `xy_spacing` / grid geometry | projected metric spacing required | `m` | reproject geographic degrees to a metric CRS before loading | `convert_dem_to_grid.py`; dt_004 |
| `water__unit_flux_in` | rainfall / runoff from forcing data | `m/s` | mm/hr to m/s: divide by 3,600,000; mm/day to m/s: divide by 86,400,000; m/yr to m/s: divide by 31,557,600 | `FlowAccumulator`; dt_001 |
| `surface_water__discharge (steady-state)` | derived steady-state accumulated runoff | `m^3/s` | Q = drainage_area * runoff_rate | `FlowAccumulator`; DAG rank-3 output |
| `soil__depth` | HWSD-style soil depth commonly in cm | `m` | cm to m: divide by 100 | `convert_soil_params.py`; dt_006 |
| `K_sp` | literature or calibrated stream-power erodibility | `length^(1-2*m_sp) / time` | units depend on `m_sp` and `n_sp`; at `m_sp=0.5`, `n_sp=1.0`, K has units `1/yr` when dt is in years | `StreamPowerEroder`; dt_003 |
| `linear_diffusivity` | literature or calibrated diffusivity | `m^2/yr` when dt is in years | m^2/s to m^2/yr: multiply by 3.156e7 | `LinearDiffuser`; dt_002 |
| `dt` | user-chosen timestep | `yr` for landscape evolution, `s` for overland flow | dt must share the time unit of K, diffusivity, uplift, and other rates | `run_landlab.py`; dt_005 |
| `uplift_rate` | scalar or field forcing | `m/yr` | applied as `z[core_nodes] += uplift_rate * dt` when dt is in years | user execution loop |
| `sediment__outflux (SPACE Qs)` | SPACE / ErosionDeposition node flux | `m^3/s (or kg/s)` | compare Qs directly against Q; do not fit SSC = Qs/Q when testing SPACE Qs-Q | DAG rank-2 output; dt_020 |
| `rainfall__daily_depth` | SoilMoisture rainfall can arrive in m or m/s | `mm` at cell | m to mm: multiply by 1000; m/s to mm/day: multiply by 86,400,000 | SoilMoisture; dt_008 |

---

## 11. Validation Status

Run results (site, period, scores) are not kept here. They live in the run records
(`test_runs` in hydrocraft.db and each run's own result files).

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

| DAG variable | Headline metric | Direction | Convention bar, cited |
|--------------|-----------------|-----------|------------------------|
| `sediment_yield` | `pbias` | zero_centered | very good: `|PBIAS| <= 10` (moriasi2015); good: `|PBIAS| <= 15` (moriasi2015); satisfactory: `|PBIAS| <= 20` (moriasi2015) |
| `topographic__elevation` | `csi` | maximize | very good: no cited threshold; good: no cited threshold; satisfactory: no cited threshold; cites: `[]` |
| `concavity (theta)` | `pbias` | zero_centered | very good: no cited threshold; good: no cited threshold; satisfactory: no cited threshold; cites: `[]` |
| `steepness (ks)` | `pbias` | zero_centered | very good: no cited threshold; good: no cited threshold; satisfactory: no cited threshold; cites: `[]` |

### Data replacement tracking

| Component | Source | Status | Notes |
|-----------|--------|--------|-------|
| Forcing | CMFD/MSWX/NASA POWER via `ki_tools_common.load_forcing.load_daily_forcing` where meteorological forcing is needed | Available | Landscape-evolution sediment-yield runs commonly use uniform runoff_rate, not a daily hydrograph. |
| Soil | HWSD or other soil data via `convert_soil_params.py` | Available | Convert depth from cm to m before Landlab fields. |
| DEM / topography | SRTM/ASTER/LiDAR through `convert_dem_to_grid.py` | Available | DEM must be projected in meters and NoData must be closed. |
| Initial conditions | Landlab grid fields (`topographic__elevation`, `soil__depth`, `bedrock__elevation`) | Available | Field locations and units must match component expectations. |

---

## Grid Field System

Fields are stored at 7 locations on the grid:

| Location | Description | Access Pattern |
|----------|-------------|----------------|
| `node` | Grid vertices (most common) | `grid.at_node["field_name"]` |
| `link` | Edges connecting nodes | `grid.at_link["field_name"]` |
| `patch` | Polygons in primal mesh | `grid.at_patch["field_name"]` |
| `cell` | Polygons in dual mesh | `grid.at_cell["field_name"]` |
| `corner` | Vertices in dual mesh | `grid.at_corner["field_name"]` |
| `face` | Edges in dual mesh | `grid.at_face["field_name"]` |
| `grid` | Scalar (whole-grid) values | `grid.at_grid["field_name"]` |

### Standard Variable Names and Units

| Variable | Units | Location | Description |
|----------|-------|----------|-------------|
| `topographic__elevation` | m | node | Land surface elevation |
| `drainage_area` | m^2 | node | Contributing upstream area |
| `surface_water__discharge` | m^3/s | node | Volumetric water discharge |
| `water__unit_flux_in` | m/s | node | Rainfall rate (depth/time) |
| `topographic__steepest_slope` | - | node | Max downhill gradient |
| `flow__receiver_node` | - | node | Downstream neighbor ID |
| `flow__upstream_node_order` | - | node | Topological sort order |
| `flow__link_to_receiver_node` | - | node | Link ID to receiver |
| `soil__depth` | m | node | Regolith/soil thickness |
| `soil_production__rate` | m/yr | node | Bedrock-to-soil conversion rate |
| `bedrock__elevation` | m | node | Bedrock surface elevation |
| `hillslope_sediment__unit_volume_flux` | m^2/s | link | Diffusive sediment flux |
| `topographic__gradient` | - | link | Surface slope on links |
| `sediment__influx` | m^3/s | node | Incoming sediment rate |
| `sediment__outflux` | m^3/s | node | Outgoing sediment rate |
| `surface_water__depth` | m | node | Water depth on surface |
| `rainfall__daily_depth` | mm | cell | Daily precipitation |
| `vegetation__cover_fraction` | - | cell | Fractional veg cover |

---

## Unit Trap Table

These are the most dangerous unit mismatches in Landlab workflows:

| # | Variable | Expected | Common Mistake | Effect | Triplet |
|---|----------|----------|----------------|--------|---------|
| 1 | `water__unit_flux_in` | m/s | mm/hr or m/yr | Discharge 1000x–31.5M× wrong | dt_001 |
| 2 | `linear_diffusivity` | m^2/yr (if dt in yr) | m^2/s | Hillslopes flatten instantly | dt_002 |
| 3 | `K_sp` (erodibility) | depends on m,n | wrong exponent combo | Erosion rate orders of magnitude off | dt_003 |
| 4 | `xy_spacing` | m | km or degrees | Area, slope, flux all wrong | dt_004 |
| 5 | `dt` (timestep) | must match K units | yr vs s mismatch | All rates wrong | dt_005 |
| 6 | `soil__depth` | m | cm | Weathering rate, SPACE output wrong | dt_006 |
| 7 | `topographic__elevation` | m | ft or cm | Slope calc wrong → erosion wrong | dt_007 |
| 8 | `rainfall__daily_depth` | mm | m | Soil moisture 1000× off | dt_008 |
| 9 | `roughness` (Manning n) | s/m^(1/3) | 1/n (Chezy) | Velocity inverted | dt_009 |
| 10 | DEM nodata value | -9999 | 0 or NaN | Flow routes into nodata cells | dt_010 |

---

## Key Components Reference

### Flow Routing
| Component | Method | Key Params |
|-----------|--------|------------|
| `FlowAccumulator` | Wraps director + accumulation | `flow_director`, `runoff_rate` |
| `FlowDirectorD8` | Steepest single-direction (D8) | — |
| `FlowDirectorDINF` | D-infinity (two receivers) | — |
| `FlowDirectorMFD` | Multiple flow directions | `partition_method` |
| `DepressionFinderAndRouter` | Fill/route pits | `routing` |
| `PriorityFloodFlowRouter` | Priority-flood (Barnes 2014) | `flow_metric`, `suppress_out` |

### Erosion & Sediment Transport
| Component | Equation | Key Params |
|-----------|----------|------------|
| `StreamPowerEroder` | E = K A^m S^n | `K_sp`, `m_sp`, `n_sp`, `threshold_sp` |
| `ErosionDeposition` | Davy & Lague (2009) ξ-q | `K`, `v_s`, `m_sp`, `n_sp` |
| `Space` | Bedrock + alluvium layers | `K_sed`, `K_br`, `v_s`, `H_star` |
| `LinearDiffuser` | q = -D ∇z | `linear_diffusivity` |
| `DepthDependentDiffuser` | Nonlinear depth-dep. | `linear_diffusivity`, `soil_transport_decay_depth` |
| `TaylorNonlinearHillslopeFlux` | Taylor expansion nonlinear | `critical_slope`, `nterms` |

### Weathering & Soil Production
| Component | Equation | Key Params |
|-----------|----------|------------|
| `ExponentialWeatherer` | P = P0 exp(-H/H*) | `soil_production_maximum_rate`, `soil_production_decay_depth` |
| `DepthDependentTaylorDiffuser` | Coupled soil transport | `soil_transport_velocity`, `soil_transport_decay_depth` |

### Hydrology
| Component | Method | Key Params |
|-----------|--------|------------|
| `KinwaveOverlandFlowModel` | Kinematic wave + Manning | `precip_rate`, `roughness` |
| `OverlandFlow` | de Almeida et al. (2012) | `h_init`, `mannings_n` |
| `SoilMoisture` | Laio et al. (2001) bucket | `soil_porosity`, `soil_field_capacity` |
| `GroundwaterDupuitPercolator` | Dupuit approx. | `hydraulic_conductivity`, `porosity` |

### Tectonics & Lithology
| Component | Description | Key Params |
|-----------|-------------|------------|
| `NormalFault` | Vertical displacement on fault | `faulted_surface`, `fault_throw_rate_through_time` |
| `Flexure` | Lithospheric flexure (Airy/elastic) | `eet`, `youngs` |
| `Lithology` | Track rock layers through erosion | `layer_ids`, `layer_thicknesses` |

---

## I/O Formats

| Format | Read | Write | Function |
|--------|------|-------|----------|
| ESRI ASCII (.asc) | Yes | Yes | `landlab.io.esri_ascii.load()` / `dump()` |
| NetCDF (.nc) | Yes | Yes | `to_netcdf()` / `from_netcdf()` |
| Shapefile (.shp) | Yes | No | `read_shapefile()` |
| VTK Legacy (.vtk) | No | Yes | `write_legacy_vtk()` |
| OBJ (.obj) | No | Yes | `write_obj()` |
| Native Landlab | Yes | Yes | `save()` / `load()` |

---

## 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

---

## Synthetic test setups (how-to)

### SPACE Qs-Q test (quasi-1D channel)
- 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

### SPACE steady-state concavity test
- Grid: 50×50 RasterModelGrid, dx=200 m, open south boundary
- Duration: 4000 × 500 yr = 2 Myr; uplift U=1e-3 m/yr

### Whipple & Tucker (1999) steady-state test
- 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

## Lessons from past runs (rules)

1. Fit Qs vs Q directly, never SSC = Qs/Q vs Q: SSC vs Q always gives a negative exponent,
   whatever the model does (`dt_020`).
2. If the SPACE Qs-Q exponent b is outside [0.35, 1.20], check whether SPACE runs in the
   detachment-limited or transport-limited regime (H, H*, v_s) (`dt_021`).
3. For slope-area on a real DEM, use the 2D spatial gradient |∇z| with log-spaced area bins,
   not the D8 receiver slope (`topographic__steepest_slope`) (`dt_018`).
4. For slope-area on a DEM clip, use `open_all` boundaries; a single-outlet boundary
   (`closed_all_but_outlet`) pushes θ toward 0 (`dt_019`).
5. Do not use SPACE for annual sediment yield on a real DEM; use detachment-limited stream
   power, set K_sp analytically from the target SY, and use `closed_all_but_outlet` for
   watershed-scale SY (`dt_022`).
6. P export from TP = SY × 10 × P_soil_ppm × ER × 1e-6 is a gross potential; the amount at the
   watershed outlet is lower because of P retention inside the basin.

---

## Tools Reference

| Tool | Stage | Path | Lines | Purpose |
|------|-------|------|-------|---------|
| `convert_dem_to_grid.py` | s2 | `tools/convert_dem_to_grid.py` | ~250 | Load DEM → RasterModelGrid with fields |
| `convert_soil_params.py` | s3 | `tools/convert_soil_params.py` | ~220 | HWSD soil → grid fields (depth, K, porosity) |
| `run_landlab.py` | s5 | `tools/run_landlab.py` | ~280 | Execute Landlab simulation from YAML config |
| `parse_landlab_output.py` | s6 | `tools/parse_landlab_output.py` | ~250 | Extract grid fields → CSV/NetCDF + metrics |
| `dissect_atchafalaya_ssc_q_surrogate.py` | validation | `tools/dissect_atchafalaya_ssc_q_surrogate.py` | ~280 | SPACE Qs-Q rating and steady-state concavity check |
| `dissect_loess_plateau_slope_area.py` | validation | `tools/dissect_loess_plateau_slope_area.py` | ~430 | Real-DEM slope-area (spatial gradient, log-binned) |
| `dissect_loess_plateau_sediment_yield.py` | validation | `tools/dissect_loess_plateau_sediment_yield.py` | ~280 | DL stream power SY + particulate-P export |

---

## Calibration Parameters (Priority Order)

| Parameter | Component | Range | Controls | Sensitivity |
|-----------|-----------|-------|----------|-------------|
| K_sp | StreamPowerEroder | 1e-7–1e-3 | Channel erosion rate | Very High |
| linear_diffusivity | LinearDiffuser | 0.001–1.0 m²/yr | Hillslope curvature | High |
| m_sp | StreamPowerEroder | 0.3–0.7 | Area-discharge scaling | High |
| n_sp | StreamPowerEroder | 0.7–2.0 | Slope sensitivity | High |
| uplift_rate | (user loop) | 1e-5–1e-2 m/yr | Relief, channel gradient | Very High |
| v_s | Space/ErosionDeposition | 0.01–10 m/s | Transport vs detachment | High |
| K_sed / K_br | Space | 1e-7–1e-3 | Alluvium vs bedrock erosion | High |
| H_star | Space | 0.1–2.0 m | Alluvial cover effect | Medium |
| soil_production_maximum_rate | ExponentialWeatherer | 1e-5–1e-3 m/yr | Soil thickness | Medium |
| soil_production_decay_depth | ExponentialWeatherer | 0.2–2.0 m | Weathering decay | Medium |

---

## Data Requirements

| Data | Source | Format | Notes |
|------|--------|--------|-------|
| DEM | SRTM/ASTER/LiDAR | GeoTIFF or ESRI ASCII | Must be projected (meters) |
| Soil depth | HWSD/SoilGrids | Raster/CSV | Convert cm → m |
| Rainfall | CMFD/MSWX/ERA5 | NetCDF | Convert mm/hr → m/s for `water__unit_flux_in` |
| Uplift rate | Literature/GPS | Scalar or field | Must match dt units |
| Erodibility (K) | Literature/calibration | Scalar or field | Units depend on m, n |
| Diffusivity (D) | Literature/calibration | Scalar | m²/yr typical |

---

## Quick Start

```bash
# 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

| ID | Severity | Domain | Short Description |
|----|----------|--------|-------------------|
| dt_001 | silent | unit_conversion | Rainfall in mm/hr instead of m/s |
| dt_002 | silent | unit_conversion | Diffusivity time-unit mismatch |
| dt_003 | silent | unit_conversion | K_sp units wrong for given m,n |
| dt_004 | silent | unit_conversion | Grid spacing in degrees not meters |
| dt_005 | silent | unit_conversion | Timestep unit mismatch with parameters |
| dt_006 | silent | unit_conversion | Soil depth in cm instead of m |
| dt_007 | silent | unit_conversion | Elevation in feet instead of meters |
| dt_008 | silent | unit_conversion | Rainfall mm vs m for soil moisture |
| dt_009 | silent | parameter_format | Manning n vs Chezy C confusion |
| dt_010 | degraded | input_format | DEM nodata not handled |
| dt_011 | fatal | boundary_condition | All boundaries closed |
| dt_012 | degraded | missing_component | No depression routing |
| dt_013 | degraded | numerical_stability | CFL violation in diffusion |
| dt_014 | fatal | component_order | Erosion before flow routing |
| dt_015 | silent | field_location | Field at wrong grid element |
| dt_016 | silent | ssc_q_rating_curve | Wrong runoff_rate breaks SPACE SSC-Q correlation |
| dt_017 | silent | ssc_q_rating_curve | m_sp/n_sp wrong → concavity outside [0.40, 0.55] |
| dt_018 | silent | slope_area_analysis | D8 receiver slope polluted by depression-filling → use spatial gradient |
| dt_019 | silent | slope_area_analysis | Single-outlet boundary depresses θ toward 0 → use open_all for DEM clips |
| dt_020 | silent | ssc_q_rating_curve | Fitting SSC vs Q gives negative exponent — always fit Qs vs Q |
| dt_021 | silent | ssc_q_rating_curve | SPACE b outside [0.35,1.20] — check DL vs TL regime (H, H*, v_s) |
| dt_022 | silent | model_scope_mismatch | SPACE near-zero SY on real DEMs — use DL stream power for annual SY |

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                 # diagnostic triplets (count: KI map)
```