Combining SfM Photogrammetry and Terrestrial Laser Scanning to Assess Event-Scale Sediment Budgets along a Gravel-Bed Ephemeral Stream

Stream power represents the rate of energy expenditure along a stream reach and can be calculated using topographic data acquired via structure-from-motion (SfM) photogrammetry and terrestrial laser scanning (TLS). This study sought to quantitatively relate morphological adjustments in the Azohía Rambla, a gravel-bed ephemeral stream in southeastern Spain, to stream power (ω), critical power (ωc), and energy gradients (∂ω/∂s), along different reference channel reaches of 200 to 300 m in length. High-resolution digital terrain models (HRDTMs), combined with ortophotographs and point clouds from 2018, 2019, and 2020, and ground-based surveys, were used to estimate the spatial variability of morphological sediment budgets and to assess channel bed mobility during the study period at different spatial scales: reference channel reaches (RCRs), pilot bed survey areas (PBSAs), and representative geomorphic units (RGUs). The optimized complementary role of the SfM technique and terrestrial laser scanning allowed the generation of accurate and reliable HRDTMs, upon which a 1-D hydrodynamic model was calibrated and sediment budgets calculated. The resulting high-resolution maps allowed a spatially explicit analysis of stream power and transport efficiency in relation to volumes of erosion and deposition in the RCR and PBSA. In addition, net incision or downcutting and vertical sedimentary accretion were monitored for each flood event in relation to bedforms and hydraulic variables. Sediment sources and sinks and bed armoring processes showed different trends according to the critical energy and stream power gradient, which were verified from field observations. During flows exceeding bankfull discharges (between 18 and 24 m3 s−1 according to channel reach), significant variations in ∂ω/∂s values andω/ωc ratios (e.g., −15 < ∂ω/∂s < 15 Wm−3; ω/ωc > 2 for a peak discharge of 31 m3 s−1) were associated with a large amount of bedload mobilized upstream and vertical accretion along the middle reach (average rise height of 0.20 to 0.35 m for the same event). By contrast, more moderate peak flows (≤10 m3 s−1) only produced minor changes resulting in surface washing, selective transport, and local bed scouring.


Introduction
Ephemeral streams are dry, irregular, and often unpredictable watercourses, particularly sensitive to short-term climatic changes and to human impacts. Consequently, their response can be highly variable, ranging from slight river-bed adjustments [1] to overall changes in channel geometry after flash floods [2][3][4]. This type of stream is found in large areas around the world, in arid and semi-arid environments, where quite fragile hydromorphological systems often form, subjected to a metastable dynamic equilibrium that strongly depends on the frequency and magnitude of floods. During periods of relatively more frequent and lower magnitude peak flows, surface bed washing and downcutting phenomena can produce secondary channels, lateral erosion, and partial destruction of low bars [5]. Infrequent and large floods affect the entire main channel, modifying its global geometry, and create prominent features of aggradation and incision, depending on the bed materials [2].
Although ephemeral streams tend to maintain similar flow characteristics over long periods of time [6], there are semi-arid Mediterranean areas whose dry and intermittent streams are being subjected to the effects of the current climate change and have begun to show signs of relatively rapid geomorphic adjustment. A clear example of this is the coastal drainage system in southeastern Spain, mainly composed of ephemeral gravel-bed streams. These ephemeral channels usually transport a high tractive load from heavily weathered headwater areas and alluvial channel banks, producing an abundant amount of coarse sediment (gravel and pebbles). In general, they are channels with a high width-to-depth ratio, which, except in sections of structural geological control, run through poorly consolidated alluvial formations, showing great dynamism and instability. Their upper reaches dissect steep alluvial fans, giving rise to entrenched and deep channels subjected to downcutting and basal erosion, which cause frequent breaks and collapses of material. In addition, downstream, these gravel-bed streams experience channel widening, aggradation, and bed armoring [7]. As a result, the streambed often is very permeable, so flood peak discharges decrease sharply downstream due to transmission losses [8,9]. Due to high infiltration rates, incision, and transport capacity decrease, while vegetation tends to grow locally on stable bars and subsurface moisture storage increases, promoting positive feedback for instream aggradation [10][11][12]. Downstream, alluvium within the channel increases in thickness and smooths the bed slope, creating a channel profile that is straighter than concave, unlike in perennial rivers [13,14].
The scientific literature on spatial morphological variability in this type of ephemeral streams in relation to erosion volumes, deposition, and sediment transport (sediment budget) on a detailed scale is scarce. Recently, high-resolution techniques, such as structure-from-motion (SfM) photogrammetry and terrestrial laser scanning (TLS), have been tested to evaluate morphological adjustments in ephemeral streams [15][16][17][18][19].
In this paper, we propose an approach to assess the event-scale sediment budget along an ephemeral gravel-bed channel that couples the results of the structure-from-motion multi-view stereo (SfM-MVS) and TLS techniques. The use of SfM-MVS in a low-altitude unmanned aerial vehicle (UAV) allows the development of orthomosaics and high-resolution digital elevation models (<5 cm per pixel) [20,21]. This survey product offers significantly better resolution and much lower processing costs than others generated by traditional imaging platforms, such as satellites or aircraft [22][23][24][25]. In geomorphological studies, UAV technology has already been applied to monitor fluvial system changes and dynamics [26][27][28]. The high temporal resolution offered by UAV allows assessment of the rapid evolution of river ecosystems. The UAV-SfM technique, when compared to traditional photogrammetry and other topographic surveying approaches (e.g., total station and DGPS), has proven to be particularly effective for the detection of channel morphological changes [29][30][31].
Laser Imaging Detection and Ranging (LiDAR) is a land surveying technique that accurately measures and collects millions of points in a few minutes, providing geometrical information from a scene by precisely measuring the time-of-flight of laser beams. A laser scanner device can be mounted on moving vehicles, such as UAVs or ground-based vehicles (for TLS). The 3-D point cloud (3DPC) datasets obtained from these techniques have usually exploited geometrical features' identification and extraction, providing unlimited possibilities in geomorphology and, more specifically, in fluvial geomorphology [20,[32][33][34].
In order to estimate the morphological sediment budgets in the study of the gravel-bed stream, three analyses were performed using SfM-MVS and TLS, as follows: (1) An analysis of hydraulic variables at the cell scale was performed for the flash floods under study, based on rasterized information from cross-section data and using 1-D hydrodynamic modeling and high-resolution digital terrain models (HRDTMs); (2) erosion and deposition values and net sediment budgets were extracted from datasets of SfM and TLS for periods with different numbers and magnitudes of events; (3) the 3DPC datasets were used to calculate a set of statistical parameters, regarding elevation and volume differences after each event in reference channel reaches (RCRs), pilot bed survey areas (PBSAs), and representative geomorphic units (RGUs).

