# A Kalman Filter Approach for Estimating Daily Discharge Using Space-Based Discharge Estimates

Siqi Ke, Mohammad J Tourian, Nico Sneeuw, Renato Prata de Moraes
Frasson, Rodrigo C D Paiva, Michael Durand, Colin Gleason, Omid Elmi,
Pierre-olivier Malaterre, Cédric David

► To cite this version:

Siqi Ke, Mohammad J Tourian, Nico Sneeuw, Renato Prata de Moraes Frasson, Rodrigo C D Paiva, et al.. A Kalman Filter Approach for Estimating Daily Discharge Using Space-Based Discharge Estimates. Water Resources Research, 2025, 61 (8), pp.e2024WR037667. 10.1029/2024WR037667. hal-05247098

HAL Id: hal-05247098

https://hal.inrae.fr/hal-05247098v1

Submitted on 9 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2024WR037667

### Special Collection:

Science from the Surface Water and Ocean Topography Satellite Mission

### Key Points:

* A Kalman filter-based method is developed for estimating daily discharge over continuous river reaches in a single-branch river

* Our method enables data assimilation along the river using a physically based correlation model and SWOT observations

* Our results show that estimating daily discharge from SWOT data is feasible

### Supporting Information:

Supporting Information may be found in the online version of this article.

### Correspondence to:

S. Ke,

siqi.ke@gis.uni-stuttgart.de

### Citation:

Ke, S., Tourian, M. J., Sneeuw, N., Frasson, R. P. D. M., Paiva, R. C. D., Durand, M., et al. (2025). A Kalman filter approach for estimating daily discharge using space-based discharge estimates. *Water Resources Research*, *61*, e2024WR037667. https://doi.org/10.1029/2024WR037667

Received 2 APR 2024

Accepted 11 JUL 2025

© 2025. The Author(s).

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## A Kalman Filter Approach for Estimating Daily Discharge Using Space-Based Discharge Estimates

Siqi Ke¹ ID, Mohammad J. Tourian¹ ID, Nico Sneeuw¹ ID, Renato Prata de Moraes Frasson² ID, Rodrigo C. D. Paiva³ ID, Michael Durand⁴ ID, Colin Gleason⁵ ID, Omid Elmi¹ ID, Pierre-Olivier Malaterre⁶ ID, and Cédric David² ID

¹Institute of Geodesy, University of Stuttgart, Stuttgart, Germany, ²Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA, ³Instituto de Pesquisas Hidráulicas, Universidade Federal do Rio Grande do Sul, Porto Alegre, Brazil, ⁴School of Earth Sciences, Ohio State University, Columbus, OH, USA, ⁵Department of Civil and Environmental Engineering, University of Massachusetts Amherst, Amherst, MA, USA, ⁶INRAE, UMR G-eau, Montpellier, France

**Abstract** The SWOT satellite mission is the first to conduct a global survey of the Earth's surface waters, measuring water surface height, river width, and water surface slope, based on which river discharge is estimated. At mid-latitudes, the repeat orbit design of SWOT only allows a sampling of twice per repeat cycle, which is considered too low for most hydrological applications. To address the spatiotemporal limitations of SWOT, we develop a method that assimilates SWOT observations across continuous reaches within a single-branch river network to obtain daily discharge estimates. Our model-free assimilation method provides a linear dynamic system that includes a process model based on a physically based spatiotemporal discharge correlation model and observation equations utilizing SWOT products. We solve this dynamic system using a simple Kalman filter in the time domain, assimilating SWOT observations and incorporating the physically based prior to estimate daily discharge. Since SWOT discharge products were not yet available during the period of this research, we used synthetic SWOT data sets, introducing random and systematic errors through Monte Carlo simulation. The validation of the estimated discharge against true discharge over all test cases leads to a median correlation as high as 0.95, a median NSE for residuals (mean-removed discharge) as high as 0.81, and a median relative bias as low as 5%, respectively. These promising results suggest that daily discharge for continuous reaches in a river network can be obtained through our data assimilation framework.

## 1. Introduction

River discharge, one of the most crucial hydrological variables, serves as a major link in the global hydrologic cycles as well as being a critical climate component (Dingman & Bjerklie, 2006). Thus, measurements of river discharge are essential for flood hazard management, water resource planning, climate and ecology studies, and compliance with transboundary water agreements (Smith & Pavelsky, 2008). However, the number of gauge station records has consistently decreased over the past few years due to the lack of an adequate monitoring network and restricted access (Elmi et al., 2024; Fekete & Vörösmarty, 2007; Saemian et al., 2025; Tourian et al., 2013).

To address this problem, the interest in using spaceborne approaches for estimating discharge has significantly risen over the last 2 decades (Alsdorf et al., 2007; Domeneghetti et al., 2021; Elmi et al., 2021; Gleason & Durand, 2020; McCabe et al., 2017; Tourian, Elmi, et al., 2017; Tourian, Schwatke, & Sneeuw, 2017). Several satellite-based initiatives have been developed to support inland water monitoring. For instance, the Global Water Watch platform provides free, near-real-time satellite-based estimates of reservoir water area and storage globally (Deltares & World Resources Institute, 2023). In addition, platforms like HydroWeb (CNES, LEGOS, and Theia, 2013), DAHITI (DGFI-TUM, 2013), and HydroSat (Tourian et al., 2022) provide long-term, altimetry-based water level time series and other hydrological information for lakes, reservoirs, and rivers.

Launched on 16 December 2022, the Surface Water and Ocean Topography (SWOT) satellite allows estimates of global river discharge for the first time, vastly improving the understanding of global hydrological processes (Biancamaria et al., 2016). SWOT was collaboratively developed by NASA, the Center National d'Etudes Spatiales (CNES, the French space agency), the Canadian Space Agency (CSA), and the United Kingdom Space Agency (UKSA) to make the first global survey of Earth's surface water. SWOT simultaneously measures river water surface height, river width, and water surface slope (JPL Internal Document, 2020) to enable discharge
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Conceptual overview of this study. SWOT measurements are temporally sparse, illustrated by the mix of observed (blue boxes) and unobserved (yellow boxes) days for each reach (left). Through the Kalman filter-based assimilation process, observational information is propagated to both upstream and downstream reaches, as shown by curved arrows (center). This enables reconstruction of daily discharge across all continuous river reaches (right), where solid green boxes represent days assimilated by SWOT measurements, and hatched green boxes indicate discharge estimates based mainly on the process model of our Kalman filter.</figcaption></figure>

estimates for river reaches with width down to 100 m or maybe even 50 m (Pavelsky et al., 2014), allowing an expected discharge uncertainty of less than 30% in ungauged basins (Durand et al., 2023).

While SWOT can provide consistent and coherent information about the spatial distribution of river storage and discharge (Biancamaria et al., 2016), it still has limitations in observing the global spatiotemporal dynamics of the continental water surface (Alsdorf et al., 2007). SWOT measures most mid-latitude locations twice on average during the 21-day repeat cycle of the science orbit (~35 observations per year), with more observations at higher latitudes (Biancamaria et al., 2016; Durand et al., 2023). Although the sampling is still appropriate to monitor the global water cycle, it is inadequate for many local-scale questions on rivers where SWOT may not fully observe temporal dynamics (Durand et al., 2016). Moreover, SWOT does not observe full river networks. According to Altenau et al. (2021), a total of 8,458 river locations (3.4% of the total) are never observed due to small gaps in SWOT coverage, and discharge estimates are not produced for locations narrower than 50 m Durand et al. (2023).

To address the spatiotemporal limitations of SWOT, this study proposes a model-free data assimilation method to estimate daily river discharge for continuous river segments, that is, reaches, within a river network using SWOT measurements along the river (see Figure 1). To achieve this, SWOT measurements from different reaches, acquired at different times, must be combined. Such a combination leads to an improved temporal resolution by propagating observations upstream and downstream to a target reach. To this end, we implement a linear dynamic system that provides a scheme for data assimilation of multiple SWOT output types along a river (Tourian, Schwatke, & Sneeuw, 2017). Within the linear dynamic system, a discharge state vector evolves in the time domain through a stochastic process model. This stochastic process model is established with the help of the auto- and cross-covariance information of discharge across the catchment, which is computed from a physically based spatiotemporal discharge correlation model and the discharge standard deviation information from the prior SWOT database, providing hydrologic data for global rivers over 30 m wide (Altenau et al., 2021). Moreover, the linear dynamic system contains an observation model with two equations combining the evolved state vector with SWOT output elements.

We solve the linear dynamic system using a Kalman filter, which can ultimately provide daily discharge estimates over all reaches in a single-branch river network. Kalman filter is a sequential data assimilation algorithm that provides estimates of unknown variables using measurements observed over time. Kalman filter has been demonstrating its usefulness in various applications and has a relatively simple form, which requires small computational power (Kim & Bang, 2018). The Kalman filter algorithm is mainly divided into two main steps: prediction and update. The prediction step uses the process model to advance the state temporally until the next scheduled time step, resulting in a predicted (a priori) state. Once the observation is available at a specific measurement epoch, the update step incorporates it to refine the current predicted state through the observation model. The improved state is referred to as the updated (a posteriori) state. Due to the unavailability of some
required SWOT products at the time of this study, including the discharge product (expected to be globally available in 2025), we employed synthetic SWOT data to evaluate the feasibility of our approach.

In related studies in the recent past, different methodologies were developed to construct a spatiotemporally complete estimate of discharge from remote sensing observations. Table 1 presents a tabulated literature review of some of these methods, categorized by their application in either gauged or ungauged basins.

In line with the previous developments, here, for the first time, we aim to achieve a daily estimation of river discharge across a single-branch river network by assimilating SWOT products using the Kalman filter, thereby improving the spatiotemporal resolution of discharge estimates.

The paper is organized as follows. We first introduce the test data sets used in this study (see Section 2). Then we describe the structure of our linear dynamic system, including both the process and observation model, as well as the Kalman filter and the validation method (see Section 3). The estimated daily discharge is presented subsequently, accompanied by their validation results (see Section 4).

## 2. Data

To evaluate the performance of our proposed algorithm, we use a perturbed version of synthetic SWOT data sets, hereafter referred to as the *perturbed data sets* (Frasson et al., 2020). This choice is due to the unavailability of certain SWOT products during the study period, particularly the discharge product, which is expected to be released globally around mid-2025. These data sets consist of 15 river sections (see Table 2) and are commonly used to evaluate the performance of SWOT discharge algorithms. They are generated using the SWOT high-fidelity simulator with a parameterized error model originally developed by Durand et al. (2020) and subsequently adopted in several peer-reviewed studies (Frasson et al., 2021, 2023), thus providing a widely accepted proxy in the field. Due to the large computational demand, the Ohio (Adams et al., 2010) and the Missouri River (USACE, 2018) were split into seven and two test cases featuring areas with distinct hydraulic characteristics (Frasson et al., 2023). Figure 2 visualizes the study area, showing river width information in cyan bands for each test case and additionally providing a global overview of river locations.

In the data sets, the measurements are in days at two spatial scales: river cross-section and reach. A river cross-section is defined along the river centerline, usually with less than 1 km distance, and serves as an approximation of SWOT nodes, typically spaced 200 m apart (Altenau et al., 2021). A reach is defined as a consecutive river segment between cross-sections, each spanning up to 20 km or more, whereas a SWOT-defined reach is approximately 10 km in length. The test cases have 4–14 reaches without tributaries, with their boundaries delineated by red lines in Figure 2. It is important to note that discharge estimation presents challenges in small reaches around low-head dams and other obstacles, which are typically classified as Type 4 reaches (dam and waterfall) in the SWOT river database (Altenau et al., 2021). Consequently, these reaches have been identified as *bad reaches*, where the perturbed data sets contain no measurements, while other reaches are referred to as *good reaches* (Frasson et al., 2017). Table 2 lists the total number of reaches and days, as well as the indices of bad reaches.

In this study, we predominantly focus on reach-scale data, as the discharge estimation is at the reach scale, consistent with the SWOT discharge product. Our approach incorporates two reach-scale measurements: change in cross-sectional area and discharge. The change in cross-sectional area, one of the SWOT outputs and an input for our observation model, represents the anomaly between the total cross-sectional area and a baseline reference.

In the perturbed data sets, 10,000 realizations of the daily reach-scale measurements are generated via a Monte Carlo simulation (Frasson et al., 2023). For the computation of perturbed discharge, the flow law parameter inversion algorithm Metropolis-Manning (MetroMan) is employed (Durand et al., 2014), which is one of the SWOT discharge inversion algorithms. It uses a Metropolis algorithm (a Bayesian sampling approach) to derive posterior distributions of the flow law parameters and the modified Gauckler-Manning-Strickler equation (Durand et al., 2023) to compute discharge. Both the perturbed change in cross-sectional area and perturbed discharge are subject to bias. For the change in cross-sectional area, the bias arises from layover effects in height measurements (Durand et al., 2020). For the discharge, the bias originates from two sources: the biased height, as in the change in cross-sectional area, and biased flow law parameter estimates resulting from the Bayesian sampling approach. The detailed procedure for generating perturbed reach-scale measurements, along with the associated biases, is introduced in Appendix A.
<caption>**Table 1**</caption>

Tabulated Literature Review for Discharge Estimation From Remote Sensing Observations, Categorized by Their Application in Gauged and Ungauged Basins

