PIHM · diff

git:20260817.960ae14 to git:20260817.6859167

11 added, 0 removed. Audit C to A.

+ ---

+ name: pihm

+ description: >-

+ MM-PIHM v1.0.0 — semi-discrete finite-volume coupled surface-subsurface watershed

+ formulation. Covers Spatially-distributed physically-based watershed hydrology on an

+ unstructured triangular (TIN) mesh; Overland (surface) flow per element; Unsaturated

+ (vadose) zone soil water and infiltration; Saturated (groundwater) zone lateral and

+ vertical flow, recharge. Use when the task involves running, configuring, calibrating or

+ interpreting PIHM.

+ ---

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

# MM-PIHM Knowledge Infrastructure

**Package**: `pihm-ki` | **Version**: 1.0.0 | **Model**: MM-PIHM v1.0.0

**Domain**: Distributed Hydrology | **Language**: C (gcc) | **Solver**: SUNDIALS CVODE v7.3.0

**License**: MIT | **Repository**: PSUmodeling/MM-PIHM

| Metric | Value |

|--------|-------|

| Pipeline stages | 9 |

| Python tools | 4 |

| Skill documents | 5 |

| Diagnostic triplets | 20 |

| Failure domains | 5 |

---

## 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 CMFD unit documentation and known traps.

See `data_ki/HWSD/SKILL.md` for soil property documentation.

See `data_ki/ObservedQ/SKILL.md` for observed discharge data.

## 1. Overview

MM-PIHM (Multi-Modular Penn State Integrated Hydrologic Model) is a physically-based,

spatially-distributed watershed model that solves coupled surface–subsurface hydrology

on unstructured triangular meshes using finite-volume methods and the CVODE implicit

ODE solver from SUNDIALS.

**Model Variants:**

- **PIHM** — core hydrology (surface water, unsaturated zone, groundwater, channel flow)

- **Flux-PIHM** — adds Noah Land Surface Model for energy balance, snow, soil temperature

- **Flux-PIHM-BGC** — adds Biome-BGC for carbon/nitrogen cycling

- **Cycles-L** — adds Cycles agroecosystem model for crop growth

**Key Capabilities:**

- Triangular Irregular Network (TIN) mesh for irregular terrain

- Coupled overland flow, vadose zone, saturated zone, and channel routing

- van Genuchten soil hydraulics with macropore flow

- Adaptive implicit time stepping via CVODE

- OpenMP parallelization for large domains

- Spin-up mode for steady-state initialization

**State Variables (ODE system):**

- Surface water depth per element (m)

- Unsaturated zone storage per element (m)

- Groundwater head per element (m)

- River stage per segment (m)

- Total ODEs: `3 * nelem + nriver`

---

## 2. Installation

### 2.1 Dependencies

- GCC with OpenMP support

- CMake >= 3.18 (auto-downloaded if missing)

- SUNDIALS CVODE v7.3.0 (bundled in `cvode/` directory)

### 2.2 Build Steps

```bash

cd MM-PIHM/

make cvode # Install CVODE library

make pihm # Compile PIHM (core hydrology)

# OR

make flux-pihm # Compile Flux-PIHM (with Noah LSM)

# OR

make flux-pihm-bgc # Compile Flux-PIHM-BGC (with BGC)

```

### 2.3 OpenMP Configuration

```bash

export OMP_NUM_THREADS=12

```

### 2.4 Test Run

```bash

make test # Compiles Flux-PIHM and runs ShaleHills example

```

---

## 3. Pipeline Stages

|---|-------|-------------|-------|--------|------|

**Stage Dependencies:**

- S1 → S2 → S3,S4,S5,S6 (parallel) → S7 → S8 → S9 → S7 (iterate)

---

## ✅ NEW BASINS ARE SUPPORTED HEADLESSLY (since 2026-08-02)

Stages **S1 (TIN mesh)** and **S2 (attribute assignment)** used to require

**PIHMgis**, an interactive QGIS-plugin GUI that is not installed and is not a

batch executable (so the WINE pattern does NOT apply). That block is CLOSED:

tools/mesh_builder.py

generates `<project>.mesh` / `.att` / `.riv` for an arbitrary watershed from a