Study Area and Field Surveys
The SfM and TLS surveys were carried out in the Rambla de la Azohía, located in southeastern Spain (Region of Murcia) ( Figure 1). The area drained by this ephemeral stream and its tributaries is relatively small (13 km 2 ), to the point that it could be completely covered by a single storm. The climate of the study area presents strong seasonal contrasts, with extreme droughts and monthly potential evapotranspiration values higher than 140 mm in the summer season. These conditions, together with a scarce vegetation cover and very irregular and intense rains, with records above 50 mm/h, often cause large flash floods [35]. This watercourse is part of a littoral hydrological system made up of numerous ephemeral gravel-bed channels, draining directly into the Mediterranean Sea. Steep slopes modeled on metamorphic poorly permeable materials (mainly phyllites, schists, and quartzites), subjected to intense weathering, constitute an important source of sediments and ensure the supply of coarse particles to the channel when heavy rains occur. As a result of this, the Rambla de la Azohía exhibits a particularly active morphodynamics, ranging from minor bed alterations to large morphological channel changes. Instead of using a monthly or annual time scale, singular extreme hydrological events, isolated in time, were analyzed to determine overall and local morphological channel changes and sediment budgets at different spatial scales: RCR, PBSA, and RGU (Figures 1 and 2). Specifically, morphological changes were monitored in three stages: (1)   Chronogram of the flow events and of the field surveys (FSs) performed using either UAV-SfM (Structure from Motion based on information collected with UAVs), Terrestrial Laser Scanning (TLS), or the combination of the two techniques (UAV-SfM + TLS). P = precipitation (mm); Qp-MDR = peak discharge (m 3 s −1 ) at the Middle Reach; Qp-UPR = peak discharge (m 3 s −1 ) at the Upper Reach. In parentheses, the maximum rainfall intensity (mm h −1 ).
The first two flow events of the analysis period (18 November 2018 and 19-20 April 2019) were monitored by UAV and TLS, covering survey areas with different extensions (Figure 2). The UAV-SfM of September 2018 covered a wide area (>600 m in length) along the upper, middle, and lower sections of the study stream, while the TLS only comprised the RCR, upper (UPR) and middle (MDR) shortened reaches (<200 m long), and their respective PBSA within them. After the aforementioned peak flows, a date and common study area were chosen to carry out both SfM-MVS and TLS, in order to obtain the same georeferencing, compare the associated errors, and analyze the possibility of the complementary optimal use in time of both techniques. In fact, the morphological changes produced after the peak flows of 12-13 September and 2 December 2019 were only detected by TLS ( Figure 2).

Material and Methods
The affected areas, elevation differences, volumes of erosion and deposition, and sediment budgets were obtained for the two RCR (upper (UPR) and middle (MDR) stretches), using data derived from UAV-SfM and TLS. To detect changes in bedload budgets at the bed level exclusively, an analysis was performed using TLS in pilot survey bed areas (PSBAs) with higher point density. The RCRs were delineated using the limits of the water sheet generated by the hydraulic model for each flood event. In addition, the values of the total area of interest (m 2 ), total volume difference average (m 3 ), net thickness difference (m) for the area of interest, and percent imbalance (departure from equilibrium), among other variables, were provided by both techniques for each RGU.

