ConvLSTM rainfall–runoff: gridded daily meteorology → river discharge

A compact, reproducible deep-learning pipeline that predicts daily discharge at a single gauge from gridded daily meteorological fields over the contributing catchment. Fields (ERA5-Land reanalysis) go into a ConvLSTM encoder; the recent discharge history enters through a separate branch; the two are fused into a single next-day prediction.

The repository contains the method only — no data. The discharge series used during development is proprietary (see Data availability), so everything here is written to run on your basin, your rasters and your gauge.

---

Overview

Classic lumped rainfall–runoff models collapse the catchment into basin-average inputs, discarding where the rain fell and where the snow lies. In a mountainous basin that spatial detail is exactly what drives the hydrograph: a storm over the lower foothills and the same storm volume over the glaciated headwaters produce very different responses, and the snowpack melts from the bottom of the elevation range upward.

A ConvLSTM keeps both dimensions at once:

• convolution over the H × W grid learns where in the basin a signal matters;
• recurrence over the T-day window learns how long it takes that signal to reach the gauge — i.e. it learns catchment lag from data rather than from a prescribed unit hydrograph.

The output is one scalar per window: discharge at the target date.

Why not just an LSTM on basin-mean series? Because averaging is a lossy step taken before the model sees the data. Why not a plain CNN? Because the response is not instantaneous — the same field means different things depending on what preceded it.

---

Method

Input channels and their hydrological meaning

Each day becomes one tile of shape (H, W, 9):

#	Channel	Source	Hydrological role
0	t2m	ERA5-Land 2 m temperature, daily mean (converted to °C)	thermal state; drives melt and evaporative demand
1	swe	ERA5-Land snow water equivalent	how much water is stored in the snowpack, spatially
2	pr	ERA5-Land total precipitation, daily sum	today's input to the system
3–5	pr3, pr7, pr14	rolling sums of pr over 3 / 7 / 14 days	antecedent precipitation index (API) — a proxy for catchment wetness, i.e. how much of the next rainfall becomes runoff instead of storage
6–8	dd3, dd7, dd14	rolling sums of positive degree-days, max(t2m − T₀, 0)	degree-day melt index — the classical temperature-index snowmelt predictor, accumulated over the same three windows

The rolling windows are computed per pixel over the whole time series before windowing, so the network receives wetness and melt-energy state as maps, not as basin averages. T₀ (prep.degday_base_c) and the window lengths (prep.roll_windows) are configurable.

Discharge is deliberately not rasterised into a channel. It is a point quantity of the outlet, and broadcasting it as a constant map only pushes it through the spatial normalisation and convolution filters before global pooling collapses it back to a scalar. It enters as a side input instead.

Architecture

flowchart LR
    A["img (T, H, W, 9)<br/>daily meteo fields"] --> B["ConvLSTM2D 3x3<br/>return_sequences"]
    B --> C[BatchNorm]
    C --> D["ConvLSTM2D 3x3<br/>last state"]
    D --> E[BatchNorm]
    E --> F["Conv2D 3x3 + ReLU"]
    F --> G[GlobalAvgPool2D]
    H["qlag (T, 3)<br/>Q(t-1), Q(t-2), Q(t-3)<br/>for every day of the window"] --> I[Flatten]
    I --> J["Dense 32 + ReLU"]
    G --> K((concat))
    J --> K
    K --> L["Dense 64 + ReLU"]
    L --> M["Dense 1 linear<br/>Q at target date"]

Loading

Plain-text version of the same graph:

img  (T,H,W,9) -> ConvLSTM2D -> BN -> ConvLSTM2D -> BN -> Conv2D -> GAP --+
                                                                          +-> concat -> Dense(64) -> Q
qlag (T,3)     -> Flatten -> Dense(32) ----------------------------------+

Training details: Huber loss (robust to flood-peak outliers), Adam with gradient clipping, ReduceLROnPlateau + EarlyStopping on validation RMSE, dropout and recurrent dropout in both ConvLSTM layers. Discharge targets and the lag branch are scaled with the same StandardScaler, fitted on the training split only; meteorological channels are standardised per channel with training-split statistics.

Pipeline

step 0  clip.py        ERA5-Land GeoTIFFs  ->  clipped to basin polygon
step 1  prep_tiles.py  clipped rasters     ->  daily .npz tiles (H, W, 9)
step 2  dataset.py     tiles + discharge   ->  windows (N,T,H,W,9) + lags (N,T,3) + y
step 3  train.py       windows             ->  trained model, metrics, exports
step 4  inference.py   trained model       ->  forecasts, recursive lag chaining

---

Repository layout

