cs.CVSep 26, 2026

Distributed Hydrological Modeling in the Feature Space

Authors: Mohamad Hakam Shams Eddin, Maria Luisa Taccari, Yikui Zhang, Shijie Jiang, Juergen Gall, Markus Reichstein

Organizations: Max Planck Institute for Biogeochemistry, Germany · European Centre for Medium-Range Weather Forecasts (ECMWF), United Kingdom · Institute of Bio- and Geosciences Agrosphere (IBG-3), Research Centre Jülich, Germany · University of Bonn, Germany · Lamarr Institute for Machine Learning and Artificial Intelligence, Germany

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Accurate forecasting of river discharge and floods is very challenging. River dynamics are affected by storage, meteorological forcing, and flow propagation at different spatial and temporal scales. Forecasting thus requires a framework that considers the upstream-to-downstream flow through river networks across grid cells and catchments. This modeling is known in hydrology as distributed modeling and routing. Existing deep learning approaches either ignore this topology, operate on lumped catchments, or route predicted physical quantities through a separate graph or physical routing model. We instead introduce feature-space routing: a topology-aware state-space operator embedded directly in the forecasting dynamics. At every forecast step, the operator gathers latent states from upstream grid cells and causally updates the downstream state according to the known river network. This preserves the physical connectivity of the river system while allowing the propagated state itself to be learned end-to-end and allows the model to predict river discharge considering both local dynamics and neighboring upstream contributions. To address uncertainty and provide probabilistic forecasts, we minimize the fair continuous ranked probability score (fCRPS) as a training objective. Our experiments on the European Flood Awareness System (EFAS) and observational data for river discharge forecasting demonstrate that encoding the physical structure of river networks explicitly in the feature space substantially improves the forecasting skill, particularly in an ungauged setting. Our approach achieves state-of-the-art results on both reanalysis and observational data and is able to forecast maps of river discharge at 1 arcminute and 6-hourly resolution up to 10 days lead time.

Figures & tables

Figure 1: Example river discharge and flood forecasting over mid Europe in Sep. 2024. In the middle is the forecast produced by our model with routing at 4-day lead time, compared to the same model but without routing mechanism on the left and EFAS reanalysis data as a reference on the right side. The flood categories are shown below, e.g., high river discharge with 2020 -year return period. Without explicit distributed modeling, the routing-free model produces heterogeneous predictions in which neighboring pixels can have inconsistent discharge, leading to many false positives. Our routing module instead yields more spatially coherent predictions that better follow the river network.
Figure 2: Overview of the routing module. (a) Shows the Local Drain Direction (LDD) which is used to generate the graph of river network (b). Each downstream point forms a sequence with its upstream neighbors and routing is conducted via SSD on concatenated point sequences at once (c).
Figure 3: An overview of the proposed model for river discharge forecasting. The model initializes the forecast at time tt and forecasts high-resolution river discharge maps Xdis06t:t+L\textbf{X}_{dis06}^{t:t+L} from the initial conditions ( XANAt−4(T−1):t\textbf{X}_{ANA}^{t-4(T-1):t} , XIMERGt−4(T−1):t\textbf{X}_{IMERG}^{t-4(T-1):t} , XCPCt−4T:t−4\textbf{X}_{CPC}^{t-4T:t-4} , XEMOt−4(T+1):t−8\textbf{X}_{EMO}^{t-4(T+1):t-8} ), AlphaEarth embedding ( XAlphaE\textbf{X}_{AlphaE} ), and weather forecasts ( XHRESt+1:t+L\textbf{X}_{HRES}^{t+1:t+L} ), conditioned on a noise vector ϵz\epsilon_{z} to generate probabilistic forecast.
Figure 4: The structure of the forecast layer. It has a similar structure as the encoder layer, but it additionally incorporates meteorological forecasts (HRES) by concatenation, the noise ϵz\epsilon_{z} , and the routing module. The bi-directional SSD block is shown in more details.
   Validation (2022-2023) Test (2024-2025)
 Model R2 r ( ↑\uparrow )KGE r ( ↑\uparrow )CRPS r ( ↓\downarrow ) R2 r ( ↑\uparrow )KGE r ( ↑\uparrow )CRPS r ( ↓\downarrow )
ProbabilisticDRUM (mean) 0.14960.15052.0872 0.09840.13522.3581
ED-LSTM (mean) 0.17260.24311.9121 0.13470.23332.1971
ME-LSTM (mean) 0.20930.30711.8395 0.12100.26402.2249
Ours (mean) 0.45150.57511.4720 0.39070.54741.7619
Table 1: Results on EFAS Reanalysis (CAMELS split).
Figure 5: Results on EFAS reanalysis with lead time.
(a) Routing and Warmup (b) Routing Module (c) Upstream Aggregation
RoutingWarmupKGE r | KGE o | F1 r ( ↑\uparrow ) ModuleKGE r | KGE o | F1 r ( ↑\uparrow ) TypeKGE r | KGE o | F1 r ( ↑\uparrow )
✗✗0.666 | 0.692 | 0.096 GS-mean0.680 | 0.721 | 0.310 Summation0.738 | 0.682 | 0.335
✓✗0.672 | 0.739 | 0.114 GS-sum0.702 | 0.702 | 0.285 Mean0.730 | 0.739 | 0.332
✗✓0.687 | 0.721 | 0.162 GAT0.714 | 0.633 | 0.268 GraphMamba0.722 | 0.720 | 0.276
✓✓0.740 | 0.747 | 0.339 SSM0.740 | 0.747 | 0.339 Ours0.740 | 0.747 | 0.339
Table 2: Ablation studies on the validation set over Mid Europe. Shown are KGE on reanalysis (KGE r ) and observational (KGE o ) CAMELS data along with F1 r @(1.5-500) for reanalysis over the full domain ( 281281 K points).
 ValTest
RiverMamba0.141 ±\pm 0.010.138 ±\pm 0.01
Ours0.334 ±\pm 0.020.314 ±\pm 0.01
Table 3: F1 comparison to RiverMamba.
   Validation (2022-2023) Test (Jun 2024- Dec 2025)
 Model R2 o ( ↑\uparrow )KGE o ( ↑\uparrow )CRPS o ( ↓\downarrow ) R2 o ( ↑\uparrow )KGE o ( ↑\uparrow )CRPS o ( ↓\downarrow )
 EFAS (control run) --- 0.04100.3964-
ProbabilisticDRUM (mean) 0.17500.28344.5375 -0.00890.23554.3474
ED-LSTM (mean) 0.28410.39573.9024 0.12140.37883.6818
ME-LSTM (mean) 0.33390.47713.4928 0.09130.38433.7962
Ours (mean) 0.60420.66892.4592 0.36930.55202.9708
Table 5: Results on gauged GRDC observations.
Figure 6: Results on GRDC observations with lead time.
Appendix figures & tables109 assets

Supplementary material from the paper’s appendix.

Appendix

