Docs

# Noah-MP Runoff Options — Skill Document

> **Stage ID**: s9_runoff_options
> **Pipeline order**: 9 of 12 (configured in namelist.hrldas, Stage s9_namelists)
> **Depends on**: s5_soil_properties, s9_namelists

## Purpose

Noah-MP's `RUNOFF_OPTION` (also called `opt_run` in the source code) controls how precipitation is partitioned into surface runoff, infiltration, and subsurface drainage. This single integer switch selects one of five fundamentally different runoff generation physics, each suited to different basin types. Choosing the wrong option can produce 24x less surface runoff than observations (dt_v018), crash the model (dt_v019), or invert the seasonal discharge cycle (dt_v017). This document covers all five options, their physics, parameter sensitivities, and calibration strategies.

## Prerequisites

Before selecting a runoff option, verify:

- [ ] `soil_properties.nc` exists in DOMAIN/ with per-cell soil hydraulic parameters (SMCMAX, DKSAT, BEXP)
- [ ] `namelist.hrldas` exists and is editable (generated by `generate_namelists.py`)
- [ ] Basin characteristics are known: climate zone, terrain type, typical soil depth, dominant runoff mechanism
- [ ] For Option 5 only: `MMF_RUNOFF_FILE` is available (not generated by the standard pipeline)

## Inputs

| Input | Type | Source | Description |
|-------|------|--------|-------------|
| `RUNOFF_OPTION` | integer (1-5) | namelist.hrldas | Selects the runoff generation scheme |
| `REFKDT` | float | soil_properties.nc or GENPARM.TBL | Infiltration scaling factor (Options 1, 3) |
| `BEXP` | float array | soil_properties.nc | Clapp-Hornberger b-exponent per soil layer |
| `DKSAT` | float array | soil_properties.nc | Saturated hydraulic conductivity per soil layer |
| `SMCMAX` | float array | soil_properties.nc | Porosity per soil layer |
| `SMCWLT` | float array | soil_properties.nc | Wilting point per soil layer |
| `SLOPE` | float | GENPARM.TBL | Slope index for subsurface drainage (0-1) |
| `MMF_RUNOFF_FILE` | file path | namelist.hrldas | Required ONLY for Option 5 |

## Procedure

### Step 1: Assess Basin Characteristics

Determine the dominant runoff mechanism for your basin:

| Basin Type | Dominant Mechanism | Recommended Option |
|------------|-------------------|-------------------|
| Humid mountainous, thin soils | Saturation excess | 1 (SIMGM/TOPMODEL) |
| General purpose, mixed terrain | Simplified saturation excess | 2 (SIMTOP) |
| Flat terrain, deep soils, alluvial | Infiltration excess | 3 (Schaake96) — **default** |
| Simple/quick applications | Bucket drainage | 4 (BATS) |
| VIC-like behavior, variable capacity | Variable infiltration | 5 (Xinanjiang/VIC) |

### Step 2: Understand Each Option's Physics

---

#### Option 1: SIMGM (TOPMODEL-based with groundwater)

**How it works**: Combines a TOPMODEL-style saturation-excess surface runoff scheme with a simple unconfined aquifer groundwater model. The topographic wetness index determines which fraction of the basin is saturated. Water table depth controls the saturated fraction: as the water table rises (wetter conditions), more area becomes saturated and generates runoff from any additional rainfall.

**Key formula (plain language)**:
- Compute water table depth from soil moisture profile
- Saturated fraction = f(topographic index, water table depth)
- Surface runoff = precipitation x saturated fraction
- Subsurface drainage = f(water table depth, SLOPE)

**Key parameters and their effects**:
| Parameter | Effect | Typical Range |
|-----------|--------|---------------|
| REFKDT | Controls infiltration rate vs surface runoff split | 0.1-10.0 |
| SLOPE | Controls subsurface lateral drainage speed | 0.01-1.0 |
| BEXP | Soil moisture retention curve shape | soil-dependent |

**When to use**:
- Mountain basins with variable topography and thin soils
- Basins where saturation-excess runoff dominates (humid, high water table)
- When you need physically-based groundwater interaction
- **Best for**: Humid mountainous basins (Chaohe, Tibetan Plateau tributaries)

**Known issues**:
- Requires no extra files beyond standard pipeline output
- More computationally expensive than Options 2-4 (groundwater table tracking)
- Can produce too much baseflow in deep-soil basins if SLOPE is too high

**Calibration strategy**:
1. Start with REFKDT=1.0 (moderate infiltration)
2. Tune SLOPE for baseflow recession match (lower = slower recession)
3. Fine-tune REFKDT for peak magnitude (lower REFKDT = more surface runoff)

---