SfM-MVS Photogrammetry
The remote information was collected with a Phantom 4 Pro, with a 20 Mpx camera and a 1-inch sensor. The photogrammetric parameters of the flights were set before the development of the first field survey and the inspection of the study area, using the software tool DJI GS Pro© with equivalent values for all the surveys; this procedure ensured optimum comparability among the digital elevation models generated at different sites and on different dates. The average flight height was c. 50 m above ground level, resulting in images with a resolution of c. 1 cm.
Standard targets extracted from Agisoft PhotoScan Pro 1.2.2© (Agisoft, St. Petersburg, Russia) were used as ground control points (GCPs) for georeferencing and analyzing the quality of the reconstructed HRDTM. These targets were randomly distributed throughout each of the stream segments and were surveyed using a Prexiso G5© (Leica, Wetzlar, Germany) GPS-RTK instrument. The GPS-RTK instrument was connected to the regional network of GPS GNSS (Network of Reference Stations in the Autonomous Community of Murcia "Meristemum") via mobile signal, to obtain differential corrections in real time. In each survey, c. 66% of the measured points were subsequently integrated into the model for georeferencing purposes and the remaining 34% were used as checkpoints for the validation of the HRDTM. Permanent GCPs were also installed, using FENO survey markers in each of the GRSs to define common reference points for all surveys and methodologies. All the GPS points were recorded in the WGS84 global reference system.
We selected images with a high level of overlap (80-90%) to ensure successful subsequent image matching [36]. These images and the coordinates of the GCPs' were uploaded into structure-from-motion photogrammetry software, Agisoft PhotoScan Pro v.1.2.2© (Agisoft, Russia), and aligned using high-precision settings and generic pair preselection to create a sparse point cloud [37,38]. We optimized the image alignment using the GCPs' data and regenerated the model using the same settings as previously; consequently, the HRDTM, orthophoto, and dense point cloud were exported for further analysis in a global georeferencing system (WGS84).
The difference between the UAVs-collected HRDTM (DoD) in different field surveys, as well as between the TLS-generated HRDTM, was calculated in ArcGis 10.5 © (ESRI, Redlands, CA, USA), by subtracting the final topography from the previous topography for the same area [17]. The same polygon was used for the comparison of the complete RCR, and similarly for the RGU, to maintain the boundaries for change detection, and the HRDTMs were created with equal pixel size (0.02 × 0.02 m), starting from a specific point, so that the grid cells were fully coincident. To assess the sediment budgets and morphological bed changes, we used statistical data on differences in area, volume, height, and percent imbalance for each of the RGUs. In addition, histograms of the volume and bar graphs were used to represent the trends of elevations.

Terrestrial Laser Scanning (TLS)
In this study, a terrestrial 3-D laser scanner (Leica ScanStation C10 model) was used. This instrument captures the coordinates of points up to 200 m away and the returned energy, and estimates the RBG color using a digital camera. The precision of the instrument is 2 mm. The TLS captures the evenly spaced points when the surface is vertical, but when it is sub-horizontal, the scanned points are sparse. Consequently, to scan the bedforms, it is necessary to choose elevated positions for the scanner, when possible. Another option is to perform multiple scanning stations and to register them, merging sparse point clouds and reducing the point spacing and the shadow areas.
Several scans were made from different stations placed strategically along the bedforms for detailed geometrical definition using TLS. For each date, several scanning stations were made and registered using HDS targets. In addition, the FENO survey markers were also scanned in order to register the TLS point clouds described in the previous section and the SfM datasets. The scan performed on 5 September 2019 was used as the benchmark for all the TLS scans. On that date, the SfM and TLS data were acquired. The coordinates of the FENO survey markers were extracted from that SfM dataset and the TLS dataset was registered using those coordinates. Consequently, the 2018 and 2020 scans were registered to that benchmark using existing buildings as the reference. The registration was performed using the ICP plugin of CloudCompare [39] and was applied to existing buildings, which were considered as stable areas for the whole study period. Then, the registration was evaluated using the M3C2 plugin of CloudCompare. Both registrations had a mean error of 2 mm.
To study the elevation differences and volume changes, two pilot bed survey areas (PBSAs) with higher point density, within the UPR and MDR reference channel reaches, were selected: the upper area (8.55 × 17.7 m 2 ), and the middle area (16.3 × 27.55 m 2 ). The existing vegetation (scrub, bushes, and trees) was removed using the Canupo (CAractérisation de NUages de POints) CloudCompare plugin with specific trainings for both areas [40] (Figure 3). Then, each area was rasterized using a 5-cm grid and empty areas were interpolated.

Criteria for the Selection of Reference Channel Reaches (RCRs), Pilot Bed Survey Areas (PBSAs), and Representative Geomorphic Units (RGUs)
Reference upper and middle reaches (upper and middle RCR) were chosen depending on the accessibility and the presence of potentially active surface areas. In addition, these stream stretches better represent the variations in the sediment budget and morphological adjustments within the main channel, as opposed to the lower reach, whose global changes are the product of sediment-laden flood water spills, affecting crevasse splays and a wide active flood bed. In fact, the choice of the upper and middle RCR in this ephemeral gravel-bed stream is intended to reflect the clear influence that coarse grain-sizes and bedload, as the dominant transport modality, have on the morphological sediment budget within the channel itself [41]. Downstream, along the lower reach to the mouth of this stream, a distributary drainage system exists, which implies that the budget of more complex sediments needs to be evaluated, accounting for mixed sources (hillslopes and secondary channels), for the temporary storage of sediment in bars and alluvial plains [42].
The upper RCR is an entrenched channel stretch that enabled one scanning station at an elevated position. The rest of the scans were performed on the ground close to the bedforms. In the middle RCR, where bedforms arises and the vegetation cover is scarce, only scanning stations on the ground were performed, reducing the density of the points. The choice of PBSA is intended to detect and quantify precise changes in bed sediment budgets and bed load mobility, while the RGU monitoring allows us to analyze behavior patterns and spatial variability of forms associated with scouring and deposition areas within the channel.

