Performance of an Unmanned Aerial Vehicle (UAV) in Calculating the Flood Peak Discharge of Ephemeral Rivers Combined with the Incipient Motion of Moving Stones in Arid Ungauged Regions

Ephemeral rivers are vital to ecosystem balance and human activities as essential surface runoff, while convenient and effective ways of calculating the peak discharge of ephemeral rivers are scarce, especially in ungauged areas. In this study, a new method was proposed using an unmanned aerial vehicle (UAV) combined with the incipient motion of stones to calculate the peak discharge of ephemeral rivers in northwestern China, a typical arid ungauged region. Two field surveys were conducted in dry seasons of 2017 and 2018. Both the logarithmic and the exponential velocity distribution methods were examined when estimating critical initial velocities of moving stones. Results reveal that centimeter-level orthoimages derived from UAV data can demonstrate the movement of stones in the ephemeral river channel throughout one year. Validations with peak discharge through downstream culverts confirmed the effectiveness of the method. The exponential velocity distribution method performs better than the logarithmic method regardless of the amount of water through the two channels. The proposed method performs best in the combination of the exponential method and the river channel with evident flooding (>20 m3/s), with the relative accuracy within 10%. In contrast, in the river channel with a little flow (around 1 m3/s), the accuracies are weak because of the limited number of small moving stones found due to the current resolution of UAV data. The poor performance in the river channel with a little flow could further be improved by identifying smaller moving stones, especially using UAV data with better spatial resolution. The presented method is easy and flexible to apply with appropriate accuracy. It also has great potential for extensive applications in obtaining runoff information of ephemeral rivers in ungauged regions, especially with the quick advance of UAV technology.


Introduction
Water tensions are particularly acute in arid regions globally due to the infrequency and short periods of rainfall. The urbanization process of expanding populations has exacerbated water scarcity method of estimating the flood peak discharge of ephemeral rivers is presented to address the conflict between the shortage of hydrological gauge data and the data needs for research and management of ephemeral rivers.

