Hydro revamp: inflow-driven representation via atlite Cutout.hydro with PUDL plant data

[ktehranchi]

Motivation

Conventional hydro in PyPSA-USA currently comes from the 2016-vintage Breakthrough/TAMU dataset, attached as non-extendable Generators with the historical hourly dispatch as p_max_pu (add_electricity.py, attach_breakthrough_renewable_plants, with the in-code note "temporarily adding hydro with breakthrough only data until I can correctly import hydro_data"). Consequences of the current representation:

• PUDL/EIA-860 "Conventional Hydroelectric" generators are not attached: hydro is listed in renewable_carriers, so attach_conventional_generators filters them out, and no other path picks them up.
• There is no reservoir representation: no StorageUnit, no inflow (storage_units_t.inflow is never populated), no energy budget, no run-of-river vs storage distinction. The optimizer sees a fixed 2016 dispatch shape regardless of weather year.
• Existing pumped storage (via prime_mover_code == 'PS') uses a uniform hardcoded max_hours = 24/0.894 for all plants; EIA-860 energy_storage_capacity_mwh is NULL for all hydro/PHS generators, so no data-driven duration is currently available.
• The renewable: hydro: config block (config.common.yaml) is not wired to any code path (build_renewable_profiles excludes hydro by wildcard constraint), and the docs describe an atlite-inflow implementation that does not exist.

This issue proposes replacing the Breakthrough hydro pathway with a weather-year-consistent, inflow-driven representation built on atlite.Cutout.hydro, following the approach used in pypsa-eur / pypsa-earth / pypsa-northamerica (ERA5 runoff routed through HydroBASINS catchments, normalized to reported generation), adapted to PyPSA-USA's PUDL-based plant data and US-specific datasets.

Prior art

pypsa-earth (and pypsa-northamerica, which inherits its implementation) builds per-plant inflow with atlite.Cutout.hydro: plants are point-joined to HydroBASINS level-6 basins, ERA5 runoff is aggregated over each plant's upstream basin set with a travel-time lag, and the result is normalized to national annual generation statistics (EIA/IRENA). Plants attach as ror (Generator), hydro (StorageUnit with inflow), and PHS (StorageUnit). This issue adopts that basic architecture and proposes US-specific choices for the plant data, normalization targets, and storage parameters.

Feasibility (verified)

• atlite==0.3.0 (our pin) includes Cutout.hydro (atlite/convert.py, atlite/hydro.py). Requirements: a plants frame with lon/lat (index becomes the output plant coordinate), a HydroBASINS layer with HYBAS_ID/NEXT_DOWN/DIST_MAIN/geometry, and a cutout with the runoff feature.
• The published PyPSA-USA ERA5 cutouts (Zenodo 14611937) already contain the runoff feature (verified from the file headers); build_cutout.py prepares all ERA5 features.
• PUDL v2025.5.0 (out_eia__yearly_generators, which build_powerplants.py already queries): 4,177 HY generators / 78.0 GW → 1,481 plants after grouping to plant_id_eia; 145 PS generators / 20.8 GW → 38 plants; lat/lon present for nearly all.
• 100% of CONUS HY plants fall inside a HydroBASINS NA level-6 polygon (9 plants / 38 MW need a snap-or-drop rule); Alaska plants fall in the ar tile.
• EIA-923 via PUDL (out_eia923__monthly_generation_fuel_combined, prime_mover_code in ('HY','PS')) covers 98.1% of hydro capacity; 97.9% of capacity has a complete 12-month series for 2019. National 2019 total (287.8 TWh) matches EIA's published figure.
• PyPSA 1.3's default one-port clustering strategies aggregate StorageUnit inflow by sum and max_hours capacity-weighted, so plant-level inflow survives cluster_network without new configuration.

Design decisions specific to PyPSA-USA

1. Site-level inflow, plant-level normalization. Group generators to plant_id_eia before calling Cutout.hydro (one point per site), and calibrate each plant's inflow so its annual sum matches its own EIA-923 net generation, with a fallback hierarchy plant → state → national for the ~2% of capacity without a usable target. Rationale: 200 level-6 basins host more than one plant (one hosts 51, on the mid-Columbia cascade); per-plant calibration gives cascaded plants the same inflow shape with individually correct levels, and makes the co-located-unit question moot.

2. Annual calibration by default; monthly only for monthly reporters. ~90% of US hydro plants report EIA-923 generation annually; EIA imputes their monthly split from census-division allocation factors (Turner, Voisin & Nelson 2022, Sci Data 9:675, doi:10.1038/s41597-022-01748-x). Measured for 2019: monthly (M) reporters = 39.7 GW (50.6% of capacity), annual (A) = 38.2 GW. Monthly calibration should therefore be gated on reporting_frequency_code in ('M','AM'); applying it to A-coded plants would imprint the imputation climatology rather than hydrology.