#### Option 2: SIMTOP (Simplified TOPMODEL)

**How it works**: A simplified version of Option 1 without the explicit groundwater model. Uses a statistical representation of the topographic index distribution to compute the saturated area fraction. Subsurface drainage is controlled by the SLOPE parameter and soil hydraulic properties.

**Key formula (plain language)**:
- Saturated fraction = exponential function of mean soil moisture deficit
- Surface runoff = precipitation x saturated fraction (saturation excess)
- Subsurface drainage = DKSAT x SLOPE x f(soil moisture)

**Key parameters and their effects**:
| Parameter | Effect | Typical Range |
|-----------|--------|---------------|
| SLOPE | Dominant control on drainage rate | 0.01-1.0 |
| BEXP | Soil water retention, affects saturation fraction | soil-dependent |
| DKSAT | Base drainage rate | soil-dependent |

**When to use**:
- General-purpose applications where TOPMODEL behavior is desired without groundwater complexity
- Medium-sized basins with moderate terrain
- When Option 1 is too slow but you want saturation-excess physics

**Known issues**:
- Less physical than Option 1 (no explicit water table)
- Sensitive to SLOPE parameter — too high causes excessive drainage
- May underestimate runoff in very steep terrain

**Calibration strategy**:
1. SLOPE is the primary knob (0.1 for slow drainage, 0.5 for moderate, 1.0 for fast)
2. BEXP multiplier for fine-tuning seasonal soil moisture dynamics
3. Check that total ET + Q approximately equals P (water balance closure)

---

#### Option 3: Schaake96 (Infiltration Excess) — DEFAULT

**How it works**: Uses the Schaake (1996) infiltration partitioning formula derived from Philip's equation. Computes maximum infiltration capacity as a function of soil moisture deficit and saturated hydraulic conductivity, then any rainfall exceeding this capacity becomes surface runoff.

**Key formula (plain language)**:
```
KDT = REFKDT * DKSAT(layer1) / REFDK
VAL = 1 - exp(-KDT * dt / 86400)
DD  = sum(layer_depth * (SMCMAX - SMCWLT) * (1 - (SH2O - SMCWLT) / (SMCMAX - SMCWLT)))
INFMAX = (P * DD * VAL) / (P + DD * VAL) / dt
RUNSRF = max(0, QINSUR - INFMAX)
```

Where:
- `REFKDT` scales infiltration capacity (higher = more infiltration = less runoff)
- `DD` is the soil moisture deficit (mm of available storage)
- `REFDK` is a reference conductivity constant (2.0e-6 m/s)
- `QINSUR` is the total water input rate (rain + snowmelt)

**Key parameters and their effects**:
| Parameter | Effect | Typical Range |
|-----------|--------|---------------|
| **REFKDT** | **#1 most important**: controls 80% of peak magnitude | 0.1-10.0 (default 3.0) |
| DKSAT | Base infiltration rate (soil-dependent) | x0.1 to x10 multiplier |
| BEXP | Affects DD computation (soil moisture deficit) | soil-dependent |
| SLOPE | Subsurface drainage rate | 0.01-1.0 |

**When to use**:
- Flat terrain with deep soils (alluvial plains, lowlands)
- Basins where Hortonian (infiltration-excess) overland flow dominates
- Default choice for WRF-Hydro — most extensively tested
- **Good for**: Bengbu/Huai River (r=0.84 vs VIC at 0.25deg)

**Known issues** (verified in Chaohe investigation):
- **REFKDT=3.0 (default) is too high for mountain basins** — absorbs nearly all rainfall into soil. Chaohe test: 0.158mm surface runoff from 23.8mm rain, while VIC produced 3.78mm (24x more). See dt_v018.
- Sensitivity: REFKDT=3.0 gives flat baseflow; REFKDT=0.1 gives 5x too-high peaks (1750 vs 340 m3/s). Optimal for Chaohe: 0.3-0.5 with routing attenuation.
- Does NOT account for saturation from below (no TOPMODEL concept)
- In dry soils (large DD), almost all rain infiltrates regardless of intensity

**Calibration strategy** (from calibration_guide.md):
1. **REFKDT first** — start at 3.0, reduce to 0.5-1.0 for mountain basins, increase to 5-10 for swamp/paddy
2. **MannN + OVROUGHRTFAC** — attenuate remaining peak overshoot after REFKDT adjustment
3. **SLOPE + GW params** — match baseflow recession curve
4. **CZIL** — adjust total water balance (ET) if runoff volume is correct but ET is wrong

---

#### Option 4: BATS (Bucket-style)

**How it works**: Uses the BATS (Biosphere-Atmosphere Transfer Scheme) approach — a simple bucket model where surface runoff is generated when soil moisture exceeds field capacity. Drainage is proportional to the excess above field capacity.

