Satellite Video Remote Sensing for Estimation of River Discharge

We demonstrate that river discharge can be estimated by deriving water surface velocity estimates from satellite‐derived video imagery when combined with high‐resolution topography of channel geometry. Large Scale Particle Image Velocimetry (LSPIV) was used to map surface velocity from 28 s of 5 Hz satellite video acquired at a 1.2 m nominal ground spacing over the Darling River, Tilpa, Australia, during a 1‐in‐5‐year flood. We stabilized and assessed the uncertainty of the residual motion induced by the satellite platform, enhancing our sub‐pixel motion analysis, and quantified the sensitivity of image extraction rates on computed velocities. In the absence of in situ observations, LSPIV velocity estimates were validated against predictions from a calibrated 2D hydrodynamic model. Despite the confounding influence of selecting a surface velocity depth‐averaging coefficient, inference of discharge was within 0.3%–15% compared with gauging station measurements. These results provide a valuable foundation for refining satellite video LSPIV techniques.


Introduction
Globally, 29% of the world's population is exposed to flood risk and insecure water supplies, yet knowledge of the river discharges upon which flood and water resource management depends remains inadequate (Rentschler & Salhab, 2020).Global monitoring networks for quantifying river discharge are in decline, gauging remains logistically difficult and there are political influences on data sharing (King et al., 2018;Lins, 2008;Zakharova et al., 2020).However, satellite-based remote sensing approaches to monitoring discharge are helping to alleviate these issues (e.g., Sichangi et al., 2016).
Approaches to satellite-based river discharge monitoring typically rely upon statistical and hydraulic approximations to make indirect estimates of river discharge.Widely applied satellite radar altimetry measures water elevations at virtual river cross-sections (Revel et al., 2023;Tarpanelli et al., 2013;Zakharova et al., 2020) and near-simultaneous optical imagery can be used to infer water surface flow velocity from space (Kääb et al., 2019).Other satellite approaches have relied on remote sensing of discharge (RSQ) algorithms, which retrieve hydraulic variables from remotely sensed data and then relate these quantities to river discharge (Q) (e.g., Gleason & Durand, 2020;Riggs et al., 2022).These techniques are limited by relatively coarse spatial resolution and the requirement for near-simultaneous satellite swath overlaps, constraining global coverage.
High resolution commercial satellite video sensors can record the dynamics of river flow and floods.Optical flow measurement algorithms can estimate velocity by tracking the movement of visible features between frames (e.g., Eltner et al., 2020;Perks et al., 2020).Currently, optical satellite video acquired by low earth orbiting sensors offer spatial resolutions (pixel sizes) ranging from 0.9-1.2m at frame rates up to 30 Hz [for example, SkySat (Bhushan et al., 2021) and Jilin-1 (European Space Agency, 2022) constellations].Inference of flow velocities using satellite video has previously been demonstrated by Legleiter and Kinzel (2021a), who used 17 frames of cloud-free satellite video acquired by Planet Labs SkySat constellation of the Tanana River in central Alaska.Surface flow velocities were estimated to within 8.65% of radar gauging measurements and were further assessed using asynchronously acquired acoustic Doppler current profiler (aDcp) velocity data.
We apply and test the use of satellite video-based velocities for estimating discharge.We couple freely available, high-resolution topographic data with velocity estimates derived from satellite video Large-Scale Image Velocimetry (LSPIV) (Lewis et al., 2018;Muste et al., 2008) and some critical assumptions regarding channel hydraulics to estimate flood discharge following monsoonal rainfall in Darling River at Tilpa, Australia.The accuracy of satellite video-derived velocity estimates was assessed via comparison to hydraulic model simulations; and discharge estimates were compared with in situ gauging station observations.

Study Area
The River Darling at Tilpa (Figure 1) is located within the 502,500 km 2 Murray-Darling basin (Matheson & Thoms, 2018;Murray-Darling Basin Authority, 2010).The basin has a strongly episodic climate, with large floods followed by lengthy dry spells due to the influence of the El Niño Southern Oscillation (Grimaldi et al., 2019).Extensive and prolonged rainfall from late February to early April 2022 led to a flood with a 5-year return period (Q = 722 m 3 s −1 ).This location is ideally suited to testing our ability to measure discharge using non-contact, image-based velocity calculation techniques due to the availability of: (a) cloud-free satellite video; (b) a high-resolution LiDAR digital elevation model (DEM) acquired during dry river bed conditions; and (c) gauged in situ discharge observations at Tilpa (Water NSW, 2023).