DEM + outlet coordinate, using

`ki_tools_common.terrain_ops.delineate_basin` (whitebox breach → D8 →

accumulation → snap → watershed) +

`extract_river_centerline_from_streams` for the main stem, then a

**constrained-Delaunay TIN** (meshpy/Triangle) in which the basin boundary and

the river centreline are CONSTRAINED segments — so every river reach lands on

element edges, which is what `InitRiver` requires.

```bash

python tools/mesh_builder.py --dem KISSPATH_DATA/MERIT_DEM/n40e005_dem.tif \

--outlet-lon 8.6124 --outlet-lat 42.1771 \

--out-dir input/ChiuniFR --project ChiuniFR \

--target-elem-area-km2 0.12 --aquifer-thickness-m 15

```

Verified 2026-08-02: 49.5 km² Chiuni catchment (Corsica, FR) → 752 elements,

446 nodes, 14 river segments; `flux-pihm` runs it end to end.

**Conventions the generator encodes** (all verified against `input/ShaleHills`):

| Rule | Why | Failure if broken |

|---|---|---|

| Elements must be **counter-clockwise** | `InitTopo` uses a SIGNED area `0.5*((x1-x0)(y2-y0)-(y1-y0)(x2-x0))` | negative element area → negative storage, silent |

| `NABR[j]` = neighbour across the edge **OPPOSITE node j**, 0 = none | verified 1544/1544 ShaleHills records | fluxes routed to the wrong neighbour |

| River `LEFT`/`RIGHT` must be **mesh neighbours of each other** | `InitRiver` looks up `elem[left].nabr[j] == right` | river never couples to the hillslope |

| Node elevations must **decrease** along the river | outlet `DOWN = -3` (ZERO_DPTH_GRAD) computes `sqrt(grad_h)` | NaN discharge → CVODE dies on step 1 |

| Drainage area for channel width must be a **running maximum** | a segment midpoint can miss the 1-cell stream raster | 1 m wide channel mid-mainstem |

**Still true (dt_023):** a headless TIN gives a SINGLE-THREAD main stem, so

**basin-outlet DISCHARGE on a large basin is still not trustworthy** — tributary

inflow is not accumulated. Use the headless path for small/medium basins and for

**state** variables (GW head, soil moisture, ET), and validate discharge only

where the delineated main stem carries the basin.

**S5/S6** (`.para` / `.calib` / `.lsm` / `.lai`, spin-up IC) are likewise no

longer "Manual": see `tools/config_writer.py`.

## 4. Input Files Reference

All input files share a project prefix: `<project>.<ext>`.

For project "ShaleHills", files are in `input/ShaleHills/`.

### 4.1 Mesh File (`.mesh`)

```

NUMELE 535

INDEX NODE1 NODE2 NODE3 NABR1 NABR2 NABR3

1 1 2 3 0 2 5

...

NUMNODE 293

INDEX X(m) Y(m) ZMIN(m) ZMAX(m)

1 634595.29 4520225.42 225.13 228.08

```

### 4.2 Attribute File (`.att`)

```

INDEX SOIL GEOL LC METEO LAI BC1 BC2 BC3

1 1 1 28 1 0 0 0 0

```

- **There is NO `NUMATT` record.** `src/read_att.c` CheckHeaders the FIRST line

against `INDEX SOIL GEOL LC METEO LAI BC1 BC2 BC3` and takes the row count

from `nelem` in the `.mesh`. Writing a leading `NUMATT <n>` aborts with

"input file format error ... at Line 1" (dt_024).

- BC values: 0 = no flow, >0 = boundary condition series index

- `GEOL 0` is valid (no `.geol` file); the soil column is then used throughout.

### 4.3 Soil File (`.soil`)

`src/read_soil.c` calls `CheckHeader(16, INDEX SILT CLAY OM BD KINF KSATV KSATH

MAXSMC MINSMC ALPHA BETA MACHF MACVF DMAC QTZ)` on the line after `NUMSOIL` —

the column header is MANDATORY (dt_024).

