Discharge Estimation Using Video Recordings from Small Unoccupied Aircraft Systems

Measuring discharge during flashy flood events is difficult and dangerous in the Southwest US, owing to high flow velocity and suspended sediment loads. Flowmeters based on acoustic principles often perform poorly in the hyper-concentrated sediment laden flows typical of flash floods in this region. Small unoccupied aircraft systems (sUAS) have been increasingly used in the study of hydraulics and hydrology because they are affordable, stable, reliable, and can be easily and safely operated by a trained crew. With sUAS positioned directly above water surface, onboard cameras can capture high-resolution video, which can be used to calculate the two-dimensional (2D) distribution of surface velocities within the camera’s field of view computed by image velocimetry techniques including large-scale particle image velocimetry (LSPIV) (Fujita et al. 1998). Discharge can then be estimated by converting LSPIV-obtained surface velocities occurring along a cross section into mean cross-sectional velocities and multiplying by measured cross-sectional area (e.g., from a stage-area relation).

Techniques for the acquisition and processing of LSPIV surface velocities are well documented and described in the literature, including (1) how to collect oblique (Fujita et al. 1998; Bradley et al. 2002) and/or nadir (Perks et al. 2016; Legleiter et al. 2017; Lewis et al. 2018) imagery from airborne (Bandini et al. 2021) or ground-mounted (Creutin et al. 2003), handheld (Le Coz et al. 2016), and citizen science collected (Guillén et al. 2017) cameras; (2) how to geo-rectify imagery using a pinhole camera model (Jodeau et al. 2008); and (3) the comparison and evaluation of results among several image velocimetry methods (Pearce et al. 2020; Perks et al. 2020) and/or seeding approaches (Tauro et al. 2016; Pizarro et al. 2020).

To convert the surface velocities to depth-averaged velocity, Rantz (1982) employed a simple correction coefficient between the surface- and depth-averaged velocities occurring in a cross section. Rantz’s velocity index was found to be 0.85–0.86 for natural channels and 0.90 for smooth artificial channels based on the field measurements of discharges in streams and canals through the conterminous US (Hulsing et al. 1966). For flows that occur in smooth channels with vertical walls, the velocity indexes are lowest in the middle of a cross section and may reach over 1.05 toward the sidewalls due to the presence of secondary currents. This finding suggests that the velocity index can vary significantly over a given cross section, depending on the characteristics of secondary currents generated from the sidewalls, vegetation, and other local disturbances, such as bed forms.

Hauet et al. (2018) evaluated how the velocity index correlated with depth, roughness, and other flow characteristics through review of many discharge measurements collected at gauging locations in France. Through sensitivity analysis, they concluded that a velocity index of 0.80 should be used for rivers less than 2.0 m deep, and a velocity index of 0.90 should be used for water depth over 2.0 m (or concrete channels, regardless of flow depth). Similarly, Polatel (2005) also found that the velocity index was closely related to channel boundary material type and flow depth, concluding that the velocity indices ranged from 0.789 to 0.928, with larger indices for smooth-bed and large-flow-depth conditions, and smaller indices for rough-bed, low-flow-depth conditions.

A short video (e.g., 1–3 min) can be decomposed into a series of consecutive images and used to derive a 2D instantaneous surface velocity field in the streamwise and transverse directions. If time between image frame is small enough, turbulence characteristics can be calculated from the time series of instantaneous velocities. Lewis and Rhoads (2015) used the 2D surface velocity distribution obtained from LSPIV to investigate the turbulence characteristics at river confluences. Johnson and Cowen (2017a, b) computed friction velocity from the surface turbulence dissipation rate (TDR) directly, and Johnson and Cowen (2017b) used the ratio of surface velocity and friction velocity to estimate the depth-averaged velocity based on the power-law velocity profile. In a like manner, Jin and Liao (2019) also used TDR to estimate friction velocity and then flow depth, assuming the Manning’s equation is valid on the Milwaukee River.