Results on EFAS Reanalysis
StaticCPCIMERGEMOAnalysisIFS-HRESR2 cr{}_{c}^{r} | KGE cr{}_{c}^{r} | R2 gr{}_{g}^{r} | KGE gr{}_{g}^{r} | F1 r ( ↑\uparrow )
✗✓✓✓✓✓0.6963 | 0.7338 | 0.6919 | 0.7443 | 0.2734
✓✗✓✓✓✓0.6741 | 0.6972 | 0.6674 | 0.7145 | 0.3053
✓✓✗✓✓✓0.7021 | 0.7303 | 0.7041 | 0.7411 | 0.3100
✓✓✓✗✓✓0.6863 | 0.7474 | 0.6637 | 0.7116 | 0.3101
✓✓✓✓✗✓0.6869 | 0.6979 | 0.6835 | 0.7107 | 0.3069
Appendix
Table 6: Ablation studies for feature importance on the validation set over Europe. Shown are KGE for CAEMLS (KGE cr{}_{c}^{r} ) and GRDC (KGE gr{}_{g}^{r} ) splits on reanalysis. along with the F1 r @(1.5-500) for reanalysis over the full domain ( 281281 K points). In the second half of the table, we report the results of R2 co{}_{c}^{o} , KGE co{}_{c}^{o} , and F1 co{}_{c}^{o} @(1.5-20) on the fine-tunning on CAMELS observations.
Results on EFAS Reanalysis
LISFLOOD StaticAlphaEarthR2 cr{}_{c}^{r} | KGE cr{}_{c}^{r} | R2 gr{}_{g}^{r} | KGE gr{}_{g}^{r} | F1 r ( ↑\uparrow )
✗✗0.6963 | 0.7338 | 0.6919 | 0.7443 | 0.2734
✓✗0.7045 | 0.7404 | 0.6943 | 0.7462 | 0.3360
✗✓0.7061 | 0.7396 | 0.7043 | 0.7493 | 0.3394
Results on CAMELS observational data
LISFLOOD StaticAlphaEarthR2 co{}_{c}^{o} | KGE co{}_{c}^{o} | F1 co{}_{c}^{o} ( ↑\uparrow )
Appendix
Table 7: Ablation studies for static feature importance on the validation set over Europe. Shown are KGE for CAEMLS (KGE cr{}_{c}^{r} ) and GRDC (KGE gr{}_{g}^{r} ) splits on reanalysis. along with the F1 r @(1.5-500) for reanalysis over the full domain ( 281281 K points). In the second half of the table, we report the results of R2 co{}_{c}^{o} , KGE co{}_{c}^{o} , and F1 co{}_{c}^{o} @(1.5-20) on the fine-tunning on CAMELS observations.
 EFAS ReanalysisCAMELS Observations
Bi-SSDR2 cr{}_{c}^{r} | KGE cr{}_{c}^{r} | R2 gr{}_{g}^{r} | KGE gr{}_{g}^{r} | F1 r ( ↑\uparrow )R2 co{}_{c}^{o} | KGE co{}_{c}^{o} | F1 co{}_{c}^{o} ( ↑\uparrow )
✗0.6715 | 0.7057 | 0.6307 | 0.6915 | 0.32250.6793 | 0.7324 | 0.2947
✓0.7061 | 0.7396 | 0.7043 | 0.7493 | 0.33940.6990 | 0.7472 | 0.3109
Appendix
Table 8: Ablation studies for the importance of Bi-SSD in the forecasting block. Shown are KGE for CAEMLS (KGE cr{}_{c}^{r} ) and GRDC (KGE gr{}_{g}^{r} ) splits on reanalysis. along with the F1 r @(1.5-500) for reanalysis over the full domain ( 281281 K points). We also report the results of R2 co{}_{c}^{o} , KGE co{}_{c}^{o} , and F1 co{}_{c}^{o} @(1.5-20) on the fine-tunning on CAMELS observations.
Results on EFAS Reanalysis
# Routing modulesR2 cr{}_{c}^{r} | KGE cr{}_{c}^{r} | R2 gr{}_{g}^{r} | KGE gr{}_{g}^{r} | F1 r ( ↑\uparrow )
×1\times 10.7061 | 0.7396 | 0.7043 | 0.7493 | 0.3394
×2\times 20.6547 | 0.7057 | 0.6289 | 0.6650 | 0.2839
Appendix
Table 9: Ablation studies for the number of routing module/distributed steps in the forecasting block. Shown are KGE for CAEMLS (KGE cr{}_{c}^{r} ) and GRDC (KGE gr{}_{g}^{r} ) splits on reanalysis. along with the F1 r @(1.5-500) for reanalysis over the full domain ( 281281 K points).
 EFAS ReanalysisCAMELS Observations
TypeR2 cr{}_{c}^{r} | KGE cr{}_{c}^{r} | R2 gr{}_{g}^{r} | KGE gr{}_{g}^{r} | F1 r ( ↑\uparrow )R2 co{}_{c}^{o} | KGE co{}_{c}^{o} | F1 co{}_{c}^{o} ( ↑\uparrow )
Connection per Neighbor0.7027 | 0.7264 | 0.7064 | 0.7456 | 0.29500.6936 | 0.7318 | 0.2461
Binary Mask0.6622 | 0.6816 | 0.6404 | 0.6751 | 0.11220.6245 | 0.7158 | 0.2954
All ++ Sort(area)0.7061 | 0.7396 | 0.7043 | 0.7493 | 0.33940.6990 | 0.7472 | 0.3109
Appendix
Table 10: Ablation studies for different upstream aggregation techniques as a supplementary to Table 2 (c). Shown are KGE for CAEMLS (KGE cr{}_{c}^{r} ) and GRDC (KGE gr{}_{g}^{r} ) splits on reanalysis. along with the F1 r @(1.5-500) for reanalysis over the full domain ( 281281 K points). We also report the results of R2 co{}_{c}^{o} , KGE co{}_{c}^{o} , and F1 co{}_{c}^{o} @(1.5-20) on the fine-tunning on CAMELS observations.
Figure 7: Comparison between EFAS reanalysis (top, used as a reference), a model without routing (middle row), a model with routing (bottom row) during an extreme flood in Saarland Germany , mid May 2024. Shown are flood severity maps at 1 arcminute resolution where each panel shows flood extent at different lead times from nowcasting to 4 days lead time. Without explicit distributed modeling, the routing-free model produces heterogeneous predictions in which neighboring pixels can have inconsistent discharge, leading to many false positives. Our routing module instead yields more homogeneous and spatially clustered predictions. In addition, the routing module makes the predicted river discharge consistent with topology across the river skeletons (Sec. B ).
Results on EFAS Reanalysis
Bi-directionalityTraining time* ( ↓\downarrow )Inference time* ( ↓\downarrow )
 1 epoch 16 steps280,820 points
×2\times 2 Mamba (RiverMamba)∼\sim 159 min6.20 ±\pm 0.005 sec
×1\times 1 Mamba2 (Ours)∼\sim 49 min4.10 ±\pm 0.008 sec
∗ tested on NVIDIA A40 GPU.
Appendix
Table 11: Ablation studies for bi-directional Mamba in comparison to RiverMamba. Shown are training time for in epoch in minutes and inference time in seconds.
 OverallLead Time tt