<table><thead><tr><th>Reference</th><th>Methodology</th></tr></thead><tbody><tr><td colspan="2">in ungauged basins</td></tr><tr><td>Yoon et al. (2012)</td><td>Estimating river bathymetry for retrieving river discharge from SWOT using an ensemble Kalman filter-based data assimilation algorithm coupled with a hydrodynamic model.</td></tr><tr><td>Durand et al. (2014)</td><td>Using a Bayesian algorithm to estimate discharge from water surface elevation and slope measurements based on the diffusive approximation of the shallow water equations.</td></tr><tr><td>Sun et al. (2015)</td><td>Estimating daily time series of streamflow using a hydrological model calibrated based on satellite observations of river water surface width.</td></tr><tr><td>Oubanas et al. (2018)</td><td>Using a variant of the classical variational data assimilation method to estimate discharge from simulated observations of SWOT.</td></tr><tr><td>Brinkerhoff et al. (2020)</td><td>Improving discharge estimation by classifying and parameterizing rivers based on their unique geomorphology and hydraulics through an update to the Bayesian At-many-stations hydraulic geometry-Manning's (BAM) algorithm (Hagemann et al., 2017).</td></tr><tr><td>Andreadis et al. (2020)</td><td>Developing an ensemble Kalman filter-based data assimilation discharge algorithm by integrating a data-driven approach to estimate priors with a model informed by hydraulic geometry relations.</td></tr><tr><td>Larnier et al. (2020)</td><td>Using an inversion method to estimate discharge by combining the Saint-Venant equations with a dedicated algebraic system.</td></tr><tr><td>Frasson et al. (2021)</td><td>Estimating discharge relying on extrapolating geometric and hydraulic characteristics derived from large data sets in conjunction with SWOT data.</td></tr><tr><td>Ishitsuka et al. (2021)</td><td>Providing river discharge estimation by integrating Mass-conserved Flow Law Inversion (McFLI) estimated by satellite imagery data with hydrological models within a globally applicable data assimilation framework.</td></tr><tr><td>Gejadze et al. (2022)</td><td>Developing a robust discharge estimation method using Variational Data Assimilation based on the full Saint Venant hydraulic model.</td></tr><tr><td colspan="2">in gauged basins (incl. Sparsely gauged basins)</td></tr><tr><td>Neal et al. (2009)</td><td>Using an ensemble Kalman filter to assimilate water level information with simulations from a coupled hydrological and hydrodynamic model to estimate discharge, given information about the channel.</td></tr><tr><td>Nathanson et al. (2012)</td><td>Using a fluid mechanics-based model, constrained by airborne LiDAR-derived topographic data, to define stage-discharge rating curves for stream discharge estimation in channels in remote locations.</td></tr><tr><td>Brakenridge et al. (2012)</td><td>Applying the calibration-measurement ratio approach (Brakenridge et al., 2007) to thousands of global stations, using a previously in situ calibrated hydrology model to provide the training data for discharge predictions.</td></tr><tr><td>Paiva et al. (2013)</td><td>Introducing a data assimilation framework for gauged and radar altimetry-based discharge and water levels applied to a large-scale hydrologic-hydrodynamic model for stream flow forecasts over the Amazon River basin.</td></tr><tr><td>Lorenz et al. (2015)</td><td>Presenting an ensemble Kalman filter to estimate river discharge through least squares prediction, leveraging the joint temporal and spatial covariance of precipitation, evapotranspiration, water storage changes, and river discharge.</td></tr><tr><td>Huang et al. (2018)</td><td>Demonstrating the feasibility of using satellite imagery and altimetry data with a gauge on the information-sparse Tibetan plateau to estimate river discharge by employing power function equations.</td></tr><tr><td>Emery et al. (2018)</td><td>Introducing an ensemble Kalman filter-based data assimilation platform that calibrates the hydrological model of river storage and discharge using satellite altimetry-based discharge products and in situ data.</td></tr><tr><td>Lin et al. (2019)</td><td>Integrating advanced ground and remote sensing precipitation data, machine learning-derived runoff maps, and precise river flowlines from a detailed digital elevation model to produce high-resolution global river discharge.</td></tr></tbody></table>
<table><caption>**Table 1**
*Continued*</caption><thead><tr><td>Reference</td><td>Methodology</td></tr></thead><tbody><tr><td>Emery et al. (2020)</td><td>Assimilating daily averaged in situ discharge measurements through a classical Kalman filter employing a hydrodynamic model based on David et al. (2011) to correct daily averaged runoff.</td></tr></tbody></table>

To summarize, Table 3 specifies the variables in the perturbed data sets used in this study. We separate data into two categories: river information and time series data. *River information* provides fundamental and general river characteristics, whereas *time series data* offers daily reach-scale data. The daily data includes perturbed measurements, available only at good reaches, and the true discharge, which later serves as a key benchmark for validating our results. In the Monte Carlo simulation, there are 100 sets of flow law parameters generated by MetroMan with varying parameter errors; each set of parameters is associated with 100 distinct realizations of the measurement process (Frasson et al., 2023). Hence, there are a total of 10,000 draws of perturbed change in cross-sectional area and perturbed discharge, respectively.

For example, Figures A1 and A2 in Appendix A present histogram plots, including the realizations of flow law parameters and measurements for a reach in the test case Missouri Upstream, highlighting the existence and sources of bias in both perturbed change in cross-sectional area and perturbed discharge. Moreover, Table A1 in the same appendix summarizes the median relative bias of these two measurements across all reaches for each test case. The relative bias of the change in cross-sectional area ranges from 3% to 72%, while that of the discharge spans from 33% to 494%, with most cases falling below 100%.

Note that, although the perturbed data sets can offer daily measurements, our goal is to estimate discharge benefiting from SWOT data and imitate the actual SWOT observation sampling. Thus, we employ the CNES Large Scale SWOT Hydrology Simulator (Elmer et al., 2020) to simulate the SWOT orbit pass plans for each test case and then sample the measurements based on them to serve as SWOT data used in this study. In the simulation, we set 1 January 2023, as the simulation start day and the simulation end day depending on the number of days in each test case. Table 2 additionally lists the total number of measured days across all reaches during the simulation period.

<table><caption>**Table 2**
*List of the Test Cases in This Study Published by Frasson et al. (2021)*</caption><thead><tr><td>No.</td><td>Test case name</td><td>Number of reaches/days</td><td>Total number of meas.</td><td>Indices of bad reaches</td><td>Reference</td></tr></thead><tbody><tr><td>1</td><td>Brahmaputra</td><td>5/107</td><td>89</td><td>–</td><td>Maswood and Hossain (2016)</td></tr><tr><td>2</td><td>Iowa River</td><td>10/366</td><td>380</td><td>Reach 2,4</td><td>Gilles et al. (2012)</td></tr><tr><td>3</td><td>Jamuna</td><td>4/250</td><td>168</td><td>–</td><td>Maswood and Hossain (2016)</td></tr><tr><td>4</td><td>Kushiyara</td><td>5/366</td><td>227</td><td>–</td><td>Maswood and Hossain (2016)</td></tr><tr><td>5</td><td>Missouri Downstream</td><td>5/595</td><td>256</td><td>–</td><td>USACE (2018)</td></tr><tr><td>6</td><td>Missouri Upstream</td><td>4/595</td><td>116</td><td>–</td><td>USACE (2018)</td></tr><tr><td>7</td><td>Ohio Section 1</td><td>10/220</td><td>299</td><td>Reach 6</td><td>Adams et al. (2010)</td></tr><tr><td>8</td><td>Ohio Section 2</td><td>8/220</td><td>197</td><td>–</td><td>Adams et al. (2010)</td></tr><tr><td>9</td><td>Ohio Section 3</td><td>14/220</td><td>528</td><td>–</td><td>Adams et al. (2010)</td></tr><tr><td>10</td><td>Ohio Section 4</td><td>5/220</td><td>99</td><td>–</td><td>Adams et al. (2010)</td></tr><tr><td>11</td><td>Ohio Section 5</td><td>8/220</td><td>295</td><td>–</td><td>Adams et al. (2010)</td></tr><tr><td>12</td><td>Ohio Section 7</td><td>7/220</td><td>234</td><td>–</td><td>Adams et al. (2010)</td></tr><tr><td>13</td><td>Ohio Section 8</td><td>8/220</td><td>172</td><td>Reach 2,6</td><td>Adams et al. (2010)</td></tr><tr><td>14</td><td>Padma</td><td>5/366</td><td>260</td><td>–</td><td>Maswood and Hossain (2016)</td></tr><tr><td>15</td><td>Seine Downstream</td><td>6/365</td><td>208</td><td>Reach 3,5</td><td>Even et al. (2007)</td></tr></tbody></table>

*Note.* The total number of measurements (meas.) across all observed reaches is according to the pass plans generated by the CNES Large Scale SWOT Hydrology Simulator (Elmer et al., 2020). Bad reaches refer to areas around low-head dams and other obstacles and thus have limitations in discharge estimation. They are enumerated separately and are helpful for our result assessment later on.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Study area maps of the 15 data sets. River locations are shown individually in each subplot, with width information (cyan bands) and reach boundaries (red lines). The final subplot provides a global overview of the river distribution, with marker numbers corresponding to the river indices in the preceding subplots.</figcaption></figure>

### 3. Methodology

To estimate daily discharge for all the reaches in a single-branch river network, we develop a linear dynamic system featuring a stochastic process model and observation equations. The system is then solved using a Kalman filter. A summary of our methodology is provided in Figure 3. In the following, we describe the construction of

<table><caption>**Table 3**
*List of Variables Used in This Study*</caption><thead><tr><th>Variable</th><th>Unit</th><th>Dimension</th><th>Description</th></tr></thead><tbody><tr><td colspan="4">River information</td></tr><tr><td>name</td><td>—</td><td>1</td><td>River name</td></tr><tr><td>rch_bnd</td><td>m</td><td>R+1,1</td><td>Reach boundaries</td></tr><tr><td>gdrch</td><td>—</td><td>GR, 1</td><td>Indices of good reaches</td></tr><tr><td colspan="4">Time series data</td></tr><tr><td>t</td><td>day</td><td>1, T</td><td>Time measured in days</td></tr><tr><td>Q<sub>true</sub></td><td>m<sup>3</sup>/s</td><td>R, T</td><td>True discharge</td></tr><tr><td>Q<sub>ptb</sub></td><td>m<sup>3</sup>/s</td><td>GR, T, 10 000</td><td>Perturbed discharge</td></tr><tr><td>δA<sub>ptb</sub></td><td>m<sup>2</sup></td><td>GR, T, 10 000</td><td>Perturbed change in cross-sectional area</td></tr><tr><td colspan="4">T: Number of time steps, R: Number of reaches, GR: Number of good reaches</td></tr></tbody></table>

*Note.* The data sets are generally divided into two categories: river information and time series data. In *time series data*, data mainly contains two components: true discharge and perturbed measurements. For each reach, there are 10,000 realizations of perturbed measurements per time step.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Flowchart of the proposed methodology. First, we develop a linear dynamic system that incorporates a stochastic process model and two observation equations, utilizing the perturbed data sets. Then, this system is solved using the Kalman filter, which enables the estimation of daily discharge time series across continuous river reaches.</figcaption></figure>

the process model (see Section 3.1) and observation equations (see Section 3.2), as well as the implementation of the Kalman filter to obtain sequential solutions (see Section 3.3).

### 3.1. Process Model

We conceptualize the linear process model as follows (Kalman, 1960):

$$r_t = \Phi r_{t-1} + e \quad e \sim N(0, W) \quad (1)$$

where

$$r_t = q_t - \bar{q}.$$

The state vector $r_t$ incorporates discharge anomalies, defined as the difference between the discharge time series of each reach within a river network $q_t$ and its respective long-term mean discharge $\bar{q}$ (Hofmann-Wellenhof & Moritz, 2005; Tourian, Schwatke, & Sneeuw, 2017). This linear process model defines the evolution of the state vector from one time step $t-1$ to the next at time $t$, that is, $r_{t-1}$ to $r_t$, by applying the state transition matrix $\Phi$ and adding the process noise $e$. The process noise is assumed to be normally distributed with a mean of 0 and a covariance matrix $W$. The reason for using discharge anomalies instead of discharge is that the process model setup in this study relies on the auto- and cross-covariance matrices of the state vector, which requires the expectation of the states to be zero as a prerequisite (Tourian, Schwatke, & Sneeuw, 2017), as detailed in Appendix B. Additionally, the process model consistently predicts the states at the same reaches on the following day from $r_{t-1}$ to $r_t$, that is, the assimilation cycle length is 1 day. The specific details of the state vector will be presented later in Section 3.3 (see Equation 10), following the introduction of the observation model, as it accommodates both the process and the observation model setups. In this section, we focus on the definition of the transition matrix and its process noise covariance matrix.
By incorporating the prerequisite that the state vector has an expected value of zero, the transition matrix $\Phi$ and its process noise covariance matrix $W$ can be derived by minimizing the process noise $e$ (Gelb & Corporation, 1974; Hofmann-Wellenhof & Moritz, 2005; Kurtenbach, 2020), as detailed in Appendix B, and are expressed as

$$ \Phi = \Sigma_{\Delta} \Sigma^{-1} \qquad (2) $$

and

$$ W = \Sigma - \Sigma_{\Delta} \Sigma^{-1} \Sigma_{\Delta}^{\top}, \qquad (3) $$

where $\Sigma$ represents the auto-covariance matrix of the state vector $r_t$, and $\Sigma_{\Delta}$ represents the cross-covariance matrix between the state vectors $r_{t-1}$ and $r_t$, both of which can be interpreted separately as (Johnson & Wichern, 2007)

$$ \begin{aligned} \Sigma &= C\{r_t\} = S_t P_{t,t} S_t \\ \Sigma_{\Delta} &= C\{r_{t-1}, r_t\} = S_{t-1} P_{t-1,t} S_t \end{aligned} \qquad (4) $$

where $S_{t-1}$ and $S_t$ are diagonal matrices containing the standard deviations of the state vectors $r_{t-1}$ and $r_t$, respectively. $P_{t,t}$ and $P_{t-1,t}$ denote the auto- and cross-correlation matrices of the state vectors, with correlation coefficients computed by a physically based spatiotemporal discharge correlation model developed by Paiva et al. (2015) (see Section S1 in Supporting Information S1). This model describes how discharge at a given river location relates to discharge at a downstream point after a time lag. The model output is determined by several hydrodynamic parameters, including diffusive wave celerity, hydraulic diffusivity, and temporal decorrelation length parameter, together with the spatial and temporal offsets.

Moreover, by substituting the definitions of the auto- and cross-covariance matrices (see Equation 4) into the expression for the transition matrix (see Equation 2), the transition matrix can be reformulated as follows:

$$ \Phi = S_{t-1} P_{t-1,t} P_{t,t}^{-1} S_t^{-1}. \qquad (5) $$

It is evident that the establishment of the transition matrix can be achieved by performing basic matrix operations using the standard deviation and correlation matrices, significantly enhancing computational efficiency.

## 3.2. Observation Model

The discharge estimation can be improved by combining the predicted state from the process model with the measured state in a discrete spatiotemporal field while considering their respective uncertainties. In our setup, there are two linear observation equations within the observation model. One is the rearranged discretized mass conservation equation, and the other directly incorporates the SWOT discharge outputs as observations.

The first observation equation is obtained from the discretized form of the mass conservation equation (Durand et al., 2014), which expresses the linear relationship among consecutive reaches of a river (see discretization Section S2 in Supporting Information S1). As a result, for a given reach $s$, the mass conservation equation can be written as follows:

$$ \Delta_s \begin{bmatrix} Q(s-1, t-\Delta t) \\ Q(s, t-\Delta t) \\ Q(s+1, t-\Delta t) \\ Q(s-1, t) \\ Q(s, t) \\ Q(s+1, t) \end{bmatrix} + \frac{\delta A(s, t) - \delta A(s, t-\Delta t)}{\Delta t} = 0 \qquad (6) $$
where the first term on the left side is the flow imbalance, representing the reach average of the spatial derivative of discharge, and the second term is the mass imbalance, representing the temporal derivative of the change in cross-sectional area. $\Delta_s$ is a matrix composed of reach-length-dependent discretization coefficients (see Equation S11 in Supporting Information S1), and $f$ is a vector containing discharge values at time steps $t$ – $\Delta t$ and $t$ for the reaches adjacent to reach $s$.

