Add dynamic tailwater elevation with outflow-dependent TW rating curves

**Source Structure**:

- `src/lib.rs` - Main module: data structures (`ReservoirData`, `RiverData`, `ConfluenceData`), `simulate_cascade()` orchestration, per-timestep reservoir simulation

- `src/helpers.rs` - Interpolation functions: `interpolate_elevation()`, `interpolate_storage()`, bilinear HPF interpolation for power/release computation

- `src/helpers.rs` - Interpolation functions: `interpolate_elevation()`, `interpolate_storage()`, `interpolate_tailwater()`, bilinear HPF interpolation for power/release computation

- `python/powersheds/__init__.py` - Package init, re-exports from `_lib`

- `python/powersheds/diagnostics.py` - Interactive diagnostic visualization (see Diagnostics section)

**Three Object Types**:

- **Reservoir**: Pool elevation from storage-elevation curves, head calculation, release/power via bilinear HPF interpolation, storage/capacity constraints

- **Reservoir**: Pool elevation from storage-elevation curves, head calculation (with optional dynamic tailwater), release/power via bilinear HPF interpolation, storage/capacity constraints

- **River**: Lag-based flow routing with legacy flows for initialization

- **Confluence**: Instantaneous flow merging (no lag)

**Key Algorithm Flow** (per timestep):

1. Sort objects by simulation order

2. For each reservoir: interpolate pool elevation → compute head → bilinear HPF lookup for target release → constrain release → update storage → compute actual power

2. For each reservoir: interpolate pool elevation → compute head (with dynamic tailwater iteration if TW table provided) → bilinear HPF lookup for target release → constrain release → update storage → compute actual power

3. Rivers delay flows by `lag` hours; confluences pass through immediately

4. Upstream releases/spills become downstream inflows

**Dynamic Tailwater**: When `set_tw_outflow`/`set_tw_elevation` are provided, tailwater elevation is computed from total outflow via a rating curve, using fixed-point iteration (tolerance=0.001m, max_iter=10, damping=0.5). When empty, the static `tailwater_elevation` scalar is used (legacy behavior).

## Units

- Storage: Mm³ (Million Cubic Meters)

- Head/Elevations: meters

- Power: MW

- HPF flows (`hpf_q`): m³/s (cumecs)

- TW outflow (`set_tw_outflow`): m³/s (cumecs)

## Diagnostics

The `powersheds.diagnostics` module provides interactive visualization for simulation results. Requires the `viz` optional dependency: `uv pip install -e ".[viz]"`

**Public API** (`from powersheds.diagnostics import ...`):

- `cascade_diagnostics(results, ...)` → Plotly `go.Figure` with dropdown reservoir selector, three linked time series panels (flows, elevation, power), and infeasible power shading

- `cascade_diagnostics(results, cascade_data=None, ...)` → Plotly `go.Figure` with dropdown reservoir selector, three linked time series panels (flows with release constraint band, elevation with pool/tailwater/head, power with infeasible shading). Pass `cascade_data` to enable min/max release bands on the flows panel.

- `network_diagram(cascade_data, ...)` → matplotlib `Figure` with publication-quality static network diagram (trapezoid reservoir nodes, Bezier edges, lag labels)

- `save_diagnostics_html(results, output_path, ...)` → writes standalone HTML file

# Rebuild the module first

uv pip install -e .

# Run all tests (148 tests)

# Run all tests (158 tests)

uv run python -m pytest tests/ -v

# Run a specific test file

uv run python -m pytest tests/test_cascade.py -v