SourceMean24 h24\,\mathrm{h}48 h48\,\mathrm{h}72 h72\,\mathrm{h}96 h96\,\mathrm{h}
EFAS Reanalysis (GT)99.97%----
✗ routing99.09%99.19%99.34%99.37%99.35%
✓ routing99.62%99.73%99.87%99.89%99.89%
Δ\Delta Benefit of adding routing+0.53%+0.54%+0.53%+0.52%+0.54%
Appendix
Table 12: LDD consistency rate (ratio of pixel pairs satisfying downstream ≥\geq upstream) for the EFAS reanalysis as a reference and the two forecast models. Evaluation is averaged for the years 20242024 - 20252025 over the Mid European domain for 1616 time steps and 11 nowcasting step.
Figure 8: Local Drain Direction (LDD) consistency check (downstream ≥\geq upstream) for the forecast issued on 2024-02-09T00:00:00 UTC, at four lead times ( 2424 , 4848 , 7272 , and 96 h96\,\mathrm{h} . The magenta (GT) or red (model) markers flag pixels where the LDD upstream point carries more discharge than its downstream neighbour.
Figure 9: Upstream drainage area distribution over the Mid European domain with 280,820280,820 EFAS grid points and the 66 classification classes we used in Sec. C , Table 13 , and Fig. 10 . On the left side is the histogram of upstream area. On the right is the spatial map of the assigned class at every grid point starting from small headwater streams (class 1) to the largest rivers (class 6).
Class (Uparea km 2 )MAE ( ↓\downarrow )RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-500) ( ↑\uparrow )
1: 1.95-3.040.0143 →\rightarrow 0.00780.0213 →\rightarrow 0.01330.016 →\rightarrow 0.4170.236 →\rightarrow 0.5320.1218 →\rightarrow 0.1279
2: 3.04-5.280.0180 →\rightarrow 0.01440.0306 →\rightarrow 0.02440.287 →\rightarrow 0.4470.429 →\rightarrow 0.5540.0687 →\rightarrow 0.1249
3: 5.28-12.880.0349 →\rightarrow 0.02570.0589 →\rightarrow 0.04370.243 →\rightarrow 0.4760.399 →\rightarrow 0.5810.0471 →\rightarrow 0.1215
4: 12.88-84.870.0960 →\rightarrow 0.08600.1681 →\rightarrow 0.15420.422 →\rightarrow 0.5210.540 →\rightarrow 0.6260.0578 →\rightarrow 0.1071
5: 84.87-559.260.5460 →\rightarrow 0.53271.0069 →\rightarrow 0.97180.559 →\rightarrow 0.5880.630 →\rightarrow 0.6710.0503 →\rightarrow 0.0851
6: 559.26-1600027.9160 →\rightarrow 7.930014.408 →\rightarrow 14.0030.557 →\rightarrow 0.6170.641 →\rightarrow 0.6930.0410 →\rightarrow 0.0840
Appendix
Table 13: Gauged performance of the model with routing against a model without it by upstream-area class. Evaluation is done over 280,820280,820 EFAS grid points for the validation and test sets combined ( 20222022 - 20252025 ). For each class and metric, values are the class-median of the per-pixel MAE/RMSE/R2/KGE and averaged across the 1717 lead times. Cyan marks cells where routing improves the metric and Purple marks cells where it is worse.
Figure 10: Per pixel percentage improvement from the model with routing over a model without routing, for RMSE r , R2 r , KGE r and F1 r @(1.5-500). For each grid point, the percentage improvement is computed as ΔRMSE=RMSE(w/o routing)−RMSE(routing)RMSE(w/o routing)×100, ΔR2=R2(routing)−R2(w/o routing)∣R2(w/o routing)∣×100, ΔKGE=KGE(routing)−KGE(w/o routing)∣KGE(w/o routing)∣×100 ,\Delta{\text{RMSE}}=\frac{\text{RMSE}_{\text{(w/o routing)}}-\text{RMSE}_{\text{(routing)}}}{\text{RMSE}_{\text{(w/o routing)}}}\times 100,~\Delta{\text{R2}}=\frac{\text{R2}_{\text{(routing)}}-\text{R2}_{\text{(w/o routing)}}}{|\text{R2}_{\text{(w/o routing)}}|}\times 100,~\Delta{\text{KGE}}=\frac{\text{KGE}_{\text{(routing)}}-\text{KGE}_{\text{(w/o routing)}}}{|\text{KGE}_{\text{(w/o routing)}}|}\times 100\,, where subscripts "routing" and "w/o routing" denote a model with routing and without it. F1 is computed per pixel and the plotted quantity is the plain difference ΔF1=F1(routing)−F1(w/o routing).\Delta{\text{F1}}=\text{F1}_{\text{(routing)}}-\text{F1}_{\text{(w/o routing)}}. The values are bounded in [−1,1][-1,1] . Positive (green) values indicate routing improves the metric at that pixel and negative (magenta) values indicate it is worse.
Figure 11: The 44 spatial-holdout regions used to construct the ungauged evaluation in Table 14 and Fig. 12 . The Mid European domain is split into 44 equal quadrants at its geometric midpoint (51 ∘ N, 10.5 ∘ E). Region 11 (top-left), Region 22 (top-right), Region 33 (bottom-left) and Region 44 (bottom-right).
Class (Uparea km 2 )MAE ( ↓\downarrow )RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-500) ( ↑\uparrow )
1: 1.95-3.040.0346 →\rightarrow 0.01260.0536 →\rightarrow 0.0183-6.758 →\rightarrow 0.089-1.159 →\rightarrow 0.2370.0539 →\rightarrow 0.0543
2: 3.04-5.280.0502 →\rightarrow 0.02280.0835 →\rightarrow 0.0338-2.085 →\rightarrow 0.105-0.349 →\rightarrow 0.2680.0391 →\rightarrow 0.0519
3: 5.28-12.880.0843 →\rightarrow 0.04090.1348 →\rightarrow 0.0616-1.034 →\rightarrow 0.109-0.098 →\rightarrow 0.2910.0334 →\rightarrow 0.0531
4: 12.88-84.870.2624 →\rightarrow 0.14690.4032 →\rightarrow 0.2264-0.633 →\rightarrow 0.119-0.031 →\rightarrow 0.3000.0288 →\rightarrow 0.0454
5: 84.87-559.261.4830 →\rightarrow 0.97202.1890 →\rightarrow 1.4710-0.602 →\rightarrow 0.143-0.044 →\rightarrow 0.3190.0266 →\rightarrow 0.0390
6: 559.26-16000220.564 →\rightarrow 14.65828.943 →\rightarrow 22.410-0.918 →\rightarrow -0.144-0.036 →\rightarrow 0.2390.0186 →\rightarrow 0.0223
Appendix
Table 14: Ungauged (spatial-holdout) performance of routing against without routing by upstream-area class. The evaluation is at all 280,820280,820 EFAS grid points over the years 20222022 - 20252025 . For each class and metric, values are the class-median of the per pixel metric and averaged across the 1717 lead times. Cyan marks points where routing improves the metric and Purple marks points where it is worse.
Figure 12: Per pixel percentage improvement for routing over without routing in the ungauged (spatial-holdout) setting. Metrics RMSE r , R2 r , KGE r , and F1 r @(1.5-500) are averaged over lead time for the years 20222022 - 20252025 . For each grid point, the percentage improvement is computed as ΔRMSE=RMSE(w/o routing)−RMSE(routing)RMSE(w/o routing)×100, ΔR2=R2(routing)−R2(w/o routing)∣R2(w/o routing)∣×100, ΔKGE=KGE(routing)−KGE(w/o routing)∣KGE(w/o routing)∣×100 ,\Delta{\text{RMSE}}=\frac{\text{RMSE}_{\text{(w/o routing)}}-\text{RMSE}_{\text{(routing)}}}{\text{RMSE}_{\text{(w/o routing)}}}\times 100,~\Delta{\text{R2}}=\frac{\text{R2}_{\text{(routing)}}-\text{R2}_{\text{(w/o routing)}}}{|\text{R2}_{\text{(w/o routing)}}|}\times 100,~\Delta{\text{KGE}}=\frac{\text{KGE}_{\text{(routing)}}-\text{KGE}_{\text{(w/o routing)}}}{|\text{KGE}_{\text{(w/o routing)}}|}\times 100\,, where subscripts "routing" and "w/o routing" denote a model with routing and without it. F1 is computed per pixel and the plotted quantity is the plain difference ΔF1=F1(routing)−F1(w/o routing).\Delta{\text{F1}}=\text{F1}_{\text{(routing)}}-\text{F1}_{\text{(w/o routing)}}. The F1 values are bounded in [−1,1][-1,1] . Positive (green) values indicate routing improves the metric at that pixel and negative (magenta) values indicate it is worse.
Figure 13: The two random holdout splits of the 892892 CAMELS stations ( 446446 each) used for the fully ungauged evaluation.
 KGE ( ↑\uparrow )R2 ( ↑\uparrow )F1@(1.5-500) ( ↑\uparrow )