In our case, the mass imbalance term serves as the observation for the first observation equation, where the change in cross-sectional area is one of the specific outputs derived from SWOT measurements, denoted as $\delta A_{Sat}$. If we tentatively use $f$ to build the state vector $r_t$, that is, $r_t = f - \bar{f}$ (see Equation 1), we can establish the observation equation by reformulating the discretized mass conservation equation as follows:

$$z_1 = \mathbf{H}_1 r_t + v_1, \qquad (7)$$

with

$$z_1 = -\frac{\delta A_{Sat}(s,t) - \delta A_{Sat}(s,t - \Delta t)}{\Delta t} - \mathbf{H}_1 \bar{f} \\ \mathbf{H}_1 = \Delta_s \\ v_1 \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_1) \text{ with } \mathbf{R}_1 = \frac{\sigma_{\delta A_{Sat}(s,t)}^2 + \sigma_{\delta A_{Sat}(s,t - \Delta t)}^2}{\Delta t^2}$$

where the observation vector $z_1$ contains the mass imbalance term, along with an extra term $\mathbf{H}_1 \bar{f}$ to balance the average term $\bar{f}$ in the state vector definition. The observation matrix $\mathbf{H}_1$ comprises reach-length-dependent discretization coefficients, derived from reach boundaries in the perturbed data sets (see Table 3). The observation noise $v_1$ is assumed to be zero-mean Gaussian white noise, with its covariance matrix $\mathbf{R}_1$ obtained from the error propagation of the uncertainty in $\delta A_{Sat}$, denoted by $\sigma_{\delta A_{Sat}}$. It should be noted that the formal state vector has a higher dimensionality than $f$, and the observation vector $z_1$ typically consists of multiple observations, necessitating a corresponding expansion of the observation matrix $\mathbf{H}_1$ and the observation covariance matrix $\mathbf{R}_1$ (see Equations C4 and C5).

For the second observation equation, we ingest SWOT discharge outputs themselves as the observations. Then the linear observation equation reads

$$z_2 = \mathbf{H}_2 r_t + v_2 \quad v_2 \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_2) \qquad (8)$$

where the observation vector $z_2$ incorporates all available SWOT discharge measurements $Q_{Sat}$ during this assimilation window with the subtraction of their respective long-term averages and the observation matrix $\mathbf{H}_2$ is a selection matrix. Accordingly, the observation covariance matrix $\mathbf{R}_2$ is a diagonal matrix with the error variances of SWOT discharge measurements, denoted as $\sigma_{Q_{Sat}}^2$, on the diagonal.

In this study, the SWOT-measured change in cross-sectional area $\delta A_{Sat}$ and discharge $Q_{Sat}$ are defined as the medians of 10,000 Monte Carlo realizations of the perturbed measurements ($\delta A_{ptb}$ and $Q_{ptb}$), respectively (see Table 3). Their corresponding uncertainties, $\sigma_{\delta A_{Sat}}$ and $\sigma_{Q_{Sat}}$, are determined based on the interquartile range (IQR) of the realizations:

$$\sigma_{\delta A_{Sat}(s,t)} \approx \text{IQR}[\delta A_{ptb}=1,...,10\,000(s,t)]/1.349 \\ \sigma_{Q_{Sat}(s,t)} \approx \text{IQR}[Q_{ptb}=1,...,10\,000(s,t)]/1.349, \qquad (9)$$

where IQR is a measure of statistical dispersion, calculated as the difference between the 25th and 75th percentiles of the data. In a normal distribution, dividing the IQR by 1.349 provides a robust estimate of the standard deviation, effectively accounting for all measurement realizations while mitigating the influence of outliers (Fox, 2015). Moreover, as stated in Section 2, our created SWOT measurements conform to the SWOT sampling, under the support of pass plans generated by the CNES Large Scale SWOT Hydrology Simulator.
### 3.3. Kalman Filtering

With the process model (see Equation 1) and observation equations (see Equations 7 and 8) established, we have effectively outlined the components of a linear dynamic system. The next critical step involves solving this system to estimate daily discharge. To achieve this, we employ a Kalman filter, which recursively integrates the model-predicted state with the observation data to produce a weighted average estimate. The construction of a classical Kalman filter is grounded on two fundamental assumptions (Kalman, 1960): first, that the dynamic systems are linear; and second, that all errors, including those from the process model and observations, follow a normal distribution with a zero mean, indicating unbiased errors. Our model configuration sufficiently satisfies both assumptions, allowing us to formally construct our algorithm.

We formally establish the state vector $\mathbf{r}_t$ as follows:

$$ \mathbf{r}_t = \mathbf{q}_t - \bar{\mathbf{q}} = \underbrace{\begin{bmatrix} Q(1, t-21) \\ \vdots \\ Q(N, t-21) \\ \vdots \\ Q(1, t-20) \\ \vdots \\ Q(N, t-20) \\ \vdots \\ Q(1, t) \\ \vdots \\ Q(N, t) \end{bmatrix}}_{\mathbf{q}_t} - \underbrace{\begin{bmatrix} \bar{Q}(1) \\ \vdots \\ \bar{Q}(N) \\ \vdots \\ \bar{Q}(1) \\ \vdots \\ \bar{Q}(N) \\ \vdots \\ \bar{Q}(1) \\ \vdots \\ \bar{Q}(N) \end{bmatrix}}_{\bar{\mathbf{q}}} \quad (10) $$

$N$: Total number of river reaches

which consists of discharge states $\mathbf{q}_t$ across all reaches (from reach 1 to $N$) of a river for 22 days, with each value subtracted by its respective long-term mean daily discharge $\bar{\mathbf{q}}$. The 22-day duration was chosen to ensure that two distinct SWOT-measured changes in cross-sectional area are available at the same reach on different days, as required by the first observation equation (see Equation 7). Given the 21-day SWOT repeat orbit, the 22-day interval allows for at least two measurements of the SWOT-measured change in cross-sectional area at the same reach. With the setup of the state vector, the established dynamic system should be adapted accordingly (see Appendix C).

To implement the Kalman filter, Table 4 summarizes its general steps. As outlined in Section 1, the Kalman filter includes two main steps: prediction and update. The prediction step uses the process model to transit the state temporally, obtaining a predicted (a priori) state $\mathbf{r}_t^-$, the corresponding state covariance $\Sigma_t^-$ is derived by propagating the state covariance from the previous time step, that is, $\Sigma_{t-1}^+$, via the process model and adding it to the covariance of process noise $\mathbf{W}$.

The update step uses the observation model to refine the predicted state with the observation $\mathbf{z}_t$ once it is available at a specific time step. The improved state is referred to as the updated (a posteriori) state $\mathbf{r}_t^+$. It is computed by multiplying the innovation $\mathbf{y}_t$, that is, the residual between the current prediction and the current observation information, with the optimal Kalman gain $\mathbf{K}_t$, and then integrating this result with the predicted state. The optimal Kalman gain represents the weighting assigned to the observation and the current predicted state when the state errors are minimized. One of the terms to calculate $\mathbf{K}_t$ is the innovation covariance $\mathbf{V}_t$, which is derived by the error propagation of the innovation equation. The updated state covariance $\Sigma_t^+$ is obtained by multiplying the predicted state covariance $\Sigma_t^-$ with the term $(\mathbf{I} - \mathbf{K}_t\mathbf{H}_t)$, which reduces the uncertainty of the state estimate by incorporating the influence of new observations. It is noteworthy that the update step is not necessary for every scheduled time point, and in case the observation is unavailable for any reason, the update can be skipped and multiple prediction steps can be performed until a new observation becomes accessible.
<table><caption>**Table 4**
Kalman Filter Procedure</caption><thead><tr><th>Prediction</th><th>Equation</th></tr></thead><tbody><tr><td>Predicted state</td><td>r<sub>t</sub><sup>-</sup> = Φ<sub>t</sub>r<sub>t-1</sub><sup>+</sup></td></tr><tr><td>Predicted state covariance</td><td>Σ<sub>t</sub><sup>-</sup> = ΦΣ<sub>t-1</sub><sup>+</sup>Φ<sup>T</sup> + W</td></tr><tr><td>Update</td><td>Equation</td></tr><tr><td>Observation innovation</td><td>y<sub>t</sub> = z<sub>t</sub> - H<sub>t</sub>r<sub>t</sub><sup>-</sup></td></tr><tr><td>Innovation covariance</td><td>V<sub>t</sub> = H<sub>t</sub>Σ<sub>t</sub><sup>-</sup>H<sub>t</sub><sup>T</sup> + R<sub>t</sub></td></tr><tr><td>Optimal Kalman gain</td><td>K<sub>t</sub> = Σ<sub>t</sub><sup>-</sup>H<sub>t</sub><sup>T</sup>V<sub>t</sub><sup>-1</sup></td></tr><tr><td>Updated state</td><td>r<sub>t</sub><sup>+</sup> = r<sub>t</sub><sup>-</sup> + K<sub>t</sub>y<sub>t</sub></td></tr><tr><td>Updated state covariance</td><td>Σ<sub>t</sub><sup>+</sup> = (I - K<sub>t</sub>H<sub>t</sub>)Σ<sub>t</sub><sup>-</sup></td></tr></tbody></table>

Note. The prediction step estimates the prior state and its covariance based on the previous updated state. The update step incorporates the new observation to correct the prediction.

Given that the assimilation cycle length is 2 day and the state vector includes 22-day discharge across all reaches of a river, running the Kalman filter can typically yield up to 22 daily discharge estimates for each reach at a single epoch day. We thus need to decide how to combine these data. For this purpose, we present the Kalman filter output diagram on the left side of Figure 4, where the estimated discharge after each filtering cycle day is arranged in columns. One alternative for the combination is to extract the nth ($n = 1, ..., 22$) set of discharge estimates from each cycle day and aggregate them into an ensemble. Another alternative is to summarize all available discharge estimates for each epoch $t$ across cycle days (in the row direction) and calculate their weighted mean, arithmetic mean, and median, respectively. Figure 4 also illustrates both alternatives.

## 3.4. Implementation

To implement the methodology, we first set up the linear dynamic system. Within the process model, to compute the transition matrix and its process noise covariance matrix, we begin by constructing the auto- and cross-

correlation matrices $P_{t,t}$ and $P_{t-1,t}$ (see Equation C2) by calculating their correlation coefficients using the spatiotemporal correlation model (see Equation S2 in Supporting Information S1), which relies on the input variables, that is, diffusive wave celerity, hydraulic diffusivity, and temporal decorrelation length parameter, together with the spatial and temporal offsets, derived from data in perturbed data sets (see Appendix C). Next, we calculate the standard deviations of the true discharge for each reach over time series as the elements in the standard deviation matrices $S_t$ and $S_{t-1}$ (see Equation C1), and further compute the auto- and cross-covariance matrices $\Sigma$ and $\Sigma_\Delta$ (see Equation 4). Finally, the transition matrix $\Phi$ and its process noise covariance matrix $W$ are obtained using Equations 3 and 5, respectively.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Aggregation of the estimated discharge. The left side displays the Kalman filter output diagram. One way to combine the estimates is to extract the nth ($n = 1, ..., 22$) set of discharge estimates from each cycle day in columns and aggregate them as an ensemble, denoted by $Q_{est(n)}$. Another way is to calculate the weighted mean, arithmetic mean, and median of the discharge sets for each epoch $t$ (in the row direction) and summarize them into ensembles, denoted by $Q_{est(wm)}$, $Q_{est(am)}$, and $Q_{est(med)}$, as depicted in the diagram on the right side.</figcaption></figure>
The transition matrix and its process noise covariance matrix exhibit consistent structural patterns across test cases. As a representative example, Figures C1 and C2 in Appendix C show these matrices for the test case Missouri Downstream, respectively. Both matrices characterize how the discharge anomaly state vector $r_t$ (see Equation 10) evolves from $r_{t-1}$ within the process model (see Equation 1). In the visual layout of the transition matrix, the row and column indices of the matrices represent the elements of $r_t$ and $r_{t-1}$, respectively, consistent with the state prediction operation $r_t = \Phi r_{t-1}$. In contrast, the process noise covariance matrix has both axes corresponding to $r_t$, as it represents the covariance structure of the process noise affecting $r_t$ during the prediction step.

In the transition matrix in Figure C1, the red elements with value 1 indicate that discharge anomaly states from days $t-21$ to $t-1$ in $r_{t-1}$ are directly projected to the same temporal positions in $r_t$. These entries imply no change in the states across time, corresponding to identity mappings in the transition matrix. In contrast, the discharge anomaly states on day $t$, corresponding to the bottom rows of the transition matrix, are estimated through a linear combination of multiple past states from $r_{t-1}$. This is evident in the zoomed-in region at the bottom of Figure C1, where no single entry dominates, and the values are generally dispersed across several columns.

In the process noise covariance matrix in Figure C2, elevated noise levels are concentrated in the block related to states on day $t$ in $r_t$, visible as a distinct yellow region in the bottom-right corner (see the zoomed-in area). This increased noise arises because these states are predicted through a linear combination of states in $r_{t-1}$, in contrast to the direct projection applied for days $t-21$ to $t-1$, which results in near-zero noise elsewhere. According to the Kalman filter procedure (see Table 4), this pattern increases the predicted state uncertainty for day $t$, leading to a greater weight on observations during the update step.

Moreover, Figure C3 in Appendix C provides a visual representation of the matrices used to construct the transition matrix, as defined in Equation 5, for the same test case. We observe that both the auto- and cross-correlation matrices exhibit high correlation values, indicating strong correlations within $r_t$ and between $r_{t-1}$ and $r_t$. Notably, the auto-correlation matrix is highly sensitive to matrix inversion, suggesting the importance of carefully selecting the input variables related to river characteristics, that is, diffusive wave celerity, hydraulic diffusivity, and temporal decorrelation length parameter, when computing the correlation coefficients via the physically based model. Otherwise, inappropriate input choices may lead to a poorly conditioned transition matrix and, subsequently, significant errors in discharge estimation.

The two observation equations are established based on the number of available perturbed measurements within a filtering cycle. Hence, the dimensions and composition of the observation matrices (see Equations C5 and C8) are dynamically adjusted to accommodate changes in the observation vectors (see Equations C3 and C6), ensuring consistency with the corresponding measured states. Furthermore, the observation covariance matrices (see Equations C4 and C7) are derived from the uncertainties of the perturbed measurements (see Equation 9), respectively.

To initiate the Kalman filtering, we employ the first 22 days of perturbed discharge values at all reaches, along with their associated uncertainties derived from Equation 9, to construct the initial state vector and its state covariance matrix. In our approach, the mean daily discharge term in the state vector is calculated by averaging the SWOT-measured discharge over the available time series. However, this is not an ideal option, as the SWOT-measured discharge is not in daily sampling. Alternatively, one could utilize the daily true discharge or estimates from a hydrologic model, for example, the WBMsed-model (Cohen et al., 2014). Nevertheless, given the biases in the SWOT-measured discharge as described in Section 2, using the bias-free true discharge would offer an overly optimistic scenario.