The “probability concept” approach (Chiu 1988; Moramarco et al. 2004; Fulton and Ostrowski 2008; Moramarco and Singh 2010; Moramarco et al. 2017; Fulton et al. 2020a, b) uses Shannon’s entropy theory to correlate cross-sectional maximum velocity and depth-averaged velocity along each vertical. By maximizing entropy to estimate a probability distribution of velocities in the cross section, a relationship between the maximum and mean velocities can be developed (Chiu 1988; Farina et al. 2014). Chiu and Tung (2002) and Bjerklie et al. (2020) posited that the depth of the maximum velocity relative to the depth in the vertical containing the maximum velocity does not change with flow discharge, water level, or velocity for simple channel geometries where flows are within their banks. However, further work is necessary to verify this hypothesis.

The accuracy computing TDR from LSPIV using sUAS imagery is dependent upon the quality of surface velocity seeding and texture. Natural seeding includes wave crests, vortices, bubbles, or natural debris (Dal Sasso et al. 2020). It can difficult to achieve high levels of seeding in natural environments (Dal Sasso et al. 2020); therefore, it may be necessary to introduce artificial seeding. Dal Sasso et al. (2020), Huang et al. (1997), and Keane and Adrian (1992) all stressed that the quality of the seeding distribution is the most critical factor in using PIV algorithms.

Although approaches leveraging TDR have produced promising results in the laboratory and at least two river systems (Johnson and Cowen 2017b; Jin and Liao 2019), systematic evaluation of the use of TDR derived from surface velocities and applied to computing stream discharge in rivers and channels of various sizes is lacking. Because each river differs in geometry and flow characteristics, the velocity index is likely not constant, and the TDR may not obey the standard linear distribution expected from theory and laboratory experiments.

The purpose of this study is to evaluate the accuracy of discharge estimation methods using TDR derived from imagery taken from sUAS in field applications. Specifically, we evaluate four methods of computing depth-averaged velocities from TDR derived from LSPIV analysis of videos: (1) constant-velocity index, (2) logarithmic velocity profile, (3) power-law velocity profile, and (4) entropy method at each of eight field sites representing engineered canals and natural rivers. We describe how to determine surface TDR from LSPIV analysis and apply each method to compute friction velocity, depth-averaged velocity, and in conjunction with channel area estimate discharge. We present the results of comparisons of each method against truth measurements of discharge captured with the conventional means at each site. We conclude by discussing of the implications of using TDR to estimate discharge as well as considerations and limitations for field application of the methods.

Eight field sites located in Arizona, New Mexico, California, and Maine were studied. A summary of site characteristics is presented in Table 1. Seven of the field sites were colocated with established USGS gauging stations. The Androscoggin River field site was located 2.8 km upstream of the nearest USGS gauging station. All but the Androscoggin River site are located in regions with typically flashy flood responses (Duan et al. 2017). The Wellton-Mohawk Main Outlet Drain (Fig. 1) is representative of what the other canal sites typically look like. The sUAS images in Fig. 2 show flow conditions during this study for the three natural river sites on the Agua Fria and Androscoggin Rivers.

At each site, ground control points (GCPs) (e.g., solid circles in Fig. 1) were placed to calibrate the pixel ground scale resolution, which is required to translate pixel displacements computed from LSPIV into velocity units. The distances between the GCPs on each side of the channel were measured and recorded using either a tape measure or a laser rangefinder accurate to $\pm 0.15\mathrm{m}$.

All sUAS video footage was collected with the camera pointing perpendicular to the water surface (nadir) with at least a 1 m of bank visible on each side. The sUAS pilot also ensured that at least two GCPs were visible and centered in the field of view of the camera. The flight altitude varied depending upon the channel width; for example, wider channels required the pilot to fly at a higher elevation so that the entire channel width including 1 m of bank on each side was visible in the field of view.

Seeding was necessary for applying the LSPIV algorithm in all but the Agua Fria and Androscoggin River sites. Both rice cereal and dry wheat straw were evaluated as seeding material. To uniformly distribute the rice cereal across the channel, a leaf blower was modified by cutting the blower hose and adding a valve to regulate the rate at which seeding material was dispersed (Fig. 3). Above the valve, a container to hold seeding material was installed. The blower functioned well, distributing seed material up to 15 m. Because the straw could not fit in the blower, it had to be spread by hand, and often formed clumps and dead areas where no seeding was visible. Seeding was best distributed across the channel when there was a bridge directly upstream of the measurement location, which allowed the leaf blower operator to stand directly over the middle of the channel and spray seeding back and forth across the entire channel width. Seeding from the channel banks led to inhomogeneous distribution of seeding across the channel.