```

NUMSOIL 5

INDEX SILT(%) CLAY(%) OM(%) BD(g/cm3) KINF(m/s) KSATV(m/s) KSATH(m/s)

MAXSMC(m3/m3) MINSMC(m3/m3) ALPHA(1/m) BETA(-) MACHF(-) MACVF(-)

DMAC(m) QTZ(-)

DINF 0.10

KMACV_RO 100.0

KMACH_RO 1000.0

```

### 4.4 Meteorological Forcing (`.meteo`)

```

METEO_TS 1 WIND_LVL 10.0

TIME PRCP SFCTMP RH SFCSPD SOLAR LONGWV PRES

#TS kg/m2/s K % m/s W/m2 W/m2 Pa

2009-01-01 00:00 0.0 271.55 65.0 2.3 0.0 280.5 101325.0

```

- **`METEO_TS` and `WIND_LVL` are ONE line.** `src/read_forc.c` does

`sscanf(cmdstr, "%s %d %s %lf", ...)` and requires `match == 4`, then

`CheckHeader(8, TIME PRCP SFCTMP RH SFCSPD SOLAR LONGWV PRES)`. Two separate

lines, or a missing column header, is ERR_WRONG_FORMAT (dt_024).

- `WIND_LVL` must match the forcing's wind measurement height: **2.0 for NASA

POWER** (`ki_tools_common` requests `WS2M`), 10.0 for MSWX/CMFD.

### 4.5 River File (`.riv`)

```

NUMRIV 20

INDEX FROM TO DOWN LEFT RIGHT SHAPE MATL BC RES

SHAPE section: INDEX DPTH(m) OINT CWID

MATERIAL section: INDEX ROUGH(s/m^1/3) CWR(-) KH(m/s)

```

### 4.6 Parameter File (`.para`)

Controls simulation mode, time window, time steps, solver tolerances, and output intervals.

- `MODEL_STEPSIZE` — hydrology step (seconds, typically 60)

- `ABSTOL` — CVode absolute tolerance (m)

- `RELTOL` — CVode relative tolerance (dimensionless)

### 4.7 Calibration File (`.calib`)

Global multipliers applied to soil/river/vegetation properties:

- Hydraulic conductivity: `KSATH`, `KSATV`, `KINF` (multiplier, default 1.0)

- Macropore: `KMACSATH`, `KMACSATV`, `DMAC`, `MACVF`, `MACHF`

- Vegetation: `VEGFRAC`, `ALBEDO`, `ROUGH`, `DROOT`

- River: `ROUGH_RIV`, `KRIVH`, `RIV_DPTH`, `RIV_WDTH`

- Forcing: `PRCP` (multiplier), `SFCTMP` (offset in K)

### 4.8 Initial Conditions (`.ic`)

Per-element values for: surface water (m), unsaturated storage (m), groundwater (m), snow (m), canopy storage (m), soil temperature (K).

Per-river: stage (m).

### 4.9 LAI Forcing (`.lai`)

```

LAI_TS 1

TIME LAI(m2/m2)

2009-01-01 0.5

2009-04-01 3.2

```

### 4.10 Boundary Conditions (`.bc`)

```

BC_TS 1

TYPE 1 (1=Dirichlet head, 2=Neumann flux)

TIME VALUE(m or m3/s)

```

---

## 5. Unit Trap Table

|----------|--------------------|---------------------|------------|---------|

| Temperature | K | °C | +273.15 | dt_002 |

| Pressure | Pa | hPa/mbar | ×100 | dt_003 |

| Wind speed | m/s | km/hr | ÷ 3.6 | dt_006 |

| van Genuchten alpha | 1/m | 1/cm | ×100 | dt_008 |

| Soil depth/elevation | m | cm or ft | ÷ 100 or ÷ 3.281 | dt_009 |

| Bulk density | g/cm3 | kg/m3 | ÷ 1000 | dt_010 |

---

## 6. Tools Reference

|---|------|-------|---------|-------|

| 0 | `mesh_builder.py` | S1+S2 | **Headless** DEM+outlet → `.mesh`/`.att`/`.riv` | ~470 |

| 0b | `config_writer.py` | S5+S6 | `.para`/`.calib`/`.lsm`/`.lai` (+ spin-up modes) | ~180 |