3. Explicit efficiency accounting. EIA-923 reports metered net generation, while a StorageUnit discharges through efficiency_dispatch. Reservoir inflow must be inflow_gen / efficiency_dispatch, with the divisor taken from the same cost-table variable used for the component (no literals). Run-of-river (Generator) takes p_max_pu = inflow_gen / p_nom clipped to 1.0, with no division. Invariant to enforce in a test: with unconstrained dispatch and no spill, annual reservoir output equals the EIA-923 target.

4. Reservoir energy capacity from ORNL HESC v2. Hydropower Energy Storage Capacity Dataset v2 (doi:10.21951/HESC/1972462) provides per-plant storage energy (MWh) and discharge duration (hours), keyed on plant_id_eia/EHA_PtID/NIDID, in four confidence tiers. Use the best available tier per plant; fall back to NID storage-volume × head (via the HILARRI crosswalk) and then to a state/mode median — never a single global default. Record the source tier per plant.

5. Run-of-river vs storage classification from ORNL EHA. The Existing Hydropower Assets plant database (doi:10.21951/EHA_FY2024/2344934) carries an operating-mode classification (ROR/storage/pumped) and EIA_PtID, replacing any inference from EIA-860 (which has no such field). HESC duration (<~2 h ⇒ effectively run-of-river) is the secondary signal.

6. PHS durations from NID via HILARRI. Derive per-plant max_hours from upper-reservoir normal storage and head where available (HILARRI v4, doi:10.21951/HILARRI/3362570, crosswalks plant_id_eia ↔ NIDID ↔ GRanD ↔ NHD COMID); calibrate the fleet total against ORNL Hydropower Market Report figures (~21.9 GW / ~553 GWh ≈ 25 h average); replace the uniform hardcoded value with a config default plus per-plant overrides. PHS receives no inflow and is excluded from all normalization denominators.

7. A dedicated runoff cutout with an extended bounding box. The usa cutout's northern edge (50°N) truncates the upstream basins of 73 plants / 21.0 GW (26.9% of conventional hydro capacity) — Grand Coulee alone has 26.6% of its upstream area in Canada — and atlite returns zero (not NaN) for basins outside the cutout. Rather than growing the main cutouts (invalidating ~40 GB of published artifacts), build a runoff-only cutout at x:[-127,-64], y:[23,56], dx=dy=0.5° (runoff is basin-averaged, so 0.5° is sufficient; ~1/20 the size of the existing cutouts). The build must assert that every plant's upstream basin set lies within the cutout bounds and fail otherwise. This also decouples hydro from renewable.dataset, under which (godeeep, the default) no cutout is currently in the DAG at all.

8. HydroBASINS retrieved at build time. HydroBASINS v1c is under the HydroSHEDS v1 license (free use, redistribution restricted), so retrieve hybas_na_lev06_v1c (and ar for Alaska; deduplicate on HYBAS_ID in the 51–63°N overlap band) directly from data.hydrosheds.org rather than re-hosting. HESC/EHA/HILARRI are DOE open-access with a citation requirement and can be staged into the existing Zenodo data bundle.

Implementation plan

Phase 1 — inflow pipeline and attachment (replaces the Breakthrough hydro dependency)

1. rule retrieve_hydrobasins (retrieve.smk) + additions to the Zenodo bundle for HESC v2, EHA, HILARRI v4.
2. Runoff-only cutout (build or publish alongside the existing Zenodo cutouts).
3. rule build_hydro_inflow (build_electricity.smk), national scope, output routed through PROFILES:
   powerplants.csv → filter HY → dedup to plant_id_eia → sjoin to HydroBASINS (snap-or-drop with logged MW) → bounds assertion → Cutout.hydro(flowspeed=1.0) → per-plant EIA-923 annual calibration with plant→state→national fallback → reservoir inflow ÷ efficiency_dispatch → clip_min_inflow (post-normalization) → hydro_inflow.nc with inflow(plant, time) [MW] plus per-plant norm_source, norm_year, basin_id, n_upstream_basins and provenance attributes.
4. Attachment in add_electricity.py: remove attach_breakthrough_renewable_plants(..., ["hydro"]) and the associated bus2sub remap; use the existing match_plant_to_bus path.
   • ror → Generator, ATB ror costs (mapping already exists in constants.py).
   • hydro → StorageUnit with inflow, per-plant max_hours (HESC → NID → median, tier logged), efficiency_store=0, cyclic SOC per period.
   • PHS → StorageUnit as today with per-plant max_hours (item 6).
   • Split EIA_TECH_MAP so Conventional Hydroelectric → hydro and Hydroelectric Pumped Storage → PHS (currently both map to hydro); update RPS_CARRIERS, cluster-weighting carrier sets, and the plotting maps (colors/nice-names for hydro/ror/PHS already exist).