After executing the Kalman filter and aggregating its results, we obtain a series of estimated discharge, including the 1st to 22nd sets of estimates, as well as the weighted mean, arithmetic mean, and median estimated discharge. To summarize, Figure C4 in Appendix C provides a flowchart illustrating the detailed implementation of the methodology using perturbed data sets. Moreover, an annotated example code is provided to facilitate reproducibility (Ke et al., 2025).

### 3.5. Validation Method

We employ the following validation metrics to evaluate the performance of the estimated daily discharge $Q_{est}$ against the true daily discharge $Q_{true}$ provided by the perturbed data sets (see Table 3):
* Correlation coefficient

* Nash-Sutcliffe Efficiency (NSE)

$$NSE = 1 - \frac{\sum_{t=1}^{T} [\epsilon_{\text{true}}(t) - \epsilon_{\text{est}}(t)]^2}{\sum_{t=1}^{T} [\epsilon_{\text{true}}(t) - \bar{\epsilon}_{\text{true}}]^2} \quad (11)$$

where

$$\epsilon_{\text{est}}(t) = |Q_{\text{est}}(t) - \bar{Q}_{\text{est}}| \\ \epsilon_{\text{true}}(t) = |Q_{\text{true}}(t) - \bar{Q}_{\text{true}}|$$

NSE ranges from -infinity to 1. An efficiency of 1 corresponds to a perfect match between the estimated and true discharge, whereas an efficiency of less than 0 suggests that the true mean discharge is a better predictor than our model.

* Normalized Root Mean Squared Error (nRMSE)

$$\text{nRMSE}[\%] = \frac{\text{RMSE}}{(\bar{\epsilon}_{\text{est}} + \bar{\epsilon}_{\text{true}})/2} \times 100\% = \frac{\sqrt{\frac{1}{T} \sum_{t=1}^{T} [\epsilon_{\text{est}}(t) - \epsilon_{\text{true}}(t)]^2}}{\frac{1}{2}(\frac{1}{T} \sum_{t=1}^{T} \epsilon_{\text{est}}(t) + \frac{1}{T} \sum_{t=1}^{T} \epsilon_{\text{true}}(t))} \times 100\% \quad (12)$$

* Relative bias (rBias)

$$\text{rBias}[\%] = \frac{|\bar{Q}_{\text{est}} - \bar{Q}_{\text{true}}|}{\bar{Q}_{\text{true}}} \times 100\% \quad (13)$$

* Error

$$\text{Error} = |Q_{\text{est}} - Q_{\text{true}}| \quad (14)$$

* Relative error

$$\text{Relative error}[\%] = \frac{\text{Error}}{\bar{Q}_{\text{true}}} \times 100\% \quad (15)$$

In calculating NSE and nRMSE, we use the residuals of the estimated discharge and the true discharge, denoted as $\epsilon_{est}$ and $\epsilon_{true}$, to exclude the influence of bias, particularly the average term in the state vector. This approach ensures that both metrics assess only the variability of the estimates relative to the reference. As complementary measures to evaluate the overall discrepancy between the estimated and true discharge, the metrics rBias, error, and relative error are introduced.

Serving as a critical reference for evaluating the estimated discharge, we also validate the SWOT-measured discharge against the accordingly sampled true discharge using the above metrics. It is important to note that, although the measured discharge is typically considered the true value in the context of validation, in this study, it is identified as data to refine the model predictions within the Kalman filter framework rather than as a benchmark. If the performance of the estimated discharge is close to or even better than that of the measured discharge, it indicates that our linear dynamic system is not only capable of enhancing the temporal resolution of discharge estimation but also effectively captures the evolution of discharge time series well. As an additional reference, we implement a simple linear interpolation of the discrete measured discharge at each reach to generate interpolated
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Estimated daily discharge time series ($Q_{est(med)}$, red curves) for the reaches in Kushiyara, compared to the SWOT-measured discharge ($Q_{Sat}$, blue asterisks) and the true discharge ($Q_{true}$, black curves).</figcaption></figure>

daily discharge, denoted as $Q_{interp}$, which is also validated against the daily true discharge employing the aforementioned metrics to evaluate whether our method demonstrates advantages over this simple interpolation method.

# 4. Results and Validation

In the following, for clarity and intuitiveness, we select only the median estimated discharge, denoted by $Q_{est(med)}$, to present results (hereafter referred to simply as *the estimated discharge*), as it leverages discharge estimates from each cycle day and consistently demonstrates relatively reliable performance compared to the other sets. For completeness, a full comparison of all estimated discharge sets is provided in Supporting Information S1, including hydrographs (Section S3 in Supporting Information S1) and corresponding performance metrics (Section S4 in Supporting Information S1) across all test cases. As an additional test, we also use the true daily discharge to build the mean daily discharge term in the state vector, as mentioned in Section 3.4, resulting in an overly optimistic estimation (see Section S5 in Supporting Information S1).

We proceed by first inspecting and analyzing the hydrographs in Section 4.1, then presenting an overall validation using various performance metrics in Section 4.2, and finally assessing our method's performance under extreme conditions through the Iowa River test case in Section 4.3.

## 4.1. Hydrograph Analysis

For the sake of better readability and discussion, we select two typical test cases for demonstration: Kushiyara, without bad reaches, and Ohio Section 1, with bad reaches (see Table 2).

Figure 5 shows the hydrographs of estimated daily discharge for reaches in Kushiyara compared to the SWOT-measured discharge and the true discharge. Overall, the estimated curves capture the main trends of the true discharge, demonstrating the method's effectiveness in resolving daily variations. Nevertheless, there is a substantial deviation between the estimates and true values, which stems from the construction of the state vector. In our approach, the state vector is composed of discharge anomalies (see Equation 10), defined as deviations from the mean of the SWOT-measured discharge in the framework. However, due to significant bias in the measurements, the computed mean itself is overestimated, effectively raising the baseline for all anomaly based
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Estimated daily discharge time series ($Q_{est(med)}$, red curves) for the reaches in Ohio Section 1 against the SWOT-measured discharge ($Q_{Sat}$, blue asterisks) and the true discharge ($Q_{true}$, black curves). Reach 6, framed by a pink square, is a bad reach without observations.</figcaption></figure>

estimates. The perturbed discharge utilized to create the SWOT-measured values has a median relative bias as high as 353%, according to Table A1. These biases primarily arise from biased flow law parameters within MetroMan (Durand et al., 2014) and the layover bias in interferometric SAR of SWOT (Durand et al., 2020). If we use the true discharge to calculate the average term (see Figure S52 in Supporting Information S1), the estimates then lie at the same level as the true discharge.

Figure 6 shows the estimated daily discharge in Ohio Section 1, compared against both the SWOT-measured and the true discharge. Compared to Kushiyara (see Figure 5), the SWOT-measured discharge here exhibits less bias, with the corresponding perturbed discharge showing a median relative bias of 76% (see Table A1), resulting in a closer match between the estimates and the true values. Notably, reach 6 is a bad reach with a lack of measurements. Despite this, the method still produces a reasonable approximation of the discharge, underscoring its potential to improve the spatial resolution of discharge estimation, especially in unobserved areas.

## 4.2. Overall Validation

As an initial evaluation of the estimated discharge performance, we calculate its median NSE values (see Equation 11) across all reaches in each test case and present them on a spatial map in Figure 7. Overall, more than half of the rivers achieve NSE values above 0.6, with the Ohio River basin exhibiting notably strong performance, underscoring the ability of our approach to provide accurate discharge estimates. However, certain rivers, such as the Iowa River and the Padma, exhibit considerably lower NSE values, necessitating further investigation.

To further assess the estimation performance, we compile the validation results, including correlation, NSE, nRMSE (see Equation 12), and rBias (see Equation 13), into bar charts, as shown in Figure 8. For each test case, the median value of each metric is calculated across all reaches. In addition to the estimates, the figure also presents the performance of the SWOT-measured and interpolated discharge, both evaluated against the true discharge and serving as important references. Moreover, test cases with bad reaches are highlighted with pink frames to facilitate subsequent analysis. To complement the bar charts, Table 5 summarizes the exact metric
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spatial map of the median NSE values for the estimated discharge. Three zoomed-in panels, that is, North America, Western Europe, and South Asia, highlight the regional details, with their locations shown in the global map below. Each numbered marker denotes a test case (listed on the right), with the color indicating the median NSE across all river reaches within that case. The classification intervals reflect different levels of discharge estimation performance, as defined in the legend.</figcaption></figure>

values for each test case, while Figure 9 presents scatter plots comparing the metrics for the three types of discharge.

As shown in Figures 8 and 9 and detailed in Table 5, the estimated discharge generally shows high correlation with the true discharge, although slightly lower than that of the SWOT-measured discharge. In most cases, it still outperforms or closely matches the interpolated discharge. In terms of NSE and nRMSE, the estimated discharge performs comparably to, or even better than, the SWOT-measured discharge, while clearly outperforming the interpolated discharge across most test cases. Regarding rBias, it typically exhibits smaller deviations from the true discharge compared to both references, with pronounced exceptions observed in Brahmaputra and Ohio Section 8. Notably, four test cases contain bad reaches where SWOT measurements are unavailable. Despite this, the estimated discharge still demonstrates satisfactory performance.

In addition to the aggregated metrics, the cumulative distribution functions (CDFs) of relative errors (see Equation 15) for each test case are plotted in Figure 10 to illustrate the distribution of estimation errors across all reaches. These curves offer a complementary perspective by highlighting how the magnitude and spread of errors differ among the estimated discharge, the SWOT-measured discharge, and the interpolated discharge.

Across the majority of test cases, the CDF curves of the estimated discharge tend to be positioned to the left of those corresponding to the SWOT-measured and interpolated discharge, indicating a higher proportion of small relative errors. This trend suggests that the proposed approach not only improves the spatiotemporal resolution of the discharge series but also generally preserves or enhances estimation accuracy. In particular, it can provide estimates for bad reaches, where linear interpolation is not applicable. However, cases such as Kushiyara and Missouri Downstream exhibit mild S-shaped patterns across all three CDF curves, indicating the presence of systematic errors, which are also evident in their corresponding time series plots (see Figure 5 and Figure S6 in Supporting Information S1). Furthermore, the Jamuna, Kushiyara, and Padma display heavier-tailed distributions, reflecting a substantial portion of large relative errors. These features are consistent with the elevated rBias shown in the corresponding bar plots (see Figure 8). Collectively, the above information reinforces the notion that both
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Median of the performance metrics across all reaches for each test case in bar charts. Bars in red represent the validation results for the estimated discharge *Q*<sub>est(med)</sub>, while blue and yellow bars correspond to the SWOT-measured *Q*<sub>Sat</sub> and linearly interpolated daily discharge *Q*<sub>interp</sub>, respectively, serving as crucial references for assessing estimation quality. Test cases with bad reaches are framed by a pink frame.</figcaption></figure>

systematic and large-magnitude errors contribute to the estimation challenges in these basins, likely due to the propagation of bias and uncertainty from the SWOT-measured discharge within the Kalman filter framework.

Notably, in cases such as Ohio Sections 2 and 5, interpolation outperforms the SWOT-measured discharge over a range of percentiles. This can occur when the SWOT measurements capture the overall trend of the true discharge, and interpolation, through its smoothing effect, leverages this trend to generate more stable daily estimates.

### 4.3. Evaluation at Extremes

Furthermore, our method demonstrates robustness under flood events. The Iowa River, which features a 500-year return period flood (Gilles et al., 2012), provides a particularly relevant example. As shown in the hydrograph plot (see Figure 11), the true discharge exhibits sharp peaks exceeding 1,000 m³/s in several reaches, in contrast to near-zero flows during normal conditions. Despite the poor performance of the SWOT-measured discharge, reflected by a median NSE as low as -1.4, a median nRMSE of 111%, and a median rBias of 280%, our estimated discharge is able to maintain, and even slightly improve, these metrics, with corresponding values of -1.11%, 110%, and 261% (see Table 5). As shown in the CDF plot (see Figure 10), the relative errors of the estimated discharge lie to the left of those of the SWOT-measured discharge across most percentiles, indicating improved accuracy in most reaches. This improvement is also evident in the hydrograph plot. For example, at reach 8, the estimated discharge mitigates the impact of large errors present in the SWOT-measured discharge. More
<table><caption>**Table 5**
Median of the Performance Metrics Across All Reaches for Each Test Case</caption><thead><tr><th rowspan="2">Test case name</th><th colspan="3">Correlation</th><th colspan="3">NSE</th><th colspan="3">nRMSE (%)</th><th colspan="3">rBias (%)</th></tr><tr><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th></tr></thead><tbody><tr><td>Brahmaputra</td><td>0.81</td><td>0.99</td><td>0.78</td><td>0.33</td><td>0.64</td><td>0.35</td><td>221</td><td>104</td><td>213</td><td>86</td><td>62</td><td>62</td></tr><tr><td>Iowa River</td><td>0.68</td><td>0.76</td><td>0.69</td><td>-1.11</td><td>-1.40</td><td>-1.25</td><td>110</td><td>111</td><td>115</td><td>261</td><td>280</td><td>281</td></tr><tr><td>Jamuna</td><td>-0.34</td><td>0.22</td><td>0.16</td><td>-0.97</td><td>-0.08</td><td>-0.10</td><td>217</td><td>159</td><td>171</td><td>485</td><td>506</td><td>508</td></tr><tr><td>Kushiyara</td><td>0.83</td><td>0.82</td><td>0.81</td><td>0.61</td><td>-0.75</td><td>-0.54</td><td>73</td><td>107</td><td>98</td><td>368</td><td>398</td><td>399</td></tr><tr><td>Missouri Downstream</td><td>0.76</td><td>0.91</td><td>0.82</td><td>0.58</td><td>0.81</td><td>0.60</td><td>108</td><td>75</td><td>107</td><td>88</td><td>92</td><td>92</td></tr><tr><td>Missouri Upstream</td><td>0.78</td><td>0.95</td><td>0.88</td><td>0.56</td><td>0.72</td><td>0.65</td><td>135</td><td>132</td><td>125</td><td>25</td><td>31</td><td>31</td></tr><tr><td>Ohio Section 1</td><td>0.88</td><td>0.99</td><td>0.54</td><td>0.77</td><td>0.89</td><td>0.00</td><td>60</td><td>41</td><td>133</td><td>38</td><td>46</td><td>45</td></tr><tr><td>Ohio Section 2</td><td>0.92</td><td>0.99</td><td>0.63</td><td>0.76</td><td>0.81</td><td>0.14</td><td>53</td><td>45</td><td>112</td><td>53</td><td>56</td><td>56</td></tr><tr><td>Ohio Section 3</td><td>0.91</td><td>0.99</td><td>0.89</td><td>0.81</td><td>0.91</td><td>0.73</td><td>61</td><td>41</td><td>69</td><td>15</td><td>11</td><td>12</td></tr><tr><td>Ohio Section 4</td><td>0.86</td><td>0.98</td><td>0.82</td><td>0.61</td><td>0.70</td><td>0.59</td><td>101</td><td>86</td><td>103</td><td>5</td><td>13</td><td>14</td></tr><tr><td>Ohio Section 5</td><td>0.93</td><td>0.99</td><td>0.92</td><td>0.78</td><td>0.87</td><td>0.75</td><td>69</td><td>55</td><td>75</td><td>25</td><td>29</td><td>29</td></tr><tr><td>Ohio Section 7</td><td>0.95</td><td>0.99</td><td>0.96</td><td>0.80</td><td>0.89</td><td>0.83</td><td>65</td><td>48</td><td>59</td><td>19</td><td>20</td><td>20</td></tr><tr><td>Ohio Section 8</td><td>0.93</td><td>0.98</td><td>0.92</td><td>0.55</td><td>0.63</td><td>0.53</td><td>120</td><td>103</td><td>120</td><td>72</td><td>69</td><td>69</td></tr><tr><td>Padma</td><td>0.63</td><td>0.63</td><td>0.64</td><td>-1.96</td><td>-7.20</td><td>-5.69</td><td>148</td><td>168</td><td>167</td><td>130</td><td>177</td><td>177</td></tr><tr><td>Seine Downstream</td><td>0.93</td><td>1.00</td><td>0.97</td><td>0.80</td><td>0.96</td><td>0.88</td><td>65</td><td>28</td><td>50</td><td>12</td><td>12</td><td>12</td></tr></tbody></table>