Hydrometeorological and Hydraulic Data
A rainfall-runoff model was used to simulate runoff on the basis of rainfall events during the budget period. Five pre-simulated runoff events were defined on the basis of the rainfall amount, initial wetness conditions, and initial discharge and were subsequently adjusted according to the real flow heights measured by water pressure sensors. The HEC-HMS 3.5 program [43] was combined with a DTM in grid format and 4 × 4 m cell resolution [44] in order to generate the drainage patterns and catchment limits. Then, in accordance with [45], we assembled the headwater sub-basins that are homogeneous in texture and roughness and individualized those affected by confluences with important breaks in the bed slope. The main abstraction used was the infiltration curve number (NC) of each sub-basin, following the methodology of the Soil Conservation Service of the US Department of Agriculture [46]. The geometries and data tables, referring to the vector lithological layer of the Geological and Mining Institute of Spain (IGME) at a scale of 1: 50,000, and the land use map of the Corine Land Cover (CLC), Level 3 (2012) were combined using GIS to estimate this parameter, according to the different hydrological soil groups and substrate textures. Using the antecedent conditions of dry soil and pluviometric data for intervals of 5 min at the Cuesta del Cedacero station, we calculated the dimensionless unit hydrograph of the SCS and its peak unit flow (q p ) (m 3 s −1 ) (Equation (1)): where A is the area of the basin (km 2 ) and Tp the time to peak (h), estimated according to Equation (2): where D is the duration of the effective rainfall (h) and tp is the lag time, in our case, 0.6 times the concentration time of the catchment [47]. The output data ( Figure 4) were then transferred to the 1-D hydrodynamic model HEC-RAS [48], supported by HRDTM, to obtain the channel profile and cross-sections and hydraulic variables. A total of 212 cross-sections were used, with an average separation of 2.0-3.5 between them in both channel reaches. The theoretical hydrometeorological discharge, estimated for the conditions of a subcritical flow regime, was calibrated using that from the real water height data. Variations in the water level during each flood event were recorded with levelogger sensors and corrected using barometric calibration devices (barologger sensors), located in a station intermediate between the upper and middle reaches. For all the cross-sections, information concerning the flow surface area, hydraulic radius, velocity, shear stress, and stream power was obtained. The power per unit length of stream (Ω) and mean stream power (ω) at peak flood discharge (Qp) were calculated, according to Equations (3) and (4), respectively [49], for each cell i between each pair of cross-sections: where γ is the specific weight of water (Nm −3 ); Q is the discharge (m 3 s −1 ); S w is the water surface slope (m m −1 ), which is used to estimate the energy gradient; w is the water-surface width (m); Ω represents the energy dissipation per unit channel length (Wm −1 ); and ω the energy expenditure per unit bed area (Wm −2 ). The mean stream power gradient (∂ω/∂s) was calculated by subtracting ω in cross-sectional cell i from ω in cell i−1 and dividing the difference by the distance between the centroid of each pair of consecutive cells along the channel centerline. In this way, positive ∂ω/∂s values indicate downstream increases in ω, while a negative ∂ω/∂s indicates that ω is decreasing from one cell to the next [50]. The energy expenditure beyond the critical mean stream power (ω c ) and the stream power (ω)/resisting power (ω c ) ratios (ω/ω c ) were estimated in each cross-sectional cell to determine the spatial variability of the energy available for sediment transport and explain morphological changes during peak flows.
The ω c was obtained using Equations (10) and (16) from Parker et al. [51], which relate ω c to slope and grain size. A representative reach pebble count was made in the RCR along transect stretches from one stream bank to the other, while the RGU were sampled by volumetric extraction of bed-material to differentiate surface and subsurface particle-size distributions. The median grain size (D 50 ) and 84th percentile (D 84 ) [49] were calculated for each transect and set of cross-sections.

Sediment Budget Calculation and Detection of RGU Adjustments
The sediment budgets in an EBGS accounts for the rates and processes of erosion and bedload transport within the channel, and for the temporary storage of coarse grain-size material in bars and pools along bends and runs. Significant rates of weathering and breakdown of sediment, during transport or storage, are also frequent in semi-arid environments, such as this one [42]. Variations in bedload transport downstream lead to adjustments in the erosion and deposition volumes, affecting the sediment budget, and vice versa, changes in this budget could condition the amount of material entrained [52]. This premise makes the estimation of sediment budgets a useful and necessary task for the determination of the magnitude of morphological bed changes. An analysis relating event-scale bedload transport to changes in channel morphology is generally difficult, because this requires a detailed topographic survey of the channel immediately before and after each flow event [53]. This consideration is particularly relevant in ephemeral streams, where the channel remains dry most of the time and only carries water in flash floods that last a few hours. Consequently, the monitoring of morphological bed adjustments and the calculation of sediment budgets at the event scale can provide more realistic results than the analysis of long study periods. In our case, a strategy based on the combined use of high-definition topographic products was adopted to reveal the geomorphological effects of isolated events or, at most, two events of different magnitudes.
High-density 3DPC and HRDTM, generated from SfM-MVS and TLS, were overlain to obtain good accuracy in the sediment budget calculation. At the event level, this evaluation went further, because it allowed us to relate the net bed variations to the characteristics of the hydrograph in question. The errors associated with the use of two surface models to determine volumetric sediment budgets were described and assumed for each comparative survey analysis, according to Brasington et al. [54]. The same procedure was applied to detect morphological adjustments attributed to erosion or deposition in the RGU ( Figure 5), overlaying the 3DPC on the raster image derived from a vector layer that contained the polygons of these geomorphic units within the wetted channel perimeter in each flow event.

