Skip to content

Repository files navigation

HydroLegion logo

HydroLegion-HLM

CI License: MIT DOI

A from-scratch Python reimplementation of the Hillslope Link Model (HLM) — GPU-accelerated, ensemble-native solving of network rainfall–runoff ODEs, from a single basin to the globe.

The name pairs Hydro (hydrology) with Legion for its ensemble-native, massively-parallel core: a legion of parameter sets, initial conditions, or storm fields marched together on the GPU.

The Hillslope Link Model (HLM)

The Hillslope Link Model is a distributed rainfall–runoff model developed at the Iowa Flood Center (IFC), University of Iowa. The landscape is decomposed into hillslope–link pairs: each hillslope drains a cascade of soil-moisture buckets into its channel link, and the links are coupled along the tree structure (a directed acyclic graph) of the river network — producing streamflow at every link in the network, not just at the outlet.

The original HLM is implemented in ASYNCH, a numerical solver written in C that integrates the resulting system of ODEs asynchronously over the tree structure of the river network (Small et al., 2013, An asynchronous solver for systems of ordinary differential equations linked by a directed tree structure, Advances in Water Resources). The HLM model formulations are documented in the ASYNCH built-in models reference.

What is HydroLegion-HLM?

HydroLegion-HLM is a refactoring of the original HLM, rewritten from scratch and fully in Python — model physics, numerical solvers, forcing I/O, calibration, and CLI, end-to-end, with no C core. The rewrite includes the ASYNCH tree-structured numerical solver among its solver strategies. It differs from the original HLM in three key ways:

  1. Fully implemented in Python end-to-end on the numpy / scipy / numba / CuPy stack — readable, hackable, pip-installable.
  2. The same physics under multiple, interchangeable numerical solver structures: per-link loops, vectorized and numba-parallel CPU solvers, a decoupled hillslope/channel solver, the ASYNCH-style tree-structured solver, and — importantly — a GPU implementation with CUDA.
  3. A fully-GPU solve of both runoff and routing (see below).

⚡ Fully on the GPU — runoff and routing. In HydroLegion-HLM the entire system is solved on the GPU: hillslope runoff generation is data-parallel over links, channel routing is a GPU sparse matrix multiplication over the network, and an ensemble batch axis marches hundreds of members together — with no host↔device round-trips inside the time loop.

The same physics is implemented at six numerical tiers so you can trade accuracy, speed, and hardware:

Tier What Use
1XX per-link loop (the physics benchmark) reference / teaching
2XX numpy-vectorized (sparse adjacency) the CPU default
3XX numba-parallel upstream sum large-N CPU
4XX decoupled (parallel hillslope + global channel) continental CPU
5XX ASYNCH-style tree routing (vectorized + parallel wavefront) continental CPU
6XX GPU (CuPy / CUDA) — data-parallel march, ensemble batch axis, exp/IMEX integrator calibration / UQ / ensembles

Model families 101–107 span two- and three-bucket hillslopes, nonlinear infiltration, segmented tile drainage (106), and flow-regulating reservoirs (107).

Why? Ensemble-native by design

This redevelopment exists to make massive simulation campaigns practical — thousands of runs instead of one. Three concrete examples, each shipped as a worked recipe:

  1. Ensembles of thousands of simulations — e.g. stochastic storm/forcing ensembles for the probabilistic flood response of a basin (recipe R4).
  2. Calibration — thousands of parameter evaluations against observed streamflow, driven by Optuna or SPOTPY optimizers (recipe R3).
  3. Ensembles of initial conditions — the spread of outcomes from uncertain antecedent states, with per-member initial conditions on the GPU batch axis (recipe R2).

Preliminary benchmarks

Preliminary, single-machine results on the bundled real basin (USGS 07196500, 1126 links):

Configuration Runtime
CPU vectorized (2XX) minutes
CPU decoupled (4XX), 30-day run ~51 s
GPU (6XX), same run ~1–3 s

Test machine (laptop):

Component Spec
CPU Intel Core Ultra 9 185H — 16 cores / 22 threads
RAM 32 GB
GPU NVIDIA GeForce RTX 4060 Laptop GPU — 8 GB VRAM
Software Windows 11 Pro · Python 3.10 · CUDA 12.9 · CuPy 14.x

That is roughly a 50–100× GPU speedup over the CPU tiers — and because the GPU tier carries an ensemble batch axis, it marches hundreds of members in one run (256 members of full model state for this basin is ~13 MB), which is what makes the three scenarios above tractable. Formal, reproducible benchmark suites will be released with HydroLegion-Lab.

Install

pip install hydrolegion-hlm            # CPU tiers + forcing read (numpy/scipy/numba/networkx/joblib/xarray/zarr)
pip install hydrolegion-hlm[gpu]       # + cupy-cuda12x for the GPU tier 6 (NVIDIA / CUDA 12)
pip install hydrolegion-hlm[calib]     # + optuna / spotpy / cmaes / pathos for calibration
pip install hydrolegion-hlm[examples]  # + matplotlib / jupyter / nbmake to run the recipe notebooks
pip install hydrolegion-hlm[docs]      # + jupyter-book to build the documentation site
pip install hydrolegion-hlm[dev]       # + pytest / ruff / build

