Quantification of bedform dynamics and bedload sediment flux in sandy braided rivers from airborne and satellite imagery

Images from specially‐commissioned aeroplane sorties (manned aerial vehicle, MAV), repeat unmanned aerial vehicle (UAV) surveys, and Planet CubeSat satellites are used to quantify dune and bar dynamics in the sandy braided South Saskatchewan River, Canada. Structure‐from‐Motion (SfM) techniques and application of a depth‐brightness model are used to produce a series of Digital Surface Models (DSMs) at low and near‐bankfull flows. A number of technical and image processing challenges are described that arise from the application of SfM in dry and submerged environments. A model for best practice is presented and analysis suggests a depth‐brightness model approach can represent the different scales of bedforms present in sandy braided rivers with low‐turbidity and shallow (< 2 m deep) water.


Introduction
Over the past decade, considerable improvements in the quality, spatial scale, frequency of capture, and resolution of remote sensing imagery, have enabled new opportunities for investigating river morphodynamics (Lane et al., 2010;Javernick et al., 2014;Ishiguro et al., 2016;Vázquez-Tarrío et al., 2017). In particular, the use of unmanned aerial vehicles (UAVs) and Structure-from-Motion (SfM) photogrammetry has become increasingly popular (Carrivick et al., 2016;Kelleher et al., 2018). These techniques have been facilitated by advances in computer vision software (Fonstad et al., 2013), alongside the availability and low-cost of UAV systems. However, significant variations in data quality, both between and within surveys (Smith and Vericat, 2015) have been reported, stimulating a pressing need for more rigorous and confidencebounded data analysis methodologies (James et al., 2017).
The use of UAVs and aerial images has improved our ability to quantify the spatial organisation of river relief and bed roughness (Williams et al., 2014;Dietrich, 2016;James et al., 2017;Carbonneau et al., 2018), to track morphological change over short time periods (Lane et al., 1996(Lane et al., , 2010Palmsten et al., 2015) and to infer sediment transport rates (Lane et al., 1995;Brasington et al., 2003;Vericat et al., 2017). However, a number of outstanding challenges remain. One of the main issues with surveying rivers using photogrammetry is that unless the river is ephemeral, it is a two-media environment (Lane, 2000, Lane et al., 2010. Generating accurate elevations for inundated river beds has proved problematic using SfM photogrammetric software (Javernick et al., 2014;Woodget et al., 2015), a difficulty exacerbated in sandy river environments where well-sorted silts and sands may give poor image texture (Lane et al., 2010).
However, if these technological and analytical challenges can be overcome, these methods offer great potential for quantifying river process-form relationships. For example, the field measurement of bedload flux (defined here as bed material in continuous or intermittent contact with the bed) is notoriously difficult due to the logistical difficulties of collecting data at high stage and the inherent spatial variability in transport rates (Turowski et al., 2010;Frings and Vollmer, 2017). This represents a significant issue in terms of limiting the availability of data required to understand key links between flow, morphology and sediment transport rates, which are essential for parameterising and validating numerical models. To date, many estimates of bedload sediment transport are thus restricted to relatively small and shallow rivers where point samplers (i.e. Helley-Smith type samplers) or fixed location pit-type bedload traps have been used. In addition, more spatially extensive estimates of sediment fluxes have been generated using before and after flood topographic and lidar surveys (e.g., Goff and Ashmore, 1994;Brasington et al., 2003;Anderson and Pitlick, 2014). However, such repeat surveys may miss the record of successive erosion and deposition events and largely lack detail on how sediment transport varies through a flood (cf. Fuller et al., 2003), or how it might link to the short-term evolution of bed and barforms. The most detailed work on sediment flux has thus necessarily been limited to physical experiments where detailed topographic surveys of the entire bed are possible at high temporal and spatial resolution. While such experiments have demonstrated that sediment flux in braided rivers may correlate with the migration rate of bars within the channel (e.g., Wickert et al., 2013) or the cyclic erosion and filling of pools (Dhont and Ancey, 2018), this observation remains to be properly validated in the field. The new technologies described above now provide an opportunity to quantify the spatial distributions of bar and bedform migration rates in the field and thus advance knowledge of how flow, morphology and sediment transport are linked. This paper comprises two sections. The first details a number of methodological procedures for analysing imagery in a sandy braided river environment, while the second part illustrates how the resultant data can be used to monitor and to interpret sand-bed river dynamics over short time periods (hours to days). The objectives of this paper are thus to: (i) describe the processing steps and challenges associated with using aerial imagery and SfM photogrammetry to produce digital surface models (DSMs) in a sandy braided river with complex submergent and emergent relief; (ii) develop and refine a depth-brightness model to allow quantification of submerged elevations; (iii) use the resulting data to quantify the spatial variability of dune and bar migration rates; and (iv) assess the potential for using repeat aerial and satellite imagery of bed morphological change to map bedform and reach-scale sediment transport rates in a sand-bed braided river.

Study Site
The field-site detailed herein is the South Saskatchewan River, Canada, which is the location where the 'classic' facies model for sandy, braided rivers was developed as first proposed by Cant (1978) and Cant andWalker (1976, 1978). The river has a range of scales of different sand bedforms Sambrook Smith et al., 2006;Lunt et al., 2013), active bedload transport through most anabranches at a range of flow stages, and has been the subject of previous analysis both of planform change (Lane et al., 2010;Parker et al., 2013) but also sediment preservation and channel belt sedimentology (Woodward et al., 2003;Ashworth et al., 2011;Lunt et al., 2013). The low suspended sediment concentrations within the water column provide optimal conditions for assessing the potential for using aerial imagery to quantify bedform dynamics and morphological change.
The South Saskatchewan River flows from the Rocky Mountains in Alberta, Canada, into Lake Diefenbaker, 25 km upstream of the study site at Outlook, Saskatchewan ( Figure 1(A), (B)). Figure 1(C) shows the extent of the study site reported in this paper, which is divided into two main sections, upstream and downstream of the town of Outlook. Three smaller study reaches are reported in this paper termed 'SS1', 'SS2' and 'SS3' (see labels in Figure 1(C)). The Gardiner Dam traps much of the very fine sediment so that the downstream river flow is clear and therefore the river bed is entirely visible at even moderate flows. Previous work has shown that bed degradation downstream of the Gardiner Dam is minimal at the site of the study reaches (Galay et al., 1985;Phillips, 2003). The study reaches are ∽600 m wide, have an average bed slope of 0.0003 and possess a very well sorted medium sand bed (D 50 = 0.3 mm) with negligible clay. The channels are dominated by ripples and dunes with lobate unit bars, typically ∽1.5 m in height ).