Stream Power Maps
The mean stream power (ω), energy gradient (∂ω/∂s), and mean stream power (ω)/resisting power (ω c ) ratio (ω/ω c ) were estimated discreetly along each of the RCR for the most significant flood events: 19-20 April 2019 and 12-13 September 2019. Figure 6 shows the grids of these variables for the UPR and MDR during these events. The resulting maps were then related to the spatial patterns of sediment budgets in order to explain better the processes that control the morphological channel adjustments.  Figure 6 for these variables encompass 95% of the data; that is, 338.3 Wm −2 and 37.7 Wm −3 , respectively). Note, however, that more than 75% of the set of values of ∂ω/∂s in this event fell within the range 15 to −15 Wm −3 and that the extreme data reflect sharp local changes in ω. As in the case of the flow velocity and the shear stress, the highest values of ω were concentrated in entrenched channel cross-sections, along straight to slightly sinuous stretches, and locally over steep riffles characterized by high near-bed velocity gradients. This pattern is linked to incision processes in several Mediterranean ephemeral streams [35,55,56]. However, we observed a significantly different behavior in the gravel beds studied here, where important transitory erosion was often accompanied by extraordinary bedload mobility and high deposition rates, which tended to cause channel aggradation.
The distribution of the ω/ω c ratios was also skewed for both RCR and the flood events, with maximum values very distant from the mean and median. In the events of April and September, maximum ratios around 10 and 7 were reached, respectively, in the two reference stretches, with much lower and different medians depending on the section. Specifically, during the peak flood of 19-20 April 2019, the ω/ω c medians ranged from 2.7 in RCR (UPR) to 1.2 in RCR (MDR).

Sediment Budgets in RCR Determined from SfMData
Gravel-bed erosion and aggradation are recurrent geomorphological processes that discretely affect the global evolution of ephemeral gravel-bed channels over time. Changes at the event scale are usually very disparate; in some cases, scour and downcutting phenomena predominate, and in others deposition. The present study contributes to a better understanding of the current behavior of this type of channel and the subsequent trend of its minor adjustments, but it contributes little in relation to possible long-term morphodynamic changes. A good indicator of these adjustments is the sediment budget in the reaches with the highest bedload. In particular, in the RCR analyzed here, the bedload experienced important spatial variations at the event scale, significantly affecting the sediment sources (areas of erosion) and sinks (areas of deposition) in the short period of analysis. In the September 2018-September 2019 stage, lateral erosion from steep alluvial banks, active low bars, partially destroyed coarse bar heads, and finer-grained bar tails provided a large bedload in the downstream direction. As a result, the greatest deposition thicknesses were recorded in the flanks of the longitudinal and medial alluvial bars, in both RCR (see Figure 7, drawn up on the basis of SfM-MVS data).
Previous field observations suggest that during low-water stages, as in the event of November 2018, when the top of the bar emerged, vertical accretion of these bars ceases and new secondary channels form, causing small island bars to migrate. The flash flood of 19-20 April 2019 resumed the aggradation process, with very widespread increases in bed height. The SfM statistical data in Table 1 show significant cumulative changes in ground surface elevation after the flow events of 19 November 2018 and 19-20 April 2019. An average net thickness difference of around 22 and 21 cm was found in the upper and middle reaches, respectively. The sedimentary balance was positive in both RCR, verifying a clear dominance of deposition over erosion. The proximity of both reaches to abundant sources of coarse sediment, mainly due to locally strong connectivity between the channel bed and active alluvial banks, together with a similar bed slope, promoting a notable mobility of the bedload, explains the scarce variation between the RCR in the average sediment budget. The high average deposition rate observed in the middle reach (unit TVSL = 0.218 m 3 m −2 ) was not accompanied by equally significant unit volumes of erosion upstream (unit TVSR = 0.128 m 3 m −2 ). This suggests that the sources of sediments include not only the upper RCR but also other more entrenched sections close to the headwater area and intermediate sections between the two RCR, where bank breaking and gravel bar removal are especially active processes. In addition, the minor unit net erosion recorded here during this event decreased even more downstream, reaching 0.086 m 3 m −2 in the middle reach (Table 1), so the bed cutting could not progress in a downstream direction.   Heterogeneities due to differences in grain sizes and bedforms may create substantial velocity and shear stress variations across the channel or downstream during a single discharge [57]. This can explain the patchwork of ω values in both reaches (Figure 7) and, thus, the spatial changes in sediment budgets displayed in Table 1. The upper RCR registered the highest average unit volume of surface raising and lowering (0.13 and 0.23 m 3 m −2 , respectively), in line with the highest mean values of ω (72.4 Wm −2 ), but also the highest SD in the net thickness differences and in the ∂ω/∂s gradients (16.4 Wm −3 ).

