A Dataset and Model for Imputing Water Surface Elevation on a Large and Extremely Sparse Spatiotemporal Graph Ruben Cartuyvels1∗ , Karim Douch2, 7 , Gabriele Bertoli3, 4 , Mounia El Baz1 , Artemis Vrettou7 , Sébastien Lefèvre5, 6 , Diego Fernandez Prieto2
arXiv:2609.11580v1 [cs.LG] 10 Sep 2026
1
European Space Agency, Φ-lab, 2 European Space Agency, Science Hub, 3 University of Florence, 4 Imperial College London, 5 Université Bretagne Sud, IRISA, 6 University of Tromsø, 7 Serco Italia SpA, Rome, Italy
Abstract Continuous monitoring of water surface elevation across river networks is critical for flood forecasting, water resource management, and understanding the global water cycle. Yet, the scarcity of in situ gauges across much of the globe constrains the development of reliable modeling frameworks. Satellite altimetry has the potential to alleviate this problem but its use is currently hindered by sparse temporal coverage. To this end, we introduce AmazonWSE, a dataset for training and evaluating large-scale spatiotemporal graph imputation methods that integrates processed satellite altimetry measurements from a range of sources, including the recent wide-swath SWOT sensor. The dataset covers approximately 19K river sections in the Amazon river basin over 10 years (2016–2026), with in situ gauges held out for evaluation. Besides contributing a novel real-world use case with the potential for societal impact, AmazonWSE introduces significant technical challenges: with fewer than 1% of sections observed per day, the dataset is far sparser than existing imputation benchmarks, and its directed acyclic river topology is both structurally different from and larger than graphs in existing datasets. We show that prior spatiotemporal graph imputation methods are not adapted to this topology, scale and sparsity, and propose a simple bidirectional selective state space model that outperforms them by sampling connected subgraphs and flattening space and time into a single token sequence with topology-aware positional encodings. Compared to the state-of-the-art published method for SWOT-based WSE densification, which integrates statistics with physical modeling, our model reduces RMSE against in situ gauges by 18-39%, while producing predictions for every river section rather than only those with sufficient nearby satellite coverage.
Introduction Monitoring water surface elevation (WSE) across river networks is essential for flood forecasting, water resource management, and understanding the global water cycle. The Amazon basin, the largest river system on Earth, is of particular importance because of its role in global climate regulation and biodiversity (Fassoni-Andrade et al. 2021). Observations of WSE come from diverse sensors that each cover the spatiotemporal domain only partially. In situ flow gauges, such as those operated by Brazil’s National Water and Sanitation agency (ANA), provide sub-daily measurements at fixed locations but cover unevenly a small fraction ∗
[email protected] This is a preprint.
Figure 1: Example time series and local graph topology in AmazonWSE. Locations can be observed by zero, one (such as the location observed by SWOT) or more (such as the location observed by ICESat-2 and S3A) sources. The satellite sources measure with reference to the approximate sea level, while the in situ gauges measure deviations from a chosen location-specific reference level.
of the river network. Satellite altimetry missions have measured WSE globally since the 1990s, but these instruments observe a given location only every 10–91 days along narrow ground tracks (Abdalati et al. 2010; Normandin et al. 2018). The Surface Water and Ocean Topography (SWOT) mission, launched in 2022, provides for the first time wideswath observations of river WSE with a 21-day repeat cycle (Biancamaria, Lettenmaier, and Pavelsky 2016), yielding unprecedented spatial coverage but a short temporal record. We synthesize these sensor observations into a benchmark dataset for spatiotemporal graph imputation of daily WSE that observes the Amazon river network. The spatiotemporal graph formed by the Amazon and its tributaries is extremely large with over 19K nodes (or reaches; river sections of ∼812 km as defined in (Altenau et al. 2021)) and highly sparse: on any given day, fewer than 1% of nodes carry an observation. Spatiotemporal graph neural networks (STGNN) for imputation often take full graphs as input, which for large graphs is resource intensive. Consequently, current methods are evaluated on small graphs of ∼1,000 nodes or less. Fur-
Dataset
Nodes
Obs.
Missing
Graph
Source
Nodes
# obs.
Period
% obs.
AQI-36 METR-LA PEMS-BAY CausalRivers LamaH-CE LargeST-CA
36 207 325 666 859 8.6K
274K 6.52M 16.94M 107M 11B 4.52B
13.24% 8.11% 0.003% ∼8% <20% n.r.
Cities Road Road River River Road
SWORD reaches SWOT RiverSP HydroWeb ICESat-2
19,172 10K 3.8K 18K
– (static) 353K 391K 208K
– 2023–26 2016–26 2018–26
1.69% 0.50% 0.32%
Total observed In situ (eval)
18.5K 375
1.9M 955K
2016–26 2016–26
0.97% 1.32%
AmazonWSE
19.2K
1.9M
99%
River
Table 1: Spatiotemporal dataset comparison. Obs.: measurements, reported or estimated from time steps and sparsity.
thermore, under extreme sparsity, same-day spatial neighborhoods contain almost no observed nodes, either starving the message-passing mechanism or introducing prohibitive amounts of virtual tokens. Finally, we incorporate multiple sources that may observe the same node on the same day, yet due to different sensor characteristics, yield a different measurement. Figure 1 shows example time series and a fragment of the underlying topology. We show that recent transductive and inductive graph or recurrence based imputation methods do not perform well in our setting. We introduce a simple sequence model which flattens space and time into a single token sequence of observed measurements. This model handles the sparsity more naturally than grid-based STGNNs and its inductive bias fits the sequential upstream→downstream and temporal order inherent to river systems. Our results support recent findings that GNNs applied to river networks for discharge prediction show limited benefit from graph topology (Kirschstein and Sun 2024). Importantly, we find that river topology does inform better predictions, but is better exploited through subgraph sampling and metadata encodings than through explicit graph convolutions. In summary, this paper makes the following contributions:
Table 2: Dataset statistics for the Amazon basin. Percentage observed (% obs.) is calculated over the respective periods of operation of the sources.
gates measurements from ERS-1/2, Topex/Poseidon, Jason1/2/3, Envisat, Saral/AltiKa, Sentinel-3A/B, and Sentinel-6A into time series at locations where the satellite ground track crosses a river (Santos da Silva et al. 2010; Normandin et al. 2018). While other altimeters observe only along their precise ground track, the novel wide-swath InSAR sensor of the SWOT satellite observes entire contiguous river segments rather than point crossings (Biancamaria, Lettenmaier, and Pavelsky 2016). Reach-Reg (Halicki et al. 2026) exploits the SWOT measurement geometry for spatiotemporal WSE densification by chaining linear regressions between simultaneously observed reaches and by modeling water velocity. They achieve the best accuracy with re-processed SWOT data which is not available at scale (Schwatke et al. 2015) but show that their method also works on public SWOT RiverSP data. However, like other works (Tourian et al. 2016; Nielsen et al. 2022), they consider only large and well observed rivers without complex topology.
Related Work
Spatiotemporal Imputation. Early work focused on modeling the temporal dimension (Yi et al. 2016; Cao et al. 2018). Subsequent works model spatial interactions more explicitly and can be categorized into transductive methods, that assume to see every to-be-predicted node during training (Cini, Marisca, and Alippi 2022; Marisca, Cini, and Alippi 2022; Liu et al. 2023a; Cheng et al. 2024; De Felice et al. 2024; Nie et al. 2024; Yang et al. 2025b) and inductive methods, that transfer to entirely unseen nodes (Wu et al. 2021; Zheng et al. 2023; Li et al. 2025; Xu et al. 2025; Ren et al. 2026; Liang et al. 2026). Even though some methods introduce mitigations for the cost of processing large graphs, none of these has been evaluated on graphs with more than 1,200 nodes. SPIN (Marisca, Cini, and Alippi 2022) and IGNNK (Wu et al. 2021) support subgraph sampling, similar to the model proposed here. Only SPIN (Marisca, Cini, and Alippi 2022) and ImputeFormer (Nie et al. 2024) have been shown to work (transductively) at sparsity levels exceeding 90%, but not up to 99% (Marisca, Cini, and Alippi 2022; Nie et al. 2024). Forecasting methods exist for larger but denser graphs (Cini et al. 2023; Liu et al. 2023b).
Satellite Altimetry for River Monitoring. Classical altimetry satellites provide an estimate of the WSE along their ground-track by measuring the distance between the satellite and the water surface with a radar, and have provided observations since 1991. The HydroWeb database aggre-
Existing Datasets. Benchmarks used for spatiotemporal imputation, though having more observations along their temporal axis, are substantially smaller in their spatial extent than AmazonWSE (Table 1). AQI-36 contains 36 hourly airquality series, METR-LA and PEMS-BAY contain 207 and
1. We introduce a dataset of Amazon WSE time series on an underlying river graph that observes 19K river sections with a daily frequency over 2016–2026 and serves as a benchmark for imputation under extreme sparsity. 2. We propose a Mamba-based sequence model as baseline that we train for masked reconstruction on the heterogeneous observations with topology-aware subgraph samples (Gu and Dao 2023). 3. We show that this model outperforms interpolation and existing transductive and inductive neural methods, and provide ablation studies. We demonstrate improvements in coverage and accuracy against Reach-Reg (Halicki et al. 2026), a current state-of-the-art method for densification of altimetry-derived WSE from SWOT data.
Figure 2: Spatial distribution of data sources in AmazonWSE, after quality filters, with SWORD river topology background. Left to right: SWORD reaches, SWOT-observed reaches, HydroWeb virtual stations, ICESat-2 transects, ANA in situ gauges. 325 five-minute traffic series, respectively (Zheng, Liu, and Hsieh 2013; Li et al. 2018). LargeST scales traffic forecasting to 8,600 sensors over five years, and CER-E and PV-US contain energy grid time series over ∼5,000-6,000 nodes, but none of these has been used for imputation (Commission for Energy Regulation 2012; Hummon et al. 2012; Liu et al. 2023b). Table 1 shows that these datasets offer denser data, which hampers the comparability of imputation studies that report incompatible artificial sparsification scenarios, from random missing entries (Nie et al. 2024) and temporal/spatial blocks (Liu et al. 2023a) to entirely withheld target nodes (Yang et al. 2025a), removed context nodes (Zheng et al. 2023), and held-out variable channels (De Felice et al. 2024). AmazonWSE offers a consistent benchmark for extreme sparsity; by combining 1.9M irregular observations at 19K nodes over 2016–2026, fewer than 1% of reaches are observed on any day. Its directed acyclic topology is also novel with respect to the cyclical graphs of existing datasets. Hydrology datasets with discharge or water level time series aggregate data spatially into catchments and thereby erase the topology (Caravan; Kratzert et al. 2023), or only make available in situ gauge data, no satellite altimetry: LamaH-CE; (Klingler, Schulz, and Herrnegger 2021) and CausalRivers (Stein et al. 2025). Rivers and GNNs. Acosta et al. (2025); Taghizadeh et al. (2025) propose GNNs for spatial flood forecasting. Dufourg et al. (2024) find that a GNN slightly outperforms (Conv)LSTMs when forecasting a water index from satellite image time series. Kirschstein and Sun (2024), on the other hand, report that reflecting river network topology through node adjacency in GNNs does not improve discharge forecasting over using isolated node-level MLPs. Here, we find that topology information does improve WSE reconstruction if used to sample context and encoded through positional encodings rather than through GNN message passing.
Dataset We construct a multi-source dataset of water surface elevation time series for the Amazon river by integrating five data sources, each with different spatial and temporal characteristics. Figure 2 shows the spatial distribution of the data sources. Table 2 summarizes the resulting dataset.
SWORD and Topology The SWORD database (v17b) was designed as a complementary static dataset to SWOT products and provides a
topological graph of 19,172 river reaches in the Amazon basin (Altenau et al. 2021). Each reach is characterized by a unique ID, geographic coordinates, and upstream/downstream connectivity. SWORD provides the graph structure, shown in Figures 1-3, on which our model operates and the static metadata used for spatial encodings. The river topology defines a directed graph G = (V, E), where V is the set of SWORD reaches and (u, v) ∈ E indicates that v is immediately downstream of u. The graph is acyclic, v ̸⇝ v, and approximately but not strictly a tree, since reaches sometimes have multiple parents.
Observation Sources SWOT RiverSP. The SWOT RiverSP product provides WSE observations mapped to SWORD reaches, from the wide-swath InSAR sensor (SWOT 2025). We use data from the science cycle that started in July, 2023. Of the 19,172 SWORD reaches, approximately 17K are observed by SWOT, and 10K are retained by quality filtering. The repeat cycle takes 21 days, but due to the wide-swath sensor some reaches are observed several times during this period. HydroWeb.next. The HydroWeb database provides WSE time series at 7.3K river crossings in the Amazon basin, derived from altimetry missions (including S3A, shown in Figure 1) spanning 1993–2026 (Crétaux and Calmant 2015). We use both the operational and research collections and retain 3.7K locations with data between 2016–2026 matched to a nearby reach. Each virtual station has a measurement every 10–35 days depending on the altimetry mission. ANA In Situ Gauges. The Brazilian National Water Agency (ANA) operates flow gauges across the Amazon basin. We retrieve WSE records via their API (Agencia Nacional de Aguas e Saneamento Basico (ANA) 2026). We retain 375 gauges with data between 2016–2026 after quality filtering, manual inspection and matching to SWORD reaches. Following standard practice in hydrology, we provide the gauge data as ground truth. They provide dense temporal sampling (with an aggregated daily observation during periods of operation) but are geographically sparse. ICESat-2. The NASA Ice, Cloud and land Elevation Satellite-2 (ICESat-2) mission carries a laser altimetry instrument, which measures along three pairs of narrow laser beams (Abdalati et al. 2010). We use the mean along-track surface elevation for each beam for each transect across a water body from the Level3B ATL22 Mean Inland Surface
Processing and Quality Filtering HydroWeb, ANA, and ICESat-2 locations are assigned to the nearest SWORD reach, retaining only matches within 10 km (unmatched locations are discarded, imposing a stricter limit than the 20 km limit used by Halicki et al. (2026)). SWOT RiverSP observations are already mapped to SWORD reaches. HydroWeb, SWOT, and ICESat-2 report elevations referenced to the EGM2008 geoid (Pavlis et al. 2012), which provides a gravity-adjusted estimate of “elevation above sea level”. ANA records are in a local gauge datum; for evaluation, each gauge is therefore mapped to the prediction datum by a fitted linear transformation, following Halicki et al. (2026). Every source is placed on the common daily grid from 2016-01-01 to 2026-05-01 (3,774 days). Further processing is source-specific and according to expert hydrology standards. For SWOT, we screen product quality fields including reach quality, width, crosstrack distance, crossover calibration, random WSE uncertainty, and severe bit flags. These criteria combine filters adapted from Andreadis et al. (2025) and Halicki et al. (2026) with additional harmonic-residual, observed-pixel precision and ∆width/∆WSE consistency checks we introduce. HydroWeb text products are parsed from both operational and research collections; invalid fill values and the short Jason-2 interleaved (J2N) record are removed before daily alignment. ICESat-2 transects over reservoirs and transects shorter than 50 m are discarded, followed by a per-reach harmonic-residual filter. ANA sub-daily measurements are reduced to a daily median, and filtered using seasonal-trend residuals, followed by manual inspection by an expert. More details are given in Appendix A.2. All data (filtered and original measurements, uncertainties, quality flags, metadata) is provided in h5netcdf format and can easily be read by xarray. Each dataset may be used freely for any purpose under its provider’s terms; some providers require attribution. Appendix A.3 explains the storage format. All data export, processing, loading and model code will be made public upon acceptance.1
Task: Spatiotemporal Imputation Given whichever altimetry measurements are available, the task is to reconstruct a daily WSE at each SWORD reach. The evaluation is broken down in two tracks: reconstruction during the SWOT era (2023-07 to 2026-05), and a pre-SWOT hindcast. ANA gauges are not used as model inputs or as training targets but are reserved only for evaluation.
Model We formulate WSE densification as masked reconstruction on the directed SWORD graph G = (V, E). An observation (r) is a tuple (j, i, t, r, hj,t ), where source location j is mapped 1
https://github.com/rubencart/AmazonWSE
−0.5 −1.0 −1.5 Latitude
Water Data product (Jasinski et al. 2025). Observations span Oct. 2018–Mar. 2026; with a nominal repeat cycle of 91 days and a large spatial coverage. After filtering and matching, time series remain for >90% of SWORD reaches.
−2.0 −2.5 −3.0 −58
−57
−56
Longitude
−55
−54
−53
Figure 3: Example of sampled subgraph (local edges in grey) with river topology (SWORD reaches in blue). to SWORD reach i ∈ V by i = ν(j), t is a daily date bin, (r) r identifies the observing source, and hj,t is the measured WSE. Multiple sources may therefore produce distinct tokens for the same reach and day.
Sample Construction A sample is defined by an anchor location a, a contiguous time interval I, and a connected local neighborhood Va ⊂ V obtained by following the SWORD topology upstream and downstream from the anchor (Example in Figure 3). Let Ja be the selected source locations mapped to these reaches. The observed part of the sample is Da,I = (j, i, t, r) : j ∈ Ja , i = ν(j), t ∈ I, (1) (r)
hj,t is observed . We create tokens only for available measurements rather than constructing a dense reach–time grid. Thus, sparsity reduces sequence length instead of filling the input with missingvalue tokens. Token Representation. For token ℓ, let zℓ be the WSE normalized by the mean/std of the entire time series at source location j (different sources observing the same location are normalized separately). The embedded token eℓ that is input to the model is a linear projection of zℓ if the token is not masked and a learned mask token embedding emask if it is masked (for training) or a query token (for inference). Node-Independent Metadata Encodings. Each dynamic token (masked and non-masked) receives temporal, source, geographic, and topological metadata: pℓ = Emonth [mℓ ] + PEday (tℓ − tmin ) + Esrc [rℓ ] + Wa caiℓ + Wg cgiℓ + TreePE(biℓ ).
(2)
Here, mℓ is the month, Emonth and Esrc are month and source embedding tables, PEday is a sine/cosine position encoding, tℓ − tmin is the token’s date index offset to the earliest token in the sample, cai contains normalized coordinates relative to the most-downstream sampled reach, cgi contains absolute geographic coordinates, and bi is the local branch path from reach i to the sampled downstream root. The source embedding distinguishes measurements produced by different altimeters. Metadata are added to model inputs before every layer: eℓ ← eℓ + pℓ .
...
Dim = 1
Model input
Linear Projection
Mamba Backward Scan
xL
Add query tokens
... Model input Forward Scan
Dim = 192 Linear Projection + Metadata Encoding Add query tokens
Dim = 1
Masking
Figure 4: Model architecture . During training, sparse WSE observations (colored) are sampled and part of these are masked. Learned query tokens (grey) and temporal, satellite, coordinate, and river-tree metadata are added. Tokens are processed by N bidirectional selective-SSM layers, and mapped to normalized WSE by satellite-specific heads. In inference, the masking is replaced by the introduction of query tokens carrying the metadata of what is to be predicted. Each sampled location additionally contributes a static token encoding its mean WSE w.r.t. the most-downstream reach in the sample: vimean = Normβ (µi − µi∗ ), where i∗ is the most-downstream sampled reach. Static metadata tokens form a prefix. Dynamic tokens are then sorted so that they are ordered first by day and then along the river topology according to downstream flow rank. We approximate a local subgraph with a tree and encode the position of reach i by the sequence of branch choices bi = (bi,0 , . . . , bi,K−1 ) that separate i from the local downstream root following Shiv and Quirk (2019). Let ui,k ∈ {0, 1}B be the one-hot encoding of branch choice bi,k , where B = 3 is the branching factor and ui,k = 0 for padded path positions. A decay rate applied to depth k is computed from learned weights w ∈ RF , for channel f = 1, . . . , F : r F k 1 − ρ2f , ρf = tanh(wf ) (3) [gk ]f = ρf 2 The tree encoding is obtained by assigning this weight vector to the branch selected at each depth: ⊤ TreePE(bi ) = Wtree vec concatK−1 (4) k=0 ui,k gk . Each path depth therefore activates one branch-specific block of F components. The learned geometric decay rates allow different channels to emphasize different parts of the path, while shared path segments receive the same encoding. Learned weights Wtree project the encoding from K · B · F dimensions to the model dimension. These encodings contain no learned lookup indexed by node index (SWORD reach). The same functions encode coordinates, relative branch paths, source identities, and scalar
Figure 5: Construction of inference sample that is input to our model (top) vs. to graph based methods for imputation like SPIN (Marisca, Cini, and Alippi 2022) or KITS (Xu et al. 2025). Learned query tokens are added to sparse conditioning inputs to inform the model of what to decode. Graph based methods add large amounts of query tokens (grey) to obtain a full spatiotemporal grid, while we add tokens only for a single location, but still using context from different locations. WSE metadata at every location. The model is therefore inductive to reaches not observed during training, provided their SWORD topology and static metadata are available.
Bidirectional Mamba The ordered sequence is processed using standard Mamba blocks (Gu and Dao 2023). At its core, Mamba applies an input-dependent state-space recurrence sℓ = Aℓ sℓ−1 + Bℓ xℓ , yℓ = Cℓ sℓ + Dxℓ ,
(5) (6)
where the discretization and the maps Bℓ and Cℓ depend on the current input xℓ . This selectivity allows the model to retain or discard information as it scans the sequence. Each residual block has the form MBlock(X) = X + Mamba(LN(X)) .
(7)
The standard formulation of Mamba is 1-directional, but we combine a forward and backward scan, as proposed by Zhu et al. (2024), but in our case to integrate past information from upstream nodes with future information from downstream nodes. After restoring the reverse output to the original order, both directions are fused: h i (k) H(k) = Wbi H(k) . (8) → ∥ revdyn H← Only the dynamic sequence is reversed; the static metadata prefix remains fixed. The resulting model uses observations on both sides of a query and is therefore a non-causal reconstruction model rather than a forecaster. A source-specific linear head maps each final token representation to normalized WSE prediction zbℓ .
RMSE ↓ >2023/07 <2022/06
Inductive kNN IGNNK KITS Ours
5 25 012
20 -1
Avg nearby observations / month (binned)
Figure 6: Ours vs. baseline RMSE as a function of the number of observations in the subgraph neighborhood per month.
Table 3: Baseline comparison on held-out in situ gauges for the 2.7% (>2023/07) and 0.6% (<2022/06) sparsity regimes. sub: subgraph sampling and full: full graph in each sample.
Training Training combines location masking, which hides complete source-location time series to teach spatial reconstruction at unseen locations, and random masking, which hides individual observations to teach temporal densification. Static metadata remain visible. If Ω is the set of masked, non-padding observations in a minibatch, we minimize 1 X 2 L(θ) = (b z ℓ − zℓ ) . (9) |Ω| ℓ∈Ω
Inference To predict reach i over interval I, we use it as anchor a = i to construct a subgraph sample as described above, and we add i as virtual location containing exactly one masked query token (i, t) for each day t ∈ I. No daily queries are introduced for the other reaches in the sampled neighborhood: their sparse satellite measurements remain the observed context, inf Di,I = Di,I ∪ {(i, t, masked) : t ∈ I} .
83
0.84
3
0.62
-8
sub
0.4
49
2.24 2.38 1.34 2.40 1.61
0.6
9
1.77 2.33 1.11 2.38 1.23
0.8
-4
full full sub full sub
RMSE
2.47 1.28 2.25 1.58 0.94 1.56 0.91
30
ImputeFormer
1.47 0.88 2.52 0.93 0.74 1.69 0.67
0
SPIN-H
sub full sub full sub full sub
-3
Transductive Temp. LSTM GRIN
Ours ImputeFormer SPIN-H
1.0
14
Model
(10)
Predicting only the target reach prevents the sequence from being dominated by virtual tokens and maintains a favorable ratio of observed context to queries even under extreme graph sparsity (illustrated in Figure 5). Each window is reconstructed independently, without feeding predictions back into the model, and we average overlapping windows.
Experiments Hyperparameters are listed in Appendix A.4. Our model has 3.3M parameters and takes 2 h to train with less than 4 GB of GPU memory on one H100, predicting all 19K nodes takes another 1,5 h. Training data spans 2023-07–2026-05 and combines SWOT with HydroWeb and ICESat-2 observations. For our training and for baselines, time series with less
than 10 (for ICESat-2) or 20 (for the remaining sources) values are discarded, reducing the number of observed ICESat-2 locations from 18.5K to 10K (the data are kept in the published dataset for subsequent works to use as they see fit). With ∼75% of reaches observed by the remaining time series, our task involves both transduction and induction, even though 95% of the evaluation gauges are located on a reach that is observed (this can be expected to slightly disadvantage transductive methods). We evaluate in two temporal settings: 1. SWOT period (2023-07 to 2026-05): The sparsity of this period is lower due to the presence of SWOT measurements (in both training and inference): 2.7% observed. 2. Hindcast (2016-01 to 2022-06-30): the model must extrapolate backward in time using only classical altimetry from HydroWeb and ICESat-2 observing 0.6% of day– reach slots, without any SWOT data. We keep gauge data from the year 2022-07 to 2023-07 as evaluation setting for ablations (not for early stopping). The model never sees in situ data during training, preventing leakage also when the training and inference periods coincide. Baselines. We select baselines to cover complementary settings and model families. A per-reach bidirectional LSTM inspired by BRITS-I (Cao et al. 2018) and spatial kNN isolate temporal modeling and non-parametric spatiotemporal interpolation. Among transductive imputers, GRIN (Cini, Marisca, and Alippi 2022) represents recurrent graph message passing, while SPIN-H (Marisca, Cini, and Alippi 2022) and ImputeFormer (Nie et al. 2024) represent attention-based approaches and demonstrated high (95%) point sparsity. IGNNK (Wu et al. 2021) and KITS (Xu et al. 2025) test inductive reconstruction at unseen reaches; KITS is especially relevant because it addresses the gap between sparse training graphs and virtual targets at inference. Some baselines natively support multivariate inputs; for those that do not, we add separate source channels, compare them with merging observations into one WSE series, and report whichever performs best. We run graph baselines both on the full graph and on the same sampled subgraphs as our model. We also compare with Reach-Reg (Halicki et al. 2026), the domainspecific baseline for SWOT WSE densification, which we fit on SWOT RiverSP data to make predictions for our gauges. Appendix A.5 describes baseline details.
RMSE ↓
Configuration Base
RMSE ↓
0.56
Period
Cov. RR
KGE ↑
IF
RR
IF
Metadata enc
No tree encoding No satellite enc No mean WSE tokens
0.58 0.58 0.58
SWOT era 184/284 0.90 0.61 0.55 0.85 0.92 0.94 Hindcast 149/250 0.85 0.80 0.70 0.88 0.88 0.91
Data sources
No ICESat-2 No SWOT & ICEsat-2
0.59 0.60
Token order
Flow Random
0.57 0.64
Architecture & inputs
1-directional Isolated
0.72 2.58
Table 5: Comparison with Reach-Reg (RR) and ImputeFormer (IF). Mean metrics are computed on gauges overlapping with Reach-Reg; Cov. reports RR/total gauge coverage. Scores differ from other tables because we evaluate only on gauges for which Reach-Reg makes a prediction.
7 14 28
1.0 RMSE
Table 4: Ablation study (validation period, mean over heldout gauges). Base has all metadata encodings and tokens, uses all 3 data sources in training and inference (when available), sees input sequences sorted by time, and is 2-directional.
Validation RMSE
1.1 0.9
Days
91 182 365
0.8 0.7
Results Baseline Comparison. Table 3 compares our model with existing ML methods; we outperform all methods in all settings. It also shows that baselines have a larger gap between performance on SWOT period imputation where 2.7% of days have an input measurement and on the Hindcast setting where only 0.6% is observed. SPIN-H and ImputeFormer with subgraphs reach good accuracy, which is in line with their reported improved results on sparse settings compared to other models. Except for GRIN, subgraph sampling outperforms handling the full graph at once, we attribute this to our larger graph compared to those of existing datasets. The transductive baselines do better than inductive baselines (but not than our model which is also inductive). Figure 6 shows performance of our model vs. ImputeFormer and SPIN-H, in function of number of the average amount of context tokens available in the subgraph (colored tokens in Figure 5), where baseline performance deteriorates faster when less conditioning observations are available. These results support the hypothesis that sequence models outperform GNN-based approaches under extreme spatiotemporal sparsity. Ablations. Table 4 presents ablation experiments on metadata encoding and data source inclusion. In contrast to Kirschstein and Sun (2024), we find that graph topology and spatial context do help: the metadata encodings as well as tokens being ordered by time or flow improve predictions slightly and the subgraph model (Base) significantly outperforms the Isolated model (that sees only 1 node timeseries). All data sources (SWOT, ICESat-2 and HydroWeb) contribute to better performance. Figure 7 shows a comparison of time window size and the spatial neighborhood size that the subgraph is sampled from: a neighborhood size of 3001000 km works best and performance is robust to temporal window size, with 3 months giving the lowest RMSE. Comparison with Reach-Reg. Table 5 compares our model against Reach-Reg and ImputeFormer on overlapping in situ gauges across both evaluation periods. Kling– Gupta Efficiency (Gupta et al. 2009) is defined as KGE =
0.6 0
1000
2000
Km
3000
4000
5000
Figure 7: Ablation of sampled subgraph size limits. p 1 − (r − 1)2 + (α − 1)2 + (β − 1)2 , where r is the Pearson correlation between predictions and observations, α = σpred /σobs measures variability, and β = µpred /µobs measures bias. It is widely used in Hydrology and jointly evaluates correlation, variability, and mean agreement, with a perfect score of 1 and higher values indicating better performance. Our model outperforms both methods on both periods and provides greater coverage than Reach-Reg. We evaluate Reach-Reg for making hindcasts, even though it was originally developed only to make predictions for the SWOT-era (details in Appendix A.5). As a state-of-the-art hydrology method, it sets a level of accuracy along with a coverage baseline that new methods should aim to meet or surpass on the AmazonWSE benchmark. Visualizations and comparisons of predictions are provided in Appendix A.6.
Conclusion We introduced AmazonWSE, a benchmark dataset for spatiotemporal graph imputation on time series of water levels observed satellite altimetry, covering the Amazon river basin, with 19K river reaches as nodes, 99% sparsity, and a directed acyclic graph topology. We showed that prior spatiotemporal graph imputation methods are not adapted to this scale and sparsity, and proposed a bidirectional selective state space model that outperforms them by sampling connected subgraphs and flattening space and time into a single token sequence. Our model also outperforms a state-of-theart non-neural baseline for SWOT-based WSE densification and produces denser reconstructions. We hope that these results inspire the community to develop improved methods for spatiotemporal imputation in extremely sparse settings such as river monitoring from satellite altimetry.
References Abdalati, W.; Zwally, H. J.; Bindschadler, R.; Csatho, B.; Farrell, S. L.; Fricker, H. A.; Harding, D.; Kwok, R.; Lefsky, M.; Markus, T.; Marshak, A.; Neumann, T.; Palm, S.; Schutz, B.; Smith, B.; Spinhirne, J.; and Webb, C. 2010. The ICESat2 Laser Altimetry Mission. Proceedings of the IEEE, 98(5): 735–751. Acosta, C. M.; Herath, H. M. V. V.; Lim, J. Y.; Saha, A.; Rasnayaka, S.; and Marshall, L. 2025. DUALFloodGNN: Physics-informed Graph Neural Network for Operational Flood Modeling. arXiv:2512.23964. Agencia Nacional de Aguas e Saneamento Basico (ANA). 2026. HidroWeb Web Service API. https://www.ana.gov.br/ hidrowebservice/swagger-ui. Accessed: February 5, 2026. Altenau, E. H.; Pavelsky, T. M.; Durand, M. T.; Yang, X.; Frasson, R. P. D. M.; and Bendezu, L. 2021. The Surface Water and Ocean Topography (SWOT) Mission River Database (SWORD): A Global River Network for Satellite Data Products. Water Resources Research, 57(7): e2021WR030054. Andreadis, K. M.; Coss, S. P.; Durand, M.; Gleason, C. J.; Simmons, T. T.; Tebaldi, N.; et al. 2025. A First Look at River Discharge Estimation from SWOT Satellite Observations. Geophysical Research Letters, 52(9): e2024GL114185. Biancamaria, S.; Lettenmaier, D. P.; and Pavelsky, T. M. 2016. The SWOT Mission and Its Capabilities for Land Hydrology. Surveys in Geophysics, 37: 307–337. Cao, W.; Wang, D.; Li, J.; Zhou, H.; Li, L.; and Li, Y. 2018. BRITS: Bidirectional Recurrent Imputation for Time Series. In Advances in Neural Information Processing Systems, volume 31, 6775–6785. Cheng, S.; Osman, N.; Qu, S.; and Ballan, L. 2024. FastSTI: A fast conditional pseudo numerical diffusion model for spatio-temporal traffic data imputation. IEEE Transactions on Intelligent Transportation Systems, 25(12): 20547–20560. Cini, A.; Marisca, I.; and Alippi, C. 2022. Filling the G_ap_s: Multivariate Time Series Imputation by Graph Neural Networks. In International Conference on Learning Representations. Cini, A.; Marisca, I.; Bianchi, F. M.; and Alippi, C. 2023. Scalable spatiotemporal graph neural networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, 7218–7226. Commission for Energy Regulation. 2012. CER smart metering project-electricity customer behaviour trial, 2009-2010 [Dataset]. https://github.com/wwzjustin/CER-Smart-MeterProject-by-Irish-Social-Science-Data-Archive. Accessed: February 5, 2026. Crétaux, J.-F.; and Calmant, S. 2015. Hauteur des lacs et des rivières (HYDROWEB) [Dataset]. https://doi.org/ 10.24400/329360/HYDROWEB\_WATER\_LEVEL. Accessed: February 5, 2026. De Felice, G.; Cini, A.; Zambon, D.; Gusev, V.; and Alippi, C. 2024. Graph-based virtual sensing from sparse and partial multivariate observations. In International Conference on Learning Representations, 17111–17132.
Dufourg, C.; Pelletier, C.; May, S.; and Lefèvre, S. 2024. Forecasting water resources from satellite image time series using a graph-based learning strategy. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 48: 81–88. Fassoni-Andrade, A. C.; Fleischmann, A. S.; Papa, F.; Paiva, R. C. D.; Wongchuig, S.; Melack, J. M.; et al. 2021. Amazon Hydrology from Space: Scientific Advances and Future Challenges. Reviews of Geophysics, 59(4): e2020RG000728. Gu, A.; and Dao, T. 2023. Mamba: Linear-Time Sequence Modeling with Selective State Spaces. arXiv:2312.00752. Gupta, H. V.; Kling, H.; Yilmaz, K. K.; and Martinez, G. F. 2009. Decomposition of the Mean Squared Error and NSE Performance Criteria: Implications for Improving Hydrological Modelling. Journal of Hydrology, 377(1–2): 80–91. Halicki, M.; Niedzielski, T.; Schwatke, C.; Scherer, D.; and Dettmering, D. 2026. Daily River Water Levels from MultiMission Altimetry: A Reach-Based Regression Method Using the Unique SWOT Data Geometry. Journal of Hydrology, 673: 135367. Hummon, M.; Ibanez, E.; Brinkman, G.; and Lew, D. 2012. Sub-hour solar data for power system modeling from static spatial variability analysis. Technical report, National Renewable Energy Laboratory (NREL), Golden, CO (United States). Jasinski, M.; Stoll, J.; Hancock, D.; Robbins, J.; and Nattala, J. 2025. ATLAS/ICESat-2 L3B Mean Inland Surface Water Data, Version 4 [Dataset]. https://doi.org/10.5067/ATLAS/ ATL22.004. Accessed: February 5, 2026. Kirschstein, N.; and Sun, Y. 2024. The Merit of River Network Topology for Neural Flood Forecasting. In Proceedings of the 41st International Conference on Machine Learning, 24713–24725. Klingler, C.; Schulz, K.; and Herrnegger, M. 2021. LamaHCE: LArge-SaMple DAta for hydrology and environmental sciences for central Europe. Earth System Science Data, 13(9): 4529–4565. Kratzert, F.; Nearing, G.; Addor, N.; Erickson, T.; Gauch, M.; Gilon, O.; Gudmundsson, L.; Hassidim, A.; Klotz, D.; Nevo, S.; et al. 2023. Caravan-A global community dataset for large-sample hydrology. Scientific Data, 10(1): 61. Li, Y.; Yu, R.; Shahabi, C.; and Liu, Y. 2018. Diffusion Convolutional Recurrent Neural Network: Data-Driven Traffic Forecasting. In International Conference on Learning Representations. Li, Y.; Zezhi, S.; Yu, C.; Qian, T.; Zhang, Z.; Du, Y.; He, S.; Wang, F.; and Xu, Y. 2025. Sta-gann: A valid and generalizable spatio-temporal kriging approach. In Proceedings of the 34th ACM International Conference on Information and Knowledge Management, 1726–1736. Liang, Z.; Li, W.; Zhang, D.; Jia, Z.; Chen, Y.; Wang, Z.; Zheng, X.; and Youssef, M. 2026. DarkFarseer: Robust Spatio-Temporal Kriging Under Graph Sparsity and Noise. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 40, 23451–23459.
Liu, M.; Huang, H.; Feng, H.; Sun, L.; Du, B.; and Fu, Y. 2023a. Pristi: A conditional diffusion framework for spatiotemporal imputation. In Proceedings of the IEEE International Conference on Data Engineering (ICDE), 1927–1939. Liu, X.; Xia, Y.; Liang, Y.; Hu, J.; Wang, Y.; Bai, L.; Huang, C.; Liu, Z.; Hooi, B.; and Zimmermann, R. 2023b. LargeST: A Benchmark Dataset for Large-Scale Traffic Forecasting. In Advances in Neural Information Processing Systems, volume 36, 75354–75371. Marisca, I.; Cini, A.; and Alippi, C. 2022. Learning to Reconstruct Missing Data from Spatiotemporal Graphs with Sparse Observations. In Advances in Neural Information Processing Systems, volume 35. Nie, T.; Qin, G.; Ma, W.; Mei, Y.; and Sun, J. 2024. ImputeFormer: Low Rankness-Induced Transformers for Generalizable Spatiotemporal Imputation. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, 2260–2271. Nielsen, K.; Zakharova, E.; Tarpanelli, A.; Andersen, O. B.; and Benveniste, J. 2022. River levels from multi mission altimetry, a statistical approach. Remote Sensing of Environment, 270: 112876. Normandin, C.; Frappart, F.; Diepkilé, A. T.; Marieu, V.; Mougin, E.; Blarel, F.; et al. 2018. Evolution of the Performances of Radar Altimetry Missions from ERS-2 to Sentinel-3A over the Inner Niger Delta. Remote Sensing, 10(6): 833. Pavlis, N. K.; Holmes, S. A.; Kenyon, S. C.; and Factor, J. K. 2012. The Development and Evaluation of the Earth Gravitational Model 2008 (EGM2008). Journal of Geophysical Research: Solid Earth, 117(B4): B04406. Ren, X.; Zhao, K.; Taškova, K.; and Riddle, P. 2026. AnchorGK: Anchor-based Incremental and Stratified Graph Learning Framework for Inductive Spatio-Temporal Kriging. In Proceedings of the ACM SIGKDD Conference on Knowledge Discovery and Data Mining, 1251–1262. Santos da Silva, J.; Calmant, S.; Seyler, F.; Rotunno Filho, O. C.; Cochonneau, G.; and Mansur, W. J. 2010. Water Levels in the Amazon Basin Derived from the ERS-2 and ENVISAT Radar Altimetry Missions. Remote Sensing of Environment, 114(10): 2160–2181. Schwatke, C.; Dettmering, D.; Bosch, W.; and Seitz, F. 2015. DAHITI–an innovative approach for estimating water level time series over inland waters using multi-mission satellite altimetry. Hydrology and Earth System Sciences, 19(10): 4345–4364. Shiv, V.; and Quirk, C. 2019. Novel Positional Encodings to Enable Tree-Based Transformers. In Advances in Neural Information Processing Systems, volume 32. Stein, G.; Shadaydeh, M.; Blunk, J.; Penzel, N.; and Denzler, J. 2025. CausalRivers – Scaling Up Benchmarking of Causal Discovery for Real-World Time-Series. In International Conference on Learning Representations. SWOT. 2025. SWOT Level 2 River Single-Pass Vector Data Product [Dataset]. https://doi.org/10.5067/SWOTRIVERSP-D. Accessed: February 5, 2026.
Taghizadeh, M.; Zandsalimi, Z.; Nabian, M. A.; ShafieeJood, M.; and Alemazkoor, N. 2025. Interpretable physics-informed graph neural networks for flood forecasting. Computer-Aided Civil and Infrastructure Engineering, 40(18): 2629–2649. Tourian, M.; Tarpanelli, A.; Elmi, O.; Qin, T.; Brocca, L.; Moramarco, T.; and Sneeuw, N. 2016. Spatiotemporal densification of river water level time series by multimission satellite altimetry. Water Resources Research, 52(2): 1140– 1159. Wu, Y.; Zhuang, D.; Labbe, A.; and Sun, L. 2021. Inductive graph neural networks for spatiotemporal kriging. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, 4478–4485. Xu, Q.; Long, C.; Li, Z.; Ruan, S.; Zhao, R.; and Li, Z. 2025. Kits: Inductive spatio-temporal kriging with increment training strategy. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, 12945–12953. Yang, C.; Zhao, C.; Wang, C.; and Fan, J. 2025a. DRIK: Distribution-Robust Inductive Kriging without Information Leakage. arXiv:2509.23631. Yang, X.; Sun, Y.; Chen, X.; Zhang, Y.; and Yuan, X. 2025b. Graph structure learning for spatial-temporal imputation: Adapting to node and feature scales. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, 959– 967. Yi, X.; Zheng, Y.; Zhang, J.; and Li, T. 2016. ST-MVL: Filling Missing Values in Geo-Sensory Time Series Data. In Proceedings of the International Joint Conference on Artificial Intelligence, 2704–2710. Zheng, C.; Fan, X.; Wang, C.; Qi, J.; Chen, C.; and Chen, L. 2023. Increase: Inductive graph representation learning for spatio-temporal kriging. In Proceedings of the ACM Web Conference, 673–683. Zheng, Y.; Liu, F.; and Hsieh, H.-P. 2013. U-air: When urban air quality inference meets big data. In Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 1436–1444. Zhu, L.; Liao, B.; Zhang, Q.; Wang, X.; Liu, W.; and Wang, X. 2024. Vision Mamba: Efficient Visual Representation Learning with Bidirectional State Space Model. In International Conference on Machine Learning, 62429–62442. PMLR.
A
Appendix
Overview of Supplementary Material • This Appendix (Technical Supplement): 1. Dataset Sources 2. Data Preprocessing Details 3. Storage Format 4. Hyperparameters and Model Details 5. Baselines 6. Predictions • Media Supplement: animated .gif files of 1) sparse ICESat-2 and HydroWeb inputs and 2) imputation prediction of daily WSE for all SWORD reaches by our model. • Code and Data Supplement: – Netcdf files for each source in the AmazonWSE dataset with a tiny amount of measurements: as an example of the storage format. – Code of data export and preprocessing scripts. – Code to run training and evaluation of our model and the baselines on the tiny included dataset.
A.1
Dataset Sources
Figure 8 shows the temporal coverage of different sources in AmazonWSE and in HydroWeb. The next paragraphs provide more details about the sources when relevant. SWORD. The Surface Water and Ocean Topography (SWOT) Mission River Database (v17b) was designed as a complementary static dataset to SWOT products and provides a topological graph of 19,172 river reaches in the Amazon basin (Altenau et al. 2021).2 ICESat-2. The NASA Ice, Cloud and land Elevation Satellite-2 (ICESat-2) mission carries the ATLAS (Advance Topographic Laser Altimeter System) instrument. Data is collected along three pairs of narrow laser beams that operate at 532nm wavelength. The beams in the pair are separated by a distance of around 90m, while each pair is separated by the next by a distance of 3 km. The data in AmazonWSE are from the Level3B ATL22 Mean Inland Surface Water Data product version 4 (Jasinski et al. 2025), which reports mean along-track surface elevation for each beam for each transect across a water body. Our observations span Oct. 2018–Mar. 2026 because the product has been released up until March 2026 at the time of writing. SWOT RiverSP. We use data from the “SWOT Level 2 River Single-Pass Vector Reach Data Product, Version D”, and of the science cycle that started in July, 2023 (SWOT 2025). HydroWeb.next. We use the operational and research collections (Crétaux and Calmant 2015) and include data from the following satellites: the SWOT nadir altimeter (i.e., a different sensor on the same satellite than the InSAR swath altimetry sensor used for the SWOT data source in AmazonWSE), Sentinel-3A/B, Sentinel-6A, JASON-2/3, SARAL. 2
https://zenodo.org/records/15299138.
A.2
Data Preprocessing Details
This section describes the data processing in more detail. All timestamps are first converted from local time to UTC (if not already in UTC). Subsequently, all measurements are assigned to daily bins on the shared temporal grid. A quality flag in the shared data files records missing, rejected (by quality filters), and retained observations. SWOT. Our product-level screening draws on the criteria used by Andreadis et al. (2025) and Halicki et al. (2026), with additional criteria added by ourselves. RiverSP records are aggregated into 24-hour bins. WSE and width use the median within each reach–day, provided measurement standard deviations are aggregated accordingly, and quality bit fields use the maximum so that severe flags are preserved. Many filters use variables shared as part of the SWOT RiverSP products, the names of those variables are marked here by the monospaced font, like wse_r_u. Altogether filters take away about 50% of measurements. This may seem like a lot, but examples shown in Figure 11 of unfiltered RiverSP SWOT time series along with collocated HydroWeb (“classical” altimetry) and in situ gauge time series can hopefully motivate the need for this extensive filtering. Filters that Apply to Reaches. • River Width. We remove reaches with a static (mean) width in SWORD lower than 50 m. Andreadis et al. (2025) used an 80 m cutoff for their early discharge validation; we instead use SWOT’s 50 m river-observation goal (Biancamaria, Lettenmaier, and Pavelsky 2016; Andreadis et al. 2025) to retain the mission’s intended lower-width range. • Reach Identifier. SWORD assigns identifiers to each reach of the form “CBBBBBRRRRT”, where the ending “T” encodes a reach type. The entire identifier follows the Pfafstetter coding system where the characters encode the topology of the river network (Verdin and Verdin 1999). We remove reaches with T ∈ {4, 5, 6}, since they encode “dam or waterfall”, “unreliable topology”, or “ghost reach/node”. Filters that Apply to Individual Measurements. • Cross-track Distance. Observations within 15 km crosstrack distance (x_trk_dist) of nadir ground track are removed, following Andreadis et al. (2025). • Crossover Calibration and WSE Uncertainty. Observations with crossover-calibration quality flag xovr_cal_q > 1 or with the random component of their reported WSE uncertainty wse_r_u > 0.5 m are removed, following Andreadis et al. (2025). The wse_u field provided in the RiverSP products (which estimates the standard deviation of the WSE measurement error of the SWOT sensor) is stored in our dataset’s netcdf files (after root-sum-square aggregation if a reach is observed more than once on a single day). • Reach Quality Flags. Observations with reach quality flag reach_q above 2 are removed, following Halicki et al. (2026). Unlike Halicki et al. (2026), who remove
Temporal Coverage of Hydroweb Data Sources SWT
Temporal Coverage of Amazon Data Sources ICESat-2
S6A
2026-03-10
2018-10-14
Source
2026-05-01
In-Situ
Source
2026-05-01
2023-03-29
2026-05-01
2016
2017
2018
2019
2020
2021
Date
2022
2023
2024
2025
2026
2026-05-01
2018-12-17
S3A
2026-05-01
2016-07-26
J3
2022-03-29
2016-02-26 2016-06-17 2016-01-01
SRL HydroWeb
2026-05-01
2021-09-21
S3B SWOT
2026-05-01
2023-07-21
J2
2016-09-25 2016-01-01
2016
2017
2018
2019
2020
2021
Date
2022
2023
2024
2025
2026
Figure 8: Left: temporal coverage of sources in AmazonWSE. Right: temporal coverage of satellite sensors included in HydroWeb: the SWOT nadir altimeter (i.e., a different sensor than the swath altimetry sensor used for the SWOT data source in AmazonWSE), Sentinel-3A/B, Sentinel-6A, JASON-2/3, SARAL.
Figure 9: Examples of SWOT observations retained and rejected by the harmonic-residual (left). Observations retained and discarded by the rating curve filter that considers ∆width/∆WSE (right). bitwise quality flags reach_q_b > 2097152, we remove only reach_q_b > 8388608 which additionally allows observations flagged as coming from a lake, but still removes severe flags like outliers.
tion is retained only if: ui,t 2q ≤p obs Ni,t
ui,t p p = 0.1 m.
obs Ni,t ≥4
• Observed Precision. Instead of applying a fixed threshold on the fraction of dark pixels (pixels for which no or not enough backscattered signal was measured and for which hence no measurement could be derived that is fed into the RiverSP measurement calculations) like (Andreadis et al. 2025; Halicki et al. 2026), we estimate whether the observed water pixels support a target precision p = 0.1 m. We approximate observed-pixel errors to be unbiased, independent, and approximately equal-variance, and norobs mally distributed. If Ni,t is the number of observed water pixels and ui,t is reported WSE uncertainty, an observa-
2 ,
Thus a larger reported uncertainty requires support from more observed pixels, where the factor of 2 approximates a two-sided 95% normal confidence interval. We crudely obs estimate Ni,t using the resolution of the InSAR sensor (∼7.5 m × 20 m and the static mean river width and reach length reported in SWORD (thereby assuming that the entire reach falls within one of the two 50 km wide observed swaths by the sensor). • Deviation from Harmonics. For reaches with sufficient observations, a two-harmonic seasonal model is fitted and WSE residuals beyond 4.5 standard deviations from the mean (residual) are removed (Figure 9, left). • Rating Curve. A robust delta rating-curve filter addi-
Figure 10: Interpolated SWOT mean and std fields that are used to de-normalize predictions. tionally identifies inconsistent changes in river width and WSE as observed by SWOT (Figure 9, right; river width and the water level are bound by a monotonously increasing relationship determined physically by the shape of the river bed). We fit a regression line on pairs of ∆width, ∆WSE between subsequent observations of the same reach, then remove observations of which the ∆width, ∆WSE deviates more than 4.5 standard deviations from the regressed line, or observations of which the ∆width and ∆WSE have non-matching signs (allowing a noise buffer). • Pass Correction. We investigate the bias between overpasses with different IDs that observe the same reach and found that the bias between passes is minimal. We try to correct it by shifting the mean WSE observed by the same pass ID towards each other, but found that that it makes negligible difference. • Mean/Std Field Interpolation. Because SWOT provides the densest coverage, we compute its per-location mean and standard deviation, and interpolate this linearly and spatially to obtain dense fields as shown in Figure 10. The interpolation uses inverse distance weighting and averages upstream branch contributions if there are several. We use these fields for de-normalization of predictions (since predictions may have to be made for locations where no measurements are available). HydroWeb. The HydroWeb database already cleans and homogenizes measurements as described in (Santos da Silva et al. 2010; Normandin et al. 2018). • Text files returned from the HydroWeb.next API3 (for the operational and research collections, or HYDROWEB_RIVERS_OPE and 3
https://github.com/CNES/py-hydroweb
HYDROWEB_RIVERS_RESEARCH) with time series are parsed, heights encoded as missing or fill values are converted to consistent missing data indicators; otherwise, HydroWeb’s supplied measurements and uncertainties are preserved without an additional value-level outlier filter. • Jason-2 Interleaved (J2N) is excluded because its Oct. 2016–May 2017 record is too short for stable statistics. SARAL and Jason-2 are retained, with normalization statistics estimated from 2012 onward rather than 2016 to use their longer records. ICESat-2. The filters remove about 4% of measurements (of the subset of measurements that has been matched to a SWORD reach). • Transect Matching. The ATL22 contains measurements per transect, which is a contiguous portion of one ICESat2 beam crossing a water body; land interruptions such as islands can split one beam pass into several transects. Transects are joined to their nearest SWORD reach. Reservoir observations and transects shorter than 50 m are excluded. • Deviation from Harmonics. For reaches with at least 10 transects, a two-harmonic model is fitted to transect orthometric heights and observations with residual modified Z-scores above 4.5 are rejected. • Transect Aggregation. Remaining transects are first aggregated within each reach-day and beam pass (source file and beam ID) using their median orthometric height. Reach-day WSE is the median of these beam-pass values, so a beam split into several transects does not receive extra weight. • Uncertainties. We compute wse_u as the leave-onebeam-pass-out jackknife standard uncertainty of the
Figure 11: Unfiltered SWOT RiverSP measurements (in red) with time series of collocated in situ gauges (top) and HydroWeb locations (bottom). While for some reaches SWOT measurements follow HydroWeb/in situ gauges closely (left), for other reaches, SWOT measurements are very noisy or low quality (right).
Figure 12: ICESat-2 harmonics filter (left) and aggregation of transect measurements per reach per day with accepted and rejected measurements (right). reach-day median. The ATL22 ht_stdev field is provided as the median within-transect standard deviation of filtered short-segment heights, but it should not be interpreted as a standard error of the transect mean. The output also retains the slope-fit RMSE, transect counts, beam-pass counts, and distinct-beam counts. ANA in situ gauges Figure 13 show examples or in situ gauges before and after quality filters. 21% of the 1000 timeseries are rejected automatically (10% of all measurements), of the remaining 810, 589 have a SWORD reach match within 10 km, and of those, we retain 375 series. • ANA Timestamps. are converted from Brazilian standard time (UTC-3), to UTC and binned to calendar days. Multiple observations within a day are represented by their median and converted from centimetres to meters. • Uncertainty. As ANA provides no per-observation measurement uncertainty, we use the within-day population standard deviation of the contributing readings as a dataderived uncertainty measure. • Minimum Amount and Quality. Stations must span at least two annual cycles (730 days) and contain at least 100 finite daily observations before filtering. Repeated-
value runs longer than three observations are removed. Negative observations are retained because gauge stage is referenced to a local datum. • Seasonality. Daily series are linearly interpolated only to fit an STL seasonal–trend decomposition; the original observations are retained for filtering. Outliers are identified iteratively, for up to five STL fits, using a modified z-score of the residuals with a threshold of 4.5. • Manual Inspection. An expert hydrology scientist manually inspects the remaining 581 time series and discards those that contain clear faulty measurements: 375 time series remain. • Datum Alignment. Satellite sources report WSE w.r.t. a reference geoid (Pavlis et al. 2012), i.e., the elevation w.r.t. an approximate sea level. On the other hand, the in situ time series are measured in reference to a local gauge datum, for instance w.r.t. the local river bed bottom, which does not include the terrain elevation. Hence, before comparison, in all evaluations, we follow common practice in hydrology (Halicki et al. 2026) to project the in situ measurements with a per-location fitted linear regression to the reference scale of the satellite data and/or predictions. This operation is considered part of our eval-
uation benchmark, will be part of published evaluation code, and should be replicated by future authors that use this benchmark dataset. Matching to SWORD Reaches. Figure 14 shows matched and unmatched locations per source. • HydroWeb virtual stations, ICESat-2 transects, and ANA gauges are matched to their nearest SWORD reach using a 10 km distance threshold (adapted from the 20 km threshold Halicki et al. (2026) uses for the Solimoes, which is part of the Amazon basin). Time series that are further than 10 km from the nearest SWORD reach are not used (though they are still included in the data files). • SWOT is linked directly by the RiverSP reach identifier.
A.3
Storage Format
Each source is stored in a netcdf file with h5netcdf as a CF 1.13 indexed-ragged time series (Eaton et al. 2025), which means they are stored along a location dimension and an observation dimension (rather than a time dimension), where time is an additional variable. This means only measured location–time pairs are stored, and missing (unobserved) location–time pairs do not take up storage space. Table 6 summarizes the shared schema. An example netcdf file per source is included in the Code and Data Supplement along with the submission. Indexing and Validity. The index variable declares instance_dimension=location. Locations are unique and sorted by location_id; observations are sorted first by location and then by time. Duplicate (location_id, time) pairs are forbidden. A quality flag of 0 marks a recorded but rejected value, a flag of 1 marks an accepted value, and a missing observation has no corresponding row. Global Metadata. Each file records schema_version=1.0, featureType=timeSeries, a non-empty source, and quality_convention="quality_flag: 0=rejected, 1=accepted". Source-specific extensions. Sources may add variables without changing the two-dimensional organization: • SWOT adds observation-level width, width_u, slope, and reach_q_b, together with int32 cycle_id and pass_id. • ICESat-2 adds float32 ht_stdev and slope_cross_error (m), as well as int32 counts of transects, beams, and their quality categories. • HydroWeb adds observation-level satellite, orbit, and retracking_algorithm metadata. • ANA adds the location-level string river_name.
A.4
Hyperparameters and Model Details
Table 7 gives an overview of the used hyperparameters. Importantly, since none of the models sees in situ data during training, and since our model makes source-aware predictions (i.e., with a separate linear decoding head per source),
we make predictions as if for one of the training sources, also when we compare to in situ data. We de-normalize the prediction with the interpolated statistics from Figure 10: since the source that is being decoded might not have measurements hence also no statistics on the reach that is being predicted and statistics of the in situ gauge series should not be used. We find that decoding as HydroWeb works best when comparing to in situ gauge data, presumably since this source has long records of clean data. This creates a discrepancy: we decode normalized HydroWeb predictions, and de-normalize them with interpolated SWOT statistics, but of all the tested combinations, this gave best results.
A.5
Baselines
All baselines use the splits, masking conventions, and sampled-subgraph settings, and hyperparameters from Table 7 unless stated otherwise. Subgraphs vs full graph. We distinguish fixed-graph and sampled-subgraph settings. Fixed graphs use a basin-wide time–node grid with stable reach identities, following the transductive protocols of GRIN, SPIN-H, and ImputeFormer (Cini, Marisca, and Alippi 2022; Marisca, Cini, and Alippi 2022; Nie et al. 2024). Sampled sub-graph windows instead contain the same local river neighborhoods as our method (cfr. Table 7 and Figure 3). SPIN-H and IGNNK explicitly support subgraph sampling (Marisca, Cini, and Alippi 2022; Wu et al. 2021), while for GRIN and ImputeFormer we run the subgraph version as adaptations rather than replications of their published protocols. KITS is closest to its published increment-training setting when run on a fixed graph with virtual nodes (Xu et al. 2025). Source representation. We consider two representations of co-located sources (SWOT, HydroWeb, ICESat-2): the channels mode assigns each source a separate value (in a separate input channel), input mask, and target mask, preserving source identity and allowing one source to remain visible while another is held out. The merged mode first normalizes sources with common SWOT-derived reach statistics; after masking, visible source values are averaged into one scalar. Thus, held-out observations never enter the aggregate, but source identity and source-specific biases are removed. Models are never trained on in situ data. At in situ gauges, we use the training source predictions and dense SWOTderived reach statistics. IGNNK and KITS use merged inputs because they are scalar-WSE models (Wu et al. 2021; Xu et al. 2025). For ImputeFormer, GRIN, and SPIN-H we implement channel variants, even though these methods were primarily evaluated with scalar targets (Nie et al. 2024; Cini, Marisca, and Alippi 2022; Marisca, Cini, and Alippi 2022). We also tested representing each source–location pair as a separate node (allowing thus several nodes for the same reach), but it gave worse performance. Non-graph controls. The non-parametric kNN baseline uses 91-day windows, ten spatial neighbors, and equal spatial and temporal weights. It fits a dense field separately per input source and outputs the mean of the per-source predictions.
Figure 13: Filtered (blue) and unfiltered (orange) in situ gauge time series: a clean time series (top, left), a faulty time series that shows sudden jumps that clearly do not represent natural variations in river WSE (top, right), and two time series where automatic filters removed outliers (bottom, left and right). A temporal-only, one-layer bidirectional LSTM inspired by BRITS-I (Cao et al. 2018) uses the channels mode and 365day per-reach windows, without graph information. Graph-imputation models. For fixed-graph runs we compared 7-, 14-, and 28-day fixed windows (smaller than subgraph runs to fit in memory), as well as merged vs. channels mode, and chose the setting that gave best RMSE on validation samples. We had to reduce the hidden dimension size of some baselines to fit their training in our GPU memory. We reuse settings from the original publications and/or default settings in the published repositories where applicable. GRIN (Cini, Marisca, and Alippi 2022) uses one recurrent graph-imputation layer, hidden size 64, kernel size 2, and decoder order 1. Its fixed and sampled variants use channels and windows of 28 and 91 days and feed-forward widths of 128 and 64, respectively. SPIN-H (Marisca, Cini, and Alippi 2022) uses five layers, hidden size 32, η = 3, and two message-passing layers. Its fixed variant uses 28-day merged inputs, latent size 128, and four heads; its sampled variant uses 91-day channels, latent size 64, and two heads. IGNNK (Wu et al. 2021) uses hidden size 192, first-order diffusion, and merged windows of 28 days (fixed graph) or 91 days (sampled subgraph). Its graph convolutions reconstruct each time step from visible neighbors without explicit temporal dynamics. KITS (Xu et al. 2025) uses hidden size 64, a 0.3 virtual-node ratio, unit cycle-consistency weight, and merged windows of 7 days (fixed) or 91 days (sampled). Training-time virtual nodes are connected around random one-hop neighborhoods and removed before evaluation. ImputeFormer (Nie et al. 2024) uses input, node-embedding,
and feed-forward dimensions of 64, 96, and 256; three layers; four temporal heads; rank 8; dropout 0.1; and Fourier-loss weight 0.01. Its fixed variant uses a 14-day window and its sampled subgraph variant uses a 91-day window, both with source channels. Repositories. For GRIN, SPIN-H and the BRITS-I inspired bidirectional LSTM we use the implementations of torch-spatiotemporal.4 ImputeFormer5 , IGNNK6 and KITS7 are adapted from their published repositories. Baselines that were not included. GgNet (De Felice et al. 2024) is excluded because it assumes dynamic covariates available at target reaches, unlike our sparse observations of the WSE target itself. Diffusion imputers such as PriSTI and FastSTI (Liu et al. 2023a; Cheng et al. 2024) are deferred pending a matched, deterministic extreme-sparsity protocol. Forecasting-only methods are excluded because they provide neither an imputation objective nor a corresponding masking protocol. Reach-Reg. Reach-Reg is a recent method that leverages the SWOT measurement geometry to fit a chain of linear regressions (orthogonal distance regressions) to propagate measurements between consecutive reaches. Subsequently, a time lag is estimated with the Manning formula for river velocity, and aggregation and interpolation then yield a daily 4
https://github.com/TorchSpatiotemporal/tsl https://github.com/tongnie/ImputeFormer 6 https://github.com/Kaimaoge/IGNNK 7 https://github.com/Sam1224/KITS 5
Figure 14: Spatial distribution of data sources in AmazonWSE, matched (kept) and unmatched (discarded) to SWORD reaches within 10 km (before quality filtering). From left to right: HydroWeb virtual stations, ICESat-2 transects, and ANA in situ gauges. Variable
Type
Description
Location variables (location) location_id location_quality_flag latitude longitude sword_reach_id reach_distance_m
int64 int8 float32 float32 int64 float32
Unique site identifier (cf_role=timeseries_id) Site validity: 0 rejected, 1 accepted Latitude (degrees_north) Longitude (degrees_east) Associated SWORD reach identifier Distance to the associated reach (m)
Observation variables (observation) observation_location_index int32 time datetime64[ns] wse float32 wse_u float32 quality_flag int8
Index into the location dimension Observation time (CF convention) Water-surface elevation (m) Water-surface elevation uncertainty (m) Observation validity: 0 rejected, 1 accepted
Table 6: NetCDF schema shared by all data sources. Variables are grouped by their associated dimension. WSE. The approach is applied either to combined SWOT, S3 and S6 measurements from the DAHITI dataset (Schwatke et al. 2015) or to SWOT RiverSP measurements only. We use the official Reach-Reg implementation.8 We group the 375 gauges into 107 unique river segments (sequences of river reaches formed by traversing downstream and upstream from the target gauge, choosing the upstream branch with highest accumulated flow). We run Reach-Reg with the settings of Halicki et al. (2026) on SWOT RiverSP data (version D, as our model and baselines) retrieved and filtered from the HydroCron API with the settings and filters proposed by Reach-Reg. We compare linear with Akima interpolation and find that it does not make a significant difference. A prediction is not made for all gauges, since Reach-Reg does not make a prediction when quality filters leave insufficient neighboring stations for propagation or when propagation and path-error filtering leave no valid predictions. The authors achieve a better accuracy when running Reach-Reg on timeseries from the DAHITI dataset, combining S3, S6 and SWOT (Schwatke et al. 2015), rather 8 https://github.com/MichalHalicki4/Reach-Reg/, commit 59aa780.
fetched
at
than on SWOT only, but since SWOT timeseries for only a small part of reaches in the Amazon have been included in DAHITI as of yet, this further reduces coverage by over %50. We next adapt Reach-Reg to make hindcast predictions for the available gauges (i.e., for the period from 2016-01 to 2022-06-30 for which no SWOT measurements are available). We reuse the fitted orthogonal distance regression weights and Manning parameters that resulted from the prediction that was just described for SWOT RiverSP data. We then inject historical DAHITI observations into the matching (and fitted) SWOT-era station slots and propagate them according the fitted and chosen parameters. We increase the maximum cumulative ODR path-error, defined as the sum of the ODR RMSEs along a propagation path, from 10 to 50 m. We also increase the per-link ODR-RMSE threshold used to choose a direct link versus an alternative path from 0.2 to 1.0 m. We retain a propagated observation only when its cumulative path-error is below max(0.2 · WSE amplitude, 1.5 · median(cumulative path-error)). These relaxations accommodate the longer regression chains needed in the much sparser hindcast setting. It should be noted that this evaluation takes Reach-Reg out of the context it was developed for, and that results are hence only indicative.
Group
Hyperparameter
Value
Architecture
Temporal model Bidirectional; decoding order Embedding / SSM state dimension Bidirectional SSM layers Step rank; convolution width Dropout rate Metadata combination; reinjection Output heads Tree embedding dimension F (Shiv and Quirk 2019)
Mamba-1 true; time then flow 192 / 16 3 12; 4 0.5 elementwise addition; true satellite-specific 4
Optimization
Loss Optimizer Learning rate; weight decay Batch size; maximum steps Warmup steps; gradient clipping Learning-rate schedule Validation interval; early-stopping patience Early-stopping metric Seed
MSE AdamW (Loshchilov and Hutter 2019) 10−4 ; 0.05 16; 100,000 1,000; 1.0 reduce on plateau (factor 0.2, patience 5) 500 steps; 20 checks RMSE on location masked HydroWeb 43
Input subgraphs
Bin size; temporal bins Minimum / maximum measurements Maximum spatial distance Upstream / downstream sampling Main-upstream trunk proportion Maximum upstream / downstream hops WSE normalization
24 hours; 91 (days) 15 (training), 0 (eval) / 500 300 km 0.75 / 0.25 0.33 30 / 30 per-location mean and standard deviation
Masking
Training strategies (probabilities) Mask ratio
location (0.9), random (0.1) 0.66
Table 7: Hyperparameters used for the SSM model.
A.6
Predictions
Figure 15 shows examples of predicted time series.
Figure 15: Example predictions of our method compared to Reach-Reg on SWOT-era (top) and on hindcast (bottom, where Reach-Reg does not have enough input observations to predict the whole series). The dotted lines represent the in situ gauge series linearly projected to the datum of our predictions and of Reach-Reg.
Appendix References Eaton, B.; Gregory, J.; Drach, B.; Taylor, K.; Hankin, S.; Caron, J.; Signell, R.; Bentley, P.; Rappa, G.; Höck, H.; Pamment, A.; Juckes, M.; Raspaud, M.; Blower, J.; Horne, R.; Whiteaker, T.; Blodgett, D.; Zender, C.; Lee, D.; Hassell, D.; Snow, A. D.; Kölling, T.; Allured, D.; Jelenak, A.; Soerensen, A. M.; Gaultier, L.; Herlédan, S.; Manzano, F.; Bärring, L.; Barker, C.; Bartholomew, S. L.; Lavergne, T.; Lawrence, B.; Massey, N.; Cofiño, A. S.; McGinnis, S.; and Laake, P. V. 2025. NetCDF Climate and Forecast (CF) Metadata Conventions. https://doi.org/10.5281/zenodo.14274886. Loshchilov, I.; and Hutter, F. 2019. Decoupled Weight Decay Regularization. In International Conference on Learning Representations. Verdin, K. L.; and Verdin, J. P. 1999. A topological system for delineation and codification of the Earth’s river basins. Journal of Hydrology, 218(1-2): 1–12.