Landlab git:20260414.2b184d0C

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

> **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.

# 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-03-25
**Stats**: 4 tools | 5 skill documents | 15+ diagnostic triplets | ~3,400 lines of validated code
**Validation status**: production_validated (Whipple & Tucker 1999 steady-state test)

---

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


## 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()
```

---

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

---

## 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
| Metric | Simulated | Expected | Status |
|--------|-----------|----------|--------|
| Slope-area r² | 0.9978 | ~1.0 | PASS |
| Concavity (θ) | 0.488 | 0.500 | PASS (2.4% error) |
| Relief (m) | 810.2 | ~800 | PASS |
| Wall time | 2.2 s | — | — |

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

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

---

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

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
├── 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                 # 15 symptom→diagnosis→remedy
```