mean0.5630.3160.185
median0.6020.5810.159
Appendix
Table 15: Per station skill of the ungauged CAMELS prediction against observations.
Figure 14: Per station KGE, R2, and F1 of the ungauged CAMELS prediction. Evaluation is done for the years 20102010 to 20242024 and for 44 days lead time + 11 nowcasting.
Figure 15: The two random holdout splits of the 480480 GRDC stations ( 240240 each) used for the fully ungauged evaluation.
 KGE ( ↑\uparrow )R2 ( ↑\uparrow )F1@(1.5-500) ( ↑\uparrow )
mean0.5420.2020.200
median0.5950.5450.175
Appendix
Table 16: Per station skill of the ungauged GRDC prediction against observations.
Figure 16: Per station KGE, R2, and F1 of the ungauged GRDC prediction. Evaluation is done for the years 20102010 to 20252025 and for 44 days lead time + 11 nowcasting.
Model# Param.auto-regressiveProbabilisticNowcasting
RiverMamba4.4 million✗✗✗
Ours2.3 million✓✓✓
Appendix
Table 17: Comparison to RiverMamba [ Shams Eddin et al., 2025 ] .
   Validation (2022-2023) Test (2024-2025)
 Model R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-20) ( ↑\uparrow ) R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-20) ( ↑\uparrow )
DetermRiverMamba 0.7075 ±\pm 0.01380.7255 ±\pm 0.00770.1242 ±\pm 0.0160 0.5806 ±\pm 0.02210.6580 ±\pm 0.01180.0611 ±\pm 0.0104
Ours 0.7075 ±\pm 0.00560.7410 ±\pm 0.01340.1627 ±\pm 0.0143 0.5676 ±\pm 0.00970.6589 ±\pm 0.01620.0841 ±\pm 0.0037
ProbOurs (mean) 0.72470.76950.2129 0.62800.72840.2266
Appendix
Table 18: Results on EFAS Reanalysis (CAMELS split). ( ±\pm ) denotes the standard deviation for 3 runs.
   Validation (2022-2023) Test (2024-2025)
 Model R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-20) ( ↑\uparrow ) R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-20) ( ↑\uparrow )
DetermRiverMamba 0.7326 ±\pm 0.00990.7579 ±\pm 0.00520.1106 ±\pm 0.0125 0.6161 ±\pm 0.02220.6927 ±\pm 0.01410.0990 ±\pm 0.0151
Ours 0.6957 ±\pm 0.02730.7507 ±\pm 0.02730.1326 ±\pm 0.0272 0.5631 ±\pm 0.02370.6645 ±\pm 0.02410.1132 ±\pm 0.0136
ProbOurs (mean) 0.73590.78760.1998 0.65050.75520.2125
Appendix
Table 19: Results on EFAS Reanalysis (GRDC split). ( ±\pm ) denotes the standard deviation for 3 runs.
   Validation (2022-2023) Test (2024-2025)
 Model F1@(1.5-20) ( ↑\uparrow )F1@(1.5-500) ( ↑\uparrow ) F1@(1.5-20) ( ↑\uparrow )F1@(1.5-500) ( ↑\uparrow )
DetermRiverMamba 0.1670 ±\pm 0.01220.1407 ±\pm 0.0099 0.1713 ±\pm 0.01340.1379 ±\pm 0.0108
Ours 0.3542 ±\pm 0.01990.3337 ±\pm 0.0176 0.3183 ±\pm 0.01910.3143 ±\pm 0.0129
Appendix
Table 20: Results on EFAS Reanalysis (Full mid Europe domain). ( ±\pm ) denotes the standard deviation for 3 runs.
Figure 17: KGE of the river discharge forecasting with different lead time on CAMELS reanalysis split (test set 20242024 - 20252025 temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 18: Averaged F1-score for RPs 1.5-20 with lead time on CAMELS reanalysis split (test set 20242024 - 20252025 temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 19: Cumulative distribution function (CDF) of station-level KGE on the CAMELS reanalysis split (N=892), across lead times. Dashed vertical lines shows each model’s media.
Figure 20: Cumulative distribution function (CDF) of station-level F1@(1.5-10) on the CAMELS reanalysis split (N=892), across lead times. Dashed vertical lines shows each model’s media.
Figure 21: Comparison of KGE between Ours (mean) and RiverMamba on CAMELS EFAS reanalysis split across lead times (test set 20242024 - 20252025 temporally out-of-sample).
Figure 22: Comparison of F1-score between Ours (mean) and RiverMamba on CAMELS reanalysis EFAS split for the 1.5-year return period events across lead times (test set 20242024 - 20252025 temporally out-of-sample).
   Validation (2022-2023) Test (2024)
 Model R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-20) ( ↑\uparrow ) R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-20) ( ↑\uparrow )
