*3.1. Data*

#### 3.1.1. UAV Data Collection

Aerial images were acquired via a DJI Phantom 3 Professional UAV (DJI, Shenzhen, China), manually flown by sight and equipped with an independent camera set to take photographs at a fixed interval (Table 1). The DJI Phantom 3 Professional UAV has a vertical takeoff weight of 1280 g and can reach a maximum speed of 16 m/s. Fully charged batteries provide a maximum flight time of 23 min. The camera was a standard lens mounted on a Cardan suspension and based on complementary metal-oxide

semiconductor (CMOS) image sensor technology. The FC300X sensor provides images with 12.4 million effective pixels and has separate channels for red, green, and blue light. The spatial and elevation resolution of the acquired UAV data was confirmed to be centimeter-level in our previous study [33].


**Table 1.** Key parameters of the unmanned aerial vehicle (UAV).

UAV images were acquired twice, in August of 2017 and 2018, the dry season of the study area with no flow through the channel. In each survey, the same study area was surveyed in two flight missions pre-programmed with Pix4D Capture software for several waypoints to achieve a 90% image overlap along the track. The flight altitude was on average 70 m and was a trade-o ff between the size of the studied channel and the desired spatial detail of approximately 3 cm per pixel. The two in situ surveys produced 380 and 384 UAV images separately. A total of 8 Ground Control Points (GCPs), 4 for each river channel, were measured by Real-Time Kinematic in situ surveys. The root-mean-square error (RMSE) between the ground measurements and UAV data was ±2.79 cm in the horizontal direction, which is feasible to identify the stone movement in the exposed channel.

#### 3.1.2. UAV Data Processing

The stereoscopic images were processed using the rapid and automatic professional processing software Pix4Dmapper. The software supports not only UAV data but also aerial photography, oblique photography, and close-range photogrammetry. Through data importing, initial processing, and point cloud encryption, Pix4D mapper provides a high-resolution Digital Ortho Map (DOM) and a Digital Surface Model (DSM) [34,35] (Figure 3). The processing steps include (1) determining the internal orientation element and calibrating the camera, (2) searching and fast matching the same name point employing the scale-invariant feature transform (SIFT) algorithm [36], (3) using the position information of the same name point and the image to adjust the regional network and to restore the position and posture of the image, (4) performing aerial triangulation encryption and generating an image point cloud using the control point or the matching point, and (5) generating the DSM from an image point cloud and then generating the DOM.

**Figure 3.** (**a**) The derived Digital Surface Model (DSM), (**b**) the three-dimensional triangulation process by UAV and the flight path, and (**c**) the derived Digital Orthophoto Map (DOM).

#### 3.1.3. In Situ Survey Data

In situ surveys were conducted in August of 2017 and 2018, the same time as UAV data collection. During the surveyed dry season, river channel H and S were completely exposed. Evident watermarks were found both on the bank of some reaches and the downstream culvert of the studied river channels. There is also a clear vegetation growth boundary to help determine the water level. The watermark conditions of the two years on the riverbank were compared and measured by tape to determine the maximum water level of the peak flood throughout the year. Cross sections with the determined maximum water levels were selected from the measured river reaches for subsequent calculation. The watermarks on the culvert were also measured to determine the maximum water level when the peak flood passed through the downstream culvert in each study site for validation.

#### *3.2. Critical Velocity of Particle Transport*

In arid ungauged regions, upstream ephemeral rivers rarely flow, and river channels are scarcely affected by human activities. Limited stones scattered in the river channel, especially identifiable moving stones, are regarded as only moved by water flow. When the flood flow promotes a stone to move, the movement of the stone is less likely a ffected by the surrounding stones and bed materials. Thus, the description of a moving stone targeted at the incipience by a hydraulically rough wall-shear flow can be illustrated in Figure 4. The catalyst forces for the moving stone are the hydrodynamic drag force F D and the lift force FL. The submerged weight W brings the stabilizing condition to the stone. The three forces reach equilibrium and have the relationship [37]:

$$F\_D = \tan \mathcal{Q} (\mathcal{W} - F\_L),\tag{1}$$

where ∅ is the angle of repose.

The submerged weight W, drag F D, and lift FL can be expressed as

$$F\_D = C\_D \frac{\pi d^2}{4} \frac{\rho Ul\_0^2}{2} ,\tag{2}$$

$$F\_L = C\_L \frac{\pi d^2}{4} \frac{\rho Ul\_0^2}{2} \tag{3}$$

$$\mathcal{W} = \frac{1}{6}(\rho\_{\text{s}} - \rho)\pi d^3 \,, \tag{4}$$

where U0 is the fluid velocity at the bottom of the channel, C D and CL are respectively the drag and lift coe fficients, d is particle nominal diameter, ρ*s* is stone density, and ρ is liquid density.

In the studied river channels, size compositions of the surface sediment mainly include sand fraction (<2 mm), and gravel fraction (2–20 mm). Cobbles (20–200 mm) and boulders (>200 mm) are scattered on the riverbed, and only some could be identified from UAV images due to the resolution. Most of the stones in the channel are granite with an average empirical density ρ*s* = 2.80 g/cm<sup>3</sup> [38,39].

In the studied river channels H and S, four cross sections (sections A, B, C, and D) were analyzed and measured in situ individually (Figure 5). The research reach refers to the range of 5 meters upstream and downstream of each section. In each research reach, moved and unmoved stones were identified by comparing the high-resolution orthoimages taken over two years. For each moved stone, the average flow rate when it starts to move can express the water flow intensity. The power law and Prandt–von Karman logarithmic velocity were used to simulate the velocity profile and to calculate the critical velocity of the incipient motion of each moving stone separately.

