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

# Hapi: A Multivariable Land-Surface Transformer for Medium-Range Hydrological Forecasting at Continental Scale
Source: [https://arxiv.org/html/2609.22702](https://arxiv.org/html/2609.22702)
John K\. Hutchison[https://orcid.org/0000-0002-7419-6611](https://orcid.org/0000-0002-7419-6611)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USARao Kotamarthi[https://orcid.org/0000-0002-2612-7590](https://orcid.org/0000-0002-2612-7590)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USAJeremy Feinstein[https://orcid.org/0009-0009-5362-6860](https://orcid.org/0009-0009-5362-6860)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USAHaiwen Guan[https://orcid.org/0009-0008-3005-1497](https://orcid.org/0009-0008-3005-1497)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USARomit Maulik[https://orcid.org/0000-0002-4370-0416](https://orcid.org/0000-0002-4370-0416)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USAVijay P\. Ramalingam[https://orcid.org/0009-0008-1599-5676](https://orcid.org/0009-0008-1599-5676)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USAJason Stock[https://orcid.org/0009-0006-2354-2373](https://orcid.org/0009-0006-2354-2373)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USATom Wall[https://orcid.org/0000-0002-4070-8623](https://orcid.org/0000-0002-4070-8623)Affiliation:Argonne National Laboratory, Lemont, IL 60439, USAAffiliation:Purdue University, 610 Purdue Mall, West Lafayette, IN 47907\*Corresponding author:[hongzhang@anl\.gov](mailto:[email protected])

###### Abstract

Accurate flood forecasts several days in advance are essential for flood control, water\-resource management, and emergency response\. Producing them at high resolution over a continental domain calls for local hydrological detail together with spatial context extending from river basins to synoptic weather systems\. We developed Hapi, a U\-Net Swin Transformer that uses fine three\-dimensional patches and hierarchical shifted\-window attention to forecast discharge, surface runoff, snow water equivalent, and soil wetness across the contiguous United States\. The model produces 24–72\-hour forecasts at0\.05∘0\.05^\{\\circ\}resolution, with learned Laplacian task weights adjusting each variable’s contribution to training\. On 2024 test data using reconstructed weather and land\-surface inputs from ERA5\-Land, Hapi outperformed an operational physics\-based model and a state\-of\-the\-art AI model in flood detection\. Independent validation against 3,881 U\.S\. Geological Survey gauges and a Hurricane Helene case study supported its advantage over the physics\-based model in reproducing daily discharge\. Controlled experiments showed that learned task weighting strengthens rare\-flood detection, which is particularly sensitive to changes in precipitation inputs\. Hapi produced a four\-variable, 72\-hour forecast across the contiguous United States with an average inference time of 0\.11 seconds on a single A100 GPU\.

Keywords:deep learning, hydrological forecasting, flood detection, land\-surface modeling, Swin Transformer, continental\-scale forecasting

## 1Introduction

Accurate flood forecasts several days in advance are essential for flood control, water\-resource management, and emergency response\([Tellman et al\., 2021](https://arxiv.org/html/2609.22702#bib.bib1);[Alfieri et al\., 2013](https://arxiv.org/html/2609.22702#bib.bib2)\)\. A central challenge is to produce these forecasts at high resolution over continental domains while representing local hydrological variation and basin\-to\-synoptic\-scale context\. Operational systems address this challenge with process\-based simulators\. The Global Flood Awareness System \(GloFAS\)\([Alfieri et al\., 2013](https://arxiv.org/html/2609.22702#bib.bib2)\), and the U\.S\. National Water Model\([Cosgrove et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib6)\), for example, solve land\-surface and river\-routing equations at continental scale and produce a full land\-surface state\. Their ensembles consume millions of CPU\-hours per year, depend on hand\-tuned regional parameters, and take months to adapt to new variables or regions\.

![Refer to caption](https://arxiv.org/html/2609.22702v1/conus_forecast_demo.png)Figure 1:Continental illustration of a Hapi discharge forecast\.a, Discharge at initialization on 20 June 2024, represented by the field carried forward in the persistence forecast\.b, Hapi’s 72\-h discharge forecast valid on 23 June\.c, GloFAS reanalysis for the same valid date\.d, Correctly predicted exceedances \(hits\), missed exceedances, and false alarms at each cell’s local 2\-year return level; cells with no forecast or reference exceedance are shown in grey\. The valid date was selected independently of Hapi as the 2024 test date with the largest number of reference 2\-year\-threshold exceedances over CONUS\. The discharge maps share a logarithmic colour scale\.Deep learning provides a complementary approach, but most hydrological applications have represented rivers as collections of individual basins or gauges\. Within this setting, long short\-term memory networks \(LSTMs\), are widely used to predict streamflow\([Nearing et al\., 2021](https://arxiv.org/html/2609.22702#bib.bib33)\)\. A single model trained across many catchments can outperform individually calibrated process models when conditioned on basin attributes\([Kratzert et al\., 2019](https://arxiv.org/html/2609.22702#bib.bib11);[Newman et al\., 2015](https://arxiv.org/html/2609.22702#bib.bib28);[Addor et al\., 2017](https://arxiv.org/html/2609.22702#bib.bib29)\)\. This approach now spans multiple timescales, global regionalization, and operational flood warning\([Gauch et al\., 2021](https://arxiv.org/html/2609.22702#bib.bib26);[Beck et al\., 2020](https://arxiv.org/html/2609.22702#bib.bib34);[Nearing et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib12);[Nevo et al\., 2022](https://arxiv.org/html/2609.22702#bib.bib13)\)\. Caravan\([Kratzert et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib30)\)and NeuralHydrology\([Kratzert et al\., 2022](https://arxiv.org/html/2609.22702#bib.bib27)\)support standardized development and evaluation of these basin\-based models\. Transformers and time\-series foundation models have also been applied to 120\-h streamflow prediction and zero\-shot forecasting at unseen stations\([Demiray and Demir, 2026](https://arxiv.org/html/2609.22702#bib.bib44);[Sun and Sun, 2026](https://arxiv.org/html/2609.22702#bib.bib45)\)\. Despite their broader reach, these models produce basin or gauge sequences without explicit spatial interaction across a gridded domain\.

Gridded models extend deep\-learning discharge prediction beyond this pointwise representation\. RiverMamba\([Shams Eddin et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib15)\)brings spatial context to global0\.05∘0\.05^\{\\circ\}discharge forecasting and outperforms an operational LSTM baseline, particularly beyond 48 h\. SwinFlood\([Song et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib14)\)combines a convolutional neural network with Swin Transformer blocks within a single urban–rural catchment on the Shenzhen River\. These models extend hydrological deep learning from independent basin sequences to spatially distributed discharge forecasts\.

Related developments in weather and Earth\-system forecasting offer strategies for handling the spatial and multivariable dimensions of this problem\. Hierarchical and graph\-based models build broad spatial context without applying global attention to every grid point\([Bi et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib7);[Lam et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib8);[Chen et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib9);[Nguyen et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib35)\)\. Generalist Earth\-system models forecast heterogeneous fields through shared architectures\([Bodnar et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib23);[Nguyen et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib10);[Schmude et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib37);[Vaughan et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib38)\)\. Similarly, ecLand emulators show that several prognostic land\-surface states can be advanced together\([Wesselkamp et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib43)\)\. A complementary line of hydrological research encodes conservation or process structure through mass\-conserving recurrent networks, differentiable process models, or differentiable routing\([Hoedt et al\., 2021](https://arxiv.org/html/2609.22702#bib.bib25);[Feng et al\., 2022](https://arxiv.org/html/2609.22702#bib.bib31);[Bindas et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib32)\)\. Together, these studies provide ingredients for spatial or multivariable prediction, but they do not jointly establish high\-resolution continental forecasting of multiple hydrological states, threshold\-based flood skill, independent gauge validation, and computationally efficient inference in one model\.

To address this gap, we introduceHapi, a hierarchical transformer that forecasts river discharge, surface runoff, snow water equivalent, and soil wetness over the contiguous United States \(CONUS\) out to 72 h while representing local hydrological variation and basin\-to\-synoptic\-scale context\. Figure[1](https://arxiv.org/html/2609.22702#S1.F1)shows an example of Hapi’s 72\-h discharge forecast across CONUS\. Hapi combines fine three\-dimensional patches, which retain local detail, with hierarchical shifted\-window attention that expands spatial context without global attention over the full grid\. Learned task weights account for the different scales and distributions of the four targets\. During autoregressive rollout, the four predicted fields feed into the next step, while meteorological forcing is updated externally\. On a 2024 CONUS holdout, we evaluate reanalysis\-conditioned flood detection and continuous discharge skill against persistence, the operational GloFAS forecast, and RiverMamba\. ERA5\-Land fields supply the meteorological forcing in this experiment, isolating hydrological\-state prediction from errors in numerical weather prediction \(NWP\); a precipitation\-perturbation experiment probes sensitivity to departures from this forcing\. Our contributions are

- •A hierarchical vision architecture for continental hydrology\.Hapi preserves the two\-dimensional grid, uses fine three\-dimensional convolutional patches to limit early spatial aggregation, and expands context through hierarchical Swin V2 stages, as described in Section[3](https://arxiv.org/html/2609.22702#S3)\. A mirrored decoder and multiscale skip connections then reconstruct fine spatial structure while the hierarchy builds broader meteorological context within one continental model\.
- •Automatic balancing of heterogeneous hydrological targets\.Hapi forecasts discharge, runoff, snow water equivalent, and soil wetness without manually tuning their relative loss weights, despite their different scales and error distributions, as described in Section[3\.3](https://arxiv.org/html/2609.22702#S3.SS3)\. Under matched training conditions, learned task weighting improves discharge accuracy and flood\-detection skill relative to uniform or fixed weighting, with the largest gains at rarer thresholds, as shown in Section[4\.5](https://arxiv.org/html/2609.22702#S4.SS5)\.
- •Leading flood detection with independent observational validation\.Hapi achieves the highest flood\-detection F1 in 20 of 21 lead–threshold comparisons in Section[4\.2](https://arxiv.org/html/2609.22702#S4.SS2), including every threshold from 1\.5 to 100 yr at 48 and 72 h\. All 21 F1 advantages over RiverMamba remain significant after multiplicity correction, as shown in Appendix[N](https://arxiv.org/html/2609.22702#A14)\. Hapi also has the highest base\-rate\-invariant symmetric extremal dependence index, SEDI, in 18 of 21 comparisons with GloFAS and persistence in Appendix[O](https://arxiv.org/html/2609.22702#A15)\. Against 3,881 U\.S\. Geological Survey gauges, hereafter USGS gauges, Hapi also achieves higher median Nash–Sutcliffe efficiency, NSE, than GloFAS and persistence at every lead, as shown in Section[4\.3](https://arxiv.org/html/2609.22702#S4.SS3)\.

## 2Data

We use a daily gridded dataset assembled from GloFAS reanalysis\([Harrigan et al\., 2020](https://arxiv.org/html/2609.22702#bib.bib3)\), ERA5\-Land\([Copernicus Climate Change Service, 2019](https://arxiv.org/html/2609.22702#bib.bib51)\), and static terrain and upstream drainage fields\. The dataset spans twenty years, 2005–2024, over CONUS on a0\.05∘0\.05^\{\\circ\}latitude–longitude grid, approximately 5 km, at a 24\-h time step\. ERA5\-Land data are regridded from their native0\.1∘0\.1^\{\\circ\}resolution to the GloFAS grid\. The GloFAS grid consists of512×1152=589,824512\\times 1152=589\{,\}824cells\. Ocean cells and cells outside CONUS are zero\-padded to preserve the rectangular grid shape and masked out when computing the loss and evaluation metrics, leaving 416,298 active land cells, as described in Appendix[B](https://arxiv.org/html/2609.22702#A2)\.

### 2\.1Variables

Each sample from the dataset contains 12 input variables grouped by their role in the autoregressive rollout described in Section[3](https://arxiv.org/html/2609.22702#S3)\. Table[1](https://arxiv.org/html/2609.22702#S2.T1)lists their sources, units, and definitions\.

Table 1:The 12 model inputs, grouped by their role in the autoregressive rollout\. ERA5\-Land fields come from the Copernicus Climate Change Service Climate Data Store\([Copernicus Climate Change Service, 2019](https://arxiv.org/html/2609.22702#bib.bib51)\); upstream drainage area comes from the GloFAS v4\.0 auxiliary data\([ECMWF, 2023](https://arxiv.org/html/2609.22702#bib.bib4)\)\. DEM denotes digital elevation model\.ERA5\-Land supplies the six forcing variables at each rollout step, making the prediction of hydrological states reanalysis\-conditioned\. The two runoff variables differ because GloFAS runoff is routed and prognostic, whereas ERA5\-Landrois an unrouted forcing\. Similarly,swvl1describes the top 7 cm, while the prognostic soil wetness index aggregates the top three layers with vegetation\-dependent weights\. Appendix[A](https://arxiv.org/html/2609.22702#A1)reports the variable\-specific normalization choices and the distributions that motivate them\.

### 2\.2Data Splits

We follow the WeatherBench convention and split the data by calendar year\. Models are trained on 2005–2022, comprising eighteen years and approximately 6,575 daily samples, while the 365 samples from 2023 are used only for validation\. The held\-out test set covers the first ten months of 2024, from 2024\-01\-01 through 2024\-11\-01, after accounting for theT=2T\{=\}2history buffer at the start of the window\. It is used only for the results in Section[4](https://arxiv.org/html/2609.22702#S4)\. The three calendar ranges are disjoint and ordered in time\.

### 2\.3Flood Detection Thresholds

For flood detection, we use the per\-cell return levels distributed as the GloFAS v4\.0 auxiliary data to define flood events\([ECMWF, 2023](https://arxiv.org/html/2609.22702#bib.bib4)\)\. GloFAS derives these fields by fitting a Gumbel distribution, theξ=0\\xi=0member of the generalized extreme value family, to annual maxima using L\-moments\. The annual maxima come from LISFLOOD\([Burek et al\., 2013](https://arxiv.org/html/2609.22702#bib.bib5)\)physical simulations forced with ERA5 in the GloFAS reanalysis\([Harrigan et al\., 2020](https://arxiv.org/html/2609.22702#bib.bib3)\)for 1979–2022\. The fitted distributions provide return levels for 1\.5, 2, 5, 10, 20, 50, 100, 200, and 500 years\. At each grid cell𝐱\\mathbf\{x\}, therr\-year return levelqr​\(𝐱\)q\_\{r\}\(\\mathbf\{x\}\)is the discharge threshold satisfyingPr\[Qmax\(𝐱\)\>qr\(𝐱\)\]=1/r\\Pr\[Q\_\{\\max\}\(\\mathbf\{x\}\)\>q\_\{r\}\(\\mathbf\{x\}\)\]=1/r, whereQmax​\(𝐱\)Q\_\{\\max\}\(\\mathbf\{x\}\)is the annual maximum discharge\. Thus, a 100\-year return level has a 1% chance of being exceeded in any given year, and exceedances can occur in consecutive years\. Using these published fields keeps the event definition consistent with the one used by GloFAS for flood warnings\. The fitting window ends in 2022, before the 2024 test window, so the thresholds carry no information from the validation or test periods\.

The return\-level fields align exactly with our grid, allowing us to select a CONUS subset without interpolation, as described in Appendix[B](https://arxiv.org/html/2609.22702#A2)\. Cells that are undefined in the source fields are excluded, leaving414,795414\{,\}795valid cells within the416,298416\{,\}298\-cell CONUS land mask\.

For each return periodrr, a cell is labeled as flooded when the discharge exceeds therr\-year threshold\. The corresponding predicted flood mask is obtained by applying the same return level to the model’s predicted discharge\.

## 3Method

Hapi is a U\-Net Swin Transformer V2 for gridded multivariable hydrology at0\.05∘0\.05^\{\\circ\}resolution\. Its backbone integrates a U\-Net\([Ronneberger et al\., 2015](https://arxiv.org/html/2609.22702#bib.bib17)\)encoder–decoder, Swin Transformer V2 blocks\([Liu et al\., 2022](https://arxiv.org/html/2609.22702#bib.bib16)\), 2D rotary position embeddings\([Su et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib18);[Fang et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib19)\), and gated attention\([Qiu et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib39)\)\. We also use learned per\-variable Laplacian task weighting in place of the hand\-tuned weights used by prior multivariable forecasters, together with multi\-step fine\-tuning at a randomized rollout depth for the four prognostic land\-surface variables\. GivenTThistorical frames ofCinC\_\{\\mathrm\{in\}\}input variables on anH×WH\\times Wgrid, the model predictsCoutC\_\{\\mathrm\{out\}\}prognostic variables at the next time step, as shown in Figure[2](https://arxiv.org/html/2609.22702#S3.F2)\.

Existing gridded Earth\-system forecasters take different architectural routes, which typically suit the atmosphere better than the land surface\. RiverMamba\([Shams Eddin et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib15)\)traverses the grid as a 1D sequence with bidirectional Mamba blocks\([Gu and Dao, 2024](https://arxiv.org/html/2609.22702#bib.bib22)\), committing to a serial scan order on a 2D field\. Pangu\-Weather\([Bi et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib7)\)adds a 3D Earth\-specific transformer with an explicit pressure\-level axis, matched to atmospheric vertical structure but unnatural for surface\-only variables\. GraphCast\([Lam et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib8)\)uses a multi\-mesh icosahedral graph that trades resolution uniformity for global coverage\. FuXi\([Chen et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib9)\)and Aurora\([Bodnar et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib23)\)use U\-shaped vision transformers similar to our backbone but designed for atmospheric targets\. We use the Swin hierarchy because its multiscale windows span river\-channel to synoptic scales\. Hapi retains a two\-dimensional hierarchy and reconstructs full\-resolution outputs through a mirrored decoder and skip connections\.

Hapi also differs from existing models in its token representation\. Strategies include per\-variable tokenization with cross\-attention\([Nguyen et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib10)\), physically grouped branches\([Bi et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib7)\), and variable\-agnostic latent encoders\([Jaegle et al\., 2022](https://arxiv.org/html/2609.22702#bib.bib41);[Bodnar et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib23)\)\. Hapi keeps a monolithic patch embedding but distinguishes prognostic, forcing, and static inputs at the rollout level, a separation described in Section[3\.2](https://arxiv.org/html/2609.22702#S3.SS2)\. The hydrological targets span six orders of magnitude, and up to 78% of their values are exactly zero, as shown in Table[12](https://arxiv.org/html/2609.22702#A1.T12), making uniform standardization unsuitable\. We combine per\-variable log1p reparameterization with learnable Laplacian task weighting in Section[3\.3](https://arxiv.org/html/2609.22702#S3.SS3), adapting the likelihood\-based weighting form of[Kendall et al\. \(2018\)](https://arxiv.org/html/2609.22702#bib.bib20)to a mean absolute error \(MAE\), objective for heavy\-tailed hydrological targets\. For autoregressive rollout, AI weather models stabilize long horizons with curriculum fine\-tuning\([Lam et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib8)\), randomized forecast intervals\([Nguyen et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib35)\), or diffusion\-based forecasting\([Price et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib24);[Hatanpää et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib36)\); we fine\-tune on rollouts of randomized depth, in the manner of[Brandstetter et al\. \(2022\)](https://arxiv.org/html/2609.22702#bib.bib21)\.

### 3\.1Architecture

To limit early aggregation of local hydrological variation, Hapi uses smallp×pp\\times pspatial patches withp=2p=2\. A three\-dimensional convolutional \(Conv3D\) layer, with kernel and stride\(T,p,p\)\(T,p,p\)maps the input tensor𝐗∈ℝB×Cin×T×H×W\\mathbf\{X\}\\in\\mathbb\{R\}^\{B\\times C\_\{\\mathrm\{in\}\}\\times T\\times H\\times W\}toHp×Wp\\frac\{H\}\{p\}\\times\\frac\{W\}\{p\}tokens of dimensiondd, whereBBis the batch size\. This projection combines the variables andTThistorical frames within each patch, leaving a two\-dimensional token grid for subsequent attention\.

Hapi builds spatial context throughNNstages of Swin Transformer blocks to reduce attention costs\. Attention operates within local windows, and shifted windows allow information to pass across window boundaries\. Between stages, patch merging halves each spatial dimension and doubles the channel width\. This encoder expands the area covered by each window while reducing the number of tokens at the next stage\.

The decoder mirrors the encoder’s structure by upsampling through patch unmerging, which uses a learned linear expansion and pixel\-shuffle rearrangement to double each spatial dimension\. Skip connections concatenate encoder and decoder features at corresponding resolutions, giving the decoder direct access to encoder representations that bypass the bottleneck\. A final linear projection producesCoutC\_\{\\mathrm\{out\}\}channels at the originalH×WH\\times Wresolution, one predicted field per prognostic variable\.

![Refer to caption](https://arxiv.org/html/2609.22702v1/arch.png)Figure 2:Hapi architecture\. The encoder on the left downsamples through patch merging acrossNNstages; each Swin block uses 2D rotary position embeddings and gated attention, both described below\. The decoder on the right upsamples through patch unmerging, with skip connections transferring fine\-grained features from matching encoder stages\. During autoregressive rollout, prognostic channels are predicted and fed back, forcing channels are supplied externally, and static channels remain fixed, as described in Section[3\.2](https://arxiv.org/html/2609.22702#S3.SS2)\.#### Position encoding\.

Within an attention window, token interactions depend on relative spatial position as well as the local features\. Swin V2 represents position through a learned multilayer perceptron \(MLP\), that maps coordinate offsets\(Δ​h,Δ​w\)\(\\Delta h,\\Delta w\)to an additive attention bias\. We use axial 2D rotary position embeddings \(RoPE\)\([Su et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib18);[Heo et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib40)\), to encode these offsets directly in the query–key interaction without learning a position\-bias network\. For each attention head of dimensiondhd\_\{h\}, we group vector components into even–odd pairs and alternate row and column rotations across successive pairs\. For a query vector𝐪\\mathbf\{q\}at window position\(h,w\)\(h,w\), the rotation is

RoPE⁡\(𝐪2​k,𝐪2​k\+1,h,w\)=\(cos⁡ϕk−sin⁡ϕksin⁡ϕkcos⁡ϕk\)​\(𝐪2​k𝐪2​k\+1\),\\mathrm\{RoPE\}\(\\mathbf\{q\}\_\{2k\},\\mathbf\{q\}\_\{2k\+1\};h,w\)=\\begin\{pmatrix\}\\cos\\phi\_\{k\}&\-\\sin\\phi\_\{k\}\\\\ \\sin\\phi\_\{k\}&\\cos\\phi\_\{k\}\\end\{pmatrix\}\\begin\{pmatrix\}\\mathbf\{q\}\_\{2k\}\\\\ \\mathbf\{q\}\_\{2k\+1\}\\end\{pmatrix\},\(1\)where the angles are defined by

ϕ2​j\(h,w\)=hωj,ϕ2​j\+1\(h,w\)=wωj,ωj=θ−4j/dh,\\phi\_\{2j\}\(h,w\)=h\\,\\omega\_\{j\},\\qquad\\phi\_\{2j\+1\}\(h,w\)=w\\,\\omega\_\{j\},\\qquad\\omega\_\{j\}=\\theta^\{\-4j/d\_\{h\}\},\(2\)withj=0,…,dh/4−1j=0,\\ldots,d\_\{h\}/4\-1indexing the frequency bands\. The coordinateshhandwwindex the row and column within each attention window\. Applying the same rotations to key vectors makes the query–key dot product depend on position through relative row and column offsets\. We setθ=10\\theta=10for the8×88\\times 8windows\.

#### Gated attention\.

A Swin window can contain river channels, dry land, and ocean cells, so the relevance of neighboring information can vary across locations and features\. Within each head, standard softmax attention uses one normalized set of weights to mix all features of the value vectors\. Consequently, it passes a mixture of values forward, even when none provides useful information for a particular feature of the receiving token\. We add separate, query\-dependent control using the element\-wise sigmoid gate of[Qiu et al\. \(2025\)](https://arxiv.org/html/2609.22702#bib.bib39), computed from the receiving token and applied after attention\. The gate allows each token to attenuate individual attended features before the output projection\. For input token features𝐗\\mathbf\{X\}within a window, we compute the gate logits𝐆\\mathbf\{G\}together with the query, key, and value projections:

\[𝐐,𝐊,𝐕,𝐆\]=𝐗𝐖q​k​v​g\+𝐛q​k​v​g\.\[\\mathbf\{Q\},\\mathbf\{K\},\\mathbf\{V\},\\mathbf\{G\}\]=\\mathbf\{X\}\\mathbf\{W\}\_\{qkvg\}\+\\mathbf\{b\}\_\{qkvg\}\.\(3\)The gated attention output is

𝐎=σ⁡\(𝐆\)⊙Attn⁡\(𝐐,𝐊,𝐕\),\\mathbf\{O\}=\\sigma\(\\mathbf\{G\}\)\\odot\\mathrm\{Attn\}\(\\mathbf\{Q\},\\mathbf\{K\},\\mathbf\{V\}\),\(4\)whereσ\\sigmais the sigmoid function and⊙\\odotdenotes element\-wise multiplication\. The gateσ⁡\(𝐆\)\\sigma\(\\mathbf\{G\}\)assigns a scale between zero and one to each feature in every token and attention head\.

### 3\.2Autoregressive Rollout with Prognostic–Forcing Separation

The model rolls out by feeding its output back as input\. Variables update by different rules at each step:prognosticvariables, comprising river discharge, runoff, snow water equivalent, and soil wetness index, are replaced by their predicted fields;forcingvariables, including precipitation, temperature, and radiation, are supplied from ERA5\-Land;staticdescriptors, comprising elevation and upstream drainage area, stay fixed\. This separates the prognostic hydrological state from meteorological forcing: forecast products can supply the forcing fields, while terrain remains fixed at every step\. Prognostic predictions are held within declared physical ranges before they are fed back: the nonnegative variables are snapped to the point mass at zero in their target distributions, and the soil wetness index is clipped to\[0,1\]\[0,1\], as described in Appendix[C](https://arxiv.org/html/2609.22702#A3)\. The rule is applied only at intermediate steps and is identical during training and inference; the loss at the final step is computed from the raw output\. Section[3\.3](https://arxiv.org/html/2609.22702#S3.SS3)describes the training schedule for multi\-step fine\-tuning at a randomized rollout depth\.

### 3\.3Training Objective

Foundation models such as Pangu\-Weather\([Bi et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib7)\), GraphCast\([Lam et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib8)\), FuXi\([Chen et al\., 2023](https://arxiv.org/html/2609.22702#bib.bib9)\), and GenCast\([Price et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib24)\)balance per\-variable losses with hand\-tuned scalar weights\. RiverMamba also fixes its severity and lead\-time weights by hand\([Shams Eddin et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib15)\), although its single target does not require cross\-variable balancing\. We adapt the likelihood\-based Gaussian task weighting of[Kendall et al\. \(2018\)](https://arxiv.org/html/2609.22702#bib.bib20), which uses mean squared error, to a Laplacian form with MAE for hydrology’s heavy\-tailed targets with frequent zeros\. Per\-variable log\-scales are learned end to end\. Loss is computed only on theCpC\_\{p\}prognostic outputs, soCout=CpC\_\{\\mathrm\{out\}\}=C\_\{p\}\.

Heavy\-tailed targets, comprising discharge, snow water equivalent, and runoff, are reparameterized withlog⁡\(1\+x\)\\log\(1\+x\)so that rare large values do not dominate the loss, as described in Appendix[A](https://arxiv.org/html/2609.22702#A1)\. For a prediction𝐘^∈ℝB×Cp×H×W\\hat\{\\mathbf\{Y\}\}\\in\\mathbb\{R\}^\{B\\times C\_\{p\}\\times H\\times W\}and target𝐘\\mathbf\{Y\}, define the residualδb,i,h,w=Y^b,i,h,w−Yb,i,h,w\\delta\_\{b,i,h,w\}=\\hat\{Y\}\_\{b,i,h,w\}\-Y\_\{b,i,h,w\}in normalized space\. The masked mean absolute error for variableiiis

MAEi=∑b,h,wmi,h,w​\|δb,i,h,w\|B​∑h,wmi,h,w,i=1,…,Cp,\\mathrm\{MAE\}\_\{i\}=\\frac\{\\sum\_\{b,h,w\}m\_\{i,h,w\}\\,\\lvert\\delta\_\{b,i,h,w\}\\rvert\}\{B\\sum\_\{h,w\}m\_\{i,h,w\}\},\\qquad i=1,\\dots,C\_\{p\},\(5\)wheremi∈\{0,1\}H×Wm\_\{i\}\\in\\\{0,1\\\}^\{H\\times W\}is the per\-variable validity mask of Appendix[B](https://arxiv.org/html/2609.22702#A2), so zero\-padded ocean and out\-of\-domain cells neither contribute to the numerator nor inflate the denominator\.

Each variable has a learnable log\-scalesis\_\{i\}\. Writingδi\\delta\_\{i\}for a scalar residual of variableii, the Gaussian formulation of[Kendall et al\. \(2018\)](https://arxiv.org/html/2609.22702#bib.bib20)gives the negative log\-likelihood12​e−2​si​δi2\+si\\tfrac\{1\}\{2\}e^\{\-2s\_\{i\}\}\\delta\_\{i\}^\{2\}\+s\_\{i\}up to an additive constant, withsis\_\{i\}denoting the log\-standard deviation\. For our MAE objective, we use a zero\-centered Laplacian residual model with scaleesie^\{s\_\{i\}\}:

p⁡\(δi∣si\)=12​e−si​exp⁡\[−e−si​\|δi\|\]\.p\(\\delta\_\{i\}\\mid s\_\{i\}\)=\\frac\{1\}\{2\}e^\{\-s\_\{i\}\}\\,\\exp\\\!\\left\[\-e^\{\-s\_\{i\}\}\\,\\lvert\\delta\_\{i\}\\rvert\\right\]\.\(6\)Averaging its negative log\-likelihood over valid cells and the batch, and dropping the constantlog⁡2\\log 2, givese−si​MAEi\+sie^\{\-s\_\{i\}\}\\,\\mathrm\{MAE\}\_\{i\}\+s\_\{i\}\. Averaging over variables yields

ℒ=1Cp​∑i=1Cp\(e−si​MAEi\+si\),\\mathcal\{L\}=\\frac\{1\}\{C\_\{p\}\}\\sum\_\{i=1\}^\{C\_\{p\}\}\\bigl\(e^\{\-s\_\{i\}\}\\,\\mathrm\{MAE\}\_\{i\}\+s\_\{i\}\\bigr\),\(7\)wheree−sie^\{\-s\_\{i\}\}is the inverse Laplace scale and serves as the learned task weight\. Increasingsis\_\{i\}reduces this weight, and the additivesis\_\{i\}term prevents all weights from collapsing to zero\. For fixedMAEi\>0\\mathrm\{MAE\}\_\{i\}\>0, stationarity with respect tosis\_\{i\}givese−si=1/MAEie^\{\-s\_\{i\}\}=1/\\mathrm\{MAE\}\_\{i\}\. The learned weights therefore balance the variables according to their error scales, removing the need to tune their relative loss weights by hand\. Thesis\_\{i\}are global task\-level scales shared across samples, locations, and lead times\. They adapt the relative contributions of the four prognostic variables during optimization; they are not input\-dependent predictive uncertainties and do not provide calibrated forecast intervals\. The implementation also admits fixed per\-variable prior weights multiplying each term; they are uniform in every run reported here, and Appendix[F](https://arxiv.org/html/2609.22702#A6)gives the converged\{si\}\\\{s\_\{i\}\\\}\.

Multi\-step fine\-tuning follows[Brandstetter et al\. \(2022\)](https://arxiv.org/html/2609.22702#bib.bib21)by randomizing the rollout horizon\. At each iteration, we sampleJ∼Uniform⁡\{1,…,Jmax\}J\\sim\\mathrm\{Uniform\}\\\{1,\\dots,J\_\{\\max\}\\\}, autoregressively roll out the model forJJsteps, and evaluate the loss on the prediction at stepJJ:

ℒtotal=𝔼J∼Uniform⁡\{1,…,Jmax\}​\[ℒ\(J\)\]\.\\mathcal\{L\}\_\{\\mathrm\{total\}\}=\\mathbb\{E\}\_\{J\\sim\\mathrm\{Uniform\}\\\{1,\\dots,J\_\{\\max\}\\\}\}\\left\[\\mathcal\{L\}^\{\(J\)\}\\right\]\.\(8\)The intermediate predictions are fed back into the model without detaching the computational graph, so gradients fromℒ\(J\)\\mathcal\{L\}^\{\(J\)\}propagate through the completeJJ\-step rollout\. We useJmax=3J\_\{\\max\}\{=\}3\.

### 3\.4Model Configuration

The final model uses two consecutive daily frames,T=2T\{=\}2, of 12 input variables to predict the next\-day fields of four prognostic variables\. Longer histories produced only marginal accuracy gains at greater computational cost\. Although the temporal projection can be reconfigured for any history length,T=2T\{=\}2provided the best accuracy–efficiency trade\-off among those tested\.

The encoder comprises three stages with\(1,2,4\)\(1,2,4\)Swin blocks and uses embedding dimensiond=192d\{=\}192, patch sizep=2p\{=\}2, and window sizew=8w\{=\}8\. The decoder mirrors this stage structure\. The complete model has 43\.96 million parameters\. In a Bayesian search over these four hyperparameters, shallow encoders with fine tokenization achieved the lowest validation losses\. The selected configuration follows this pattern, as shown in Appendix[E](https://arxiv.org/html/2609.22702#A5)\.

Training begins with single\-step pretraining and continues with randomized multi\-step fine\-tuning usingJmax=3J\_\{\\max\}\{=\}3\. We optimize the task log\-scales\{si\}\\\{s\_\{i\}\\\}during pretraining and freeze them during fine\-tuning\. Freezing the scales prevents them from compensating for accumulated rollout error, leaving the model parameters to reduce that error, as described in Appendix[D](https://arxiv.org/html/2609.22702#A4)\. Appendix[F](https://arxiv.org/html/2609.22702#A6)gives the full training details and learned log\-scales\.

## 4Results

We evaluated Hapi from overall forecast skill to specific capabilities and limitations\. We first compared discharge forecasts with persistence, GloFAS, and RiverMamba at continental and gauge scales\. We then examined the contributions of multivariable prediction and learned task weighting, evaluated the other prognostic variables, and tested variation across river subsets and perturbed meteorological forcing\. Finally, we used Hurricane Helene to illustrate event\-scale behavior and assessed computational performance\.

### 4\.1Baselines, Metrics, and Statistical Analysis

We evaluated gridded discharge against GloFAS reanalysis using a common verification target and compared Hapi with persistence, the operational GloFAS forecast, and RiverMamba\. The comparisons cover 24\-, 48\-, and 72\-h leads during the 2024 test period\. Flood\-threshold detection is the primary outcome, complemented by continuous discharge metrics, per\-cell analyses, and paired moving\-block bootstrap inference\. The verification target is the GloFAS reanalysis 24\-h mean river discharge, the field Hapi is trained to predict in Section[2\.1](https://arxiv.org/html/2609.22702#S2.SS1)and from which the per\-cell return\-level thresholds are derived in Section[2\.3](https://arxiv.org/html/2609.22702#S2.SS3)\. Scores therefore measure agreement with that reanalysis, using an identical target for every model\.

We compare Hapi against three baselines\.Persistencecarries forward the last discharge field observed at forecast issue time and holds it constant across all leads; the 72\-h forecast therefore uses the field from three days before the target\. The lag\-1 autocorrelation of daily streamflow exceeds 0\.95 on large catchments, making persistence a strong short\-lead baseline\.GloFAS forecastis the operational, physics\-based global flood forecasting system of the Copernicus Emergency Management Service\([Alfieri et al\., 2013](https://arxiv.org/html/2609.22702#bib.bib2)\)\. The GloFAS v4 forecasts evaluated here are generated with the distributed hydrological model LISFLOOD OS\([Burek et al\., 2013](https://arxiv.org/html/2609.22702#bib.bib5)\)\. We use the Copernicus Emergency Management Service productcems\-glofas\-forecast, with system versionoperational, hydrological model LISFLOOD, thecontrol\_forecastmember, and variableriver\_discharge\_in\_the\_last\_24\_hours, at lead times of 24, 48, and 72 h, retrieved over the same CONUS box as the training grid\. This operational forecast is evaluated against the same verification target as the other models\.

RiverMamba\([Shams Eddin et al\., 2025](https://arxiv.org/html/2609.22702#bib.bib15)\)is a recent gridded deep\-learning discharge forecaster at the same0\.05∘0\.05^\{\\circ\}resolution\. We follow its published issue\-time convention: RiverMamba uses seven historical days throught−1t\-1, issues attt, and predictst\+1t\+1throught\+7t\+7, without using a day\-ttdischarge nowcast\. Thus the tabulated 24\-, 48\-, and 72\-h leads follow each system’s stated forecast issue time, although RiverMamba’s latest discharge input is one day older than Hapi’s at the same nominal lead\. Only Hapi and GloFAS carry the four prognostic land\-surface variables; RiverMamba is discharge\-only\. RiverMamba was run on its native global domain for forecast evaluation, and its CONUS predictions were extracted for the comparisons below\.

Each system is evaluated with its respective meteorological inputs\. Hapi uses ERA5\-Land forcing over the forecast horizon, whereas the operational GloFAS product uses its native forecast forcing\. The Hapi–GloFAS comparison therefore measures reanalysis\-conditioned Hapi skill against the archived operational product, not performance under identical meteorological inputs\.

We exclude the U\.S\. National Water Model because its streamflow predictions are defined on NHDPlus river reaches rather than the regular grid used here\([Cosgrove et al\., 2024](https://arxiv.org/html/2609.22702#bib.bib6)\)\. A direct comparison would require a separately validated reach\-to\-grid mapping or an independent gauge evaluation that includes the National Water Model, both outside the scope of this study\.

Continuous performance is evaluated using MAE and root mean squared error \(RMSE\), in m3s\-1, together with the dimensionless NSE and Kling–Gupta efficiency, KGE\. The primary flood\-detection analysis uses per\-cell GloFAS v4\.0 return levels at seven return periods from 1\.5 to 100 yr and reports F1, precision, recall, and the critical success index \(CSI\)\. Exploratory results for the more uncertain 200\- and 500\-yr levels are reported in Appendix[L](https://arxiv.org/html/2609.22702#A12)\. Definitions and aggregation conventions for all evaluation metrics are given in Appendix[H](https://arxiv.org/html/2609.22702#A8)\. The primary discharge metrics use the 2024 test window at 24\-, 48\-, and 72\-h leads on the CONUS land mask, which contains 414,795 valid cells after excluding 1,503 cells without defined return levels\.

Threshold\-crossing skill is the primary outcome, and continuous skill is secondary\. The decision a flood forecast supports is binary: at a given return period, warn or do not warn\. Threshold\-based metrics directly quantify this decision, whereas continuous error statistics do not\. A forecast can have low RMSE while misclassifying threshold exceedances during flood events\. Unadjusted accuracy is uninformative at these base rates because a forecast with no warnings would score close to one\. F1 balances missed events against false alarms without requiring a cost ratio, which depends on the application\.

F1 depends on the event base rate, precluding direct comparison between values at different return\-period thresholds, such as 1\.5 and 100 yr\. We use F1 to compare models at a fixed return period and the base\-rate\-invariant SEDI to compare skill across return periods, as described in Appendix[O](https://arxiv.org/html/2609.22702#A15)\. Pooled continuous metrics characterize domain\-wide performance and emphasize high\-discharge rivers\. We therefore also compute the metrics independently at each grid cell and summarize their spatial distributions in Appendix[M](https://arxiv.org/html/2609.22702#A13)\.

Uncertainty in F1 and SEDI differences is estimated with a paired moving\-block bootstrap using 30\-day blocks andR=10,000R=10\{,\}000replicates, as described in Appendix[N](https://arxiv.org/html/2609.22702#A14)\. Within each replicate, the same days are resampled for both models\. Holm correction is applied across the 21 lead–return\-period comparisons within each model pair and metric\. Hapi–persistence and Hapi–GloFAS comparisons use 302 common dates per lead, whereas Hapi–RiverMamba comparisons use the 294 dates available for both models\. We report paired differences with percentile intervals and multiplicity\-adjusted decisions\.

### 4\.2Discharge Forecast Skill

Hapi achieves the highest flood\-detection F1 in 20 of the 21 lead–threshold comparisons in Table[2](https://arxiv.org/html/2609.22702#S4.T2), including all seven thresholds from 1\.5 to 100 yr at both 48 and 72 h\. Hapi exceeds RiverMamba throughout and falls below persistence only at the 100\-yr threshold at 24 h, by 0\.0015\. Figure[3](https://arxiv.org/html/2609.22702#S4.F3)a shows that Hapi’s margin over the strongest baseline is larger at 48 and 72 h than at 24 h across all seven thresholds\. At the 1\.5\-yr and 100\-yr thresholds shown in Figure[3](https://arxiv.org/html/2609.22702#S4.F3)b,c, Hapi has the lowest false\-alarm ratio among the four models at every lead\.

We next used the paired moving\-block bootstrap described in Appendix[N](https://arxiv.org/html/2609.22702#A14)to quantify uncertainty in the F1 comparisons reported in Table[2](https://arxiv.org/html/2609.22702#S4.T2)\. Hapi’s higher F1 than RiverMamba was significant after Holm–Bonferroni correction in all 21 lead–threshold comparisons\. Against GloFAS, Hapi’s higher F1 remained significant at every threshold at 24 h\. At the 5\-yr threshold, the Hapi–GloFAS difference was\+0\.203\+0\.203at 24 h, with a 95% interval of0\.1650\.165–0\.2270\.227\. The difference was\+0\.145\+0\.145at 72 h, with a 95% interval of0\.0870\.087–0\.1820\.182; both comparisons were significant after Holm–Bonferroni correction\.

Among Hapi, GloFAS, and persistence, Hapi has the highest SEDI in 18 of 21 lead–threshold comparisons, including all seven at both 48 and 72 h\. The SEDI analysis in Appendix[O](https://arxiv.org/html/2609.22702#A15)complements the fixed\-threshold F1 comparisons by allowing skill to be compared across event rarity\. At 72 h, Hapi’s SEDI decreases from 0\.931 at the 1\.5\-yr threshold to 0\.724 at the 100\-yr threshold, indicating lower flood\-detection skill at more extreme thresholds\. Unlike F1, this decrease is not an automatic consequence of the lower event frequency\.

In Table[3](https://arxiv.org/html/2609.22702#S4.T3), Hapi has the lowest pooled MAE and RMSE at 24 h, but GloFAS leads on RMSE, NSE, and KGE at 48 and 72 h\. At 48 h, Hapi and GloFAS have nearly identical MAE, both rounding to 2\.94 m3/s; GloFAS has the lowest MAE at 72 h, while persistence has the highest KGE at 24 h\. Pooling emphasizes high\-discharge rivers in a domain spanning five orders of magnitude of discharge, motivating a separate evaluation at individual grid cells\.

Although GloFAS has higher pooled NSE at 48 and 72 h, Hapi has higher median per\-cell NSE at both leads in Table[4](https://arxiv.org/html/2609.22702#S4.T4)\. Among Hapi, GloFAS, and persistence, Hapi has the highest median per\-cell NSE at every lead: 0\.902, 0\.823, and 0\.776, compared with 0\.698, 0\.637, and 0\.559 for GloFAS\. At 72 h, Hapi exceeds GloFAS at 79\.2% of scored cells and persistence at 92\.3%, with median NSE differences of\+0\.142\+0\.142and\+0\.363\+0\.363, respectively\. These spatial comparisons show that Hapi’s NSE advantage extends across most scored cells despite its lower pooled NSE at longer leads\.

Table 2:Flood\-detection F1 on the CONUS land mask in the 2024 test window at 24\-, 48\-, and 72\-h leads and GloFAS return\-period thresholds from 1\.5 to 100 yr\. Higher is better; bold denotes the best score for each lead and threshold\. Results at the more uncertain 200\- and 500\-yr thresholds appear in Appendix[L](https://arxiv.org/html/2609.22702#A12)\.Figure 3:Flood\-detection margin and false\-alarm ratio on the CONUS land mask in the 2024 test window\. Panel a shows Hapi’s F1 minus the highest F1 among the three baselines at each lead and return\-period threshold\. Line color denotes forecast lead; marker shape identifies the strongest baseline at each point\. Panels b and c show the false\-alarm ratio,FPTP\+FP\\frac\{\\mathrm\{FP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FP\}\}, at the 1\.5\-yr and 100\-yr thresholds, respectively\. Color, line style, and marker encode the model, with Hapi shown by thicker lines\.Table 3:Discharge forecast skill on the CONUS land mask in the 2024 test window, pooled over space and time\. MAE and RMSE are in m3/s; NSE and KGE are dimensionless\. Bold denotes the best score per column\. Per\-cell skill is reported in Table[4](https://arxiv.org/html/2609.22702#S4.T4)\.Table 4:Per\-cell continuous skill on the CONUS land mask in the 2024 test window\. NSE, KGE, and the KGE bias ratioβ\\betaare computed independently at every grid cell, and their*medians over cells*are reported, complementing the values pooled over space and time in Table[3](https://arxiv.org/html/2609.22702#S4.T3)\. The last two columns give the percentage of cells at which Hapi’s NSE exceeds that of the comparison model and the median per\-cell difference\. Scores use 374,297 cells, representing 90\.2% of the mask, where all three sources are defined at the relevant lead; the remainder have near\-zero target variance, which makes NSE undefined\. Bold marks the best NSE and KGE at each lead\.
### 4\.3Independent Validation Against USGS Gauges

The preceding grid\-based comparisons use GloFAS reanalysis as the verification target\. Against USGS daily discharge observations, Hapi also has the highest median NSE at every lead, exceeding GloFAS at about two\-thirds of gauges with finite paired NSE values in Table[5](https://arxiv.org/html/2609.22702#S4.T5)\. This comparison covers 3,881 gauges matched by drainage area to unique grid cells and retained with paired days during the 2024 test window\. The matching and filtering protocol is given in Appendix[I](https://arxiv.org/html/2609.22702#A9)\.

Hapi’s median NSE remains near 0\.16 through 72 h, while that of persistence declines from 0\.070 to−0\.108\-0\.108\. Its median percent bias is also smaller, ranging from\+0\.5\+0\.5to\+1\.2\+1\.2%, compared with\+4\.3\+4\.3to\+4\.8\+4\.8% for GloFAS and\+5\.3\+5\.3to\+5\.5\+5\.5% for persistence\. Hapi thus retains its advantage in median NSE against an observational target, although the gauge and grid\-cell samples differ\.

This aggregate advantage coexists with substantial variation among gauges\. At 24 h, Hapi’s median NSE is 0\.160, with 10th and 90th percentiles of−4\.775\-4\.775and 0\.661; its lowest median NSE occurs in basins larger than 50,000 km2\. Figure[4](https://arxiv.org/html/2609.22702#S4.F4)illustrates this variation using gauges nearest the median Hapi 72\-h NSE within each region and drainage\-area band\. The Missouri panels show missed seasonal peaks on the Two Medicine and Milk rivers and muted variability on the North Loup River\. The South Atlantic–Gulf and Pacific Northwest panels generally track the observed variation more closely, though the Owyhee forecasts contain pulses with weak observational counterparts\. These hydrographs show why a small median bias does not imply accurate discharge timing or magnitude at every gauge\.

Table 5:Independent validation against USGS daily discharge in the 2024 test window\. Of 7,881 candidate gauges, 3,881 meet the matching and coverage criteria in Appendix[I](https://arxiv.org/html/2609.22702#A9)\. NSE, KGE, RMSE, and percent bias, PBIAS, are medians across gauges\. The final columns give the percentages with NSE above 0 and 0\.5\. Bold marks the best summary in each column\.Figure 4:Discharge at 72\-h lead against USGS observations during the 2024 test window\. Rows show Missouri, South Atlantic–Gulf, and Pacific Northwest water\-resource regions; columns show drainage areas of 500–5,000, 5,000–50,000, and at least 50,000 km2\. Within each bin, the selected gauge has Hapi NSE nearest the bin median among gauges with at least 300 paired days and finite NSE\. The Alabama River bin contains only one eligible gauge\. Panel titles give Hapi NSE against observations; selection details are in Appendix[I](https://arxiv.org/html/2609.22702#A9)\.
### 4\.4Four\-Variable and Discharge\-Only Model Comparison

We compare the available four\-variable and discharge\-only configurations at 24\-h lead to characterize their discharge skill\. The discharge\-only configuration treats the other three land\-surface variables as prescribed inputs and uses weighted MAE\. Because the configurations also differ in training schedule, input\-data revision, and loss formulation, this comparison cannot isolate the effect of multivariable prediction or its propagation during longer rollouts\.

The four\-variable configuration has lower MAE and RMSE, while leading flood\-detection F1 at six of seven thresholds in Table[6](https://arxiv.org/html/2609.22702#S4.T6)\. Discharge\-only MAE is 32% higher, at 2\.07 versus 1\.57 m3/s, and RMSE is 27% higher, at 41\.14 versus 32\.30 m3/s\. The discharge\-only configuration nevertheless has higher KGE and 100\-yr F1, so the advantage depends on the metric\. The next comparison holds the four\-variable configuration fixed to examine the contribution of loss weighting\.

Table 6:Comparison of four\-variable and discharge\-only configurations at 24\-h lead over CONUS in the 2024 test window\. Training differences are described in Section[4\.4](https://arxiv.org/html/2609.22702#S4.SS4)\. Full Hapi matches Table[3](https://arxiv.org/html/2609.22702#S4.T3)\. MAE and RMSE are in m3/s; other metrics are dimensionless\. Bold marks the higher skill in each column\.
### 4\.5Loss\-Weighting Ablation

To distinguish the effect of task scaling from the effect of adapting the scales during training, we compare three loss formulations with the same four\-variable architecture, inputs, targets, data splits, batch size, and 16\-GPU training budget\. Uniform task weights minimize the unweighted mean of the four normalized\-variable MAEs\. Fixed task scales use the Laplacian objective in Equation[7](https://arxiv.org/html/2609.22702#S3.E7), with the four log\-scales set to the values learned by the adaptive control and held constant throughout training\. Learned task scales use the same Laplacian objective but optimize the log\-scales jointly with the model parameters\. Table[7](https://arxiv.org/html/2609.22702#S4.T7)compares the three single\-step models at 24 h on the 2024 test window\.

Learned task scales yield the lowest discharge errors and highest F1 at every threshold\. Uniform weights increase MAE by 24%, from 1\.59 to 1\.97 m3/s, and reduce F1 by 4\.7–11\.5 percentage points, with the largest loss at the 100\-yr threshold\. Fixed scales bring MAE to 1\.64 m3/s, only 3% above the learned\-scale result, accounting for most of the improvement in this metric\.

Learning the task scales reduces MAE only slightly relative to fixed scales but produces larger gains in F1 at rare\-event thresholds\. At the 100\-yr threshold, fixed scales raise F1 from 0\.345 to 0\.382, whereas learning the scales raises it further to 0\.460\. The additional gain from learning the scales exceeds the gain from fixed scaling at every threshold from 10 to 100 yr, indicating that adaptation particularly benefits rare\-flood detection\.

Table 7:Loss\-weighting ablation at 24\-h lead over CONUS in the 2024 test window\. The same four\-variable architecture and data were used for single\-step training on 16 GPUs, with checkpoints selected by minimum validation loss at epochs 952–953\. Each configuration has one training run\. MAE and RMSE are in m3/s; other metrics are dimensionless\. These pretraining checkpoints precede the fine\-tuning used in Table[3](https://arxiv.org/html/2609.22702#S4.T3)\.
### 4\.6Non\-Discharge Prognostic Skill

The preceding comparisons focus on discharge, one of Hapi’s four prognostic variables\. We next evaluate the three other variables against persistence to assess their forecast skill beyond carrying the initial state forward\. Against their reanalysis targets, Hapi predicts runoff and soil wetness with lower RMSE and MAE and higher NSE than persistence at all three leads in Table[8](https://arxiv.org/html/2609.22702#S4.T8)\. Runoff NSE remains positive, from 0\.328 at 24 h to 0\.298 at 72 h, whereas carrying the initial runoff field forward gives negative NSE at every lead\.

For soil wetness and screened snow, Hapi’s advantage becomes clearer at longer leads\. Persistence is already accurate at 24 h for these slowly varying stores, but its RMSE increases more rapidly during the rollout\. By 72 h, Hapi reduces RMSE relative to persistence by 36% for soil wetness and 32% for screened snow\. Hapi also has lower MAE and higher NSE for both stores at every lead, while persistence retains higher 24\-h KGE\.

The snow comparison depends on three cells with persistent target snow accumulation among the 349,402 evaluated cells; these three cells are excluded by the screening rule in Appendix[J](https://arxiv.org/html/2609.22702#A10)\. These cells contribute about 99% of Hapi’s unscreened squared error\. Retaining them raises Hapi’s 24\-h RMSE from 2\.51 to 22\.20 mm and gives persistence lower RMSE and higher NSE at every lead\. Table[8](https://arxiv.org/html/2609.22702#S4.T8)reports both masks, showing the improvement across the screened snow field alongside the concentrated errors that reverse the full\-mask ranking\.

Table 8:Skill of the non\-discharge prognostic variables against reanalysis over the 2024 test window\. Hapi and analysis\-time persistence use identical finite cell–day pairs on each variable’s channel mask\. Runoff and snow RMSE and MAE are in mm of water equivalent; soil wetness, NSE, and KGE are dimensionless\. Bold marks the better score for each variable and lead\. The screened snow rows exclude three persistent accumulation cells using the rule in Appendix[J](https://arxiv.org/html/2609.22702#A10)\.VariableLeadModelRMSEMAENSEKGERunoff water equivalent, 418,792 cells24 hHapi1\.85800\.28500\.32810\.34740\.3474persistence2\.45002\.45000\.43630\.4363−0\.1686\-0\.16860\.415648 hHapi1\.89000\.29990\.31070\.3365persistence2\.67902\.67900\.52450\.5245−0\.3862\-0\.38620\.30390\.303972 hHapi1\.90700\.31050\.29800\.3327persistence2\.71202\.71200\.56150\.5615−0\.4197\-0\.41970\.28710\.2871Snow water equivalent, 349,402 cells24 hHapi22\.200022\.20000\.60480\.60480\.91770\.91770\.87110\.8711persistence2\.73500\.58890\.99880\.996448 hHapi25\.690025\.69000\.84760\.88970\.88970\.84540\.8454persistence4\.63101\.09801\.09800\.99640\.992572 hHapi28\.170028\.17001\.05000\.86730\.86730\.82600\.8260persistence6\.23101\.56601\.56600\.99350\.9883screened, 349,399 cells24 hHapi2\.51000\.56920\.99770\.99230\.9923persistence2\.73502\.73500\.58880\.58880\.99720\.99720\.996148 hHapi3\.46100\.80400\.99560\.9947persistence4\.62904\.62901\.09801\.09800\.99210\.99210\.99150\.991572 hHapi4\.24601\.00100\.99330\.9941persistence6\.22906\.22901\.56601\.56600\.98560\.98560\.98640\.9864Soil wetness index, 426,768 cells24 hHapi0\.00910\.00370\.99800\.99650\.9965persistence0\.00930\.00930\.00420\.00420\.99800\.99800\.998648 hHapi0\.01040\.00450\.99740\.9969persistence0\.01470\.01470\.00740\.00740\.99490\.99490\.99680\.996872 hHapi0\.01220\.00570\.99650\.9957persistence0\.01890\.01890\.01010\.01010\.99150\.99150\.99500\.9950
### 4\.7Per\-Subset Skill

For discharge, the continental F1 comparison combines rivers of different sizes and flood activity, leaving open whether Hapi’s advantage holds across the river network\. We therefore compare Hapi and RiverMamba at their nominal forecast leads on AIFAS river points, flood\-active cells, and three drainage\-area groups separated at 500 and 5,000 km2\. The first two subsets overlap the drainage groups; their definitions and cell counts are given in Appendix[L](https://arxiv.org/html/2609.22702#A12)\. RiverMamba’s latest discharge input is one day older at each nominal lead, as described in Section[4\.1](https://arxiv.org/html/2609.22702#S4.SS1)\.

Hapi has higher F1 in all 75 comparisons across five subsets, three leads, and the five thresholds from 1\.5 to 20 yr\. At the 2\-yr threshold, the largest gaps occur in the small\-river subset, where Hapi leads by 0\.185–0\.212 across the three leads in Table[9](https://arxiv.org/html/2609.22702#S4.T9)\. The four exceptions within the full 1\.5–100\-yr comparison occur on large rivers at the 50\- and 100\-yr thresholds at 48 and 72 h\. At 72 h, RiverMamba’s F1 is 0\.187 versus Hapi’s 0\.090 at 50 yr and 0\.131 versus 0\.013 at 100 yr\. The broad advantage therefore weakens for rare exceedances on large rivers at longer leads\.

Continuous discharge metrics test whether these spatial differences extend beyond threshold crossings\. Hapi has lower MAE on every subset at 24 h and on three of the five subsets at every lead, as shown in Appendix[L](https://arxiv.org/html/2609.22702#A12)\. On large rivers, RiverMamba has lower MAE from 48 h and higher NSE at 72 h, while its KGE is already higher at 24 h\.

Table 9:Per\-subset flood\-detection F1 at the 2\-yr return period for Hapi and RiverMamba, abbreviated RM, at nominal forecast leads over CONUS in the 2024 test window\. Bold marks the higher value for each subset and lead\. Appendix[L](https://arxiv.org/html/2609.22702#A12)reports all seven thresholds from 1\.5 to 100 yr\.
### 4\.8Forcing\-Perturbation Sensitivity

The preceding forecasts use ERA5\-Land forcing, whereas operational predictions depend on uncertain weather forecasts\. We therefore test how flood\-detection skill responds when precipitation departs from reanalysis\. We focus on precipitation because errors in rainfall amount directly alter the water input available for runoff generation and can shift discharge across flood thresholds\. At rollout steps 2 and 3, precipitation is multiplied byeσ​εe^\{\\sigma\\varepsilon\}, whereε\\varepsilonis a standardized Gaussian field correlated in space and across rollout steps\. This intervention preserves dry cells and scales the perturbation with rainfall intensity\. We testσ=0\\sigma=0, 0\.2, 0\.5, and 1\.0; the 24\-h forcing is unchanged, and the zero\-amplitude run recovers the unperturbed baseline\.

The log\-normal multiplier changes both rainfall variability and mean amount\. Its median is one, but its mean iseσ2/2e^\{\\sigma^\{2\}/2\}, reaching about 1\.65 atσ=1\\sigma=1\. The multiplier ranges corresponding to Gaussian perturbations within one standard deviation of zero are 0\.82–1\.22 atσ=0\.2\\sigma=0\.2and 0\.61–1\.65 atσ=0\.5\\sigma=0\.5, representing asymmetric multiplicative departures rather than symmetric percentage errors\.

F1 and NSE decrease or remain unchanged across the four tested amplitudes, with larger relative F1 losses at rare thresholds\. At 48 h andσ=1\\sigma=1, NSE changes from 0\.990 to 0\.987, while 100\-yr F1 falls from 0\.269 to 0\.039\. At the intermediate amplitudeσ=0\.5\\sigma=0\.5, the 100\-yr F1 loss is already 38%, compared with 4% at the 1\.5\-yr threshold\. The small change in pooled NSE therefore conceals substantial sensitivity in rare\-event detection\. Because these perturbations are not calibrated to a weather forecast system, operational skill still requires evaluation with archived NWP forecasts\.

Table 10:Sensitivity to precipitation perturbations at rollout steps 2–3 over CONUS in the 2024 test window\. The zero\-preserving multipliereσ​εe^\{\\sigma\\varepsilon\}uses a standardized Gaussian field smoothed with a spatial kernel of width 3 grid cells and a first\-order autoregressive temporal coefficient of 0\.8\. The log\-noise amplitudeσ\\sigmais distinct from the smoothing width\. The multiplier is not mean\-centered; its nominal mean iseσ2/2e^\{\\sigma^\{2\}/2\}\. The 24\-h forcing is unchanged, andσ=0\\sigma=0is the unperturbed baseline\.
### 4\.9Case Study: Hurricane Helene, 27 September 2024

To complement the aggregate flood\-detection evaluation with an event\-scale assessment, we examine how the unperturbed forecasts represent discharge magnitude during Hurricane Helene\. Figure[5](https://arxiv.org/html/2609.22702#S4.F5)compares forecasts valid on 27 September 2024 with GloFAS reanalysis over the southern Appalachians, following Helene’s landfall the previous day\. The 24\-, 48\-, and 72\-h forecasts were initialized on 26, 25, and 24 September, respectively\. We quantify magnitude using the high\-flow median ratio: the median forecast discharge over the 204 cells in the highest 1% of reference discharge divided by the reanalysis median over those same cells\.

The spatial median in GloFAS reanalysis over these cells is 3,547 m3/s\. Hapi’s high\-flow median ratio is 42%, 18%, and 19% at 24, 48, and 72 h, compared with 26%, 16%, and 12% for GloFAS\. A ratio of 42%, for example, means that the forecast median over the selected cells is 42% of the reanalysis median\. Absolute magnitudes remain below half of the reference median for both systems\. Hapi nevertheless reproduces a larger fraction of the high\-discharge magnitude at every lead, with its clearest advantage at 24 h\.

![Refer to caption](https://arxiv.org/html/2609.22702v1/case_study_helene.png)Figure 5:Hurricane Helene discharge forecasts valid on 27 September 2024 over the southern Appalachians, spanning 32\.5–38\.5∘N and 87\.0–78\.5∘W, on a common logarithmic scale\. The top row shows Hapi and the bottom row GloFAS forecasts, each beside the GloFAS reanalysis reference\. The 24\-, 48\-, and 72\-h forecasts were initialized on 26, 25, and 24 September, respectively\. The high\-flow median ratio is the forecast median divided by the GloFAS reanalysis median over the same 204 cells in the highest 1% of reference discharge; 100% would indicate equal medians\. The reference median is 3,547 m3s\-1\.
### 4\.10Computational Performance

Hapi completes a three\-step, four\-variable 72\-h forecast over the512×1152512\\times 1152CONUS grid with an average inference time of 0\.11 s and 1\.69 GB of peak allocated memory, as reported in Table[11](https://arxiv.org/html/2609.22702#S4.T11)\. RiverMamba produces seven discharge leads over its native 6\.22\-million\-point global domain in 13\.75 s, with 13\.33 GB of peak allocated memory before CONUS extraction\. Both measurements exclude data loading, host\-to\-device transfer, output extraction, metric calculation, and file writing\.

Normalizing execution time by spatial units and output leads gives0\.062​μ0\.062\\,\\mus per grid cell per lead for Hapi and0\.316​μ0\.316\\,\\mus per river point per lead for RiverMamba\. RiverMamba’s normalized execution time is therefore5\.1×5\.1\\timesHapi’s\. This ratio measures implementation throughput rather than performance on equivalent workloads because the systems differ in domain, spatial representation, forecast outputs, GPU memory capacity, and software environment, as detailed in Appendix[F](https://arxiv.org/html/2609.22702#A6)\.

Hapi controls computational scaling through fixed8×88\\times 8attention windows, which keep attention cost linear in the number of spatial tokens, and hierarchical downsampling, which places the widest layers on smaller grids\. The resulting arithmetic estimate is 0\.815 trillion multiply–accumulates per 24\-h pass and 2\.45 trillion for the three\-step forecast; Appendix[G](https://arxiv.org/html/2609.22702#A7)provides the derivation\. This estimate characterizes Hapi’s computational structure and is distinct from the measured execution times above\.

Table 11:Average GPU\-resident inference time over 10 repetitions at a batch size of 1 using 32\-bit floating\-point precision on one NVIDIA A100\-SXM4 GPU, with 80 GB capacity for Hapi and 40 GB for RiverMamba\. Normalized time is execution time divided by the number of spatial units and output leads\. Memory is peak allocated memory\. Data loading, transfer, CONUS extraction, metrics, and file writing are excluded\.

## 5Discussion

Hapi provides spatially distributed, multivariable hydrological forecasts over CONUS at 24–72\-h lead times while maintaining efficient continental\-scale inference\. Its clearest advantage is flood\-threshold detection: Hapi leads 20 of 21 F1 comparisons and 18 of 21 SEDI comparisons, with the advantage most consistent beyond the first forecast day\. Its ranking on continuous metrics, by contrast, depends on the scale of aggregation\. Together, the results show that learned task weighting, spatially distributed prediction, and efficient inference can be combined in a single continental model without limiting the forecast to discharge alone\.

#### Continuous versus flood\-detection skill\.

Continuous and flood\-detection metrics rank the models differently because they emphasize different parts of the discharge distribution\. Pooled RMSE, NSE, and KGE give greater influence to large rivers, whereas median per\-cell metrics describe performance at a typical location\. Hapi’s median bias ratio in Appendix[M](https://arxiv.org/html/2609.22702#A13)remains slightly below unity, indicating slight underprediction at a typical cell rather than the domain\-wide positive bias that might be inferred from its pooled KGE\. Flood\-detection F1 measures whether forecasts correctly cross locally defined warning thresholds and is therefore more closely aligned with the warning decision motivating this study\.

This threshold\-detection advantage generally becomes more pronounced as lead time increases and persistence weakens\. Hapi has the lowest false\-alarm ratio at the illustrated common and rare thresholds at every lead, although false alarms increase with lead for both Hapi and RiverMamba\. The subset analysis shows that the advantage spans most of the river network, but also identifies its boundary: RiverMamba performs better for rare exceedances on large rivers in four longer\-lead comparisons\. Thus, the aggregate result reflects broad gains rather than uniform superiority across every river size and event rarity\.

#### Spatial structure of continuous skill\.

Agreement between the per\-cell and gauge evaluations, including higher NSE than GloFAS at about two\-thirds of gauges, indicates that Hapi’s continental advantage is not solely inherited from the gridded reanalysis target\. Performance nevertheless varies substantially among locations, with weaker relative skill clustered in the upper Midwest, Northeast, and mountain West, as shown in Figure[10](https://arxiv.org/html/2609.22702#A13.F10)\. Reservoir regulation is one plausible contributor because Hapi receives no storage state, operating rule, or release schedule, whereas LISFLOOD includes a reservoir module\. A conservative name\-based proxy supports an association but not a causal explanation: the 220 identified gauges have a lower median 24\-h NSE than the remaining gauges, 0\.030 versus 0\.168, with Mann–Whitneyp=0\.004p=0\.004\. The difference is not statistically significant among basins larger than 50,000 km2, withp=0\.36p=0\.36\. Although station names do not identify all regulated rivers or account for other basin differences, the association identifies regulated rivers as a distinct challenge for spatial hydrological forecasting\.

#### Learned weighting and multivariable prediction\.

The better discharge performance of the four\-variable configuration is consistent with positive transfer among coupled hydrological variables, as reported when streamflow and water temperature are learned jointly\([Sadler et al\., 2022](https://arxiv.org/html/2609.22702#bib.bib50)\)\. However, the two configurations also differ in training schedule, input\-data revision, and loss formulation, so their comparison does not isolate the effect of joint prediction\.

The controlled loss ablation more directly identifies adaptive task balancing as an important part of the multivariable design\. Relative to uniform weighting, learned scales reduce discharge MAE by 19% and improve F1 by 4\.7–11\.5 percentage points, with the largest gain at the rarest threshold\. Fixed scales recover most of the improvement in average discharge error, whereas adapting the scales during training contributes more strongly at rare\-flood thresholds\. This distinction suggests that balancing target magnitudes and adapting their influence serve related but different roles in multivariable training\.

Skill in the three additional prognostics establishes the complementary result that the shared representation retains useful information beyond discharge\. The snow conclusion is spatially sensitive: three persistent accumulation cells dominate Hapi’s squared error and reverse the RMSE and NSE ranking on the full mask\. Outside these isolated accumulation errors, the shared model maintains skill across runoff, soil wetness, and snow throughout the rollout\.

#### Sensitivity to meteorological forcing\.

Operational hydrological forecasts depend on uncertain weather predictions, whereas the main evaluation supplies ERA5\-Land forcing over the forecast horizon\. The precipitation perturbations show why this distinction matters: pooled NSE changes little, but relative F1 losses grow at rarer thresholds\. Because the multiplier changes both rainfall variability and mean amount and is not calibrated to a forecast system, the experiment diagnoses sensitivity rather than predicting operational degradation\. The contrast between stable pooled NSE and declining rare\-threshold F1 shows that meteorological errors can affect flood detection well before they dominate aggregate continuous metrics\.

The Hurricane Helene case study provides an event\-scale view of the same forecast system\. At 24 h, Hapi reproduces 42% of the reference high\-flow median, compared with 26% for GloFAS; the corresponding ratios remain higher for Hapi at 48 and 72 h, although both systems remain below half of the reference median\. The case study complements the aggregate threshold metrics with a spatially concrete example of both Hapi’s relative advantage and the substantial magnitude error that remains during an individual extreme event\.

#### Computational efficiency\.

Shifted\-window attention and hierarchical downsampling expand spatial context while keeping attention local, and the mirrored decoder restores full\-resolution fields without global attention\. Training the final model through its two phases used approximately 830 A100\-hours\. A three\-step CONUS forecast takes an average of 0\.11 s with 1\.69 GB of peak allocated GPU memory, showing that model inference is computationally efficient once trained\.

#### Limitations\.

The present evaluation focuses on CONUS during 2024 and uses ERA5\-Land forcing over the forecast horizon\. Extending the evaluation across multiple years with archived operational NWP forecasts would test performance under realistic meteorological forcing\. The deterministic framework can be extended with calibrated ensembles or a probabilistic prediction head to characterize forecast uncertainty\. Reservoir storage and operating information could likewise be added as inputs, particularly for the regulated rivers highlighted by the gauge analysis\.

#### Operational relevance\.

Hapi combines continental coverage, multivariable prediction, strong flood\-threshold detection, and low inference cost in one forecasting framework\. The forcing experiment shows that rare\-event detection is particularly sensitive to meteorological input errors, while the spatial analysis points to regulated rivers as a distinct forecasting challenge\. Its low execution cost makes repeated forecasts with updated meteorological inputs and ensemble prediction computationally practical\.

## 6Conclusion

Hapi, a hierarchical Swin Transformer V2–based U\-Net, forecasts discharge and evolving land\-surface states together across CONUS at0\.05∘0\.05^\{\\circ\}resolution with efficient continental\-scale inference\. The encoder–decoder combines fine three\-dimensional patches, shifted\-window attention, and multiscale skip connections, while learned task weights balance discharge, runoff, snow water equivalent, and soil wetness over 24–72\-h forecasts\. Hapi leads 20 of 21 F1 lead–threshold comparisons across the four models, and its F1 advantage over RiverMamba remains significant after multiplicity correction in all 21 comparisons\. Among Hapi, GloFAS, and persistence, Hapi leads 18 of 21 SEDI comparisons, including every threshold at 48 and 72 h\. Although GloFAS leads several pooled metrics at longer leads, Hapi has the highest median per\-cell NSE at every lead\. Independent validation at 3,881 USGS gauges confirms that this local\-skill advantage persists against observations\.

The controlled ablation shows that learned task weighting reduces 24\-h discharge error and contributes most strongly to F1 at rare\-flood thresholds\. Beyond discharge, Hapi outperforms persistence in RMSE, MAE, and NSE for runoff, soil wetness, and screened snow throughout the rollout\. The sensitivity of rare\-event detection to precipitation perturbations makes forcing quality an important consideration for operational application\. On one A100 GPU, the average runtime for a four\-variable 72\-h CONUS forecast is just 0\.11 s, with 1\.69 GB of peak allocated memory\. Hapi provides an efficient framework for high\-resolution, multivariable hydrological forecasting across CONUS\.

## Acknowledgments

We gratefully acknowledge the computing resources provided on Swing, a high\-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory\.

## Funding

This work has been supported by the U\.S Department of Energy \(DOE\) Office of Cybersecurity, Energy Security, and Emergency Response \(CESER\) through the U\.S\. Department of Energy \(DOE\) contract DE\-AC02\-06CH11357\. Computing resources come from the Argonne Leadership Computing Facility, a U\.S\. Department of Energy \(DOE\) Office of Science user facility at Argonne National Laboratory\.

## Data availability

The persistent repository identifier, archived version, license, and public\-release date will be inserted here before publication\. Ground\-truth ERA5\-Land and GloFAS reanalysis inputs, needed only for retraining or a fresh inference pass, are publicly available from the Copernicus Climate Change Service at[https://cds\.climate\.copernicus\.eu](https://cds.climate.copernicus.eu/)and the GloFAS portal at[https://global\-flood\.emergency\.copernicus\.eu](https://global-flood.emergency.copernicus.eu/)\. The upstream\-drainage\-area field is distributed with the GloFAS v4\.0 auxiliary data\([ECMWF, 2023](https://arxiv.org/html/2609.22702#bib.bib4)\)\. Representative\-gauge hydrographs are given in Figure[4](https://arxiv.org/html/2609.22702#S4.F4)and per\-cell skill maps in Appendix[M](https://arxiv.org/html/2609.22702#A13); spatial flood\-risk maps at the 100\-yr return period are deferred to a future release\.

## References

- Addoret al\.\(2017\)N\. Addor, A\. J\. Newman, N\. Mizukami, and M\. P\. ClarkThe CAMELS data set: catchment attributes and meteorology for large\-sample studies\.Hydrology and Earth System Sciences21\(10\),pp\. 5293–5313\.External Links:[Document](https://dx.doi.org/10.5194/hess-21-5293-2017)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Alfieriet al\.\(2013\)L\. Alfieri, P\. Burek, E\. Dutra, B\. Krzeminski, D\. Muraro, J\. Thielen, and F\. PappenbergerGloFAS – global ensemble streamflow forecasting and flood early warning\.Hydrology and Earth System Sciences17\(3\),pp\. 1161–1175\.External Links:[Document](https://dx.doi.org/10.5194/hess-17-1161-2013)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p1.1),[§4\.1](https://arxiv.org/html/2609.22702#S4.SS1.p2.1)\.
- Becket al\.\(2020\)H\. E\. Beck, M\. Pan, P\. Lin, J\. Seibert, A\. I\. J\. M\. van Dijk, and E\. F\. WoodGlobal fully distributed parameter regionalization based on observed streamflow from 4,229 headwater catchments\.Journal of Geophysical Research: Atmospheres125\(17\),pp\. e2019JD031485\.External Links:[Document](https://dx.doi.org/10.1029/2019JD031485)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Biet al\.\(2023\)K\. Bi, L\. Xie, H\. Zhang, X\. Chen, X\. Gu, and Q\. TianAccurate medium\-range global weather forecasting with 3D neural networks\.Nature619\(7970\),pp\. 533–538\.External Links:[Document](https://dx.doi.org/10.1038/s41586-023-06185-3)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1),[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p1.1),[§3](https://arxiv.org/html/2609.22702#S3.p2.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Bindaset al\.\(2024\)T\. Bindas, W\. Tsai, J\. Liu, F\. Rahmani, D\. Feng, Y\. Bian, K\. Lawson, and C\. ShenImproving river routing using a differentiable Muskingum–Cunge model and physics\-informed machine learning\.Water Resources Research60\(1\),pp\. e2023WR035337\.External Links:[Document](https://dx.doi.org/10.1029/2023WR035337)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1)\.
- Bodnaret al\.\(2025\)C\. Bodnar, W\. P\. Bruinsma, A\. Lucic, M\. Stanley, A\. Allen, J\. Brandstetter, P\. Garvan, M\. Riechert, J\. Weyn, H\. Dong, J\. K\. Gupta, K\. Tambiratnam, A\. Archibald, C\. Wu, E\. Heider, M\. Welling, R\. E\. Turner, and P\. PerdikarisA foundation model for the Earth system\.Nature641,pp\. 1180–1187\.External Links:[Document](https://dx.doi.org/10.1038/s41586-025-09005-y)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1),[§3](https://arxiv.org/html/2609.22702#S3.p2.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Brandstetteret al\.\(2022\)J\. Brandstetter, D\. Worrall, and M\. WellingMessage passing neural PDE solvers\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p4.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Bureket al\.\(2013\)P\. Burek, J\. van der Knijff, and A\. de RooLISFLOOD – distributed water balance and flood simulation model – revised user manual 2013\.Technical reportTechnical ReportEUR 26162 EN,Publications Office of the European Union,Joint Research Centre, European Commission,Luxembourg\.Cited by:[§2\.3](https://arxiv.org/html/2609.22702#S2.SS3.p1.1),[§4\.1](https://arxiv.org/html/2609.22702#S4.SS1.p2.1)\.
- Chenet al\.\(2023\)L\. Chen, X\. Zhong, F\. Zhang, Y\. Cheng, Y\. Xu, Y\. Qi, and H\. LiFuXi: a cascade machine learning forecasting system for 15\-day global weather forecast\.npj Climate and Atmospheric Science6,pp\. 190\.External Links:[Document](https://dx.doi.org/10.1038/s41612-023-00512-1)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1),[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p1.1),[§3](https://arxiv.org/html/2609.22702#S3.p2.1)\.
- Copernicus Climate Change Service \(2019\)Copernicus Climate Change ServiceERA5\-Land hourly data from 1950 to present\.Note:Copernicus Climate Change Service \(C3S\) Climate Data Store \(CDS\)External Links:[Document](https://dx.doi.org/10.24381/cds.e2161bac),[Link](https://doi.org/10.24381/cds.e2161bac)Cited by:[Table 1](https://arxiv.org/html/2609.22702#S2.T1),[§2](https://arxiv.org/html/2609.22702#S2.p1.1)\.
- Cosgroveet al\.\(2024\)B\. Cosgrove, D\. Gochis, T\. Flowers, A\. Dugger, F\. Ogden, T\. Graziano, E\. Clark, R\. Cabell, N\. Casiday, Z\. Cui, K\. Eicher, G\. Fall, X\. Feng, K\. Fitzgerald, N\. Frazier, C\. George, R\. Gibbs, L\. Hernandez, D\. Johnson, R\. Jones, L\. Karsten, H\. Kefelegn, D\. Kitzmiller, H\. Lee, Y\. Liu, H\. Mashriqui, D\. Mattern, A\. McCluskey, J\. L\. McCreight, R\. McDaniel, A\. Midekisa, A\. Newman, L\. Pan, C\. Pham, A\. RafieeiNasab, R\. Rasmussen, L\. Read, M\. Rezaeianzadeh, F\. Salas, D\. Sang, K\. Sampson, T\. Schneider, Q\. Shi, G\. Sood, A\. Wood, W\. Wu, D\. Yates, W\. Yu, and Y\. ZhangNOAA’s National Water Model: advancing operational hydrology through continental\-scale modeling\.Journal of the American Water Resources Association60\(2\),pp\. 247–272\.External Links:[Document](https://dx.doi.org/10.1111/1752-1688.13184)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p1.1),[§4\.1](https://arxiv.org/html/2609.22702#S4.SS1.p5.1)\.
- Davison and Hinkley \(1997\)A\. C\. Davison and D\. V\. HinkleyBootstrap methods and their application\.Cambridge Series in Statistical and Probabilistic Mathematics,Cambridge University Press\.Cited by:[Appendix N](https://arxiv.org/html/2609.22702#A14.SS0.SSS0.Px5.p1.1)\.
- Demiray and Demir \(2026\)B\. Z\. Demiray and I\. DemirTowards generalized hydrological forecasting using transformer models for 120 h streamflow prediction\.Machine Learning: Earth2\(2\),pp\. 025006\.External Links:[Document](https://dx.doi.org/10.1088/3049-4753/ae81b8)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- ECMWF \(2023\)ECMWFGloFAS v4\.0 auxiliary data\.Note:[https://confluence\.ecmwf\.int/pages/viewpage\.action?pageId=242067380](https://confluence.ecmwf.int/pages/viewpage.action?pageId=242067380)Upstream area, elevation, local drainage direction, and flood\-threshold fields for GloFAS v4\.0\. Accessed 2026\-09\-03Cited by:[Appendix B](https://arxiv.org/html/2609.22702#A2.SS0.SSS0.Px1.p1.1),[§2\.3](https://arxiv.org/html/2609.22702#S2.SS3.p1.1),[Table 1](https://arxiv.org/html/2609.22702#S2.T1),[Data availability](https://arxiv.org/html/2609.22702#Sx3.p1.1)\.
- Fanget al\.\(2024\)Y\. Fang, Q\. Sun, X\. Wang, T\. Huang, X\. Wang, and Y\. CaoEVA\-02: a visual representation for neon genesis\.Image and Vision Computing149,pp\. 105171\.External Links:[Document](https://dx.doi.org/10.1016/j.imavis.2024.105171)Cited by:[§3](https://arxiv.org/html/2609.22702#S3.p1.1)\.
- Fenget al\.\(2022\)D\. Feng, J\. Liu, K\. Lawson, and C\. ShenDifferentiable, learnable, regionalized process\-based models with multiphysical outputs can approach state\-of\-the\-art hydrologic prediction accuracy\.Water Resources Research58\(10\),pp\. e2022WR032404\.External Links:[Document](https://dx.doi.org/10.1029/2022WR032404)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1)\.
- Ferro and Stephenson \(2011\)C\. A\. T\. Ferro and D\. B\. StephensonExtremal dependence indices: improved verification measures for deterministic forecasts of rare binary events\.Weather and Forecasting26\(5\),pp\. 699–713\.Cited by:[Appendix H](https://arxiv.org/html/2609.22702#A8.p2.4)\.
- Gauchet al\.\(2021\)M\. Gauch, F\. Kratzert, D\. Klotz, G\. Nearing, J\. Lin, and S\. HochreiterRainfall–runoff prediction at multiple timescales with a single Long Short\-Term Memory network\.Hydrology and Earth System Sciences25\(4\),pp\. 2045–2062\.External Links:[Document](https://dx.doi.org/10.5194/hess-25-2045-2021)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Gu and Dao \(2024\)A\. Gu and T\. DaoMamba: linear\-time sequence modeling with selective state spaces\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[§3](https://arxiv.org/html/2609.22702#S3.p2.1)\.
- Guptaet al\.\(2009\)H\. V\. Gupta, H\. Kling, K\. K\. Yilmaz, and G\. F\. MartinezDecomposition of the mean squared error and NSE performance criteria: implications for improving hydrological modelling\.Journal of Hydrology377\(1–2\),pp\. 80–91\.External Links:[Document](https://dx.doi.org/10.1016/j.jhydrol.2009.08.003)Cited by:[Appendix H](https://arxiv.org/html/2609.22702#A8.p1.2)\.
- Harriganet al\.\(2020\)S\. Harrigan, E\. Zsoter, L\. Alfieri, C\. Prudhomme, P\. Salamon, F\. Wetterhall, C\. Barnard, H\. Cloke, and F\. PappenbergerGloFAS\-ERA5 operational global river discharge reanalysis 1979–present\.Earth System Science Data12\(3\),pp\. 2043–2060\.External Links:[Document](https://dx.doi.org/10.5194/essd-12-2043-2020)Cited by:[§2\.3](https://arxiv.org/html/2609.22702#S2.SS3.p1.1),[§2](https://arxiv.org/html/2609.22702#S2.p1.1)\.
- Hatanpääet al\.\(2025\)V\. Hatanpää, E\. Ku, J\. Stock, M\. Emani, S\. Foreman, C\. Jung, S\. Madireddy, T\. Nguyen, V\. Sastry, R\. A\. O\. Sinurat, H\. Zheng, S\. Wheeler, T\. Arcomano, V\. Vishwanath, and R\. KotamarthiAERIS: Argonne earth systems model for reliable and skillful predictions\.InProceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis,SC ’25,pp\. 72–85\.External Links:[Document](https://dx.doi.org/10.1145/3712285.3772094),2509\.13523Cited by:[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Heoet al\.\(2024\)B\. Heo, S\. Park, D\. Han, and S\. YunRotary position embedding for vision transformer\.InEuropean Conference on Computer Vision \(ECCV\),External Links:2403\.13298Cited by:[§3\.1](https://arxiv.org/html/2609.22702#S3.SS1.SSS0.Px1.p1.1)\.
- Hoedtet al\.\(2021\)P\. Hoedt, F\. Kratzert, D\. Klotz, C\. Halmich, M\. Holzleitner, G\. S\. Nearing, S\. Hochreiter, and G\. KlambauerMC\-LSTM: mass\-conserving LSTM\.InProceedings of the 38th International Conference on Machine Learning \(ICML\),PMLR, Vol\.139,pp\. 4275–4286\.Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1)\.
- Holm \(1979\)S\. HolmA simple sequentially rejective multiple test procedure\.Scandinavian Journal of Statistics6\(2\),pp\. 65–70\.Cited by:[Appendix N](https://arxiv.org/html/2609.22702#A14.SS0.SSS0.Px5.p1.1)\.
- Jaegleet al\.\(2022\)A\. Jaegle, S\. Borgeaud, J\. Alayrac, C\. Doersch, C\. Ionescu, D\. Ding, S\. Koppula, D\. Zoran, A\. Brock, E\. Shelhamer, O\. Hénaff, M\. M\. Botvinick, A\. Zisserman, O\. Vinyals, and J\. CarreiraPerceiver IO: a general architecture for structured inputs and outputs\.InInternational Conference on Learning Representations \(ICLR\),External Links:2107\.14795Cited by:[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Kendallet al\.\(2018\)A\. Kendall, Y\. Gal, and R\. CipollaMulti\-task learning using uncertainty to weigh losses for scene geometry and semantics\.InProceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition \(CVPR\),pp\. 7482–7491\.External Links:[Document](https://dx.doi.org/10.1109/CVPR.2018.00781)Cited by:[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p1.1),[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p3.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Kratzertet al\.\(2022\)F\. Kratzert, M\. Gauch, G\. Nearing, and D\. KlotzNeuralHydrology — A Python library for deep learning research in hydrology\.Journal of Open Source Software7\(71\),pp\. 4050\.External Links:[Document](https://dx.doi.org/10.21105/joss.04050)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Kratzertet al\.\(2019\)F\. Kratzert, D\. Klotz, G\. Shalev, G\. Klambauer, S\. Hochreiter, and G\. NearingTowards learning universal, regional, and local hydrological behaviors via machine learning applied to large\-sample datasets\.Hydrology and Earth System Sciences23\(12\),pp\. 5089–5110\.External Links:[Document](https://dx.doi.org/10.5194/hess-23-5089-2019)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Kratzertet al\.\(2023\)F\. Kratzert, G\. Nearing, N\. Addor, T\. Erickson, M\. Gauch, O\. Gilon, L\. Gudmundsson, A\. Hassidim, D\. Klotz, S\. Nevo, G\. Shalev, and Y\. MatiasCaravan – A global community dataset for large\-sample hydrology\.Scientific Data10,pp\. 61\.External Links:[Document](https://dx.doi.org/10.1038/s41597-023-01975-w)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Künsch \(1989\)H\. R\. KünschThe jackknife and the bootstrap for general stationary observations\.The Annals of Statistics17\(3\),pp\. 1217–1241\.Cited by:[Appendix N](https://arxiv.org/html/2609.22702#A14.p1.1)\.
- Lamet al\.\(2023\)R\. Lam, A\. Sanchez\-Gonzalez, M\. Willson, P\. Wirnsberger, M\. Fortunato, F\. Alet, S\. Ravuri, T\. Ewalds, Z\. Eaton\-Rosen, W\. Hu, A\. Merose, S\. Hoyer, G\. Holland, O\. Vinyals, J\. Stott, A\. Pritzel, S\. Mohamed, and P\. BattagliaLearning skillful medium\-range global weather forecasting\.Science382\(6677\),pp\. 1416–1421\.External Links:[Document](https://dx.doi.org/10.1126/science.adi2336)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1),[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p1.1),[§3](https://arxiv.org/html/2609.22702#S3.p2.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Liuet al\.\(2022\)Z\. Liu, H\. Hu, Y\. Lin, Z\. Yao, Z\. Xie, Y\. Wei, J\. Ning, Y\. Cao, Z\. Zhang, L\. Dong, F\. Wei, and B\. GuoSwin Transformer V2: scaling up capacity and resolution\.InProceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition \(CVPR\),pp\. 11999–12009\.External Links:[Document](https://dx.doi.org/10.1109/CVPR52688.2022.01170)Cited by:[§3](https://arxiv.org/html/2609.22702#S3.p1.1)\.
- Nearinget al\.\(2024\)G\. Nearing, D\. Cohen, V\. Dube, M\. Gauch, O\. Gilon, S\. Harrigan, A\. Hassidim, D\. Klotz, F\. Kratzert, A\. Metzger, S\. Nevo, F\. Pappenberger, C\. Prudhomme, G\. Shalev, S\. Shenzis, T\. Y\. Tekalign, D\. Weitzner, and Y\. MatiasGlobal prediction of extreme floods in ungauged watersheds\.Nature627\(8004\),pp\. 559–563\.External Links:[Document](https://dx.doi.org/10.1038/s41586-024-07145-1)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Nearinget al\.\(2021\)G\. S\. Nearing, F\. Kratzert, A\. K\. Sampson, C\. S\. Pelissier, D\. Klotz, J\. M\. Frame, C\. Prieto, and H\. V\. GuptaWhat role does hydrological science play in the age of machine learning?\.Water Resources Research57\(3\),pp\. e2020WR028091\.External Links:[Document](https://dx.doi.org/10.1029/2020WR028091)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Nevoet al\.\(2022\)S\. Nevo, E\. Morin, A\. G\. Rosenthal, A\. Metzger, C\. Barshai, D\. Weitzner, D\. Voloshin, F\. Kratzert, G\. Elidan, G\. Dror, G\. Begelman, G\. Nearing, G\. Shalev, H\. Noga, I\. Shavitt, L\. Yuklea, M\. Royz, N\. Giladi, N\. Peled Levi, O\. Reich, O\. Gilon, R\. Maor, S\. Timnat, T\. Shechter, V\. Anisimov, Y\. Gigi, Y\. Levin, Z\. Moshe, Z\. Ben\-Haim, A\. Hassidim, and Y\. MatiasFlood forecasting with machine learning models in an operational framework\.Hydrology and Earth System Sciences26\(15\),pp\. 4013–4032\.External Links:[Document](https://dx.doi.org/10.5194/hess-26-4013-2022)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Newmanet al\.\(2015\)A\. J\. Newman, M\. P\. Clark, K\. Sampson, A\. Wood, L\. E\. Hay, A\. Bock, R\. J\. Viger, D\. Blodgett, L\. Brekke, J\. R\. Arnold, T\. Hopson, and Q\. DuanDevelopment of a large\-sample watershed\-scale hydrometeorological data set for the contiguous USA: data set characteristics and assessment of regional variability in hydrologic model performance\.Hydrology and Earth System Sciences19\(1\),pp\. 209–223\.External Links:[Document](https://dx.doi.org/10.5194/hess-19-209-2015)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Nguyenet al\.\(2023\)T\. Nguyen, J\. Brandstetter, A\. Kapoor, J\. K\. Gupta, and A\. GroverClimaX: a foundation model for weather and climate\.InProceedings of the 40th International Conference on Machine Learning \(ICML\),PMLR, Vol\.202,pp\. 25904–25938\.Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Nguyenet al\.\(2024\)T\. Nguyen, R\. Shah, H\. Bansal, T\. Arcomano, R\. Maulik, V\. Kotamarthi, I\. Foster, S\. Madireddy, and A\. GroverScaling transformer neural networks for skillful and reliable medium\-range weather forecasting\.InAdvances in Neural Information Processing Systems \(NeurIPS\),External Links:2312\.03876Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Priceet al\.\(2025\)I\. Price, A\. Sanchez\-Gonzalez, F\. Alet, T\. R\. Andersson, A\. El\-Kadi, D\. Masters, T\. Ewalds, J\. Stott, S\. Mohamed, P\. Battaglia, R\. Lam, and M\. WillsonProbabilistic weather forecasting with machine learning\.Nature637,pp\. 84–90\.External Links:[Document](https://dx.doi.org/10.1038/s41586-024-08252-9)Cited by:[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p1.1),[§3](https://arxiv.org/html/2609.22702#S3.p3.1)\.
- Qiuet al\.\(2025\)Z\. Qiu, Z\. Wang, B\. Zheng, Z\. Huang, K\. Wen, S\. Yang, R\. Men, L\. Yu, F\. Huang, S\. Huang, D\. Liu, J\. Zhou, and J\. LinGated attention for large language models: non\-linearity, sparsity, and attention\-sink\-free\.InAdvances in Neural Information Processing Systems,External Links:2505\.06708Cited by:[§3\.1](https://arxiv.org/html/2609.22702#S3.SS1.SSS0.Px2.p1.1),[§3](https://arxiv.org/html/2609.22702#S3.p1.1)\.
- Ronnebergeret al\.\(2015\)O\. Ronneberger, P\. Fischer, and T\. BroxU\-Net: convolutional networks for biomedical image segmentation\.InMedical Image Computing and Computer\-Assisted Intervention – MICCAI 2015,pp\. 234–241\.External Links:[Document](https://dx.doi.org/10.1007/978-3-319-24574-4%5F28)Cited by:[§3](https://arxiv.org/html/2609.22702#S3.p1.1)\.
- Sadleret al\.\(2022\)J\. M\. Sadler, A\. P\. Appling, J\. S\. Read, S\. K\. Oliver, X\. Jia, J\. A\. Zwart, and V\. KumarMulti\-task deep learning of daily streamflow and water temperature\.Water Resources Research58\(4\),pp\. e2021WR030138\.External Links:[Document](https://dx.doi.org/10.1029/2021WR030138)Cited by:[§5](https://arxiv.org/html/2609.22702#S5.SS0.SSS0.Px3.p1.1)\.
- Schmudeet al\.\(2024\)J\. Schmude, S\. Roy, W\. Trojak, J\. Jakubik, D\. S\. Civitarese, S\. Singh, J\. Kuehnert, K\. Ankur, A\. Gupta, C\. E\. Phillips, R\. Kienzler, D\. Szwarcman, V\. Gaur, R\. Shinde, R\. Lal, A\. Da Silva, J\. L\. G\. Diaz, A\. Jones, S\. Pfreundschuh, A\. Lin, A\. Sheshadri, U\. Nair, V\. Anantharaj, H\. Hamann, C\. Watson, M\. Maskey, T\. J\. Lee, J\. B\. Moreno, and R\. RamachandranPrithvi WxC: foundation model for weather and climate\.arXiv preprint arXiv:2409\.13598\.Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1)\.
- Shams Eddinet al\.\(2025\)M\. H\. Shams Eddin, Y\. Zhang, S\. Kollet, and J\. GallRiverMamba: a state space model for global river discharge and flood forecasting\.InAdvances in Neural Information Processing Systems \(NeurIPS\),External Links:2505\.22535Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p3.1),[§3\.3](https://arxiv.org/html/2609.22702#S3.SS3.p1.1),[§3](https://arxiv.org/html/2609.22702#S3.p2.1),[§4\.1](https://arxiv.org/html/2609.22702#S4.SS1.p3.1)\.
- Songet al\.\(2025\)W\. Song, M\. Guan, and D\. YuSwinFlood: a hybrid CNN–Swin Transformer model for rapid spatiotemporal flood simulation\.Journal of Hydrology,pp\. 133280\.External Links:[Document](https://dx.doi.org/10.1016/j.jhydrol.2025.133280)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p3.1)\.
- Suet al\.\(2024\)J\. Su, M\. Ahmed, Y\. Lu, S\. Pan, W\. Bo, and Y\. LiuRoFormer: enhanced transformer with rotary position embedding\.Neurocomputing568,pp\. 127063\.External Links:[Document](https://dx.doi.org/10.1016/j.neucom.2023.127063)Cited by:[§3\.1](https://arxiv.org/html/2609.22702#S3.SS1.SSS0.Px1.p1.1),[§3](https://arxiv.org/html/2609.22702#S3.p1.1)\.
- Sun and Sun \(2026\)A\. Y\. Sun and A\. A\. SunZero\-shot forecasting of streamflow using time series foundation models: are we there yet?\.Machine Learning: Earth2\(1\),pp\. 010501\.External Links:[Document](https://dx.doi.org/10.1088/3049-4753/ae4982)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p2.1)\.
- Tellmanet al\.\(2021\)B\. Tellman, J\. A\. Sullivan, C\. Kuhn, A\. J\. Kettner, C\. S\. Doyle, G\. R\. Brakenridge, T\. A\. Erickson, and D\. A\. SlaybackSatellite imaging reveals increased proportion of population exposed to floods\.Nature596\(7870\),pp\. 80–86\.External Links:[Document](https://dx.doi.org/10.1038/s41586-021-03695-w)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p1.1)\.
- Vaughanet al\.\(2025\)A\. Vaughan, S\. Markou, W\. Tebbutt, J\. Requeima, W\. P\. Bruinsma, T\. R\. Andersson, M\. Herzog, N\. D\. Lane, M\. Chantry, J\. S\. Hosking, and R\. E\. TurnerEnd\-to\-end data\-driven weather prediction\.Nature\.External Links:[Document](https://dx.doi.org/10.1038/s41586-025-08897-0)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1)\.
- Wesselkampet al\.\(2025\)M\. Wesselkamp, M\. Chantry, E\. Pinnington, M\. Choulga, S\. Boussetta, M\. Kalweit, J\. Boedecker, C\. F\. Dormann, F\. Pappenberger, and G\. BalsamoAdvances in land surface forecasting: a comparison of LSTM, gradient boosting, and feed\-forward neural networks as prognostic state emulators in a case study with ecLand\.Geoscientific Model Development18,pp\. 921–949\.External Links:[Document](https://dx.doi.org/10.5194/gmd-18-921-2025)Cited by:[§1](https://arxiv.org/html/2609.22702#S1.p4.1)\.

## Appendix APer\-Variable Normalization

The 12 input variables have very different scales and distributions, so preprocessing is specified separately for each variable\. We use the following transformations for the distributions shown in Figure[6](https://arxiv.org/html/2609.22702#A1.F6):

- •log1p:x↦log⁡\(1\+x\)x\\mapsto\\log\(1\+x\), without subsequent standardization\. Used for variables with strong right skew and dynamic ranges spanning several orders of magnitude \(e\.g\., river discharge has skewness\>30\>30and a max/median ratio above10510^\{5\}\)\.
- •min\-max:x↦\(x−xmin\)/\(xmax−xmin\)x\\mapsto\(x\-x\_\{\\min\}\)/\(x\_\{\\max\}\-x\_\{\\min\}\)\. Used for volumetric soil water\.
- •z\-score:x↦\(x−μ\)/σx\\mapsto\(x\-\\mu\)/\\sigma\. Used for variables that are close to Gaussian \(e\.g\., 2 m temperature, solar radiation\)\.
- •identity:x↦xx\\mapsto x\. Used for the soil wetness index, which already lies in\[0,1\]\[0,1\]\.

The normalization choice and key statistics for each variable are provided in Table[12](https://arxiv.org/html/2609.22702#A1.T12)\. Total evaporation retains ERA5\-Land’s downward\-positive flux convention and is standardized without a logarithmic transform\.

Table 12:Per\-variable normalization choices and descriptive distributional statistics\. The model’s normalization parameters were estimated from the 2005–2022 training period only\. The descriptive statistics were computed over the full 2005–2024 record \(every tenth day\) on the0\.05∘0\.05^\{\\circ\}CONUS grid\. The basis is the 419,682 cells where 2 m temperature is defined, the widest of the per\-variable masks\. Flood detection is scored on the narrower 414,795\-cell threshold\-valid mask described in Appendix[B](https://arxiv.org/html/2609.22702#A2)\. Skewness is the standard biased estimator\. Zero % counts exact zeros, and Max/Med is the maximum divided by the absolute median\.†Snow depth has a zero median, so its Max/Med ratio uses the median of nonzero values\.Figure 6:Raw distributions of the 12 sampled input variables over the temperature\-defined CONUS land mask\.\. This illustration uses a shorter window than the descriptive statistics in Table[12](https://arxiv.org/html/2609.22702#A1.T12)\. Each panel shows skewness, zero fraction, and the maximum divided by the absolute median\. Med0denotes the median of nonzero values when the overall median is zero\. Variables normalized with log1p use a logarithmicxx\-axis to show the tail structure\. Red and orange lines mark the median and mean\. Large separation indicates high skewness\.
## Appendix BPhysical Masking

The CONUS grid at0\.05∘0\.05^\{\\circ\}resolution contains ocean pixels, the Great Lakes, and zero\-padded regions outside the domain boundary\. We compute a separate mask𝐦c∈\{0,1\}H×W\\mathbf\{m\}\_\{c\}\\in\\\{0,1\\\}^\{H\\times W\}for each channelcc:

mc,h,w=𝟏\[∃ts\.t\.xc,t,h,w≠0\]\.m\_\{c,h,w\}=\\mathbf\{1\}\\bigl\[\\exists\\,t\\text\{ s\.t\. \}x\_\{c,t,h,w\}\\neq 0\\bigr\]\.\(9\)A pixel is included if it is nonzero at any time step in the training period\. Different variables produce different masks\. River discharge, for example, is nonzero only along river channels, while 2 m temperature covers all land pixels\. Static fields \(DEM, upstream area\) use a simpler rule: the mask is 1 wherever the field value is nonzero\.

These masks serve two purposes:

1. 1\.Loss masking: masked\-out pixels do not contribute to the training loss, so the model does not spend capacity on zero\-padded regions\.
2. 2\.Output zeroing: predictions are multiplied by the mask before being fed back during autoregressive rollout, keeping zero\-padded regions exactly at zero\.

#### Preparation of the flood\-detection thresholds\.

The GloFAS v4\.0 return\-level fields\[[ECMWF, 2023](https://arxiv.org/html/2609.22702#bib.bib4)\]are distributed globally at0\.05∘0\.05^\{\\circ\}on a3000×72003000\\times 7200grid whose cell centers coincide with those of our512×1152512\\times 1152CONUS grid\. Preparation therefore requires only coordinate subsetting: we select the CONUS cells by exact coordinate match\. We verified that the maximum coordinate offset between target and source is00, that the subset values are bit\-identical to the source, and that the missing\-data pattern is preserved exactly\. We avoid interpolation because it would smooth the return\-level field, including the local maxima that define the scored events\. Of the589,824589\{,\}824cells in the CONUS box,426,768426\{,\}768\(72\.4%72\.4\\%\) have a defined return level\. Intersecting these cells with the416,298416\{,\}298\-cell land mask leaves the414,795414\{,\}795cells used throughout\. Cells left undefined in the source are assigned\+∞\+\\inftyat load time so that no event can ever be triggered on them\.

## Appendix CPhysical Bounds on Prognostic Feedback

Several prognostic variables \(river discharge, runoff, and snow water equivalent\) are physically nonnegative and carry large fractions of exact zeros\. The soil wetness index is bounded in\[0,1\]\[0,1\]by construction\. During autoregressive rollout, predictions that stray outside those ranges are fed back as input and accumulate across steps\. Each prognostic channel therefore declares optional physical bounds in the configuration, and both sides are enforced on the prediction before it is fed back\.

The constants are computed in physical units and then mapped through the channel’s own normalizationfcf\_\{c\}from Appendix[A](https://arxiv.org/html/2609.22702#A1), usinglog⁡\(1\+x\)\\log\(1\+x\)for discharge and runoff\. Using the standardization formula\(x−μc\)/σc\(x\-\\mu\_\{c\}\)/\\sigma\_\{c\}would place the thresholds incorrectly for these variables\.

The floor uses snapping because the targets have real mass at the floor\. For a channel with declared minimumlcl\_\{c\}and smallest training valuezcz\_\{c\}strictly above it, the cutoff is the midpoint between the two,

τc=fc​\(12​\(lc\+zc\)\),\\tau\_\{c\}=f\_\{c\}\\\!\\left\(\\tfrac\{1\}\{2\}\\left\(l\_\{c\}\+z\_\{c\}\\right\)\\right\),\(10\)and any prediction belowτc\\tau\_\{c\}is replaced byfc​\(lc\)f\_\{c\}\(l\_\{c\}\)\. For the nonnegative variableslc=0l\_\{c\}=0, so the cutoff in physical units is half the smallest nonzero value observed in training\. Predictions below this cutoff are set to zero\. For a channel with no observed value above its lower bound, the rule reduces to standard clipping\.

The upper bound uses standard clipping: a prediction above a declared maximumucu\_\{c\}is set tofc​\(uc\)f\_\{c\}\(u\_\{c\}\), while all other predictions remain unchanged\. Snapping up would be wrong wherever no target reaches the ceiling, as with the soil wetness index, which is bounded at 1 but attains 0\.988 at most in the data\.

Both rules are applied at intermediate rollout steps only, in training and at inference alike, so the two see the same feedback\. Neither is applied at the final step, where the loss is taken, so the loss sees raw output\. At intermediate steps the clip does zero the gradient on the cells it moves, which is a cost we accept in exchange for training and inference following the same rollout\.

## Appendix DTwo\-Phase Training

Training uses the objective in Section[3\.3](https://arxiv.org/html/2609.22702#S3.SS3)under the following two\-phase schedule\.

#### Phase 1: single\-step pretraining\.

The model is trained at a single autoregressive step \(J=1J=1\) under Eq\.[7](https://arxiv.org/html/2609.22702#S3.E7), on a cosine learning\-rate schedule with warm restarts\. This fixes the one\-step operator before the compounding errors of rollout are introduced\.

#### Phase 2: multi\-step fine\-tuning\.

The pretrained weights are then fine\-tuned at randomized rollout depthJ∼Uniform⁡\{1,…,Jmax\}J\\sim\\mathrm\{Uniform\}\\\{1,\\dots,J\_\{\\max\}\\\}withJmax=3J\_\{\\max\}=3, backpropagating through the whole rollout, as described in Section[3\.3](https://arxiv.org/html/2609.22702#S3.SS3)\. Intermediate predictions are held to their declared physical bounds before feedback, as specified in Appendix[C](https://arxiv.org/html/2609.22702#A3), so the model sees its own accumulated errors under the same rule used at inference\.

Phase 2 uses two constraints\. Phase 2 uses a base learning rate of1×10−51\{\\times\}10^\{\-5\}, compared with5×10−45\{\\times\}10^\{\-4\}in Phase 1, and keeps the learned log\-scalessis\_\{i\}fixed\. If left trainable, the log\-scales can reduce the rollout objective by increasingsis\_\{i\}for variables with accumulating multi\-step error and incurring the corresponding regularization penalty\. In Phase 1, trainable log\-scales adjust the relative contributions of the four variables to the objective\. During Phase 2, the same adjustment could lower the objective without improving rollout predictions, so the log\-scales are held fixed\.

## Appendix EArchitecture Hyperparameter Sweep

We ran a Bayesian hyperparameter sweep over 224 configurations \(206 completed\) using Weights & Biases, with a budget of 20 epochs per run\. The search space is specified in Table[13](https://arxiv.org/html/2609.22702#A5.T13)\.

Table 13:Hyperparameter search space\.Two architectural regimes stood out\.Shallow encoders with fine\-grained tokenization\(p=2p\{=\}2, depths\(1,1,4\)\(1,1,4\)or\(1,2,4\)\(1,2,4\),d=192d\{=\}192–240240\) dominated the top 10 runs by validation loss and converged fastest\.Deep encoders with coarser tokenization\(p=4p\{=\}4, depths\(3,3,8\)\(3,3,8\),d=192d\{=\}192–240240\) formed a competitive second tier that converged more slowly but used less compute per step\. Window size had little effect on final performance\. Embedding dimensions of 192 and 240 consistently beat 48 and 96, showing that model capacity matters at this spatial resolution\.

## Appendix FImplementation Details

#### Optimizer\.

We train on 16 NVIDIA A100 GPUs distributed across two eight\-GPU nodes\. We use AdamW withβ1=0\.9\\beta\_\{1\}\{=\}0\.9,β2=0\.999\\beta\_\{2\}\{=\}0\.999, and two parameter groups\. Model parameters use base learning rate5×10−45\{\\times\}10^\{\-4\}during single\-step pretraining and1×10−51\{\\times\}10^\{\-5\}during multi\-step fine\-tuning\. A 2\-epoch linear warmup is used for pretraining\. Cosine cycles span 20 epochs in pretraining and 500 epochs in multi\-step fine\-tuning, with weight decay1×10−81\{\\times\}10^\{\-8\}\. Pretraining uses batch size 6 per GPU\. Multi\-step fine\-tuning uses batch size 4 per GPU andJmax=3J\_\{\\max\}\{=\}3\. The learning rate is scaled byB×world\_size\\sqrt\{B\\times\\text\{world\\\_size\}\}\(square\-root scaling for AdamW\)\. No gradient clipping is applied\. The log\-scale parameters\{sc\}\\\{s\_\{c\}\\\}in the adaptive task\-weighting loss use a fixed learning rate of10−310^\{\-3\}with no weight decay, keeping them tied to per\-variable error magnitudes throughout the model’s annealing schedule\. Under multi\-step fine\-tuning the log\-scales are frozen, requiring the model to reduce the multi\-step rollout error directly\. Per\-variable prior weights are uniform \(wi=1w\_\{i\}\{=\}1\) in every reported run\. The objective supports non\-uniform priors, but they are not used\. Initial log\-scales aresi\(0\)=\(−4\.24,−3\.17,−3\.00,−4\.48\)s\_\{i\}^\{\(0\)\}\{=\}\(\-4\.24,\\,\-3\.17,\\,\-3\.00,\\,\-4\.48\)for \(discharge, surface runoff, snow water equivalent, soil wetness index\), seeded from per\-variable error magnitudes observed in a short warmup run\. The evaluated checkpoint retainssi=\(−4\.4134,−3\.7854,−4\.4622,−6\.2489\)s\_\{i\}\{=\}\(\-4\.4134,\\,\-3\.7854,\\,\-4\.4622,\\,\-6\.2489\), giving inverse scalesexp⁡\(−si\)\\exp\(\-s\_\{i\}\)of\(82\.6,44\.1,86\.7,517\.4\)\(82\.6,\\,44\.1,\\,86\.7,\\,517\.4\)\. These log\-scales are inherited unchanged from the pretraining checkpoint\. The learned weighting is far from uniform and would be hard to guess\. Relative to discharge, the soil wetness index is weighted6\.3×6\.3\\timesand surface runoff0\.53×0\.53\\times\. All reported configurations use one training realization and identical evaluation procedures\. The codebase includes an optional Binned Spectral Power auxiliary loss for spectral\-fidelity fine\-tuning, but it is disabled for all analyzed runs\.

#### Data loading\.

Each input variable is stored as a separate NumPy memory\-mapped file with shape\(Ttotal,H,W\)\(T\_\{\\text\{total\}\},H,W\)\. DataLoaders use 2 workers with a prefetch factor of 1 and persistent workers to reduce startup cost\. Per\-variable normalization statistics are computed once and cached to disk\.

#### Distributed training\.

The training script auto\-detects the distributed environment \(SLURM, torchrun, or JSM\) and falls back to single\-GPU training if initialization fails\. Logging and checkpointing run only on the main process\. Per\-node file operations, such as memory\-mapped file copies to local scratch, run once per node\.

#### Hapi benchmark protocol\.

We benchmarked the model of record on one NVIDIA A100\-SXM4\-80GB with driver 580\.95\.05 using PyTorch 2\.10\.0\+cu128\. Inference used FP32 with autocast disabled, matching the production configuration\. We enabledtorch\.compile\(mode="default"\), which optimizes the model’s computation graph before repeated execution\. CUDA events measured the complete rollout, including forward evaluation, masking, physical bounds, prognostic feedback, and forcing injection\. Measurements began after data loading and host\-to\-device transfer and ended before metric computation and file writing\. Each configuration used 3 warm\-up rollouts and 10 timed repetitions\. Peak allocated memory fromtorch\.cuda\.max\_memory\_allocatedincludes the model and device\-resident inputs\.

#### Hapi benchmark results\.

The 43\.96 M\-parameter model requires one forward pass for a 24\-h forecast and three autoregressive passes for a 72\-h forecast\. At batch 1, average inference time increases from 36\.6 ms at 24 h to 109\.5 ms at 72 h, with interquartile ranges of 36\.5–36\.6 and 109\.4–109\.5 ms, respectively\. Peak allocated memory increases from 1\.63 to 1\.69 GB\. At the production batch size of 4, the corresponding average inference times are 139\.2 and 418\.5 ms, with peak allocated memory of 4\.49 and 4\.74 GB\. The interquartile range remained below 0\.3% of the average in all four configurations\.

#### RiverMamba timing protocol\.

We measured the public RiverMamba checkpoint with 4\.38 M parameters in FP32 on one NVIDIA A100\-SXM4\-40GB using PyTorch 2\.4\.1\+cu124\. As for Hapi, graph optimization was enabled before timing:torch\.compile\(mode="default"\)transforms supported operations into optimized GPU graphs, while RiverMamba’s custom selective\-scan and causal\-convolution operators remain eager operations\. RiverMamba produces seven leads jointly for 6,221,926 global river points, processed in the same 20 chunks used for its native global predictions\. We summed CUDA\-event execution times across all chunks for each complete global forecast\. After 3 warm\-up forecasts, 10 repetitions gave an average of 13\.75 s, with an interquartile range of 13\.748–13\.749 s, and peak allocated memory of 13\.33 GB\. Input preparation, host\-to\-device transfer, CONUS extraction, metric calculation, and file writing were excluded from both systems’ measurements\.

The two runs used the same A100\-SXM4 compute architecture but different memory capacities and software versions\. Because Hapi predicts four variables on a regular CONUS grid while RiverMamba predicts discharge on global river points, the raw times do not represent equal workloads\. Table[11](https://arxiv.org/html/2609.22702#S4.T11)therefore also reports execution time divided by spatial units and output leads\. This normalization measures implementation throughput\. It does not remove differences in model structure, spatial representation, or forecast construction\.

## Appendix GHapi Arithmetic Estimate

We count one multiplication and its accumulation as one multiply–accumulate \(MAC\)\. Counting them separately gives twice as many floating\-point operations\. ForNNtokens, channel widthCC, MLP expansionr=4r=4, and window widthw=8w=8, a Hapi block requires

MACblock=13​N​C2\+2​N​w2​C\.\\mathrm\{MAC\}\_\{\\mathrm\{block\}\}=13NC^\{2\}\+2Nw^\{2\}C\.\(11\)The first term covers the query–key–value–gate, output, and MLP projections\. The second covers the two attention matrix products\. We omit normalization, activation, softmax, rotary embedding, masking, and data movement\. Patch size 2 converts the512×1152512\\times 1152input into 147,456 tokens\. Two merging steps reduce this count to 36,864 and 9,216\. Table[14](https://arxiv.org/html/2609.22702#A7.T14)gives the resulting component totals\. The complete encoder–decoder requires 814\.9 billion MACs per 24\-h pass and 2\.45 trillion MACs for the three\-pass forecast\. The two attention matrix products contribute 2\.5% of the estimated MACs within Hapi’s transformer blocks\. Their cost is linear in the number of spatial tokens because the window size is fixed\.

Table 14:Estimated Hapi arithmetic for one 24\-h forward pass\. The three backbone rows include both encoder and decoder blocks at each resolution\.ComponentTokensWidthMACs, billionBackbone, 2 blocks147,456192148\.6Backbone, 4 blocks36,864384289\.9Backbone, 4 blocks9,216768286\.3Patch embedding147,4561922\.7Merge, unmerge, and skip projections——65\.2Output expansion and projection——22\.2One 24\-h pass814\.9Three\-pass 72\-h forecast2,444\.8
## Appendix HEvaluation Metric Definitions

Letyiy\_\{i\}andy^i\\hat\{y\}\_\{i\}denote the reference and forecast values on theNNpaired finite samples used in a calculation, and let an overbar denote the sample mean\. The continuous metrics are

MAE\\displaystyle\\mathrm\{MAE\}=1N​∑i=1N\|y^i−yi\|,\\displaystyle=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\\lvert\\hat\{y\}\_\{i\}\-y\_\{i\}\\rvert,\(12\)RMSE\\displaystyle\\mathrm\{RMSE\}=\[1N​∑i=1N\(y^i−yi\)2\]1/2,\\displaystyle=\\left\[\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\(\\hat\{y\}\_\{i\}\-y\_\{i\}\)^\{2\}\\right\]^\{1/2\},\(13\)NSE\\displaystyle\\mathrm\{NSE\}=1−∑i=1N\(y^i−yi\)2∑i=1N\(yi−y¯\)2,\\displaystyle=1\-\\frac\{\\sum\_\{i=1\}^\{N\}\(\\hat\{y\}\_\{i\}\-y\_\{i\}\)^\{2\}\}\{\\sum\_\{i=1\}^\{N\}\(y\_\{i\}\-\\bar\{y\}\)^\{2\}\},\(14\)KGE\\displaystyle\\mathrm\{KGE\}=1−\[\(r−1\)2\+\(α−1\)2\+\(β−1\)2\]1/2,\\displaystyle=1\-\\left\[\(r\-1\)^\{2\}\+\(\\alpha\-1\)^\{2\}\+\(\\beta\-1\)^\{2\}\\right\]^\{1/2\},\(15\)whererris the Pearson correlation betweeny^\\hat\{y\}andyy,α=sy^/sy\\alpha=s\_\{\\hat\{y\}\}/s\_\{y\}is the variability ratio, andβ=y^¯/y¯\\beta=\\bar\{\\hat\{y\}\}/\\bar\{y\}is the bias ratio\[[Gupta et al\., 2009](https://arxiv.org/html/2609.22702#bib.bib42)\]\. For the gauge evaluation, percent bias is

PBIAS=100​∑i=1N\(y^i−yi\)∑i=1Nyi\.\\mathrm\{PBIAS\}=100\\frac\{\\sum\_\{i=1\}^\{N\}\(\\hat\{y\}\_\{i\}\-y\_\{i\}\)\}\{\\sum\_\{i=1\}^\{N\}y\_\{i\}\}\.\(16\)NSE is undefined when the reference variance is zero\. KGE is undefined when the reference mean or either standard deviation needed by its ratios and correlation is zero\. PBIAS is undefined when the summed reference discharge is zero\.

For flood detection, a forecast and reference are each classified by whether they exceed the same local return\-level threshold\. The resulting contingency table contains true positives \(TP\), false positives \(FP\), false negatives \(FN\), and true negatives \(TN\)\. We define

precision\\displaystyle\\mathrm\{precision\}=TPTP\+FP,\\displaystyle=\\frac\{\\mathrm\{TP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FP\}\},recall=H\\displaystyle\\mathrm\{recall\}=H=TPTP\+FN,\\displaystyle=\\frac\{\\mathrm\{TP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FN\}\},\(17\)F1\\displaystyle\\mathrm\{F1\}=2​TP2​TP\+FP\+FN,\\displaystyle=\\frac\{2\\,\\mathrm\{TP\}\}\{2\\,\\mathrm\{TP\}\+\\mathrm\{FP\}\+\\mathrm\{FN\}\},CSI\\displaystyle\\mathrm\{CSI\}=TPTP\+FP\+FN,\\displaystyle=\\frac\{\\mathrm\{TP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FP\}\+\\mathrm\{FN\}\},\(18\)FAR\\displaystyle\\mathrm\{FAR\}=FPTP\+FP,\\displaystyle=\\frac\{\\mathrm\{FP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FP\}\},F\\displaystyle F=FPFP\+TN\.\\displaystyle=\\frac\{\\mathrm\{FP\}\}\{\\mathrm\{FP\}\+\\mathrm\{TN\}\}\.\(19\)Here,HHis the hit rate, also called the probability of detection \(or recall\), FAR is the false\-alarm ratio, andFFis the false\-alarm rate \(also called the probability of false detection\)\. The frequency bias shown by the rays in Figure[13](https://arxiv.org/html/2609.22702#A15.F13)is

B=TP\+FPTP\+FN=recallprecision\.B=\\frac\{\\mathrm\{TP\}\+\\mathrm\{FP\}\}\{\\mathrm\{TP\}\+\\mathrm\{FN\}\}=\\frac\{\\mathrm\{recall\}\}\{\\mathrm\{precision\}\}\.\(20\)The symmetric extremal dependence index is

SEDI=ln⁡F−ln⁡H−ln⁡\(1−F\)\+ln⁡\(1−H\)ln⁡F\+ln⁡H\+ln⁡\(1−F\)\+ln⁡\(1−H\)\\mathrm\{SEDI\}=\\frac\{\\ln F\-\\ln H\-\\ln\(1\-F\)\+\\ln\(1\-H\)\}\{\\ln F\+\\ln H\+\\ln\(1\-F\)\+\\ln\(1\-H\)\}\(21\)for0<H<10<H<1and0<F<10<F<1\[[Ferro and Stephenson, 2011](https://arxiv.org/html/2609.22702#bib.bib49)\]\. We report SEDI as undefined when either rate is 0 or 1 rather than clipping the rate before evaluation\.

Unless stated otherwise, the main discharge metrics pool all valid cell–day pairs over the CONUS mask before each score is formed\. Likewise, pooled flood\-detection scores sum TP, FP, FN, and TN over valid cells and days before applying the equations above\. The per\-cell analysis in Appendix[M](https://arxiv.org/html/2609.22702#A13)instead computes a score independently at each cell and summarizes the resulting spatial distribution\. The gauge analysis computes each score independently at each retained gauge and reports medians across gauges\.

## Appendix IGauge Matching and Hydrograph Selection

We match each USGS gauge to the cell within±2\\pm 2grid cells whose GloFAS v4\.0 upstream drainage area best matches the reported gauge basin area\. This search reduces mismatches between adjacent tributaries and main\-stem reaches that coordinate\-only matching can introduce\. Matches differing by more than a factor of two in drainage area are rejected, leaving 7,838 of 7,881 candidates with a median area ratio of 1\.01\. Where multiple gauges claim one cell, we retain the closest drainage\-area match before checking observation availability, leaving 7,256 unique cells\. Observations are available for 3,976 retained gauges\. Requiring at least 200 paired days leaves 3,881 scored gauges\. Paired finite NSE values are available at 3,869 gauges at each lead\.

Hydrograph selection uses Hapi’s 72\-h NSE and a stricter requirement of at least 300 paired days, leaving 3,710 gauges with finite NSE before regional and drainage\-area selection\. We take the three most populated eligible water\-resource regions and choose the gauge nearest each region–area bin’s median NSE\. The drainage\-area bands are 500–5,000, 5,000–50,000, and at least 50,000 km2\. Eligible bin counts from smallest to largest are 267, 149, and 60 in Missouri, 304, 117, and 1 in South Atlantic–Gulf, and 247, 89, and 27 in the Pacific Northwest\. The selected gauge approximates a bin’s median skill\. The single\-gauge Alabama River bin cannot represent within\-bin variability\.

## Appendix JNon\-Discharge Evaluation Masks and Snow Screening

Each variable is evaluated on its own prediction mask, reconstructed from cells with nonzero temporal variance in the 24\-h prediction\. Masked outputs remain constant at the inverse transform of zero, which is zero for the three non\-discharge prognostic variables\. Hapi and analysis\-time persistence are scored on identical finite cell–day pairs at each lead\. The original GloFAS GRIB headers specify kg m\-2for both runoff and snow water equivalent\. Preprocessing retains these numerical values, which are equivalent to millimetres of water equivalent\. Runoff and snow errors are therefore reported in millimetres without numerical rescaling\. Soil wetness is dimensionless\.

The snow sensitivity analysis excludes cells whose minimum finite target snow water equivalent exceeds 100 mm over the available evaluation window\. This target\-based rule removes three of 349,402 snow cells and does not use prediction error to select them\. The scoring script applies the rule to the evaluation target, so the screened scores are a sensitivity analysis conditional on that mask\. The excluded cells retain large snow accumulations in the Cascades and account for about 99% of unscreened squared error\. They remain in training and in the full\-mask evaluation\. Table[8](https://arxiv.org/html/2609.22702#S4.T8)reports both masks\.

## Appendix KExtended\-Lead RiverMamba Diagnostics

Hapi is trained and evaluated at lead times of 24, 48, and 72 hours, so the head\-to\-head comparisons in Section[4](https://arxiv.org/html/2609.22702#S4)stop at 72 hours\. Table[15](https://arxiv.org/html/2609.22702#A11.T15)reports RiverMamba’s flood\-detection metrics at 4–7\-day leads on the same CONUS test window and return\-period thresholds\. The numbers come from the public RiverMamba checkpoint evaluated under our pipeline\. They are*not*a head\-to\-head comparison with Hapi at those leads\.

F1 generally declines as lead time increases beyond the leads in Table[2](https://arxiv.org/html/2609.22702#S4.T2)\. At the 1\.5\-yr threshold the CONUS\-pooled series runs0\.728→0\.697→0\.5920\.728\\to 0\.697\\to 0\.592over the 24\-, 48\-, and 72\-h leads of the main text, then0\.580→0\.566→0\.571→0\.4770\.580\\to 0\.566\\to 0\.571\\to 0\.477at 96, 120, 144, and 168 hours\. F1 declines monotonically with return period at every lead\.

Table 15:RiverMamba flood\-detection F1 at 4–7\-day leads on the CONUS land mask, 2024 test window\. Hapi is not evaluated at these leads\. CONUS\-pooled from the per\-subset TP/FP/FN of the same evaluation used throughout the paper\.
## Appendix LPer\-Subset Flood Detection Skill

The per\-lead grids for all seven primary thresholds extend the per\-subset comparison in Section[4\.7](https://arxiv.org/html/2609.22702#S4.SS7)\. Tables[16](https://arxiv.org/html/2609.22702#A12.T16)and[18](https://arxiv.org/html/2609.22702#A12.T18)report flood\-detection F1 and continuous discharge skill for every subset at the nominal 24\-, 48\-, and 72\-h leads reported by each system\.

Both models’ pooled F1 aggregates TP, FP, and FN across the cells of a subset before the ratio is formed\. Hapi’s numbers come from the per\-leaddaily\_discharge\.nc, and RiverMamba’s come from its per\-lead evaluation logs\. The AIFAS subset contains 146,265 cells co\-located with operational AIFAS diagnostic river points\. The flood\-active subset contains 135,523 cells exceeding the 2\-yr threshold at least once in 2024\. Upstream drainage area partitions all 414,795 valid CONUS cells into small rivers below 500 km2, medium rivers from 500 to below 5,000 km2, and large rivers of at least 5,000 km2, containing 338,893, 51,744, and 24,158 cells, respectively\. The AIFAS and flood\-active subsets overlap these drainage groups\.

Table 16:Per\-subset flood\-detection F1 at all seven primary return\-period thresholds over CONUS in the 2024 test window\. Hapi and RiverMamba, abbreviated RM, use their nominal forecast leads\. RiverMamba’s latest discharge input is one day older\. Bold marks the higher value using unrounded scores\.#### Extreme\-tail return periods \(200\- and 500\-yr\)\.

The GloFAS v4\.0 return levels are obtained from a Gumbel fit to a 44\-year record of annual maxima from 1979–2022, as described in Section[2\.3](https://arxiv.org/html/2609.22702#S2.SS3)\. The 100\-yr level therefore extrapolates to2\.3×2\.3\\timesthe record length, which is routine in design\-flood practice\. The 200\- and 500\-yr levels reach4\.5×4\.5\\timesand11×11\\timesthe fitted record length and require substantially longer extrapolation\. The main\-text F1 results in Tables[2](https://arxiv.org/html/2609.22702#S4.T2)and[6](https://arxiv.org/html/2609.22702#S4.T6)are therefore restricted to return periods≤100\\leq 100yr\. Table[17](https://arxiv.org/html/2609.22702#A12.T17)provides CONUS\-aggregated results at the more uncertain 200\- and 500\-yr levels for all four sources at 24, 48, and 72 h leads\. Persistence leads both extreme rows at 24\-h lead, driven by the tiny event count per cell \(median≤2\\leq 2over the test window on the CONUS land mask\)\. Hapi leads at both levels at 48\- and 72\-h lead, although the difference at the 500\-yr level and 48 h is only0\.0010\.001\. Hapi outperforms GloFAS at every lead\. We report these extreme\-tail rankings descriptively without tests of statistical significance\.

Table 17:CONUS\-aggregated flood\-detection F1 at the 200\- and 500\-yr return levels\. These levels extrapolate the GloFAS Gumbel fit to4\.5×4\.5\\timesand11×11\\timesits 44\-year record\. The rankings are descriptive and are not accompanied by significance tests\.Table 18:Per\-subset continuous discharge skill, pooled over paired cell–day values within each subset, at nominal forecast leads over CONUS in the 2024 test window\. MAE is in m3/s\. NSE and KGE are dimensionless\. Bold marks lower MAE or higher NSE and KGE for each subset and lead, using values before rounding to two decimal places\.Per\-cell distributions\.The pooled F1 numbers in Table[16](https://arxiv.org/html/2609.22702#A12.T16)aggregate TP, FP, and FN over every cell in a subset and report a single number, which the heaviest\-event cells dominate\. Figures[7](https://arxiv.org/html/2609.22702#A12.F7)and[8](https://arxiv.org/html/2609.22702#A12.F8)show the corresponding per\-cell distributions for Hapi, persistence, GloFAS, and RiverMamba\. At each skill threshold, the vertical axis gives the fraction of cells attaining that value or higher\. The NSE figures display thresholds from−1\-1to 1 without clipping the underlying values\. The F1 figures use the 2\-yr return\-period threshold and cells with≥\\geq1 event over the test window\. The curves use every finite cell value and compare Hapi and RiverMamba at their nominal leads\. RiverMamba provides leads 1–7, and the figures show leads 1–3\. The figures report per\-cell distributions, whereas Table[16](https://arxiv.org/html/2609.22702#A12.T16)reports ratios pooled across cells\.

The pattern is consistent across all three subset groups: Hapi’s curve generally lies above the three baselines and declines least over 24–72 h, while persistence declines most rapidly and GloFAS and RiverMamba are intermediate over most of the skill range\. Per\-cell F1 retains a visible mass at 0 \(cells where the model never aligns with a target event\) and at the small\-integer rationals\{1/3,1/2,2/3,1\}\\\{1/3,1/2,2/3,1\\\}\(cells with 1–3 events where partial agreement is possible\)\. The steps are intrinsic to event sparsity at fine spatial resolution\.

Per\-cell F1 is sensitive to timing at this resolution: for a single\-event cell, a predicted exceedance displaced by one day gives F1=0=0, since TP=0=0and FP=FN=11\. Most cells at the 2\-yr threshold see only one or two events over the test window, so the initial drop at F1=0=0reflects timing sensitivity at sparse event counts\. The resulting mass at zero occurs for all four sources and does not by itself imply complete model failure\. The pooled F1 in Table[16](https://arxiv.org/html/2609.22702#A12.T16)avoids this floor because TP, FP, and FN are accumulated across cells before the F1 ratio is computed\. Per\-cell F1 at the other eight return periods \(1\.5, 5, 10, 20, 50, 100, 200, 500 yr\) is in the supplementary file set and shows the same qualitative ordering across models\.

Figure 7:Complementary empirical distributions of per\-cell NSE for four forecast sources across five spatial subsets \(columns\) and nominal 24\-, 48\-, and 72\-h leads \(rows\)\. Each curve gives the fraction of cells attaining at least the indicated NSE and uses every finite cell value\. The horizontal axis is displayed from−1\-1to 1\. Values below−1\-1remain in the empirical sample and contribute to the fraction below the plotted range\.Figure 8:Complementary empirical distributions of per\-cell F1 at the 2\-yr return\-period threshold for four forecast sources across five spatial subsets \(columns\) and nominal 24\-, 48\-, and 72\-h leads \(rows\)\. Each curve gives the fraction of event\-containing cells attaining at least the indicated F1 and uses every finite cell value\. Hapi and RiverMamba are compared at their nominal reported leads\.

## Appendix MPer\-Cell Continuous Skill

The continuous metrics in Table[3](https://arxiv.org/html/2609.22702#S4.T3)are pooled over space and time, using all paired cell–day values in each calculation\. On a domain spanning five orders of magnitude of discharge, pooling gives high\-discharge rivers greater influence, particularly on RMSE and NSE through their squared\-error terms\. Typical\-cell skill is quantified by computing NSE, KGE, and the KGE bias ratioβ\\betaindependently at every cell and taking the median, as reported in Table[4](https://arxiv.org/html/2609.22702#S4.T4)\. A cell is included when all three sources are defined at that lead and the cell lies inside the CONUS mask: 374,297 of 414,795 cells, 90\.2%\. The excluded 9\.8% have near\-zero target variance over the test window, which makes the NSE denominator collapse\.

Pooled, GloFAS overtakes Hapi from 48\-h lead onward on NSE and KGE\. Per cell, Hapi has the highest median NSE at every lead \(0\.902 / 0\.823 / 0\.776 against 0\.698 / 0\.637 / 0\.559 for GloFAS\) and the highest median KGE at 48 and 72 h\. At 72\-h lead Hapi’s per\-cell NSE exceeds GloFAS’s at 79\.2% of cells with a median difference of\+0\.142\+0\.142, and persistence at 92\.3% of cells with a median difference of\+0\.363\+0\.363\. Pooling conceals a broad advantage at typical cells behind the behavior of a small number of very large ones\. The pooled metrics characterize skill over all cell–day values, whereas the per\-cell medians characterize skill at a typical location\.

The per\-cell medianβ\\betais 0\.973, 0\.967 and 0\.962 at 24, 48 and 72 h, indicating slight underprediction at a typical cell despite the pooled result\. Persistence has a median bias ratio ofβ=1\.000\\beta=1\.000to three decimal places at every lead\.

Figure[9](https://arxiv.org/html/2609.22702#A13.F9)maps the per\-cell scores and Figure[10](https://arxiv.org/html/2609.22702#A13.F10)the Hapi\-minus\-GloFAS difference\. The difference map shows that the residual GloFAS advantage clusters in the upper Midwest, the Northeast and parts of the mountain West\.

![Refer to caption](https://arxiv.org/html/2609.22702v1/percell_nse_maps.png)Figure 9:Per\-cell NSE over the 2024 test window for each model and lead, CONUS land mask\. Cells where the score is undefined \(near\-zero target variance\) are left blank\. The pooled values in Table[3](https://arxiv.org/html/2609.22702#S4.T3)are dominated by high\-discharge cells and do not reflect the typical cell shown here\.![Refer to caption](https://arxiv.org/html/2609.22702v1/percell_nse_diff.png)Figure 10:Per\-cell NSE difference, Hapi minus GloFAS, by lead\. Positive \(blue\) favors Hapi\. The remaining GloFAS\-favored cells form coherent regional clusters\. Possible explanations of these regional differences are discussed in Section[5](https://arxiv.org/html/2609.22702#S5)\.
## Appendix NPaired Block Bootstrap and Multiplicity Correction

The flood\-detection numbers in Table[2](https://arxiv.org/html/2609.22702#S4.T2)come from one 2024 test window whose daily values are strongly autocorrelated, so any statement about which model ranks higher needs an uncertainty estimate attached to it\. We use a moving\-block bootstrap\[[Künsch, 1989](https://arxiv.org/html/2609.22702#bib.bib46)\]over the time axis\.

#### Resampling scheme\.

Within each ofR=10,000R=10\{,\}000replicates we draw⌈n/L⌉\\lceil n/L\\rceilblocks ofL=30L=30consecutive days with replacement, truncate tonndays, and recompute the pooled contingency table\. As specified in Section[4\.1](https://arxiv.org/html/2609.22702#S4.SS1), RiverMamba’s latest discharge input is one day older than Hapi’s at the same nominal lead\. The same block draw is applied to every model within a replicate, making the difference between two models the resampled quantity\. The scheme preserves two statistical properties\. Whole days are resampled, preserving within\-day spatial dependence without treating the414,795414\{,\}795spatially clustered cells as independent observations\. Within each replicate, F1 and SEDI are recomputed from the pooled resampled contingency counts\.

#### Paired differences and marginal intervals\.

Every model is scored on the same days, and their errors are strongly positively correlated because all models encounter the same sequence of hydrological conditions\. Pairing accounts for shared day\-to\-day variation and reduces the variance of the difference relative to treating the model estimates as independent\. Requiring the two marginal 95% intervals to be disjoint is more conservative than testing the paired difference directly\. For the persistence and GloFAS comparisons, the marginal intervals overlap in every 72\-h cell, whereas the paired test separates 20 of 28 cells\. The paired test is therefore the appropriate inference for these correlated forecasts\.

#### Block length\.

The block must be long enough to preserve temporal dependence in the resampled series, which here is the daily count of cells exceeding the threshold\. That dependence is not a single number: it falls steeply with return period, because frequent exceedances track slow catchment wetness while rare ones track individual storms\. The observed series’ lag\-1 autocorrelation and integrated autocorrelation timeτint=1\+2​∑kρk\\tau\_\{\\mathrm\{int\}\}=1\+2\\sum\_\{k\}\\rho\_\{k\}, truncated at the first non\-positiveρk\\rho\_\{k\}, appear in Table[19](https://arxiv.org/html/2609.22702#A14.T19)\.

Table 19:Temporal dependence of the observed daily exceedance count on the CONUS mask, 2024 test window, by return\-period threshold\. Computed at 24\-h lead\. The series is the common verification target, so the other leads agree to within four cells per day\.τint\\tau\_\{\\mathrm\{int\}\}is in days\.No single block length matches all seven columns\. We setL=30L=30days near the upper end of the estimated dependence range and assess sensitivity to this choice below\. At this block length, the 50\- and 100\-yr F1 and SEDI differences at 72\-h lead against persistence and GloFAS are not significant after Holm correction\. Their nonsignificance does not establish equivalence\.

#### Sensitivity to the block length\.

Re\-running the whole procedure at block lengths from 7 to 60 days leaves the conclusion in place, as shown in Table[20](https://arxiv.org/html/2609.22702#A14.T20)\. The surviving count moves over a narrow range with no trend across 7–45 days\. At a block length of 60 days, each replicate draws five blocks and six blocks respectively, with the final block truncated as needed\.

Table 20:Cells surviving Holm correction \(of 21 per family\) as the bootstrap block length varies,R=10,000R=10\{,\}000throughout\. The 30\-day column is the configuration reported in Figure[11](https://arxiv.org/html/2609.22702#A14.F11)\. Reproduced bybootstrap\_block\_sensitivity\.py, which also emits Table[19](https://arxiv.org/html/2609.22702#A14.T19)\.
#### pp\-values and multiplicity\.

Each cell’spp\-value is the two\-sided percentile bootstrappp: the proportion of replicates whose difference falls on the opposite side of zero from the point estimate, doubled, with the\(1\+k\)/\(R\+1\)\(1\+k\)/\(R\+1\)correction of[Davison and Hinkley \[1997\]](https://arxiv.org/html/2609.22702#bib.bib48)so that no cell can reportp=0p=0\. The smallest attainable raw value is therefore2/\(R\+1\)2/\(R\+1\)\. Because 21 \(lead×\\timesreturn\-period\) cells are inspected at once, we apply a Holm step\-down correction\[[Holm, 1979](https://arxiv.org/html/2609.22702#bib.bib47)\]across each 21\-cell family, separately per comparison and per metric, for six families and 126 cells in total\. We define six families by crossing the two prespecified metrics \(F1 and SEDI\) with the three model comparisons \(Hapi–persistence, Hapi–GloFAS, and Hapi–RiverMamba\)\. Within each family, Holm correction covers all leads and return periods\. Holm holds under arbitrary dependence between the tests, which is what this application needs: the seven return periods within a lead are nested, since every 100\-yr exceedance is also a 1\.5\-yr exceedance\. After correction the smallest attainable adjusted value is42/\(R\+1\)42/\(R\+1\)\. AtR=1,000R=1\{,\}000, this bound is0\.0420\.042, placing a comparison with no sign reversals just belowα=0\.05\\alpha=0\.05because of Monte Carlo resolution\. IncreasingRRto10,00010\{,\}000reduces the bound to0\.00420\.0042\. IncreasingRRfrom 1,000 to 10,000 changed few uncorrected decisions but lowered the minimum attainable adjustedpp\-value\.

#### Result\.

The underlying bootstrap output also reports the paired point difference and its 2\.5th–97\.5th percentile interval for every cell\. The intervals quantify effect\-size uncertainty, whereas Figure[11](https://arxiv.org/html/2609.22702#A14.F11)emphasizes the point differences and multiplicity\-adjusted decisions\. At 24\-h lead Hapi’s F1 advantage over the GloFAS forecast is significant after correction at*every*return period from 1\.5 to 100 yr\. Hapi’s F1 advantage over RiverMamba survives correction in all 21 cells\. The corresponding SEDI difference survives in six cells, at the 1\.5\- and 2\-yr thresholds at each lead\. Across all 126 cells, 88 survive correction, distributed 29 / 30 / 29 over the 24\-, 48\- and 72\-h leads\. None of the 126 comparisons shows a significant loss for Hapi after Holm correction\. The seven negative point differences comprise six SEDI comparisons and one F1 comparison\.

![Refer to caption](https://arxiv.org/html/2609.22702v1/bootstrap_paired_heatmap.png)Figure 11:Paired moving\-block\-bootstrap differences in flood\-detection skill over the 2024 CONUS test window\. Values are Hapi minus the indicated comparator at each nominal forecast lead and return\-period threshold\. Positive values favor Hapi\. Panel a shows F1 differences and panel b SEDI differences\. Bold values and dots indicatep<0\.05p<0\.05after Holm correction across the 21 lead–return\-period cells within each comparison and metric\. RiverMamba’s latest discharge input is one day older than Hapi’s at the same nominal lead\. The corresponding 2\.5th–97\.5th percentile intervals are retained in the source\-data file\.

## Appendix OBase\-Rate\-Invariant Skill and the Precision–Recall Decomposition

F1 depends on the event base rate, which falls by orders of magnitude as the return\-period threshold rises\. Consequently, F1 supports comparisons among models at a fixed threshold but not comparisons of skill across thresholds\. For cross\-threshold comparisons, we use SEDI because it remains non\-degenerate as the event base rate tends to zero\.

Figure[12](https://arxiv.org/html/2609.22702#A15.F12)reports SEDI at the same 21 lead–threshold combinations as F1\. Hapi’s 72\-h SEDI decreases from 0\.931 at the 1\.5\-yr threshold to 0\.724 at the 100\-yr threshold, indicating lower skill for rarer events under this base\-rate\-invariant measure\. Among Hapi, GloFAS, and persistence, Hapi has the highest SEDI in 18 of the 21 comparisons, including all seven at both 48\- and 72\-h lead\. The three exceptions are all at 24\-h lead, at the 20\-, 50\- and 100\-yr thresholds, where persistence leads\.

The same contingency tables yield the precision, recall, and critical\-success index \(CSI\) values in Table[21](https://arxiv.org/html/2609.22702#A15.T21), with their geometric relationship shown in Figure[13](https://arxiv.org/html/2609.22702#A15.F13)\. Figure[13](https://arxiv.org/html/2609.22702#A15.F13)favors Hapi most clearly at 48\- and 72\-h lead, where it lies above and to the right of both baselines at every threshold and therefore occupies higher CSI contours\. At 24 h, Hapi retains higher precision, but persistence has higher recall at the rarest thresholds, producing a trade\-off rather than uniform dominance\. Hapi’s advantage over the GloFAS forecast at 72\-h lead is larger in precision than in recall at every return period: at the 5\-yr threshold it gains\+0\.191\+0\.191in precision and\+0\.106\+0\.106in recall, and at the 100\-yr threshold\+0\.081\+0\.081and\+0\.067\+0\.067, respectively\. Hapi has higher recall than GloFAS at every lead–threshold combination\.

At 72\-h lead and the 1\.5\-yr threshold, its precision exceeds GloFAS by0\.1000\.100, its recall by0\.0610\.061, and its CSI by0\.1040\.104\. Hapi also has higher CSI at 48 h \(0\.7130\.713versus0\.5970\.597\)\. GloFAS does not exceed Hapi on CSI in any comparison in Table[21](https://arxiv.org/html/2609.22702#A15.T21)\. Persistence does so in one cell only, at 24\-h lead and the 100\-yr threshold, and by0\.0010\.001\(0\.3070\.307against0\.3060\.306\), consistent with its known short\-lead advantage\. Hapi’s precision exceeds its recall in every cell but one, corresponding to fewer predicted than observed exceedances\. The exception is 72\-h lead at the 100\-yr threshold, where recall edges ahead by0\.0010\.001\(0\.2090\.209against0\.2080\.208\)\. Persistence’s precision and recall differ by at most0\.0010\.001over the test window\.

![Refer to caption](https://arxiv.org/html/2609.22702v1/sedi_heatmap.png)Figure 12:Symmetric extremal dependence index \(SEDI\) for flood detection on the CONUS land mask over the 2024 test window\. SEDI is invariant to the event base rate and ranges from−1\-1to 1, with 0 indicating no skill\. Rows group persistence, GloFAS, and Hapi at 24\-, 48\-, and 72\-h lead\. Columns give the GloFAS return\-period thresholds\. Exact SEDI values are printed in each cell\. Bold values and dots mark the best model within each lead–return\-period comparison\.Table 21:Precision, recall, and critical success index \(CSI\) for flood detection at GloFAS return\-period thresholds, CONUS land mask, 2024 test window\. Metric definitions are given in Appendix[H](https://arxiv.org/html/2609.22702#A8)\. Recall exceeding precision indicates overforecasting of exceedances\. The reverse indicates underforecasting\. Bold marks the highest value for each lead and return\-period threshold\.Figure 13:Performance diagram for flood detection over the CONUS land mask in the 2024 test window\. The horizontal axis shows success ratio \(precision\), and the vertical axis shows probability of detection \(recall\)\. Performance improves toward the upper\-right corner, and the shaded contours show CSI\. Dashed rays show frequency bias, with values above and below one indicating overforecasting and underforecasting, respectively\. Marker size decreases from the 1\.5\-yr to the 100\-yr return\-period threshold\.