convlstm_runoff/
    config.py       YAML loading, path resolution, --set overrides
    clip.py         step 0: clip rasters to the basin polygon
    prep_tiles.py   step 1: build daily 9-channel tiles
    dataset.py      step 2: windowing, splits, normalisation, CSV reader
    model.py        two-branch ConvLSTM architecture
    metrics.py      RMSE, NSE, persistence baseline
    train.py        step 3: training, evaluation, rolling-origin CV
    interpret.py    saliency / IG / Grad-CAM / occlusion / permutation
    inference.py    step 4: forecasting with recursive discharge lags
configs/
    example.yaml    every path and hyper-parameter lives here

There are no hard-coded paths anywhere in the code — all of them come from the config file or from command-line arguments.

---

Installation

python -m venv .venv && . .venv/bin/activate    # Windows: .venv\Scripts\activate
pip install -r requirements.txt

Tested on Python 3.12 with TensorFlow 2.17 (CPU). rasterio, rioxarray and geopandas are needed only for steps 0–1; if you already have tiles, the rest of the pipeline runs without them.

---

Input data

1. Catchment polygon

An ESRI shapefile of the basin draining to your gauge (paths.basin_shapefile). Any CRS — it is reprojected to the raster CRS automatically. Attributes are irrelevant; only the geometry is used. Derive it from any DEM-based delineation (HydroSHEDS, MERIT Hydro, your own flow-accumulation run).

2. Meteorological rasters

Daily GeoTIFFs, one file per day per variable, in three separate folders. This repository does not download them: get ERA5-Land from the Copernicus Climate Data Store with the cdsapi client, aggregate the hourly product to daily, and export one GeoTIFF per day.

Expected daily variables and file names:

Variable	ERA5-Land source	Aggregation	Default file name
t2m	2 m temperature	daily mean	ERA5L_t2m_daily_mean_YYYY_MM_DD.tif
swe	snow depth water equivalent	12:00 UTC snapshot	ERA5L_swe_hourly_YYYY_MM_DD_12h.tif
tp	total precipitation	daily sum	ERA5L_tp_daily_sum_YYYY_MM_DD.tif

The names are not baked into the code: filenames.*.pattern (regex used to parse the date out of a file name) and filenames.*.template (name written after clipping) are config entries. Change them to match your own export instead of renaming thousands of files.

Grid, CRS and extent may differ between variables — the first t2m raster defines the reference grid and everything else is bilinearly reprojected onto it. Temperature in kelvin is detected and converted to °C automatically.

3. Discharge series

A single CSV (paths.discharge_csv):

date,Q
01.01.2000,100
02.01.2000,105
03.01.2000,98

(illustrative values — substitute your own gauge record)

• date column named either date or data, case-insensitive; format set by data.q_date_format (default %d.%m.%Y, i.e. DD.MM.YYYY);
• discharge column named Q (or anything starting with q);
• decimal comma is accepted (12,5 → 12.5);
• UTF-8 with or without BOM;
• duplicate dates are averaged; unparseable rows are dropped.

Units are never converted — whatever you feed in (m³/s is the natural choice) is what comes out.

Gaps in the record

Missing days are handled explicitly rather than silently: a training window is dropped if its days are not consecutive, if the target date is not exactly pred_horizon days after the window, or if any of Q(t-1..t-3) is missing. The number of windows dropped for each reason is printed per split, e.g.

[dataset] val  : окон собрано 140; отброшено — разрыв внутри окна: 13,
                 разрыв до целевой даты: 1, нет лагов Q: 3

---

Usage

cp configs/example.yaml configs/local.yaml   # then edit the paths inside

Run from the repository root (relative paths in the config resolve against the current directory, or against --base-dir):

python -m convlstm_runoff.clip        --config configs/local.yaml
python -m convlstm_runoff.prep_tiles  --config configs/local.yaml
python -m convlstm_runoff.train       --config configs/local.yaml
python -m convlstm_runoff.inference   --config configs/local.yaml \
    --range-start 2025-01-01 --range-end 2025-03-31 --out-csv out/forecast.csv

Any config key can be overridden without editing the file:

python -m convlstm_runoff.train --config configs/local.yaml \
    --set train.epochs=5 --set cv.enabled=true --set interpret.enabled=true

Forecast mode

inference.py walks forward day by day. While observed discharge is available it is used for the lag branch; past the last observation the model's own prediction is written into memory and becomes the lag for the following step. --prefer-facts keeps observed values in the output instead of overwriting them with predictions. The output CSV carries a source column (fact / forecast) so the two are never confused.