Note. This table complements the bar charts in Figure 8 by providing the exact values of correlation, NSE, nRMSE, and rBias for the estimated $Q_{est(med)}$, the SWOT-measured $Q_{Sat}$, and the interpolated discharge $Q_{interp}$, respectively.

importantly, while reaches 2 and 4 are bad reaches, our method still produces relatively stable estimates there, further demonstrating its robustness and ability to fill data gaps effectively.

From the evaluation above, it can be concluded that the proposed method effectively enhances the spatiotemporal resolution of discharge estimates while maintaining accuracy levels comparable to those of the SWOT-measured discharge. The method even demonstrates robustness under more dynamic flow conditions, such as flood events. In addition, it offers clear advantages over linear interpolation, particularly in terms of consistency and applicability to unobserved reaches. Although the estimated discharge does not consistently outperform the SWOT-measured discharge in terms of correlation, NSE, and nRMSE, the discrepancies remain within an acceptable range. This is reasonable given the higher temporal resolution of the daily estimates compared to the discrete SWOT measurements.

## 5. Discussion

Although the overall validation confirms the effectiveness of the proposed method across a wide range of test cases, the estimates still reveal underperformance in certain rivers. To better understand the underlying causes, we use Brahmaputra as a representative example and present its hydrographs (see Figure 12), performance metrics in both bar chart and tabular form (see Figure 13 and Table 6), and CDF plots of relative errors (see Figure 14) at the reach level.

As shown, while the estimated discharge captures the overall trends across most reaches in the Brahmaputra, its accuracy generally remains limited. Notably, reach 5 exhibits substantial challenges for discharge retrieval, characterized by systematic overestimation in the hydrograph, low correlation, negative NSE, and large rBias. The S-shaped patterns in the CDF plots of relative errors at this reach further underscore the presence of systematic bias in both the estimated and SWOT-measured discharge, consistent with the hydrograph plot.

Several key factors may explain the underperformance of the discharge estimates:

* Measurement quality and frequency. The accuracy of the estimated discharge depends heavily on the quality and frequency of SWOT measurements. At reach 5, the SWOT-measured discharge is significantly biased, which leads to an inflated mean used in the construction of the discharge anomaly state vector. As a result, the
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of the performance metrics between the estimated discharge *Q*<sub>est(med)</sub>, the SWOT-measured discharge *Q*<sub>Sat</sub>, and the linearly interpolated discharge *Q*<sub>interp</sub> across all test cases, as a complement to the bar charts in Figure 8. Each subplot shows a scatter plot for a specific metric: (a) correlation, (b) NSE, (c) nRMSE, and (d) rBias. In each plot, the x-axis shows the metric value for the SWOT-measured or interpolated discharge, with markers depicted as circles and crosses, respectively, while the y-axis shows the corresponding metric for the estimates. The dashed line represents the 1:1 reference line.</figcaption></figure>

baseline of the estimates is shifted upward, producing persistent overestimation (see Figure 12). On the other hand, infrequent SWOT data reduces the dynamic system's ability to reflect short-term flow fluctuations. For instance, around day 30 at reaches 1 to 4, the estimated discharge fails to capture the rapid rise in true discharge due to the lack of timely SWOT measurements incorporated into the Kalman filter process (see Figure 12).

* *Morphology-induced parameter error.* The wide, braided channels and low banks of the Brahmaputra (Pareta, 2024) complicate SWOT observations and affect the reliability of measured variables, such as river width and slope. As a result, several key input parameters for the physically based discharge correlation model, such as diffusive wave celerity and hydraulic diffusivity, are derived from these SWOT values and thus suffer from reduced accuracy. These inaccuracies can distort the correlation structure and transition matrix, thereby degrading the performance of the Kalman filtering process and ultimately the discharge estimates. Notably, although we use perturbed data sets in this study, the error characteristics follow parameterized models based on the SWOT high-fidelity simulator, which also capture morphology-induced uncertainties (Frasson et al., 2021).

* *Mismatched discharge variability assumption.* There are strong seasonal contrasts at Brahmaputra reaches across the Brahmaputra reaches. Discharge variability is high during the monsoon but remains relatively flat during other periods (see Figure 12). Applying a uniform variability assumption over time, as done in the current implementation when computing the standard deviation matrices in the transition matrix, can result in inaccurate state predictions in both high- and low-flow regimes.

* *Correlation distortion at confluences.* Our approach neglects the geomorphologic dispersion due to the superposition of flows at river confluences, which can affect the spatiotemporal dependence of discharge (Paiva et al., 2015). In the Brahmaputra (see Figure 2), the presence of tributaries disrupts downstream correlation patterns. Ignoring this effect leads to overestimated discharge correlation and reduced estimation accuracy.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Cumulative distribution function (CDF) plots of the relative errors for the estimated discharge ($Q_{est(med)}$, in red), the SWOT-measured discharge ($Q_{Sat}$, in blue), and the interpolated discharge ($Q_{interp}$, in yellow), evaluated against the true discharge across all reaches in each test case. Test cases containing bad reaches are highlighted with pink boxes.</figcaption></figure>

Based on the findings discussed above, several directions for improving the proposed approach can be identified. First, integrating complementary information from existing altimetry or imagery data sets may help mitigate the impact of biased SWOT measurements and improve the robustness of discharge estimation. Second, the observation covariance structure could be refined by accounting for spatiotemporal dependence across measurements, particularly to better handle sparse SWOT data coverage. Third, adopting seasonally varying discharge variability, such as using separate standard deviation representations for high-flow and low-flow periods, could improve the flexibility of the state transition formulation. Finally, incorporating the influence of tributaries into the correlation structure may help to better capture discharge dynamics at confluences (Paiva et al., 2015).

**6. Summary and Conclusions**

In this study, we propose a method for estimating daily discharge using SWOT data by characterizing the spatiotemporal dynamics of continuous reaches in a single-branch river network. Our model-free assimilation method is based on building a linear dynamic system that assimilates multiple SWOT outputs along a river. The linear dynamic system consists of a stochastic process model and two observation equations. The stochastic process model comprises the auto- and cross-covariance information of river discharge along the river, employing discharge correlation computed from a physically derived spatiotemporal correlation model and discharge standard deviation based on prior river information. One observation equation is the discretized mass conservation equation, using the change in cross-sectional area, one of the SWOT outputs, as observations. The other observation equation directly assimilates SWOT discharge. The linear dynamic system is effectively solved using a Kalman filter, which is executed in the time domain to obtain a sequential discharge for consecutive reaches.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Estimated daily discharge time series ($Q_{est(med)}$, red curves) for the reaches in the Iowa River, compared to the SWOT-measured discharge ($Q_{Sat}$, blue asterisks) and the true discharge ($Q_{true}$, black curves). Test cases containing bad reaches are highlighted with pink boxes.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Estimated daily discharge time series ($Q_{est(med)}$, red curves) for the reaches in Brahmaputra, compared to the SWOT-measured discharge ($Q_{Sat}$, blue asterisks) and the true discharge ($Q_{true}$, black curves).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Bar charts of the performance metrics for individual reaches in Brahmaputra. Each subplot shows the validation results for the estimated discharge ($Q_{est(med)}$, in red) alongside the SWOT-measured ($Q_{Sat}$, in blue) and interpolated discharge ($Q_{interp}$, in yellow), evaluated against the true values across five reaches. At most reaches, the estimates exhibit inferior performance relative to the references.</figcaption></figure>

Since SWOT data were unavailable during the period of this study, particularly as the discharge product is scheduled for release around mid-2025, we employ a perturbed version of synthetic SWOT data sets obtained by Monte Carlo simulation to verify the feasibility of our approach. The data sets contain 15 river test cases with a range of 4–14 reaches, where some of them have biased measurements. It should be noted that the perturbed data sets do not contain observations in reaches with low-head dams and other obstacles, which are referred to as bad reaches.

After the Kalman filtering, we obtain the daily discharge time series at each reach (including bad reaches) without data outages. Due to the configuration of the state vector and prediction scheme, we can obtain a maximum of 22 discharge realizations per day at each reach. Thus, we build 25 combination ensembles for the outcomes: 22

<table><caption>**Table 6**

Performance Metrics for Individual Reaches in Brahmaputra</caption><thead><tr><th rowspan="2">Reach No.</th><th colspan="3">Correlation</th><th colspan="3">NSE</th><th colspan="3">nRMSE (%)</th><th colspan="3">rBias (%)</th></tr><tr><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th><th>Q<sub>est(med)</sub></th><th>Q<sub>Sat</sub></th><th>Q<sub>interp</sub></th></tr></thead><tbody><tr><td>1</td><td>0.80</td><td>0.99</td><td>0.93</td><td>0.29</td><td>0.48</td><td>0.35</td><td>242</td><td>169</td><td>225</td><td>93</td><td>71</td><td>70</td></tr><tr><td>2</td><td>0.83</td><td>0.99</td><td>0.96</td><td>0.33</td><td>0.77</td><td>0.54</td><td>229</td><td>104</td><td>176</td><td>86</td><td>62</td><td>62</td></tr><tr><td>3</td><td>0.81</td><td>0.99</td><td>0.78</td><td>0.33</td><td>0.84</td><td>0.42</td><td>221</td><td>62</td><td>182</td><td>70</td><td>49</td><td>49</td></tr><tr><td>4</td><td>0.88</td><td>0.98</td><td>0.61</td><td>0.42</td><td>0.64</td><td>0.22</td><td>195</td><td>83</td><td>235</td><td>73</td><td>61</td><td>61</td></tr><tr><td>5</td><td>0.54</td><td>-0.31</td><td>-0.27</td><td>0.29</td><td>-5.47</td><td>-1.63</td><td>175</td><td>204</td><td>213</td><td>732</td><td>785</td><td>785</td></tr></tbody></table>

Note. This table complements Figure 13 by providing exact values of correlation, NSE, nRMSE, and rBias for the estimated discharge $Q_{est(med)}$, the SWOT-measured discharge $Q_{Sat}$, and the interpolated discharge $Q_{interp}$, respectively.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** CDF plots of the relative errors in Brahmaputra for the estimated discharge ($Q_{est(med)}$, in red), the SWOT-measured discharge ($Q_{Sat}$, in blue), and the interpolated discharge ($Q_{interp}$, in yellow), evaluated against the true discharge at individual reaches. The curves exhibit an S-shaped pattern at reach 5, which reflects the presence of systematic bias.</figcaption></figure>

discharge sets extracted directly from realizations, as well as the weighted mean, the arithmetic mean, and the median discharge sets. As a simple reference approach, we also calculate a linearly interpolated daily discharge time series. Then the estimated ensembles, the measured, and the interpolated daily discharge are separately validated against true discharge to assess the effectiveness of our dynamic system.

Generally, the median validation values across reaches for the median discharge sets indicate a correlation of -0.34 to 0.95, a NSE for residuals (mean-removed discharge) of -1.96 to 0.81, a nRMSE for residuals of 53%-221%, and a rBias of 5%-485%, while the measured discharge shows the median correlation ranging 0.22-0.99, the median NSE for residuals ranging between -7.20 and 0.96, the median nRMSE for residuals ranging 28%-168%, and the median rBias ranging 11%-506%, and for the interpolated daily discharge, they are 0.16-0.97, -5.69-0.88, 50%-213%, and 12%-508%, respectively. We also summarize the relative errors of the discharge across all reaches into one CDF plot for each test case. The estimated discharge generally shows smaller errors than the interpolated daily discharge and follows the measured discharge closely. By comparing the estimated and the measured discharge, we can conclude that our estimates are close to the measured discharge in terms of correlation, NSE, nRMSE, and relative errors, and typically even outperform the measured values in terms of rBias. It goes without saying that our estimated discharge outperforms the simple interpolation method. Especially, it can provide estimates for unobserved reaches, which is not possible by linear interpolation.

Validation shows that our method can successfully estimate river discharge at a daily resolution for the 15 test cases and reveals great potential for utilizing SWOT products. In particular, although observation data is
unavailable at bad reaches and significantly biased at many reaches, our estimates can nevertheless demonstrate remarkably good performance. Inevitably, the performance of the estimates depends strongly on the quality and density of the observation data, as well as a careful characterization of the transition matrix.

Compared to previous efforts, this approach provides a model-free data assimilation framework that enhances SWOT discharge products by improving temporal resolution to a daily scale and spatial resolution to a reach scale. It leverages the classical Kalman filter, one of the simplest forms of data assimilation method, ensuring computational efficiency. In addition, the method is simple to implement and comprehend, without involving any complex hydrological models, while allowing direct ingestion of SWOT measurements into the observation model.

In future work, we will investigate the applicability of our method using SWOT products in larger river basins, even globally. The algorithm presented here was only tested on a single-branch river network with a limited number of reaches. We aim to further improve the method to estimate discharge across entire river networks by considering tributaries and the spatiotemporal dependence of discharge at confluences. Additionally, an enhancement to our assimilation method can be achieved by taking the lateral inflow into account via a modification of the observation covariance within the discretized mass conservation equation. Furthermore, our algorithm can benefit from existing altimetry and imagery data as complementary sources to SWOT observations to deliver assimilated river discharge estimates at the full catchment scale.