Video footage from seven of the eight study sites was captured using a Parrot Anafi sUAS (Paris, France) shooting at $\mathrm{3,840}\times \mathrm{2,160}$ at 30 frames per second (Anafi Parrot Drone 2021). The Androscoggin River field site used a Da Jiang innovations (DJI) Mavic Pro Government Edition drone shooting at $\mathrm{4,096}\times \mathrm{2,160}$ at 24 frames per second to capture videos (DJI 2023). Video was collected for as long as wind and seeding conditions permitted (generally 20 to 30 s at a minimum).

To evaluate the sUAS derived methods in this study, flow velocities and discharge were measured at each field site at the time of sUAS data collection using either a Sontek M9 (San Diego) acoustic Doppler current profiler (ADCP) or Sontek Flowtracker operated by USGS hydrologists. Subsurface velocity data were collected either directly before or after video footage was captured with the sUAS. All ADCP data were postprocessed using Velocity Mapping Toolbox (VMT) version 4.09 (Parsons et al. 2013) to obtain the bathymetry and depth-averaged velocities. All discharge values measured by ADCP or acoustic Doppler velocimeter (ADV) were reviewed and processed using the QRev version 4.32 software application (Mueller 2020) according to current USGS policy (Turnipseed and Sauer 2010).

We treated the USGS reported discharge values as the truth reference values for evaluating the accuracies of sUAS derived discharges. The uncertainties of ADCP measured discharge at all sites were described by Engel et al. (2022) and are summarized in Table 3. In all but two cases (Table 3) the USGS reported quality was better than 9%. All input data for this study, including the sUAS videos, annotated images of GCP locations and distances, and subsurface velocity measurement files, are available (Engel et al. 2022).

This study used the Rectification of Image Velocity Results (RIVeR) toolbox version 2.5 (Patalano et al. 2017) to extract the frames from the video clips recorded from the sUAS and process LSPIV. To evaluate turbulent characteristics, all available frames were extracted because TDR relies on instantaneous velocity fluctuations. RIVeR uses an embedded version of PIVLab (Thielicke and Stamhuis 2014) to compute the time-resolved image velocimetry results. Image frame stabilization was performed using the approach of Farid and Woodward (2007), adapted to mask out the water (i.e., noncamera platform) motion from analysis in order to remove excess camera platform movement (e.g., the sUAS was not perfectly stable).

Each data set was preprocessed to enhance visualization of water surface seeding texture of the images. After conducting several trials and comparing the estimated surface velocities with the measured velocities from the ADCP, it was concluded that a contrast limited histogram equalization (CLAHE) with background subtraction produced the best particle tracking results. For all sites, a CLAHE window size of 20 pixels, corresponding to ground distances ranging from 0.04 to 0.73 m was used because the application of CLAHE was relatively insensitive to changes in window size.

A fast Fourier transform (FFT) multiwindow (three-pass) deformation PIV algorithm was utilized. Passes one, two, and three window sizes were set to ground distances ranging from 0.27 to 4.7 m (128 pixel), 0.13 to 2.35 m (64 pixel), and 0.07 to 1.17 m (32 pixel), respectively, depending on the site and pixel ground scale distances. These settings produced a displacement velocity vector at the center of each 32-pixel window in the region of interest, resulting in maximum vector spacings ranging from 0.067 m at the Sile Main Canal to 1.17 m at the Androscoggin River. Images were calibrated into real-world units using the GCPs. The output LSPIV analyses for each site comprised a time series of 2D surface velocities magnitudes ($\mathrm{m}/\mathrm{s}$) in the streamwise ($u$) and cross-stream ($v$) directions for each frame pair of the processed video. The time-series frequency was dictated by video framerate and ranged from 23.98 to 30 Hz (Engel et al. 2022).