Study Area
The study area comprised two small watersheds, namely Harulin (45 • 0'22"-45 • 8'52"N, 80 • 53'15"-80 • 59'17"E) and Sailimu (44 • 47'10"-44 • 49'2"N, 80 • 9'15"-81 • 11'26"E), located in the Bortala River Basin of northwestern China (Figure 1). The two watersheds are in the upstream mountainous area of the Bortala River, separately covering an area of 47.9 km 2 (Harulin watershed) and 4.0 km 2 (Sailimu watershed). The climate in Bortala River Basin is a typical temperate continental climate with annual precipitation less than 100 mm and annual evaporation more than 1600 mm. With an arid climate, the primary source of the water system in the basin is seasonal ice and snow meltwater from the top of the mountain and the surface runoff formed by heavy precipitation in summer [32]. Therefore, many streams in the basin are usually ephemeral, and the runoff from April to August often accounts for more than 50% of the annual runoff. Culverts with flood watermarks were only found in the main river channels of Harulin and Sailimu watersheds. Thus, the maximum flood peak flow of ephemeral rivers based on the water level can be calculated in the two watersheds. The main river channels in the two watersheds are separately referred to as river channel H and river channel S.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 23 resolution UAV data, (2) to obtain the flood peak discharge of the ephemeral river, and (3) to compare the performance of the logarithmic and the exponential velocity distribution methods. The proposed method of estimating the flood peak discharge of ephemeral rivers is presented to address the conflict between the shortage of hydrological gauge data and the data needs for research and management of ephemeral rivers.

Study Area
The study area comprised two small watersheds, namely Harulin (45°0'22"-45°8'52"N, 80°53'15"-80°59'17"E) and Sailimu (44°47'10"-44°49'2"N, 80°9'15"-81°11'26"E), located in the Bortala River Basin of northwestern China ( Figure.1). The two watersheds are in the upstream mountainous area of the Bortala River, separately covering an area of 47.9 km 2 (Harulin watershed) and 4.0 km 2 (Sailimu watershed). The climate in Bortala River Basin is a typical temperate continental climate with annual precipitation less than 100 mm and annual evaporation more than 1600 mm. With an arid climate, the primary source of the water system in the basin is seasonal ice and snow meltwater from the top of the mountain and the surface runoff formed by heavy precipitation in summer [32]. Therefore, many streams in the basin are usually ephemeral, and the runoff from April to August often accounts for more than 50% of the annual runoff. Culverts with flood watermarks were only found in the main river channels of Harulin and Sailimu watersheds. Thus, the maximum flood peak flow of ephemeral rivers based on the water level can be calculated in the two watersheds. The main river channels in the two watersheds are separately referred to as river channel H and river channel S.

Methodology
In this study, a new method is presented to estimate the peak discharge of ephemeral rivers, in which high-resolution UAV data was taken as a primary data source, and the method of initiation of motion was carried out. The detailed workflow of the study is as follows (Figure 2): Step 1, high-resolution UAV data acquired in the field were processed to obtain Digital Orthophoto Maps (DOMs) and Digital Surface Models (DSMs) of the studied channels.
Step 2, orthoimages were compared to identify moved stones around selected cross sections, then the critical velocities of these stones were calculated separately using logarithmic and exponential velocity distributions.
Step 3, the average velocity of selected cross sections was calculated from the critical velocities of moving stones, and the discharge through cross sections was then obtained with the extracted cross-sectional area taken from the DSM. Step 4, essential hydraulic variables of culverts were extracted from orthoimages and the DSM to calculate the discharge through a culvert for validation.

Methodology
In this study, a new method is presented to estimate the peak discharge of ephemeral rivers, in which high-resolution UAV data was taken as a primary data source, and the method of initiation of motion was carried out. The detailed workflow of the study is as follows (Figure 2): Step 1, highresolution UAV data acquired in the field were processed to obtain Digital Orthophoto Maps (DOMs) and Digital Surface Models (DSMs) of the studied channels.
Step 2, orthoimages were compared to identify moved stones around selected cross sections, then the critical velocities of these stones were calculated separately using logarithmic and exponential velocity distributions.
Step 3, the average velocity of selected cross sections was calculated from the critical velocities of moving stones, and the discharge through cross sections was then obtained with the extracted cross-sectional area taken from the DSM. Step 4, essential hydraulic variables of culverts were extracted from orthoimages and the DSM to calculate the discharge through a culvert for validation.

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 minutes. 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

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]. 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-off 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.

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)  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 centimeterlevel in our previous study [33]. 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-off 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 rootmean-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.

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.

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.

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 affected 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 F L . The submerged weight W brings the stabilizing condition to the stone. The three forces reach equilibrium and have the relationship [37]: where ∅ is the angle of repose.

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.

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 affected 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 FD 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]: where ∅ is the angle of repose. The submerged weight W, drag FD, and lift FL can be expressed as where U0 is the fluid velocity at the bottom of the channel, CD and CL are respectively the drag and lift coefficients, d is particle nominal diameter, is stone density, and is liquid density. The submerged weight W, drag F D , and lift F L can be expressed as where U 0 is the fluid velocity at the bottom of the channel, C D and C L are respectively the drag and lift coefficients, 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 3 [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.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 23 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 =2.80 g/cm 3 [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.
where ks is the height of the protrusion of the boundary roughness, also known as the side-wall roughness or the riverbed roughness. * is the shear velocity corresponding to the sand resistance and * = . χ is a correction factor related to ks/δ. δ is the calculated thickness of the viscous bottom layer and δ = 11.6 / * . R' is the hydraulic radius corresponding to the sand resistance.
The critical initial shear stress can be expressed as: Sand waves have not yet 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:

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]: where k s 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 * = gR J. χ is a correction factor related to k s /δ. δ is the calculated thickness of the viscous bottom layer and δ = 11.6ν/U * . R' is the hydraulic radius corresponding to the sand resistance. The critical initial shear stress can be expressed as: Sand waves have not yet formed, so let R' = R and U c 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 U c can be derived as: Remote Sens. 2020, 12, 1610 8 of 22 where Θ = f (R * ) ≈ 0.06 under a turbulent condition (R * 500). The density of water ρ = 1 g/cm 3 . g is the acceleration due to gravity and has a value of 9.8 m/s 2 . 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 k s = D.

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 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].

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 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.

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 where v c is the average velocity of the flow, R c 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. n c 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 Remote Sens. 2020, 12, 1610 9 of 22 the performance of logarithmic and exponential velocity distribution in the method. These metrics are calculated as: where Q i is the individual result of the estimated peak discharge through each cross section (i = 1, 2, 3, 4), Q e is the estimated peak discharge through each culvert, Q m is the average peak discharge through all calculated cross sections in each river channel, and n is the number of calculated cross sections.