**Key formula (plain language)**:
- If soil moisture > field capacity: runoff = excess
- Drainage = gravity drainage function of soil moisture
- Simple, computationally cheap

**Key parameters and their effects**:
| Parameter | Effect | Typical Range |
|-----------|--------|---------------|
| SMCMAX | Porosity (bucket size) | soil-dependent |
| SMCREF | Field capacity (overflow threshold) | soil-dependent |
| DKSAT | Gravity drainage rate | soil-dependent |

**When to use**:
- Quick feasibility studies
- Very simple applications where detailed runoff physics is not important
- Teaching and demonstration runs
- Basins where you have no calibration data and want a robust default

**Known issues**:
- No infiltration-excess mechanism — cannot produce Hortonian overland flow
- Tends to underestimate flood peaks (all rain infiltrates until bucket is full)
- Less physically realistic than Options 1-3

**Calibration strategy**:
1. SMCMAX controls bucket size — larger = more storage before runoff
2. DKSAT controls drainage speed — affects recession limb
3. Limited calibration potential — consider switching to Option 1 or 3 for better results

---

#### Option 5: Xinanjiang/VIC (Variable Infiltration Capacity)

**How it works**: Implements the Xinanjiang (or VIC) variable infiltration capacity curve. The basin is conceptualized as having a distribution of infiltration capacities (from zero at fully saturated areas to maximum at the driest areas). The fraction of the basin that generates runoff varies continuously with soil moisture state.

**Key formula (plain language)**:
- Infiltration capacity varies spatially following a power-law distribution
- Surface runoff = integral of the area where rainfall > local infiltration capacity
- Similar physics to VIC model's ARNO baseflow scheme

**Key parameters and their effects**:
| Parameter | Effect | Typical Range |
|-----------|--------|---------------|
| b_parameter (from MMF file) | Shape of infiltration capacity distribution | 0.01-1.0 |
| Topographic index fields | Spatial distribution of infiltration | Derived from DEM |

**When to use**:
- When you want WRF-Hydro results comparable to VIC (same runoff generation physics)
- Basins where both saturation-excess and infiltration-excess contribute
- Chinese basins where the Xinanjiang model has been widely calibrated

**Known issues** (verified — dt_v019):
- **CRASHES without MMF_RUNOFF_FILE** — error: `"Problem opening wrfinput file for read_mmf_runoff"`. The standard HydroCraft pipeline does NOT generate this file. You must either:
  1. Use NCAR's preprocessing tools to create the MMF file from topographic data
  2. Use Option 1 (TOPMODEL) instead, which provides similar saturation-excess physics without the extra file
- MMF_RUNOFF_FILE must be specified in namelist.hrldas
- Less tested in WRF-Hydro than Options 1 and 3

**Calibration strategy**:
1. First, generate the MMF_RUNOFF_FILE (beyond current pipeline scope)
2. The b_parameter in the MMF file controls runoff sensitivity — higher b = more variable infiltration = more runoff
3. Compare directly with VIC output for the same basin to validate

---

### Step 3: Set the Runoff Option

Edit `namelist.hrldas` in the `&NOAHLSM_OFFLINE` section:

```fortran
RUNOFF_OPTION = 3     ! Change to 1, 2, 3, 4, or 5
```

For Option 5 only, also add:
```fortran
MMF_RUNOFF_FILE = "./DOMAIN/MMF_runoff.nc"
```

**Command**:
```bash
# To change RUNOFF_OPTION in an existing namelist:
sed -i 's/RUNOFF_OPTION.*=.*/RUNOFF_OPTION = 1/' namelist.hrldas
```

Or regenerate namelists via `generate_namelists.py` (which defaults to Option 3).

**Expected result**: namelist.hrldas contains the desired `RUNOFF_OPTION` value.

**If this fails**: For Option 5, see dt_v019. Switch to Option 1.

### Step 4: Adjust REFKDT for Options 1 and 3

REFKDT is the single most important calibration parameter for Schaake96 (Option 3) and also affects TOPMODEL (Option 1).

```python
import netCDF4 as nc

ds = nc.Dataset("DOMAIN/soil_properties.nc", "a")
refkdt = ds["refkdt"][:]
refkdt[:] = 0.5  # Reduce from default 3.0 for mountain basins
ds["refkdt"][:] = refkdt
ds.close()
```

| Basin type | Recommended REFKDT |
|------------|-------------------|
| Mountain, thin soils | 0.3-0.5 |
| General humid | 1.0-2.0 |
| Flat, deep soils (default) | 3.0 |
| Paddy/wetland | 5.0-10.0 |