We compared seeding material influence on LSPIV velocity results to ADCP near-surface velocities at a cross section for the Wellton-Mohawk Main Outlet Drain site. The Wellton-Mohawk results shown in Fig. 4 are representatives of our findings, in which the root-mean square error (RMSE) using rice cereal, straw, and no seedings were 0.090, 0.094, and $0.39\mathrm{m}/\mathrm{s}$, and the corresponding normalized RMSE (NRMSE) were 9.26%, 9.67%, and 40.62%, respectively.

Nadir sUAS images at Reservation Main Canal illustrate the difference in surface texture between seeding material type (Fig. 5). Comparison of velocities by seeding material indicated that both straw and rice cereal are preferable with similar accuracy to no seeding material, but seeding is required in low-velocity, low-turbidity channels. For this reason, rice cereal was used because it was more economical and easier to distribute using the seeding distribution system (Fig. 3), except for the Gila and Androscoggin Rivers which were (1) too large to distribute seeding material across the measurement area, and (2) had better surface texture because of increased turbidity and velocity.

For each instantaneous velocity contained in the LSPIV-derived surface velocity time series ($u_{i,t}$, $v_{i,t}$) along a given channel cross section, the streamwise ($u_{i,t}^{\prime }$) and transverse ($v_{i,t}^{\prime }$) fluctuations can be found by using Eqs. (1) and (2)where primes = instantaneous velocity fluctuations; overbars = time-averaging; and subscripts $i$ and $t=i$th position along the cross section at the $t$th time step. Assuming the turbulence flow is homogenous and isotropic, the TDR can be computed from the one-dimensional (1D) longitudinal ($E_{11}$) or transverse ($E_{22}$) velocity spectra from the velocities fluctuation expressed in 1D wave-number space, ${\kappa }_1$, whereand $\mathrm{\Lambda }$ = wavelength. Based on Kolomogorov’s second similarity hypothesis, within the inertia subrange the statistics of motion are uniquely determined by the TDR, and the temporal energy spectra can be found as follows:where $C_1=9/55\times 1.5C$; and $C_1^{\prime }=4/3C_1$, where $C=1.5$ and $\epsilon$ is TDR (Pope 2014). TDR was estimated by plotting the velocity spectra ($E_{11}$, $E_{22}$) in the log space versus the wave number (Fig. 6), and taking the average of all spectra conforming to the best fit of a straight line with a slope of $-5/3$ in the inertial subrange (i.e., the portion of the spectra in which $E_{11}$ or $E_{22}$ conforms to a $-5/3$ slope line). Each line in Fig. 6 represents the results calculated from a time series of measured surface velocity at a transverse location in a section. The dissipation rate was averaged over several cross sections at a fixed transversal location. We found that in all sites, the streamwise velocity spectra were much wider in the inertia subrange than the transverse velocity spectra. The transverse spectra ($E_{22}$) tended to better follow the $-5/3$ slope than the streamwise spectra ($E_{11}$), the same as observed by Johnson and Cowen (2017a) in laboratory shallow open-channel flow experiments. Therefore, the transverse velocity spectra were selected for the estimation of turbulence dissipation rate for all further computations in this paper, consistent with the method used by Johnson and Cowen (2017a).

Nezu (1977) proposed a universal relation for the vertical distribution of turbulence dissipation in open-channel flows:where ${\epsilon }_{zi}$ = turbulence energy dissipation rate at a location ($z$) along the $i$th vertical; angled brackets = ensemble averaging; $u_{*,i}$ = friction velocity along the $i$th vertical; and $h_i$ = flow depth at the $i$th vertical in a cross section. Values of $E$ vary and typically range from 8.43 to as high as 13.4 (Pope 2014; Sukhodolov et al. 2006; Johnson and Cowen 2017a). We adopted the method of Johnson and Cowen (2017a) with $E=9.8$ by solving Eq. (5) for friction velocity using surface TDR (i.e., $z/h_i=1.0$)

Furthermore, Johnson and Cowen (2017a) found that a constant ratio 1.24 between the near-surface and surface dissipation rate fit experimental data and provided a means to convert surface TDR into near-surface TDR. Although they rightly suggested further investigation is needed to confirm their laboratory results, we found good agreement in our observations (“Results” section). For our study, measured dissipation rates at the surface were converted to subsurface dissipation rates for each vertical $i$ in a cross section by multiplying by 1.24, and then Eq. (7) was applied to calculate the friction velocity, in which $z/h$ is set to 0.9.