Stone Movement and Velocity Distribution in the River Channel
Comparing the derived orthoimages in two years, stone movement by water flow was identified in 10 m river reaches over each selected cross section in two river channels ( Figure 6). In the magnified view in the white rectangular area in each river reach over the two years, it is possible to distinguish the movement or non-movement of stones in each river channel. Stones found in the 2017 image but that had disappeared in the 2018 image were identified as moving with the flow throughout the year, while stable stones were regarded as unmoved stones. The nominal particle size of the identified stones was at least 8 cm due to the resolution of the UAV images, and most of these stones can be described as cobbles. The particle sizes of some large boulders in the river section are even more than the average water depth of the selected cross section during the large peak flow event. These boulders were excluded from the circle selection in Figure 6 and the following calculation. Considering that all moving cobbles were mainly moved by water flow, the initial velocities of all identified moving cobbles calculated by the logarithmic or the exponential velocity distribution are presented in Figure 7. In four river reaches of river channel H (Figure 7 a-h), the critical initial velocity of the same moving cobble calculated using the logarithmic method is appreciably larger than the velocity calculated using the exponential method. In river reach A, the critical initial velocities of moving cobbles greatly varied, 1.96-2.65 m/s for the exponential method and 3.40-4.39 m/s for the logarithmic method. In addition, there were more cobbles in the higher velocity grading (shown with the longer arrow) for the logarithmic method (Figure 7a) than for the exponential In both river channel H and S, the movement of cobbles was visible, and cobbles apparently moved by water flow were mainly distributed in the main channel. The main channel was more specific in river channel H with more moved cobbles (8)(9)(10)(11)(12)(13)(14)(15) found in red circles (Figure 6a-h). In river channel S (Figure 6i-p), cobbles of various sizes were found, and the main channel was not obviously washed by flow over one year. Consequently, more unmoved cobbles  in blue circles were identified in river channel S.
Considering that all moving cobbles were mainly moved by water flow, the initial velocities of all identified moving cobbles calculated by the logarithmic or the exponential velocity distribution are presented in Figure 7. In four river reaches of river channel H (Figure 7a-h), the critical initial velocity of the same moving cobble calculated using the logarithmic method is appreciably larger than the velocity calculated using the exponential method. In river reach A, the critical initial velocities of moving cobbles greatly varied, 1.96-2.65 m/s for the exponential method and 3.40-4.39 m/s for the logarithmic method. In addition, there were more cobbles in the higher velocity grading (shown with the longer arrow) for the logarithmic method (Figure 7a) than for the exponential method (Figure 7e). In river reach B, the critical initial velocities of moving cobbles were low, and there were more cobbles with the lowest velocity grading for the logarithmic method (Figure 7b) than for the exponential method (Figure 7f). The channel in river reach C was slowly washed by the water flow and the critical initial velocities of moving cobbles were high (above 4 m/s by the logarithmic method and above 2.4 m/s by the exponential method). There were more cobbles with the highest velocity for the logarithmic method ( Figure 7c) than for the exponential method (Figure 7g), with there being more red arrows in Figure 7c.