GPU is always optional — CPU-only users (including macOS / non-NVIDIA) run everything except tier 6.

Quickstart

import numpy as np
import hydrolegion
from hydrolegion.synthetic_net import generate_random_river_network_fast
from hydrolegion.synthetic_forcing import generate_synthetic_rainfall, generate_synthetic_pet

conn, outlet = generate_random_river_network_fast(50, seed=1)
N = len(conn)
param = {"L": np.full(N, 2.0), "Ah": np.full(N, 1.0),
         "k1": 50., "k2": 5., "k3": 20., "k4": 2., "k5": .5,
         "lambda1": .2, "lambda2": -.1, "V0": 1.}                 # family 4 (3-bucket)
ic = np.concatenate([np.zeros(N), np.full(N, 1e-8), np.full(N, 1e-8), np.zeros(N), np.zeros(N)])
P = generate_synthetic_rainfall(duration_days=7)
PET = generate_synthetic_pet(duration_days=7)

res = hydrolegion.run(204, conn, P, PET, param, ic, [0, 7*86400], 3600, max_step=600)   # 2XX vectorized
print(res["Qc"][:, outlet].max())                                # peak outlet discharge (m3/s)

Switch tiers by changing the model number: 104 (loop), 304 (numba), 404 (decoupled), 504 (tree), 604 (GPU — needs [gpu]). hydrolegion info reports versions + GPU availability.

Recipes

Four worked recipes — each on a real basin (USGS 07196500) and a synthetic network you can scale, runnable as a notebook (examples/notebooks/) or from the CLI (examples/cli/*.json):

Recipe Question Tier
R1 — Single run "Simulate this basin for this period." CPU
R2 — IC ensemble "Many initial conditions — the spread of outcomes." GPU 6XX
R3 — Calibration "Fit parameters against observed streamflow." CPU 4XX/5XX or GPU
R4 — Storm ensemble "Many stochastic rainfall fields — the probabilistic response." GPU 6XX

These are worked examples bundled with the core; the full recipe collections for ensembles, calibrations, and initial conditions will ship with HydroLegion-Lab (see Roadmap). The bundled example data (network + daily forcing + observed streamflow) was prepared with HydroLegion-Forge — see examples/data/README.md.

CLI (config-driven)

hydrolegion validate-config examples/cli/r1_single_run_real.json   # check + print the resolved plan
hydrolegion run             examples/cli/r1_single_run_real.json   # single run / GPU ensemble
hydrolegion calibrate       examples/cli/r3_calibration_real.json  # fit to observed streamflow

A single JSON config specifies the model, network, forcing, parameters, time span, solver, optional ensemble, and outputs; the notebook and the CLI share one code path (hydrolegion.config). Every run writes a reproducibility manifest (version + git SHA + hardware + config). Full schema in the docs. Build the docs site with pip install hydrolegion-hlm[docs] && jupyter-book build ..

Correctness

Every tier is held to the others by a pytest suite (tests/):

  • Tier equivalence — 2XX ≡ 3XX (bit-identical), 1XX ≈ 2XX (< 2e-5), 2XX ≈ 4XX/5XX (< 5e-3), 6XX ≈ 2XX (< 1.5e-2, explicit vs BDF).
  • Mass-balance closure — < 0.05% for every family/tier, lumped and distributed forcing.
  • Golden regression — the deterministic 2XX outlet hydrograph is pinned per family.
  • GPU features — per-member ICs / per-member forcing / snapshot-restart bit-exactness (gated on CUDA).
pytest -m "not slow"      # CPU + GPU-equivalence, fast
pytest                    # + the slow GPU feature scripts

Roadmap — the HydroLegion family

This repository is HydroLegion-HLM: the core modeling engine only. Future releases will complete the family:

  • HydroLegion-Forge — creates the input data: river-network extraction and forcing preparation (already used to produce the bundled example data).
  • HydroLegion-Lab — the recipes for ensembles, calibrations, and initial-condition experiments, plus the reproducible benchmark suites.

Author & Citation

Gabriel Perez — Oklahoma State University — Perez-HydroSystems Lab — ORCID 0000-0003-3880-0874

How to cite: see CITATION.cff (GitHub's "Cite this repository" button). DOI (all versions): 10.5281/zenodo.21401928 · DOI (v0.3.0): 10.5281/zenodo.21401929

Credits

HydroLegion-HLM is an independent reimplementation that stands on the foundation of the original Hillslope Link Model and its ASYNCH solver, created by their authors at the Iowa Flood Center, University of Iowa. Full credit for the model formulation and the original asynchronous tree-structured solver belongs to them:

License

MIT — see LICENSE.

About

Hydrological modeling framework with a GPU-accelerated solver for large-scale calibration across hundreds to thousands of parameter sets, uncertainty quantification, ensemble and realization analysis using multiple forcing or parameter scenarios, scenario exploration, and regional reanalysis.

Topics

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Contributors

Languages