```

**Test Coverage** (148 tests across 8 files):

**Test Coverage** (158 tests across 9 files):

- `tests/test_hpf.py` - HPF bilinear interpolation (74 tests): grid corner/edge/interior roundtrips, statistical robustness, performance benchmarks, edge cases

- `tests/test_diagnostics.py` - Diagnostics visualization (28 tests): topology extraction, network layout, axis ranges, network_diagram matplotlib output, cascade_diagnostics Plotly output, HTML export, Cumberland integration, edge cases

- `tests/test_diagnostics.py` - Diagnostics visualization (29 tests): topology extraction, network layout, axis ranges, network_diagram matplotlib output, cascade_diagnostics Plotly output (incl. release constraint bands), HTML export, Cumberland integration, edge cases

- `tests/test_cascade.py` - Full cascade integration (13 tests): output structure, multi-reservoir routing, simulation ordering, mass balance, physical invariants, Cumberland NaN check

- `tests/test_simulate_timestep.py` - Per-timestep reservoir constraints (10 tests): below power pool, min/max release, insufficient water, spill, mass balance, unconstrained power equality

- `tests/test_interpolation.py` - Storage-elevation interpolation (10 tests): breakpoints, midpoints, clamping, roundtrip, monotonicity, Cumberland validation

- `tests/test_dynamic_tailwater.py` - Dynamic tailwater (9 tests): flat TW matches static, higher TW reduces head, TW output in range, convergence on steep curve, mass balance, Cumberland integration (bounds, Barkley flat, no NaN, head positive)

- `tests/test_river_confluence.py` - Flow routing (6 tests): river lag, legacy flows, lag=0 passthrough, legacy length mismatch, confluence summation

- `tests/test_regression.py` - Regression tests (2 tests): Cumberland golden baseline, CLAUDE.md smoke test

- `tests/test_unit_conversion.py` - Unit conversion (2 tests): cumecs to Mm3/hr at known grid points

%% Cell type:markdown id: tags:

# Powersheds Cumberland Cascade Diagnostics

This notebook loads the 8-reservoir Cumberland River Basin cascade,

runs the Powersheds simulation, and produces:

1. A **static network diagram** showing the cascade topology

2. An **interactive time series diagnostic** (Plotly) with dropdown reservoir selection

%% Cell type:code id: tags:

``` python

from dataclasses import dataclass

from pathlib import Path

import pandas as pd

import yaml

import powersheds

from powersheds.diagnostics import cascade_diagnostics, network_diagram, save_diagnostics_html

```

%% Cell type:markdown id: tags:

## 1. Load cascade configuration and time series

%% Cell type:code id: tags:

``` python

from dataclasses import dataclass, field

@dataclass

class ReservoirData:

    object_type: str

    capacity: float

    initial_pool_elevation: float

    min_power_pool: float

    set_storage: list[float]

    set_elevation: list[float]

    hpf_h: list[float]

    hpf_p: list[float]

    hpf_q: list[float]

    tailwater_elevation: float

    max_release: float

    min_release: float

    catchment_inflow: list[float]

    target_power: list[float]

    simulation_order: int

    downstream_object: str

    set_tw_outflow: list[float] = field(default_factory=list)

    set_tw_elevation: list[float] = field(default_factory=list)

@dataclass

class RiverData:

    object_type: str

    simulation_order: int

    downstream_object: str

    lag: int

    legacy_flows: list[float]

@dataclass

class ConfluenceData:

    object_type: str

    simulation_order: int

    downstream_object: str

@dataclass

class CascadeData:

    reservoirs: dict[str, ReservoirData]

    rivers: dict[str, RiverData]

    confluences: dict[str, ConfluenceData]

```

%% Cell type:code id: tags:

``` python

CUMBERLAND_DIR = Path(".").resolve()

with open(CUMBERLAND_DIR / "cascade_config.yaml") as f:

    config_dict = yaml.safe_load(f)

reservoir_dict = {}

river_dict = {}

confluence_dict = {}

for name, specs in config_dict.items():

    if specs["object_type"] == "reservoir":

        ts = pd.read_csv(CUMBERLAND_DIR / "time_series" / f"{name}.csv")

        se = pd.read_csv(CUMBERLAND_DIR / "storage_HWelevation_tables" / f"{name}.csv")

        hpf = pd.read_parquet(CUMBERLAND_DIR / "head_power_flow_tables" / name)

        tw_path = CUMBERLAND_DIR / "outflow_TWelevation_tables" / f"{name}.csv"

        if tw_path.exists():

            tw = pd.read_csv(tw_path)

            set_tw_outflow = tw["outflow_cumecs"].tolist()

            set_tw_elevation = tw["elevation_m"].tolist()

        else:

            set_tw_outflow, set_tw_elevation = [], []

        reservoir_dict[name] = ReservoirData(

            **specs,

            catchment_inflow=ts["catchment_inflow"].tolist(),

            target_power=ts["target_power"].tolist(),

            set_storage=se["storage_Mm3"].tolist(),

            set_elevation=se["elevation_m"].tolist(),

            hpf_h=hpf["H_m"].astype(float).tolist(),

            hpf_p=hpf["P_MW"].astype(float).tolist(),

            hpf_q=hpf["Q_cumecs"].astype(float).tolist(),

            set_tw_outflow=set_tw_outflow,

            set_tw_elevation=set_tw_elevation,

        )

    elif specs["object_type"] == "river":

        legacy = pd.read_csv(CUMBERLAND_DIR / "time_series" / f"{name}.csv")

        river_dict[name] = RiverData(

            **specs,

            legacy_flows=legacy["legacy_flow"].tolist(),

        )

    elif specs["object_type"] == "confluence":

        confluence_dict[name] = ConfluenceData(**specs)