Error handling is deliberate: dates with a missing tile or missing lags are skipped with a [SKIP] message, while genuine failures (channel mismatch, corrupt file, model error) are raised with the offending date attached rather than swallowed.

---

Evaluation

Three things are reported, and the first one matters most:

1. Persistence baseline. "Tomorrow equals today" — Q(t) = Q(t−1). On a daily step in a large basin this baseline is strong, and a model that fails to beat it has learned nothing beyond autocorrelation. It is computed on exactly the same target dates as the model, and it tolerates gaps in the record (missing dates are skipped with a warning instead of aborting the run).

2. RMSE and NSE, both in the original discharge units after inverse-transform. NSE = 1 is perfect; NSE = 0 means "no better than predicting the long-term mean".

3. Rolling-origin cross-validation (cv.enabled: true). For each year y the model is retrained on everything up to y−1 and validated on y, producing out/cv/cv_summary.csv plus per-year prediction files. This is the honest way to estimate skill on a short hydrological record: a single train/val/test split can be lucky or unlucky depending on whether a wet or dry year landed in the test window.

---

Interpretability

A discharge model that cannot be questioned is not usable in an operational hydrological setting — you have to know why it raised the hydrograph before you act on it. interpret.py implements five complementary attribution methods, all operating on the meteorological branch (the lag branch is left at its observed values, so what you see is what the network extracts from meteorology beyond trivial discharge inertia).

Enable with interpret.enabled: true; everything is written to out/explain/.

Gradient-based — what the model is sensitive to right now

• Saliency — |∂Q/∂x| for a single window, reduced to a map (H, W) and to a profile over the T days of the window. Fast, but noisy and prone to gradient saturation. Read it as "where a small perturbation would move the prediction".
• Integrated Gradients — the same derivative integrated along a straight path from a zero baseline to the actual input (interpret.ig_steps steps). More stable than raw saliency and approximately conserving: attributions sum toward the difference in output between baseline and input.
• Grad-CAM — gradients of the output with respect to the last Conv2D feature maps, pooled into channel weights and projected back onto the grid. Coarser than saliency and therefore easier to overlay on a basin map: it answers "which part of the catchment carries this forecast".

Perturbation-based — model-agnostic and free of gradient artefacts

• Occlusion in time — zero out one day of the window and measure the shift in the prediction. The resulting profile is effectively a learned response-time curve: if the peak sits at lag 3–5 days, the network has discovered the catchment's concentration time on its own, and you can compare that against what you know about the basin.
• Occlusion by channel — zero out one channel across the whole window. This is the check that reveals whether the model is actually using snow (swe, dd*) in the melt season or leaning entirely on recent rainfall.
• Occlusion in space — slide a interpret.occlusion_patch-sized patch across the grid, zeroing it and recording the change. Produces an importance map over the catchment. Expensive: H × W forward passes per window.

Dataset-level

• Permutation importance (over time steps and over channels) — shuffle one slice across the batch and measure the increase in RMSE over a subset of the test set (interpret.permutation_subset). Unlike the single-window methods above, this measures importance on the sample as a whole, which is what you want before claiming "SWE matters in this basin". Computed in normalised units, consistently with the model output.

Reading them together is the point: saliency and Grad-CAM say where, occlusion in time says when, permutation says how much it matters overall. Agreement between independent methods is what makes an attribution trustworthy; disagreement is a signal to distrust the model, not the method.

---

Data availability

The discharge series used to develop and validate this pipeline is proprietary and not publicly available; neither it, nor any derived artefact (trained weights, fitted scalers, normalisation statistics, prediction tables, hydrograph plots) is included in this repository, and .gitignore is configured to keep it that way.

The method itself is basin-agnostic. It runs on any catchment for which you have: a basin polygon, daily ERA5-Land fields (globally available, free) and a discharge series at the outlet. Open discharge datasets such as GRDC or the CAMELS family are suitable substitutes for reproducing the workflow.

---

Known limitations

• The full training set is materialised in memory as one (N, T, H, W, C) float32 array. Fine for a small basin at ERA5-Land resolution; for a large grid or a long record, convert dataset.py to a tf.data generator.
• recurrent_dropout > 0 disables the cuDNN fast path for ConvLSTM — expect slow epochs on GPU. Set it to 0 if throughput matters more than regularisation.
• NaNs in a raster are filled with the scene mean, and no valid-pixel mask is propagated to the model.
• In recursive forecast mode, error accumulates with distance from the last observation, and the normalisation statistics come from the training period only — predictions far outside the training range are extrapolation.

---

Note

Developed with AI-assisted coding. The method design, the choice of predictors and the hydrological validation are the author's.

License

MIT — see LICENSE.