We employed four methods to compute depth-averaged velocities: (1) constant-velocity index, (2) logarithmic velocity profile, (3) power-law profile, and (4) entropy method. Except for the entropy method, all the methods assumed a 1D velocity profile is valid for each vertical $i$ in each cross section. The following are detailed descriptions of each method.

When applying the constant-velocity index, the logarithmic law and the power-law methods, the depth-averaged velocity was calculated for each vertical and the total discharge is the summation of all the subdischarges according to the mean-section methodwhere $n$ = number of subsections in a cross section; and $b_i$ = distance between the verticals $i$ and $i+1$. Discharge was computed from the entropy method as the product of cross-sectional mean velocity determined by Eq. (20) and cross-sectional area

Friction velocity was calculated for each vertical using Eq. (7), and then the depth-averaged velocities were estimated according to the four methods described previously. Flow discharges were estimated from Eqs. (21) and (22).

A flowchart summarizing the process of computing discharge for each method, labeled with the key equations used in each step, is shown in Fig. 7. Table 3 summarizes the estimated discharges from the four methods and the corresponding conventionally measured discharges, which were considered truth for evaluation purposes. The normalized errors of estimated discharges were calculated from the estimated ($Q_e$) and the observed ($Q_o$) discharges as $(Q_e-Q_o)/Q_o$, and are summarized in Table 4 together with the channel types and $B/H$ ratios. For each field site, the estimated discharges are reported by method, and the methods with the least normalized error are in bold.

When applying the logarithmic law profile, a roughness height $z_0$ needs to be selected in Eqs. (9) and (10). Because bed sediment size was not available, the $z_0$ was assumed to be proportional to local flow depth as $z_0=ah$, where $a$ was optimized using field measurements and varied from 0.04 to 0.1, which is consistent with other research (Bagherimiyab and Lemmin 2013). Because Eq. (9) is theoretically only valid to $0.5h$ due to free surface effects, the log-wake law [Eq. (10)] was also evaluated. However, the log-law without wake correction performed better than the log-wake law at all sites in this study (Table 4).

Moreover, the log-law estimated discharges at the four sites with $B/H>10$—including three natural rivers and one earthen canal—performed best relative to the other sites with $B/H\le 10$. At the four sites with $B/H>10$ the values of $a$ were relatively uniform $0.063<a\le 0.10$ with an average of 0.076 ($\mathrm{SE}=0.0083$, $n=4$, and $95\%\mathrm{CI}=0.026$). However, at the sites with $B/H\le 10$, the values of $a$ were less uniform $0.007<a\le 0.04$ with an average of 0.014 ($\mathrm{SE}=0.0988$, $n=4$, and $95\%\mathrm{CI}=0.315$). This suggests the log-law method may be more accurate for estimating the depth-averaged velocity for $B/H>10$ presuming $a$ can be optimized with known discharge reference data.

The results from the power-law method were not as accurate as the either log-law method (Table 4), and values of $H$ at multiple verticals at the Androscoggin, Reservation Main and Cochiti East Side Main Canal sites were less than 1.0, leading to $m<0$, which is nonsensical. Thus, the depth-averaged velocity was assumed equal to the observed surface velocity for these cases. The power-law profile performed best at Reservation Main Canal Site (3.41% error), but all other sites had larger reported percent errors (Table 4). We suspect that the computed skin friction factor ($C_f$) based on surface velocities performed poorly owing to the assumption that there was no form drag present (Johnson and Cowen 2017b). Calibration of $G$ and further evaluation of the empirical relation between $G$ and $H$ could be re-examined as more laboratory and field data become available.

For the entropy method, we assumed the maximum velocity occurred at the free surface, and the observed maximum velocity at surface was the maximum velocity in the entire cross section. The entropy method performed best for the concrete canal sites with $4<B/H<5$. However, at the Sile Main Canal ($B/H=1.64$), flow was noted in the field to be 2D, and the maximum velocity likely occurred below the water surface, violating our assumptions and increasing the percent error of discharges at this site to a higher value (9.33%) relative to the other sites with $B/H<5$. Estimated discharges using the logarithmic law method [Eq. (9)] were much more accurate than the discharges computed from the entropy method for the sites where $B/H>10$.