Sediment Budget and Stream Power in PBSA Estimated from TLS Data
The results in Table 1 (SfM-based results) differ from those obtained using TLS for the same reference channel reach during the period of December 2018 to September 2019, for two reasons. Firstly, there were differences caused by the fact that the UAV-SfM surveys included an event (the peak flow of November 18, 2018) that was not monitored with TLS. Secondly, it is widely accepted that the evaluation of sediment budgets is sensitive to differences in DTM quality and in the process used to suppress vegetation cover, since they inherently incorporate errors into the generated terrain models [54]. Although this type of error was considered small in the two reference reaches, given the high resolution of the DTMs used, this factor cannot be completely discounted in the more vegetated areas.
In addition, the high-density 3-D point clouds generated from TLS were applied to pilot bed areas (PBSAs) in order to detect detailed changes in bedload budgets. In Figures 8 and 9 and Table 2     The analysis of budgets expressed in unit volume for 10-cm elevation intervals show that the largest volumetric variations implied morphological adjustments, ranging between −0.2 and +0.2 m in bed elevation (Table 2). Within these intervals, the greatest unit volume of deposition (41.3 dm 3 m −2 ) corresponded to PBSA (MDR) after the April 2019 event, and the greatest erosion volume (48.9 dm 3 m −2 ) to PBSA (UPR) after the two minor events of the period September 2019 to January 2020. These data are consistent with those displayed for both cases in the normalized histograms of elevation differences (Figure 9), where the interval from 0 to +0.2 m had a cumulative frequency greater than 60% and that from 0 to −0.2 m exceeded 70%.
Likewise, the variations and signs of the mean energy gradient faithfully reflected the different patterns observed in the bedload budgets of each PBSA. Assuming that ∂ω/∂s represents the clearest expression of sediment change in the nearest cell downstream, we can hypothesize that this hydraulic variable had a strong influence on the net erosion or deposition, and on the total mean bedload volume per unit area. However, while the net deposition pattern expected with a higher positive ∂ω/∂s did occur, the same did not happen with the net erosion pattern, since the most negative ∂ω/∂s registered in PBSA (UPR), during the flood peak on 19 April 2019, finally translated into a unit net volume of downcutting lower than the one accumulated by the minor events of September 2019 and January 2020. The larger flow event in April generated significantly higher mean ∂ω/∂s valuesthan the minor peak flow in September, with a negative sign in PBSA (UPR) and a positive one in PBSA (MDR), which resulted in a decrease from −3.8 to −0.4 Wm −3 for the first case and from 1.9 to 1.6 Wm −3 for the second (Table 3). This suggests that in the major event, a much larger bed load was mobilized and replaced than in the minor event, and that during the latter the energy expenditure was used for scouring and there was not enough surplus (ω/ω c = 1.6vs. 2.6 in the April flood) to drag and deposit the material transported from upstream.  Such results corroborate the hypothesis that in this ephemetal gravel-bed stream, the events of greater magnitude (Qp ≥30 m 3 s −1 ) produce vertical sedimentary accretion on the bed [41], exceeding even the critical stress of the coarsest particles (gravel and pebbles) in PBSA (UPR) (mean ω/ω c >2) and thus moving a large bedload downstream. In contrast, more moderate peak flow rates (Qp ≤10 m 3 s −1 ) were only able to exert a surface washing action, selective transport, and local scouring. In both cases, all the events analyzed tended to reinforce and arm the bed, as seen in situ, making it more resistant to erosion in future floods.

Morphological Bed Adjustments Observed in RCR from TLS Data
The detailed topographic survey performed for each channel reach prior to and following a bedload-transporting flow event made it possible to detect precise changes in channel morphology in relation to event-scale bedload fluxes. Such an approach was highly recommended by Kasprak et al. [53], when an HRDTM is available, and by Singer and Michaelides [58] for the case of ephemeral channels, where topography appears to remain unchanged.
In Figure 10, the morphological changes observed along the channel center line of both reference reaches using TLS datasets. From this, different patterns of spatial variability were inferred, depending on the reference site, the magnitude of the flood, and the textural and topographic characteristics immediately prior to the occurrence of each event.
In the upper RCR, the major morphological bed changes corresponded to two types of channel narrowing: at the exit of a meander and that due to structural control in less sinuous reaches. In the first case, the forms of excavation and vertical accretion reflected a marked bed alteration just 15 m in length, the incision prevailing in the largest event and the deposition in the least torrential (Figure 10c in UPR). In the second case, downstream narrowing, the major event caused significant variations in the topography of the bed, creating new active bars that later, in the event of September 2019, were partially or totally destroyed. These morphological adjustments must have been caused by high stream power values concentrated in short distances, as already suggested by Conesa-García et al. [59] for entrenched channel cross-sections, along stretches with sudden changes in bed roughness. This pattern has been found to be linked to channel degradation in semi-arid ephemeral streams (e.g., [35,56]). However, we observed significantly different behavior in the gravel beds studied here, where important transitory erosion was often followed by considerable replenishment of coarse-grained sediments, finally producing bed aggradation. The rest of the stretch showed slight changes in the inter-bends or runs transition zones and significant homogeneous variation in the final part, where progressive erosion downstream caused by the flash flood of April 2019 was partially offset by aggradation processes during the following minor event.
The middle RCR presents two spatial patterns of geomorphic change: one of them around an extensive central bar (65 m long), characterized by important variations of different sign, depending on whether it is the anterior sector, the bar itself, or the posterior area in a downstream direction. At the back of the bar head, an incision of 0 to −0.5 m was observed, while the head and platform of the bar experienced vertical accretion in both events, and the tail suffered a lowering in height of up to −0.3 m. The minor event in September introduced hardly any variations in this reach except in the final part of the platform, which was lowered by around −0.4 m (Figure 10c in MDR). The second pattern, represented in the rest of the reference channel reaches, arises from a certain bed stability, only interrupted by slight alternative adjustments of incision and vertical accretion, resulting in net elevation differences between +0.2 and −0.2 m.