## Appendix A: Perturbed Measurements

At the cross-section scale, the test cases include simulated measurements of height $H$ and width $W$ for the entire area, and these fine-scale measurements are aggregated into reach-averaged measurements of $H$ and $W$, along with estimations of slope $S$. In the original synthetic data sets, the measurement set accounts for daily sampling and introduces errors in the cross-sectional $H$ and $W$. The reach-averaged $H$, $W$, and $S$ are then computed from these cross-sectional measurements while considering their error corruption. The error in cross-sectional $H$ consists of layover and random errors depending on $W$, river orientation, and topographic roughness, and the error of $W$ includes mainly the random error depending on the cross-section position in the SWOT swath (Frasson et al., 2021). The perturbed version, based on the original data sets, includes measurement error models consistent with the original for $H$ and $W$, and is extended by adding the capability to generate the change in cross-sectional area ($\delta A$) containing errors. $\delta A$, one of the SWOT outputs as well as an input for our observation model, represents the anomaly between the total cross-sectional area ($A$) and the baseline cross-sectional area ($A_0$), with $A_0$ serving as one of the flow law parameters in discharge estimation.

Moreover, in the perturbed data sets, 10,000 realizations of the daily reach-averaged measurements are generated through a Monte Carlo simulation, which entails generating multiple realizations of observation errors for each of the measurements utilized in the SWOT discharge ($Q_{ptb}$) production: $H$, $W$, $S$, and $\delta A$ (Frasson et al., 2023). For the computation of $Q_{ptb}$, the flow law parameter inversion algorithm Metropolis-Manning (MetroMan) is employed (Durand et al., 2014), which is one of the SWOT discharge inversion algorithms and is also conceived as an optimization algorithm, aiming to minimize the error in mass conservation. It uses the modified Gauckler-Manning-Strickler (Durand et al., 2023) equation along with the mass conservation equation and a Metropolis algorithm (a Bayesian sampling approach) to derive posterior distributions of the flow law parameters relying on input measurements of $H$, $W$, $S$, and $\delta A$, as well as a prior estimate of the mean annual discharge at the location of interest, derived from the water balance model WBMsed in our case, denoted by $Q_{WBM}$ (Cohen et al., 2014; Frasson et al., 2023). As a constraint, MetroMan must be applied within sets of at least three consecutive reaches that contain no tributaries. In this study, the MetroMan configuration assumes the discharge of a reach $s$ at time step $t$ can be determined using the modified Gauckler-Manning-Strickler equation in this version:

$$Q_{ptb}(s, t) = \frac{1}{n_a(s) \left[ \frac{A_0(s) + \delta A_{ptb}(s, t)}{W_{ptb}(s, t)} \right]^{b(s)}} [A_0(s) + \delta A_{ptb}(s, t)]^{5/3} W_{ptb}(s, t)^{-2/3} \sqrt{S_{ptb}(s, t)} \quad (A1)$$

where $W_{ptb}$, $\delta A_{ptb}$, and $S_{ptb}$ are the perturbed measurements. The baseline cross-sectional area $A_0$, the friction coefficient $n_a$, and the exponent coefficient $b$ denote the time-invariant perturbed flow law parameters estimated by MetroMan.
The detailed procedure for generating reach-averaged perturbed measurements $H_{ptb}$, $S_{ptb}$, $W_{ptb}$, and $\delta A_{ptb}$ is as follows:

* $H_{ptb}$: As the starting point for generating perturbed measurements, the height error model initially produces cross-section-scale (denoted by c) realizations at each time step t:

$$H_c(t) = H_T + \epsilon_H, \epsilon_H \sim \mathcal{N}(H_L, \sigma_H) \quad (A2)$$

where the measured cross-sectional height $H_c(t)$ is the summation of the true cross-sectional height $H_T$ and a measurement error $\epsilon_H$, which is generated through a normal distribution with the layover bias $H_L$ as the mean value and $\sigma_H$ as the standard deviation. Both $H_L$ and $\sigma_H$ depend on time- and space-varying river width, river orientation relative to the SWOT swath, and topographic standard deviation, as specified by Durand et al. (2020). The reach-averaged height $H_{ptb}(t)$ is then calculated by fitting a line to $H_c(t)$ within a reach. The midpoint of this line is designated as $H_{ptb}(t)$.

* $S_{ptb}$: The reach-averaged slope $S_{ptb}(t)$ is also derived from the fitted line to $H_c(t)$.

* $W_{ptb}$: Similarly to the height measurements, the production of width also starts at the cross-section scale c at each time step t:

$$W_c(t) = W_T + \epsilon_W, \epsilon_W \sim \mathcal{N}(0, \sigma_W) \quad (A3)$$

where $W_{ptb}(t)$ is computed as the true cross-sectional width $W_T$ perturbed by $\epsilon_W$, which is sampled from a normal distribution with a zero mean and a standard deviation $\sigma_W$. $\sigma_W$ is influenced by the position of the cross-section inside the swath and is determined by a function fitted to several simulated SWOT overpasses in the Po (Domeneghetti et al., 2017) and Sacramento Rivers (Frasson et al., 2017). Subsequently, the reach-averaged width $W_{ptb}(t)$ is obtained by averaging $W_c(t)$ within the reach boundaries.

* $\delta A_{ptb}$: First, the cross-sectional width-height observation relationships within each reach are constructed and fitted with a linear polynomial. Then, the reach-averaged $\delta A_{ptb}$ is calculated by the integral from the reference height $H_0$ to the height at any given time. This approach is similar to the calculation method in RiverObs, the SWOT processing algorithm (SWOT Algorithm Definition Team, 2021), but with the key difference that instead of a normal linear regression, RiverObs applies a piecewise error-model-based linear fit (Fuller, 1987) with two breakpoints. Therefore, the method employed in this study provides a comparatively conservative way to calculate errors of $\delta A_{ptb}$.

To identify the sources of bias in the perturbed discharge $Q_{ptb}$, Figure A1 shows the realizations of flow law parameters in histograms for the 3rd reach of the test case Missouri Upstream (see Table 2), along with the mean

<figure><img src="image_17.png" /><figcaption>Figure A1. 100 realizations of flow law parameters $n_a$, $A_0$, and b in histograms for the 3rd reach of Missouri Upstream (see Table 2). The differences between the red lines (mean value) and dashed black lines (ensemble mean) present the ensemble range of parameter errors, which are sampled from the posterior distributions of the Bayesian estimates (MetroMan), with no bias correction.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure A2.** 10,000 realizations of measurements $Q_{ptb}$, $W_{ptb}$, $H_{ptb}$, $S_{ptb}$, and $\delta A_{ptb}$ in histograms for the 3rd reach of the Missouri Upstream (see Table 2) at the 8th time step. Discrepancies between ensemble means and true values in the histograms of $Q_{ptb}$, $H_{ptb}$, $S_{ptb}$, and $\delta A_{ptb}$ reveal the existence of bias.</figcaption></figure>

(dashed black lines) and true values (red lines) to indicate the extent of errors. The differences between them highlight biases in the Bayesian posterior estimates within MetroMan. Additionally, Figure A2 shows the realizations of measurements in histograms for the same reach at the 8th time step. The histograms exhibit normal distributions, with $Q_{ptb}$, $H_{ptb}$, $S_{ptb}$, and $\delta A_{ptb}$ showing significant discrepancies between the mean and true values. The bias in $H_{ptb}$ arises from layover effects. For $S_{ptb}$ and $\delta A_{ptb}$, the biases result from error propagation due to the biased $H_{ptb}$. The bias in $Q_{ptb}$ is attributed to its derivation using the modified Gauckler-Manning-Strickler equation (see Equation A1), which relies on biased inputs including $A_0$, $n_a$, $b$, $S_{ptb}$, and $\delta A_{ptb}$.

Furthermore, Table A1 summarizes the median relative bias of $\delta A_{ptb}$ and $Q_{ptb}$ across all reaches for each of the 15 test cases. As shown, the biases vary considerably across rivers.

<caption>**Table A1**</caption>

Median Relative Bias (rBias) of Perturbed Change in Cross-Sectional Area $\delta A_{ptb}$ and Perturbed Discharge $Q_{ptb}$ Across All Reaches for Each Test Case (rBias[%] = $\frac{|\bar{X}_{ptb} - \bar{X}_{true}|}{\bar{X}_{true}} \times 100\%$, Where X Denotes the $\delta A$ or $Q$ Time Series)

<table><thead><tr><th>No.</th><th>Test case name</th><th>Median rBias of δA<sub>ptb</sub> (%)</th><th>Median rBias of Q<sub>ptb</sub> (%)</th></tr></thead><tbody><tr><td>1</td><td>Brahmaputra</td><td>61</td><td>96</td></tr><tr><td>2</td><td>Iowa River</td><td>72</td><td>286</td></tr><tr><td>3</td><td>Jamuna</td><td>9</td><td>494</td></tr><tr><td>4</td><td>Kushiyara</td><td>14</td><td>353</td></tr><tr><td>5</td><td>Missouri Downstream</td><td>28</td><td>99</td></tr><tr><td>6</td><td>Missouri Upstream</td><td>40</td><td>49</td></tr><tr><td>7</td><td>Ohio Section 1</td><td>15</td><td>76</td></tr><tr><td>8</td><td>Ohio Section 2</td><td>9</td><td>79</td></tr><tr><td>9</td><td>Ohio Section 3</td><td>9</td><td>39</td></tr><tr><td>10</td><td>Ohio Section 4</td><td>6</td><td>42</td></tr><tr><td>11</td><td>Ohio Section 5</td><td>6</td><td>47</td></tr><tr><td>12</td><td>Ohio Section 7</td><td>6</td><td>33</td></tr><tr><td>13</td><td>Ohio Section 8</td><td>3</td><td>73</td></tr><tr><td>14</td><td>Padma</td><td>19</td><td>147</td></tr><tr><td>15</td><td>Seine Downstream</td><td>22</td><td>38</td></tr></tbody></table>

Note. The strong correlation between them indicates that the bias in $Q_{ptb}$ is largely attributable to perturbations in $\delta A_{ptb}$.
## Appendix B: Derivation of Transition Matrix and Its Process Noise Covariance

In principle, the transition matrix $\mathbf{\Phi}$ of the process model (see Equation 1) should be estimated by minimizing the process error $e$ (Gelb & Corporation, 1974; Kurtenbach, 2020). The process error's dispersion can be quantified as (Tourian, Schwatke, & Sneeuw, 2017)

$$C\{e\} = E\{ee^{\top}\}, \qquad (B1)$$

and the minimum error variance then occurs when the trace of the error covariance matrix is minimum:

$$\mathrm{tr}[C\{e\}] \rightarrow \min. \qquad (B2)$$

The error covariance matrix of the process model can be extended as

$$\begin{aligned} C\{e\} &= E\{ee^{\top}\} = E\{(\mathbf{\Phi}r_{t-1} - r_t)(\mathbf{\Phi}r_{t-1} - r_t)^{\top}\} \\ &= E\{\mathbf{\Phi}r_{t-1}\mathbf{r}_{t-1}^{\top}\mathbf{\Phi}^{\top} - r_t\mathbf{r}_{t-1}^{\top}\mathbf{\Phi}^{\top} - \mathbf{\Phi}r_{t-1}\mathbf{r}_{t-1}^{\top} + r_t\mathbf{r}_{t-1}^{\top}\}. \end{aligned} \qquad (B3)$$

If the expectation of state vector $r_t$ equals $\mathbf{0}$:

$$E\{r_t\} = \mathbf{0}, \qquad (B4)$$

we can directly define the auto-covariance $\mathbf{\Sigma}$ of the state vector $r_t$ as the expectation of $r_t\mathbf{r}_{t-1}^{\top}$:

$$\mathbf{\Sigma} = C\{r_t\} = E\{r_t\mathbf{r}_{t-1}^{\top}\}, \qquad (B5)$$

and similarly the cross-covariance $\mathbf{\Sigma}_{\Delta}$ between $r_{t-1}$ and $r_t$ as

$$\mathbf{\Sigma}_{\Delta} = C\{r_{t-1}, r_t\} = E\{r_t\mathbf{r}_{t-1}^{\top}\}. \qquad (B6)$$

Then the error covariance matrix of the process model can be rewritten as

$$C\{e\} = E\{ee^{\top}\} = \mathbf{\Phi}\mathbf{\Sigma}\mathbf{\Phi}^{\top} - \mathbf{\Sigma}_{\Delta}\mathbf{\Phi}^{\top} - \mathbf{\Phi}\mathbf{\Sigma}_{\Delta}^{\top} + \mathbf{\Sigma}, \qquad (B7)$$

which can be reformed to

$$C\{e\} = \mathbf{\Sigma} - \mathbf{\Sigma}_{\Delta}\mathbf{\Sigma}^{-1}\mathbf{\Sigma}_{\Delta}^{\top} + (\mathbf{\Phi} - \mathbf{\Sigma}_{\Delta}\mathbf{\Sigma}^{-1})(\mathbf{\Sigma}(\mathbf{\Phi} - \mathbf{\Sigma}_{\Delta}\mathbf{\Sigma}^{-1})^{\top}). \qquad (B8)$$

Because in Equation B8 the first two terms are independent of $\mathbf{\Phi}$, then $\mathrm{tr}[C\{e\}]$ becomes minimum if the third term disappears, which leads to the definitions of transition matrix $\mathbf{\Phi}$ (see Equation 2) and its process noise $\mathbf{W}$ (see Equation 3).

## Appendix C: Setup of the Linear Dynamic System and Methodology Implementation

In the formal setup, the state vectors (see Equation 10) consist of discharge states from all reaches (reach 1 to N) of a river over 22 days with a subtraction of the long-term mean daily discharge for each reach, rather than a single state. Consequently, we should adjust the transition matrix $\mathbf{\Phi}$ and its corresponding covariance matrix of process error $\mathbf{W}$ for the process model, as well as expand the components of the two observation equations.

Before computing $\mathbf{\Phi}$, it is necessary to build the standard deviation and correlation matrices. The standard deviation matrices contain the uncertainties of discharge states on the diagonal:
$$S_t = \begin{bmatrix} \sigma_{Q(1,t-21)} & \mathbf{0} \\ \vdots & \vdots \\ \mathbf{0} & \sigma_{Q(N,t)} \end{bmatrix} \qquad (C1)$$

$$S_{t-1} = \begin{bmatrix} \sigma_{Q(1,t-22)} & \mathbf{0} \\ \vdots & \vdots \\ 0 & \sigma_{Q(N,t-1)} \end{bmatrix}$$

where $\sigma_Q$ are computed as the standard deviations of $Q_{true}$ for each reach over the time series. The correlation matrices should have a consistent structure with the state vectors and read