| 1 | `forcing_converter.py` | S3 | Convert ERA5/CSV met data → PIHM `.meteo` format | ~280 |

| 2 | `soil_converter.py` | S4 | Convert HWSD/SSURGO → PIHM `.soil` and `.geol` | ~250 |

| 3 | `run_pihm.py` | S7 | Execute PIHM binary with validation | ~200 |

| 4 | `output_parser.py` | S8 | Parse binary output → CSV with headers | ~260 |

---

## 7. Critical Domain Knowledge

### 7.1 Precipitation units are kg/m2/s, NOT mm/hr (dt_001)

PIHM expects precipitation as mass flux (kg/m2/s). Since water density = 1000 kg/m3,

1 kg/m2/s = 1 mm/s = 3600 mm/hr. Passing mm/hr directly without dividing by 3600

produces rainfall 3600× too high, causing immediate flooding and solver crash.

### 7.2 Temperature must be in Kelvin (dt_002)

All temperature fields (SFCTMP, TBOT) use Kelvin. Passing Celsius creates ~273 K offset

that makes snow/ice calculations fail silently — winter snow accumulates all year.

### 7.3 Pressure in Pascals, not hectopascals (dt_003)

Surface pressure uses Pa (e.g., 101325). Passing hPa (1013.25) makes pressure 100× too

low, causing incorrect vapor pressure calculations and unrealistic ET.

### 7.4 van Genuchten alpha is in 1/m, not 1/cm (dt_008)

HWSD/Rosetta typically provide alpha in 1/cm. PIHM needs 1/m. Forgetting the ×100

conversion makes soils appear 100× more permeable, draining unrealistically fast.

### 7.5 Calibration file uses multipliers, not absolute values (dt_011)

KSATH=10 means multiply base Ksat by 10, not set Ksat to 10 m/s. Setting absolute

values (e.g., 1e-5 for Ksat) when the field is a multiplier produces near-zero flow.

### 7.6 Spin-up is essential for groundwater initialization (dt_012)

Running without spin-up (INIT_MODE=0, SIMULATION_MODE=0) starts from arbitrary IC.

Groundwater may take years to equilibrate, producing unrealistic baseflow for the

first several years. Use SIMULATION_MODE=1 for spin-up, then restart with the `.ic` output.

### 7.7 CVODE tolerance controls mass conservation (dt_013)

Setting ABSTOL too loose (>0.01 m) causes mass balance errors that accumulate.

Too tight (<1e-6) causes solver failures. Default ABSTOL=1e-4 m is usually safe.

### 7.8 Output intervals use special codes (dt_014)

Para file output codes: -1=yearly, -2=monthly, -3=daily, -4=hourly, 0=off, or

positive integer for seconds. Confusing -3 (daily) with 3 (every 3 seconds) produces

massive output files that fill disk.

---

## 8. Output Variables

### 8.1 Element-based outputs (per triangular element)

|----------|-------------|------|-------------|

| GW | `.gw.dat` | m | Groundwater head **above the element's aquifer bottom (`zmin`)**, NOT depth-to-water. Water-table elevation = `elem.topo.zmin + GW`; depth below ground = `zmax - zmin - GW` (`lat_flow.c`: `elem.ws.gw + elem.topo.zmin`) |

| CMC | `.cmc.dat` | m | Canopy moisture storage |

| EC | `.ec.dat` | m/s | Canopy evaporation |

| ETT | `.ett.dat` | m/s | Transpiration |

### 8.2 River-based outputs (per river segment)

|----------|-------------|------|-------------|

### 8.3 Output File Format

Binary files contain double-precision values. Each record:

`[time_stamp (int)] [value_elem1 (double)] [value_elem2 (double)] ... [value_elemN (double)]`

Use `output_parser.py` or the `PIHM-utils` Python package to convert to CSV/DataFrame.

---

## 9. Execution

### 9.1 Command Line

```bash

./pihm [-b] [-c] [-d] [-f] [-s] [-t] [-V] [-v] [-o output_dir] <project_name>

```

| Flag | Description |

|------|-------------|

| `-b` | Brief mode (minimal screen output) |

| `-c` | Elevation correction (fix surface sinks) |

| `-d` | Debug mode (CVODE log file) |