**Figure 5.** (**b**) and (**c**) present the location of cross sections and the downstream culvert in river channels H and S. (**a**) and (**d**) are fieldwork photos of measuring flood max water level.

#### 3.2.1. Logarithmic Distribution of the Flow Velocity

Assuming that the flow velocity profile in the channel follows the power law representing a logarithmic velocity distribution, the average velocity through a cross section of the open channel with turbulence can be expressed using the Einstein unified formula [40]:

$$\frac{\dot{M}}{\dot{M}\_\*} = 5.75 \lg \left( \frac{12.27 R' \chi}{k\_s} \right) \tag{5}$$

where ks is the height of the protrusion of the boundary roughness, also known as the side-wall roughness or the riverbed roughness. *U*∗ is the shear velocity corresponding to the sand resistance and *U*∗ = *gRJ*. χ is a correction factor related to ks/δ. δ is the calculated thickness of the viscous bottom layer and δ = 11.6ν/*U*<sup>∗</sup>. R' is the hydraulic radius corresponding to the sand resistance.

The critical initial shear stress can be expressed as:

$$
\Theta = \frac{\tau\_c}{(\rho\_s - \rho)D} = f\left(\frac{\mathcal{U}\_\* D}{\nu}\right). \tag{6}
$$

Sand waves have not ye<sup>t</sup> formed, so let R' = R and Uc represent the average vertical velocity of the critical condition of incipient motion. The critical initial condition expressed by the average critical velocity through the cross section Uc can be derived as:

$$\frac{U\_{\rm c}}{\sqrt{\frac{\rho\_{\rm s} - \rho}{\rho} gD}} = 5.75 \sqrt{\Theta} \lg \left(\frac{12.75 R \chi}{k\_{\rm s}}\right) \tag{7}$$

where Θ = *f*(*<sup>R</sup>*∗) ≈ 0.06 under a turbulent condition (*R*∗ 500). The density of water ρ = 1 g/cm3. g is the acceleration due to gravity and has a value of 9.8 m/s2. D is the nominal diameter of stones. In the case of a wide and shallow channel, the hydraulic radius R equals the water depth, the Einstein correction factor χ = 1.0, and the roughness height ks = D.

#### 3.2.2. Exponential Distribution of the Flow Velocity

Assuming the average velocity on the cross section is exponentially distributed on the vertical line, the critical initial velocity can be written as

$$
\hbar L\_{\varepsilon} = K \sqrt{\frac{\rho\_s - \rho}{\rho} gD(\frac{h}{D})} \,\,^{1/6}\prime\,\,\,\,\tag{8}
$$

where K is a comprehensive coefficient that considers the submerged weight of the stones, lifting force, drag force, and effect of the irregular shape of stones on the three force coefficients. Shamov combined measurement and laboratory data of a river, determining K as 1.14, which has been widely used with good effect [41].

#### *3.3. Discharge Calculation of River Sections*

It is assumed that the hydraulic conditions in a relatively short reach of a river do not vary obviously, and the average velocity through each section can then be expressed as the mean of the critical velocity of all scattered moving stones. The DSM acquired in the dry season and the identified riverbank watermarks are used to obtain the profile in each section when the flood peak flows through the river channel. The peak discharge through each section is then calculated as

$$
\mathbb{Q} = \upsilon \rtimes A,\tag{9}
$$

where Q is the peak flood discharge through each section, v is the average velocity of the flow through each section, and A is the cross-sectional area when the flood peak flows through culverts.

## *3.4. Performance Evaluation*

Studies have long been conducted to estimate the discharge through culverts on various hydraulic conditions, and flood watermarks have also been used to identify the maximum water level during a specific period [42]. In each river channel, the peak discharge through the downstream culvert (Figure 5) was calculated by Equations (9) and Equation (10) (the Manning formula) to validate the effectiveness of the proposed method and to evaluate the performance of two velocity distributions during the calculation process. The underwater cross-sectional area when the flood peak flows through the culverts is derived from the culvert profile obtained from the DSM and the measured maximum water depth. The Manning formula is normally used to calculate the flow velocity of open channels [43], but its reliability to estimate a peak discharge from the channel cross section has also been confirmed [44]. The general form of the Manning formula is written as

$$w\_c = \frac{k}{n\_c} R\_c^{2/3} f^{1/2},\tag{10}$$

where vc is the average velocity of the flow, Rc is the hydraulic radius, and J is the hydraulic gradient. k is the factor of conversion between the International System of Units (SI) and English units, normally regarded as 1 for SI units. nc is the roughness coefficient of the culvert, determined by the riverbed materials, water bank, and growing plants [45].

The relative accuracy (RA) was calculated to indicate the accuracy of individual calculation results, and the threshold for the relative accuracy of discharge was set to 20% according to empirical studies [46]. Root-mean-square error and mean absolute percentage error (MAPE) were used to evaluate the performance of logarithmic and exponential velocity distribution in the method. These metrics are calculated as:

$$\text{Relative accuracy} = \frac{|Q\_{\bar{i}} - Q\_{\mathcal{W}}|}{Q\_{\mathcal{W}}},\tag{11}$$

$$RMSE = \sqrt{\left(\frac{\sum \left(Q\_{\varepsilon} - Q\_{m}\right)^{2}}{n}\right)} \,\tag{12}$$

$$MAPE = \frac{1}{n} \sum \frac{|Q\_c - Q\_m|}{Q\_m},\tag{13}$$

where Qi is the individual result of the estimated peak discharge through each cross section (*i* = 1, 2, 3, 4), Qe is the estimated peak discharge through each culvert, Qm is the average peak discharge through all calculated cross sections in each river channel, and *n* is the number of calculated cross sections.