$$\mathbf{P}_{t,t} = \begin{bmatrix} \rho_{Q(1,t-21),Q(1,t-21)} & \cdots & \rho_{Q(1,t-21),Q(N,t)} \\ \vdots & \ddots & \vdots \\ \rho_{Q(N,t),Q(1,t-21)} & \cdots & \rho_{Q(N,t),Q(N,t)} \end{bmatrix} \qquad (C2)$$

$$\mathbf{P}_{t-1,t} = \begin{bmatrix} \rho_{Q(1,t-21),Q(1,t-22)} & \cdots & \rho_{Q(1,t-21),Q(N,t-1)} \\ \vdots & \ddots & \vdots \\ \rho_{Q(N,t),Q(1,t-22)} & \cdots & \rho_{Q(N,t),Q(N,t-1)} \end{bmatrix}$$

where the correlation elements are computed by using the analytical solution of the spatiotemporal correlation model (see Equation S2 in Supporting Information S1), which requires the inputs of distance $\Delta x$, time lag $\Delta t$, diffusive wave celerity $c_0$, hydraulic diffusivity $D_0$, and temporal decorrelation length $\tau$. For $\Delta x$, reach boundaries rch_bnd are already defined (see Table 3), thus the distances between reaches are easily accessible. $\Delta t$ depends on the period between states and can range from 0 to 21 days. The value of $c_0$ is obtained by fitting a power-law model $Q_{true} = i(A - A_0)^j$ to the variables total cross-sectional area A and $Q_{true}$, where the coefficients i and j and the baseline cross-sectional area $A_0$ are constants. Then the mean derivative (see Equation S3 in Supporting Information S1) of this power-law function within the data range at each reach is taken as $c_0$ (Paiva et al., 2015). $D_0$ is defined as the median of the values computed across all days using Equation S4 in Supporting Information S1, where we employ $W_{ptb}$ and $S_{ptb}$ in time series and a from WBMsed-model derived constant $Q_{WBM}$ (Cohen et al., 2014) for all reaches within a single test case (see Table 3). Note that if in a test case, there is a considerable disparity between $Q_{WBM}$ and the median of $Q_{true}$ over time series and reaches, we use the median of $Q_{true}$ instead of $Q_{WBM}$. $\tau$ is estimated by fitting Equation S5 in Supporting Information S1 to the empirical auto-correlation function using $Q_{true}$.

After forming the standard deviation and correlation matrices, we can compute the transition matrix $\Phi$ (see Equation 5) and afterward its corresponding process noise covariance matrix $\mathbf{W}$ (see Equation 3) by building the covariance matrices of the state vector (see Equation 4). Figures C1 and C2 show the matrices $\Phi$ and $\mathbf{W}$ for the test case Missouri Downstream, respectively. Additionally, Figure C3 visualizes the matrices used to construct $\Phi$, as defined in Equation 5, for the same test case.

Regarding the two observation equations, we assume that within the 22-day time range in the state vector, there are k measurements of $\delta A_{Sat}$ and $Q_{Sat}$, respectively. For the first observation equation, because the observation includes two sequential $\delta A_{Sat}$ measurements, the size of $z_1$ is determined by the combination formula: $m_1 = C_k^2 = \frac{k!}{2!(k-2)!}$, and the observation matrix $\mathbf{H}_1$ is thus a $m_1 \times (22 \times N)$-matrix. The observation vector $z_1$ can be expanded to
<figure><img src="image_19.png" /><figcaption>**Figure C1.** Transition matrix for Missouri Downstream, with matrix entries representing the state prediction from $r_{t-1}$ (columns) to $r_t$ (rows). This test case includes five reaches over 22 days, resulting in a matrix size of $(5 \times 22 \text{ days}) \times (5 \times 22 \text{ days}) = 110 \times 110$ for the transition matrix. The red diagonal elements with a value of 1 indicate that the discharge anomaly states from day $t-1$ to $t-21$ in $r_{t-1}$ are directly projected to the same temporal positions in $r_t$ during the transition. In contrast, the states on day $t$ in $r_t$, associated with the bottom rows of the matrix (see the zoomed-in part at the bottom), are computed as a linear combination of multiple states from $r_{t-1}$.</figcaption></figure>

$$z_1 = \begin{bmatrix} -\frac{\delta A_{\text{Sat}}(s_1, t_1) - \delta A_{\text{Sat}}(s_1, t_1 - \Delta t_1)}{\Delta t_1} \\ -\frac{\delta A_{\text{Sat}}(s_2, t_2) - \delta A_{\text{Sat}}(s_2, t_2 - \Delta t_2)}{\Delta t_2} \\ \vdots \\ -\frac{\delta A_{\text{Sat}}(s_{m_1}, t_{m_1}) - \delta A_{\text{Sat}}(s_{m_1}, t_{m_1} - \Delta t_{m_1})}{\Delta t_{m_1}} \end{bmatrix} - \mathbf{H}_1 \bar{q} \quad (C3)$$

where $z_1$ includes an extra term $\mathbf{H}_1 \bar{q}$, which aligns with the subtraction term $\mathbf{H}_1 \bar{f}$ in Equation 7. The corresponding observation covariance $\mathbf{R}_1$ is a diagonal matrix:

$$\mathbf{R}_1 = \begin{bmatrix} \frac{\sigma_{\delta A(s_1, t_1)}^2 + \sigma_{\delta A(s_1, t_1 - \Delta t_1)}^2}{\Delta t_1^2} & \mathbf{0} \\ \mathbf{0} & \frac{\sigma_{\delta A(s_{m_1}, t_{m_1})}^2 + \sigma_{\delta A(s_{m_1}, t_{m_1} - \Delta t_{m_1})}^2}{\Delta t_{m_1}^2} \end{bmatrix}. \quad (C4)$$

Additionally, $\mathbf{H}_1$ reads
<figure><img src="image_20.png" /><figcaption>**Figure C2.** Process noise covariance matrix for Missouri Downstream, with matrix entries representing the uncertainty associated with $r_t$ (rows and columns). It has the same dimensions as the transition matrix (110 × 110). Most areas exhibit values close to zero, except for the bottom-right area (see the zoomed-in part on the right), which corresponds to states on day t in $r_t$. According to the Kalman filter procedure (see Table 4), this pattern leads to a significant increase in the predicted state covariance for day t, resulting in a greater reliance on observations during the update step. The elevated noise arises from the fact that these states are estimated through a linear combination of multiple states in $r_{t-1}$, introducing greater uncertainty compared to the direct projections applied for days $t-21$ to $t-1$.</figcaption></figure>

$$ \mathbf{H}_1 = \begin{bmatrix} 0 & \cdots & \delta_{s_1} & \cdots & \delta_{s_1} & \cdots & 0 & \cdots & 0 & \cdots & 0 & \cdots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & \delta_{s_2} & \cdots & \delta_{s_2} & \cdots & 0 & \cdots & 0 & \cdots \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \ddots & \ddots & \vdots & \ddots & \vdots & \ddots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 0 & \cdots & 0 & \cdots & \delta_{s_{m_1}} & \cdots & \delta_{s_{m_1}} & \cdots \end{bmatrix} \quad (C5) $$

which is a sparse matrix, containing only the reach-length-dependent elements $\delta$ (see Equation S11 in Supporting Information S1) at the appropriate positions associated with the relevant state elements.

For the second observation equation, the dimension of the observation vector $z_2$ equals the number of $Q_{Sat}$, denoted by $m_2 = k$. Consequently, $\mathbf{H}_1$ is a $k \times (22 \times N)$-matrix. $z_2$ is interpreted as

$$ z_2 = \begin{bmatrix} Q_{Sat}(s_1, t_1) \\ Q_{Sat}(s_2, t_2) \\ \vdots \\ Q_{Sat}(s_{m_2}, t_{m_2}) \end{bmatrix} - \begin{bmatrix} \bar{Q}_{Sat}(s_1) \\ \bar{Q}_{Sat}(s_2) \\ \vdots \\ \bar{Q}_{Sat}(s_{m_2}) \end{bmatrix} \quad (C6) $$

which also involves a subtraction of the average values to balance the state vector. Correspondingly, $\mathbf{R}_2$ is written as
<figure><img src="image_21.png" /><figcaption>**Figure C3.** Computation of the transition matrix $\Phi_{s,t}$ in Missouri Downstream utilizing the correlation ($P_{t,t}$ and $P_{t-1,t}$) and standard deviation matrices ($S_t$ and $S_{t-1}$), as defined in Equation 5.</figcaption></figure>

$$ \mathbf{R}_2 = \begin{bmatrix} \sigma_{Q\text{Sat}}^2(s_1, t_1) & \cdots & \mathbf{0} \\ \mathbf{0} & \sigma_{Q\text{Sat}}^2(s_{m_2}, t_{m_2}) \end{bmatrix}. \qquad (C7) $$

The observation matrix $\mathbf{H}_2$ is expressed as

$$ \mathbf{H}_2 = \begin{bmatrix} 0 & \cdots & 1 & \cdots & 0 & \cdots \\ 0 & \cdots & 0 & \cdots & 1 & \cdots \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1 & \cdots \end{bmatrix} \qquad (C8) $$

which is a sparse matrix composed of elements 1 to match the measured states.

If we aggregate the two observation equations together, we have

$$ \underbrace{\begin{bmatrix} z_1 \\ z_2 \end{bmatrix}}_{z_t} = \underbrace{\begin{bmatrix} \mathbf{H}_1 \\ \mathbf{H}_2 \end{bmatrix}}_{\mathbf{H}_t} r_t + \underbrace{\begin{bmatrix} v_1 \\ v_2 \end{bmatrix}}_{v_t} \qquad (C9) $$

as the observation model in the implementation of the Kalman filter. Furthermore, the aggregated observation covariance $\mathbf{R}_t$ reads

$$ \mathbf{R}_t = \begin{bmatrix} \mathbf{R}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{R}_2 \end{bmatrix}. \qquad (C10) $$

Figure C4 presents a flowchart detailing the implementation of the entire methodology using perturbed data sets.
<figure><img src="image_22.png" /><figcaption>**Figure C4.** Flowchart illustrating the detailed implementation of the methodology, including all necessary inputs for constructing a linear dynamic system using perturbed data sets.</figcaption></figure>

## Data Availability Statement

The full data sets needed for this research can be freely downloaded from Zenodo (Ke et al., 2025) under the Creative Commons Attribution 4.0 International license.

## Acknowledgments

This study is carried out within the ERDSWOT project as part of the SWOT-science team for the period 2020-2024. Siqi Ke gratefully acknowledges financial support from China Scholarship Council. Omid Elmi is supported by the German Research Foundation (Deutsche Forschungsgemeinschaft, grant 324641997) and contributes to the Research Units 2630, GlobalCDA: understanding the global freshwater system by combining geodetic and remote sensing information with modeling using a calibration/data assimilation approach (www.globalcda.de). Renato Frasson acknowledges financial support from the SWOT Science Team Grant 80NSSC20K1141. A portion of this research was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). Pierre-Olivier Malaterre acknowledges support from the CNES-TOSCA HYDROS project. Open Access funding enabled and organized by Projekt DEAL.

## References

Adams, T., Chen, S., Davis, R., Schade, T., & Lee, D. (2010). The Ohio River community HEC-RAS model. In *World environmental and water resources Congress 2010* (pp. 1512–1523). https://doi.org/10.1061/41114(371)160

Alsdorf, D. E., Rodríguez, E., & Lettenmaier, D. P. (2007). Measuring surface water from space. *Reviews of Geophysics*, *45*(2). https://doi.org/10.1029/2006RG000197

Altenau, E. H., Pavelsky, T. M., Durand, M. T., Yang, X., Frasson, R. P. d. M., & 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. https://doi.org/10.1029/2021WR030054

Andreadis, K. M., Brinkerhoff, C. B., & Gleason, C. J. (2020). Constraining the assimilation of SWOT observations with hydraulic geometry relations. *Water Resources Research*, *56*(5), e2019WR026611. https://doi.org/10.1029/2019WR026611

Biancamaria, S., Lettenmaier, D. P., & Pavelsky, T. M. (2016). The SWOT mission and its capabilities for land hydrology. *Surveys in Geophysics*, *37*(2), 307–337. https://doi.org/10.1007/s10712-015-9346-y

Brakenridge, R., Cohen, S., Kettner, A., De Groeve, T., Syvitski, J., & Fekete, B. (2012). Calibration of satellite measurements of river discharge using a global hydrology model. *Journal of Hydrology*, *475*, 123–136. https://doi.org/10.1016/j.jhydrol.2012.09.035

Brakenridge, R, Nghiem, S. V., Anderson, E., & Mic, R. (2007). Orbital microwave measurement of river discharge and ice status. *Water Resources Research*, *43*(4). https://doi.org/10.1029/2006WR005238

Brinkerhoff, C. B., Gleason, C. J., Feng, D., & Lin, P. (2020). Constraining remote river discharge estimation using reach-scale geomorphology. *Water Resources Research*, *56*(11), e2020WR027949. https://doi.org/10.1029/2020WR027949

CNES, LEGOS, & Theia. (2013). HydroWeb [Web Platform]. Retrieved from https://hydroweb.next.theia-land.fr

Cohen, S., Kettner, A. J., & Syvitski, J. P. (2014). Global suspended sediment and water discharge dynamics between 1960 and 2010: Continental trends and intra-basin sensitivity. *Global and Planetary Change*, *115*, 44–58. https://doi.org/10.1016/j.gloplacha.2014.01.011

David, C. H., Maidment, D. R., Niu, G.-Y., Yang, Z.-L., Habets, F., & Eijkhout, V. (2011). River network routing on the NHDPlus dataset. *Journal of Hydrometeorology*, *12*(5), 913–934. https://doi.org/10.1175/2011JHM1345.1

Deltares, W. W. F., & World Resources Institute. (2023). Global Water Watch (GWW). [Web Platform]. Retrieved from https://globalwaterwatch. earth

DGFI-TUM. (2013). Database for hydrological time series of inland waters (DAHITI). [Web Platform]. Retrieved from https://dahiti.dgfi.tum. de/en/

Dingman, S. L., & Bjerklie, D. M. (2006). Estimation of river discharge. In *Encyclopedia of hydrological sciences* (Vol. 5, p. chap. 61). John Wiley & Sons, Ltd. https://doi.org/10.1002/0470848944.hsa069

Domeneghetti, A., Molari, G., Tourian, M. J., Tarpanelli, A., Behnia, S., Moramarco, T., et al. (2021). Testing the use of Single- and multi-mission satellite altimetry for the calibration of hydraulic models. *Advances in Water Resources*, *151*, 103887. https://doi.org/10.1016/j.advwatres. 2021.103887
Domeneghetti, A., Schumann, G. J. P., Wei, R., Frasson, R. P. M., Durand, M., Pavelsky, T., et al. (2017). Water surface elevation from the upcoming SWOT mission under different flows conditions. Geophysical Research Abstracts, 19, 6551. Retrieved from https://meetingorganizer.copernicus.org/EGU2017/EGU2017-6551.pdf