| `-f` | Fixed-length spin-up (don't stop at equilibrium) |

| `-s` | Silent mode |

| `-t` | Tecplot output format |

| `-v` | Verbose mode |

| `-V` | Print version and exit |

| `-o dir` | Custom output directory name |

### 9.2 Directory Layout

```

MM-PIHM/

input/<project>/ # All input files

<project>.mesh

<project>.att

<project>.soil

<project>.geol

<project>.lc

<project>.meteo

<project>.lai

<project>.bc

<project>.para

<project>.calib

<project>.ic

<project>.riv

output/<run_name>/ # Output directory

<project>.surf.dat

<project>.gw.dat

...

```

### 9.3 Spin-up Workflow

```bash

# Step 1: Run spin-up

# Set SIMULATION_MODE = 1 in .para

./pihm -o spinup_run ShaleHills

# Step 2: Copy IC file

cp output/spinup_run/ShaleHills.ic input/ShaleHills/ShaleHills.ic

# Step 3: Run with restart

# Set SIMULATION_MODE = 0, INIT_MODE = 1 in .para

./pihm -o production_run ShaleHills

```

---

## 10. Calibration Guide

### 10.1 Priority Parameters (sensitivity order)

|----------|-----------|-------|----------|

| 1 | KSATH | 0.01–100 | Lateral groundwater flow, baseflow |

| 2 | KSATV | 0.01–100 | Vertical percolation, recharge |

| 3 | KINF | 0.01–100 | Infiltration rate, surface runoff |

| 4 | KMACSATH | 0.01–100 | Macropore lateral flow |

| 5 | DMAC | 0.1–10 | Macropore depth, preferential flow |

| 6 | ROUGH_RIV | 0.1–10 | Channel flow velocity |

| 7 | PRCP | 0.8–1.2 | Precipitation bias correction |

| 8 | ALPHA | 0.1–10 | Soil water retention curve |

| 9 | POROSITY | 0.5–2.0 | Soil storage capacity |

| 10 | DROOT | 0.1–5.0 | Root zone depth, ET partitioning |

### 10.2 Calibration Strategy

1. Start with spin-up to get realistic IC

2. Adjust KSATH for baseflow recession shape

3. Adjust KINF/KSATV for storm peak magnitude

4. Adjust DMAC/KMACSATH for quick-flow response

5. Adjust PRCP if systematic volume bias exists

6. Fine-tune ROUGH_RIV for peak timing

---

## 11. Validation Reference

**Test Basin:** Shale Hills, Pennsylvania, USA

**Coordinates:** 40.66°N, 77.91°W

**Area:** ~8 ha (0.08 km2)

**Period:** 2009-01-01 to 2009-12-31

**Elements:** 535 triangular elements, 20 river segments

**Data Source:** Bundled ShaleHills example

---

## 12. Quick Start

```bash

# 1. Clone and build

git clone https://github.com/PSUmodeling/MM-PIHM.git

cd MM-PIHM

make cvode

make pihm

# 2. Set OpenMP threads

export OMP_NUM_THREADS=4

# 3. Run ShaleHills example

./pihm -o test_run ShaleHills

# 4. Parse output

python3 ki/tools/output_parser.py \

--input output/test_run/ \

--project ShaleHills \

--variables gw,surf,rivflx1 \

--output results.csv

```

---

## 13. Diagnostic Triplets Summary

|----|-------|---------------|----------|-------|

| dt_022 | — | domain_validity | n/a | Pumped-aquifer water-table obs are out of PIHM's domain |

---

## 14. File Structure

```

ki/

SKILL.md # This file

tools/

forcing_converter.py # S3: Met data → PIHM .meteo

soil_converter.py # S4: Soil DB → PIHM .soil/.geol

run_pihm.py # S7: Execute PIHM binary

output_parser.py # S8: Binary output → CSV

docs/

s3_forcing_preparation.md # Forcing preparation skill

s4_soil_parameterization.md # Soil parameter skill

s5_configuration.md # Parameter configuration skill

s7_model_execution.md # Model execution skill

s8_output_analysis.md # Output analysis skill

diagnostics/

triplets.yaml # 20 diagnostic triplets

```