DetermRiverMamba 0.7236 ±\pm 0.00000.7563 ±\pm 0.01130.3181 ±\pm 0.0041 0.5428 ±\pm 0.01510.6257 ±\pm 0.02350.2864 ±\pm 0.0675
Ours 0.7286 ±\pm 0.00110.7716 ±\pm 0.00650.2883 ±\pm 0.0008 0.6216 ±\pm 0.03130.6761 ±\pm 0.02560.2805 ±\pm 0.0107
ProbOurs (mean) 0.71750.77440.3151 0.61340.69410.3580
Appendix
Table 21: Results on gauged CAMELS Observations. ( ±\pm ) denotes the standard deviation for 3 runs.
Figure 23: KGE of the river discharge forecasting with different lead time on CAMELS observations (test set 20242024 temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 24: Averaged F1-score for 1.5-20 year return periods with lead time on CAMELS observations (test set 20242024 temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
   Test set (2024-05-30 to 2025-12-21)
 Model RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )F1@(1.5-20) ( ↑\uparrow )F1@(1.5-500) ( ↑\uparrow )
DetermEFAS (control run) 9.83940.04810.39670.09450.0525
Ours 7.2889 ±\pm 0.16280.3413 ±\pm 0.01350.5180 ±\pm 0.00730.0703 ±\pm 0.01590.0512 ±\pm 0.0175
ProbOurs (mean) 7.28600.36730.55150.10150.1515
Appendix
Table 22: Results on gauged GRDC Observations. ( ±\pm ) denotes the standard deviation for 3 runs.
Figure 25: KGE of the river discharge forecasting with different lead time on GRDC observations (test set (2024-05-30 to 2025-12-21) temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 26: Cumulative distribution function (CDF) of station-level KGE on the GRDC observations, across lead times. Dashed vertical lines shows each model’s media.
   Test set (2024-06-01 to 2024-08-31) 
 Model RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )CRPS ( ↓\downarrow )CRPS mae ( ↓\downarrow )CRPS spread ( ↑\uparrow )
ProbEFAS (mean of 51 ens members) 6.8122-0.85170.16673.88625.40431.5180
Ours (mean of 12 ens members) 4.56270.05190.37182.15994.12801.9682
Appendix
Table 23: Results on gauged GRDC Observations.
Figure 27: CRPS of the river discharge forecasting with different lead time on GRDC observations (test set (2024-06-01 to 2024-08-31) temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 28: Ensemble spread of the river discharge forecasting with different lead time on GRDC observations (test set (2024-06-01 to 2024-08-31) temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 29: Cumulative distribution function (CDF) of station-level CRPS on the GRDC observations, across lead times. Dashed vertical lines shows each model’s media.
Figure 30: Ensemble forecast against EFAS reanalysis at Giessen (GRDC 6335450) on the Argen River ( 47.62∘47.62^{\circ} N, 9.61∘9.61^{\circ} E, catchment area 673673 km 2 ), for lead time 2424 h. For forecast, the ensemble median (solid line) is shown with nested uncertainty bands spanning the 2525 th– 7575 th percentile (darkest), 1010 th– 9090 th percentile, and min–max range (lightest) across all 1818 ensemble members, against the EFAS reanalysis (black). Dashed olive lines are the EFAS flood return period thresholds ( 1.51.5 - 20.020.0 years).
Figure 31: Ensemble forecast against EFAS reanalysis at Buetzow (GRDC 6341400) on the Warnow River ( 53.84∘53.84^{\circ} N, 11.99∘11.99^{\circ} E, catchment area 1,4461,446 km 2 ), for lead time 4848 h. For forecast, the ensemble median (solid line) is shown with nested uncertainty bands spanning the 2525 th– 7575 th percentile (darkest), 1010 th– 9090 th percentile, and min–max range (lightest) across all 1818 ensemble members, against the EFAS reanalysis (black). Dashed olive lines are the EFAS flood return period thresholds ( 1.51.5 - 20.020.0 years).
Figure 32: Ensemble forecast against EFAS reanalysis at Achleiten (GRDC 6342900) on the Danube River ( 48.59∘48.59^{\circ} N, 13.51∘13.51^{\circ} E, catchment area 76,65276,652 km 2 ), for lead time 120120 h. For forecast, the ensemble median (solid line) is shown with nested uncertainty bands spanning the 2525 th– 7575 th percentile (darkest), 1010 th– 9090 th percentile, and min–max range (lightest) across all 1818 ensemble members, against the EFAS reanalysis (black). Dashed olive lines are the EFAS flood return period thresholds ( 1.51.5 - 20.020.0 years).
Figure 33: Ensemble forecasts against EFAS reanalysis at Guben (GRDC 6357501) on the Lech River ( 51.98∘51.98^{\circ} N, 14.71∘14.71^{\circ} E, catchment area 1,1841,184 km 2 ), for five lead times ( 4848 h, 9696 h, 120120 h, 144144 h, 168168 ). Each row shows the EFAS reanalysis discharge (black) against individual ensemble members (colored lines) for that lead time. Dashed olive lines are the EFAS flood return period thresholds ( 1.51.5 - 20.020.0 years).
Figure 34: Ensemble forecast against CAMELS-CZ in Czechia for the station U420730 on the Ohře River ( 50.17∘50.17^{\circ} N, 12.62∘12.62^{\circ} E, catchment area 1,7141,714 km 2 ), for lead time 4848 h. For forecast, the ensemble median (solid line) is shown with nested uncertainty bands spanning the 2525 th– 7575 th percentile (darkest), 1010 th– 9090 th percentile, and min–max range (lightest) across all 1212 ensemble members, against the CAMELS-CZ observations (black). Dashed olive lines are the flood return period thresholds ( 1.51.5 - 20.020.0 years) computed from observations.
Figure 35: Ensemble forecast against CAMELS-FI in Finland for the station Kärnäjärvi-Kellankoski (FI_980) at the lake Kymijoki ( 63.14∘63.14^{\circ} N, 25.94∘25.94^{\circ} E, catchment area 1,5831,583 km 2 ), for lead time 9696 h. For forecast, the ensemble median (solid line) is shown with nested uncertainty bands spanning the 2525 th– 7575 th percentile (darkest), 1010 th– 9090 th percentile, and min–max range (lightest) across all 1212 ensemble members, against the CAMELS-FI observations (black). Dashed olive lines are the flood return period thresholds ( 1.51.5 - 20.020.0 years) computed from observations.
Figure 36: Ensemble forecasts against CAMELS-FI in Finland for the station Kärnäjärvi-Kellankoski (FI_980) at the lake Kymijoki ( 63.14∘63.14^{\circ} N, 25.94∘25.94^{\circ} E, catchment area 1,5831,583 km 2 ), for five lead times ( 4848 h, 9696 h, 120120 h, 144144 h, 168168 ). Each row shows the CAMELS-CZ observational discharge (black) against individual ensemble members (colored lines) for that lead time. Dashed olive lines are the flood return period thresholds ( 1.51.5 - 20.020.0 years) computed from observations.
Figure 37: Ensemble forecast against GRDC observations at Castle Roberts (GRDC 6503120) on the Maigue River ( 52.56∘52.56^{\circ} N, 8.78∘8.78^{\circ} W, catchment area 830830 km 2 ), for lead time 2424 h. For forecast, the ensemble median (solid line) is shown with nested uncertainty bands spanning the 2525 th– 7575 th percentile (darkest), 1010 th– 9090 th percentile, and min–max range (lightest) across all 1818 ensemble members, against the EFAS reanalysis (black). Dashed olive lines are the flood return period thresholds ( 1.51.5 - 20.020.0 years) computed from observations.
Figure 38: Ensemble forecasts against GRDC observations at Castle Roberts (GRDC 6503120) on the Maigue River ( 52.56∘52.56^{\circ} N, 8.78∘8.78^{\circ} W, catchment area 830830 km 2 ), for five lead times ( 4848 h, 9696 h, 120120 h, 144144 h, 168168 ). Each row shows the GRDC observational discharge (black) against individual ensemble members (colored lines) for that lead time. Dashed olive lines are the flood return period thresholds ( 1.51.5 - 20.020.0 years) computed from observations.
Figure 39: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the Central European flood, September 2024. Shown are flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days.
Figure 40: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the Central European flood, September 2024. First row shows flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days. Rows 22 to 66 show the exceedance probabilities for flood return periods ( 1.51.5 - 20.020.0 ) computed from 1212 ensemble members.
Figure 41: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the Spanish flood, October 2024. Shown are flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days.
Figure 42: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the Spanish flood, October 2024. First row shows flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days. Rows 22 to 66 show the exceedance probabilities for flood return periods ( 1.51.5 - 20.020.0 ) computed from 1212 ensemble members.
Figure 43: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the United Kingdom flood, January 2025. Shown are flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days.
Figure 44: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the United Kingdom flood, January 2025. First row shows flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days. Rows 22 to 66 show the exceedance probabilities for flood return periods ( 1.51.5 - 20.020.0 ) computed from 1212 ensemble members.
Figure 45: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the Balkan Rivers flood, November 2025. Shown are flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days.
Figure 46: Comparison between EFAS reanalysis (first row, used as a reference) and our ensemble forecast (second to sixth rows) during the Balkan Rivers flood, November 2025. First row shows flood severity maps at 11 -arcminute resolution where each panel shows flood extent at different lead times from 11 to 55 days. Rows 22 to 66 show the exceedance probabilities for flood return periods ( 1.51.5 - 20.020.0 ) computed from 1212 ensemble members.
Figure 47: Routing state at a fixed 48 48\, h lead time across four independent forecasts, one randomly sampled per year from 2022 to 2025. Top row shows the L 22 -norm of the 192192 -dimensional routing feature vector across the years.
   Validation (2022-2023) Test (2024-2025)
 Model RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow ) RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )
DeterministicPersistence 4.9490.0860.533 6.013-0.1180.456
ED-LSTM 4.622 ±\pm 0.0160.229 ±\pm 0.0070.300 ±\pm 0.009 5.316 ±\pm 0.0100.157 ±\pm 0.0070.264 ±\pm 0.013
ME-LSTM 4.830 ±\pm 0.0340.181 ±\pm 0.0110.286 ±\pm 0.018 5.584 ±\pm 0.0420.095 ±\pm 0.0140.235 ±\pm 0.020
Ours 3.972 ±\pm 0.0520.428 ±\pm 0.0150.534 ±\pm 0.011 4.747 ±\pm 0.0400.333 ±\pm 0.0110.481 ±\pm 0.008
ProbabilisticClima (median) 5.382-0.005-0.035 6.025-0.068-0.074
Appendix
Table 24: Results on EFAS Reanalysis (CAMELS split). ( ±\pm ) denotes the standard deviation for 3 runs.
Figure 48: R2 of the river discharge forecasting with different lead time on reanalysis CAMELS split (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 49: Mean F1@(1.5-20) of the river discharge forecasting with different lead time on reanalysis CAMELS split (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 50: Cumulative distribution function (CDF) of station-level F1@(1.5-20) on the CAMELS reanalysis split (N=2206), across lead times. Dashed vertical lines shows each model’s media.
Figure 51: Cumulative distribution function (CDF) of station-level KGE on the CAMELS reanalysis split (N=2206), across lead times. Dashed vertical lines shows each model’s media.
Figure 52: Cumulative distribution function (CDF) of station-level R2 on the CAMELS reanalysis split (N=2206), across lead times. Dashed vertical lines shows each model’s media.
   Validation (2022-2023) Test (2024-2025)
 Model RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow ) RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )
DeterministicPersistence 10.4170.0960.542 12.056-0.0760.468
ED-LSTM 9.739 ±\pm 0.0650.239 ±\pm 0.0090.318 ±\pm 0.009 10.975 ±\pm 0.0050.151 ±\pm 0.0040.269 ±\pm 0.013
ME-LSTM 10.242 ±\pm 0.0670.186 ±\pm 0.0100.321 ±\pm 0.020 11.515 ±\pm 0.0620.085 ±\pm 0.0130.249 ±\pm 0.021
Ours 8.494 ±\pm 0.1380.420 ±\pm 0.0190.531 ±\pm 0.014 9.851 ±\pm 0.0740.314 ±\pm 0.0100.474 ±\pm 0.006
ProbabilisticClima (median) 11.0140.002-0.019 12.169-0.062-0.058
Appendix
Table 25: Results on EFAS Reanalysis (GRDC split). ( ±\pm ) denotes the standard deviation for 3 runs.
Results on the Validation set (2022-2023)
  ClimaDRUMED-LSTMME-LSTMOurs
CRPS ( ↓\downarrow ) 4.33284.51044.13793.97763.1837
Results on the Test set (2024-2025)
  ClimaDRUMED-LSTMME-LSTMOurs
CRPS ( ↓\downarrow ) 4.66674.93764.63844.64923.7604
Appendix
Table 26: Results for probabilistic models on EFAS Reanalysis (GRDC split).
Figure 53: CRPS of the river discharge forecasting with different lead time on reanalysis GRDC split (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 54: KGE of the river discharge forecasting with different lead time on reanalysis GRDC split (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 55: R2 of the river discharge forecasting with different lead time on reanalysis GRDC split (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 56: Mean F1@(1.5-20) of the river discharge forecasting with different lead time on reanalysis GDRC split (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 57: Cumulative distribution function (CDF) of station-level F1@(1.5-20) on the GRDC reanalysis split (N=1971), across lead times. Dashed vertical lines shows each model’s media.
Figure 58: Cumulative distribution function (CDF) of station-level KGE on the GRDC reanalysis split (N=1971), across lead times. Dashed vertical lines shows each model’s media.
Figure 59: Cumulative distribution function (CDF) of station-level R2 on the GRDC reanalysis split (N=1971), across lead times. Dashed vertical lines shows each model’s media.
   Validation (2022-2023) Test (2024-2025)
 Model RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow ) RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )
DeterministicED-LSTM 6.0542 ±\pm 0.00820.2705 ±\pm 0.00120.3491 ±\pm 0.0214 11.3616 ±\pm 0.16320.2366 ±\pm 0.01640.4501 ±\pm 0.0097
ME-LSTM 5.8810 ±\pm 0.02060.2993 ±\pm 0.00440.3550 ±\pm 0.0083 11,4754 ±\pm 0.16830.2384 ±\pm 0.02460.4330 ±\pm 0.0138
Ours 4.6795 ±\pm 0.06810.5603 ±\pm 0.00950.5890 ±\pm 0.0038 9.0163 ±\pm 0.09520.4576 ±\pm 0.00660.5839 ±\pm 0.0076
ProbabilisticDRUM (mean) 7.02360.15140.2177 12.63580.11530.3059
ED-LSTM (mean) 6.42170.20040.2780 11.56290.23080.4266
Appendix
Table 27: Results on GRDC observations. ( ±\pm ) denotes the standard deviation for 3 runs.
Results on the Validation set (2022-2023)
  DRUMED-LSTMME-LSTMOurs
CRPS ( ↓\downarrow ) 4.53753.90243.49282.4592
Results on the Test set (2024-2025)
  DRUMED-LSTMME-LSTMOurs
CRPS ( ↓\downarrow ) 5.69345.00915.03613.5194
Appendix
Table 28: Results for probabilistic models on GRDC observations.
Figure 60: CRPS of the river discharge forecasting with different lead time on GRDC observations (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 61: KGE of the river discharge forecasting with different lead time on GRDC observations (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 62: R2 of the river discharge forecasting with different lead time on GRDC observations (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 63: mean F1@(1.5-20) of the river discharge forecasting with different lead time on GRDC observations (test set 2024-2025, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 64: Cumulative distribution function (CDF) of station-level F1@(1.5-20) on GRDC observations (N=758), across lead times. Dashed vertical lines shows each model’s media.
Figure 65: Cumulative distribution function (CDF) of station-level KGE on GRDC observations (N=758), across lead times. Dashed vertical lines shows each model’s media.
Figure 66: Cumulative distribution function (CDF) of station-level R2 on GRDC observations (N=758), across lead times. Dashed vertical lines shows each model’s media.
   Validation (2022-2023) Test (2024)
 Model RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow ) RMSE ( ↓\downarrow )R2 ( ↑\uparrow )KGE ( ↑\uparrow )
DeterministicED-LSTM 6.0542 ±\pm 0.00820.2705 ±\pm 0.00120.3401 ±\pm 0.0214 6.8658 ±\pm 0.13260.0544 ±\pm 0.00720.0906 ±\pm 0.0136
ME-LSTM 5.8810 ±\pm 0.02060.2993 ±\pm 0.00440.3550 ±\pm 0.0083 6.9621 ±\pm 0.08740.0500 ±\pm 0.01000.1186 ±\pm 0.0123
Ours 4.6795 ±\pm 0.06810.5603 ±\pm 0.00950.5890 ±\pm 0.0038 5.8406 ±\pm 0.13110.3485 ±\pm 0.02020.4162 ±\pm 0.0284
ProbabilisticDRUM (mean) 7.02360.15140.2177 7.34330.01150.0576
ED-LSTM (mean) 6.42170.20040.2780 7.08960.04520.0831
Appendix
Table 29: Results on CAMELS observations. ( ±\pm ) denotes the standard deviation for 3 runs.
Results on the Validation set (2022-2023)
  DRUMED-LSTMME-LSTMOurs
CRPS ( ↓\downarrow ) 2.72392.31102.12371.4840
Results on the Test set (2024)
  DRUMED-LSTMME-LSTMOurs
CRPS ( ↓\downarrow ) 2.36852.26832.26831.5678
Appendix
Table 30: Results for probabilistic models on CAMELS observations.
Figure 67: CRPS of the river discharge forecasting with different lead time on CAMELS observations (val set 2022-2023, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 68: KGE of the river discharge forecasting with different lead time on CAMELS observations (val set 2022-2023, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 69: R2 of the river discharge forecasting with different lead time on CAMELS observations (val set 2022-2023, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 70: Mean F1@(1.5-20) of the river discharge forecasting with different lead time on CAMELS observations (val set 2022-2023, temporally out-of-sample). Distribution quartiles are displayed in boxes, and the entire range excluding outliers is displayed in whiskers. The median score for the model is shown by the cyan line in the box.
Figure 71: Cumulative distribution function (CDF) of station-level F1@(1.5-20) on CAMELS observations (N=715), across lead times. Dashed vertical lines shows each model’s media.
Figure 72: Cumulative distribution function (CDF) of station-level KGE on CAMELS observations (N=715), across lead times. Dashed vertical lines shows each model’s media.
Figure 73: Cumulative distribution function (CDF) of station-level R2 on CAMELS observations (N=715), across lead times. Dashed vertical lines shows each model’s media.
Figure 74: Overview of the input dynamic data to the model. We use 44 daily time steps as input for the initial condition ( T=4T=4 ). For example, CPC data has 11 day latency ( t−4t-4 ) so we take 4 data samples at t−4t-4 , t−8t-8 , t−12t-12 , and t−16t-16 , where tt is the 6-hourly time step. Note that we ignore here the 44 hours latency for IMERG since we use the daily mean. In an operational setting, step tt daily averaged for IMERG can be obtained using the mean over 00:00 UTC →∼\rightarrow\sim 20:00 UTC.
Figure 75: Timeframe and coverage of input dynamic datasets and target streamflow data. Unavailable sources are replaced with ERA5-Land.
ConfigurationPre-trainingFine-tuning
(EFAS Reanalysis)(GRDC/CAMELS)
OptimizerAdamWAdamW
Learning rate0.00030.0001
Minimum learning rate0.00010.0001
Batch size ( BB )11
Learning rate schedulerCosine annealingCosine annealing
Appendix
Table 31: Implementation details.
Figure 76: A simplified view of splitting along Gilbert space-filling curve over the EFAS domain. Each split has 200,000200,000 points. The model runs over these splits and the final results are stitch together to form the prediction over EFAS domain with 5.25.2 million points.
10K20K40K80K160K300K600K
0.01540.3010.5911.1722.3474.3958.913
Appendix
Table 32: Inference time in seconds.
Figure 77: A spatio-temporal scans along 44 frames in the encoder block using Gilbert serialization.
Figure 78: Overview of the EFAS domain.
Figure 79: Overview of the Mid European domain as a subset of EFAS domain. We used this domain for the ablation studies.
Figure 80: Locations of the 1,9711{,}971 selected GRDC stations used to train and evaluate the model. Color encodes the KGE of the EFAS reanalysis discharge against the GRDC observations at each station.
Figure 81: Hydrograph for EFAS reanalysis (purple line) from 1 January 2000 to 31 December 2022 and observations (gold line), for a gauging station on the Saale River (Rudolstadt, GRDC-6340302). The top-right box displays summary statistics from the reanalysis’s evaluation against the observations.
DatasetSourceCountryStartEnd# Stations
CAMELS-CHHöge et al. [2023]Switzerland19812020136
CAMELS-CZInstitute et al. [2025]Czechia19782024173
CAMELS-DELoritz et al. [2024]Germany19512020605
CAMELS-DKLiu et al. [2025]Denmark1989201954
CAMELS-ESCasado-Rodríguez et al. [2026]Spain19912020244
CAMELS-FISeppä et al. [2026]Finland19612023234
Appendix
Table 33: CAMELS datasets used in this work. Shown are the source of data, observation period, and the number of stations retained after matching to the EFAS drainage network.
Figure 82: Locations of the 2,2402{,}240 selected CAMELS stations used to train and evaluate the model. Color encodes the KGE of the EFAS reanalysis discharge against the CAMELS observations at each station.
Figure 83: Hydrograph for EFAS reanalysis (purple line) from 1 January 2015 to 31 December 2024 and observations (gold line), for a gauging station in Switzerland (CAMELS-CH-2063). The top-right box displays summary statistics from the reanalysis’s evaluation against the observations.
VariableLong nameUnitHeightSurface parameters
d2m2m dewpoint temperatureK2minstantaneous
etotal evaporationm of water equivalentsurfaceaccumulated
rorunoffmsurface & subsurfaceaccumulated
sdsnow depthm of water equivalentsurfaceaccumulated
sfsnowfallm of water equivalentsurfaceaccumulated
smltsnowmeltm of water equivalentsurfaceaccumulated
Appendix
Table 34: Details about the processed variables from the ECMWF Integrated Forecast System (IFS) High Resolution (HRES) atmospheric model (Sec. O.3 ) and ERA5-Land reanalysis (Sec. O.5 ). IFS-HRES accumulated state variables start at 0000 : 0000 UTC of the forecast time and ends at the end of the forecast. ERA5-Land accumulation starts at 0000 : 0000 UTC previous day and ends at 0000 : 0000 UTC current day.
VariableLong nameUnitHeightSurface parameters
pdwater vapour pressurehPa-instantaneous
prtotal precipitationmmsurfaceaccumulated
rgdownward surface solar radiationJ/m 2surfaceaccumulated
tnminimum air temperatureC ∘2minstantaneous
txmaximum air temperatureC ∘2minstantaneous
wswind speedm/s10minstantaneous
Appendix
Table 35: Details about the processed variables from EMO data [ Thiemig et al., 2022 ] .
Figure 84: AlphaEarth Embedding over Germany for the years 20172017 , 20212021 , and 20242024 . Shown are Bands; A 0101 , A 0909 , and A 1616 as an RGB image.
Figure 85: AlphaEarth Embedding over EFAS domain. Shown are Bands; A 0101 , A 0909 , and A 1616 as an RGB image.

Explore similar work

CardsList
Jul 3, 2026cs.LG

Joint distribution of upstream runoff governs downstream river-discharge prediction uncertainty in distributed ML models

Uncertainty quantification of hydrological predictions is necessary to inform operational decisions. Recent generative machine-learning methods have advanced probabilistic streamflow prediction, but have remained confined to lumped models that predict a basin outlet directly. At the same time, deterministic LSTM runoff models are increasingly applied at grid or catchment scale and routed through river networks to produce spatially continuous, physically consistent discharge fields. This technical note argues that moving probabilistic prediction from lumped to distributed models introduces a specific new requirement: the joint distribution of upstream runoff generation must be sampled jointly. In lumped inference, the model predicts the outlet distribution directly and can modulate spread from basin attributes. In distributed inference, downstream discharge is obtained by routing many upstream runoff predictions, so independent local sampling averages uncertainty away. Using Japan as a case study, we train two probabilistic basin-scale runoff LSTMs and route their runoff through a Hayami routing scheme. Randomly matching upstream ensemble members produces severely under-dispersed downstream ensembles, whereas a simple quantile matching strategy restores much of the spread of the direct basin-scale reference. The shift from lumped to distributed probabilistic hydrology therefore requires explicit attention to the spatial joint structure of runoff uncertainty.

Karan Ruparell, Tristan Hascoet, Takemasa Miyoshi +4

University of Reading, Department of Meteorology, Reading, UK · European Centre for Medium-Range Weather Forecasts (ECMWF), Reading, UK · RIKEN Center for Computational Sciences, Kobe, Japan +3

Uncertainty QuantificationTime Series ForecastingProbabilistic Forecasting
AskPDF
Sep 19, 2026cs.AI

Hapi: A Multivariable Land-Surface Transformer for Medium-Range Hydrological Forecasting at Continental Scale

Accurate flood forecasts several days in advance are essential for flood control, water-resource management, and emergency response. A central challenge is to produce high-resolution forecasts across continental domains where hydrological behavior varies widely from place to place. We developed Hapi, a U-Net Swin Transformer that uses fine three-dimensional patches and hierarchical shifted-window attention to forecast river discharge, surface runoff, snow water equivalent, and soil wetness index across the contiguous United States. The model produces medium-range forecasts (24--72~h) at 0.05∘0.05^{\circ} resolution with adaptive task weighting and required only 0.11 seconds for a four-variable 72-h CONUS forecast on one A100 GPU. In a held-out 2024 potential-skill evaluation with ERA5-Land inputs prescribed over the forecast horizon, Hapi achieved the highest F1-score for floods in 20 of 21 comparisons across seven GloFAS return periods and three forecast leads. Independent validation against observed daily discharge at 3{,}881 U.S. Geological Survey gauges showed that Hapi achieved the highest median Nash--Sutcliffe efficiency at every lead, supported by regional-cluster bootstrap intervals. In a matched 24-h comparison of loss formulations, adaptive task balancing produced the lowest discharge errors and the highest F1-score for floods.

Hong Zhang, John K. Hutchison, Rao Kotamarthi +7

Argonne National Laboratory, Lemont, IL 60439, USA

Time Series ForecastingSpatiotemporal ForecastingFlood Forecasting
AskPDF
Oct 1, 2026cs.LG

Physics-Refined Spatiotemporal Forecasting on Open-Boundary Hydrologic Graphs

Spatiotemporal forecasting on hydrologic graphs is especially prone to instability in open-boundary systems, where the forecast domain exchanges fluxes with an unobserved exterior. In such systems, boundary nodes receive external forcing, e.g., upstream inflows in rivers or tidal signals in coastal regions, that is typically unavailable at prediction time. The absence of this information can compound errors as forecasts unfold in an autoregressive fashion, leading to inferior long-horizon performance. This paper dissects this instability issue by exploring two questions. 1) What boundary forcing enters the forecast domain when information beyond the boundary is missing? 2) How should this forcing propagate through the domain without incurring error amplification under autoregressive rollout? To address both, we propose a new computing framework comprising two key components. First, to compensate for the boundary forcing, our framework learns ghost node proxies from the boundary and interior nodes, striving to approximate unobserved external inputs. Second, to control error accumulation from these learned proxies, we leverage two physics refiners. In particular, one refiner enforces local consistency by aligning ghost proxies with their two-hop neighbors (i.e., boundary nodes and their immediate interiors). The other refiner enhances global stability by correcting the model forecasts through a physics-guided graph neural operator, reducing long-horizon numerical drift. Two real-world hydrologic graphs are employed for empirical evaluation. Comparative results show that our proposal enjoys higher prediction accuracy and long-horizon stability over both learning-based and physics-informed model competitors.

Haoyang Jiang, Zhengui Wang, Shenghan Gao +3

William & Mary, Williamsburg, VA, USA · Virginia Institute of Marine Science, William & Mary, Gloucester Point, VA, USA · Florida Atlantic University, Boca Raton, FL, USA

Spatiotemporal ForecastingPhysics-Informed GNNsSpatiotemporal GNNs
AskPDF