Satellite Video
Satellite video was acquired over our study area on 5 February 2022 at 23:12 UTC by a Jilin-1 GF-03 sensor, part of the constellation by the Chang Guang Satellite Technology Company.Video had a 1.22 m spatial resolution and native frame rate of 5 Hz for 28 s.To counter sensor platform movement and scene "morphing" due to the changing view angle of the satellite overpass, we stabilized the video using FIJI's TrakEM2 plugin (Cardona et al., 2012;FIJI-ImageJ, 2020;Schindelin et al., 2012).FIJI is an open-source image processing toolkit.TrakEM2 relies on a Scale Invariant Feature Transform (SIFT) algorithm to align image stacks based on common features.To avoid geometric distortions, and because all video frames had a similar resolution, we utilized an affine transform to register our image stacks.
Temporal displacement errors related to image stabilization can significantly influence the accuracy of LSPIV velocities.We quantified the temporal distribution of frame-by-frame residual motion, by evaluating the cumulative frame-by-frame displacement (d) of six manually selected ground control points (GCPs; Figure 1) in every frame of the stabilized frame sequence.The movement of these GCPs post-stabilization provided a clear picture of the residual motion.To analyze this residual motion, we employed the differential Root Mean Square Difference metric (d(RMSD)) (Ljubičić et al., 2021).The d(RMSD) metric quantified the magnitude of the residual displacement of static features, based on a pixel intensity RMSD.This RMSD metric operates by directly comparing several subregions within subsequent images.In each subregion, it calculates the differences in pixel intensities at corresponding locations between two images.These differences are then squared, summed over all pixels in the subregions, and averaged.The square root of this average provides the d(RMSD) value.This value quantified residual motion magnitude and aided in understanding the temporal distribution of residual displacements in the video.
Computation of surface flow velocities in PIVlab is attained by cross-correlation algorithms applied to orthorectified images recorded at a known time interval.We evaluate the accuracy of both Fast Fourier Transform window deformation (direct FFT correlation with multiple passes and deforming windows) and Ensemble correlation (Figure 2).Interrogation areas (IA), which are small windows of defined size (in pixels), are used to track the displacement of image patterns within a chosen larger search area (SA) in subsequent images.The multi-pass FFT window deformation approach allows for the spatial resolution of velocity measurements to be improved through multiple reductions in the size of the interrogation areas over which correlations are calculated.Ensemble correlation is better suited for sparsely seeded images as it relies on averaging correlation matrices followed Writing -review & editing: Christopher Masafu, Richard Williams, Martin D. Hurst by detecting a correlation peak with the resultant benefit of lower bias and displacement errors (Thielicke & Sonntag, 2021).Given the relatively coarse spatial resolution of satellite video frames for PIV, where inter-frame movement of features are small, often less than the width of a single pixel (e.g., Legleiter & Kinzel, 2021a), PIVlab's sub-pixel motion estimation functions allow for more accurate and reliable sub-pixel peak determination.PIVlab implements both 2.3-point and 9-point Gaussian functions to resolve sub-pixel displacements (see Thielicke, 2014 for detail) making sub-pixel motion estimation possible.
We focused on two cloud-free and straight river reaches A and B (Figure 1) to reduce computational cost.Image pre-processing was performed to amplify the visibility of surface tracers with respect to the background (riverbanks/static ground), applying a Contrast-limited adaptive histogram equalization (CLAHE) filter (with a window size of 8 pixels, matching our smallest IA size, see Section 3.2.1) to enhance image contrast (Li & Yan, 2022;Masafu et al., 2022).Pre-processing of images for LSPIV has a significant impact on the quality of flow velocity estimates.Although image enhancement techniques such as high pass filters and intensity capping (amongst a multitude of others) exist, CLAHE offered a balanced approach to image enhancement due to its ability to enhance local contrast and visibility of tracer particles without excessively amplifying noise.
Distinct features on the water surface were difficult to discern in the raw images, which would be expected in natural rivers observed from the height of the optical sensor.However, CLAHE contrast enhancement enabled the tracking of seeding surrogates in the image sequences, which occur when specular reflection formed by incident light interacts with free-surface deformations on the river.Image intensity variations associated with these surface deformations were visible in post-processed images.