The constant-velocity index method performed best at the Androscoggin River site (0.14% error) (Tables 4 and 5). However, the absolute errors at the four sites with $B/H>10$ ranged from 0.14% to 15.14%. At the other four sites, $B/H<5$, and the errors ranged from 12.00% to 23.25%. When comparing with the results from the entropy method, the constant-velocity index method was more accurate than the entropy method for $B/H>10$, but less than accurate for $B/H<5$.

Based on results from the eight field sites studied here, the log-law, constant-velocity index, and the entropy method performed best at estimating discharge from measured surface velocity and TDR. The most accurate results from this study were achieved using the log-law or log-wake-law approach with the $a$ value calibrated to roughness values at each site. By summarizing the results of optimization of the $a$ parameter, we find it may be possible to select $a$ values without any precalibration. The average of precalibrated values for the sites in this study having $B/H>10$ (Table 4) was 0.076, corresponding to a roughness height $z_0=0.076\mathrm{h}$. Similarly, for $B/H<5$ (Table 4) at the four concrete canal sites, the constant $a$ values were 0.003, 0.007, 0.04, and 0.07, with average and median values of 0.014 and 0.007, respectively.

Discharge was recomputed using the summarized values of $a$ for each site based on $B/H$ (Table 5). When choosing $a=0.076$ for $B/H>10$, the largest error was 21.11% at the Reservation Main Canal, and the errors at other three sites ranged from 0.9% to 11.1%. These errors were smaller than the results using the entropy method but larger than the errors using the constant-velocity index method at three of the four sites, indicating that for natural rivers with $B/H>10$, the constant-velocity index method is more accurate than the log-law or log-wake law methods regardless when selecting a constant $a=0.076$.

For $B/H<5$, errors using the entropy method were smaller than either the average value of $a$ or the constant-velocity index methods. Without knowing $z_0$ in log-law or log-wake law, the results suggest the constant-velocity index method performed the best for natural rivers of $B/H>10$, whereas the entropy method was the best for concrete canals with $B/H<5$. Based on our results, $a=0.076$ should be used for the log-law method for natural rivers of $B/H>10$. The log-law method may not yield accurate results when using $a=0.007$ or $a=0.014$ for concrete canals.

This study evaluated the methods of river discharge estimation using videos recorded by sUAS at eight unique field sites. Of the eight study sites, four were engineered concrete channels, three were natural channels, and one was an engineered earthen channel. Each of these sites differed in channel geometry, flow, and seeding characteristics. With the instantaneous and average surface velocity distributions computed from LSPIV at each site, the TDR was estimated from 1D transverse velocity spectra for calculating the friction velocity. Depth-averaged velocity was calculated using four methods: the constant-velocity index, the logarithmic law (with or without wake correction), the power-law, and the entropy method. From our findings, these general recommendations emerged:

One should select the best approach for velocity calculations using $B/H$ ratios and channel types. The recommended $a$ value for $B/H>10$ is based on the data from these eight sites, and the maximum error can be as high as 21.11%. Additional field data are needed to verify if $a=0.076$ is the most suitable value for using the log-law.

This study did not evaluate sites with $5<B/H<=10$, and thus does not draw conclusions about this range. More research is needed to search for and evaluate approaches to precalibrate the roughness parameter used in the logarithmic law method and to include sites of medium $B/H$ ratios. Although limited in scope, this research shows that under ideal seeding and flow conditions, discharge can be accurately calculated from TDR and surface flow velocities obtained from the images collected remotely.

The original video recordings at eight sites are published in accordance with the USGS data retention policy online (Engel et al. 2022). The Python program for calculating discharge using methods described in the paper is also published by the University of Arizona data repository (Duan et al. 2022). Other data (e.g., RIVeR software) that support the findings of this study are available from the corresponding author upon reasonable request.

This research is funded by USGS Water Resource Research Center 104B Grant to the University of Arizona. USGS personnel in the Yuma, Arizona, and Albuquerque, New Mexico, field offices have provided all the ADCP and Flowtracker velocity data for matching the data collected by the sUAS. We are grateful to the funding support and technical assistance from USGS. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government.