Data Sources
Three types of airborne imagery were used to quantify and to monitor morphological change in the South Saskatchewan River: (i) specially-commissioned airplane sorties (manned aerial vehicle, termed here 'MAV') over an 18-km-long stretch of the river; (ii) UAV surveys over defined sub-reaches up to 1 km long; and (iii) satellite data from the Planet CubeSat constellation (Planet Team, 2017; https://www.planet.com). The following sections describe the process of image acquisition and data processing.

Aerial (MAV) photography to generate whole reach DSMs
Conventional aerial plane images (~0.06 m ground resolution) were captured at a height of~1500 m from a fixed-wing aeroplane with an UltraCamXp sensor for 2015, 2016 and two dates in 2017 (Table I). The flight lines for the 2015 images were single straight corridor flights down the river valley (similar to Lane et al., 2010), whereas those in 2016 and 2017 possessed 'back and forth' flight lines, which resulted in greater image overlap (see Table I). The greater overlap allows for better image matching in the initial processing of the SfM program Pix4D (https://pix4d.com/) that was employed throughout the analysis reported here.
Images are supplied with GPS information from the aircraft, along with Omega, Phi, and Kappa (rotations in the XYZ-axis respectively) to aid in the initial processing. In theory, ground control points (GCPs) are not necessary with such information, 954 R. J. P. STRICK ET AL. but analyses suggest that the SfM process tends to produce tilt and/or doming unless additional effort is made to constrain the SfM solution even with high-grade photogrammetric cameras (Bakker and Lane, 2017). A series of GCPs (crosses on 1 × 1 m targets with the centre of each cross occupying between 1 and 2 pixels) were therefore laid out on the floodplain and bar surfaces immediately prior to the overhead flight ( Figure 2). GCP locations were measured with a Leica 1230 real-time kinematic differential GPS, precise to ±0.02 m horizontally and ±0.03 m vertically. Fewer GCPs were surveyed in 2017 because the water level was much higher, resulting in less exposed bar surfaces on which to place targets.
Verifying bedform morphology from water depth measurements Water depth and bed height were quantified in 2015, 2016 and 2017, using a NAVISOUND 215 single beam echosonar (SBES) unit with a 200 kHz transducer deployed from a small rib boat.  The SBES operates at a vertical resolution of 0.01 m and a ping rate of 5 Hz, with the SBES surveys being located using the Leica 1230 dGPS system. The echosonar had a blanking distance of 0.15 m and was typically set 0.1 m below the water surface resulting in the minimum water depth that could be recorded as 0.25 m. SBES data were corrected for any dGPS errors associated with lost radio links and spikes in the depth data series, using a range of <0.3 m and >4 m. The SBES transects were conducted concurrent to the aerial sorties as well as at other times, across a variety of locations and depths in order to capture a range of different dune and bar morphologies present in the river.
Fixed-wing UAV image collection to quantify bedform dynamics Two fixed-wing unmanned aerial vehicles (UAVs) were used to collect repeat aerial imagery to quantify short-term bed evolution. Initial surveys in 2016 used the eBee RTK but this was replaced in 2017 by the eBee Plus. Flight lines were completed on a grid pattern. Trials demonstrated that an 80% lateral and longitudinal overlap of images was optimal for creating orthomosaics. Such large overlap values were required due to the extensive areas of water within the study reach (Javernick et al., 2014;Dietrich, 2016). Multiple flights were flown at 90 m altitude and 1 h intervals (maximum of five epochs in a day) to allow tracking of bedforms that could move at over 1 m h -1 (Table II). The overall quality of the images collected was affected by the presence of large areas of open water within the image capture zone and weather conditions at the time of flight. Excessive wind will both impact the ability of the UAV to stay on course (both eBee RTK and eBee Plus can operate in wind speeds up to 10 m s -1 ) and will also generate water surface waves that reduce visibility through the water column. The time of day and overall ambient light levels also affect the quality of the images due to surface shimmer and illumination. Table II summarises the UAV sorties used in this paper and the quality and number of epochs collected during the 2016 and 2017 field seasons. Ground Sampling Distance (GSD) was typically 0.025-0.03 m (Table II) depending on the flight parameters and subsequent processing steps. Figure 3 shows the spatial extent of the various UAV sorties in the 'downstream' reach ( Figure 1(C)).

CubeSat images
In order to compare the results obtained from the MAV and UAV tracking of barforms with independent satellite data, we conducted preliminary analysis of high-resolution (~3.5 m) 4-band imagery available from the Planet CubeSat constellation (Planet Team, 2017). This data source provides frequent (weekly to daily) imaging of much of the Earth's surface, and allows the tracking of features that show change over such timescales. Three epochs of satellite imagery were obtained from Planet for the period of the MAV and UAV surveys and where cloud cover was minimalon June 7th, June 12th and June 28th 2017. The images for the entire survey area were imported into Global Mapper where rectification was checked and adjusted where required using identical features within the image area. The position of unit bar fronts for the 18 km study reach was traced manually for each epoch. This yielded a sequence of 125 unit bars whose migration could be traced, and allowed the method of bar migration applied to the MAV and UAV imagery (see below) to also be applied to these CubeSat images.

Water level measurement
Water level in the study reaches was measured at 1 min intervals using multiple Solinst© pressure transducers (precise to ± 3 mm) mounted on dexion frames within the downstream reach and geolocated with dGPS. Corrections were made for atmospheric pressure via a Solinst© barometer. Figure 4 shows the water level variation during the period of data collection in the study reaches in June 2017. Daily water releases from the upstream Gardiner Dam for power generation during peak demand resulted in a diurnal water level fluctuation with maximum fluctuations of~0.2 m typically in the early morning hours each day.

Point Cloud and DSM Generation and Orthomosaic Production from Aerial Imagery
Images collected from the MAV sorties were used to produce DSMs (see Figure 5 for workflow). The processing protocols are described below together with an overview of the technical challenges and procedures adopted to overcome them.

Point cloud generation
Point clouds and orthomosaics were created in the Structurefrom-Motion (SfM) program Pix4D with the 'large frame' supplementary software package. In the first stage of processing, key points (tie points in classical photogrammetry) were extracted from across the images available. As the images already contained some geolocation information, the ground control point data were not added initially. A first bundle adjustment was undertaken using Pix4D's key point matching system optimised for aerial nadir imagery. Camera properties were specified using calibration certificates and were not allowed to vary during calculation. Once a solution had been obtained, the GCPs were added, using the WGS84/UTM zone 13 N (egm96) coordinate system. At this point, the precision of the bundle adjustments was compared with the theoretical precision defined by the image scale, and as the two were found to be of the same order of magnitude, the bundle adjustments were accepted. Then, the point cloud was derived. Matching was set to be to a lower density than the image resolution (25%) to avoid redundancy in determined elevations.

Tilt correction
Even with imagery acquired with a photogrammetric standard camera with precisely known focal distance, principal point offsets and lens distortion, and a reliable bundle adjustment, random error in the bundle adjustment can translate into systematic error, or tilt, in the DSM surface (Lane et al., 2004;Bakker and Lane, 2017). To check whether the DSMs contained residual tilt, elevation error values were derived by comparing the dGPS elevation values from the GCPs with the DSM elevation estimates at the GCP locations. Figure 6(A), (B) show the relationship between elevation error and Northing (approximately equivalent to distance along the valley axis) for two DSMs. Figure 6(A) shows a DSM that is free from systematic tilt error. Figure 6(B) shows a DSM with a systematic tilt error of c. 0.04 m km -1 in the Northing direction (assuming no significant tilt in the Eastward direction). To correct the DSM for the occurrence of systematic tilt, a manual and iterative approach is used to fit a 3D (X, Y, Z) modelled error plane (represented by red dots in Figure 6(A), (B)) to the elevation error values (blue dots in Figure 6(A), (B)). For the case where the DSM contains no systematic tilt, for example a fitted trend surface with no Northing dependence ( Figure 6(A)), then no correction is necessary. For the case where the DSM contains tilt such as a Northing dependence ( Figure 6(B)), the associated trend surface is used to correct the DSM. The example in Figure 6(A), (B) shows only the case of a potential Northing dependency, but it should be noted that the error surface is a 3D plane that could also be titled along the Easting.  Low image texture Bakker and Lane (2017) demonstrate that image texture is a critical control on the quality of SfM photogrammetry. Low image texture occurs when there is a homogenous distribution of pixel values with poorly identifiable distinct pixel regions. Low image texture can be a significant problem in sandy braided rivers where large areas may comprise water and low-relief sand bars (Lane et al., 2010). Some filtering of erroneous points is thus required to create a more reliable reconstruction of the mean elevation.
Points were filtered in MATLAB (point density was first reduced to 0.5 m) by applying a Chauvenet-type criterion (Chauvenet, 1960;Lane et al., 2004). Individual elevation (z) values were treated as noise and removed if where z and σ z are the mean and standard deviation of the point elevation values within a 35 pixel filtering window (17.5 m with the 0.5 m resampled grid) and m is equal to 2. Elevation points that were removed were replaced by the mean elevation value in a 3 × 3 pixel window centred on the point that was removed. The filtered point cloud was interpolated via kriging in ArcGIS to produce the initial DSM. The DSM was then classified into wet and dry areas so the former could be corrected for the refraction of water. A water surface was created by interpolating points from along bar and bank edges in the DSM. This water surface was used to identify (1) data points that are submerged and so need refraction correction; and (2) the water depth of those points. The refraction caused by water was then corrected by multiplying the newly created depths for each pixel by a refraction index of 1.34 (Westaway et al., 2000(Westaway et al., , 2001. These new depths were then subtracted from the water surface to create the corrected river bed elevations. The wet bed elevations were then merged with the dry DSM to create the overall initial DSM for the reach. Depth-Brightness Model in Aerial (MAV) Imagery The above method relies on the principle that it is possible to obtain a dense number of correct matches in inundated zones. As already shown for this kind of environment, this may be a challenge (Figure 7, see also Lane et al., 2010) even though the images themselves show a rich level of morphological detail the images themselves show a rich level of morphological detail (e.g. Figure 1(D)). For this reason, a second approach to depth determination was also employed, based upon depthbrightness relationships, that is used quite widely for deriving stream bed bathymetry from aerial imagery (Gilvear et al., 1995;Winterbottom and Gilvear, 1997;Westaway et al., 2003;Legleiter et al., 2004;Carbonneau et al., 2006;Marcus and Fonstad, 2008;Legleiter, 2016).
Application A depth-brightness model predicts flow depth from image pixel brightness (based on the absorption of light as it passes through the water; Gilvear et al., 1995). If light with an incoming intensity l in passes through a water depth of x, the remaining outgoing intensity l out can be estimated as where c is the rate of light absorption by water, which varies with turbidity and frequency of the incident light (Carbonneau et al., 2006). Using calibration data from the SBES and manual depth measurements, which were taken as close as possible to the time of image acquisition, the relationship between water depth and pixel log brightness was modelled using linear regression. Changes in water level between the times of image capture and the depth measurement surveys were accounted for by applying a linear depth correction, based on data from the water level recorders. Depth values were plotted in a GIS superimposed on a log brightness image, so that depth measurement locations could be matched precisely to the corresponding pixel value. Once the pixel values had been extracted, log brightness values were plotted against the corresponding depth values to create the depth-brightness relationship (Table III and Figure 8). The linear equation from this relationship was then used to convert the log brightness values from the orthomosaic into water depths.
The depth-brightness derived depths were subtracted from the water surface (created earlier) to produce river bed elevations, which were then merged with the DSM. This creates a combined DSM of SfM emergent elevations and depthbrightness submerged bed elevations.

Pixel colour saturation
One of the limitations of using a depth-brightness model is that for depths > 2 m colour saturation occurs in the image and depths can no longer be estimated accurately (Lane et al., 2010). Analysis of the extensive SBES data set from 2016 (747 756 points distributed over 18 km) shows that only 7.5% of values surveyed for the South Saskatchewan River were greater than 2 m depth ( Figure 9). However, these SBES data are biased towards deeper areas because of boat access and the need for rapid surveys of active dune fields. A better comparison is probably the proportion of saturated pixels in the total wetted area of the DSMs which is 1% in 2016 and 1.5% for 2017. This confirms that although this method of DSM production may not be suitable for deeper rivers, it works well to represent the majority of river bed elevations for the South Saskatchewan River.

Water surface reflection
An additional factor that can affect the quality of the depthbrightness model is surface reflection (Overstreet and Legleiter, 2017) that produces speckles of white on aerial imagery (see the white dots in Figure 10). These specks are converted into  large spikes in the depth-brightness model and subsequently derived DSMs. Problems can exist due to DSM-related mismatches that then cause incorrect local orthoimage production, typically associated with areas of expansive water.
This can leave poorly reconstructed areas in the orthomosaic, or even gaps (Figure 10), caused by Pix4D having difficulties in identifying matching points in inundated zones. The constant movement of the water surface also causes slight differences in the pixel values, as well as minor changes in illumination.
Removal of these artefacts from the bed elevations produced by the depth-brightness model was achieved by converting the bed elevation raster DSM into points. These points were then subjected to universal kriging with a constant order of trend removal, exponential Kernel, and a lag size of 1 with 31 lags. These parameters maintained the bed morphology, removed the large spikes caused by surface reflection and also filled any holes in the orthomosaics.
Accounting for pixel discolouration An additional challenge with using a depth-brightness model occurs where the spectral properties of the bed vary significantly, leading to a breakdown of the assumption that the optical signature of a point is only a function of water depth. For example, where there is a change in sediment type from sand to gravel, or there is a build-up of biofilm, there will be a change in spectral properties. This can result in darker pixels (associated with biofilms, for example) at similar depths to lighter pixels (associated with sand). The consequence of this is that the brightness model predicts greater depths for these darker areas than is the case in reality. To help overcome this problem, a difference map was generated between the depths from the DSM created using Pix4D and those from the depth- brightness model and any zones with a consistent difference of more than 1 m were masked. This resulted in 12% of masking in the 2016 wetted area portions of the DSM. Generally, these zones were in shallow water, which is where the original Pix4D DSM performed best (Figure 7). Therefore, the Pix4D DSM was laid on top of the brightness model in these zones. This issue was particularly prevalent in the 2016 images when biofilm growth was extensive following a warm spring period.

Shadow effects
Depending on the time of image acquisition, there is a possibility that shadows behind topography (e.g. dune crests, bar fronts, boulders) can cause pixel discolouration (darkening) that then translates into a perceived 'over deepening' of the bed. There was thus a trade-off between flying early or late in the day (low sun angle but more shadows) versus near midday (high sun angle but more surface glare). All aeroplane (MAV) images (that were used for DSM construction) were acquired around 10.00 and 15.00 with a low sun angle. Comparison of echosonar surveys taken immediately after the MAV image acquisition with DSMs from the MAV aerial orthomosaics shows the shadowing, and thus potential overdeepening, is a maximum of 30% (worse-case scenario) and averages 10-20%, but it only impacts 1-2 pixels adjacent to the crest of topographic highs.
Error analysis and validation for the plane imagery DSM error analysis The following error analysis is for the DSMs created from the MAV imagery. It compares the z values from field dGPS surveys from exposed/emergent topography, on or 24 hours after the aerial photography flight, with the z values for the same location as estimated from the DSMs (n = 13-76 for the different reaches and epochs, see Table IV). Although the root mean square error (RMSE) is reported for comparison with other studies, the standard deviation of error (SDE) is used for the main error analysis as RMSE conflates systematic (mean) and random (standard deviation) errors, which are different parameters. Table IV shows that the SDE ranged from ±0.08 m to ±0.19 m for the different DSMs and SDE/GSD ranged from ±1.33 to ±3.17, which is very similar to the range of values reported for other studies (see Table V).

Depth-brightness model evaluation
The performance of the depth-brightness model was evaluated by comparison with equivalent SBES and dGPS data. Comparing the average dune heights measured from the SBES data with those measured from the DSM, showed no significant difference (Paired t-test, P = 0.078) between the two different methods, giving confidence in the ability of the brightness model to recreate accurately the bedform height. Additionally, dGPS measured depth points were compared with those estimated from the brightness-model DSM in the downstream reach for June 8th 2017 ( Figure 11). The SDE for the brightness DSM points was ±0.17 m, which is similar to the overall DSM SDE value of ±0.12 m (Table IV).

DOD error analysis
When comparing two different DSMs, it is critical to understand the survey uncertainties and potential error propagation between the two surfaces (James et al., 2017). Additionally, when using DSMs of difference (DoDs), changes smaller than a specified 'level of detection' (LoD) should be omitted from analysis. The LoD can be defined using : where σ Z1 and σ Z2 are the vertical standard deviations of error of the two DSMs used to generate the DoD, and t is an appropriate value for the required confidence level, typically 95% (1.96). A single LoD value (Table IV) can then be estimated from SDE. However, depending on the distribution of points used in the SDE, significant smaller changes that are spatially coherent over large areas may be neglected (Brasington et al., 2003). The errors reported for the different DSMs used in the present paper (see SDE/GSD, Table IV) compare favourably with other studies (Table V), bearing in mind the different image scales.

Measurement of Unit Bar and Dune Migration and Estimates of Sediment Transport Rates from MAV and UAV Imagery
Unit bar sediment transport rates from MAV imagery The orthomosaics created from the MAV images were used to trace the migration of all unit bar avalanche faces for the entire reach between June 8th and June 12th 2017 (n = 295). Each visible unit bar avalanche face was digitised manually to produce a bar front length (L) for each flight date. A polygon was then created by joining the bar front positions from each date together and a planform area (A b ) of migration was calculated. The migration rate for each bar between epochs t1 and t2 (M b in m h -1 ) was then calculated by: This methodology was also applied to the CubeSat satellite imagery.
In addition to quantifying unit bar migration rates, the average unit bar migration direction was computed as the bearing (from north) of the centroids of each bar crest as positioned in the images from 8th and 12th June 2017.
For the MAV imagery sediment transport rates of individual unit bars were calculated with reference to the DoD produced for 8-12th of June 2017. The DoD provided a volume of difference (V) between the two surfaces contained within each polygon area. Sediment transport rate per unit width (Q sb , m 2 h -1 ) was calculated from unit bar migration as: where P is the sediment porosity and is taken as 0.4 for well sorted sands (van Rijn, 1993).
Dune migration rates and estimated sediment transport rates from UAV imagery Dune migration rates were quantified from repeat UAV flights in reaches SS1 and SS2 (2016) and SS3 (2017) (see Figure 3 (A), (B)). Dune migration was quantified for 'patches' in each reach that varied in size from 340 to 5198 m 2 (see Table VI). Patch locations were selected to provide a spatially-distributed data set that represented the different morphological elements and a range of local flow depths. Individual dune migration rates were calculated using the dune crestlines identified in each georectified orthomosaic and the same calculation procedure described above for unit bars. Dune sediment transport rates per unit width (Q sd , m 2 h -1 ) were calculated for individual dunes in each of the 20 patches of the SS3 2017 reach only (Figure 3(B), Table VI) as: where b is dune shape factor calculated as A d /(H x λ), A d is dune area calculated between the two successive troughs, H is dune height (from trough to crest), λ is dune wavelength (from trough to trough), and M d is dune migration rate calculated as for unit bars.
Values of H and b were obtained for each dune from downstream profiles taken through individual dunes using the brightness model DSM from the MAV images taken on 12th June 2017. Dunes identified in the UAV images could be tied to the same dunes in the MAV orthophoto/DSM because the MAV sortie was flown only 2 hours earlier. Typically, the morphology of between 9 and 20 dunes was measured in each patch depending on patch size and dune prevalence.

Using Aerial Imagery to Understand the Dynamics of Sand-bed Braided Rivers
Sand-bed braided river morphology and change Despite the challenges of using airborne imagery to visualise and quantify river bed topography, Figure 12 shows that the resulting DSMs are capable of representing well the complex and multiple scales of bed topography that characterise sandy braided rivers (Cant, 1978;Skelly et al., 2003;Sambrook Smith et al., 2006;Horn et al., 2012). For example, the DSM ( Figure 12) displays many small (label A), medium (label B) and large (label C) 2D and 3D dunes, incipient unit bars (label D), exposed unit bars with elongated bar tails (label E), and early-stage bar top hollows (labels F). With a GSD of 0.06 m for the aerial imagery it is impossible to resolve bedforms at the scale of ripples but it is possible to observe small superimposed dunes (label G). Even though Figure 12 is just one snapshot in time, the resolution and therefore clarity of the DSM, provides valuable insight into the interrelationship between different bedforms including the spatial transformation from 2D to 3D dunes and the progressive response of dune planform shape and wavelength to a change in water depth and unit bar emergence. The DSM (Figure 12) also illustrates that the South Saskatchewan River can have some reaches that contain few dominant main channels (although see discussion below). Except for the 20-40-m-wide channel on the right of the reach (Figure 12, near label B), where there is ongoing bank erosion, the DSM is characterised by a series of migrating and stacked unit bars. Figure 13(A) shows the morphology of the South Saskatchewan River at a different scale to that shown in Figure 12. Unlike Figure 12, the~6-km-long reach at low flow shows a dominant main channel that is mostly pinned against the west (left) bank. The dominant thalweg splits and re-joins both around exposed unit bars and major compound bars. The DoD over a period of 15 months (Figure 13(B)) shows examples of bar head erosion (labelled H), lateral bar accretion of up to 2 m (labelled M), thalweg scour up to 4 m (labelled T), and anabranch abandonment and filling (labelled F) (cf. Lane et al., 2010;Parker et al., 2013). The majority of the compound bars, and particularly those that are attached to the right bank (east), have not aggraded appreciably and have been simply eroded at their upstream ends, with deposition at their tail as the bars migrate slowly downstream (cf. Bristow, 1987;Ashworth et al., 2000;Best et al., 2003;Nicholas, 2013;Schuurman and Kleinhans, 2015).
The average distribution of bed heights from the DSMs for the upstream reach (Figure 13(A)) for the 3 years of aerial imagery (2015)(2016)(2017) shows negligible change (Figure 14) with an overall near-normal distribution and annual bed height fluctuations between lower (2015) and higher (2017) average bed elevations. The annual reach channel change (Figure 13(B)) can be contrasted with the DoD for the 4-day period in June 2017 when discharge was relatively high for the year (Figure 13(C) and 13(D)). There is little topographic change within this short time period (represented by white colour shading in Figure 13(D)), but the DoD does clearly pick out the    Figures 13(B), 13(D) and 14 suggest the aerial imagery and processing is capturing and resolving both the short and longer-term morphological changes in the dynamic South Saskatchewan River.

Spatial variability in dune and unit bar migration rates
Based on the UAV imagery from three study reaches ( Figure 15, Table VI) dune migration rates range over an order of magnitude from 0.1 m h -1 to 1.0 m h -1 . These rates are towards the lower end of rates reported in previous work, for example 0.8 m h -1 for the Fraser River, Canada (Villard and Church, 2005), 0.8-1.6 m h -1 for the Calamus River, USA (Gabel, 1993), 1-2.7 m h -1 in the Gels A, Denmark (Kisling-Moller, 1992), but up to~100 m h -1 for highly mobile gravel dunes in the Toutle River, USA (Dinehart, 1989). The primary reachscale morphological control on dune migration rate appears to be the extent to which flow is either largely contained within a single channel or more widely spread across the braidplain. This is most clearly illustrated by comparison between reaches SS2 and SS3 (Figure 15(B), (C)). For reach SS2, flow was largely constrained within a single~150-m-wide channel. In particular, the downstream section of this channel was uniform in width and cross-sectional morphology (encompassing the last four measurement patches) resulting in very similar dune migration rates that only varied between 0.43 m h -1 and 0.53 m h -1 . In contrast, reach SS3 was constrained within one channel at its most upstream point before changing to a braidplain > 500 m wide with multiple unit and compound bars. As the flow is distributed across a much wider area, dune migration rates diminished from 1.15 m h -1 at the head of the reach to a minimum of 0.13 m h -1 downstream of a sheltered bar tail. Dune fields that approach actively accreting bar heads were also zones of high dune migration rates as illustrated by reaches SS1 (Figure 15(A)) and SS3 (Figure 15(C)).
There was a difference in reach flow discharge when the dune migration data were measured for the reaches shown in Figure 15 (A), (C). For reaches SS1 and SS2, the flow discharge through the reach was~65 m 3 s -1 , whereas for reach SS3 it was~275 m 3 s -1 . Despite this four-fold difference in local discharge, the range in dune migration rate was broadly similar for all reaches, suggesting that: (a) increases in discharge simply spread water further across the braidplain, rather than further concentrating flow in the deeper thalwegs, which might increase dune migration rates; and (b) the main driver of the spatial changes in dune migration rates is the local topography and multiple scales of roughness. Dune size has a finite range of variation in these limited water depths that are similar between reaches, and that do not change appreciably once the flow has reached bar-top level. Consequently, this limits the range of dune sizes, and thus migration velocities, that are broadly similar between reaches.
The migration rates for unit bars (n = 295) are typically an order of magnitude lower than those of dunes, with an overall average of 0.06 m h -1 . However, as with the dune migration rates discussed above, variability exists between the different bars studied, with a maximum of 0.32 m h -1 . The spatial pattern of migration rates, and therefore sediment transport rates for a given bar height, throughout the river is characterised by high variability (Figure 16), although some trends are apparent. The lowest migration rates are either largely found on the braidplain margins away from the main channels or in compound bar tail areas, similar to the dunes discussed above. The migration rates of unit bars also vary within different channels that flow around a compound bar, presumably in response to flow steering around the bar. For example, a central vegetated compound bar in Figure 16(B) has greater rates of unit bar migration down the east (right) compared with the west (left) channel. Likewise, unit bar migration rates are often lower within cross-bar channels compared with the main thalwegs (e.g. lower part of Figure 16(C)). The migration direction of unit bars is also highly variable (Figure 16(D)). Unit bars typically migrate in directions ±90°to the 'downstream' direction of the valley (shown as a shaded zone on Figure 16(D)) and the main channels. There is some indication (Figure 16(D)) that these migration rates diminish when the unit bar orientation deviates more than 30°from the overall downstream direction, perhaps reflecting the lower competence of cross-braidplain channels. These quantitative measurements highlight the spatial complexity of sediment transport within sandy braided rivers, and suggest that local factors such as changing flow depths and flow steering, are likely critical in determining the spatial distribution of sediment transport rates.
Bedload sediment transport rates calculated from dune and unit bar migration For reach SS3, 20 patches ( Figure 3B), comprising a total of 255 dunes were monitored to calculate local bedload sediment transport rates using the morphological measurements of each individual dune (see Equation (6)). Average bed sediment transport rate per unit width based on migration of these dunes was 0.016 m 2 h -1 , with rates varying between 0.001 m 2 h -1 and 0.134 m 2 h -1 . Within the same reach, and over approximately the same time period, 19 unit bars were also used to estimate the bed sediment transport rate (see Equation (5)), yielding an average of 0.024 m 2 h -1 , with rates varying between <0.001 m 2 h -1 and 0.100 m 2 h -1 . These dune and unit bar derived rates are thus very similar and demonstrate that both methods return a bed sediment transport rate per unit width of 0.02 m 2 h -1 . In addition, the unit bar bedload transport rates were also calculated over the entire study reach, totalling 295 bars, yielding a value of 0.02 m 2 h -1 , which is identical to that for the 19 unit bars used in the analysis of reach SS3. However, a key difference between the dune and unit bar bedload transport rates is in their distributions ( Figure 17) that are positively skewed but bimodal for dune and unit bar bedload transport rates respectively. Unit bar bedload transport rates have~29% of values <0.0025 m 2 h -1 and~17% of values >0.05 m 2 h -1 .
There is also significant variability in bedload transport rates within reach SS3, with a tendency for transport rates (calculated from dune migration) to be lower downstream of unit bar fronts compared with rates calculated from dunes migrating on the stoss slope of unit bars. This attribute is most clearly illustrated by comparing patches 5, 17, 18 and 8 (see labels in Figure 15(C)) that are on the stoss side of a classic lobate fronted unit bar, and patch 10 ( Figure 15(C)) that is downstream of the unit bar front. The stoss side patch of dunes yields an average bed sediment transport rate of 0.015 m 2 h -1 , while the leeside value is only 0.007 m 2 h -1 . Likewise patch 9 (Figure 15(C)), also downstream of a lobate fronted unit bar,

The Potential of Using Aerial Imagery to Quantify the Dynamics of Sand-bed Braided Rivers
Previous studies that have used UAV imagery coupled with SfM have highlighted the potential problems in generating DSMs from inundated parts of the channel and highlighted water depth, turbidity and turbulence as key limiting factors. For example, Woodget et al. (2017) studied a small 3-13-m-wide meandering gravel-bed river, and suggested a limit of 1.4 m depth for the methodology. Based on the addition of using a depth-brightness model coupled with SfM, for low-turbidity waters, we suggest herein that this upper limit can be extended to 2 m depth and for conventional MAV imagery as well as UAV drone imagery. It should be noted that application of SfM, together with a depth-brightness model, works well in the South Saskatchewan River where suspended sediment loads are low due to the trapping of fine sediment by the upstream Gardiner Dam. Bedform tracking in much deeper channels and higher turbidity (e.g. in many of the world's largest rivers), may be more challenging unless there are future significant technological advances (e.g. UAVs tethered to miniaturised, boat-mounted, echo-sounders as illustrated by Alvarez et al. (2018);Bandini et al. (2018)).
One advantage of using UAVs that has been demonstrated herein is that because they are relatively simple and inexpensive to deploy, they provide an excellent tool for elucidating the temporal morphodynamics of rivers. To date, most studies have highlighted the improvement in spatial resolution that UAV imagery provides, such as for habitat mapping (Woodget et al., 2017) or detection of grain size variation (Carbonneau et al., 2018). The results described herein for a sand-bed river, where sediment is being transported at most flow stages, demonstrates unambiguously that in optimal conditions, bedform dynamics can be quantified at hourly timescales and at both low and near-bankfull stages.
While the concept of dune tracking to estimate bedload transport rates has a long history, this has traditionally been achieved using repeat bed profiles derived from boat-based echosounder surveys (van den Berg, 1987;Ten Brinke et al., 1999;Dinehart, 2002;Gaeuman and Jacobson, 2007;Claude et al., 2012) or at-a-point with fixed depth profilers (Dinehart, 1989). The results described herein demonstrate that in shallow, low-turbidity rivers, this approach also lends itself to image analysis methods without the need for boat-based surveys and the logistical challenges that involves. The ease with which drones can be deployed also means that even in environments where dunes migrate rapidly, such as the South Saskatchewan River, repeat surveys can be undertaken with sufficient frequency to ensure dunes can be identified and tracked from one aerial survey epoch to the next, before they lose their morphological coherence and visually recognisable form.
Although bedform tracking using repeat airborne imagery can yield insightful quantification of the bedload sediment transport rate it should be noted these estimates represent a minimum value. This is because not all sediment in transport contributes to the migration of the dune or unit bar and some sediment is likely steered around the bedform or temporarily passes into suspension over the lee face (see discussion for the case of dunes in Mohrig and Smith (1996)). However, in bedload dominated streams, such as gravel-bed rivers and the South Saskatchewan River, this underestimation of bedload sediment transport rate (particles in continuous or intermittent contact with the bed) is probably minimal.
The similarity between the mean bedload transport rates derived from both dune and unit bar migration suggests that most of the sediment being transported by unit bars is being supplied by dune migration over the stoss sides of the unit bars. A similar conclusion was reached by Villard and Church (2005) in their study of bar and dune dynamics in the Fraser River, Canada. This suggests that consideration of just unit bar migration rates may enable estimation of bedload flux at the scale of the channel width. In this case, a key consideration is the minimum area of study required to generate a robust sediment transport estimate. Reach SS3 in the present study, which was 685 m × 1190 m and contained 19 unit bars with measurable migration rates, provides an indication of the likely minimum area required. For example, the bedload sediment transport rate derived from the 19 unit bars in reach SS3 was the same as for the 295 bars of the entire study area. This suggests that the measurement reach was sufficiently large to capture the range of transport rates found throughout the 18 km study area. As a guide, a reach length of~3 x the braidplain width may be a good approximation for the monitoring area required for a first-order estimate of bedload sediment transport rate in such rivers. Whether this rule of thumb of optimal reach size is transferable to other rivers requires further testing and may depend on the local complexity of braidplain topography and prevalence of unit bars.
A final question is whether high-resolution DSMs and DoDs are actually required to generate a robust average bedload sediment transport rate for a reach in a river like the sandy braided South Saskatchewan. Generating simple rectified orthomosaics from drone imagery is much quicker than generating DSMs, especially for the submerged regions. In essence, can the sediment transport rate be generated from just 2D imagery using the unit bar migration rate? To explore this question, migration rates for the 19 unit bars of SS3 reach were all assigned the same shape factor (0.5, although it is recognised that bars may have slightly higher values depending on their aggradational history and migration rate), the same porosity (0.4 for sands) and given a bar height value of 1.25 m (an average value reported for the South Saskatchewan River by Lunt et al. (2013)). These bar attributes together with the measured migration rate then allowed calculation of the bedload sediment transport rate. The results reveal the average for the measurement reach was exactly the same as that for where the shape and height of each unit bar was individually established based on the DSMs. While this observation clearly needs additional testing to assess its utility, it could revolutionise our ability to estimate reach-based bedload sediment transport rates.
To provide a preliminary test of this idea, we used three epochs for June 7th, 12th and 28th 2017 from the 3.5 m resolution satellite imagery of the South Saskatchewan River provided by the CubeSat constellation of Planet.com (Planet Team, 2017). This imagery ( Figure 18) clearly shows the position of unit bars within the river and allows their crestlines to be tracked between epochs (see colour lines on Figure 18(A), (B) and (C)), although dunes cannot be discerned. The migration rates obtained from these images (Figure 18(D)) over a 21-day time period show a range of values comparable with the UAV and MAV imagery, with mean migration rates of 0.062 and 0.088 m h -1 for the UAV/MAV and Planet images respectively (both with standard deviations of 0.055 m h -1 ). The distribution of migration rates is more skewed to smaller values for the UAV/MAV estimates (Figure 18(D)), probably due to the difficulty of delineating smaller unit bars in the coarser-resolution satellite imagery. Although dunes in the South Saskatchewan River cannot be resolved in these CubeSat images, such bedforms may be resolved in larger rivers (so long as turbidity levels are low), where dune size may range up to many metres in height and tens to hundreds of metres in wavelength (see Galeazzi et al., 2018 for details of such dunes in the Amazon River). In this case, dune migration rates may also be quantifiable, although the results herein suggest that estimation of solely the unit bar migration rates may yield comparable results.
In addition, the pixel resolution of the satellite imagery also dictates the period between epochs that will allow slower migrating bars to be quantified. To demonstrate this contention, the distribution of migration rates for short (5 days) and long (21 days) temporal gaps between Planet images (Figure 18(E) inset) shows that, at this spatial resolution (3.5 m), slower migration rates can only be quantified at longer temporal spacings. The shorter temporal period of 5 days is unable to capture migration rates of less than~0.06 m h -1 (see grey shaded area in Figure 18(E)).
The utility of using such frequent satellite imagery thus relies on: (i) cloud-free skies; (ii) the spatial resolution of the sensor with respect to both the epoch frequency and bar migration rates, and (iii) detection of unit bar fronts in flows that may be deep and turbid. However, if these conditions can be addressed satisfactorily, or technological advances can solve Figure 18. (A)-(C) Example Planet.com images of the South Saskatchewan River acquired in the same period as the specially-commissioned aerial (MAV) flight sorties together with tracing of unit bar crests; (D) calculated bar migration rates (see text for methodology) for unit bars using the Planet. com satellite imagery compared with measurements from MAV aerial imagery. (E) Inset shows bar migration rates derived from the Planet.com images for two different temporal gaps between images, illustrating the shorter (5-day) interval is unable to resolve migration rates less than~0.06 m h -1 .
Determination of migration rates must thus be considered with respect to pixel resolution and temporal frequency of images. [Colour figure can be viewed at wileyonlinelibrary.com] 969 BEDFORM DYNAMICS AND BEDLOAD SEDIMENT FLUX IN SANDY BRAIDED RIVERS these challenges, such frequent satellite imagery opens up the possibility of estimating bedload sediment transport rates for inaccessible rivers, both large and small, where there is currently a dearth of data. This may improve our ability to estimate bedload sediment flux in rivers across the globe, and its spatial and temporal variation. Such a data requirement is central to addressing issues of channel change, anthropogenic influences on sediment transport and channel morphology, validating and improving numerical models of longer-term channel change and helping guide river channel management (Best, 2018).

Conclusions
Four principal conclusions can be drawn from the work reported herein. First, airplane (MAV) and drone (UAV) aerial imagery can be used with Structure-from-Motion (SfM) processing techniques, and application of a depth-brightness model, to generate robust DSMs for sandy braided rivers at flows near bankfull and with depths up to 2 m. The importance of following best practice in the application of SfM techniques is pivotal here, including collecting ground control data, checking for residual tilt in the analysis, and suitable filtering of acquired data. In addition, the present study confirms the value of using the optical richness of underwater zones to produce extremely detailed topographic data, and proposes a series of solutions to address the problems that can arise in doing so (e.g. reflection at the water surface, spatial variability in bottom reflection). Second, dune and bar migration rates and directions in sandy braided rivers are highly variable spatially, with evidence for local planform morphology determining some of this variability. Third, bedload sediment transport estimated from unit bar migration rates broadly integrates the sediment supplied by the dunes migrating over their stoss sides, such that reach-based estimates of bedload sediment transport rates can be estimated from unit bar data alone. Fourth, reasonable estimates of bedload sediment transport may be possible from just unit bar migration as measured from 2D imagery, including from satellites, and assumptions regarding average bar shape and height. With further technological advances, such a methodology could revolutionise our ability to provide global estimates of bedload sediment flux from large and/or inaccessible rivers that are currently very poorly quantified.