5. Wire the renewable: hydro: config block for real (via config_provider), update config-configuration.md and the # docs : HYDRO marker together (docs config test enforces this), and keep workflow/repo_data/config in sync.

Config sketch:

renewable:
  hydro:
    cutout: era5_runoff
    carriers: [ror, PHS, hydro]
    resource: {method: hydro, hydrobasins_level: 6, flowspeed: 1.0}
    normalization:
      target: eia923
      resolution: annual        # monthly applies only to M/AM reporters
      year: 2019                # default = weather year; warn if overridden
      fallback: [plant, state, national]
    max_hours: {source: hesc, phs_default: 25.0}
    clip_min_inflow: 1.0

Phase 2 — operational bounds and validation harness

• Minimum-generation / monthly power-target bounds from HydroWIRES B1 (Zenodo 13351949, CC0; monthly and weekly max/min generation for 1,492 plants, 2001–2022) as an opts/ constraint module alongside reserves.py/policy.py.
• Optional smooth (168 h rolling mean, as in pypsa-eur's build_hydro_profile) settled empirically rather than inherited (see validation below).

Phase 3 — extensions

• Cascade coupling (upper-plant release as lower-plant inflow) using HILARRI NHD COMIDs + HydroBASINS NEXT_DOWN topology.
• GloFAS routed discharge in place of raw ERA5 runoff (Cutout.hydro(module="glofas"); requires an atlite upgrade beyond 0.3.0). ERA5 runoff is not routed between grid cells, and published evaluations report degraded skill in snowmelt-dominated and heavily regulated basins, so this is the main accuracy upgrade path.
• GODEEEP-hydro for future-climate scenarios (a GODEEEP ingestion path already exists for wind/solar); Canadian hydro via the PNNL Western Canada monthly dataset (Zenodo 14984725, 110 BC/AB facilities, 1981–2019, same schema as B1).

Validation plan

Targets chosen to be independent of the calibration input:

1. Unit invariant (test): unconstrained annual reservoir output equals the plant's EIA-923 target; RoR annual p_max_pu · p_nom sum equals its target. Guards the efficiency accounting.
2. Monthly shape vs RectifHyd/RectifHydPlus (Zenodo 11584567 / ORNL doi:10.21951/rectifhydplus/2440153): compare modeled monthly inflow energy per plant against the observation-derived series (built from reservoir releases and USGS gauges; median KGE 0.78 vs. observations). Report KGE/NSE distributions by region; the Columbia, Sierra snowmelt, and Southeast fleets are the diagnostic cases. Raw EIA-923 monthlies are not used as a shape target because of the imputation issue above.
3. Hourly/weekly aggregate shape vs EIA-930: BA-level hydro (core_eia930__hourly_net_generation_by_energy_source, hydro_excluding_pumped_storage, already in PUDL) compared against solved dispatch and against inflow aggregates, extending the existing gridemissions validation plots.
4. Spot checks vs USGS streamflow: for a sample of large plants, compare inflow timing against downstream gauge discharge via HILARRI COMID → NLDI navigation → NWIS daily values (new api.waterdata.usgs.gov OGC API; the legacy WaterServices endpoint is scheduled for decommissioning in Q1 2027).
5. Held-out-year test for the calibration design: calibrate on year A, evaluate shape/level on year B, to quantify how much single-year plant-level scaling bakes in operational (non-hydrological) variance; informs whether state/basin-level scaling should be preferred where plant-level residuals are large.
6. Clustering invariance (test): total inflow energy and total max_hours-weighted storage are preserved through simplify/cluster_network.
7. System-level A/B: solved runs (current Breakthrough representation vs. new) at the standard validation configuration, comparing hydro generation by month/region against EIA-930 and emissions/production validation figures.

Open questions

1. Weather year vs normalization year when they differ (future horizons): calibrate on a multi-year EIA-923 mean and preserve the ERA5 interannual signal, or pin to a reference year?
2. Monthly calibration policy: mixed by reporting_frequency_code (recommended above), annual-only, or adopting B1/RectifHyd as the shape target for all plants (best shape, adds a dependency)?
3. Distribution of HESC/EHA/HILARRI: stage into the project Zenodo bundle vs. retrieve from ORNL HydroSource (registration-gated) at build time.
4. Whether to also extend the main usa cutout bbox eventually, vs. keeping the dedicated runoff cutout indefinitely.

Related cleanups noticed during this investigation (separate PRs)

• The PHS ITC in apply_itc matches carrier PHS against keys 8hr_PHS/10hr_PHS/12hr_PHS and therefore never applies.
• build_powerplants's PUDL S3 read is not declared as a Snakemake input.
• docs/source/config-configuration.md describes an atlite hydro implementation ahead of this issue being implemented.