cascade_data = CascadeData(

    reservoirs=reservoir_dict,

    rivers=river_dict,

    confluences=confluence_dict,

)

print(f"Loaded {len(reservoir_dict)} reservoirs, {len(river_dict)} rivers, {len(confluence_dict)} confluences")

```

%% Cell type:markdown id: tags:

## 2. Run the simulation

%% Cell type:code id: tags:

``` python

results = powersheds.simulate_cascade(cascade_data)

print(f"Simulation complete: {len(results)} objects")

for name, data in results.items():

    if isinstance(data, dict) and "actual_power" in data:

        n = len(data["actual_power"])

        print(f"  {name}: {n} hours")

```

%% Cell type:markdown id: tags:

## 3. Network diagram

Static visualization of the cascade topology. Trapezoid nodes represent

reservoirs (amber bar = hydropower), circles represent confluences, and

edges show river connections with lag times in hours.

%% Cell type:code id: tags:

``` python

network_diagram(cascade_data, title="Cumberland River Basin Cascade");

```

%% Cell type:markdown id: tags:

## 4. Interactive time series diagnostics

Three linked panels (Flows, Elevation, Power) with shared x-axis zoom/pan.

Use the **dropdown** to switch between reservoirs.

Red shading in the Power panel highlights **infeasible generation** where

actual power falls below the target.

- **Flows panel**: release, target release, inflow, spill, with grey band showing min/max release constraints

- **Elevation panel**: pool elevation, tailwater elevation, and head

- **Power panel**: target vs actual power, with red shading for **infeasible generation**

%% Cell type:code id: tags:

``` python

fig = cascade_diagnostics(

    results,

    cascade_data=cascade_data,

    title="Cumberland River Basin \u2014 Cascade Diagnostics",

    height=900,

    width=1100,

)

fig.show()

```

%% Cell type:markdown id: tags:

## 5. Export to standalone HTML

Save the interactive time series as a self-contained HTML file.

%% Cell type:code id: tags:

``` python

output_path = save_diagnostics_html(

    results,

    output_path="cumberland_diagnostics.html",

    title="Cumberland River Basin \u2014 Cascade Diagnostics",

)

print(f"Saved: {output_path.resolve()}")

```

outflow_cumecs,elevation_m

0,10.973

566.34,10.973

1132.68,10.973

1699.02,10.973

2265.36,10.973

2831.7,10.973

3398.04,10.973

3964.38,10.973

4530.72,10.973

5097.06,10.973

5663.4,10.973

6796.08,10.973

7928.76,10.973

9061.44,10.973

10194.12,10.973

11326.8,10.973

12459.48,10.973

13592.16,10.973

outflow_cumecs,elevation_m

0,143.256

141.585,146.304

283.17,147.676

424.755,148.895

566.34,149.962

707.925,150.724

849.51,151.486

1132.68,153.01

1415.85,154.381

1699.02,155.448

1982.19,156.515

2265.36,157.277

2548.53,158.039

2831.7,158.648

4247.55,161.087

5663.4,163.068

7079.25,164.897

8495.1,166.421

9910.95,167.945

11326.8,169.164

12742.65,170.383

outflow_cumecs,elevation_m

0,110.642

566.34,111.692

1132.68,113.229

1699.02,114.3

2265.36,115.473

2831.7,116.386

3398.04,117.409

3964.38,118.181

7390.737,121.107