Peak Discharge of the River Cross Section
The profiles of each selected cross section at the highest water level are presented in Figure 8. In two river channels, the elevation of the riverbed and riverbank in each cross section continuously decreased downstream. The underwater area of each cross section of river channel H was relatively large, and the section shape was likely rectangular. In contrast, the underwater area of river channel S was relatively small, and the section shape was irregular. Regarding river channel S (Figure 7i-p), velocities obtained by the logarithmic method and exponential method were similar, but the former method produced a slightly broader range of velocities (1.66-2.31 m/s) than the latter method (1.68-1.90 m/s). The logarithmic velocities of moving cobbles in river reach A (Figure 7i) were the lowest among four reaches, all graded at the minimum velocity level, whereas the calculated exponential velocities (Figure 7m) were graded at the average level for the four sections. The critical initial velocities of the moving cobbles in river reach C were similar for the two methods. The condition of river reaches B and D was alike. The maximum velocities calculated using the logarithmic method were concentrated in river reaches B and D (Figure 7j,l) while the maximum and minimum velocities obtained using the exponential method were both in sections B and D (Figure 7n,p).

Peak Discharge of the River Cross Section
The profiles of each selected cross section at the highest water level are presented in Figure 8. In two river channels, the elevation of the riverbed and riverbank in each cross section continuously decreased downstream. The underwater area of each cross section of river channel H was relatively large, and the section shape was likely rectangular. In contrast, the underwater area of river channel S was relatively small, and the section shape was irregular. When the maximum peak flood flowed through cross sections of the two rivers, the maximum water depth and average water depth of the section greatly varied, while the submerged underwater area was similar (Table 2). At the highest peak water level, the average critical initial velocity of all identified moving cobbles in each river reach was regarded as the cross-sectional flow velocity, and the cross-sectional flow rate and peak discharge calculated using the logarithmic and exponential methods were obtained, as shown in Table 2. The cross-section velocities of the river channel H are different using the two methods, with all logarithmic velocities exceeding 3 m/s and all exponential velocities being around 2 m/s. The difference in the cross-sectional velocity in river channel S is not apparent between the two calculation methods, with all velocities being about 2 m/s. Additionally, When the maximum peak flood flowed through cross sections of the two rivers, the maximum water depth and average water depth of the section greatly varied, while the submerged underwater area was similar (Table 2). At the highest peak water level, the average critical initial velocity of all identified moving cobbles in each river reach was regarded as the cross-sectional flow velocity, and the cross-sectional flow rate and peak discharge calculated using the logarithmic and exponential methods were obtained, as shown in Table 2. The cross-section velocities of the river channel H are different using the two methods, with all logarithmic velocities exceeding 3 m/s and all exponential velocities being around 2 m/s. The difference in the cross-sectional velocity in river channel S is not apparent between the two calculation methods, with all velocities being about 2 m/s. Additionally, the peak discharge calculated by the logarithmic method is larger than that by the exponential method in the case of river channel H. For river channel S, the calculated peak discharges are also similar in four cross sections using the logarithmic and exponential method, with the biggest peak discharges both in section C and the smallest both in section D.

Validation of the Estimated River Discharges
The specific profile and key hydraulic variables of each culvert at the time of the maximum peak flood within one year were determined from the topographic information from the DSM and the field measurement ( Figure 9). Culvert H is a wide and shallow box-type culvert with a maximum water depth of 1.25 m when the flood peak passed through. The abundant water in river channel H was almost close to the maximum design discharge of the culvert, while culvert S was square with a smaller water flow. The maximum water depth of culvert S was only 0.3 m. The bases of the two culverts were relatively flat laying a large amount of sediment mixed with sporadically distributed gravel stones. The culvert roughness n c was confirmed as 0.033, a practical value from local hydrological work experience according to conditions of the culvert.
The validation of the peak flood in each river channel is given in Table 3. The flooding in river channel H is evident with a larger peak discharge through culvert H than culvert S, while in river channel S there is little water with larger relative accuracies. The relative accuracies of the logarithmic method and in river channel S all exceed the threshold of 20%. Only the exponential method applied in river channel H indicates an accurate estimation, with the relative accuracies within 10%. The RMSE and MAPE help to further evaluate the performance of different velocity distributions used in the proposed method. RMSE and MAPE are both lower for the exponential method than for the logarithmic method in two river channels, indicating that the former method is more accurate in terms of the incipient motion of moving cobbles regardless of the amount of water in the river channel. depth of 1.25 m when the flood peak passed through. The abundant water in river channel H was almost close to the maximum design discharge of the culvert, while culvert S was square with a smaller water flow. The maximum water depth of culvert S was only 0.3 m. The bases of the two culverts were relatively flat laying a large amount of sediment mixed with sporadically distributed gravel stones. The culvert roughness nc was confirmed as 0.033, a practical value from local hydrological work experience according to conditions of the culvert.  ) and (d) present crosssection profiles and hydraulic parameters used to calculate average velocity and peak discharge through culverts H and S. Ac is the underwater cross-sectional area of the culvert, Xc is the wetted perimeter, J is the hydraulic gradient, Rc is the hydraulic radius, and nc is the roughness coefficient.  ) and (d) present cross-section profiles and hydraulic parameters used to calculate average velocity and peak discharge through culverts H and S. A c is the underwater cross-sectional area of the culvert, X c is the wetted perimeter, J is the hydraulic gradient, R c is the hydraulic radius, and n c is the roughness coefficient.