Sediment Budgets for Different RGUs in Relation to Stream Power Data and Field Surveys
In Table 4, standard variables commonly used in sediment budgeting (sediment area and volume) for each representative geomorphic unit (RGU) and RCR are shown, derived from SfM and TLS data, concerning the peak flows produced during their respective field surveys. As mentioned above, the SfM monitoring periods included one more event than those of the TLS; therefore, despite offering fairly similar surface area data for each RGU, they provided different volume data and errors. Although the time interval considered (September 2018 to September 2019) is relatively short, and only covers one or two events, depending on the type of survey performed, the fact that almost all RGUs registered a positive total average balance in favor of deposition is striking (Table 4 and Figure 11). Only the inter-bar active bed in the upper RCR scored a slightly negative budget, according to the TLS data, after the flash flood of 19-20 April 2019. This result, added to the fact that the difference in average total thickness or the average unit volume difference for each RGU was always lower with TLS than with SfM, undoubtedly shows that the November 2018 peak flow, although lower than that of April 2019, also contributed to the deposition. In summary, it can be considered that almost all the bed forms experienced over-sedimentation and, therefore, vertical accretion. The maximum accumulation produced after both events (0.18 to 0.23 m 3 /m 2 ) took place on the active bed zones and on the scarcely vegetated gravel bars, which behaved as somewhat unstable forms (Table 4 and Figure 11). Even the main active bed in run stretches within the middle reach rose. The positive sediment budget in lateral tail and high bars could be misleading (average total thickness difference around 0.20 m), but, nevertheless, corresponds to local collapses of material caused by basal undermining in high detrital banks (Figure 12(1,2)) and lateral accretion by the adjoining of finer sediments (sand and gravel) in alternate and central bars (Figure 12(4)). However, this contribution was not sufficient to explain the large imbalance between the average volume of sediment returned to the channel through bank collapse (equivalent to 11% of the total flux) and the net volume of sand and gravel deposited on the bars (66% of the total budget). Consequently, most of the bedload mobilized and deposited in both stream stretches could come from upper sediment sources, associated with intense downcutting processes in deep gullies that dissect the apical zone of thick alluvial fans. López-Bermúdez et al. [60], Harvey [61], and Aguilar et al. [62] found similar results when studying steep alluvial fans dissected by gullies in semi-arid and arid environments. Simultaneously, low-order gravel bed tributaries at the headwaters currently develop high geomorphic activity, leading to considerable variations in downstream sediment yield, which is consistent with the pattern described by Lisle et al. [63] and Yuill et al. [64] for this type of ephemeral gravel-bed channel.
By comparing the sediment budgets obtained using TLS after the flash flood of April 2019 with the stream power indicators referring to this event (see Figure 6), we observed that the ω/ω c ratio exerted a greater influence than ∂ω/∂s on the incision and vertical accretion processes, depending on the predominant type of RGU in each set of cells. Specifically, in sites with dominant changing RGU, such as active channel inter-bars or runs (a), low, active non-vegetated bars (b), and low, scarcely vegetated bars (c), strong relationships were found between the average total thickness difference and the ω/ω c ratio (r 2 ≈ 0.89 in both channel reaches). The relationship of surface height change with ∂ω/∂s was also statistically significant but not as strong as that with the average relative excess energy (r 2 = 0.83 and 0.66 in UPR and MDR respectively, p-value < 0.05 for the 95%confidence interval). In any case, a different behavior could be intuited in the most unstable bars of each of the RCR in relation to the ratios ∂ω/∂s and ω/ω c estimated for each channel stretch. In the upper RCR cell sets, where these RGUs predominated, average values of ∂ω/∂s and ω/ω c of −4.4 Wm −3 and 2.3, respectively, were found, which led to a slight lowering in the active bed (−0.024 m on overage) and a general increase in height of only 5 cm in longitudinal non-vegetated bars (Figure 12(4)). In contrast, smaller values of ∂ω/∂s and ω/ω c (−1.6 Wm −3 and 1.5, respectively) at the front and head of a mid-channel bar (central island), crossing the middle RCR, caused a thicker deposit (average bed rise of 11 to 15 cm) and greater morphological adjustments (Figure 12(5)).
Field observations immediately after this event showed that such sediment budgets were accompanied by a relocation of the bed material and changes in bed texture ( Figure 12). Along all the upper RCR, riffles proved to be the most stable bedforms and preserved the coarsest particles. Additionally, a large amount of gravel was deposited on runs or pool exit-point bars in bed stretches, in agreement with the results obtained by Thompson et al. [65] and Wohl [57] in steep bed-gravel streams.
In the ascending part of the hydrograph, the steepened water-surface gradients over alternate bars and runs in the upper RCR must have promoted bed dissection and removal of finer particles, which were then stored in pools downstream. This geomorphic process was described by Lisle and Hilton [66] when studying the fine bed material in pools of natural gravel-bed channels.
Downstream, in the middle RCR, new deposits at the head of the mid-bar referred to earlier contributed to raising its surface, mainly with sand and fine gravel (Figure 12(5)). The somital platform of the central bar was not totally submerged, and acted as a roughness element that must have locally increased the vertical velocity of the flow on the flanks, in turn driving the lateral vortex cells, with the consequent start or reactivation of entrainment of the less coarse particles. As the flow cells enlarged and moved locally, incision occurred in secondary active channels around the bar, increasing the bedload pulses downstream.
A similar contrast is provided by the comparison of the differences in absolute net volumes with the corresponding variations in elevation (negative or positive) detected with SfM and TLS. The total net deposition volumes were considerably higher than the erosion in all RGUs and in both RCR, according to the SfM data. Among these bedforms, within the upper RCR, the active bed registered the highest absolute volume increase (12 to 14 m 3 ), associated with bed raising of 0.2 to 0.3 m ( Figure 13). Active bars and unstable scarcely vegetated low bars also experienced a significant increase in volume, but with different ranges of heights. The largest net deposition volumes (from 4 to 8 m 3 ) meant increases in height from 0.15 to 0.30 m in the first case (leptokurtic distribution), and from 0.1 to 0.4 m in the second (dominant mesokurtic distribution). The high sparsely vegetated bars and vegetated tail only exhibited occasional variations that had little impact on the overall channel adjustments. The middle RCR showed a somewhat different pattern, characterized by Gaussian distributions and a single mode in most RGUs. Here, the active bed once again scored by far the highest volume charge, with maximums of 30 to 37 m 3 accompanied by the same bed rises (0.20 to 0.35 m) as in the upper RCR. The rest of the bedforms received, with respect to the upstream reference reach, a greater sedimentary contribution that did not translate into substantial morphological changes since the accretion height hardly varied. The reason for this could be the greater wetted channel width and the presence of transverse and central bars in the middle RCR, compared to the relative channel narrowness and development of longitudinal bars in the upper reach. The TLS results for the April event provided quite different patterns, with single-mode leptokurtic distributions in the most unstable bed forms: in a, b, and c of the upper RCR and in a and c of the middle RCR. The only exception was in the behavior of the active bars without vegetation, whose distribution showed two pointed curves in the central part, representing modes of erosion (average depth of lowering of −0.1 to −0.2 m) and deposition (average rise height of 0.20 to 0.35 m), respectively. This was probably due to unequal morphological adjustments, after transitory erosion occurred in both cases. Field observations verified lower rates of transitory erosion (only 5 cm for this event) in this morphosedimentary unit than in the active bed itself, where depths of 20 to 25 cm were reached (Figure 12(7)). In the rest of the RGUs, the net volumes of deposition exceeded those of erosion.
On a local scale, it is worth noting the bimodal distribution observed in the scarcely vegetated banks or tail, especially in the middle RCR, where the maximum volume of erosion supposed an average lowering of the bed surface of around 5 cm and that of deposition an average raising of 7 cm.
Consequently, all the patterns of spatial variability in the sediment budgets, described with SfM and TLS data, showed a clear tendency to increase bedload and granular bed armoring. The supplies of gravel in both RCR exceeded the volume of scour, causing an increasing alluvium that tends to make it more and more difficult to lower the bed surface by downcutting, increases the hydraulic resistance of the bed materials, and redirects the energy expenditure towards the banks. The result is a gradual widening of the channel and an increase in the width-to-depth ratio, particularly significant in sites of the middle and lower reaches, not subject to geological constrictions.