### Step 5: Verify Runoff Generation After First Run

After a short test run (7-14 days), check surface runoff magnitude:

```bash
# Check SFCRNOFF in LDASOUT (should be non-zero for rain events)
python -c "
import netCDF4 as nc, numpy as np
ds = nc.Dataset('2001010200.LDASOUT_DOMAIN1')
sfcrnoff = ds['SFCRNOFF'][0]
print(f'SFCRNOFF: mean={np.nanmean(sfcrnoff):.3f}, max={np.nanmax(sfcrnoff):.3f} mm')
ds.close()
"
```

If mean SFCRNOFF is near zero during rain events, the runoff option or REFKDT needs adjustment.

## Expected Outputs

After successful completion, the following should exist:

| Output | Path | Verification |
|--------|------|--------------|
| Updated namelist.hrldas | `<run_dir>/namelist.hrldas` | Contains `RUNOFF_OPTION = <N>` |
| LDASOUT with non-zero SFCRNOFF | `<run_dir>/*.LDASOUT_DOMAIN1` | `SFCRNOFF > 0` during rain events |
| Reasonable runoff ratio | Computed from output | Q/P between 0.1 and 0.8 for most basins |

## Validation Checks

Run these checks before proceeding to output analysis:

1. **Surface runoff magnitude check**: Verify SFCRNOFF is physically reasonable
   - Command: `python -c "import netCDF4 as nc; ds=nc.Dataset('<latest_LDASOUT>'); print(ds['SFCRNOFF'][0].mean())"`
   - Expected: 0.01-10 mm per output timestep during rain events
   - If near zero: REFKDT too high (dt_v018) or wrong runoff option for basin type
   - If unexpected: See diagnostic triplet dt_v018

2. **Runoff ratio check**: Compare total runoff to total precipitation
   - Command: Sum SFCRNOFF + UGDRNOFF over entire run, divide by total RAINRATE
   - Expected: 0.1-0.8 depending on climate (arid=0.1, humid=0.5-0.8)
   - If > 0.9: Possible cold-start artifact (dt_v014) — discard first 6 months

3. **Seasonal cycle check**: Verify runoff timing matches precipitation timing
   - Expected: Summer peaks for monsoon basins, winter/spring for snowmelt-dominated
   - If inverted (winter Q > summer Q): GW bucket dominance due to insufficient surface runoff (dt_v017)

## Common Pitfalls

> **PITFALL**: Option 5 (Xinanjiang) crashes without MMF_RUNOFF_FILE
> This happens when RUNOFF_OPTION=5 is set but no MMF_RUNOFF_FILE is provided in namelist.hrldas. The symptom is a fatal error: `"Problem opening wrfinput file for read_mmf_runoff"`.
> **Do this instead**: Use RUNOFF_OPTION=1 (TOPMODEL), which provides similar saturation-excess physics without requiring extra files.
> See diagnostic triplet dt_v019 for full details.

> **PITFALL**: REFKDT=3.0 (default) produces near-zero surface runoff in mountain basins
> This happens when Option 3 (Schaake96) is used with default REFKDT=3.0 in mountainous terrain. The symptom is flat baseflow-dominated hydrographs with no storm peaks, and correlation r near 0 with observations.
> **Do this instead**: Reduce REFKDT to 0.3-0.5 for mountain basins. Alternatively, switch to RUNOFF_OPTION=1 (TOPMODEL) which is better suited for saturation-excess terrain.
> See diagnostic triplets dt_v018, dt_v020 for full details.

> **PITFALL**: Choosing runoff option based on VIC comparison without understanding physics differences
> This happens when users expect WRF-Hydro to match VIC exactly by switching to Option 5. The symptom is a crash (MMF file missing) or unexpected results because the two models have different soil layer structures and process representations.
> **Do this instead**: Accept that WRF-Hydro and VIC will differ by 20-40% in magnitude (Bengbu: WRF-Hydro/VIC ratio = 0.71) even with similar physics. Calibrate each model independently against observations.

> **PITFALL**: Not accounting for cold-start spinup artifact when evaluating runoff option performance
> This happens when users compare first-month output against observations. The symptom is 237 mm/day surface runoff on day 1, making any runoff option appear to produce too much water.
> **Do this instead**: Always discard the first 6-12 months of a cold-start simulation as spinup. Evaluate runoff option performance only on the post-spinup period.
> See diagnostic triplet dt_v014 for full details.

---

*This skill document is part of the `hydrocraft-wrfhydro-standalone` knowledge infrastructure.*
*Stage s9 of 12 | Tools used: generate_namelists, build_soil_properties | Related triplets: dt_v017, dt_v018, dt_v019, dt_v020*