Value and Extension of the Proposed Method
The present paper provides a convenient, flexible, and reliable method of obtaining the annual peak discharge of an ephemeral river in arid ungauged regions. In empirical studies, statistical methods, hydrological models, and multi-source remote sensing methods are most commonly used to estimate peak discharge (Table 4). In the application of these methods, a variety of ways have been conducted for verification. Most studies used records from nearby gauging stations, such as the water level, the flow rate, and the discharge, to validate discharge estimations [47][48][49][50][51][52][53]. A few studies were verified directly by field observations [16,46,54]. Furthermore, flood marks on the deposits, places, cliffs and other facilities were also commonly used to calculate the flood peak discharge with reliable accuracy in regions lacking flood gauging measurements or to reconstruct historical floods [48,[55][56][57][58]. In arid ungauged regions, flood processes of ephemeral rivers are few and continue to be difficult to directly monitor. Therefore, the identification of clear flood watermarks on construction facilities can be used to verify the existence of flood processes effectively and to validate the maximum peak discharge.  [54] Gauging reports and field data The grid-based rainfall-runoff model (GRM), using the relationship between the runoff coefficient, intensity of rainfall, and curve number and the rational method Multiple satellite altimetry, Landsat series, Sentinel-1/2, and Google Earth Engine (GEE) Using river width and water depth derived from the water surface and water level Obtained high-spatial-resolution images with a UAV Daily precipitation data of the study area were also collected from the nearest meteorological station during the study period to reveal the characteristics of the precipitation in the study area. According to empirical studies, the average rainfall intensity in the study area is around 22.14 mm/h [60]. Generally, it will take more than one hour to generate surface runoff for continuous precipitation with this rainfall intensity in arid regions [61]. In the study area, considerable precipitation was rare, and the top 10 precipitations (>30 mm) during the study period are presented ( Figure 10). Only during concentrated precipitations would the surface flow probably yield and then lead to flooding in the channel. Stones in the ephemeral river channel could have been moved by only one or two flood peaks in September of 2017 or July of 2018. Therefore, the stone movement by flooding in the channel could be used to estimate the discharge of the flood peak.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Daily precipitation data of the study area were also collected from the nearest meteorological station during the study period to reveal the characteristics of the precipitation in the study area. According to empirical studies, the average rainfall intensity in the study area is around 22.14 mm/h [60]. Generally, it will take more than one hour to generate surface runoff for continuous precipitation with this rainfall intensity in arid regions [61]. In the study area, considerable precipitation was rare, and the top 10 precipitations (>30 mm) during the study period are presented ( Figure 10). Only during concentrated precipitations would the surface flow probably yield and then lead to flooding in the channel. Stones in the ephemeral river channel could have been moved by only one or two flood peaks in September of 2017 or July of 2018. Therefore, the stone movement by flooding in the channel could be used to estimate the discharge of the flood peak. The proposed method takes full advantage of the UAV data and the long-term water-free characteristic of ephemeral rivers to estimate peak discharge of ephemeral rivers in ungauged regions. Considering both geographical variables and hydraulic variables, the proposed method is similar to the classic slope-area method, but with the high-resolution UAV data playing a dominant role in the process. On the one hand, the UAV data help to identify terrain information and stone movement. On the other hand, the UAV data are used to derive significant hydraulic variables, such as hydraulic slope, hydraulic radius, and underwater cross-sectional area. The proposed method is vital to the effective long-term monitoring of ephemeral rivers in vast ungauged regions and has great potential to promote at spatial and temporal scales. The primary data required for this method is high-resolution imagery in the dry season. In addition to the UAV data used in this study, other centimeter-level high-resolution images can also be used as data sources. Selecting a typical cross section during the dry season is only needed for in situ measurements, and the water level at the time of the maximum peak flood flow can be determined by comparing the distribution of the soil layer, vegetation, and cobbles. The method established in this paper can also be extended to calculate the maximum peak floods of ephemeral rivers in other typical arid and semi-arid regions with determined dry seasons in these regions.

Performance Evaluation of the Estimated Velocities
Due to the short and variable flowing time of ephemeral rivers during the year, it remains hard to capture the flooding process in time and to obtain field measurements for validation. According to empirical studies, the culvert peak discharge estimated by flood watermarks was regarded as reliable and could be used to evaluate the estimated peak discharge. In addition to the validation of the peak discharge, we also examined the performance of cross-sectional velocity estimations by the incipient motion of stones. The culvert peak discharge was regarded as the real flood peak discharge and was used to calculate cross-sectional velocities (Equation (15)) for evaluating the estimated velocities by logarithmic and exponential velocity distribution methods: The proposed method takes full advantage of the UAV data and the long-term water-free characteristic of ephemeral rivers to estimate peak discharge of ephemeral rivers in ungauged regions. Considering both geographical variables and hydraulic variables, the proposed method is similar to the classic slope-area method, but with the high-resolution UAV data playing a dominant role in the process. On the one hand, the UAV data help to identify terrain information and stone movement. On the other hand, the UAV data are used to derive significant hydraulic variables, such as hydraulic slope, hydraulic radius, and underwater cross-sectional area. The proposed method is vital to the effective long-term monitoring of ephemeral rivers in vast ungauged regions and has great potential to promote at spatial and temporal scales. The primary data required for this method is high-resolution imagery in the dry season. In addition to the UAV data used in this study, other centimeter-level high-resolution images can also be used as data sources. Selecting a typical cross section during the dry season is only needed for in situ measurements, and the water level at the time of the maximum peak flood flow can be determined by comparing the distribution of the soil layer, vegetation, and cobbles. The method established in this paper can also be extended to calculate the maximum peak floods of ephemeral rivers in other typical arid and semi-arid regions with determined dry seasons in these regions.

Performance Evaluation of the Estimated Velocities
Due to the short and variable flowing time of ephemeral rivers during the year, it remains hard to capture the flooding process in time and to obtain field measurements for validation. According to empirical studies, the culvert peak discharge estimated by flood watermarks was regarded as reliable and could be used to evaluate the estimated peak discharge. In addition to the validation of the peak discharge, we also examined the performance of cross-sectional velocity estimations by the incipient motion of stones. The culvert peak discharge was regarded as the real flood peak discharge and was used to calculate cross-sectional velocities (Equation (15)) for evaluating the estimated velocities by logarithmic and exponential velocity distribution methods: where v i represents the cross-sectional velocities calculated by the culvert peak discharge (i = 1,2,3,4,5), Q c is the culvert peak discharge through each river channel, and A i is the cross-sectional area (i = 1,2, 3,4,5). Results indicate that the performance of velocity estimation was best using the exponential method in channel H (Figure 11). The overestimation by the logarithmic method was mostly due to the poor performance of original parameters of large grain sediment [27,40]. In contrast, the exponential velocity equation used in this study was the classic formula proposed by Shamov. The formula is concise and generalized, only using K as a comprehensive coefficient to represent three major forces exerted on the stones and the effects of the irregular shape of stones on the incipient motion. The K value used in this study, also proposed by Shamov, has been verified by many experiments and popularization, and is most suitable for inland arid areas [35]. Furthermore, the river channel S only has a small peak discharge of 0.74 m 3 /s, which could move only small stones during the flood. Due to the current spatial resolution of UAV data, enough small stones moved by water could not be identified. Therefore, the estimated velocities by both methods in channel S significantly deviated from the velocities calculated by culvert peak discharge.
Results indicate that the performance of velocity estimation was best using the exponential method in channel H ( Figure 11). The overestimation by the logarithmic method was mostly due to the poor performance of original parameters of large grain sediment [27,40]. In contrast, the exponential velocity equation used in this study was the classic formula proposed by Shamov. The formula is concise and generalized, only using K as a comprehensive coefficient to represent three major forces exerted on the stones and the effects of the irregular shape of stones on the incipient motion. The K value used in this study, also proposed by Shamov, has been verified by many experiments and popularization, and is most suitable for inland arid areas [35]. Furthermore, the river channel S only has a small peak discharge of 0.74 m 3 /s, which could move only small stones during the flood. Due to the current spatial resolution of UAV data, enough small stones moved by water could not be identified. Therefore, the estimated velocities by both methods in channel S significantly deviated from the velocities calculated by culvert peak discharge. We also added one more cross section (cross section E) to each river channel to evaluate the performance of velocity estimation by the proposed method when smaller stones were used. The identified stone movements of cross section E are presented in Figure 12. The main difference of the new cross section E is that moving stones were much smaller than the other four cross sections, with an average nominal diameter of around 8 cm (Table 5). When smaller stones were identified and used, estimated velocities by the exponential method decreased. In channel S, using smaller stones significantly lead to a decrease of the error of cross-sectional velocity, and the error of cross section E is the minimum in channel S. In channel H, with evident flooding, only using small stones in the calculation could not reveal the real peak flood and led to underestimations. Regarding the logarithmic method, using smaller stones also reduced the error of velocity estimations in cross section E of channel H. Regardless of using smaller stones, the logarithmic method still performed poorly in channel S with only a little flow. We also added one more cross section (cross section E) to each river channel to evaluate the performance of velocity estimation by the proposed method when smaller stones were used. The identified stone movements of cross section E are presented in Figure 12. The main difference of the new cross section E is that moving stones were much smaller than the other four cross sections, with an average nominal diameter of around 8 cm (Table 5). When smaller stones were identified and used, estimated velocities by the exponential method decreased. In channel S, using smaller stones significantly lead to a decrease of the error of cross-sectional velocity, and the error of cross section E is the minimum in channel S. In channel H, with evident flooding, only using small stones in the calculation could not reveal the real peak flood and led to underestimations. Regarding the logarithmic method, using smaller stones also reduced the error of velocity estimations in cross section E of channel H. Regardless of using smaller stones, the logarithmic method still performed poorly in channel S with only a little flow.

The Effects of the Selection of Large Boulders on the Estimation of Peak Discharge
Some of the identified moving stones in the river channel are irregular and appreciably larger than other moving stones in the riverbed. They cannot be entirely submerged during the flood peak events, and major forces exerted on them are different when they start to move [62]. Therefore, whether large boulders are selected in the calculation may affect the estimation of peak discharge and is worth exploring [63,64]. Moved boulders were only identified in cross sections A and D of river channel H, and were found in all four cross sections of river channel S. The results of the critical initial velocity and peak discharge in two river channels using two velocity distribution methods are compared in Figure 13.

The Effects of the Selection of Large Boulders on the Estimation of Peak Discharge
Some of the identified moving stones in the river channel are irregular and appreciably larger than other moving stones in the riverbed. They cannot be entirely submerged during the flood peak events, and major forces exerted on them are different when they start to move [62]. Therefore, whether large boulders are selected in the calculation may affect the estimation of peak discharge and is worth exploring [63,64]. Moved boulders were only identified in cross sections A and D of river channel H, and were found in all four cross sections of river channel S. The results of the critical initial velocity and peak discharge in two river channels using two velocity distribution methods are compared in Figure 13. In the case of river channel H with a large amount of water flow, the selection of large boulders apparently increases the cross-sectional velocity and discharge for both logarithmic and exponential methods (Figure 13a and 13c). The peak discharge in river channel H obtained using the logarithmic method is initially greater than the culvert flow, and the peak discharge thus deviates further from the culvert peak discharge considering the large boulders. In the case of the exponential method, the calculated peak discharge is similar to the culvert peak discharge, irrespective of whether large In the case of river channel H with a large amount of water flow, the selection of large boulders apparently increases the cross-sectional velocity and discharge for both logarithmic and exponential methods (Figure 13a,c). The peak discharge in river channel H obtained using the logarithmic method is initially greater than the culvert flow, and the peak discharge thus deviates further from the culvert peak discharge considering the large boulders. In the case of the exponential method, the calculated peak discharge is similar to the culvert peak discharge, irrespective of whether large boulders are included in the calculation. In the case of river channel S with only a small amount of water, there are many stones with an equal nominal diameter more substantial than the average water depth (around 0.1 m). When large boulders are included, the velocity and peak discharge calculated by the logarithmic method are relatively stable in three cross sections A, C, and D, and only decrease in cross section B with a corresponding decline in the deviation of peak discharge (Figure 13b). Instead, the velocity and peak discharge obtained by the exponential method in four sections of channel S increase, and the deviation from the peak discharge of the culvert further increases (Figure 13d). In all, whether large boulders are considered does affect the discharge estimation, and the effect greatly varies in two river channels with different scales of water flow. With a large amount of water flow, the inclusion of large moved boulders in the calculation significantly increases the estimation of peak discharge in the river channel. Nevertheless, the inclusion of large moved boulders does not significantly affect the results of peak discharge if there is only a small amount of water in the river channel.

Limitations and Uncertainties of the Present Research
The proposed method of calculating the peak discharge of ephemeral rivers using the critical initial velocity still has some weaknesses and uncertainties. First of all, the resolution of the UAV dramatically affects the number of stones identified moving, and the measurement accuracy of the length and width of the cobbles. With the current image resolution, we failed to distinguish enough small moving stones, which leads to less accurate estimations in river channel S with a small discharge. Topographic data and orthoimages with higher resolution could effectively enhance the performance of the method in the future, especially in ephemeral rivers with a small discharge [65]. Secondly, stone properties, such as density and shape, would also generate uncertainties [66][67][68]. In the study area, most of the river stones are irregular granite stones. The constant value of density and the generalized size of stones used in critical initial velocity estimations influenced the accuracies of the peak flood discharge. In addition, only classic incipient motion theories were adopted in this study, while the equation to calculate the critical initial velocity of each moved cobble could be modified with local empirical studies for future extension.

Conclusions
Ephemeral rivers, as an essential part of surface runoff, maintain the stability of the watershed ecosystem and fulfill the water requirements of humans living in arid ungauged regions. With the advantage of high-resolution data and flexible application, the UAV technique has provided more perspectives for runoff monitoring in addition to the traditional hydrological model and satellite remote sensing methods. In this study, a method was proposed to estimate the flood peak discharge of ephemeral rivers using the sediment mobility and UAV data derived in dry seasons. The main conclusions are as follows: 1.
The proposed method performs best in the combination of the exponential method and the river channel with evident flooding (>20 m 3 /s), with relative accuracy within 10%. In the river channel with a little flow (around 1 m 3 /s), the accuracies are weak because of the limited number of small moving stones found due to the current resolution of UAV data.

2.
The exponential velocity distribution method performs better regardless of the amount of water through the two channels, because of the reliable comprehensive coefficient used in the generalized formula.

3.
The effects of using small moving stones or large boulders in the proposed method depend on the discharge in the ephemeral river. In the river with a little flow, identifying smaller moving stones would increase the estimation accuracy. In large ephemeral rivers, estimation results are greatly influenced by using smaller stones or large boulders.
The proposed methodology provides a fast and flexible way of estimating the peak discharge of ephemeral rivers in arid ungauged regions with appropriate accuracy. The method is worthy of further development and local calibration for extensive applications in obtaining runoff information of ephemeral rivers in ungauged regions in the future, especially with the quick advance of UAV technology.