Conclusions
An excess or deficit of sediment in a sediment budget implies different responses of the channel morphology and bed forms. In the case of ephemeral gravel-bed streams, where the morphological adjustments are the product of complex flow dynamics at the event scale, the approach proposed here, based on the integrated use of UAV-SfM and TLS, has allowed this budget assessment to be satisfactorily addressed. Additionally, such an approach means economic savings and optimization of monitoring resources, through the synchronized application of both techniques, and the use of one or the other, depending on the magnitude of the event and the type of geomorphic change (overall channel adjustments or changes in bedforms). The SfM-MVS proved to be an appropriate technique for estimating sediment budgets in longer stream stretches, such as RCR over 100 m in length, while TLS provided excellent results for the determination of changes in bedload budgets at more detailed spatial scales (PBSA and RGU). Different complex morphodynamic processes, including total scouring-and-bedload transport and transitory erosion, which must have taken place simultaneously during each flood, could not be directly quantified. Nevertheless, the sediment budgets expressed in terms of net bed elevations and volumes were calculated with high accuracy for the wetted channel perimeter and the RGU during the monitored peak flows. This information was very useful, especially when coupled with in situ observations made after each event and the stream power parameters, obtained from the calibrated 1-D hydrodynamic model. The results from the study sites provided detailed information on the bed loading budgets, location of sediment sources and sinks, and changes in bedforms but also contributed to a better knowledge of the current morphodynamics of ephemeral gravel-bed channels in semi-arid Mediterranean coastal areas. With the support of field research, several hypotheses were verified, the most significant of which relates the magnitude of the events with the geomorphic process and morphological adjustments. Specifically, we found that flows exceeding bankfull discharge tend to produce vertical sedimentary accretion (0.20 to 0.35 m for a peak flow of 31 m 3 s −1 ), after having mobilized a large amount of bedload upstream. Great variations in the stream power gradient and high average excess energy (−15< ∂ω/∂s <15 Wm −3 ; mean ω/ω c ratios>2 for the same event) occurred during this process. By contrast, more moderate peak discharges (≤10 m 3 s −1 ) were only capable of exerting surface washing, selective transport, and local scouring, mainly affecting active low bars. The lowest net thickness of deposition detected by SfM and TLS in the upper stretch (less than 0.06 m) often masked an extraordinary mobility of the gravel deposits, crossing riffles and inter-bar runs, which could only be proved by measuring transitory erosion in the field (0.20 to 0.25 m at the upper RCR). Regardless of the magnitude of the monitored events, bed armoring was identified in situ, which tends to make the bed more resistant to erosion in future floods. If an increase in the frequency of flash floods due to climate change is confirmed in southeastern Spain, it is foreseeable that in the short and medium term, the trend of bed aggradation will continue in these ephemeral gravel streams, fed by the increasing inputs of coarse particles from the slopes, thus promoting lateral erosion and widening the channel, as has already occurred in more arid environments.

Conflicts of Interest:
The authors declare no conflict of interest. D 50 median grain size (m) D 84 particle size corresponding to the 84% of the sample weight (