Durand, M., Chen, C., de Moraes Frasson, R. P., Pavelsky, T. M., Williams, B., Yang, X., & Fore, A. (2020). How will Radar layover impact SWOT measurements of water surface elevation and slope, and estimates of river discharge? Remote Sensing of Environment, 247, 111883. https://doi.org/10.1016/j.rse.2020.111883

Durand, M., Gleason, C. J., Garambois, P. A., Bjerklie, D., Smith, L. C., Roux, H., et al. (2016). An intercomparison of remote sensing river discharge estimation algorithms from measurements of river height, width, and slope. Water Resources Research, 52(6), 4527–4549. https://doi.org/10.1002/2015WR018434

Durand, M., Gleason, C. J., Pavelsky, T. M., Frasson, R. P. D. M., Turmon, M., David, C. H., et al. (2023). A framework for estimating global river discharge from the surface water and ocean topography satellite mission. Water Resources Research, 59(4), e2021WR031614. https://doi.org/10.1029/2021WR031614

Durand, M., Neal, J., Rodríguez, E., Andreadis, K. M., Smith, L. C., & Yoon, Y. (2014). Estimating reach-averaged discharge for the river severn from measurements of river water surface elevation and slope. Journal of Hydrology, 511, 92–104. https://doi.org/10.1016/j.jhydrol.2013.12.050

Elmer, N. J., Hain, C., Hossain, F., Desroches, D., & Pottier, C. (2020). Generating proxy SWOT water surface elevations using WRF-Hydro and the CNES SWOT Hydrology simulator. Water Resources Research, 56(8), e2020WR027464. https://doi.org/10.1029/2020WR027464

Elmi, O., Tourian, M. J., Bárlossy, A., & Sneeuw, N. (2021). Spaceborne river discharge from a nonparametric stochastic quantile mapping function. Water Resources Research, 57(12), e2021WR030277. https://doi.org/10.1029/2021WR030277

Elmi, O., Tourian, M. J., Saemian, P., & Sneeuw, N. (2024). Remote sensing-based extension of GRDC discharge time series—A monthly product with uncertainty estimates. Scientific Data, 11(1), 240. https://doi.org/10.1038/s41597-024-03078-6

Emery, C., David, C., Andreadis, K., Turmon, M., Reager, J., Hobbs, J., et al. (2020). Underlying fundamentals of Kalman filtering for river network modeling. Journal of Hydrometeorology, 21(3), 453–474. https://doi.org/10.1175/JHM-D-19-0084.1

Emery, C., Paris, A., Biancamaria, S., Boone, A., Calmant, S., Garambois, P.-A., & Da Silva, J. (2018). Large-scale hydrological model river storage and discharge correction using a satellite altimetry-based discharge product. Hydrology and Earth System Sciences, 22(4), 2135–2162. https://doi.org/10.5194/hess-22-2135-2018

Even, S., Mouchel, J.-M., Servais, P., Flipo, N., Poulin, M., Blanc, S., et al. (2007). Modelling the impacts of combined sewer overflows on the river seine water quality. Science of the Total Environment, 375(1), 140–151. https://doi.org/10.1016/j.scitotenv.2006.12.007

Fekete, B., & Vörösmarty, C. (2007). The current status of global river discharge monitoring and potential new technologies complementing traditional discharge measurements [proceedings of the PUB kick-off meeting held in Brasília, 20–22 November 2002]. IAHS Publ. Predictions in Ungauged Basins: Pub Kick-Off, 309, 129–136.

Fox, J. (2015). Applied regression analysis and generalized Linear models. In (p. 36). SAGE Publications. Retrieved from https://books.google.de/books?id=3wrwCQAAQBAJ

Frasson, R. P. D. M., Durand, M. T., Larnier, K., Gleason, C., Andreadis, K. M., Hagemann, M., et al. (2020). Synthetic river datasets built for testing and development of the surface water and ocean topography mission discharge algorithms [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.3941890

Frasson, R. P. D. M., Durand, M. T., Larnier, K., Gleason, C., Andreadis, K. M., Hagemann, M., et al. (2021). Exploring the factors controlling the error characteristics of the surface water and ocean topography mission discharge estimates. Water Resources Research, 57(6), e2020WR028519. https://doi.org/10.1029/2020WR028519

Frasson, R. P. D. M., Turmon, M. J., Durand, M. T., & David, C. H. (2023). Estimating the relative impact of measurement, parameter, and flow law errors on discharge from the surface water and ocean topography mission. Journal of Hydrometeorology, 24(3), 425–443. https://doi.org/10.1175/JHM-D-22-0078.1

Frasson, R. P. D. M., Wei, R., Durand, M., Minear, J., Domeneghetti, A., Schumann, G., et al. (2017). Automated river reach definition strategies: Applications for the surface water and ocean topography mission. Water Resources Research, 53(10), 8164–8186. https://doi.org/10.1002/2017WR020887

Fuller, W. (1987). Measurement error models. In (p. 474). John Wiley & Sons, Ltd.

Gejadze, I., Malaterre, P.-O., Oubanas, H., & Shutyavae, V. (2022). A new robust discharge estimation method applied in the context of SWOT satellite data processing. Journal of Hydrology, 610, 127909. https://doi.org/10.1016/j.jhydrol.2022.127909

Gelb, A., & Corporation, T. A. S. (1974). Applied optimal estimation. MIT Press.

Gilles, D., Young, N., Schroeder, H., Piotrowski, J., & Chang, Y.-J. (2012). Inundation mapping initiatives of the Iowa flood center: Statewide coverage and detailed urban flooding analysis. Water, 4(1), 85–106. https://doi.org/10.3390/w4010085

Gleason, C. J., & Durand, M. T. (2020). Remote sensing of river discharge: A review and a framing for the discipline. Remote Sensing, 12(7), 1107. https://doi.org/10.3390/rs12071107

Hagemann, M. W., Gleason, C. J., & Durand, M. T. (2017). Bam: Bayesian AMHG-manning inference of discharge using remotely sensed stream width, slope, and height. Water Resources Research, 53(11), 9692–9707. https://doi.org/10.1002/2017WR021626

Hofmann-Wellenhof, B., & Moritz, H. (2005). Physical geodesy. SpringerWienNewYork, (pp. 361–363). https://doi.org/10.1007/978-3-211-33545-1

Huang, Q., Long, D., Du, M., Zeng, C., Qiao, G., Li, X., et al. (2018). Discharge estimation in high-mountain regions with improved methods using multisource remote sensing: A case study of the upper Brahmaputra River. Remote Sensing of Environment, 219, 115–134. https://doi.org/10.1016/j.rse.2018.10.008

Ishitsuka, Y., Gleason, C. J., Hagemann, M. W., Beighley, E., Allen, G. H., Feng, D., et al. (2021). Combining optical remote sensing, McFLI discharge estimation, global hydrologic modeling, and data assimilation to improve daily discharge estimates across an entire large watershed. Water Resources Research, 57(3), e2020WR027794. https://doi.org/10.1029/2020WR027794

Johnson, R. A., & Wichern, D. W. (2007). Applied multivariate statistical analysis. (p. 139). Pearson Prentice Hall.

JPL Internal Document. (2020). SWOT product description: Level 2 KaRIN high rate river single pass vector product. [Document]. California Institute of Technology. Retrieved from https://podaac-tools.jpl.nasa.gov/drive/files/misc/web/misc/swot_mission_docs/pdd/D-56413_SWOT_Product Description_L2_HR_RiverSP_20200825a.pdf

Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35–45. https://doi.org/10.1115/1.3662552

Ke, S., Tourian, M. J., & Frasson, R. P. d. M. (2025). Perturbed synthetic SWOT datasets for testing and development of a Kalman filter approach to estimate daily discharge [Dataset]. Zenodo. https://doi.org/10.5281/ZENODO.15482735
Kim, Y., & Bang, H. (2018). Introduction to Kalman filter and its applications. In F. Govaers (Ed.), *Introduction and implementations of the Kalman filter* (chap. 2). *Rijeka: IntechOpen*. https://doi.org/10.5772/intechopen.80600

Kurtenbach, E. (2020). Entwicklung eines Kalman-Filters zur Bestimmung kurzzeitiger Variationen des Erdschwerefeldes aus Daten der Satellitenmission GRACE. (Ph.D. thesis). University of Bonn.

Larnier, K., Monnier, J., Garambois, P.-A., & Verley, J. (2020). River discharge and bathymetry estimation from SWOT altimetry measurements. *Inverse Problems in Science and Engineering*, 29(6), 1–31. https://doi.org/10.1080/17415977.2020.1803858

Lin, P., Pan, M., Beck, H. E., Yang, Y., Yamazaki, D., Frasson, R., et al. (2019). Global reconstruction of naturalized river flows at 2.94 million reaches. *Water Resources Research*, 55(8), 6499–6516. https://doi.org/10.1029/2019WR025287

Lorenz, C., Tourian, M. J., Devaraju, B., Sneeuw, N., & Kunstmann, H. (2015). Basin-scale runoff prediction: An ensemble Kalman filter framework based on global hydrometeorological data sets. *Water Resources Research*, 51(10), 8450–8475. https://doi.org/10.1002/2014WR016794

Maswood, M., & Hossain, F. (2016). Advancing river modelling in ungauged basins using satellite remote sensing: The case of the Ganges-Brahmaputra–Meghna Basin. *International Journal of River Basin Management*, 14(1), 103–117. https://doi.org/10.1080/15715124.2015.1089250

McCabe, M. F., Rodell, M., Alsdorf, D. E., Miralles, D. G., Uijlenhoet, R., Wagner, W., et al. (2017). The future of Earth observation in hydrology. *Hydrology and Earth System Sciences*, 21(7), 3879–3914. https://doi.org/10.5194/hess-21-3879-2017

Nathanson, M., Kean, J. W., Grabs, T. J., Seibert, J., Laudon, H., & Lyon, S. W. (2012). Modelling rating curves using remotely sensed LiDAR data. *Hydrological Processes*, 26(9), 1427–1434. https://doi.org/10.1002/hyp.9225

Neal, J., Schumann, G., Bates, P., Buytaert, W., Matgen, P., & Pappenberger, F. (2009). A data assimilation approach to discharge estimation from space. *Hydrological Processes*, 23(25), 3641–3649. https://doi.org/10.1002/hyp.7518

Oubanas, H., Gejadze, I., Malaterre, P.-O., Durand, M., Wei, R., Frasson, R. P. M., & Domeneghetti, A. (2018). Discharge estimation in ungauged basins through variational data assimilation: The potential of the SWOT mission. *Water Resources Research*, 54(3), 2405–2423. https://doi.org/10.1002/2017WR021735

Paiva, R. C. D., Collischonn, W., Bonnet, M.-P., de Gonçalves, L. G. G., Calmant, S., Getirana, A., & Santos da Silva, J. (2013). Assimilating in situ and Radar altimetry data into a large-scale hydrologic-hydrodynamic model for streamflow forecast in the amazon. *Hydrology and Earth System Sciences*, 17(7), 2929–2946. https://doi.org/10.5194/hess-17-2929-2013

Paiva, R. C. D., Durand, M. T., & Hossain, F. (2015). Spatiotemporal interpolation of discharge across a river network by using synthetic SWOT satellite data. *Water Resources Research*, 51(1), 430–449. https://doi.org/10.1002/2014WR015618

Pareta, K. (2024). Morphological model for erosion prediction of India’s largest braided River using MIKE 21C model. *Earth Science, Systems and Society*, 4(1), 10075. https://doi.org/10.3389/esss.2024.10075

Pavelsky, T. M., Durand, M. T., Andreadis, K. M., Beighley, R. E., Paiva, R. C., Allen, G. H., & Miller, Z. F. (2014). Assessing the potential global extent of SWOT River discharge observations. *Journal of Hydrology*, 519, 1516–1525. https://doi.org/10.1016/j.jhydrol.2014.08.044

Saemian, P., Elmi, O., Stroud, M., Riggs, R., Kitambo, B. M., Papa, F., et al. (2025). Satellite Altimetry-based Extension of global-scale in situ river discharge Measurements (SAEM). *Earth System Science Data*, 17(5), 2063–2085. https://doi.org/10.5194/essd-17-2063-2025

Smith, L. C., & Pavelsky, T. M. (2008). Estimation of river discharge, propagation speed, and hydraulic geometry from space: Lena river, Siberia. *Water Resources Research*, 44(3). https://doi.org/10.1029/2007WR006133

Sun, W., Ishidaira, H., Bastola, S., & Jingshan, Y. (2015). Estimating daily time series of streamflow using hydrological model calibrated based on satellite observations of river water surface width: Toward real world applications. *Environmental Research*, 139, 36–45. https://doi.org/10.1016/j.envres.2015.01.002

SWOT Algorithm Definition Team. (2021). RiverObs – Processing algorithm used to transform raw pixel clouds obtained by the surface water and ocean topography mission into river products [Software]. Retrieved from https://github.com/SWOTAlgorithms/RiverObs

Tourian, M. J., Elmi, O., Mohammadnejad, A., & Sneeuw, N. (2017a). Estimating river depth from SWOT-type observables obtained by satellite altimetry and imagery. *Water*, 9(10), 753. https://doi.org/10.3390/w9100753

Tourian, M. J., Elmi, O., Shafaghi, Y., Behnia, S., Saemian, P., Schlesinger, R., & Sneeuw, N. (2022). HydroSat: Geometric quantities of the global water cycle from geodetic satellites. *Earth System Science Data*, 14(5), 2463–2486. https://doi.org/10.5194/essd-14-2463-2022

Tourian, M. J., Schwatke, C., & Sneeuw, N. (2017b). River discharge estimation at daily resolution from satellite altimetry over an entire river Basin. *Journal of Hydrology*, 546, 230–247. https://doi.org/10.1016/j.jhydrol.2017.01.009

Tourian, M. J., Sneeuw, N., & Bárdossy, A. (2013). A quantile function approach to discharge Estimation from Satellite Altimetry (ENVISAT). *Water Resources Research*, 49, 1–13. https://doi.org/10.1002/wrcr.20348

USACE. (2018). Missouri river recovery management plan and environmental impact statement: Existing conditions, unsteady HEC-RAS model calibration report. (p. 52). Retrieved from https://usace/contentdm.oclc.org/digital/collection/p16021coll17/id/7886

Yoon, Y., Durand, M., Merry, C. J., Clark, E. A., Andreadis, K. M., & Alsdorf, D. E. (2012). Estimating river bathymetry from data assimilation of synthetic SWOT measurements. *Journal of Hydrology*, 464–465, 363–375. https://doi.org/10.1016/j.jhydrol.2012.07.028