Sensitivity to Image Frame Rate and PIV Algorithm
The primary free parameters in LSPIV are the sampling frequency (frame extraction rate), interrogation (IA) and search (SA) areas; optimal configurations vary significantly (Kim et al., 2008;Legleiter & Kinzel, 2020;Sharif, 2022).IA should be small enough to eliminate spurious velocities whilst being large enough to accommodate an adequate window for surface pattern tracking (Tauro et al., 2018;Zhu & Lipeme Kouyi, 2019).Sampling frequency (frame extraction rate) and the IA are closely coupled and must be considered in tandem, with frame-to-frame displacement influencing the accuracy of pattern/particle detection on images.
FFT window deformation and Ensemble correlation algorithms were utilized with the maximum allowable number of PIV algorithm passes within PIVlab (four) for our sensitivity analysis (see Zhu & Lipeme Kouyi, 2019).We processed images using an IA of 64 × 64 pixels, with successive passes of 32 × 32, 16 × 16, and 8 × 8 pixels, all with 50% overlaps, corresponding to a minimum spatial distance of 9.8 m.SA sizes for our analyses were 128, 64, 32, and 16 pixels.For the ∼70 m wide river, this was sufficient to allow the detection of displaced surface features.Whilst smaller IAs would allow for higher-resolution vector maps, this would also significantly increase noise and thus the number of erroneous correlations.
We processed two configurations based on FFT window deformation and ensemble correlation algorithms at three sampling rates (1, 0.5, and 0.25 Hz), resulting in 6 different LSPIV runs for each scenario.These sampling frequencies resulted in image sequences consisting of 28, 14, and 7 frames which enabled us to experiment with varied frame extraction rates for image-based velocity analysis.Subsampling our original 5 Hz video to lower frame rates (similar to the approach taken by Legleiter & Kinzel, 2021a) was beneficial for detecting velocities, especially for slower-moving phenomena.At a lower frame rate, features in our video had more time to move between frames, resulting in larger displacements that are easier to detect and measure, particularly when dealing with a nominal ground resolution of 1.2 m.Following LSPIV cross-correlation, we post-processed the resultant velocity fields to filter out spurious velocities.Specifically, we utilized filters that removed velocity vectors that differed by 8 × (PIVLab's default threshold) the standard deviation from the mean velocity, and further applied a local median filter threshold of 3 × 3 pixels to remove outliers.Velocity vectors were georeferenced within PIVlab from an image coordinate system back into a projected coordinate reference system (GDA 1994 MGA Zone 55).We used GCP coordinates to assess the accuracy of our georeferencing against actual locations, using 1 m Maxar satellite imagery.

Validation of PIV Velocity Vectors
Validation of LSPIV velocities with 2D hydrodynamic models offers an alternative where the deployment of velocity sensors (e.g., aDcps) is complex, time-consuming, or hazardous.We used velocity predictions from a calibrated 2D HEC-RAS (Brunner et al., 2020) hydraulic model (see Supporting Information S1), which solved the full-momentum Saint-Venant equations.A discharge hydrograph from the Tilpa gauging station (Figure 1) was used as the upstream boundary condition.

Discharge Estimation Using LSPIV Velocities
The standard velocity-area method was used to calculate discharge (Q) (Turnipseed & Sauer, 2010) (see Text S2 in Supporting Information S1, for detail).Water depths were estimated by intersecting the flood extent limits in the satellite imagery with a DEM.Depths at each vertical are computed by subtracting the local bed elevation from the maximum water elevation along a cross-section from a 1 m resolution LiDAR DEM with a vertical and horizontal accuracy of 0.3 and 0.8 m respectively (Geoscience Australia, 2022).The DEM was acquired when the river channel was dry, thus incorporating bathymetry.Discrete PIVLab velocity measurements were interpolated by inverse distance weighting to obtain continuous velocity maps.LSPIV-derived surface velocities were converted to depth-averaged using a specified coefficient α.Hauet et al. (2018) and Le Coz et al. (2010) constrain α between 0.8-1 for deep natural channels experiencing flood discharges.

LSPIV Velocity Accuracy
Stabilization and georectification of frames used in PIV are subject to errors that propagated uncertainty to computed velocity estimates.Maximum, minimum, and mean displacement errors associated with stabilization of extracted frame sequences were 0.420, 0.055, 0.237 and 0.442, 0.15, 0.261 pixels for reach A and B respectively, all less than a single pixel width.Total georectification root mean square error (RMSE) was 0.50 and 0.77 m at reach A and B, respectively.Since our smallest search area was 8 pixels, equivalent to a distance of 9.8 m, our residual georeferencing errors were 5.1% and 7.9% of the spacing between our PIV velocity vectors.
Results of the frame-by-frame analysis of residual motion showed that our GCP locations had a high R 2 of the d(RMSD) metric at both reach A and B (Figures 3a and 3b).The displacement of our GCPs in the stabilized frame sequence further confirmed that all our residual motion at both reaches was within the subpixel range, with average displacements of 0.475 and 0.462 pixels at reach A and B respectively (Figures 3c and 3d).
Figure 4 summarizes the results of the quantitative velocity accuracy assessment of LSPIV (processed using both FFT and Ensemble correlation PIV algorithms at frame rates of 0.25, 0.5, and 1 Hz) against calibrated HEC-RAS 2D model predictions.Regression analysis results in R 2 values of 0.32-0.51(p < 0.001) between LSPIV estimates and HEC-RAS 2D model velocities.To contextualize these results, Legleiter and Kinzel (2021a) attained R 2 values of between 0.34-0.39when comparing aDcp versus satellite video-based PIV velocities across their study area.Root Mean Square Error (RMSE) values at the 0.25 Hz frame rate were 0.18 m s −1 for FFT and 0.22 m s −1 for EC in Reach A, and 0.16 m s −1 for FFT and 0.20 m s −1 for EC in Reach B. Our results indicated a tendency of the FFT algorithm to underestimate flow velocities in both study reaches.Specifically, for Reach A, the Mean Error observed was −0.071 m s −1 suggesting that the FFT algorithm, on average, estimated velocities lower than those predicted by the HEC-RAS model.Similarly, in Reach B, this underestimation persisted, albeit to a slightly lesser degree, with a Mean Error of −0.041 m s −1 .

Discharge Accuracy
The measured discharge was 582.01 m 3 s −1 at the Tilpa gauge at 23:12 UTC on 5/2/2022.LSPIV-based discharge estimates were computed at three cross-sections located in each reach (Figure 1) and ranged from 429.7 m 3 s −1 to 710.1 m 3 s −1 , with median discharges of 536.2 m 3 s −1 (reach A) and 483.4 m 3 s −1 (reach B).Mean absolute percentage error for LSPIV-based discharge estimate was 10% (reach A) and 19.7% (reach B), with significant sensitivity to α.At both reaches, we experimented with α values between 0.7-1.0;previous studies have found that α values of between 0.8-1 are appropriate for computing depth-averaged velocities in natural rivers with a depth of greater than 2 m (Hauet et al., 2018;Vigoureux et al., 2022).At reach A, α values in the range 0.8-0.9minimize the difference between PIV-derived discharge and gauged discharge to within 15%.At reach B a narrow band of α values in the range 0.94-0.97minimize the error, and values in the range 0.9-1.0result in MAE < 10%.

LSPIV Velocity Estimation
We quantified video stabilization uncertainties using the d(RMSD) metric.Our mean values of displacement following stabilization were both within the subpixel range, with mean d(RMSD) being slightly higher than the mean displacement values we obtained following initial stabilization using SIFT (Section 4.1).The presence of a lower mean displacement alongside higher d(RMSD) values highlighted the complex nature of surface flow dynamics and the challenges in capturing these using satellite based LSPIV.It underscored the importance of considering not just the average movement but also the distribution and variability of movement across the video frames.Given our discharge analysis was conducted using velocities derived from a 0.25 Hz sampling rate (i.e., displacements >1 pixel), stabilization errors did not significantly impact the accuracy of our computed velocities (at the 0.25 Hz sampling rate, which provided best correspondence to modeled velocities).
Our sensitivity analysis (Section 3.2.1)highlighted the significance of frame sampling frequency when computing LSPIV velocities, similar to other investigations (e.g., Legleiter & Kinzel, 2021a;Muste et al., 2008; Pearce  et al., 2020).In lieu of reference field velocity measurements, we conducted a direct comparison to 2D model velocity predictions.Statistical analysis of LSPIV velocity deviations, using our best-case scenario of 0.25 Hz processed using the FFT algorithm, showed that LSPIV tended to underestimate velocities compared to 2D model predictions.Our analysis showed variability in the performance metrics (R 2 , RMSE, ME) across different algorithms and settings, indicating the accuracy and reliability of satellite PIV can be context-dependent: R 2 values for both algorithms ranged from 0.29 to 0.51.Potential reasons for this variability include limitations due to satellite image resolution, atmospheric interference, and the inherent limitations of PIV in capturing the complex flow dynamics.Both algorithms generally performed better at Reach B than at Reach A, indicating that channel geometry and flow conditions could impact PIV accuracy.Additionally, some variability may be associated with 2D hydrodynamic model uncertainty (Bates, 2022;Dewals et al., 2023;Pasternack, 2011).Whilst calibrated 2D models are a viable means to assess PIV velocities for cases where flows exceed the safe operating ranges of conventional sensors, we recommend PIV velocity assessment with aDcp measurements.

Discharge Accuracy Assessment
LSPIV-based surface velocities, combined with preexisting, independent information on channel bathymetry, have been successfully used to obtain river discharge estimates in previous studies (e.g., Le Coz et al., 2010;Lewis et al., 2018).Using the velocity-area technique, we estimate discharge with a maximum mean absolute error of 35% which could be reduced to 0.3% and 3.78% at reaches A and B, respectively, by tuning α.The accuracy and precision of our reported discharge estimates compare favorably with Sun et al. (2010) and Lewis et al. (2018) who computed river discharges using LSPIV-based measurements to within −5 to 7% and <20% respectively.The ephemeral nature of the River Darling at Tilpa is advantageous for acquiring high accuracy bathymetric topography, here using airborne LiDAR.In other ephemeral locations, lower resolution data sets with near-global coverage could be used for large rivers.In temperate and tropical locations, direct bathymetric surveys or bathymetry derived from multispectral satellite imagery, altimetry (e.g., Liu et al., 2020;Moramarco et al., 2019) or the inference of depths from PIV-derived velocities using a flow resistance equation-based framework (Legleiter & Kinzel, 2021b) would be required.Despite these additional data demands, our results demonstrate that satellite-based optical video sensors could be deployed for near real-time estimation of riverine velocity and discharge, within tolerable uncertainties common to traditional discharge estimation techniques.

Variability of Surface Coefficient Values, α
Our satellite-video based LSPIV discharge estimation procedure yielded promising results, in terms of absolute flow magnitude, but the selection of the coefficient (α), used to convert surface to depth-averaged velocities, remains a source of uncertainty in discharge estimation.Fulton et al. (2020), Moramarco et al. (2017), andWelber et al. (2016) all observed local variability of α (0.52-0.78; 0.85-1.05,and 0.71-0.92respectively) when estimating discharge using non-contact techniques, attributable to variations in stage (especially during higher flows due to changes in wetted channel perimeter), channel geometry, slope, and channel alignment.Significant shifts in the error of LSPIV-based discharges due to variations in α indicated that sufficient cross-section specificity in defining α is critical to our technique.
Higher α values generally led to less error in our discharge estimates, particularly in areas where LSPIV velocities differed substantially from velocities in the hydraulic model sued for benchmarking.Hauet et al., 's. (2018) recommendation of α values based on a river's hydraulic radius and depth, with a noted uncertainty of ±15% at a 90% confidence level, are consistent with the optimal values reported here.However, without an empirically formulated, river-specific α based on in situ measurements, the appropriate values of α remains largely unclear (Legleiter et al., 2023).
When estimating flood flows in remote locations where remote sensing instruments are the sole source of depths (i.e., derived from a DEM), experimenting with values provided by Rantz (1982) (α = 0.85 or 0.86), Turnipseed and Sauer (2010) (α = 0.84-0.90),and, in extreme cases, α > 1 due to non-standard velocity distributions (see, e.g., Moramarco et al., 2017) is a sensible approach to improve the precision of flow measurements from surface velocimetry techniques.On average, in our study the α values that led to the closest approximations of observed discharge were all less than unity, indicating our velocity-depth distributions could be well approximated using logarithmic or power laws.The variability of our best fitting, cross-section averaged α at our reaches implies that the commonly used default value of 0.85 is not always appropriate in field conditions where spatial heterogeneities in channel beds have a significant impact on velocity profiles.Although we provide a method for assessing the variability of α, calibration of site-specific α values based on traditional contact measurements remains the preferred solution for accurate discharge estimation.

Conclusion
Satellite-based PIV presents a promising tool for estimating river discharges during cloud-free conditions.Key to constraining uncertainty and enhancing the accuracy and reliability of LSPIV-derived velocity estimates is the stabilization of satellite video frames and the independent assessment of residual error, particularly for sub-pixel displacements.Performance metrics from the comparison of PIV velocity magnitude vectors against 2D model predictions of surface velocity exhibited reasonable correspondence.The FFT algorithm at a frame rate of 0.25 Hz, revealed best correspondence, but differences between study reaches highlight how site-specific characteristics can influence LSPIV performance.The observed R 2 values (0.3-0.5) highlight the need for careful consideration in the application of PIV techniques, particularly for low-frame-rate satellite videos.LSPIV accuracy also depends on α.Using realistic α values (0.7-1.0) from literature, our resulting errors were −6.9 m 3 s −1 and −85.6 m 3 s −1 , and biases were −0.01 and −0.15, at our two study reaches, respectively.Despite these uncertainties, when combined with high-resolution topographic data, the ability of satellite-based LSPIV to provide large-scale, non-intrusive river surface discharge measurements in inaccessible or dangerous areas remains a compelling advantage.The level of accuracy offers a promising foundation for enhancing LSPIV methodologies; uncertainties are comparable to traditional methods and avoid the need for extrapolation of rating curves during high flow conditions.While acknowledging the necessity for further ground truthing to assess uncertainty, there is considerable potential for satellite video to be used to estimate discharge.We extend our sincere gratitude to the editor, Valeriy Ivanov, and the insightful reviewers, Carl Legleiter and Matthew Perks for their constructive feedback and valuable suggestions that significantly enhanced the quality and clarity of this manuscript.We would also like to thank Robert Ljubičić for providing the script used in our stabilization analysis.Christopher Masafu was funded by the Lord Kelvin Adam Smith (LKAS) PhD scholarship at the University of Glasgow.

Figure 1 .
Figure 1.Study area indicating investigated reaches A and B.

Figure 2 .
Figure 2. Discharge estimation and validation workflow.Green shaded boxes show the required input data.

Figure 3 .
Figure 3. Displacement versus RMSD for Ground Control Points in reach A (a) and B (b).Each color-coded scatter point corresponds to a different GCP, showing how the displacement within a predefined area around each point affects the RMSD.Stabilization effectiveness across Ground Control Points (GCPs) in reach A (c) and B (d).These plots present the calculated RMSD values for each GCP across different frames, showing the stabilization performance.Each point represents the distance (distortion measure) at a GCP for a specific frame, offering insights into the temporal consistency of stabilization accuracy.

Figure 4 .
Figure 4. Comparative visualization of river flow velocities.(a, d): HEC-RAS 2D model-derived velocity fields at reaches A and B. (b, c, e, and f): Surface velocity vectors derived from satellite-based LSPIV, computed using different algorithms.The LSPIV velocity vectors are positioned to correspond precisely with the locations used in the HEC-RAS 2D model.(g-h): scatter plots showing correlations between LSPIV velocities and HEC-RAS 2D model predictions.